infer

Inference operations: asymptotic, bootstrap, permutation, cross-validation, profile, and simulation.

Call chain:

model.infer() -> dispatch_infer() -> dispatch_{params,mee,prediction}_inference()
Each dispatcher routes to the method-specific backend (asymptotic, bootstrap, etc.)

Classes:

Functions:

Modules:

Classes¶

InferResult¶

Immutable result of inference dispatch.

Each field corresponds to a model attribute that may be updated by .infer(). Fields left as None are not modified.

Attributes:

Attributes¶

cv¶

cv: CVState | None = None

last_op¶

last_op: str | None = None

mee¶

mee: MeeState | JointTestState | None = None

params_inference¶

params_inference: InferenceState | None = None

pred¶

pred: PredictionState | None = None

profile¶

profile: ProfileState | None = None

resamples¶

resamples: ResamplesState | None = None

sim_inference¶

sim_inference: SimulationInferenceState | None = None

varying_spread¶

varying_spread: VaryingSpreadState | None = None

Functions¶

augment_spread_with_profile_ci¶

augment_spread_with_profile_ci(spread: 'VaryingSpreadState', profile: 'ProfileState', conf_level: float) -> 'VaryingSpreadState'

Augment variance components with profile likelihood confidence intervals.

Maps profile CIs (on the SD scale) to variance component CIs by squaring, matching component names in spread.components to keys in profile.ci_sd.

Parameters:

Returns:

build_emm_reference_grid¶

build_emm_reference_grid(data: pl.DataFrame, bundle: 'DataBundle', focal_var: str) -> np.ndarray

Build reference grid X matrix for EMM computation.

Parameters:

Returns:

compute_bootstrap_params¶

compute_bootstrap_params(spec: 'ModelSpec', bundle: 'DataBundle', n_boot: int = 1000, seed: int | None = None, n_jobs: int = 1) -> np.ndarray

Generate bootstrap distribution of coefficient estimates.

Uses case resampling: resample observations with replacement, refit the model, and collect coefficient estimates.

Parameters:

Returns:

compute_bootstrap_pvalue¶

compute_bootstrap_pvalue(observed: np.ndarray, boot_samples: np.ndarray, null: float = 0.0, alternative: str = 'two-sided') -> np.ndarray

Compute bootstrap p-values.

Uses the proportion of bootstrap samples as extreme or more extreme than the observed value (relative to null).

Parameters:

Returns:

compute_cv_metrics¶

compute_cv_metrics(spec: 'ModelSpec', bundle: 'DataBundle', k: int | str = 10, seed: int | None = None, holdout_group_ids: np.ndarray | None = None) -> 'CVState'

Compute k-fold or leave-one-out cross-validation metrics.

Splits the data into k folds, trains on k-1 folds, and predicts on the held-out fold. Aggregates out-of-sample prediction metrics.

Parameters:

Returns:

compute_jackknife_coefs¶

compute_jackknife_coefs(spec: 'ModelSpec', bundle: 'DataBundle') -> np.ndarray

Compute leave-one-out jackknife coefficient estimates.

Used for BCa acceleration factor calculation.

Parameters:

Returns:

compute_mee_bootstrap¶

compute_mee_bootstrap(spec: 'ModelSpec', bundle: 'DataBundle', mee: 'MeeState', data: pl.DataFrame, conf_level: float = 0.95, n_boot: int = 1000, ci_type: str = 'bca', seed: int | None = None, null: float = 0.0, alternative: str = 'two-sided', save_resamples: bool = True) -> 'tuple[MeeState, np.ndarray | None]'

Compute bootstrap inference for marginal effects.

Parameters:

Returns:

compute_mee_permutation¶

compute_mee_permutation(spec: 'ModelSpec', bundle: 'DataBundle', fit: 'FitState', mee: 'MeeState', data: 'pl.DataFrame', conf_level: float, se_obs: np.ndarray, n_perm: int = 1000, seed: int | None = None, null: float = 0.0, alternative: str = 'two-sided', save_resamples: bool = True) -> 'tuple[MeeState, np.ndarray | None]'

Compute permutation-based inference for marginal effects.

Permutes the response variable, refits and re-explores on each permutation, and computes p-values based on the null distribution.

Parameters:

Returns:

compute_params_asymptotic¶

compute_params_asymptotic(spec: 'ModelSpec', bundle: 'DataBundle', fit: 'FitState', conf_level: float, errors: str | None = None, data: 'pl.DataFrame | None' = None, null: float = 0.0, alternative: str = 'two-sided') -> 'InferenceState'

Compute asymptotic (Wald) inference for model parameters.

Resolves vcov (robust via sandwich or standard from fit), resolves df (Welch, t, or z), and delegates to compute_coefficient_inference.

Parameters:

Returns:

compute_params_bootstrap¶

compute_params_bootstrap(spec: 'ModelSpec', bundle: 'DataBundle', fit: 'FitState', conf_level: float = 0.95, n_boot: int = 1000, ci_type: str = 'bca', seed: int | None = None, save_resamples: bool = True, n_jobs: int = 1, null: float = 0.0, alternative: str = 'two-sided') -> 'InferenceState'

Compute bootstrap inference for parameters.

Parameters:

Returns:

compute_params_bootstrap_mixed¶

compute_params_bootstrap_mixed(spec: 'ModelSpec', bundle: 'DataBundle', fit: 'FitState', conf_level: float = 0.95, n_boot: int = 1000, ci_type: str = 'bca', seed: int | None = None, save_resamples: bool = True, n_jobs: int = 1, null: float = 0.0, alternative: str = 'two-sided') -> 'InferenceState'

Compute bootstrap inference for mixed model parameters.

