model

Unified statistical model with simulation-first support.

Provides a single interface for lm, glm, lmer, and glmer. Model type is inferred from formula structure and family parameter.

Parameters:

Notes: Formula Operators

In-Formula Transforms

All stateful transforms (center, norm, zscore, scale, rank, signed_rank) learn parameters from training data and reapply them on new data during .predict(). Nesting is supported: zscore(rank(x)).

Contrast Coding (In-Formula)

Set a custom reference level or ordering with keyword args: treatment(group, ref=B), sum(group, omit=A), helmert(group, [low, med, high]). ::

model("y ~ sum(group) + x", data)                    # sum coding
model("y ~ treatment(group, ref=B) + x", data)       # treatment with ref="B"
model("y ~ group + x", data, contrasts={"group": M}) # custom ndarray matrix

Attributes:

Methods:

Attributes¶

contrasts¶

contrasts: dict | None = field(default=None, repr=False)

Custom contrast matrices mapping column names to ndarray matrices.

data¶

data: pl.DataFrame | str | None = field(default=None, repr=False, converter=_coerce_data, validator=(validators.optional(validators.instance_of(pl.DataFrame))), on_setattr=(setters.frozen))

Input data as a Polars DataFrame, file path, or None for simulation.

designmat¶

designmat: pl.DataFrame

Fixed-effects design matrix as a named DataFrame.

Returns:

diagnostics¶

diagnostics: pl.DataFrame

Model-level goodness-of-fit diagnostics.

Returns a single-row DataFrame with metrics that vary by model type: df, AIC, BIC, loglik, R-squared (OLS), deviance/dispersion (GLM), sigma/ICC (mixed). After .infer(how='cv'), adds cv_error.

Returns:

effects¶

effects: pl.DataFrame

Marginal effects or estimated marginal means (EMMs) table.

Returns:

family¶

family: str = field(default='gaussian', validator=(validators.in_(FAMILIES)))

Response distribution ("gaussian", "binomial", "poisson", "gamma", "tdist").

formula¶

formula: str = field(on_setattr=(setters.frozen))

R-style model formula string (e.g., "y ~ x + (1|group)").

link¶

link: str | None = field(default=None)

Link function (e.g., "identity", "logit", "log"). None uses canonical link.

metadata¶

metadata: pl.DataFrame

Model metadata: observation counts, parameter count, group counts.

Returns a single-row DataFrame with structural info about the model and data. For mixed models, includes ngroups.

Returns:

method¶

method: str | None = field(default=None)

Estimation method ("ols", "ml", "reml"). None auto-selects.

missing¶

missing: str = field(default='drop', validator=(validators.in_(['drop', 'fail'])))

Missing value handling: "drop" (default) or "fail".

params¶

params: pl.DataFrame

Coefficient estimates table from the fitted model.

Returns:

power_results¶

power_results: pl.DataFrame

Power analysis results from simulate(power=...).

Returns:

predictions¶

predictions: pl.DataFrame

Predictions table from the fitted model.

Returns:

resamples¶

resamples: pl.DataFrame | None

Raw resampled values from bootstrap or permutation inference.

Returns a long-format DataFrame with one row per (resample, term) combination:

• resample (int): Resample index (0 to n_resamples - 1).

• term (str): Parameter or effect name.

• value (float): Resampled estimate (bootstrap) or null test statistic (permutation).

Total rows = n_resamples × n_terms.

Returns None if .infer() hasn’t been called with how="boot"|"perm", or if save_resamples=False was passed.

Returns:

simulations¶

simulations: pl.DataFrame

Simulations table from data generation or response simulation.

Pre-fit mode (n=): full generated dataset with predictors and response. Post-fit mode (nsim=): sim_1, sim_2, ... columns with simulated responses.

Returns:

varying_corr¶

varying_corr: pl.DataFrame

Random effect correlations for mixed models.

