`statspai.selection`¶

selection ¶

Variable selection tools.

stepwise: Stepwise regression with AIC/BIC/p-value criteria
lasso_select: LASSO-based variable selection with coordinate descent
best_subset: exact best-subset selection by branch and bound
glmnet: lasso / ridge / elastic net with the conventions of R glmnet (penalty path, penalty factors, lambda.min / lambda.1se)
shrinkage: ridge / lasso / principal-components prediction, tuned by cross-validation

GlmnetResult `dataclass` ¶

Bases: ResultProtocolMixin

Result of :func:glmnet.

Attributes:

Name	Type	Description
`params`	`Series`	Coefficients at the selected penalty, on the scale of the data.
`intercept`	`float`
`penalty`	`float`	The selected `lambda`.
`lambda_min, lambda_1se`	`float`	The penalty with the smallest cross-validated error, and the largest one within a standard error of it (`nan` without cross-validation).
`path`	`DataFrame`	One row per penalty: `lambda`, `df` (non-zero coefficients), `dev_ratio` (share of deviance explained).
`coefficients`	`DataFrame`	Coefficients along the path, one row per penalty, `intercept` first.
`cv`	`DataFrame or None`	`lambda`, `cvm` (mean cross-validated error), `cvsd` (its standard error), `cvup`, `cvlo`, `nzero`.

Examples:

>>> import numpy as np, pandas as pd, statspai as sp
>>> rng = np.random.default_rng(0)
>>> X = rng.normal(size=(200, 8))
>>> df = pd.DataFrame(X, columns=[f"x{j}" for j in range(8)])
>>> df["y"] = 2 * X[:, 0] - X[:, 1] + rng.normal(size=200)
>>> fit = sp.glmnet(df, "y", [f"x{j}" for j in range(8)], seed=0)
>>> isinstance(fit, sp.GlmnetResult)
True
>>> bool(fit.lambda_1se >= fit.lambda_min)
True
>>> fit.predict(df.head(3)).shape
(3,)

coef ¶

coef(s: Union[None, str, float] = None) -> Series

Coefficients (intercept first) at s: 'lambda.min', 'lambda.1se', a penalty on the path, or the selected one.

predict ¶

predict(data: DataFrame, s: Union[None, str, float] = None, type: str = 'response') -> ndarray

Predictions at s. type='link' is the linear predictor; 'response' is the mean (a probability for a binomial fit).

ShrinkageResult `dataclass` ¶

Bases: ResultProtocolMixin

Fit returned by :func:statspai.shrinkage.

Attributes:

Name	Type	Description
`method`	`str`
`params`	`Series`	Coefficients of the standardised predictors.
`intercept`	`float`	Mean of the outcome on the estimation sample.
`penalty`	`float or None`	The ridge / lasso penalty of the reported fit.
`n_components`	`int or None`	Number of principal components of the reported fit.
`cv`	`DataFrame or None`	One row per candidate tuning value: the value, the cross-validated mean squared prediction error and its square root.
`cv_rmspe`	`float or None`	Cross-validated root MSPE of the reported fit.
`rmspe_in`	`float`	In-sample root mean squared error (no degrees-of-freedom correction).
`n_obs, n_predictors, n_nonzero`	`int`

Examples:

>>> import numpy as np, pandas as pd, statspai as sp
>>> rng = np.random.default_rng(0)
>>> X = rng.normal(size=(120, 8))
>>> df = pd.DataFrame(X, columns=[f"x{j}" for j in range(8)])
>>> df["y"] = X[:, 0] - X[:, 1] + rng.normal(size=120)
>>> fit = sp.shrinkage(df, "y", [f"x{j}" for j in range(8)], method="ridge")
>>> type(fit).__name__
'ShrinkageResult'
>>> fit.predict(df.head(2)).shape
(2,)
>>> bool(fit.cv_rmspe > 0 and fit.penalty in set(fit.cv["penalty"]))
True

predict ¶

predict(data: Optional[DataFrame] = None) -> ndarray

Predicted outcome for the rows of data (which must hold the predictors, in their original units).

rmspe ¶

rmspe(data: DataFrame) -> float

Root mean squared prediction error on data, a frame holding the outcome and the predictors (a hold-out sample).

plot ¶

plot(ax: Any = None) -> Any

Cross-validated RMSPE against the tuning parameter.

SelectionResult `dataclass` ¶

Bases: ResultProtocolMixin

Result container for variable selection procedures.

Attributes:

