Responses separate from ModelSpec; wildcard(exclude=) takes a response

CHANGELOG.md

@@ -4,11 +4,11 @@
 
 ### New features
 
-- **Responses:** `ModelSpec(terms, response=...)` and
-  `spec.build_response(X)`, the response built from the same data as the
-  design. A response is a column, a list of columns, or
+- **Responses**, built from the same data as the design, as a
+  `ModelSpec` builds it: `response.fit(D)`, then `response.transform(D)`
+  (or `fit_transform`), in the backend of `D` (e.g. polars).
   - `Surv(time, status)` or `Surv(start, stop, status)`, as R's: a
-    `pd.DataFrame` of these columns, as glmstar's `CoxFamily` or lifelines
+    dataframe of these columns, as glmstar's `CoxFamily` or lifelines
     take it;
   - `Binomial('y')` for a binary column (0/1, boolean or two levels; the
     second level, or `level=`, is a success), or `Binomial(successes,
@@ -18,20 +18,27 @@
     one-hot (`output='onehot'`), or `Multinomial([...])` for counts.
 
   Responses learn their levels at fit; a level unseen at fit raises.
-  The response's columns cannot be in a term.
-- **`wildcard`**, R's `.`: every column of `X` but the response, so
-  `ModelSpec([wildcard], response='y')` is `y ~ .`. `wildcard(pow=2)` is
-  `.^2`, `wildcard(exclude=['x'])` is `. - x`, `(wildcard, 'g')` is `.:g`,
-  and `pca(wildcard, n_components=3)` a feature of every column. It is
-  expanded at fit into `expanded_terms_`; columns that are terms already
-  are not repeated, as in R's `y ~ x + .`.
+- **`wildcard`**, R's `.`: every column of `X`, but those of
+  `wildcard(exclude=...)`, so `ModelSpec([wildcard(exclude='y')])` is the
+  design of `y ~ .`. `exclude` also takes a response, leaving out its
+  columns: `wildcard(exclude=Surv('time', 'status'))`. `wildcard(pow=2)`
+  is `.^2`, `(wildcard, 'g')` is `.:g`, and `pca(wildcard(exclude='y'),
+  n_components=3)` a feature of the columns. It is expanded at fit into
+  `expanded_terms_`; columns that are terms already are not repeated, as
+  in R's `y ~ x + .`. `build_submodel` expands it too.
 - `pca()`'s `name` is optional.
 - `fit(X, y)` checks that the indices of `X` and `y` agree when both are
   pandas objects.
 
+### Bug fixes
+
+- `poly`, `ns` and `bs` are fit on the observed values of their variable,
+  and are missing only in its missing rows: before, a single missing value
+  made every value of the feature missing (or, for `ns`, arbitrary).
+
 ### Documentation
 
-- A page on responses and the wildcard.
+- A page on responses and the wildcard, and on missing values.
 
 ## 0.1.0 (2026-10-26)
 

docs/api/responses.rst

@@ -8,4 +8,4 @@ Responses
 .. autoclass:: modelmatrix.Binomial
 .. autoclass:: modelmatrix.Multinomial
 .. autoclass:: modelmatrix.Response
-   :members: fit, transform
+   :members: fit, transform, fit_transform

docs/index.md

@@ -14,7 +14,7 @@ pip install modelmatrix
 
 - [Building design matrices with `ModelSpec`](spec.md)
 - Features: [polynomials](poly.md), [splines](splines.md) and [PCA](PCA.md)
-- [Responses and the wildcard](responses.md): `y ~ .`, `Surv`, `Binomial` and `Multinomial`
+- [Responses and the wildcard](responses.md): `y ~ .`, `Surv`, `Binomial`, `Multinomial` and missing values
 - [ANOVA](anova.md) with `build_sequence`
 - [`ModelSpec` compared](comparison.md) with `ColumnTransformer`, patsy and formulaic
 - The [API reference](https://intro-stat-learning.github.io/modelmatrix/api/)

docs/responses.md

@@ -13,11 +13,12 @@ kernelspec:
 
 # Responses and the wildcard
 
-A `ModelSpec` can also describe the response of a model,
-`ModelSpec(terms, response=...)`, and `build_response` then builds it
-from the same data as the design. With a response,
-`wildcard` stands for every other column, as `.` in R's
-formulae: `ModelSpec([wildcard], response='Sales')` is
+As in scikit-learn, the design `X` and the response `y` are built
+separately: a `ModelSpec` builds the design, and a response, such as
+`Surv`, `Binomial` or `Multinomial`, builds `y` from the same data.
+`wildcard` stands for every column, as `.` in R's formulae, and
+`wildcard(exclude=...)` leaves out the response:
+`ModelSpec([wildcard(exclude='Sales')])` is the design of R's
 `Sales ~ .`.
 
 ```{code-cell} ipython3