Uses parametric bootstrap (simulate y* from fitted model, refit) via the existing lmer_bootstrap() / glmer_bootstrap() infrastructure instead of observation-level case resampling. This properly respects the cluster structure and matches R’s bootMer(type="parametric").

BCa is not supported for mixed models (no principled cluster jackknife definition exists). Use ci_type="percentile" or "basic".

Parameters:

Returns:

compute_params_cv_inference¶

compute_params_cv_inference(spec: 'ModelSpec', bundle: 'DataBundle', data: 'pl.DataFrame', conf_level: float, *, link_override: str | None = None, k: int | str = 10, seed: int | None = None, holdout_group_ids: np.ndarray | None = None) -> tuple['InferenceState', 'CVState']

Compute CV-based parameter importance via ablation.

Operates at the source term level rather than the design matrix column level. For each source term (excluding Intercept):

1. Build an ablated formula from the remaining source terms (plus any random effects). Main-effect ablation also removes interactions that contain that effect (hierarchical principle).

2. Fit the ablated model and run k-fold CV.

3. PRE = full_metric − ablated_metric per fold.

4. Broadcast the per-term PRE to every design matrix column that belongs to that source term.

A positive PRE value means the predictor reduces prediction error (improves model fit).

Parameters:

Returns:

Examples:

>>> m = model("y ~ x + z", data).fit().infer(how="cv")
>>> m.params  # Includes pre, pre_sd columns
>>> m.diagnostics  # Includes cv_rmse, cv_mae, cv_rsquared

compute_params_permutation¶

compute_params_permutation(spec: 'ModelSpec', bundle: 'DataBundle', fit: 'FitState', conf_level: float, n_perm: int = 1000, seed: int | None = None, save_resamples: bool = True, null: float = 0.0, alternative: str = 'two-sided') -> 'InferenceState'

Compute permutation-based inference for model parameters.

Permutes the response variable to break the relationship with predictors, refits the model on each permutation, and computes p-values based on the null distribution of test statistics.

Parameters:

Returns:

compute_prediction_asymptotic¶

compute_prediction_asymptotic(pred: 'PredictionState', bundle: 'DataBundle', fit: 'FitState', spec: 'ModelSpec', conf_level: float) -> 'PredictionState'

Compute asymptotic inference for predictions via delta method.

Computes standard errors using the delta method and confidence intervals using t or z critical values depending on the model family.

Parameters:

Returns:

compute_prediction_bootstrap¶

compute_prediction_bootstrap(spec: 'ModelSpec', bundle: 'DataBundle', pred: 'PredictionState', conf_level: float = 0.95, n_boot: int = 1000, ci_type: str = 'percentile', seed: int | None = None) -> 'PredictionState'

Compute bootstrap inference for predictions.

Each bootstrap replicate refits the model on resampled training data, then runs the same prediction path (inverse link, random effects, etc.) used by the original predict() call.

Parameters:

Returns:

compute_profile_inference¶

compute_profile_inference(spec: 'ModelSpec', bundle: 'DataBundle', fit: 'FitState', conf_level: float = 0.95, *, n_steps: int = 20, verbose: bool = False, threshold: float | None = None) -> 'ProfileState'

Compute profile likelihood CIs for variance components.

Builds a deviance function adapter from the unified model’s internal state and runs profile_likelihood(). Returns a ProfileState with profile traces, splines, and CIs on the theta and SD scales.

Parameters:

Returns:

compute_satterthwaite_emm_df¶

compute_satterthwaite_emm_df(bundle: DataBundle, fit: FitState, spec: ModelSpec, L: np.ndarray) -> np.ndarray

Compute Satterthwaite denominator df for EMM contrast rows.

For each row of the contrast matrix L, computes the Satterthwaite approximation to the denominator degrees of freedom.

Uses a hybrid approach: differentiates numerically w.r.t. theta only, then fills in sigma’s contributions analytically. This produces identical results to the full numerical approach but with ~52% fewer function evaluations for single random-intercept models.

Parameters:

Returns:

compute_simulation_inference¶

compute_simulation_inference(simulations: 'pl.DataFrame', conf_level: float) -> 'SimulationInferenceState'

Compute inference for simulations.

For post-fit simulations (sim_1, sim_2, ...): Computes summary statistics across simulation replicates: - Mean and SD per observation across simulations - Quantiles (e.g., 2.5%, 97.5% for 95% interval)

Returns a basic state indicating generated data. Returns a basic state indicating generated data.

Parameters:

Returns:

dispatch_infer¶

dispatch_infer(*, how: str, conf_level: float, errors: str, null: float, alternative: str, last_op: str | None, spec: ModelSpec, bundle: DataBundle, fit: FitState, data: pl.DataFrame | None, link_override: str | None, mee: MeeState | JointTestState | None, pred: PredictionState | None, simulations: pl.DataFrame | None, varying_spread: VaryingSpreadState | None, is_mixed: bool, n_boot: int = 1000, n_perm: int = 1000, ci_type: str = 'bca', seed: int | None = None, n_jobs: int = 1, save_resamples: bool = True, k: int | str = 10, n_steps: int = 20, verbose: bool = False, threshold: float | None = None, profile_auto: bool = True, holdout_group_ids: np.ndarray | None = None) -> InferResult

Dispatch inference to the correct backend based on method and last operation.

Pure function that implements the body of model.infer(). Accepts all required state as explicit parameters and returns an InferResult with the populated fields.

Parameters:

Returns:

dispatch_mee_inference¶

dispatch_mee_inference(*, how: str, mee: MeeState, spec: ModelSpec, bundle: DataBundle, fit: FitState, data: pl.DataFrame, conf_level: float, errors: str = 'iid', null: float = 0.0, alternative: str = 'two-sided', n_boot: int = 1000, n_perm: int = 1000, ci_type: str = 'bca', seed: int | None = None, save_resamples: bool = True) -> 'tuple[MeeState, np.ndarray | None]'

