diff --git a/CHANGELOG.md b/CHANGELOG.md
index 963bf19..c038ef1 100644
--- a/CHANGELOG.md
+++ b/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)
 
diff --git a/docs/api/responses.rst b/docs/api/responses.rst
index 735ac71..5a42aff 100644
--- a/docs/api/responses.rst
+++ b/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
diff --git a/docs/index.md b/docs/index.md
index 65c27c1..967b102 100644
--- a/docs/index.md
+++ b/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/)
diff --git a/docs/responses.md b/docs/responses.md
index b742304..5331398 100644
--- a/docs/responses.md
+++ b/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,8 +118,9 @@ 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]
 ```
 
@@ -119,14 +128,13 @@ y.iloc[:5]
 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,15 +148,15 @@ 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'`
@@ -156,6 +164,25 @@ 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
 ```
diff --git a/modelmatrix/_native.py b/modelmatrix/_native.py
new file mode 100644
index 0000000..0d1231e
--- /dev/null
+++ b/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
diff --git a/modelmatrix/model_spec.py b/modelmatrix/model_spec.py
index 8dfe008..bf16c3b 100644
--- a/modelmatrix/model_spec.py
+++ b/modelmatrix/model_spec.py
@@ -31,7 +31,8 @@
                       _categorical_from_df,
                       _check_categories)
 
-from .response import _as_response
+from ._native import _from_native
+from .response import Response
 from .transforms import (Poly,
                          BSpline,
                          NaturalSpline,
@@ -306,16 +307,17 @@ def transform(self, X):
         value[missing] = np.nan
         return value
 
-#### the wildcard: all the columns of X but the response, as R's `.`
+#### the wildcard: all the columns of X, as R's `.`
 
 @dataclass(frozen=True, repr=False)
 class Wildcard:
 
     """
-    A term standing for every column of `X` but those of the response,
-    as `.` in R's formulae: `ModelSpec([wildcard], response='y')` is
-    `y ~ .`. It is expanded when the `ModelSpec` is fit, into the terms
-    of `terms_`, and can also be
+    A term standing for every column of `X` (but those excluded), as
+    `.` in R's formulae: with the response `y` a column of `X`,
+    `ModelSpec([wildcard(exclude='y')])` is the design of R's `y ~ .`.
+    It is expanded when the `ModelSpec` is fit, into the terms of
+    `expanded_terms_`, and can also be
 
     - a part of an interaction: `(wildcard, 'g')` is R's `.:g`;
     - a variable of a feature: `pca(wildcard, n_components=3)`.
@@ -330,8 +332,10 @@ class Wildcard:
         Also the interactions of up to `pow` columns: `wildcard(pow=2)`
         is R's `.^2`. Only as a term.
 
-    exclude : column identifier, or sequence of them (optional)
-        Columns left out: `wildcard(exclude='x')` is R's `. - x`.
+    exclude : column identifier, Response, or sequence of them (optional)
+        Columns left out: `wildcard(exclude='x')` is R's `. - x`. A
+        :py:class:`Response`, e.g. `Surv('time', 'status')`, leaves out
+        its columns, so that the response is not a term of its own design.
     """
 
     pow: int = 1
@@ -342,9 +346,13 @@ def __post_init__(self):
             or self.pow < 1):
             raise ValueError(f'pow should be a positive integer, got {self.pow!r}')
         exclude = self.exclude
-        if isinstance(exclude, (str, int, np.integer)) or not np.iterable(exclude):
+        if (isinstance(exclude, (str, int, np.integer, Response)) or
+            not np.iterable(exclude)):
             exclude = (exclude,)
-        object.__setattr__(self, 'exclude', tuple(exclude))
+        columns = []   # a Response by its columns
+        for e in exclude:
+            columns.extend(e.columns if isinstance(e, Response) else (e,))
+        object.__setattr__(self, 'exclude', tuple(columns))
 
     def __call__(self, pow=1, exclude=()):
         return Wildcard(pow=pow, exclude=exclude)
@@ -498,69 +506,6 @@ def _native_columns(X):
         return None
     return list(X_nw.columns)
 
-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):
-        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
-
 class ModelSpec(TransformerMixin, BaseEstimator):
 
     '''Parameters
@@ -571,8 +516,8 @@ class ModelSpec(TransformerMixin, BaseEstimator):
         Sequence of sets whose
         elements are columns of *X* when fit.
         For :py:class:`pd.DataFrame` these can be column
-        names. :py:data:`wildcard` stands for every column of *X*
-        but the response, as R's `.` (see :py:class:`Wildcard`).
+        names. :py:data:`wildcard` stands for every column of *X*,
+        as R's `.` (see :py:class:`Wildcard`).
 
     intercept : bool, default=False
 
@@ -613,14 +558,6 @@ class ModelSpec(TransformerMixin, BaseEstimator):
         For an encoder of a single term, see :py:func:`contrast`
         and :py:func:`ordinal`.
 
-    response : column identifier, list of them, or Response (optional)
-
-        The response of the model, built by :py:meth:`build_response`:
-        a column of *X*, a list of columns, or a :py:class:`Surv`,
-        :py:class:`Binomial` or :py:class:`Multinomial`. Its columns are
-        not in :py:data:`wildcard`, and cannot be in a term.
-        `ModelSpec([wildcard], response='y')` is R's `y ~ .`.
-
     sparse_output : bool (optional)
 
         If True, the model matrix is a `scipy.sparse.csr_matrix`.
@@ -637,7 +574,6 @@ def __init__(self,
                  intercept=False,
                  categorical_features=None,
                  categorical_encoders=None,
-                 response=None,
                  sparse_output=False
                  ):
 
@@ -645,7 +581,6 @@ def __init__(self,
         self.terms = terms
         self.categorical_features = categorical_features
         self.categorical_encoders = categorical_encoders
-        self.response = response
         self.sparse_output = sparse_output
 
     def _check_X(self, X, reset):
@@ -775,7 +710,7 @@ def fit(self, X, y=None):
             Not used to fit, but if both *X* and *y* are pandas objects,
             their indices must agree. This parameter exists for
             compatibility with :py:class:`sklearn.pipeline.Pipeline`;
-            see :py:meth:`build_response` for a response from *X*.
+            see :py:class:`Response` for a response from *X*.
 
         """
         
@@ -786,14 +721,8 @@ def fit(self, X, y=None):
             not X.index.equals(y.index)):
             raise ValueError('the indices of X and y do not agree: '
                              'y should be the response of the rows of X')
-        # a copy of the response is fit, so that `self.response` is not modified
-        self.response_ = deepcopy(_as_response(self.response))
-        response_columns = () if self.response_ is None else self.response_.columns
         native_columns = _native_columns(X)
-        columns = _term_columns(self.terms)
-        if columns is not None:
-            columns |= set(response_columns)
-        X, _ = _from_native(X, columns)
+        X, _ = _from_native(X, _term_columns(self.terms))
         X = self._check_X(X, reset=True)
         if native_columns is not None:   # describe all of X, not only the columns used
             self.n_features_in_ = len(native_columns)
@@ -829,19 +758,8 @@ def fit(self, X, y=None):
                                         bool)
             self.columns_ = np.arange(X.shape[1])
 
-        if self.response_ is not None:
-            self.response_.fit(X)
-
-        # the terms, with each wildcard expanded: the columns of X but
-        # those of the response
-        self.expanded_terms_ = _expand_terms(list(self.terms),
-                                             [c.item() if isinstance(c, np.integer) else c
-                                              for c in self.columns_
-                                              if c not in response_columns])
-        used = _term_columns(self.expanded_terms_) & set(response_columns)
-        if used:
-            raise ValueError(f'the response columns {sorted(used, key=str)} are used '
-                             'in the terms: a response cannot be a term')
+        # the terms, with each wildcard expanded
+        self.expanded_terms_ = _expand_terms(list(self.terms), self._wildcard_columns())
 
         self.features_ = {}
         self.encoders_ = {}
@@ -953,6 +871,11 @@ def fit(self, X, y=None):
     
         return self
     
+    def _wildcard_columns(self):
+        # the columns a wildcard stands for: those of X
+        return [c.item() if isinstance(c, np.integer) else c
+                for c in self.columns_]
+
     def transform(self, X, y=None):
         """
         Build design on X after fitting.
@@ -1009,29 +932,6 @@ def get_feature_names_out(self, input_features=None):
 
     # ModelSpec specific methods
 
-    def build_response(self, X):
-        """
-        Build the response on X after fitting.
-
-        Parameters
-        ----------
-        X : array-like
-            X holding the columns of the response.
-
-        Returns
-        -------
-        y : pd.Series, pd.DataFrame or np.ndarray
-            The response, e.g. for `fit(X, y)` of an estimator: a
-            :py:class:`pd.Series` or :py:class:`pd.DataFrame` for
-            a dataframe *X*, an :py:class:`np.ndarray` for an array.
-        """
-        check_is_fitted(self)
-        if self.response_ is None:
-            raise ValueError('this ModelSpec has no response: give one with ModelSpec(..., response=...)')
-        X, _ = _from_native(X, set(self.response_.columns))
-        X = self._check_X(X, reset=False)
-        return self.response_.transform(X)
-
     @property
     def names(self, help='Name for each term in model specification.'):
         names = []
@@ -1053,7 +953,8 @@ def build_submodel(self,
 
         terms : [Feature]
             Sequence of features. Elements of `self.terms`
-            are replaced by their fitted versions in `self.terms_`.
+            are replaced by their fitted versions in `self.terms_`,
+            and a `wildcard` is expanded as at fit.
 
         Returns
         -------
@@ -1062,6 +963,7 @@ def build_submodel(self,
         """
 
         check_is_fitted(self)
+        terms = _expand_terms(list(terms), self._wildcard_columns())
         X, to_native = _from_native(X, _term_columns(terms))
         X = self._check_X(X, reset=False)
         terms = [self.features_.get(t, t) for t in terms]
@@ -1523,8 +1425,7 @@ def derived_feature(variables, encoder=None, name=None, use_transform=True):
         Variables to apply transform to. Could be
         column identifiers or variables: all columns
         will be stacked before encoding. A :py:data:`wildcard`
-        (alone, or in the sequence) stands for the columns of *X*
-        but the response.
+        (alone, or in the sequence) stands for the columns of *X*.
 
     name : str (optional)
         Defaults to `str(encoder)`.
@@ -1788,8 +1689,8 @@ def pca(variables, name=None, scale=False, **pca_args):
 
     variables : [column identifier, Column, Feature or Wildcard]
         Sequence whose columns will be encoded by PCA, or
-        :py:data:`wildcard` for the columns of *X* but the response:
-        `pca(wildcard, n_components=3)`.
+        :py:data:`wildcard` for the columns of *X*:
+        `pca(wildcard(exclude='y'), n_components=3)`.
 
     name : str (optional)
         Defaults to one derived from the variables.
diff --git a/modelmatrix/response.py b/modelmatrix/response.py
index 1a6fd84..29e143e 100644
--- a/modelmatrix/response.py
+++ b/modelmatrix/response.py
@@ -2,10 +2,8 @@
 Responses
 =========
 
-The response of a model, built from the columns of `X` as a
-`ModelSpec` builds its design: `ModelSpec(terms, response=...)`,
-then `spec.build_response(X)`. A response is a column (or a list of
-columns) of `X`, or one of
+The response of a model, built from the columns of a dataframe as a
+`ModelSpec` builds its design from them. A response is one of
 
 - `Surv`: a survival response, as R's `Surv(time, status)` or
   `Surv(start, stop, status)`;
@@ -14,12 +12,16 @@
 - `Multinomial`: the class of each observation, or counts of
   each class.
 
-A response is fit on the `X` the `ModelSpec` is fit on, e.g. to
-learn the levels of a categorical, so that it is coded the same way
-on new data.
+A response is fit, e.g. to learn the levels of a categorical, so that
+it is coded the same way on new data: `y = Surv('time',
+'status').fit_transform(D)`. To leave its columns out of a design,
+`ModelSpec([wildcard(exclude=response)])` is the design of R's `y ~ .`.
 """
 
 import numpy as np, pandas as pd
+from sklearn.exceptions import NotFittedError
+
+from ._native import _from_native
 
 
 class Response(object):
@@ -27,16 +29,66 @@ class Response(object):
     """
     Base class of the responses: the columns of `X` used, `fit`
     and `transform`.
+
+    `X` is a :py:class:`pd.DataFrame`, a dataframe of another backend
+    supported by narwhals (e.g. polars), whose response is then in that
+    backend, or an :py:class:`np.ndarray`, whose columns are integers.
     """
 
     columns = ()
 
-    def fit(self, X):
-        for col in self.columns:
-            _values(X, col)   # raises if the column is not in X
+    def fit(self, X, y=None):
+        """
+        Fit the response on `X`, e.g. its levels.
+
+        Parameters
+        ----------
+        X : dataframe or np.ndarray
+            Data holding the columns of the response.
+
+        y : None
+            Ignored.
+        """
+        X, _ = _from_native(X, set(self.columns))
+        self._fit(X)
+        self.columns_ = self.columns
         return self
 
     def transform(self, X):
+        """
+        Build the response from `X`, after fitting.
+
+        Parameters
+        ----------
+        X : dataframe or np.ndarray
+            Data holding the columns of the response.
+
+        Returns
+        -------
+        y : series, dataframe or np.ndarray
+            The response, e.g. for `fit(X, y)` of an estimator: of the
+            backend of *X* (e.g. a :py:class:`pd.Series` for a
+            :py:class:`pd.DataFrame`), an :py:class:`np.ndarray` for
+            an array.
+        """
+        if not hasattr(self, 'columns_'):
+            raise NotFittedError(f'this {type(self).__name__} is not fitted yet: '
+                                 'call fit (or fit_transform) first')
+        X, to_native = _from_native(X, set(self.columns))
+        y = self._transform(X)
+        return y if to_native is None else to_native(y)
+
+    def fit_transform(self, X, y=None):
+        """
+        Fit the response on `X` and build it.
+        """
+        return self.fit(X).transform(X)
+
+    def _fit(self, X):
+        for col in self.columns:
+            _values(X, col)   # raises if the column is not in X
+
+    def _transform(self, X):
         raise NotImplementedError
 
     def __repr__(self):
@@ -72,39 +124,6 @@ def _output(D, X):
     return D.to_numpy()
 
 
-class _Columns(Response):
-
-    # a column, or a list of columns, of X as they are
-
-    def __init__(self, columns, single):
-        self.columns = tuple(columns)
-        self.single = single
-
-    def transform(self, X):
-        if isinstance(X, pd.DataFrame):
-            if self.single:
-                return X[self.columns[0]].copy()
-            return X[list(self.columns)].copy()
-        X = np.asarray(X)
-        if self.single:
-            return X[:, self.columns[0]]
-        return X[:, list(self.columns)]
-
-    def __repr__(self):
-        if self.single:
-            return repr(self.columns[0])
-        return repr(list(self.columns))
-
-
-def _as_response(response):
-    # a Response for the `response` parameter of ModelSpec
-    if response is None or isinstance(response, Response):
-        return response
-    if isinstance(response, list):
-        return _Columns(response, single=False)
-    return _Columns((response,), single=True)
-
-
 class Surv(Response):
 
     """
@@ -155,11 +174,10 @@ def columns(self):
             return (self.stop, self.status)
         return (self.start, self.stop, self.status)
 
-    def fit(self, X):
-        self.transform(X)   # checks the columns and their values
-        return self
+    def _fit(self, X):
+        self._transform(X)   # checks the columns and their values
 
-    def transform(self, X):
+    def _transform(self, X):
         status = _numeric(X, self.status)   # booleans as 0 and 1
         observed = status[status.notna()]
         if not observed.isin([0, 1]).all():
@@ -231,12 +249,12 @@ def columns(self):
             return (self.successes,)
         return (self.successes, self.trials)
 
-    def fit(self, X):
+    def _fit(self, X):
         if self.trials is not None:
             if self.output not in ('counts', 'proportion'):
                 raise ValueError(f"output should be 'counts' or 'proportion', got {self.output!r}")
-            self.transform(X)   # checks the values
-            return self
+            self._transform(X)   # checks the values
+            return
         values = _values(X, self.successes)
         observed = values[values.notna()]
         if isinstance(values.dtype, pd.CategoricalDtype):
@@ -261,9 +279,8 @@ def fit(self, X):
         else:
             self.level_ = levels[1]
         self.levels_ = levels
-        return self
 
-    def transform(self, X):
+    def _transform(self, X):
         index = X.index if isinstance(X, pd.DataFrame) else None
         if self.trials is None:
             values = _values(X, self.successes)
@@ -339,10 +356,10 @@ def columns(self):
             return tuple(self.classes)
         return (self.classes,)
 
-    def fit(self, X):
+    def _fit(self, X):
         if isinstance(self.classes, list):
-            self.transform(X)   # checks the counts
-            return self
+            self._transform(X)   # checks the counts
+            return
         if self.output not in ('labels', 'onehot'):
             raise ValueError(f"output should be 'labels' or 'onehot', got {self.output!r}")
         values = _values(X, self.classes)
@@ -350,9 +367,8 @@ def fit(self, X):
             self.levels_ = list(values.cat.categories)
         else:
             self.levels_ = sorted(values[values.notna()].unique())
-        return self
 
-    def transform(self, X):
+    def _transform(self, X):
         index = X.index if isinstance(X, pd.DataFrame) else None
         if isinstance(self.classes, list):
             D = pd.DataFrame({c: _numeric(X, c).to_numpy() for c in self.classes},
diff --git a/modelmatrix/transforms.py b/modelmatrix/transforms.py
index 4c8cf60..9742bfd 100644
--- a/modelmatrix/transforms.py
+++ b/modelmatrix/transforms.py
@@ -38,6 +38,24 @@ def _single_column(X):
         raise ValueError('expecting a single column feature')
     return X
 
+def _observed(X):
+    # the observed (not missing) values of a single feature, and where
+    # the missing ones are: as scikit-learn's transformers, a feature is
+    # fit on the observed values and is missing where the variable is
+    missing = np.isnan(X)
+    if missing.all():
+        raise ValueError('all the values of the feature are missing')
+    return X[~missing], missing
+
+def _fill_missing(value, missing):
+    # the rows of `value` for the observed values, rows of NaN
+    # for the missing ones
+    if not missing.any():
+        return value
+    full = np.full((missing.shape[0],) + value.shape[1:], np.nan)
+    full[~missing] = value
+    return full
+
 class Poly(TransformerMixin, BaseEstimator):
 
     '''
@@ -85,7 +103,7 @@ def fit(self,
         """
         
         X_orig = X
-        X = _single_column(X)
+        X, _ = _observed(_single_column(X))
         n = X.shape[0]
 
         self.mean_ = X.mean()
@@ -134,8 +152,7 @@ def transform(self, X):
         check_is_fitted(self)
 
         X_orig = X
-        
-        X = _single_column(X)
+        X, missing = _observed(_single_column(X))
         n = X.shape[0]
 
         if not self.raw:
@@ -153,6 +170,7 @@ def transform(self, X):
 
         if not self.intercept:
             powX = powX[:,1:]
+        powX = _fill_missing(powX, missing)
 
         if isinstance(X_orig, (pd.Series, pd.DataFrame)):
             df = pd.DataFrame(powX,
@@ -395,7 +413,7 @@ def fit(self,
         """
         
         X_orig = X
-        X = _single_column(X)
+        X, _ = _observed(_single_column(X))
         n = X.shape[0]
 
         order = self.degree + 1
@@ -472,7 +490,7 @@ def transform(self, X):
         check_is_fitted(self)
 
         X_orig = X
-        X = _single_column(X)
+        X, missing = _observed(_single_column(X))
         n = X.shape[0]
 
         value = _splevf(X,
@@ -481,6 +499,7 @@ def transform(self, X):
                         ext=self.ext)
         if not self.intercept:
             value = value[:,1:]
+        value = _fill_missing(value, missing)
         columns_ = self.columns_
 
         if isinstance(X_orig, (pd.Series, pd.DataFrame)):
@@ -559,7 +578,7 @@ def fit(self,
         """
         
         X_orig = X
-        X = _single_column(X)
+        X, _ = _observed(_single_column(X))
         n = X.shape[0]
 
         order = 4
@@ -651,7 +670,7 @@ def transform(self, X):
         lower_bound, upper_bound = self.boundary_knots_
 
         X_orig = X
-        X = _single_column(X)
+        X, missing = _observed(_single_column(X))
         n = X.shape[0]
 
         value = np.empty((X.shape[0], self.nbasis_))
@@ -683,7 +702,7 @@ def transform(self, X):
 
         # map onto the natural splines
 
-        value = value.dot(self.Qproj_)
+        value = _fill_missing(value.dot(self.Qproj_), missing)
         if isinstance(X_orig, (pd.Series, pd.DataFrame)):
             df = pd.DataFrame(value,
                               columns=self.columns_,
diff --git a/tests/test_model_matrix.py b/tests/test_model_matrix.py
index 21990ad..4188a28 100644
--- a/tests/test_model_matrix.py
+++ b/tests/test_model_matrix.py
@@ -1,4 +1,5 @@
 import numpy as np, pandas as pd
+import pytest
 from sklearn.base import clone
 
 from modelmatrix.transforms import Poly, NaturalSpline, BSpline, Interaction
@@ -368,3 +369,19 @@ def test_pca():
 
     
     
+
+
+@pytest.mark.parametrize('feature', [poly('x', 3), ns('x', df=4), bs('x', df=5)],
+                         ids=['poly', 'ns', 'bs'])
+def test_missing_values_pass_through(feature):
+    # as scikit-learn's transformers: fit on the observed values,
+    # missing in the rows where the variable is missing
+    rng = np.random.default_rng(0)
+    D = pd.DataFrame({'x': rng.standard_normal(50)})
+    D.loc[[3, 17], 'x'] = np.nan
+    X = ModelSpec([feature]).fit(D).transform(D)
+    assert list(np.where(X.isna().any(axis=1))[0]) == [3, 17]
+    assert X.drop(index=[3, 17]).notna().all().all()
+    observed = D.dropna()
+    np.testing.assert_allclose(X.drop(index=[3, 17]),
+                               ModelSpec([feature]).fit(observed).transform(observed))
diff --git a/tests/test_model_spec_hard.py b/tests/test_model_spec_hard.py
index 479ca61..a6853b3 100644
--- a/tests/test_model_spec_hard.py
+++ b/tests/test_model_spec_hard.py
@@ -384,7 +384,6 @@ def test_default_params():
     assert MS().get_params() == {'terms': (), 'intercept': False,
                                  'categorical_features': None,
                                  'categorical_encoders': None,
-                                 'response': None,
                                  'sparse_output': False}
 
 def test_unfitted_raises_notfitted(D):
diff --git a/tests/test_response.py b/tests/test_response.py
index 28f9c78..0ceae5d 100644
--- a/tests/test_response.py
+++ b/tests/test_response.py
@@ -1,9 +1,8 @@
-# Responses: ModelSpec(..., response=...) and build_response
+# Responses: Surv, Binomial and Multinomial, fit and built from a dataframe
 
 import numpy as np, pandas as pd
 import pytest
 import statsmodels.api as sm
-from sklearn.base import clone
 from sklearn.exceptions import NotFittedError
 from sklearn.linear_model import LinearRegression
 from sklearn.model_selection import cross_val_score
@@ -35,35 +34,25 @@ def S():
                                                 categories=['lo', 'mid', 'hi'])})
 
 
-# columns
-
-def test_column(Carseats):
-    spec = MS(['Price'], response='Sales').fit(Carseats)
-    y = spec.build_response(Carseats)
-    pd.testing.assert_series_equal(y, Carseats['Sales'])
-    Y = MS(['Price'], response=['Sales', 'Income']).fit(Carseats).build_response(Carseats)
-    pd.testing.assert_frame_equal(Y, Carseats[['Sales', 'Income']])
-
-
-def test_no_response_or_not_fit(Carseats):
+def test_not_fit_or_missing_column(S):
     with pytest.raises(NotFittedError):
-        MS(['Price'], response='Sales').build_response(Carseats)
-    with pytest.raises(ValueError, match='has no response'):
-        MS(['Price']).fit(Carseats).build_response(Carseats)
+        Surv('time', 'status').transform(S)
     with pytest.raises(ValueError, match="'nope' is not in X"):
-        MS(['Price'], response='nope').fit(Carseats)
-
-
-def test_pipeline_with_response(Carseats):
-    spec = MS([wildcard], response='Sales')
-    y = clone(spec).fit(Carseats).build_response(Carseats)
-    scores = cross_val_score(make_pipeline(spec, LinearRegression()), Carseats, y,
-                             cv=3, error_score='raise')
+        Binomial('nope').fit(S)
+
+
+def test_with_a_design(Carseats):
+    # R's glm(US ~ ., family=binomial, data=Carseats)
+    response = Binomial('US')
+    y = response.fit_transform(Carseats)
+    spec = MS([wildcard(exclude=response)], intercept=True)
+    X = spec.fit_transform(Carseats)
+    assert 'US[Yes]' not in X.columns
+    fit = sm.GLM(y, X, family=sm.families.Binomial()).fit()
+    assert np.all(np.isfinite(fit.params))
+    scores = cross_val_score(make_pipeline(MS([wildcard(exclude='Sales')]), LinearRegression()),
+                             Carseats, Carseats['Sales'], cv=3, error_score='raise')
     assert np.all(np.isfinite(scores))
-    # the response object is not modified by fit
-    surv = Surv('a', 'b')
-    MS(['x'], response=surv)
-    assert not hasattr(surv, 'levels_')
 
 
 def test_index_of_y(Carseats):
@@ -78,27 +67,23 @@ def test_index_of_y(Carseats):
 # Surv
 
 def test_surv(S):
-    y = MS(['x'], response=Surv('time', 'status')).fit(S).build_response(S)
+    y = Surv('time', 'status').fit_transform(S)
     pd.testing.assert_frame_equal(y, S[['time', 'status']].astype(float))
-    y3 = MS(['x'], response=Surv('start', 'time', 'status')).fit(S).build_response(S)
+    y3 = Surv('start', 'time', 'status').fit_transform(S)
     assert list(y3.columns) == ['start', 'time', 'status']
-    y_kw = MS(['x'], response=Surv(stop='time', status='status', start='start')).fit(S).build_response(S)
+    y_kw = Surv(stop='time', status='status', start='start').fit_transform(S)
     pd.testing.assert_frame_equal(y3, y_kw)
     assert repr(Surv('start', 'time', 'status')) == "Surv('start', 'time', 'status')"
     # a boolean status
-    yb = MS(['x'], response=Surv('time', 'flag')).fit(S).build_response(S)
+    yb = Surv('time', 'flag').fit_transform(S)
     np.testing.assert_array_equal(yb['flag'], S['flag'].astype(float))
 
 
-def test_surv_wildcard_excludes_its_columns(S):
-    spec = MS([wildcard(exclude=['cases', 'size', 'flag', 'kind'])],
-              response=Surv('start', 'time', 'status')).fit(S)
-    assert spec.expanded_terms_ == ['x']
-
-
 def test_surv_matches_statsmodels_PHReg(S):
-    spec = MS(['x'], response=Surv('time', 'status')).fit(S)
-    X, y = spec.transform(S), spec.build_response(S)
+    response = Surv('time', 'status')
+    y = response.fit_transform(S)
+    X = MS([wildcard(exclude=[response, 'start', 'cases', 'size', 'flag', 'kind'])]).fit_transform(S)
+    assert list(X.columns) == ['x']
     fit = sm.PHReg(y['time'], X, status=y['status']).fit()
     ref = sm.PHReg(S['time'], S[['x']], status=S['status']).fit()
     np.testing.assert_allclose(fit.params, ref.params)
@@ -110,87 +95,102 @@ def test_surv_errors(S):
     with pytest.raises(TypeError, match='takes'):
         Surv('time')
     with pytest.raises(ValueError, match='should be 0 or 1'):
-        MS(['x'], response=Surv('time', 'cases')).fit(S)
+        Surv('time', 'cases').fit(S)
     bad = S.assign(start=S['time'] + 1)
     with pytest.raises(ValueError, match='less than'):
-        MS(['x'], response=Surv('start', 'time', 'status')).fit(bad)
+        Surv('start', 'time', 'status').fit(bad)
 
 
 # Binomial
 
 def test_binomial_binary(Carseats):
-    y = MS(['Price'], response=Binomial('US')).fit(Carseats).build_response(Carseats)
+    y = Binomial('US').fit_transform(Carseats)
     # the second level, as in R: No < Yes
     np.testing.assert_array_equal(y, (Carseats['US'] == 'Yes').astype(float))
-    y_no = MS(['Price'], response=Binomial('US', level='No')).fit(Carseats).build_response(Carseats)
+    y_no = Binomial('US', level='No').fit_transform(Carseats)
     np.testing.assert_array_equal(y_no, 1 - y)
 
 
 def test_binomial_binary_values(S):
-    y = MS(['x'], response=Binomial('status')).fit(S).build_response(S)
-    np.testing.assert_array_equal(y, S['status'])
-    y = MS(['x'], response=Binomial('flag')).fit(S).build_response(S)
-    np.testing.assert_array_equal(y, S['flag'].astype(float))
+    np.testing.assert_array_equal(Binomial('status').fit_transform(S), S['status'])
+    np.testing.assert_array_equal(Binomial('flag').fit_transform(S), S['flag'].astype(float))
     with pytest.raises(ValueError, match='should be 0 or 1'):
-        MS(['x'], response=Binomial('cases')).fit(S)
+        Binomial('cases').fit(S)
     with pytest.raises(ValueError, match='should have two levels'):
-        MS(['x'], response=Binomial('kind')).fit(S)
-    y = MS(['x'], response=Binomial('kind', level='hi')).fit(S).build_response(S)
+        Binomial('kind').fit(S)
+    y = Binomial('kind', level='hi').fit_transform(S)
     np.testing.assert_array_equal(y, (S['kind'] == 'hi').astype(float))
 
 
 def test_binomial_missing_and_unseen():
     D = pd.DataFrame({'x': [1., 2, 3, 4], 'y': ['a', 'b', None, 'a']})
-    spec = MS(['x'], response=Binomial('y')).fit(D)
-    y = spec.build_response(D)
-    np.testing.assert_array_equal(y, [0, 1, np.nan, 0])
+    response = Binomial('y').fit(D)
+    np.testing.assert_array_equal(response.transform(D), [0, 1, np.nan, 0])
     with pytest.raises(ValueError, match='not seen at fit'):
-        spec.build_response(D.assign(y=['a', 'c', 'b', 'a']))
+        response.transform(D.assign(y=['a', 'c', 'b', 'a']))
 
 
 def test_binomial_counts_and_proportion_agree(S):
     X = MS(['x'], intercept=True).fit_transform(S)
-    counts = MS(['x'], response=Binomial('cases', trials='size')).fit(S).build_response(S)
+    counts = Binomial('cases', trials='size').fit_transform(S)
     assert list(counts.columns) == ['cases', 'size-cases']
     np.testing.assert_array_equal(counts.sum(axis=1), S['size'])
-    prop = MS(['x'], response=Binomial('cases', trials='size', output='proportion')).fit(S).build_response(S)
+    prop = Binomial('cases', trials='size', output='proportion').fit_transform(S)
     assert list(prop.columns) == ['cases/size', 'size']
     family = sm.families.Binomial()
     fit_counts = sm.GLM(counts, X, family=family).fit()
     fit_prop = sm.GLM(prop['cases/size'], X, family=family, var_weights=prop['size']).fit()
     np.testing.assert_allclose(fit_counts.params, fit_prop.params)
     with pytest.raises(ValueError, match='between 0'):
-        MS(['x'], response=Binomial('size', trials='cases')).fit(S)
+        Binomial('size', trials='cases').fit(S)
     with pytest.raises(ValueError, match="output should be"):
-        MS(['x'], response=Binomial('cases', trials='size', output='odds')).fit(S)
+        Binomial('cases', trials='size', output='odds').fit(S)
 
 
 # Multinomial
 
 def test_multinomial_labels(S):
-    spec = MS(['x'], response=Multinomial('kind')).fit(S)
-    y = spec.build_response(S)
+    response = Multinomial('kind').fit(S)
+    y = response.transform(S)
     assert list(y.cat.categories) == ['lo', 'mid', 'hi']
     np.testing.assert_array_equal(y.astype(str), S['kind'].astype(str))
     with pytest.raises(ValueError, match='not seen at fit'):
-        spec.build_response(S.assign(kind=S['kind'].cat.add_categories('new').fillna('new')
-                                     .where(S.index != 0, 'new')))
+        response.transform(S.assign(kind=S['kind'].cat.add_categories('new')
+                                    .where(S.index != 0, 'new')))
 
 
 def test_multinomial_onehot_and_counts(S):
-    Y = MS(['x'], response=Multinomial('kind', output='onehot')).fit(S).build_response(S)
+    Y = Multinomial('kind', output='onehot').fit_transform(S)
     assert list(Y.columns) == ['kind[lo]', 'kind[mid]', 'kind[hi]']
     np.testing.assert_array_equal(Y.sum(axis=1), 1)
     np.testing.assert_array_equal(Y['kind[hi]'], (S['kind'] == 'hi').astype(float))
-    C = MS(['x'], response=Multinomial(['cases', 'size'])).fit(S).build_response(S)
+    C = Multinomial(['cases', 'size']).fit_transform(S)
     np.testing.assert_array_equal(C, S[['cases', 'size']].astype(float))
     strings = S.assign(kind=S['kind'].astype(str))
-    y = MS(['x'], response=Multinomial('kind')).fit(strings).build_response(strings)
+    y = Multinomial('kind').fit_transform(strings)
     assert list(y.cat.categories) == ['hi', 'lo', 'mid']   # sorted
 
 
 def test_ndarray_output(S):
     A = np.asarray(S[['time', 'status', 'x']])
-    y = MS([2], response=Surv(0, 1)).fit(A).build_response(A)
+    y = Surv(0, 1).fit_transform(A)
     assert isinstance(y, np.ndarray) and y.shape == (40, 2)
     np.testing.assert_allclose(y, A[:, :2])
+
+
+@pytest.mark.parametrize('backend', ['polars', 'pyarrow'])
+def test_native_output(S, backend):
+    lib = pytest.importorskip(backend)
+    nw = pytest.importorskip('narwhals')
+    P = lib.from_pandas(S) if backend == 'polars' else lib.Table.from_pandas(S)
+    Y = Surv('time', 'status').fit_transform(P)
+    assert nw.get_native_namespace(nw.from_native(Y)) is lib
+    np.testing.assert_allclose(nw.from_native(Y).to_numpy(), S[['time', 'status']])
+    y = Binomial('flag').fit_transform(P)
+    assert nw.get_native_namespace(nw.from_native(y, series_only=True)) is lib
+    np.testing.assert_allclose(np.asarray(y), S['flag'])
+    labels = Multinomial('kind').fit_transform(P)
+    assert nw.from_native(labels, series_only=True).to_list() == list(S['kind'])
+    if backend == 'polars':
+        # the levels of polars' (unordered) Categorical are sorted
+        assert labels.dtype == lib.Enum(['hi', 'lo', 'mid'])
diff --git a/tests/test_wildcard.py b/tests/test_wildcard.py
index 2e7c565..014bc7d 100644
--- a/tests/test_wildcard.py
+++ b/tests/test_wildcard.py
@@ -1,4 +1,4 @@
-# The wildcard: all the columns of X but the response, as R's `.`
+# The wildcard: all the columns of X (but those excluded), as R's `.`
 
 from itertools import combinations
 
@@ -9,7 +9,8 @@
 from sklearn.model_selection import GridSearchCV, cross_val_score
 from sklearn.pipeline import make_pipeline
 
-from modelmatrix import ModelSpec as MS, wildcard, Wildcard, pca, poly, derived_feature
+from modelmatrix import (ModelSpec as MS, wildcard, Wildcard, pca, poly, derived_feature,
+                         Surv, Binomial)
 
 from _datasets import load_data
 
@@ -36,85 +37,86 @@ def _same(A, B):
 
 def test_all_but_response(Carseats):
     others = list(Carseats.columns.drop('Sales'))
-    spec = MS([wildcard], response='Sales').fit(Carseats)
+    spec = MS([wildcard(exclude='Sales')]).fit(Carseats)
     _same(spec.transform(Carseats), MS(others).fit_transform(Carseats))
     assert spec.expanded_terms_ == others
     assert spec.names == others
-    assert list(spec.terms) == [wildcard]   # terms are not modified
+    assert list(spec.terms) == [wildcard(exclude='Sales')]   # terms are not modified
 
 
-def test_no_response_all_columns(D):
+def test_all_columns(D):
     _same(MS([wildcard]).fit_transform(D), MS(list(D.columns)).fit_transform(D))
 
 
 def test_explicit_terms_not_repeated(D):
     # R's y ~ x + .
-    spec = MS(['x', wildcard], response='y').fit(D)
+    dot = wildcard(exclude='y')
+    spec = MS(['x', dot]).fit(D)
     assert spec.expanded_terms_ == ['x', 'z', 'g']
-    spec = MS([wildcard, 'x', ('x', 'z')], response='y').fit(D)
+    spec = MS([dot, 'x', ('x', 'z')]).fit(D)
     assert spec.expanded_terms_ == ['z', 'g', 'x', ('x', 'z')]
     # two wildcards: their columns once
-    assert MS([wildcard, wildcard], response='y').fit(D).expanded_terms_ == ['x', 'z', 'g']
+    assert MS([dot, dot]).fit(D).expanded_terms_ == ['x', 'z', 'g']
 
 
 def test_pow(D):
     # R's y ~ .^2
-    spec = MS([wildcard(pow=2)], response='y').fit(D)
+    spec = MS([wildcard(pow=2, exclude='y')]).fit(D)
     cols = ['x', 'z', 'g']
     assert spec.expanded_terms_ == cols + list(combinations(cols, 2))
     _same(spec.transform(D), MS(cols + list(combinations(cols, 2))).fit_transform(D))
-    spec3 = MS([wildcard(pow=3)], response='y').fit(D)
+    spec3 = MS([wildcard(pow=3, exclude='y')]).fit(D)
     assert len(spec3.expanded_terms_) == 3 + 3 + 1
     # an interaction given explicitly, in another order, is not repeated
-    spec = MS([('z', 'x'), wildcard(pow=2)], response='y').fit(D)
+    spec = MS([('z', 'x'), wildcard(pow=2, exclude='y')]).fit(D)
     assert spec.expanded_terms_ == [('z', 'x'), 'x', 'z', 'g', ('x', 'g'), ('z', 'g')]
 
 
 def test_exclude(D):
     # R's y ~ . - z
-    assert MS([wildcard(exclude='z')], response='y').fit(D).expanded_terms_ == ['x', 'g']
-    assert MS([wildcard(exclude=['z', 'g'])], response='y').fit(D).expanded_terms_ == ['x']
+    assert MS([wildcard(exclude='z')]).fit(D).expanded_terms_ == ['y', 'x', 'g']
+    assert MS([wildcard(exclude=['y', 'z', 'g'])]).fit(D).expanded_terms_ == ['x']
     assert repr(wildcard(pow=2, exclude=['z'])) == '.^2 - z'
     assert repr(wildcard) == '.'
 
 
 def test_interaction_with_wildcard(D):
     # R's y ~ .:g
-    spec = MS([(wildcard, 'g')], response='y').fit(D)
+    spec = MS([(wildcard(exclude='y'), 'g')]).fit(D)
     assert spec.expanded_terms_ == [('x', 'g'), ('z', 'g')]
     p = poly('x', 2)
-    spec = MS([(p, wildcard(exclude='x'))], response='y').fit(D)
+    spec = MS([(p, wildcard(exclude=['y', 'x']))]).fit(D)
     assert spec.expanded_terms_ == [(p, 'z'), (p, 'g')]
     with pytest.raises(ValueError, match='cannot be a part of an interaction'):
-        MS([(wildcard(pow=2), 'g')], response='y').fit(D)
+        MS([(wildcard(pow=2), 'g')]).fit(D)
 
 
 def test_pca_of_wildcard(Carseats):
     numeric = [c for c in Carseats.columns.drop('Sales')
                if pd.api.types.is_numeric_dtype(Carseats[c])]
-    pc = pca(wildcard(exclude=['ShelveLoc', 'Urban', 'US']), n_components=3)
-    assert pc.name == 'pca(. - ShelveLoc - Urban - US, n_components=3)'
-    spec = MS([pc], response='Sales').fit(Carseats)
+    pc = pca(wildcard(exclude=['Sales', 'ShelveLoc', 'Urban', 'US']), n_components=3)
+    assert pc.name == 'pca(. - Sales - ShelveLoc - Urban - US, n_components=3)'
+    spec = MS([pc]).fit(Carseats)
     X = spec.transform(Carseats)
     X_explicit = MS([pca(numeric, name='pc', n_components=3)]).fit_transform(Carseats)
     np.testing.assert_allclose(X, X_explicit)
     assert X.shape == (400, 3)
-    # the feature's variables are the columns, but the response
+    # the feature's variables are the columns, but those excluded
     assert spec.terms_[0].variables[:2] == ('CompPrice', 'Income')
     assert 'Sales' not in spec.terms_[0].variables
 
 
 def test_wildcard_in_derived_feature(D):
-    f = derived_feature(wildcard(exclude='g'), name='xz')
-    X = MS([f], response='y').fit_transform(D)
+    f = derived_feature(wildcard(exclude=['y', 'g']), name='xz')
+    X = MS([f]).fit_transform(D)
     np.testing.assert_allclose(X, D[['x', 'z']])
     with pytest.raises(ValueError, match='cannot be a variable'):
         MS([derived_feature([wildcard(pow=2)], name='bad')]).fit(D)
 
 
 def test_tuning_pca_of_wildcard(Carseats):
-    pc = pca(wildcard(exclude=['ShelveLoc', 'Urban', 'US']), n_components=2)
-    pipe = make_pipeline(MS([pc, 'ShelveLoc'], response='Sales'), LinearRegression())
+    pc = pca(wildcard(exclude=['Sales', 'ShelveLoc', 'Urban', 'US']), n_components=2)
+    pipe = make_pipeline(MS([pc, 'ShelveLoc']), LinearRegression())
     grid = GridSearchCV(pipe, {'modelspec__terms__0__n_components': [1, 3]}, cv=3)
     grid.fit(Carseats, Carseats['Sales'])
     best = grid.best_estimator_[0]
@@ -122,34 +124,39 @@ def test_tuning_pca_of_wildcard(Carseats):
 
 
 def test_transform_without_response(D):
-    spec = MS([wildcard, ('x', 'g')], response='y').fit(D)
+    spec = MS([wildcard(exclude='y'), ('x', 'g')]).fit(D)
     _same(spec.transform(D.drop(columns='y')), spec.transform(D))
 
 
 def test_clone_and_cross_validate(D):
-    spec = MS([wildcard], response='y')
+    spec = MS([wildcard(exclude='y')])
     cloned = clone(spec)
-    assert cloned.terms == [wildcard]
+    assert cloned.terms == [wildcard(exclude='y')]
     scores = cross_val_score(make_pipeline(spec, LinearRegression()), D, D['y'],
                              cv=3, error_score='raise')
     assert np.all(np.isfinite(scores))
 
 
-def test_response_in_term_raises(D):
-    with pytest.raises(ValueError, match='a response cannot be a term'):
-        MS(['y', 'x'], response='y').fit(D)
-    with pytest.raises(ValueError, match='a response cannot be a term'):
-        MS([poly('y', 2)], response='y').fit(D)
+def test_exclude_response(D, Carseats):
+    # a Response is excluded by its columns
+    S = D.assign(time=np.arange(1, 31.), status=np.tile([0, 1], 15))
+    surv = Surv('time', 'status')
+    dot = wildcard(exclude=[surv, 'y'])
+    assert dot == wildcard(exclude=['time', 'status', 'y'])
+    assert repr(dot) == '. - time - status - y'
+    assert MS([dot]).fit(S).expanded_terms_ == ['x', 'z', 'g']
+    US = Binomial('US')
+    spec = MS([wildcard(exclude=US)]).fit(Carseats)
+    assert 'US' not in spec.expanded_terms_ and len(spec.expanded_terms_) == 10
 
 
 def test_ndarray(D):
     A = np.asarray(D[['y', 'x', 'z']])
-    spec = MS([wildcard], response=0).fit(A)
+    spec = MS([wildcard(exclude=0)]).fit(A)
     assert spec.expanded_terms_ == [1, 2]
     np.testing.assert_allclose(spec.transform(A), A[:, 1:])
-    np.testing.assert_allclose(spec.build_response(A), A[:, 0])
     # an explicit column is not repeated
-    spec = MS([2, wildcard], response=0).fit(A)
+    spec = MS([2, wildcard(exclude=0)]).fit(A)
     assert spec.expanded_terms_ == [2, 1]
     assert all(type(t) is int for t in spec.expanded_terms_)
 
@@ -157,12 +164,10 @@ def test_ndarray(D):
 def test_polars(D):
     pl = pytest.importorskip('polars')
     P = pl.from_pandas(D)
-    spec = MS([wildcard], response='y').fit(P)
+    spec = MS([wildcard(exclude='y')]).fit(P)
     X = spec.transform(P)
     assert isinstance(X, pl.DataFrame)
     np.testing.assert_allclose(X.to_numpy(), MS(['x', 'z', 'g']).fit_transform(D))
-    y = spec.build_response(P)
-    np.testing.assert_allclose(y, D['y'])
 
 
 def test_wildcard_validation():
@@ -172,3 +177,10 @@ def test_wildcard_validation():
         wildcard(pow=1.5)
     assert wildcard(pow=2) == Wildcard(pow=2)
     assert hash(wildcard(exclude=['x'])) == hash(Wildcard(exclude=('x',)))
+
+
+def test_build_submodel_with_wildcard(D):
+    spec = MS([wildcard(exclude='y')]).fit(D)
+    _same(spec.build_submodel(D, [wildcard(exclude='y')]), spec.transform(D))
+    _same(spec.build_submodel(D, ['x', wildcard(exclude=['y', 'g'])]),
+          spec.build_submodel(D, ['x', 'z']))