Name	Type	Description
`selected`	`list[str]`	Variable names retained in the final model.
`dropped`	`list[str]`	Variable names excluded from the final model.
`history`	`DataFrame`	Step-by-step log of the selection procedure.
`final_model`	`dict`	Summary statistics for the final model (R², adj-R², AIC, BIC, n, k).
`method`	`str`	Selection method used (`"forward"`, `"backward"`, `"both"`, `"lasso_cv"`, `"lasso_bic"`, `"lasso_aic"`).
`coefficients`	`dict \| None`	Coefficient estimates from the final model (name -> value).
`lasso_path`	`dict \| None`	For LASSO: `{"lambdas": array, "coef_paths": dict}` storing the regularisation path.

Examples:

>>> import statspai as sp
>>> df = sp.cps_wage()
>>> res = sp.stepwise(df, y="log_wage",
...                   x=["female", "education", "experience"],
...                   method="forward", criterion="aic", verbose=False)
>>> isinstance(res, sp.SelectionResult)
True
>>> sorted(res.final_model.keys())
['adj_r_squared', 'aic', 'bic', 'k', 'n', 'r_squared']

summary ¶

summary() -> str

Return a human-readable summary and print it.

plot ¶

plot(figsize: Tuple[float, float] = (10, 6)) -> Tuple[Any, Any]

Plot selection diagnostics.

For stepwise: criterion value at each step. For LASSO: coefficient path across lambda values.

Returns:

Type	Description
`(fig, ax) : matplotlib Figure and Axes`

glmnet ¶

glmnet(data: DataFrame, y: str, x: Sequence[str], *, alpha: float = 1.0, lambda_: Union[None, float, Sequence[float]] = None, family: str = 'gaussian', penalty_factor: Optional[Sequence[float]] = None, cv: bool = True, n_folds: int = 10, foldid: Union[None, str, Sequence[int]] = None, rule: str = 'min', n_lambda: int = 100, lambda_min_ratio: Optional[float] = None, standardize: bool = True, intercept: bool = True, seed: Optional[int] = None, tol: float = 1e-18, max_iter: int = 100000) -> GlmnetResult

Lasso, ridge and elastic net with the conventions of R glmnet.

Minimises, over the intercept and the coefficients,

-loglik / n + lambda * sum_j v_j (alpha |b_j| + (1 - alpha) b_j^2 / 2)