Returns a tidy DataFrame of pairwise correlations between random effects. Empty (zero rows) for models with only random intercepts or uncorrelated (diagonal) RE structures.

Returns:

varying_offsets¶

varying_offsets: pl.DataFrame

Random effects (BLUPs) -- group-level deviations from population parameters.

Returns:

varying_params¶

varying_params: pl.DataFrame

Group-specific coefficients (population params + varying offsets).

Returns:

varying_spread¶

varying_spread: pl.DataFrame

Variance components for mixed models (sigma2, tau2, rho, ICC).

Returns:

Methods¶

explore¶

explore(formula: str, *, inverse_transforms: bool = True, effect_scale: Literal['link', 'response'] = 'link', varying: Literal['exclude', 'include'] = 'exclude', how: Literal['auto', 'mem', 'ame'] = 'auto', by: str | None = None) -> ModelResult

Compute marginal effects or estimated marginal means (EMMs).

Parameters:

Returns:

Examples:

EMMs and contrasts::

m = model("y ~ treatment", data).fit()
m.explore("treatment").infer()
m.explore("pairwise(treatment)").infer()

With transforms (raw variable names resolve automatically)::

m = model("y ~ center(x) + group", data).fit()
m.explore("group ~ x@[50]")  # auto-centers 50

Notes: Terminology: The focal variable (LHS) is the variable whose effect is estimated. Conditions (RHS, after ~) control the reference grid — pin covariates at values or cross the focal with levels of another variable. The @ operator sets specific values; [A - B] brackets specify contrasts.

Formula Patterns

Contrast Functions

Keyword Arguments for Scaling

filter_effects¶

filter_effects(*, terms: list[str] | str | None = None, levels: list[str] | str | None = None, contrasts: list[str] | str | None = None) -> pl.DataFrame

Filter effects DataFrame by term, level, or contrast.

Parameters:

Returns:

Examples:

>>> m.explore("treatment").infer()
>>> m.filter_effects(levels=["A", "B"])
>>> m.explore("pairwise(treatment)").infer()
>>> m.filter_effects(contrasts="A - B")

filter_params¶

filter_params(terms: list[str] | str) -> pl.DataFrame

Filter params to specific coefficient terms.

Parameters:

Returns:

Examples:

>>> m.filter_params("x")
>>> m.filter_params(["x", "z"])

filter_significant¶

filter_significant(alpha: float = 0.05) -> pl.DataFrame

Filter params to statistically significant results.

Parameters:

Returns:

Examples:

>>> m.fit().infer().filter_significant()
>>> m.filter_significant(0.01)

fit¶

fit(*, weights: str | None = None, offset: str | None = None, nAGQ: int | None = None, **kwargs: object) -> ModelResult

Fit the model to data.

Parameters:

Returns:

Examples:

Fit and inspect::

m = model("y ~ x", data).fit()
m.params

Chain with inference::

m = model("y ~ x", data).fit().infer()

infer¶

infer(how: Literal['asymp', 'boot', 'perm', 'cv', 'profile', 'joint'] = 'asymp', conf_level: float | int | str = 0.95, errors: Literal['iid', 'HC0', 'HC1', 'HC2', 'HC3', 'hetero', 'unequal_var'] = 'iid', null: float = 0.0, alternative: Literal['two-sided', 'greater', 'less'] = 'two-sided', *, n_boot: int = 1000, n_perm: int = 1000, ci_type: Literal['bca', 'percentile'] = 'bca', seed: int | None = None, n_jobs: int = 1, save_resamples: bool = True, k: int = 10, n_steps: int = 20, verbose: bool = False, threshold: float | None = None, profile: bool = True, holdout_group: str | None = None) -> ModelResult

Augment current results with statistical inference.

Operates on the last operation: .fit() -> params, .explore() -> effects, .predict() -> predictions, .simulate() -> simulations.