Dispatch marginal effects inference to the appropriate method.

Parameters:

Returns:

dispatch_params_inference¶

dispatch_params_inference(*, how: str, spec: ModelSpec, bundle: DataBundle, fit: FitState, data: pl.DataFrame | None, conf_level: float, errors: str = 'iid', null: float = 0.0, alternative: str = 'two-sided', link_override: str | None = None, n_boot: int = 1000, n_perm: int = 1000, ci_type: str = 'bca', seed: int | None = None, n_jobs: int = 1, save_resamples: bool = True, k: int | str = 10) -> InferenceState

Dispatch parameter inference to the appropriate method.

Parameters:

Returns:

dispatch_prediction_inference¶

dispatch_prediction_inference(*, how: str, pred: PredictionState, spec: ModelSpec, bundle: DataBundle, fit: FitState, conf_level: float, n_boot: int = 1000, ci_type: str = 'percentile', seed: int | None = None, k: int | str = 10, holdout_group_ids: np.ndarray | None = None) -> tuple[PredictionState, CVState | None]

Dispatch prediction inference to the appropriate method.

Parameters:

Returns:

generate_group_kfold_splits¶

generate_group_kfold_splits(group_ids: np.ndarray, k: int, seed: int | None = None) -> list[tuple[np.ndarray, np.ndarray]]

Generate group-aware k-fold cross-validation indices.

All observations from a group stay in the same fold, preventing data leakage from shared random effects. Groups are shuffled and assigned to roughly equal-sized folds.

Parameters:

Returns:

Examples:

>>> group_ids = np.array([0, 0, 1, 1, 2, 2, 3, 3])
>>> splits = generate_group_kfold_splits(group_ids, k=2, seed=42)
>>> len(splits)
2

generate_kfold_splits¶

generate_kfold_splits(n: int, k: int, seed: int | None = None) -> list[tuple[np.ndarray, np.ndarray]]

Generate k-fold cross-validation train/test indices.

Shuffles the indices and splits into k roughly equal-sized folds. Each fold is used exactly once as the test set.

Parameters:

Returns:

Examples:

>>> splits = generate_kfold_splits(100, k=10, seed=42)
>>> len(splits)
10
>>> train_idx, test_idx = splits[0]
>>> len(test_idx)  # ~10 observations per fold
10

Modules¶

asymptotic¶

Asymptotic inference for parameters and predictions.

Provides orchestration functions for asymptotic (Wald-based) inference: choosing vcov (robust vs standard), choosing df (Welch vs t vs z), and wiring results into InferenceState/PredictionState containers.

The underlying math lives in internal/maths/inference/. This module is pure orchestration — it accepts containers and returns containers.

Functions:

Attributes¶

Classes¶

Functions¶

augment_spread_with_profile_ci¶

augment_spread_with_profile_ci(spread: 'VaryingSpreadState', profile: 'ProfileState', conf_level: float) -> 'VaryingSpreadState'

Augment variance components with profile likelihood confidence intervals.

Maps profile CIs (on the SD scale) to variance component CIs by squaring, matching component names in spread.components to keys in profile.ci_sd.

Parameters:

Returns:

compute_params_asymptotic¶

compute_params_asymptotic(spec: 'ModelSpec', bundle: 'DataBundle', fit: 'FitState', conf_level: float, errors: str | None = None, data: 'pl.DataFrame | None' = None, null: float = 0.0, alternative: str = 'two-sided') -> 'InferenceState'

Compute asymptotic (Wald) inference for model parameters.

Resolves vcov (robust via sandwich or standard from fit), resolves df (Welch, t, or z), and delegates to compute_coefficient_inference.

Parameters:

Returns:

compute_prediction_asymptotic¶

compute_prediction_asymptotic(pred: 'PredictionState', bundle: 'DataBundle', fit: 'FitState', spec: 'ModelSpec', conf_level: float) -> 'PredictionState'

Compute asymptotic inference for predictions via delta method.

Computes standard errors using the delta method and confidence intervals using t or z critical values depending on the model family.

Parameters:

Returns:

resolve_error_type¶

resolve_error_type(errors: str, family: str = 'gaussian') -> str

Normalize errors parameter to canonical type string.

Parameters:

Returns:

resolve_welch¶

resolve_welch(spec: 'ModelSpec', bundle: 'DataBundle', fit: 'FitState', data: 'pl.DataFrame | None') -> tuple[np.ndarray, np.ndarray]

Compute both Welch sandwich vcov and per-coefficient Satterthwaite df.

Shares the cell_info computation between vcov and df so the factor grouping is done only once. Returns per-coefficient df via the generalized Satterthwaite decomposition of the sandwich variance.

For Gaussian mixed models (lmer), this provides a repeated-measures Welch analog: per-cell variance of marginal residuals y - X @ beta is used instead of raw response variance, accounting for the fixed-effect structure while allowing unequal variances across factor cells.

Note: This differs from R’s approach. R’s sandwich/emmeans/coeftest use residual df (n - p) with sandwich SEs, ignoring the df correction entirely. Bossanova instead computes per-coefficient Satterthwaite df, which generalizes the classic Welch t-test to a regression framework: for a two-group y ~ factor(group) model, the slope’s df exactly matches scipy.stats.ttest_ind(equal_var=False) and R’s t.test(var.equal=FALSE).

Parameters:

Returns:

bootstrap¶

Bootstrap inference methods.

Provides helper functions for bootstrap-based inference used by the model class’s infer(how=‘boot’) method. Resampling logic is shared via build_resampled_bundle().

Functions:

Classes¶

Functions¶

build_emm_reference_grid¶

build_emm_reference_grid(data: pl.DataFrame, bundle: 'DataBundle', focal_var: str) -> np.ndarray