(for a Gaussian outcome -loglik / n is RSS / (2 n)) along a decreasing path of lambda, and chooses lambda by K-fold cross-validation.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`	Rows with a missing outcome or predictor are dropped.	required
`y`	`str`	Outcome. For `family='binomial'` it must be coded 0/1.	required
`x`	`sequence of str`	Predictors, in their original units.	required
`alpha`	`float`	Mixing weight: 1 is the lasso, 0 the ridge.	`1.0`
`lambda_`	`float or sequence of float`	Penalties to fit. `None` builds `glmnet`'s path: `n_lambda` values, log-spaced from the smallest penalty that sets every penalised coefficient to zero down to `lambda_min_ratio` of it, stopped early once the deviance explained no longer improves.	`None`
`family`	`('gaussian', 'binomial')`		`'gaussian'`
`penalty_factor`	`sequence of float`	Multiplier `v_j` of the penalty of each predictor. 0 leaves a predictor unpenalised; weights `1 / \|b_j\|` from a first-stage fit give the adaptive lasso. Rescaled to sum to the number of predictors, as in `glmnet`.	`None`
`cv`	`bool`	Cross-validate. With `False` the last penalty of the path is selected and `lambda_min` / `lambda_1se` are `nan`.	`True`
`n_folds`	`int`		`10`
`foldid`	`sequence of int or str`	Fold label of each row (or a column holding it). Without it rows are assigned at random with `seed`.	`None`
`rule`	`('min', '1se')`	Which penalty `params`, `coef()` and `predict()` use by default.	`'min'`
`n_lambda`	`int`		`100`
`lambda_min_ratio`	`float`	Default `1e-4` with more rows than predictors, else `0.01`.	`None`
`standardize`	`bool`	Scale predictors to unit standard deviation (divisor `n`) before penalising. Coefficients are always returned in original units.	`True`
`intercept`	`bool`		`True`
`seed`	`int`		`None`
`tol`	`float`	Convergence threshold of coordinate descent, on the largest squared coefficient change (standardised units). Much tighter than `glmnet`'s default `thresh = 1e-7`, so that the fit is the minimiser to about nine digits. `glmnet` at its default is up to 6e-4 from that in the coefficients of the test designs; at `thresh = 1e-14` it agrees to 1e-6.	`1e-18`
`max_iter`	`int`		`100000`

Returns:

Type	Description
`GlmnetResult`

Notes

Conventions shared with glmnet:

Predictors are scaled by the standard deviation with divisor n.
A Gaussian outcome is scaled to unit variance before the penalty is applied and lambda is reported on the original scale. The lasso part does not notice; the ridge part does, so in original units the ridge term is lambda / sd(y) * (1 - alpha) / 2 * b^2. This is why glmnet(alpha = 0) is not the textbook ridge at the same lambda.
The first point of a path built here is the fit with an infinite penalty, and the path stops early once the deviance explained no longer improves. For the ridge, where no finite penalty sets the coefficients to zero, the first point is therefore the null model (coefficients of order 1e-36, all counted in df) although it carries a finite label.
Cross-validation refits the whole path, standardisation included, on each training set. With a path built here every training set builds its own and is read at the full-sample penalties by linear interpolation in lambda; with lambda_ supplied every training set is fitted at exactly those penalties. The error is the mean squared error (Gaussian) or the binomial deviance with probabilities held in [1e-5, 1 - 1e-5], averaged within fold; cvsd is the standard error of the fold means, weighted by fold size.

:func:statspai.shrinkage is the other penalised regression in the package: ridge, lasso and principal components with the penalty stated on the residual sum of squares, for comparing predictors by cross-validated error. Use this function when the penalty has to mean what it means in glmnet, or for the elastic net, a binomial outcome, penalty factors, or lambda.1se.

Inference after selection is a different problem: see :func:statspai.rlasso_effect and :func:statspai.dml.

Examples:

>>> import numpy as np, pandas as pd, statspai as sp
>>> rng = np.random.default_rng(0)
>>> X = rng.normal(size=(200, 8))
>>> cols = [f"x{j}" for j in range(8)]
>>> df = pd.DataFrame(X, columns=cols)
>>> df["y"] = 2 * X[:, 0] - X[:, 1] + rng.normal(size=200)
>>> fit = sp.glmnet(df, "y", cols, alpha=0.5, seed=0)
>>> bool(fit.params["x0"] > 1.5)
True
>>> sparse = fit.coef("lambda.1se")
>>> bool((sparse != 0).sum() <= (fit.coef("lambda.min") != 0).sum())
True

References

[@friedman2010regularization], [@zou2005regularization], [@zou2006adaptive]

lasso_select ¶

lasso_select(data: DataFrame, y: str, x: list[str], method: Literal['cv', 'bic', 'aic'] = 'cv', n_folds: int = 10, n_lambda: int = 100, eps: float = 0.001, max_iter: int = 1000, tol: float = 1e-06, verbose: bool = True, seed: int = 42) -> SelectionResult

LASSO-based variable selection.

Solves min (1/2n)||y - Xβ||² + λ||β||₁ via coordinate descent. Selects λ by K-fold cross-validation, BIC, or AIC.

With method="cv" every training fold is centred and scaled on its own rows and the error is pooled over observations. On the same folds and penalties this is the cross-validated error of sp.glmnet (lasso_path["cv_mse"]). In a small sample that error curve can be nearly flat, and the selected λ then changes with the fold assignment; compare a few values of seed before reading much into the list of selected variables. sp.glmnet also reports lambda.1se.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`		required
`y`	`str`	Dependent variable column.	required
`x`	`list[str]`	Candidate independent variables.	required
`method`	`('cv', 'bic', 'aic')`	How to choose the regularisation parameter λ.	`"cv"`
`n_folds`	`int`	Number of cross-validation folds (only for `method="cv"`).	`10`
`n_lambda`	`int`	Number of λ values in the grid.	`100`
`eps`	`float`	Ratio of lambda_min / lambda_max.	`0.001`
`max_iter`	`int`	Maximum coordinate descent iterations per λ.	`1000`
`tol`	`float`	Convergence tolerance.	`1e-06`
`verbose`	`bool`	Print progress.	`True`
`seed`	`int`	Random seed for CV fold assignment.	`42`

Returns:

Type	Description
`SelectionResult`

Examples:

>>> import statspai as sp
>>> df = sp.cps_wage()
>>> res = sp.lasso_select(
...     df, y="log_wage",
...     x=["female", "education", "experience", "tenure", "union"],
...     method="bic", verbose=False)
>>> res.method
'lasso_bic'
>>> bool(res.lasso_path["lambda_best"] > 0)
True

statspai.selection¶

selection ¶

GlmnetResult dataclass ¶

coef ¶

predict ¶

ShrinkageResult dataclass ¶

predict ¶

rmspe ¶

plot ¶

SelectionResult dataclass ¶

summary ¶

plot ¶

glmnet ¶

lasso_select ¶

`statspai.selection`¶

GlmnetResult `dataclass` ¶

ShrinkageResult `dataclass` ¶

SelectionResult `dataclass` ¶