Use how="joint" for ANOVA-style joint hypothesis tests (F or chi-squared per term).

Which how methods apply to which operations::

Operation     asymp   boot   perm   cv   profile   joint
─────────     ─────   ────   ────   ──   ───────   ─────
.fit()          ✓      ✓      ✓     ✓
.explore()      ✓      ✓      ✓                     ✓
.predict()      ✓      ✓            ✓
.simulate()     ✓

profile is a special case: it operates on mixed-model variance components (varying_spread), not on the last operation’s results.

Parameters:

Returns:

jointtest¶

jointtest(terms: list[str] | None = None, *, errors: Literal['iid', 'HC0', 'HC1', 'HC2', 'HC3', 'hetero', 'unequal_var'] = 'iid') -> pl.DataFrame

Compute ANOVA-style joint hypothesis tests for model terms.

Convenience method that auto-fits if needed, then computes joint (Type III) F-tests (gaussian) or chi-square tests (non-gaussian) for each model term.

Parameters:

Returns:

Examples:

Auto-fits and returns ANOVA table::

m = model("y ~ x1 + x2", data)
m.jointtest()

Welch ANOVA (per-term Satterthwaite df)::

m = model("y ~ factor(group)", data)
m.jointtest(errors="unequal_var")

plot_design¶

plot_design(**kwargs: object) -> object

Plot design matrix as an annotated heatmap.

Parameters:

Returns:

plot_explore¶

plot_explore(specs: str, **kwargs: object) -> object

Plot marginal effects or estimated marginal means.

Parameters:

Returns:

plot_params¶

plot_params(**kwargs: object) -> object

Plot parameter estimates as a forest plot.

Parameters:

Returns:

plot_predict¶

plot_predict(term: str, **kwargs: object) -> object

Plot marginal predictions across a predictor range.

Accepts either a bare column name or an explore-style formula::

m.plot_predict("age")              # bare column
m.plot_predict("age ~ sex")        # hue by sex
m.plot_predict("age ~ sex@Female") # pin sex=Female
m.plot_predict("age@range(5)")     # 5-point grid
m.plot_predict("age@[25,50,75]")   # explicit grid values

Parameters:

Returns:

plot_profile¶

plot_profile(**kwargs: object) -> object

Plot profile likelihood curves.

Requires .infer(how="profile") to have been called first.

Parameters:

Returns:

plot_ranef¶

plot_ranef(**kwargs: object) -> object

Plot random effect estimates.

Parameters:

Returns:

plot_relationships¶

plot_relationships(**kwargs: object) -> object

Plot pairwise relationships between response and predictors.

Parameters:

Returns:

plot_resamples¶

plot_resamples(**kwargs: object) -> object

Plot distribution of resampled statistics.

Parameters:

Returns:

plot_resid¶

plot_resid(**kwargs: object) -> object

Plot residual diagnostics (4-panel grid).

Parameters:

Returns:

plot_vif¶

plot_vif(**kwargs: object) -> object

Plot VIF diagnostics as correlation heatmap.

Parameters:

Returns:

predict¶

predict(newdata: str | pl.DataFrame | None = None, type: Literal['response', 'link'] = 'response', *, varying: Literal['exclude', 'include'] = 'exclude', allow_new_levels: bool = False, n_points: int | Literal['data'] = 50, **kwargs: object) -> ModelResult

Generate predictions from the fitted model.

Accepts a formula string, a DataFrame, or None:

• Formula mode (str): Builds a prediction grid from an explore-style formula and returns grid columns prepended to fitted values — like R’s ggpredict() / effects package.

• DataFrame mode (pl.DataFrame): Predicts on the given data.

• None: Predicts on the training data.

Parameters:

Returns:

Examples:

Formula-mode predictions (like R ggpredict)::

m = model("mpg ~ wt + hp + cyl", mtcars).fit()
m.predict("wt ~ cyl").predictions
# Returns: wt | cyl | fitted