@@ -42,8 +43,8 @@ Carseats = load_data('Carseats')
 ## `y ~ .`
 
 ```{code-cell} ipython3
-spec = MS([wildcard], response='Sales', intercept=True).fit(Carseats)
-X, y = spec.transform(Carseats), spec.build_response(Carseats)
+spec = MS([wildcard(exclude='Sales')], intercept=True).fit(Carseats)
+X, y = spec.transform(Carseats), Carseats['Sales']
 sm.OLS(y, X).fit().params
 ```
 
@@ -55,11 +56,11 @@ spec.terms, spec.expanded_terms_
 ```
 
 Columns that are terms themselves are not repeated, as in R's
-`Sales ~ Price + .`, and `wildcard(exclude=...)` is R's `. - x`:
+`Sales ~ Price + .`, and more columns can be excluded, as R's `. - x`:
 
 ```{code-cell} ipython3
-MS([poly('Price', 2), 'Price', wildcard(exclude=['Population', 'Education'])],
-   response='Sales').fit(Carseats).expanded_terms_
+MS([poly('Price', 2), 'Price',
+    wildcard(exclude=['Sales', 'Population', 'Education'])]).fit(Carseats).expanded_terms_
 ```
 
 ## `y ~ .^2` and `.:x`
@@ -69,11 +70,11 @@ and the wildcard can be a part of an interaction, as R's `.:US`:
 
 ```{code-cell} ipython3
 small = Carseats[['Sales', 'Price', 'Income', 'US']]
-MS([wildcard(pow=2)], response='Sales').fit(small).expanded_terms_
+MS([wildcard(pow=2, exclude='Sales')]).fit(small).expanded_terms_
 ```
 
 ```{code-cell} ipython3
-MS(['US', (wildcard, 'US')], response='Sales').fit(small).expanded_terms_
+MS(['US', (wildcard(exclude=['Sales', 'US']), 'US')]).fit(small).expanded_terms_
 ```
 
 ## Features of the wildcard
@@ -82,26 +83,33 @@ A wildcard can be a variable of a feature: here, the first principal
 components of the numeric columns, but the response.
 
 ```{code-cell} ipython3
-pc = pca(wildcard(exclude=['ShelveLoc', 'Urban', 'US']), n_components=3, scale=True)
-spec = MS([pc, 'ShelveLoc'], response='Sales').fit(Carseats)
+pc = pca(wildcard(exclude=['Sales', 'ShelveLoc', 'Urban', 'US']), n_components=3, scale=True)
+spec = MS([pc, 'ShelveLoc']).fit(Carseats)
 spec.transform(Carseats).iloc[:5]
 ```
 
 Its parameters are tuned as those of any term:
 
 ```{code-cell} ipython3
-y = spec.build_response(Carseats)
-pipe = make_pipeline(MS([pc, 'ShelveLoc'], response='Sales'), LinearRegression())
+pipe = make_pipeline(MS([pc, 'ShelveLoc']), LinearRegression())
 grid = GridSearchCV(pipe, {'modelspec__terms__0__n_components': [1, 3, 5, 7]}, cv=5)