Build reference grid X matrix for EMM computation.

Parameters:

Returns:

compute_bootstrap_params¶

compute_bootstrap_params(spec: 'ModelSpec', bundle: 'DataBundle', n_boot: int = 1000, seed: int | None = None, n_jobs: int = 1) -> np.ndarray

Generate bootstrap distribution of coefficient estimates.

Uses case resampling: resample observations with replacement, refit the model, and collect coefficient estimates.

Parameters:

Returns:

compute_bootstrap_pvalue¶

compute_bootstrap_pvalue(observed: np.ndarray, boot_samples: np.ndarray, null: float = 0.0, alternative: str = 'two-sided') -> np.ndarray

Compute bootstrap p-values.

Uses the proportion of bootstrap samples as extreme or more extreme than the observed value (relative to null).

Parameters:

Returns:

compute_jackknife_coefs¶

compute_jackknife_coefs(spec: 'ModelSpec', bundle: 'DataBundle') -> np.ndarray

Compute leave-one-out jackknife coefficient estimates.

Used for BCa acceleration factor calculation.

Parameters:

Returns:

compute_mee_bootstrap¶

compute_mee_bootstrap(spec: 'ModelSpec', bundle: 'DataBundle', mee: 'MeeState', data: pl.DataFrame, conf_level: float = 0.95, n_boot: int = 1000, ci_type: str = 'bca', seed: int | None = None, null: float = 0.0, alternative: str = 'two-sided', save_resamples: bool = True) -> 'tuple[MeeState, np.ndarray | None]'

Compute bootstrap inference for marginal effects.

Parameters:

Returns:

compute_params_bootstrap¶

compute_params_bootstrap(spec: 'ModelSpec', bundle: 'DataBundle', fit: 'FitState', conf_level: float = 0.95, n_boot: int = 1000, ci_type: str = 'bca', seed: int | None = None, save_resamples: bool = True, n_jobs: int = 1, null: float = 0.0, alternative: str = 'two-sided') -> 'InferenceState'

Compute bootstrap inference for parameters.

Parameters:

Returns:

compute_params_bootstrap_mixed¶

compute_params_bootstrap_mixed(spec: 'ModelSpec', bundle: 'DataBundle', fit: 'FitState', conf_level: float = 0.95, n_boot: int = 1000, ci_type: str = 'bca', seed: int | None = None, save_resamples: bool = True, n_jobs: int = 1, null: float = 0.0, alternative: str = 'two-sided') -> 'InferenceState'

Compute bootstrap inference for mixed model parameters.

Uses parametric bootstrap (simulate y* from fitted model, refit) via the existing lmer_bootstrap() / glmer_bootstrap() infrastructure instead of observation-level case resampling. This properly respects the cluster structure and matches R’s bootMer(type="parametric").

BCa is not supported for mixed models (no principled cluster jackknife definition exists). Use ci_type="percentile" or "basic".

Parameters:

Returns:

compute_prediction_bootstrap¶

compute_prediction_bootstrap(spec: 'ModelSpec', bundle: 'DataBundle', pred: 'PredictionState', conf_level: float = 0.95, n_boot: int = 1000, ci_type: str = 'percentile', seed: int | None = None) -> 'PredictionState'

Compute bootstrap inference for predictions.

Each bootstrap replicate refits the model on resampled training data, then runs the same prediction path (inverse link, random effects, etc.) used by the original predict() call.

Parameters:

Returns:

cv¶

Cross-validation inference and metrics.

Provides k-fold cross-validation for model evaluation: splits data into folds, fits on training folds, predicts on held-out folds, and computes out-of-sample metrics (RMSE, MAE, R², deviance, accuracy).

For mixed models, group-aware splitting ensures all observations from a group stay in the same fold (generate_group_kfold_splits).

Functions:

Attributes¶

Classes¶

Functions¶

compute_cv_metrics¶

compute_cv_metrics(spec: 'ModelSpec', bundle: 'DataBundle', k: int | str = 10, seed: int | None = None, holdout_group_ids: np.ndarray | None = None) -> 'CVState'

Compute k-fold or leave-one-out cross-validation metrics.

Splits the data into k folds, trains on k-1 folds, and predicts on the held-out fold. Aggregates out-of-sample prediction metrics.

Parameters:

Returns:

compute_params_cv_inference¶

compute_params_cv_inference(spec: 'ModelSpec', bundle: 'DataBundle', data: 'pl.DataFrame', conf_level: float, *, link_override: str | None = None, k: int | str = 10, seed: int | None = None, holdout_group_ids: np.ndarray | None = None) -> tuple['InferenceState', 'CVState']

Compute CV-based parameter importance via ablation.

Operates at the source term level rather than the design matrix column level. For each source term (excluding Intercept):

1. Build an ablated formula from the remaining source terms (plus any random effects). Main-effect ablation also removes interactions that contain that effect (hierarchical principle).

2. Fit the ablated model and run k-fold CV.

3. PRE = full_metric − ablated_metric per fold.

4. Broadcast the per-term PRE to every design matrix column that belongs to that source term.

A positive PRE value means the predictor reduces prediction error (improves model fit).

Parameters:

Returns:

Examples:

>>> m = model("y ~ x + z", data).fit().infer(how="cv")
>>> m.params  # Includes pre, pre_sd columns
>>> m.diagnostics  # Includes cv_rmse, cv_mae, cv_rsquared

generate_group_kfold_splits¶

generate_group_kfold_splits(group_ids: np.ndarray, k: int, seed: int | None = None) -> list[tuple[np.ndarray, np.ndarray]]

Generate group-aware k-fold cross-validation indices.

All observations from a group stay in the same fold, preventing data leakage from shared random effects. Groups are shuffled and assigned to roughly equal-sized folds.

Parameters:

Returns:

Examples:

>>> group_ids = np.array([0, 0, 1, 1, 2, 2, 3, 3])
>>> splits = generate_group_kfold_splits(group_ids, k=2, seed=42)
>>> len(splits)
2

generate_kfold_splits¶

generate_kfold_splits(n: int, k: int, seed: int | None = None) -> list[tuple[np.ndarray, np.ndarray]]

Generate k-fold cross-validation train/test indices.

Shuffles the indices and splits into k roughly equal-sized folds. Each fold is used exactly once as the test set.

Parameters:

Returns:

Examples:

>>> splits = generate_kfold_splits(100, k=10, seed=42)
>>> len(splits)
10
>>> train_idx, test_idx = splits[0]
>>> len(test_idx)  # ~10 observations per fold
10

resolve_holdout_group_ids¶

resolve_holdout_group_ids(holdout_group: str | None, data: 'pl.DataFrame | None', spec: 'ModelSpec', bundle: 'DataBundle') -> np.ndarray | None

Resolve holdout group IDs for group-aware CV splitting.

Logic:

1. If holdout_group is provided, look up the column in data.

2. If any model has random effects:

   • Single or nested RE: auto-use the primary grouping factor.

   • Crossed REs: raise with a helpful error listing available vars.

3. No RE and no holdout_group: return None (use naive splits).

Parameters:

Returns:

dispatch¶

Infer dispatch: routes .infer() calls to the correct inference backend.

Extracts the dispatch logic from model/core.py into a pure function. Accepts all state objects as explicit params, returns an InferResult container with populated fields.

Classes:

Functions:

Classes¶

InferResult¶

Immutable result of inference dispatch.

Each field corresponds to a model attribute that may be updated by .infer(). Fields left as None are not modified.

Attributes:

Attributes¶

cv¶

cv: CVState | None = None

last_op¶

last_op: str | None = None

mee¶

mee: MeeState | JointTestState | None = None

params_inference¶

params_inference: InferenceState | None = None

pred¶

pred: PredictionState | None = None

profile¶

profile: ProfileState | None = None

resamples¶

resamples: ResamplesState | None = None

sim_inference¶

sim_inference: SimulationInferenceState | None = None

varying_spread¶

varying_spread: VaryingSpreadState | None = None

Functions¶

dispatch_infer¶

dispatch_infer(*, how: str, conf_level: float, errors: str, null: float, alternative: str, last_op: str | None, spec: ModelSpec, bundle: DataBundle, fit: FitState, data: pl.DataFrame | None, link_override: str | None, mee: MeeState | JointTestState | None, pred: PredictionState | None, simulations: pl.DataFrame | None, varying_spread: VaryingSpreadState | None, is_mixed: bool, n_boot: int = 1000, n_perm: int = 1000, ci_type: str = 'bca', seed: int | None = None, n_jobs: int = 1, save_resamples: bool = True, k: int | str = 10, n_steps: int = 20, verbose: bool = False, threshold: float | None = None, profile_auto: bool = True, holdout_group_ids: np.ndarray | None = None) -> InferResult

Dispatch inference to the correct backend based on method and last operation.

Pure function that implements the body of model.infer(). Accepts all required state as explicit parameters and returns an InferResult with the populated fields.

Parameters:

Returns:

formula_utils¶

Formula utility functions shared across operations modules.

Functions:

Functions¶

build_ablated_formula¶

build_ablated_formula(response_var: str, source_terms: list[str], ablate_term: str, random_terms: tuple[str, ...] = ()) -> str

Build a formula with one source term ablated.

For main effects, also removes interactions containing that term (hierarchical principle). For interactions, removes only that interaction.

Parameters:

Returns:

Examples:

>>> build_ablated_formula("y", ["x", "z", "x:z"], "x")
'y ~ z'
>>> build_ablated_formula("y", ["x", "z", "x:z"], "x:z")
'y ~ x + z'
>>> build_ablated_formula("y", ["x"], "x")
'y ~ 1'
>>> build_ablated_formula("y", ["x", "z"], "x", ("(1|group)",))
'y ~ z + (1|group)'

parse_value¶

parse_value(s: str) -> float | str

Parse a single value as float or string.

Tries to parse as float first; falls back to unquoted or quoted string.

Parameters:

Returns:

Examples:

>>> parse_value("3.14")
3.14
>>> parse_value("hello")
'hello'
>>> parse_value("'quoted'")
'quoted'

remove_predictor_from_formula¶

remove_predictor_from_formula(formula: str, predictor: str) -> str

Remove a predictor term from a formula string.

Handles various predictor formats including raw names, factor(), C(), c().

Parameters:

Returns:

Examples:

remove_predictor_from_formula("y ~ a + b + c", "b")
# 'y ~ a + c'
remove_predictor_from_formula("y ~ factor(a) + b", "a")
# 'y ~ b'

mee¶

Marginal effects inference dispatch.

Routes MEE inference requests to the appropriate method: asymptotic (delta method), bootstrap, or permutation.

Functions:

Classes¶

Functions¶

dispatch_mee_inference¶

dispatch_mee_inference(*, how: str, mee: MeeState, spec: ModelSpec, bundle: DataBundle, fit: FitState, data: pl.DataFrame, conf_level: float, errors: str = 'iid', null: float = 0.0, alternative: str = 'two-sided', n_boot: int = 1000, n_perm: int = 1000, ci_type: str = 'bca', seed: int | None = None, save_resamples: bool = True) -> 'tuple[MeeState, np.ndarray | None]'

Dispatch marginal effects inference to the appropriate method.

Parameters:

Returns:

params¶

Parameter inference dispatch.

Routes parameter inference requests to the appropriate method: asymptotic (Wald-based), bootstrap, permutation, or cross-validation.

Functions:

Classes¶

Functions¶

dispatch_params_inference¶

