diff --git a/CHANGELOG.md b/CHANGELOG.md index 7091c85..963bf19 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,37 @@ ## Unreleased +### 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 + - `Surv(time, status)` or `Surv(start, stop, status)`, as R's: a + `pd.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, + trials=...)` for the successes and failures (`output='counts'`, R's + `cbind`) or the proportion and number of trials (`output='proportion'`); + - `Multinomial('class')` for the class of each row, as a categorical or + 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 + .`. +- `pca()`'s `name` is optional. +- `fit(X, y)` checks that the indices of `X` and `y` agree when both are + pandas objects. + +### Documentation + +- A page on responses and the wildcard. + ## 0.1.0 (2026-10-26) `ModelSpec`, from [ISLP](https://github.com/intro-stat-learning/ISLP) 0.6.0, diff --git a/docs/api/index.rst b/docs/api/index.rst index cc5d4b0..09895c2 100644 --- a/docs/api/index.rst +++ b/docs/api/index.rst @@ -9,4 +9,5 @@ For the guide and examples, see the `modelmatrix documentation <../>`_. ModelSpec terms + responses transforms diff --git a/docs/api/responses.rst b/docs/api/responses.rst new file mode 100644 index 0000000..735ac71 --- /dev/null +++ b/docs/api/responses.rst @@ -0,0 +1,11 @@ +Responses +========= + +.. automodule:: modelmatrix.response + :no-members: + +.. autoclass:: modelmatrix.Surv +.. autoclass:: modelmatrix.Binomial +.. autoclass:: modelmatrix.Multinomial +.. autoclass:: modelmatrix.Response + :members: fit, transform diff --git a/docs/api/terms.rst b/docs/api/terms.rst index a1811ff..5399b2f 100644 --- a/docs/api/terms.rst +++ b/docs/api/terms.rst @@ -13,6 +13,13 @@ Helpers that make terms .. autofunction:: modelmatrix.derived_feature .. autofunction:: modelmatrix.build_columns +The wildcard +------------ + +.. autodata:: modelmatrix.wildcard + :annotation: = Wildcard() +.. autoclass:: modelmatrix.Wildcard + Terms and columns ----------------- diff --git a/docs/index.md b/docs/index.md index 6d58d33..65c27c1 100644 --- a/docs/index.md +++ b/docs/index.md @@ -14,6 +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` - [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/myst.yml b/docs/myst.yml index d2b9026..ac8588d 100644 --- a/docs/myst.yml +++ b/docs/myst.yml @@ -18,6 +18,7 @@ project: - file: poly.md - file: splines.md - file: PCA.md + - file: responses.md - file: anova.md - file: comparison.md - title: API Reference diff --git a/docs/responses.md b/docs/responses.md new file mode 100644 index 0000000..b742304 --- /dev/null +++ b/docs/responses.md @@ -0,0 +1,161 @@ +--- +jupytext: + formats: ipynb,md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 +kernelspec: + display_name: python3 + language: python + name: python3 +--- + +# 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 +`Sales ~ .`. + +```{code-cell} ipython3 +import numpy as np, pandas as pd +import statsmodels.api as sm +from sklearn.linear_model import LinearRegression +from sklearn.model_selection import GridSearchCV +from sklearn.pipeline import make_pipeline + +from ISLP import load_data +from modelmatrix import (ModelSpec as MS, + wildcard, + pca, + poly, + Surv, + Binomial, + Multinomial) + +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) +sm.OLS(y, X).fit().params +``` + +The wildcard is expanded when the `ModelSpec` is fit: `terms` is +unchanged, and `expanded_terms_` holds what it stood for. + +```{code-cell} ipython3 +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`: + +```{code-cell} ipython3 +MS([poly('Price', 2), 'Price', wildcard(exclude=['Population', 'Education'])], + response='Sales').fit(Carseats).expanded_terms_ +``` + +## `y ~ .^2` and `.:x` + +`wildcard(pow=2)` adds every interaction of two columns, as R's `.^2`, +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_ +``` + +```{code-cell} ipython3 +MS(['US', (wildcard, 'US')], response='Sales').fit(small).expanded_terms_ +``` + +## Features of the wildcard + +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) +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()) +grid = GridSearchCV(pipe, {'modelspec__terms__0__n_components': [1, 3, 5, 7]}, cv=5) +grid.fit(Carseats, y) +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. + +## 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 +`pd.DataFrame` with these columns, named as in the data, as glmstar's +`CoxFamily(event_id=..., status_id=..., start_id=...)` or lifelines take it. + +```{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) +y.iloc[:5] +``` + +```{code-cell} ipython3 +sm.PHReg(y['time'], X, status=y['status']).fit().params +``` + +## 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) +pd.crosstab(y, Carseats['US']) +``` + +With the number of trials, `Binomial(successes, trials=...)` gives the +successes and the failures, as R's `cbind(successes, failures)` (which +statsmodels' `GLM` takes too), or with `output='proportion'` the proportion +and the number of trials, as a response and weights: + +```{code-cell} ipython3 +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) +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('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] +``` diff --git a/modelmatrix/__init__.py b/modelmatrix/__init__.py index fbe2597..c45089b 100644 --- a/modelmatrix/__init__.py +++ b/modelmatrix/__init__.py @@ -13,8 +13,11 @@ pca, contrast, ordinal, - build_columns) + build_columns, + wildcard, + Wildcard) from .columns import Column +from .response import Response, Surv, Binomial, Multinomial try: from ._version import version as __version__ @@ -23,4 +26,5 @@ __all__ = ['ModelSpec', 'Feature', 'Column', 'Contrast', 'Ordinal', 'poly', 'ns', 'bs', 'derived_feature', 'pca', 'contrast', - 'ordinal', 'build_columns'] + 'ordinal', 'build_columns', 'wildcard', 'Wildcard', + 'Response', 'Surv', 'Binomial', 'Multinomial'] diff --git a/modelmatrix/model_spec.py b/modelmatrix/model_spec.py index 82e5f6d..8dfe008 100644 --- a/modelmatrix/model_spec.py +++ b/modelmatrix/model_spec.py @@ -7,9 +7,10 @@ """ from collections import namedtuple -from itertools import product +from itertools import product, combinations +from dataclasses import dataclass from typing import NamedTuple, Any -from copy import copy +from copy import copy, deepcopy import numpy as np, pandas as pd import narwhals as nw @@ -30,6 +31,7 @@ _categorical_from_df, _check_categories) +from .response import _as_response from .transforms import (Poly, BSpline, NaturalSpline, @@ -304,13 +306,167 @@ def transform(self, X): value[missing] = np.nan return value +#### the wildcard: all the columns of X but the response, 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 part of an interaction: `(wildcard, 'g')` is R's `.:g`; + - a variable of a feature: `pca(wildcard, n_components=3)`. + + Columns that are terms themselves are not repeated, as in R's + `y ~ x + .`. + + Parameters + ---------- + + pow : int, default=1 + 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`. + """ + + pow: int = 1 + exclude: tuple = () + + def __post_init__(self): + if (not isinstance(self.pow, (int, np.integer)) or isinstance(self.pow, bool) + 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): + exclude = (exclude,) + object.__setattr__(self, 'exclude', tuple(exclude)) + + def __call__(self, pow=1, exclude=()): + return Wildcard(pow=pow, exclude=exclude) + + def __repr__(self): + name = '.' if self.pow == 1 else f'.^{self.pow}' + return ' - '.join([name] + [str(c) for c in self.exclude]) + + @property + def name(self): + return repr(self) + + def _columns(self, columns): + return [c for c in columns if c not in self.exclude] + + +wildcard = Wildcard() + + +def _has_wildcard(term): + if isinstance(term, Wildcard): + return True + if isinstance(term, Feature): + return any(_has_wildcard(v) for v in term.variables) + if isinstance(term, tuple) and not isinstance(term, Column): + return any(_has_wildcard(v) for v in term) + return False + + +def _same_term(a, b): + # an interaction is the same in any order + if (isinstance(a, tuple) and isinstance(b, tuple) and + not isinstance(a, (Feature, Column)) and not isinstance(b, (Feature, Column))): + return len(a) == len(b) and all(any(_same_term(x, y) for y in b) for x in a) + if isinstance(a, tuple) or isinstance(b, tuple): # Features and Columns are tuples + return type(a) is type(b) and a == b + try: # column identifiers, e.g. 1 and np.int64(1) + return bool(a == b) + except (TypeError, ValueError): + return False + + +def _expand_feature(feature, columns): + # a Feature with each wildcard among its variables (recursively) + # replaced by the columns it stands for + variables = [] + for v in feature.variables: + if isinstance(v, Wildcard): + if v.pow != 1: + raise ValueError(f'{v!r} cannot be a variable of the feature {feature.name!r}: ' + 'its interactions are terms, use it as a term') + variables.extend(c for c in v._columns(columns) + if not any(_same_term(c, u) for u in feature.variables)) + elif isinstance(v, Feature) and _has_wildcard(v): + variables.append(_expand_feature(v, columns)) + else: + variables.append(v) + return feature._replace(variables=tuple(variables)) + + +def _expand_term(term, columns): + # the terms a term with a wildcard stands for + if isinstance(term, Wildcard): + cols = term._columns(columns) + terms = list(cols) + for k in range(2, term.pow + 1): + terms.extend(combinations(cols, k)) + return terms + if isinstance(term, Feature): + return [_expand_feature(term, columns)] + # an interaction: one for each column a wildcard among its parts stands for + parts = [] + for part in term: + if isinstance(part, Wildcard): + if part.pow != 1: + raise ValueError(f'{part!r} cannot be a part of an interaction; ' + 'use it as a term') + parts.append(part._columns(columns)) + elif isinstance(part, Feature) and _has_wildcard(part): + parts.append([_expand_feature(part, columns)]) + else: + parts.append([part]) + terms = [] + for combo in product(*parts): + if any(_same_term(a, b) for i, a in enumerate(combo) for b in combo[i+1:]): + continue # a column interacted with itself + terms.append(tuple(combo)) + return terms + + +def _expand_terms(terms, columns): + """ + `terms` with each wildcard expanded, `columns` being the columns + it stands for. A term made by a wildcard is dropped if it is + another term, or was already made. + """ + explicit = [t for t in terms if not _has_wildcard(t)] + expanded = [] + for term in terms: + if not _has_wildcard(term): + expanded.append(term) + continue + for new in _expand_term(term, columns): + if not any(_same_term(new, t) for t in explicit + expanded): + expanded.append(new) + return expanded + + +class _AllColumns(Exception): + pass + + def _term_columns(terms): """ The columns of X that `terms` (column identifiers, Columns, - Features and tuples of these) are built from. + Features and tuples of these) are built from, or None if a + wildcard stands for all of them. """ found = set() def visit(t): + if isinstance(t, Wildcard): + raise _AllColumns if isinstance(t, Feature): for v in t.variables: visit(v) @@ -324,8 +480,11 @@ def visit(t): found.add(t) except TypeError: # not a column identifier: an error later pass - for t in terms: - visit(t) + try: + for t in terms: + visit(t) + except _AllColumns: + return None return found def _native_columns(X): @@ -412,7 +571,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. + names. :py:data:`wildcard` stands for every column of *X* + but the response, as R's `.` (see :py:class:`Wildcard`). intercept : bool, default=False @@ -453,6 +613,14 @@ 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`. @@ -469,6 +637,7 @@ def __init__(self, intercept=False, categorical_features=None, categorical_encoders=None, + response=None, sparse_output=False ): @@ -476,6 +645,7 @@ 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): @@ -601,16 +771,29 @@ def fit(self, X, y=None): converted to `pandas`; `transform` then returns the model matrix in the backend of its input. - y : None - Ignored. This parameter exists only for compatibility with - :py:class:`sklearn.pipeline.Pipeline`. + y : array-like (optional) + 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*. """ if not isinstance(self.sparse_output, (bool, np.bool_)): raise ValueError(f'sparse_output should be True or False, got {self.sparse_output!r}') + if (isinstance(X, (pd.DataFrame, pd.Series)) and + isinstance(y, (pd.DataFrame, pd.Series)) and + 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) - X, _ = _from_native(X, _term_columns(self.terms)) + columns = _term_columns(self.terms) + if columns is not None: + columns |= set(response_columns) + X, _ = _from_native(X, columns) 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) @@ -646,6 +829,20 @@ 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') + self.features_ = {} self.encoders_ = {} @@ -675,7 +872,7 @@ def fit(self, X, y=None): # fit does not build the model matrix X_head = X.iloc[:2] if isinstance(X, pd.DataFrame) else X[:2] - for term in self.terms: + for term in self.expanded_terms_: if isinstance(term, Feature): term_ = _clone_term(term, clone_memo) self.features_[term] = term_ @@ -738,14 +935,14 @@ def fit(self, X, y=None): self.column_names_ = {} self.column_map_ = {} - self.terms_ = [self.features_[t] for t in self.terms] + self.terms_ = [self.features_[t] for t in self.expanded_terms_] idx = 0 if self._include_intercept(): self.column_map_['intercept'] = slice(0, 1) idx += 1 # intercept will be first column - for term, term_ in zip(self.terms, self.terms_): + for term, term_ in zip(self.expanded_terms_, self.terms_): term_df, term_names = build_columns(self.column_info_, X_head, term_, @@ -804,7 +1001,7 @@ def get_feature_names_out(self, input_features=None): """ check_is_fitted(self) names = ['intercept'] if self._include_intercept() else [] - for term in self.terms: + for term in self.expanded_terms_: names.extend(str(n) for n in self.column_names_[term]) if not names: # build_model returns a 'zero' column names = ['zero'] @@ -812,6 +1009,29 @@ 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 = [] @@ -1299,10 +1519,12 @@ def derived_feature(variables, encoder=None, name=None, use_transform=True): Parameters ---------- - variables : [column identifier, Column, Feature] + variables : [column identifier, Column, Feature, Wildcard] Variables to apply transform to. Could be column identifiers or variables: all columns - will be stacked before encoding. + will be stacked before encoding. A :py:data:`wildcard` + (alone, or in the sequence) stands for the columns of *X* + but the response. name : str (optional) Defaults to `str(encoder)`. @@ -1318,6 +1540,8 @@ def derived_feature(variables, encoder=None, name=None, use_transform=True): if name is None: name = str(encoder) + if isinstance(variables, Wildcard): + variables = [variables] var = Feature(tuple([v for v in variables]), name, encoder, @@ -1551,19 +1775,27 @@ def bs(col, intercept=False, name=None, **spline_args): name=name, encoder=encoder)._replace(spec=spec) -def pca(variables, name, scale=False, **pca_args): +def pca(variables, name=None, scale=False, **pca_args): """ Create PCA encoding of features from a sequence of variables. Additional args and *pca_args* - are passed to :py:class:`modelmatrix.transforms.PCA`. + are passed to :py:class:`sklearn.decomposition.PCA`. Parameters ---------- - variables : [column identifier, Column or Feature] - Sequence whose columns will be encoded by PCA. + 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)`. + + name : str (optional) + Defaults to one derived from the variables. + + scale : bool, default=False + Standardize the columns before the PCA? Returns ------- @@ -1579,6 +1811,11 @@ def pca(variables, name, scale=False, **pca_args): with_std=True) encoder = make_pipeline(scaler, encoder) + if name is None: + if isinstance(variables, Wildcard): + name = repr(variables) + else: + name = ', '.join(getattr(v, 'name', None) or str(v) for v in variables) _args = _argstring(**pca_args) if _args: @@ -1624,7 +1861,7 @@ def _rebuild(feature, **changes): helper, args, kwargs = feature.spec kwargs = dict(kwargs) kwargs.update(changes) - if helper == 'pca': + if helper == 'pca' and isinstance(args[0], tuple): args = (list(args[0]),) + tuple(args[1:]) return {'poly': poly, 'ns': ns, 'bs': bs, 'pca': pca}[helper](*args, **kwargs) diff --git a/modelmatrix/response.py b/modelmatrix/response.py new file mode 100644 index 0000000..1a6fd84 --- /dev/null +++ b/modelmatrix/response.py @@ -0,0 +1,381 @@ +""" +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 + +- `Surv`: a survival response, as R's `Surv(time, status)` or + `Surv(start, stop, status)`; +- `Binomial`: a binary response, or successes out of a number of + trials; +- `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. +""" + +import numpy as np, pandas as pd + + +class Response(object): + + """ + Base class of the responses: the columns of `X` used, `fit` + and `transform`. + """ + + columns = () + + def fit(self, X): + for col in self.columns: + _values(X, col) # raises if the column is not in X + return self + + def transform(self, X): + raise NotImplementedError + + def __repr__(self): + return f'{type(self).__name__}({", ".join(map(repr, self.columns))})' + + +def _values(X, col): + # a column of X, as a pd.Series + if isinstance(X, pd.DataFrame): + if col not in X.columns: + raise ValueError(f'the response column {col!r} is not in X') + return X[col] + X = np.asarray(X) + if not isinstance(col, (int, np.integer)) or not 0 <= col < X.shape[1]: + raise ValueError(f'the response column {col!r} is not a column of X, ' + f'an array with {X.shape[1]} columns') + return pd.Series(X[:, col], name=col) + + +def _numeric(X, col): + values = _values(X, col) + try: + return values.astype(float) + except (TypeError, ValueError): + raise ValueError(f'the response column {col!r} should be numeric, ' + f'got dtype {values.dtype}') + + +def _output(D, X): + # pandas for a pd.DataFrame X, else np.ndarray + if isinstance(X, pd.DataFrame): + return D + 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): + + """ + A survival response, as R's `Surv`: `Surv(time, status)` for + right-censored data, `Surv(start, stop, status)` for (start, stop] + intervals (counting process data), or with keywords + `Surv(stop=..., status=..., start=...)`. + + Parameters + ---------- + + stop : column identifier + Times of the events, or of censoring. + + status : column identifier + 1 (or True) for an event, 0 (or False) for censoring. + + start : column identifier (optional) + Start of the interval at risk. + + Notes + ----- + + The response is a :py:class:`pd.DataFrame` with columns `start` + (if given), `stop` and `status`, named as in `X`, e.g. for glmstar's + `CoxFamily(event_id=stop, status_id=status, start_id=start)` or for + lifelines' `CoxPHFitter`. + """ + + def __init__(self, *args, start=None, stop=None, status=None): + if args: + if start is not None or stop is not None or status is not None: + raise TypeError('give the columns of Surv by position or by keyword, not both') + if len(args) == 2: + stop, status = args + elif len(args) == 3: + start, stop, status = args + else: + raise TypeError('Surv takes (time, status) or (start, stop, status), ' + f'got {len(args)} columns') + if stop is None or status is None: + raise TypeError('Surv needs the columns of the times (stop) and of the status') + self.start, self.stop, self.status = start, stop, status + + @property + def columns(self): + if self.start is None: + 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 transform(self, X): + status = _numeric(X, self.status) # booleans as 0 and 1 + observed = status[status.notna()] + if not observed.isin([0, 1]).all(): + raise ValueError(f'the status {self.status!r} of Surv should be 0 or 1 (or boolean), ' + f'got the values {sorted(set(observed) - {0, 1})[:5]}') + stop = _numeric(X, self.stop) + D = {} + if self.start is not None: + start = _numeric(X, self.start) + if (start >= stop).any(): + raise ValueError(f'the start {self.start!r} of Surv should be less than ' + f'its stop {self.stop!r}') + D[self.start] = start.to_numpy() + D[self.stop] = stop.to_numpy() + D[self.status] = status.to_numpy() + index = X.index if isinstance(X, pd.DataFrame) else None + return _output(pd.DataFrame(D, index=index), X) + + def __repr__(self): + if self.start is None: + return f'Surv({self.stop!r}, {self.status!r})' + return f'Surv({self.start!r}, {self.stop!r}, {self.status!r})' + + +class Binomial(Response): + + """ + A binomial response: a binary column, or the number of successes + out of a number of trials. + + Parameters + ---------- + + successes : column identifier + A binary column (0/1, boolean, or with two levels), or with + *trials*, the number of successes. + + trials : column identifier (optional) + The number of trials. + + level : level identifier (optional) + For a binary column: the level counted as a success, any other + level being a failure. Defaults to the second level, as in R: of + the categories of a categorical, else of the sorted values. + + output : {'counts', 'proportion'} + With *trials*: `'counts'` gives the successes and the failures, + two columns as R's `cbind(successes, failures)` (e.g. for + statsmodels' `GLM` with a `Binomial` family); `'proportion'` the + proportion of successes and the number of trials, as a response + and its weights (e.g. for glmstar's `response_id` and `weight_id`). + + Notes + ----- + + A binary column gives a :py:class:`pd.Series` of 0.0 and 1.0 + (missing values are `NaN`); with *trials*, a :py:class:`pd.DataFrame`. + """ + + def __init__(self, successes, trials=None, level=None, output='counts'): + self.successes = successes + self.trials = trials + self.level = level + self.output = output + + @property + def columns(self): + if self.trials is None: + return (self.successes,) + return (self.successes, self.trials) + + 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 + values = _values(X, self.successes) + observed = values[values.notna()] + if isinstance(values.dtype, pd.CategoricalDtype): + levels = list(values.cat.categories) + elif values.dtype == bool: + levels = [False, True] + elif pd.api.types.is_numeric_dtype(values.dtype): + levels = [0, 1] + if self.level is None and not observed.isin([0, 1]).all(): + raise ValueError(f'the binary response {self.successes!r} should be 0 or 1, ' + f'got the values {sorted(set(observed) - {0, 1})[:5]}') + else: + levels = sorted(observed.unique()) + if self.level is not None: + if self.level not in levels and self.level not in set(observed): + raise ValueError(f'the level {self.level!r} is not a level of {self.successes!r}') + self.level_ = self.level + elif len(levels) != 2: + raise ValueError(f'the binary response {self.successes!r} should have two levels, ' + f'got {len(levels)}: {levels[:5]}; give the level ' + 'of a success with Binomial(..., level=...)') + else: + self.level_ = levels[1] + self.levels_ = levels + return self + + def transform(self, X): + index = X.index if isinstance(X, pd.DataFrame) else None + if self.trials is None: + values = _values(X, self.successes) + missing = values.isna().to_numpy() + if self.level is None: + unseen = set(values[~missing]) - set(self.levels_) + if unseen: + raise ValueError(f'the binary response {self.successes!r} has levels ' + f'not seen at fit: {sorted(unseen, key=str)[:5]}') + y = (values == self.level_).to_numpy().astype(float) + y[missing] = np.nan + return _output(pd.Series(y, index=index, name=self.successes), X) + + s = _numeric(X, self.successes).to_numpy() + n = _numeric(X, self.trials).to_numpy() + with np.errstate(invalid='ignore'): + if ((s < 0) | (s > n)).any(): + raise ValueError(f'the successes {self.successes!r} should be between 0 ' + f'and the trials {self.trials!r}') + if self.output == 'counts': + D = pd.DataFrame({self.successes: s, + f'{self.trials}-{self.successes}': n - s}, index=index) + else: + with np.errstate(invalid='ignore', divide='ignore'): + D = pd.DataFrame({f'{self.successes}/{self.trials}': s / n, + self.trials: n}, index=index) + return _output(D, X) + + def __repr__(self): + args = [repr(self.successes)] + if self.trials is not None: + args.append(f'trials={self.trials!r}') + if self.level is not None: + args.append(f'level={self.level!r}') + if self.output != 'counts': + args.append(f'output={self.output!r}') + return f'Binomial({", ".join(args)})' + + +class Multinomial(Response): + + """ + A multinomial response: the class of each observation, or the + counts of each class. + + Parameters + ---------- + + classes : column identifier, or list of column identifiers + A column holding the class of each observation, or a column + of counts for each class. + + output : {'labels', 'onehot'} + For a column of classes: `'labels'` gives the classes as a + categorical, its categories the levels seen at fit (e.g. for a + scikit-learn classifier); `'onehot'` a column of indicators for + each level. + + Notes + ----- + + The levels are the categories of a categorical, else the sorted + values. A level not seen at fit raises an error. + """ + + def __init__(self, classes, output='labels'): + self.classes = classes + self.output = output + + @property + def columns(self): + if isinstance(self.classes, list): + return tuple(self.classes) + return (self.classes,) + + def fit(self, X): + if isinstance(self.classes, list): + self.transform(X) # checks the counts + return self + 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) + if isinstance(values.dtype, pd.CategoricalDtype): + self.levels_ = list(values.cat.categories) + else: + self.levels_ = sorted(values[values.notna()].unique()) + return self + + 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}, + index=index) + if (D < 0).any().any(): + raise ValueError(f'the counts {self.classes} of Multinomial should be nonnegative') + return _output(D, X) + values = _values(X, self.classes) + missing = values.isna().to_numpy() + codes = pd.Index(self.levels_).get_indexer(np.asarray(values, dtype=object)) + unseen = (codes < 0) & ~missing + if unseen.any(): + raise ValueError(f'the classes {self.classes!r} have levels not seen at fit: ' + f'{sorted(set(np.asarray(values)[unseen]), key=str)[:5]}') + if self.output == 'labels': + labels = pd.Categorical.from_codes(codes, categories=self.levels_) + return _output(pd.Series(labels, index=index, name=self.classes), X) + D = np.equal.outer(codes, np.arange(len(self.levels_))).astype(float) + D[missing] = np.nan + return _output(pd.DataFrame(D, index=index, + columns=[f'{self.classes}[{l}]' for l in self.levels_]), X) + + def __repr__(self): + if self.output != 'labels': + return f'Multinomial({self.classes!r}, output={self.output!r})' + return f'Multinomial({self.classes!r})' diff --git a/tests/test_model_spec_hard.py b/tests/test_model_spec_hard.py index a6853b3..479ca61 100644 --- a/tests/test_model_spec_hard.py +++ b/tests/test_model_spec_hard.py @@ -384,6 +384,7 @@ 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 new file mode 100644 index 0000000..28f9c78 --- /dev/null +++ b/tests/test_response.py @@ -0,0 +1,196 @@ +# Responses: ModelSpec(..., response=...) and build_response + +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 +from sklearn.pipeline import make_pipeline + +from modelmatrix import ModelSpec as MS, wildcard, Surv, Binomial, Multinomial + +from _datasets import load_data + + +@pytest.fixture(scope='module') +def Carseats(): + return load_data('Carseats') + + +@pytest.fixture +def S(): + rng = np.random.default_rng(1) + n = 40 + start = rng.uniform(size=n) + return pd.DataFrame({'start': start, + 'time': start + rng.exponential(size=n), + 'status': rng.binomial(1, 0.7, n), + 'x': rng.standard_normal(n), + 'cases': rng.binomial(10, 0.3, n), + 'size': np.full(n, 10), + 'flag': rng.binomial(1, 0.5, n).astype(bool), + 'kind': pd.Categorical(rng.choice(['lo', 'mid', 'hi'], n), + 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): + 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) + 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') + 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): + spec = MS(['Price']) + spec.fit(Carseats, Carseats['Sales']) + with pytest.raises(ValueError, match='indices of X and y do not agree'): + spec.fit(Carseats, Carseats['Sales'].reset_index(drop=True).iloc[::-1]) + # arrays are not checked + spec.fit(Carseats, np.asarray(Carseats['Sales'])) + + +# Surv + +def test_surv(S): + y = MS(['x'], response=Surv('time', 'status')).fit(S).build_response(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) + assert list(y3.columns) == ['start', 'time', 'status'] + y_kw = MS(['x'], response=Surv(stop='time', status='status', start='start')).fit(S).build_response(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) + 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) + 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) + + +def test_surv_errors(S): + with pytest.raises(TypeError, match='by position or by keyword'): + Surv('time', 'status', start='start') + 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) + bad = S.assign(start=S['time'] + 1) + with pytest.raises(ValueError, match='less than'): + MS(['x'], response=Surv('start', 'time', 'status')).fit(bad) + + +# Binomial + +def test_binomial_binary(Carseats): + y = MS(['Price'], response=Binomial('US')).fit(Carseats).build_response(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) + 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)) + with pytest.raises(ValueError, match='should be 0 or 1'): + MS(['x'], response=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) + 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]) + with pytest.raises(ValueError, match='not seen at fit'): + spec.build_response(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) + 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) + 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) + with pytest.raises(ValueError, match="output should be"): + MS(['x'], response=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) + 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'))) + + +def test_multinomial_onehot_and_counts(S): + Y = MS(['x'], response=Multinomial('kind', output='onehot')).fit(S).build_response(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) + 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) + 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) + assert isinstance(y, np.ndarray) and y.shape == (40, 2) + np.testing.assert_allclose(y, A[:, :2]) diff --git a/tests/test_wildcard.py b/tests/test_wildcard.py new file mode 100644 index 0000000..2e7c565 --- /dev/null +++ b/tests/test_wildcard.py @@ -0,0 +1,174 @@ +# The wildcard: all the columns of X but the response, as R's `.` + +from itertools import combinations + +import numpy as np, pandas as pd +import pytest +from sklearn.base import clone +from sklearn.linear_model import LinearRegression +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 _datasets import load_data + + +@pytest.fixture(scope='module') +def Carseats(): + return load_data('Carseats') + + +@pytest.fixture +def D(): + rng = np.random.default_rng(0) + n = 30 + return pd.DataFrame({'y': rng.standard_normal(n), + 'x': rng.standard_normal(n), + 'z': rng.uniform(size=n), + 'g': pd.Categorical(rng.choice(['a', 'b', 'c'], n))}) + + +def _same(A, B): + assert list(A.columns) == list(B.columns) + np.testing.assert_allclose(A, B) + + +def test_all_but_response(Carseats): + others = list(Carseats.columns.drop('Sales')) + spec = MS([wildcard], response='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 + + +def test_no_response_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) + assert spec.expanded_terms_ == ['x', 'z', 'g'] + spec = MS([wildcard, 'x', ('x', 'z')], response='y').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'] + + +def test_pow(D): + # R's y ~ .^2 + spec = MS([wildcard(pow=2)], response='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) + 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) + 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 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) + assert spec.expanded_terms_ == [('x', 'g'), ('z', 'g')] + p = poly('x', 2) + spec = MS([(p, wildcard(exclude='x'))], response='y').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) + + +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) + 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 + 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) + 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()) + grid = GridSearchCV(pipe, {'modelspec__terms__0__n_components': [1, 3]}, cv=3) + grid.fit(Carseats, Carseats['Sales']) + best = grid.best_estimator_[0] + assert best.terms_[0].encoder.n_components == grid.best_params_['modelspec__terms__0__n_components'] + + +def test_transform_without_response(D): + spec = MS([wildcard, ('x', 'g')], response='y').fit(D) + _same(spec.transform(D.drop(columns='y')), spec.transform(D)) + + +def test_clone_and_cross_validate(D): + spec = MS([wildcard], response='y') + cloned = clone(spec) + assert cloned.terms == [wildcard] + 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_ndarray(D): + A = np.asarray(D[['y', 'x', 'z']]) + spec = MS([wildcard], response=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) + assert spec.expanded_terms_ == [2, 1] + assert all(type(t) is int for t in spec.expanded_terms_) + + +def test_polars(D): + pl = pytest.importorskip('polars') + P = pl.from_pandas(D) + spec = MS([wildcard], response='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(): + with pytest.raises(ValueError, match='positive integer'): + Wildcard(pow=0) + with pytest.raises(ValueError, match='positive integer'): + wildcard(pow=1.5) + assert wildcard(pow=2) == Wildcard(pow=2) + assert hash(wildcard(exclude=['x'])) == hash(Wildcard(exclude=('x',)))