-grid.fit(Carseats, y)
+grid.fit(Carseats, Carseats['Sales'])
 grid.best_params_
 ```
 
-The response is passed to `fit` as `y` (a `Pipeline` passes the design's
-data only to the `ModelSpec`). When both are pandas objects, `ModelSpec.fit`
-checks that their indices agree.
+When `X` and `y` are both pandas objects, `ModelSpec.fit` checks that
+their indices agree.
 
-## Survival responses: `Surv`
+## Responses
+
+A response is fit, e.g. to learn the levels of a categorical, so that
+it is coded the same way on new data, and `transform` builds it:
+`response.fit_transform(D)`. For a dataframe of another backend, such
+as polars, the response is in that backend. `wildcard(exclude=response)`
+leaves out the response's columns, so that it is not a term of its own
+design.
+
+### Survival responses: `Surv`
 
 `Surv(time, status)` is R's survival response, and `Surv(start, stop,
 status)` that of (start, stop] intervals. The response is a
@@ -110,23 +118,23 @@ status)` that of (start, stop] intervals. The response is a
 
 ```{code-cell} ipython3
 BrainCancer = load_data('BrainCancer').dropna()
-spec = MS([wildcard], response=Surv('time', 'status')).fit(BrainCancer)
-X, y = spec.transform(BrainCancer), spec.build_response(BrainCancer)
+response = Surv('time', 'status')
+y = response.fit_transform(BrainCancer)
+X = MS([wildcard(exclude=response)]).fit_transform(BrainCancer)
 y.iloc[:5]
 ```
 
 ```{code-cell} ipython3
 sm.PHReg(y['time'], X, status=y['status']).fit().params
 ```
 
-## Binomial responses
+### Binomial responses
 
 `Binomial('US')` codes a binary column as 0 and 1, the second level, as in
 R, being a success (or the one given by `level=`):
 
 ```{code-cell} ipython3
-spec = MS(['Price', 'Income'], response=Binomial('US'), intercept=True).fit(Carseats)
-y = spec.build_response(Carseats)
+y = Binomial('US').fit_transform(Carseats)
 pd.crosstab(y, Carseats['US'])
 ```
 
@@ -140,22 +148,41 @@ rng = np.random.default_rng(0)
 D = pd.DataFrame({'x': rng.standard_normal(50), 'size': rng.integers(5, 20, 50)})
 D['cases'] = rng.binomial(D['size'], 1 / (1 + np.exp(-D['x'])))
 
-counts = MS(['x'], response=Binomial('cases', trials='size')).fit(D).build_response(D)
-prop = MS(['x'], response=Binomial('cases', trials='size', output='proportion')).fit(D).build_response(D)
+counts = Binomial('cases', trials='size').fit_transform(D)
+prop = Binomial('cases', trials='size', output='proportion').fit_transform(D)
 X = MS(['x'], intercept=True).fit_transform(D)
 family = sm.families.Binomial()
 (sm.GLM(counts, X, family=family).fit().params,
  sm.GLM(prop['cases/size'], X, family=family, var_weights=prop['size']).fit().params)
 ```
 
-## Multinomial responses
+### Multinomial responses
 
 `Multinomial('ShelveLoc')` gives the class of each row as a categorical
 whose categories are the levels seen at fit, and with `output='onehot'`
 an indicator column for each level; `Multinomial([...])` takes a column
 of counts for each class.
 
 ```{code-cell} ipython3