dispatch_params_inference(*, how: str, spec: ModelSpec, bundle: DataBundle, fit: FitState, data: pl.DataFrame | None, conf_level: float, errors: str = 'iid', null: float = 0.0, alternative: str = 'two-sided', link_override: str | None = None, n_boot: int = 1000, n_perm: int = 1000, ci_type: str = 'bca', seed: int | None = None, n_jobs: int = 1, save_resamples: bool = True, k: int | str = 10) -> InferenceState

Dispatch parameter inference to the appropriate method.

Parameters:

Returns:

permutation¶

Permutation inference for parameters and marginal effects.

Provides permutation test functions: permute the response variable, refit the model, collect null distributions of test statistics, and compute permutation p-values. Accepts containers and returns containers.

Functions:

Classes¶

Functions¶

compute_mee_permutation¶

compute_mee_permutation(spec: 'ModelSpec', bundle: 'DataBundle', fit: 'FitState', mee: 'MeeState', data: 'pl.DataFrame', conf_level: float, se_obs: np.ndarray, n_perm: int = 1000, seed: int | None = None, null: float = 0.0, alternative: str = 'two-sided', save_resamples: bool = True) -> 'tuple[MeeState, np.ndarray | None]'

Compute permutation-based inference for marginal effects.

Permutes the response variable, refits and re-explores on each permutation, and computes p-values based on the null distribution.

Parameters:

Returns:

compute_params_permutation¶

compute_params_permutation(spec: 'ModelSpec', bundle: 'DataBundle', fit: 'FitState', conf_level: float, n_perm: int = 1000, seed: int | None = None, save_resamples: bool = True, null: float = 0.0, alternative: str = 'two-sided') -> 'InferenceState'

Compute permutation-based inference for model parameters.

Permutes the response variable to break the relationship with predictors, refits the model on each permutation, and computes p-values based on the null distribution of test statistics.

Parameters:

Returns:

prediction¶

Prediction inference dispatch.

Routes prediction inference requests to the appropriate method: asymptotic (delta method), bootstrap, or cross-validation.

Functions:

Classes¶

Functions¶

dispatch_prediction_inference¶

dispatch_prediction_inference(*, how: str, pred: PredictionState, spec: ModelSpec, bundle: DataBundle, fit: FitState, conf_level: float, n_boot: int = 1000, ci_type: str = 'percentile', seed: int | None = None, k: int | str = 10, holdout_group_ids: np.ndarray | None = None) -> tuple[PredictionState, CVState | None]

Dispatch prediction inference to the appropriate method.

Parameters:

Returns:

profile¶

Profile likelihood inference for variance components.

Builds a deviance function adapter from the model’s internal state and runs profile_likelihood() to compute CIs for variance components.

compute_profile_inference: Profile likelihood CIs for mixed models. compute_profile_inference: Profile likelihood CIs for mixed models.

Functions:

Classes¶

Functions¶

compute_profile_inference¶

compute_profile_inference(spec: 'ModelSpec', bundle: 'DataBundle', fit: 'FitState', conf_level: float = 0.95, *, n_steps: int = 20, verbose: bool = False, threshold: float | None = None) -> 'ProfileState'

Compute profile likelihood CIs for variance components.

Builds a deviance function adapter from the unified model’s internal state and runs profile_likelihood(). Returns a ProfileState with profile traces, splines, and CIs on the theta and SD scales.

Parameters:

Returns:

resample¶

Internal resampling operations.

Provides bootstrap procedures for bossanova models.

Module layout:

• core.py — index generation, p-values, CI methods

• common.py — shared bootstrap execution loop, NaN handling

• results.py — result container (BootstrapResult)

• simulate.py — family-specific response simulation

• lm_operators.py — hat-matrix operator for LM

• lmer.py — LMER parametric bootstrap

• glmer.py — GLMER parametric bootstrap

Modules:

Classes:

Functions:

Attributes:

Attributes¶

BCa_MIN_NBOOT¶

BCa_MIN_NBOOT = 50

Classes¶

BootstrapResult¶

BootstrapResult(observed: Array, boot_samples: Array, ci_lower: Array, ci_upper: Array, param_names: list[str], n_boot: int, ci_type: Literal['percentile', 'basic', 'bca'], level: float, term_types: list[str] | None = None) -> None

Results from a bootstrap procedure.

Attributes:

Attributes¶

boot_samples¶

boot_samples: Array

ci¶

ci: tuple[Array, Array]

Return confidence intervals as a tuple (lower, upper).

ci_lower¶

ci_lower: Array

ci_type¶

ci_type: Literal['percentile', 'basic', 'bca']

ci_upper¶

ci_upper: Array

level¶

level: float

n_boot¶

n_boot: int

observed¶

observed: Array

param_names¶

param_names: list[str]

se¶

se: Array

Bootstrap standard errors (std of boot_samples).

term_types¶

term_types: list[str] | None = None

Functions¶

build_lm_coefficient_operator¶

build_lm_coefficient_operator(X: Array, method: Literal['cholesky', 'qr'] = 'cholesky') -> Callable[[Array], Array]

Create reusable linear operator Y -> coefficients.

Precomputes and caches the expensive factorization. Returns a function that maps Y -> coefficients.

This implements the “hat-matrix trick” for efficient permutation testing: For fixed design matrix X, coefficients are a linear map from Y -> beta_hat: beta_hat = (X’X)^-1 X’y. We precompute the factorization of (X’X) once and reuse it for all permuted/bootstrapped Y values.

Parameters:

Returns:

Examples:

X = np.array([[1, 2], [1, 3], [1, 4]])
operator = build_lm_coefficient_operator(X)
y = np.array([5, 8, 11])
coef = operator(y)

compute_batched_glm_bootstrap¶