m.predict("wt ~ cyl").infer().predictions
# Returns: wt | cyl | fitted | se | ci_lower | ci_upper

m.predict("wt@range(10) ~ cyl@[4, 6]").predictions
m.predict("wt", n_points="data").predictions

DataFrame predictions (unchanged)::

m.predict(new_df).predictions

Training-data predictions (unchanged)::

m.predict().predictions

reset_contrasts¶

reset_contrasts() -> Self

Reset contrast coding to defaults (treatment/dummy coding).

.. deprecated:: Create a new model instead of resetting contrasts in place.

Clears any custom contrasts set via set_contrasts(). The model must be re-fitted after resetting.

Returns:

set_contrasts¶

set_contrasts(**contrasts: object) -> Self

Set custom contrast coding for categorical predictors.

.. deprecated:: Use formula syntax (sum(x), treatment(x, ref=B)) or the contrasts= constructor kwarg for ndarray matrices instead.

Parameters:

Returns:

set_display¶

set_display(enabled: bool) -> None

Toggle automatic result display in the REPL and notebooks.

When enabled (the default), evaluating a model object after calling .fit(), .explore(), .predict(), or .simulate() automatically shows the relevant result table below the model header. Disable for a compact one-line repr.

Parameters:

show_math¶

show_math(*, explanations: bool = True) -> 'MathDisplay'

Display structural LaTeX equation with term explanations.

Works before or after .fit(). Before data, shows generic symbolic form. With data, shows specific factor levels and transform parameters.

Parameters:

Returns:

Examples:

>>> m = model("y ~ x + group", data)
>>> m.show_math()           # renders in Jupyter
>>> m.show_math().to_latex() # raw LaTeX string

simulate¶

simulate(n: int | None = None, nsim: int | None = None, seed: int | None = None, coef: dict[str, float] | None = None, sigma: float = 1.0, varying: str = 'fitted', power: int | dict[str, Any] | None = None, **var_specs: object) -> ModelResult

Generate data from scratch or simulate responses from a fitted model.

Pre-fit mode (n=): generate a dataset from formula and distributions. Post-fit mode (nsim=): simulate response vectors from the fitted model. Power mode (power=): run simulation-based power analysis.

Parameters:

Returns:

summary¶

summary(style: Literal['r', 'compact'] = 'r', digits: int = 3) -> FormattedText

Generate a formatted model summary (R-style or compact).

Parameters:

Returns:

to_effect_size¶

to_effect_size(*, include_intercept: bool = False) -> pl.DataFrame

Compute standardized effect sizes from params.

Computes Cohen’s d, semi-partial r, eta-squared, and odds ratio (for binomial models).

Parameters:

Returns:

Examples:

>>> m = model("y ~ x1 + x2", data).fit().infer()
>>> m.to_effect_size()

to_odds_ratio¶

to_odds_ratio() -> pl.DataFrame

Transform params to odds ratio scale (binomial GLM only).

Exponentiates estimate, ci_lower, and ci_upper columns.

Returns:

Examples:

>>> m = model("y ~ x", data, family="binomial").fit().infer()
>>> m.to_odds_ratio()

to_response_scale¶

to_response_scale() -> pl.DataFrame

Transform effects from link scale to response scale.

Applies the inverse link function to estimate and CI columns on the effects DataFrame. For example, converts log-odds to probabilities for logistic models.

Returns:

Examples:

>>> m = model("y ~ x", data, family="binomial").fit()
>>> m.explore("x").infer()
>>> m.to_response_scale()

vif¶

vif() -> pl.DataFrame

Compute variance inflation factors for model predictors.

VIF measures multicollinearity — how much the variance of each coefficient is inflated due to correlation with other predictors. Only requires data (not fitting), since VIF is a property of the design matrix.

Returns:

Examples:

>>> m = model("y ~ x1 + x2 + x3", data)
>>> m.vif()