-spec = MS(['Price', 'Sales'], response=Multinomial('ShelveLoc', output='onehot')).fit(Carseats)
-spec.build_response(Carseats).iloc[:5]
+Multinomial('ShelveLoc', output='onehot').fit_transform(Carseats).iloc[:5]
+```
+
+## Missing values
+
+As scikit-learn's transformers, a `ModelSpec` keeps every row: a row
+of the design is missing where a variable it uses is missing (a
+`poly` or spline is fit on the observed values), and so is the
+response. A transformer cannot drop rows of `X` without those of `y`,
+so the rows to keep, those of R's `na.omit`, are chosen outside the
+`ModelSpec` (or the missing values are imputed, e.g. with
+scikit-learn's `SimpleImputer`, for an estimator that does not take
+them).
+
+```{code-cell} ipython3
+Hitters = load_data('Hitters')
+X = MS([wildcard(exclude='Salary')]).fit_transform(Hitters)
+y = Hitters['Salary']
+keep = X.notna().all(axis=1) & y.notna()
+X, y = X[keep], y[keep]
+X.shape
 ```

modelmatrix/_native.py

@@ -0,0 +1,89 @@
+"""
+Dataframes of other backends
+============================
+
+A dataframe of any eager backend supported by narwhals, e.g. polars or
+pyarrow, is converted to pandas, and the results are converted back.
+"""
+
+import numpy as np, pandas as pd
+import narwhals as nw
+from scipy import sparse as sp_sparse
+
+
+def _series_to_native(S, backend):
+    # a categorical is an Enum, if the backend has one, of its
+    # categories as strings
+    values, dtype = S.to_numpy(), None
+    if isinstance(S.dtype, pd.CategoricalDtype):
+        values = np.asarray([None if pd.isna(v) else str(v) for v in S], dtype=object)
+        dtype = nw.Enum([str(c) for c in S.cat.categories])
+    try:
+        return nw.new_series(str(S.name), values, dtype=dtype, backend=backend).to_native()
+    except NotImplementedError:   # e.g. no Enum in pyarrow
+        return nw.new_series(str(S.name), values, backend=backend).to_native()
+
+def _from_native(X, columns=None):
+    """
+    Convert a dataframe (or series) of any eager backend supported by
+    narwhals, e.g. polars or pyarrow, to pandas. Returns the converted
+    `X` and a function converting a `pd.DataFrame` model matrix back to
+    that backend, or None if `X` was not converted.
+
+    If `columns` is given, only those columns of a dataframe are
+    converted (the copy is of the columns a model matrix uses).
+
+    Conversion through pandas keeps categoricals, including the order
+    of a polars Enum or an ordered pyarrow dictionary (as an ordered
+    `pd.Categorical`). The levels of an unordered categorical, e.g. a
+    polars Categorical, are stored in order of appearance; they are
+    sorted, as pandas does when it infers categories, so that the
+    model matrix does not depend on the order of the rows.
+    """
+    if (isinstance(X, (pd.DataFrame, pd.Series, np.ndarray)) or
+        sp_sparse.issparse(X)):
+        return X, None
+    try:
+        X_nw = nw.from_native(X, allow_series=True)
+    except TypeError:   # not a dataframe, e.g. a list: left to _check_X
+        return X, None
+    if isinstance(X_nw, nw.LazyFrame):
+        raise TypeError('Lazy frames are not supported by ModelSpec; '
+                        'collect them first, e.g. X.collect()')
+    backend = nw.get_native_namespace(X_nw)
+    if columns is not None and isinstance(X_nw, nw.DataFrame):
+        keep = [c for c in X_nw.columns if c in columns]
+        X_nw = X_nw.select(keep or X_nw.columns[:1])
+
+    def to_native(D):
+        if isinstance(D, pd.Series):   # e.g. a response
+            return _series_to_native(D, backend)
+        duplicated = list(D.columns[D.columns.duplicated()])
+        if duplicated:
+            raise ValueError(f'the model matrix has duplicated column names {duplicated}, '
+                             'which this dataframe backend cannot hold')
+        return nw.from_dict({str(c): D[c].to_numpy() for c in D.columns},
+                            backend=backend).to_native()
+
+    X = X_nw.to_pandas()
+    if isinstance(X_nw, nw.Series):
+        dtypes = {X.name: X_nw.dtype}
+    else:
+        dtypes = dict(X_nw.schema)
+    for name, dtype in dtypes.items():
+        col = X if isinstance(X, pd.Series) else X[name]
+        if not isinstance(col.dtype, pd.CategoricalDtype):
+            continue
+        # some versions of polars convert an Enum to an unordered
+        # pd.Categorical: its categories are in the declared order
+        if dtype == nw.Enum:
+            col = col.cat.as_ordered()
+        elif not col.dtype.ordered:
+            col = col.cat.reorder_categories(sorted(col.dtype.categories))
+        else:
+            continue
+        if isinstance(X, pd.Series):
+            X = col
+        else:
+            X[name] = col
+    return X, to_native