compute_batched_glm_bootstrap(boot_indices: np.ndarray, X: np.ndarray, y: np.ndarray, family: Family, *, weights: np.ndarray | None = None, max_iter: int = 25, tol: float = 1e-08, chunk_size: int | None = None) -> np.ndarray

Compute bootstrap GLM coefficients using chunked vmap.

For GLM models (non-robust, no random effects), each bootstrap replicate is an independent IRLS solve. This function batches those solves via jax.vmap — on JAX, this compiles to a single XLA program where each lane runs its own lax.while_loop.

Memory safety is ensured by chunking: compute_batch_size determines how many solves to batch at once based on available memory. The IRLS state tuple is larger than OLS, so chunk sizes are smaller.

Parameters:

Returns:

compute_batched_ols_bootstrap¶

compute_batched_ols_bootstrap(boot_indices: np.ndarray, X: np.ndarray, y: np.ndarray, chunk_size: int | None = None) -> np.ndarray

Compute bootstrap OLS coefficients using chunked vmap.

For LM models (gaussian family, identity link, no random effects), each bootstrap replicate is an independent OLS solve. This function batches those solves via ops.vmap — on JAX, this compiles to a single XLA program; on NumPy, it degrades to a loop+stack (equivalent to sequential but without Python overhead from fit_model).

Memory safety is ensured by chunking: compute_batch_size determines how many solves to batch at once based on available memory.

Parameters:

Returns:

compute_bootstrap_ci_basic¶

compute_bootstrap_ci_basic(observed: 'ArrayLike', boot_stats: 'ArrayLike', level: float = 0.95) -> tuple[np.ndarray, np.ndarray]

Compute basic (pivotal) bootstrap confidence intervals.

The basic bootstrap CI is: [2theta_hat - q_upper, 2theta_hat - q_lower] where q_lower and q_upper are the quantiles of the bootstrap distribution.

Parameters:

Returns:

Examples:

observed = jnp.array([1.0, 2.0])
boot_stats = jnp.array([[1.1, 2.1], [0.9, 1.9], [1.2, 2.2]])
lower, upper = compute_bootstrap_ci_basic(observed, boot_stats, level=0.95)

compute_bootstrap_ci_bca¶

compute_bootstrap_ci_bca(boot_stats: 'ArrayLike', observed: 'ArrayLike', jackknife_stats: 'ArrayLike | None' = None, level: float = 0.95) -> tuple[np.ndarray, np.ndarray]

BCa (bias-corrected and accelerated) bootstrap confidence intervals.

BCa intervals adjust for bias and skewness in the bootstrap distribution, providing second-order accuracy. This is R’s gold standard method.

Parameters:

Returns:

Efron, B. (1987). Better Bootstrap Confidence Intervals. Efron, B. (1987). Better Bootstrap Confidence Intervals. DiCiccio, T. J., & Efron, B. (1996). Bootstrap confidence intervals.

compute_bootstrap_ci_percentile¶

compute_bootstrap_ci_percentile(boot_stats: 'ArrayLike', level: float = 0.95) -> tuple[np.ndarray, np.ndarray]

Compute percentile bootstrap confidence intervals.

Parameters:

Returns:

Examples:

boot_stats = jnp.array([[1.0, 2.0], [1.5, 2.5], [0.8, 1.8]])
lower, upper = compute_bootstrap_ci_percentile(boot_stats, level=0.95)

compute_pvalues¶

compute_pvalues(observed_stats: 'ArrayLike', null_stats: 'ArrayLike', alternative: Literal['two-sided', 'greater', 'less'] = 'two-sided') -> np.ndarray

Compute permutation p-values.

Uses formula: p = (1 + sum(|T_null| >= |T_obs|)) / (B + 1) where B is the number of null samples.

Parameters:

Returns:

Examples:

observed = jnp.array([2.5, -1.0])
null = jnp.array([[0.1, 0.2], [-0.3, 0.5], [0.8, -0.1]])
pvals = compute_pvalues(observed, null, alternative="two-sided")

generate_bootstrap_indices¶

generate_bootstrap_indices(rng: 'RNG', n: int, n_boot: int) -> np.ndarray

Generate n_boot bootstrap index arrays (with replacement).

Parameters:

Returns:

Examples:

from rng import RNG
rng = RNG.from_seed(42)
indices = generate_bootstrap_indices(rng, n=5, n_boot=3)
indices.shape
# (3, 5)

generate_kfold_indices¶

generate_kfold_indices(rng: 'RNG', n: int, k: int, shuffle: bool = True) -> list[tuple[np.ndarray, np.ndarray]]

Generate k-fold cross-validation train/test index splits.

Parameters:

Returns:

Examples:

from rng import RNG
rng = RNG.from_seed(42)
splits = generate_kfold_indices(rng, n=10, k=5)
len(splits)
# 5
train, test = splits[0]
len(train) + len(test)
# 10

generate_loo_indices¶

generate_loo_indices(n: int) -> list[tuple[np.ndarray, np.ndarray]]

Generate leave-one-out cross-validation train/test index splits.

Parameters:

Returns:

Examples:

splits = generate_loo_indices(n=5)
len(splits)
# 5
train, test = splits[0]
len(test)
# 1

generate_permutation_indices¶

generate_permutation_indices(rng: 'RNG', n: int, n_perm: int) -> np.ndarray

Generate n_perm permutation index arrays.

Parameters:

Returns:

Examples:

from rng import RNG
rng = RNG.from_seed(42)
indices = generate_permutation_indices(rng, n=5, n_perm=3)
indices.shape
# (3, 5)

glmer_bootstrap¶

glmer_bootstrap(X: 'Array', Z: sp.csc_matrix | 'Array', y: 'Array', X_names: list[str], theta: 'Array', beta: 'Array', family: Family, n_groups_list: list[int], re_structure: str, metadata: dict, weights: 'Array | None' = None, n_boot: int = 999, seed: int | None = None, boot_type: Literal['parametric'] = 'parametric', ci_type: Literal['percentile', 'basic', 'bca'] = 'percentile', level: float = 0.95, max_iter: int = 25, tol: float = 1e-08, max_outer_iter: int = 10000, nAGQ: int = 0, verbose: bool = False, n_jobs: int = 1) -> BootstrapResult

Bootstrap inference for GLMER coefficients.

Each bootstrap iteration requires full PIRLS optimization. Supports parallel execution via joblib.

Uses parametric bootstrap: simulate b* ~ N(0, LL’), compute eta* = X beta + Z b*, then y* from the family distribution and refit. This is the standard approach for mixed models, matching R’s lme4::bootMer(type="parametric").

Note: Only parametric bootstrap is supported for mixed models, matching R/lme4. Residual and case bootstrap are not applicable because the residual structure involves multiple variance components (within-group and between-group). Parametric bootstrap properly accounts for this hierarchical structure by simulating new random effects.

Parameters:

Returns:

Examples:

from family import binomial
family = binomial()
result = glmer_bootstrap(
    X, Z, y, ["Intercept", "x"],
    theta=theta_hat, beta=beta_hat,
    family=family, n_groups_list=[15],
    re_structure=[{"type": "intercept"}],
    metadata={}, n_boot=999
)
print(result.ci)

lmer_bootstrap¶

lmer_bootstrap(X: 'Array', Z: sp.csc_matrix | 'Array', y: 'Array', X_names: list[str], theta: 'Array', beta: 'Array', sigma: float, n_groups_list: list[int], re_structure: str, metadata: dict | None = None, lower_bounds: list[float] | None = None, n_boot: int = 999, seed: int | None = None, boot_type: Literal['parametric'] = 'parametric', ci_type: Literal['percentile', 'basic', 'bca'] = 'percentile', level: float = 0.95, method: Literal['REML', 'ML'] = 'REML', max_iter: int = 10000, tol: float = 1e-08, verbose: bool = False, which: Literal['fixef', 'ranef', 'all'] = 'fixef', group_names: list[str] | None = None, random_names: list[str] | None = None, n_jobs: int = 1) -> BootstrapResult

Bootstrap inference for lmer parameters.

Supports parallel execution via joblib for significant speedups on multi-core machines. Each bootstrap iteration requires a full PLS optimization, making parallelization highly beneficial.

Uses parametric bootstrap: simulate y* from the fitted model (y* = X beta_hat + Z b* + eps*) and refit. This is the standard approach for mixed models, matching R’s lme4::bootMer(type="parametric").

Note: Only parametric bootstrap is supported for mixed models, matching R/lme4. Residual and case bootstrap are not applicable because the residual structure involves multiple variance components (within-group and between-group). Parametric bootstrap properly accounts for this hierarchical structure by simulating new random effects.

Parameters:

Returns:

Examples:

# Parametric bootstrap for fixed effects
result = lmer_bootstrap(
    X, Z, y, ["Intercept", "Days"], theta, beta, sigma,
    n_groups_list=[18], re_structure="slope",
    n_boot=999, boot_type="parametric"
)
print(result.ci)

# Bootstrap for variance components
result = lmer_bootstrap(
    X, Z, y, ["Intercept", "Days"], theta, beta, sigma,
    n_groups_list=[18], re_structure="slope",
    n_boot=999, which="ranef",
    group_names=["Subject"], random_names=["Intercept", "Days"]
)

Bootstrap refitting can fail to converge for some samples. Convergence Bootstrap refitting can fail to converge for some samples. Convergence failures > 10% will trigger a warning. Failed samples are excluded from CI computation.

simulate_response¶

simulate_response(rng: RNG, mu: Array, family: Family, dispersion: float = 1.0, weights: Array | None = None) -> Array

Simulate response from GLM fitted values.

Dispatches to family-specific simulation functions.

Parameters:

Returns:

Modules¶

common¶

Common bootstrap execution patterns shared across model types.

This module provides the shared bootstrap execution loop (sequential/parallel), NaN-aware CI computation, and failure-handling patterns that are duplicated across lm.py, glm.py, and mixed.py.

All iteration functions should be standalone (not closures) for joblib serialization compatibility.

Functions:

Attributes:

Attributes¶

Array¶

Array = Any

Classes¶

Functions¶

collect_boot_samples¶

collect_boot_samples(boot_samples_list: list[np.ndarray], n_boot: int, context: str = 'Bootstrap') -> tuple[np.ndarray, np.ndarray, int]

Stack bootstrap samples and compute valid-sample mask.

Counts NaN failures, issues a warning if > 10% failed, and returns the full array plus a boolean mask of valid (non-NaN) rows.

Parameters:

Returns:

compute_bootstrap_ci¶

compute_bootstrap_ci(coef_obs: np.ndarray, boot_samples_valid: np.ndarray, ci_type: Literal['percentile', 'basic', 'bca'], level: float, jackknife_stats: np.ndarray | None = None) -> tuple[np.ndarray, np.ndarray]

Compute bootstrap confidence intervals using the specified method.

Dispatches to percentile, basic, or BCa CI computation from core.py.

Parameters:

Returns:

run_bootstrap_iterations¶

run_bootstrap_iterations(iteration_fn: Callable[..., np.ndarray], rngs: list[RNG], n_boot: int, iter_kwargs: dict[str, Any], n_jobs: int = 1, verbose: bool = False, desc: str = 'Bootstrap') -> list[np.ndarray]

Run bootstrap iterations sequentially or in parallel via joblib.

This is the shared execution loop used by all bootstrap functions. The iteration function must be a standalone callable (not a closure) to support joblib serialization.

Parameters:

Returns: