inference

Statistical inference utilities for bossanova.

Classes:

Functions:

Modules:

Classes¶

CellInfo¶

CellInfo(cell_labels: list[tuple[str, ...]], cell_variances: NDArray[np.float64], cell_counts: NDArray[np.int64], cell_indices: NDArray[np.int64], factor_columns: list[str]) -> None

Information about factor cells for Welch-style inference.

When using errors='unequal_var', residual variances are computed within each cell defined by the crossing of all factors in the model. This implements the Welch-James approach for factorial designs.

Attributes:

Attributes¶

cell_counts¶

cell_counts: NDArray[np.int64]

cell_indices¶

cell_indices: NDArray[np.int64]

cell_labels¶

cell_labels: list[tuple[str, ...]]

cell_variances¶

cell_variances: NDArray[np.float64]

factor_columns¶

factor_columns: list[str]

Chi2TestResult¶

Chi2TestResult(num_df: int, chi2: float, p_value: float) -> None

Result container for chi-square test.

Attributes:

Attributes¶

chi2¶

chi2: float

num_df¶

num_df: int

p_value¶

p_value: float

FTestResult¶

FTestResult(num_df: int, den_df: float, F_value: float, p_value: float) -> None

Result container for F-test.

Attributes:

Attributes¶

F_value¶

F_value: float

den_df¶

den_df: float

num_df¶

num_df: int

p_value¶

p_value: float

InferenceResult¶

InferenceResult(se: np.ndarray, statistic: np.ndarray, p_values: np.ndarray, ci_lower: np.ndarray, ci_upper: np.ndarray, df_values: list[float]) -> None

Results from coefficient inference computation.

Attributes:

Attributes¶

ci_lower¶

ci_lower: np.ndarray

ci_upper¶

ci_upper: np.ndarray

df_values¶

df_values: list[float]

p_values¶

p_values: np.ndarray

se¶

se: np.ndarray

statistic¶

statistic: np.ndarray

MixedModelProtocol¶

Bases: Protocol

Protocol for mixed-effects models compatible with profile likelihood.

This protocol defines the minimal interface required for profile_likelihood(), supporting both legacy lmer/glmer classes and the new unified model() class.

The protocol is runtime checkable, allowing isinstance() checks if needed, but the function primarily uses duck typing for flexibility.

Functions:

Functions¶

get_theta_lower_bounds¶

get_theta_lower_bounds(n_theta: int) -> list[float]

Get lower bounds for theta parameters.

TTestResult¶

TTestResult(estimate: float, se: float, t_value: float, df: float, p_value: float) -> None

Result container for t-test.

Attributes:

Attributes¶

df¶

df: float

estimate¶

estimate: float

p_value¶

p_value: float

se¶

se: float

t_value¶

t_value: float

Functions¶

adjust_pvalues¶

adjust_pvalues(pvalues: np.ndarray, method: str = 'none') -> np.ndarray

Adjust p-values for multiple comparisons.

Implements standard p-value adjustment methods matching R’s p.adjust().

Parameters:

Returns:

Examples:

>>> p = np.array([0.01, 0.04, 0.03, 0.005])
>>> adjust_pvalues(p, method="bonferroni")
array([0.04, 0.16, 0.12, 0.02])
>>> adjust_pvalues(p, method="holm")
array([0.03, 0.08, 0.06, 0.02])

build_cholesky_with_derivs¶

build_cholesky_with_derivs(theta_group: np.ndarray, n_re: int, structure: str) -> tuple[np.ndarray, list[np.ndarray]]

Build Cholesky factor L and its derivatives w.r.t. theta.

Parameters:

Returns:

compute_aic¶

compute_aic(loglik: float, k: int) -> float

Compute Akaike Information Criterion.

AIC = -2 * loglik + 2 * k. Trades off model fit against complexity, penalizing by 2 per parameter.

Parameters:

Returns:

Examples:

>>> compute_aic(-100.0, k=3)
206.0
>>> compute_aic(-50.0, k=5)
110.0

compute_akaike_weights¶

compute_akaike_weights(ic_values: np.ndarray) -> np.ndarray

Compute Akaike weights from information criterion values.

Weights quantify relative probability that each model is the best model in the set. Works with AIC or BIC values.

w_i = exp(-0.5 * delta_i) / sum(exp(-0.5 * delta_j))

where delta_i = IC_i - min(IC).

Parameters:

Returns:

Examples:

>>> weights = compute_akaike_weights(np.array([100.0, 102.0, 110.0]))
>>> np.isclose(weights.sum(), 1.0)
True
>>> weights[0] > weights[1] > weights[2]
True

compute_bic¶

compute_bic(loglik: float, k: int, n: int) -> float

Compute Bayesian Information Criterion.

BIC = -2 * loglik + k * log(n). Penalizes complexity more heavily than AIC for n >= 8, preferring simpler models.

Parameters:

Returns:

Examples:

>>> compute_bic(-100.0, k=3, n=100)
213.81...
>>> compute_bic(-100.0, k=3, n=50)
211.73...

compute_cell_info¶

compute_cell_info(residuals: NDArray[np.float64], data: pl.DataFrame, factor_columns: list[str]) -> CellInfo

Compute cell-based variance information for Welch inference.

For multi-factor models, cells are defined by the crossing of all factors. This implements the Welch-James approach for factorial designs.

Parameters:

Returns:

Examples:

import numpy as np
import polars as pl
df = pl.DataFrame({
    "group": ["A", "A", "B", "B", "B"],
    "x": [1.0, 2.0, 3.0, 4.0, 5.0],
})
residuals = np.array([0.1, -0.1, 0.2, -0.1, -0.1])
info = compute_cell_info(residuals, df, ["group"])
info.cell_labels
# [('A',), ('B',)]
info.cell_counts
# array([2, 3])

compute_chi2_test¶

compute_chi2_test(L: np.ndarray, coef: np.ndarray, vcov: np.ndarray) -> Chi2TestResult

Compute Wald chi-square test for L @ β = 0.

Used for GLMs and GLMMs where we don’t have residual df. The Wald statistic is asymptotically chi-squared distributed.

Parameters:

Returns:

Examples:

>>> L = np.array([[-1, 1, 0], [-1, 0, 1]])
>>> coef = np.array([0.5, 1.2, 0.8])
>>> vcov = np.diag([0.1, 0.15, 0.12])
>>> result = compute_chi2_test(L, coef, vcov)
>>> result.num_df
2

compute_ci¶

compute_ci(estimate: np.ndarray, se: np.ndarray, critical: float | np.ndarray, alternative: str = 'two-sided') -> tuple[np.ndarray, np.ndarray]

Compute confidence interval bounds.

Parameters:

Returns:

Examples:

>>> estimate = np.array([1.5, 2.0])
>>> se = np.array([0.2, 0.3])
>>> critical = 1.96
>>> ci_lower, ci_upper = compute_ci(estimate, se, critical)

compute_coefficient_inference¶

compute_coefficient_inference(coef: np.ndarray, vcov: np.ndarray, df: float | np.ndarray | None, conf_level: float = 0.95, null: float = 0.0, alternative: str = 'two-sided') -> InferenceResult

Compute inference statistics for regression coefficients.

Orchestrates computation of standard errors, test statistics, p-values, and confidence intervals for coefficient estimates.

Parameters:

Returns:

Examples:

>>> coef = np.array([1.5, 2.0])
>>> vcov = np.array([[0.04, 0.01], [0.01, 0.09]])
>>> result = compute_coefficient_inference(coef, vcov, df=20)
>>> result.se.shape
(2,)

compute_contrast_variance¶

compute_contrast_variance(L: np.ndarray, vcov: np.ndarray) -> np.ndarray

Compute variance-covariance of linear contrasts L @ β.

For a contrast matrix L and coefficient vcov V, computes: Var(L @ β) = L @ V @ L’

Parameters:

Returns:

Examples:

>>> L = np.array([[-1, 1, 0], [-1, 0, 1]])
>>> vcov = np.diag([1.0, 1.5, 2.0])
>>> var_L = compute_contrast_variance(L, vcov)
>>> var_L.shape
(2, 2)

compute_cooks_distance¶

compute_cooks_distance(residuals: np.ndarray, hat: np.ndarray, sigma: float, p: int) -> np.ndarray

Compute Cook’s distance for influence.

Cook’s distance measures the influence of each observation on the fitted coefficients.

Parameters:

Returns:

Note: Cook’s distance formula (matching R’s cooks.distance()): D_i = (e_i² / (p * σ²)) * (h_i / (1 - h_i)²) where e_i is the raw residual, σ is sigma, h_i is leverage.

Values > 1 are traditionally considered influential. More conservatively, D_i > 4/(n-p) is used.

Examples:

>>> D = compute_cooks_distance(residuals, hat, sigma=2.5, p=3)
>>> influential = np.where(D > 1)[0]

compute_cr_vcov¶

compute_cr_vcov(X: NDArray[np.float64], residuals: NDArray[np.float64], cluster_ids: NDArray[np.intp], XtX_inv: NDArray[np.float64], cr_type: str = 'CR0') -> NDArray[np.float64]

Compute cluster-robust covariance matrix for Gaussian mixed models.

Implements the sandwich estimator clustered by random-effects grouping variable: (X'X)^{-1} M (X'X)^{-1} where M sums outer products of per-cluster score vectors.

Matches sandwich::vcovCL() in R.

Parameters:

Returns:

compute_deviance¶

compute_deviance(loglik: float) -> float

Compute deviance from log-likelihood.

Deviance = -2 * loglik. Measures lack of fit relative to a saturated model.

Parameters:

Returns:

Examples:

>>> compute_deviance(-100.0)
200.0
>>> compute_deviance(0.0)
0.0

compute_f_pvalue¶

compute_f_pvalue(f_stat: float, df1: int, df2: float) -> float

Compute p-value from F-statistic.

Uses the upper tail of the F-distribution: P(F > f_stat | df1, df2).

Parameters:

Returns:

Examples:

>>> p = compute_f_pvalue(4.0, df1=2, df2=97)
>>> 0 < p < 0.05
True
>>> compute_f_pvalue(0.0, df1=1, df2=100)
1.0

compute_f_test¶

compute_f_test(L: np.ndarray, coef: np.ndarray, vcov: np.ndarray, df_resid: float) -> FTestResult

Compute F-test for linear hypothesis L @ β = 0.

Uses the Wald statistic divided by the number of constraints to produce an F-statistic with appropriate degrees of freedom.

Parameters:

Returns:

Note: The Wald statistic is: W = (L @ β)’ @ (L @ V @ L’)⁻¹ @ (L @ β)

The F-statistic is F = W / q, distributed as F(q, df_resid).

Examples:

>>> L = np.array([[-1, 1, 0], [-1, 0, 1]])  # 2 contrasts
>>> coef = np.array([10.0, 12.0, 15.0])
>>> vcov = np.diag([1.0, 1.5, 2.0])
>>> result = compute_f_test(L, coef, vcov, df_resid=97)
>>> result.num_df
2

compute_glm_cr_vcov¶

compute_glm_cr_vcov(X: NDArray[np.float64], residuals: NDArray[np.float64], irls_weights: NDArray[np.float64], cluster_ids: NDArray[np.intp], XtWX_inv: NDArray[np.float64], cr_type: str = 'CR0') -> NDArray[np.float64]

Compute cluster-robust covariance matrix for non-Gaussian mixed models.

Implements the weighted sandwich estimator clustered by random-effects grouping variable. Only CR0 and CR1 are supported (CR2/CR3 require leverage adjustments that are not well-defined for GLM).

Parameters:

Returns:

compute_glm_hc_vcov¶

compute_glm_hc_vcov(X: NDArray[np.float64], residuals: NDArray[np.float64], irls_weights: NDArray[np.float64], XtWX_inv: NDArray[np.float64], hc_type: str = 'HC0') -> NDArray[np.float64]

Compute heteroscedasticity-consistent covariance matrix for GLM.

(X’WX)^{-1} X’W Omega W X (X’WX)^{-1} (X’WX)^{-1} X’W Omega W X (X’WX)^{-1}

where W are the IRLS weights and Omega accounts for variance misspecification.

Parameters:

Returns:

compute_hc_vcov¶

compute_hc_vcov(X: NDArray[np.float64], residuals: NDArray[np.float64], XtX_inv: NDArray[np.float64], hc_type: str = 'HC3') -> NDArray[np.float64]

Compute heteroscedasticity-consistent covariance matrix.

Implements the sandwich estimator: (X’X)^{-1} X’ Omega X (X’X)^{-1} where Omega is a diagonal matrix of squared (possibly adjusted) residuals.

Parameters:

Returns:

compute_leverage¶

compute_leverage(X: np.ndarray, weights: np.ndarray | None = None, XtWX_inv: np.ndarray | None = None) -> np.ndarray

Compute diagonal of hat matrix (leverage values).

For the unweighted case (OLS), uses efficient QR decomposition. For the weighted case (GLM), can use pre-computed (X’WX)^{-1} to avoid redundant QR decomposition when this matrix is already available from IRLS.

Parameters:

Returns:

Note: Unweighted case: H = X(X’X)^{-1} X’, computed via QR as diag(H) = ||Q[i, :]||². Weighted case: H = W^{1/2} X (X’WX)^{-1} X’ W^{1/2}. The sum of leverage values equals p (number of predictors).

Examples:

>>> X = np.array([[1, 2], [1, 3], [1, 4]])
>>> h = compute_leverage(X)
>>> np.sum(h)  # Should equal p (number of predictors)
2.0

>>> # GLM with pre-computed XtWX_inv (avoids redundant QR)
>>> weights = np.array([1.0, 0.5, 0.8])
>>> XtWX = X.T @ np.diag(weights) @ X
>>> XtWX_inv = np.linalg.inv(XtWX)
>>> h = compute_leverage(X, weights=weights, XtWX_inv=XtWX_inv)

compute_mvt_critical¶

compute_mvt_critical(conf_level: float, corr: np.ndarray, df: float, *, tol: float = 1e-06) -> float

Compute multivariate-t critical value for simultaneous inference.

Finds c such that P(-c < T_i < c for all i) = conf_level where T ~ MVT(df, corr). Uses scipy’s Genz-Bretz CDF integration with Brent root-finding, matching R’s mvtnorm::qmvt().

For df = inf (asymptotic), uses multivariate normal instead.

Parameters:

Returns:

Examples:

>>> corr = np.array([[1.0, 0.5], [0.5, 1.0]])
>>> crit = compute_mvt_critical(0.95, corr, df=50)
>>> crit > 1.96  # wider than unadjusted
True

compute_pvalue¶

compute_pvalue(statistic: np.ndarray, df: float | np.ndarray | None = None, alternative: str = 'two-sided') -> np.ndarray

Compute p-values from test statistics.

Uses t-distribution if df is provided, otherwise uses normal distribution. Supports two-sided, greater, and less alternatives.

Parameters:

Returns:

Examples:

>>> statistic = np.array([2.5, -1.8])
>>> p = compute_pvalue(statistic, df=20)
>>> all((p >= 0) & (p <= 1))
True

compute_satterthwaite_df¶

compute_satterthwaite_df(vcov_beta: np.ndarray, jacobian_vcov: np.ndarray, hessian_deviance: np.ndarray, min_df: float = 1.0, max_df: float = 1000000.0, tol: float = 1e-08) -> np.ndarray

Compute Satterthwaite degrees of freedom for each fixed effect.

Implements the Satterthwaite approximation for denominator degrees of freedom:

df_j = 2 * (SE(β̂_j)²)² / Var(SE(β̂_j)²)

Where Var(SE(β̂_j)²) is computed via the delta method using the Jacobian of Vcov(beta) w.r.t. variance parameters and the inverse Hessian of the deviance function.

Parameters:

Returns:

Notes:

• Uses Moore-Penrose pseudo-inverse for the Hessian to handle negative or near-zero eigenvalues (convergence issues)

• Warns when negative eigenvalues are detected

• Clips df to [min_df, max_df] range for numerical stability

• Returns max_df when variance of variance is near zero

Examples:

# After fitting lmer model and computing Hessian/Jacobian
df = compute_satterthwaite_df(vcov_beta, jac, hess)
# df[0] is the denominator df for testing beta[0]

compute_satterthwaite_summary_table¶

compute_satterthwaite_summary_table(beta: np.ndarray, beta_names: list[str], vcov_beta: np.ndarray, vcov_varpar: np.ndarray, jac_list: list[np.ndarray]) -> dict[str, list]

Compute full coefficient table with Satterthwaite df and p-values.

Computes t-statistics, degrees of freedom, and p-values for all fixed effects using Satterthwaite’s approximation.

This is a convenience function that combines compute_satterthwaite_df() and compute_satterthwaite_t_test() into a single summary table.

Parameters:

Returns:

Notes:

• Computes one t-test per fixed effect

• Uses standard t-distribution for p-values

• All computations vectorized for efficiency

Examples:

# After fitting lmer and computing Hessian/Jacobian
table = compute_satterthwaite_summary_table(
    beta, beta_names, vcov_beta, vcov_varpar, jac_list
)
import pandas as pd
df = pd.DataFrame(table)
print(df)

compute_satterthwaite_t_test¶

compute_satterthwaite_t_test(beta: np.ndarray, se: np.ndarray, df: np.ndarray, conf_level: float = 0.95) -> dict[str, np.ndarray]

Compute t-statistics, p-values, and confidence intervals.

Uses Satterthwaite degrees of freedom for each coefficient to compute t-test statistics and p-values.

Parameters:

Returns:

Examples:

# After computing df from compute_satterthwaite_df()
se = np.sqrt(np.diag(vcov_beta))
results = compute_satterthwaite_t_test(beta, se, df)
print(results['p_value'])  # Two-tailed p-values

compute_sd_jacobian¶

compute_sd_jacobian(theta: np.ndarray, sigma: float, group_names: list[str], random_names: list[str] | dict[str, list[str]], re_structure: str | list[str] | dict[str, str]) -> tuple[np.ndarray, np.ndarray]

Compute SDs and Jacobian of SDs w.r.t. varpar = [theta, sigma].

Sigma = sigma^2 * L @ L.T Sigma = sigma^2 * L @ L.T SD[i] = sqrt(Sigma[i,i]) = sigma * sqrt((L @ L.T)[i,i])

The Jacobian has shape (n_sds, n_theta + 1) where the last column corresponds to sigma.

Parameters:

Returns:

compute_se_from_vcov¶

compute_se_from_vcov(vcov: np.ndarray) -> np.ndarray

Compute standard errors from variance-covariance matrix.

Parameters:

Returns:

Examples:

>>> vcov = np.array([[0.04, 0.01], [0.01, 0.09]])
>>> se = compute_se_from_vcov(vcov)
>>> np.allclose(se, [0.2, 0.3])
True

compute_sigma_se_wald¶

compute_sigma_se_wald(model: MixedModelProtocol) -> float

Compute Wald standard error for sigma.

Uses the delta method approximation from the Hessian.

Parameters:

Returns:

compute_studentized_residuals¶

compute_studentized_residuals(residuals: np.ndarray, hat: np.ndarray, sigma: float) -> np.ndarray

Compute internally studentized (standardized) residuals.

Parameters:

Returns:

Note: Formula: r_i / (σ * sqrt(1 - h_i)). Also called “internally studentized” or “standardized” residuals. They have approximately unit variance under model assumptions.

Examples:

>>> residuals = np.array([0.5, -0.3, 0.8])
>>> hat = np.array([0.2, 0.3, 0.5])
>>> sigma = 1.2
>>> r_std = compute_studentized_residuals(residuals, hat, sigma)

compute_t_critical¶

compute_t_critical(conf_int: float, df: float | np.ndarray, alternative: str = 'two-sided') -> float | np.ndarray

Compute t-distribution critical value for confidence interval.

Parameters:

Returns:

Examples:

>>> crit = compute_t_critical(0.95, df=10)
>>> abs(crit - 2.228) < 0.01
True

compute_t_test¶

compute_t_test(L: np.ndarray, coef: np.ndarray, vcov: np.ndarray, df: float) -> TTestResult

Compute t-test for a single contrast L @ β = 0.

For a single contrast (q=1), returns t-statistic instead of F.

Parameters:

Returns:

Examples:

>>> L = np.array([-1, 1, 0])  # Compare coef[1] to coef[0]
>>> coef = np.array([10.0, 12.0, 15.0])
>>> vcov = np.diag([1.0, 1.5, 2.0])
>>> result = compute_t_test(L, coef, vcov, df=97)
>>> result.estimate
2.0

compute_tukey_critical¶

compute_tukey_critical(conf_level: float, k: int, df: float) -> float

Compute Tukey HSD critical value for pairwise comparisons.

Uses the studentized range distribution to compute the critical value for simultaneous pairwise comparisons of k group means. Matches R’s qtukey(conf_level, k, df) / sqrt(2).

Parameters:

Returns:

Examples:

>>> crit = compute_tukey_critical(0.95, k=3, df=28)
>>> abs(crit - 2.474) < 0.001
True

compute_vif¶

compute_vif(X: np.ndarray, X_names: list[str]) -> pl.DataFrame

Compute variance inflation factors.

Parameters:

Returns:

Note: VIF measures how much the variance of a coefficient is inflated due to multicollinearity. VIF_j = 1 / (1 - R²_j) where R²_j is from regressing predictor j on all other predictors.

Excludes intercept column. Rules of thumb: VIF > 10: severe multicollinearity, VIF > 5: moderate multicollinearity, VIF < 5: acceptable.

Examples:

>>> X = np.array([[1, 2, 3], [1, 3, 4], [1, 4, 5]])
>>> names = ['Intercept', 'x1', 'x2']
>>> vif_df = compute_vif(X, names)

compute_wald_ci_varying¶

compute_wald_ci_varying(theta: np.ndarray, sigma: float, vcov_varpar: np.ndarray, group_names: list[str], random_names: list[str] | dict[str, list[str]], re_structure: str | list[str] | dict[str, str], conf_level: float = 0.95) -> tuple[np.ndarray, np.ndarray]

Compute Wald CIs for variance components on SD scale.

Uses the delta method to transform asymptotic variance from varpar (theta, sigma) scale to the standard deviation scale.

Parameters:

Returns:

Notes: The delta method approximation assumes:

1. Asymptotic normality of theta and sigma estimates

2. Sample size large enough for the approximation to hold

For variance components near zero (boundary), profile likelihood CIs are more accurate as they respect the parameter bounds.

compute_wald_statistic¶

compute_wald_statistic(contrast_values: np.ndarray, contrast_vcov: np.ndarray) -> float

Compute Wald statistic for testing L @ β = 0.

Computes: W = (L @ β)’ @ (L @ V @ L’)⁻¹ @ (L @ β)

Uses solve for efficiency, with pinv fallback for singular matrices.

Parameters:

Returns:

Examples:

>>> contrast_values = np.array([2.0, 5.0])
>>> contrast_vcov = np.array([[2.5, 0.5], [0.5, 3.0]])
>>> W = compute_wald_statistic(contrast_values, contrast_vcov)
>>> W > 0
True

compute_welch_satterthwaite_df_per_coef¶

compute_welch_satterthwaite_df_per_coef(X: NDArray[np.float64], cell_info: CellInfo) -> NDArray[np.float64]

Compute per-coefficient Welch-Satterthwaite degrees of freedom.

For each coefficient j, decomposes its sandwich variance into per-cell contributions and applies the Satterthwaite approximation independently. This gives genuinely different df values for different coefficients.

For a two-group treatment-coded design y ~ factor(group):

• Intercept (reference group mean): df = n_ref - 1

• Slope (group difference): classic two-group Welch formula

Algorithm for each coefficient j::

W = X @ (X'X)^{-1}                        # (n, p) projection weights
cell_W_sq[k, j] = sum_{i in cell_k} W[i,j]^2
v[k, j] = s_k^2 * cell_W_sq[k, j]        # cell k's contribution
df_j = (sum_k v[k,j])^2 / sum_k [v[k,j]^2 / (n_k - 1)]

Parameters:

Returns:

Notes: This is a bossanova extension with no direct R equivalent. R’s sandwich/emmeans/coeftest use residual df with HC SEs. Bossanova computes proper Satterthwaite df per coefficient, which for a two-group model exactly matches the classic Welch t-test (t.test(var.equal=FALSE) / scipy.stats.ttest_ind).

• For symmetric designs (equal n and variance), all df are equal

• The total variance sum_k v[k,j] equals vcov[j,j] from compute_welch_vcov

• Loops over cells (typically 2-5), vectorized over coefficients

Examples:

import numpy as np
# Two-group treatment-coded design
X = np.column_stack([np.ones(20), np.array([0]*10 + [1]*10)])
# ... with cell_info from compute_cell_info ...
df = compute_welch_satterthwaite_df_per_coef(X, cell_info)
# df[0] ≈ n_ref - 1 (intercept)
# df[1] ≈ classic Welch df (slope)

compute_z_critical¶

compute_z_critical(conf_int: float, alternative: str = 'two-sided') -> float

Compute z-distribution critical value for confidence interval.

Parameters:

Returns:

Examples:

>>> crit = compute_z_critical(0.95)
>>> abs(crit - 1.96) < 0.01
True

convert_theta_ci_to_sd¶

convert_theta_ci_to_sd(ci_theta: dict[str, tuple[float, float]], theta_opt: np.ndarray, sigma_opt: float, group_names: list[str], random_names: list[str] | dict[str, list[str]], re_structure: str | list[str] | dict[str, str]) -> tuple[dict[str, tuple[float, float]], np.ndarray, np.ndarray]

Convert theta-scale CIs to SD-scale CIs.

For variance components, SD = sigma * |theta| for simple cases. For slope models with Cholesky factors, we use the diagonal elements of sigma^2 * L @ L.T.

Parameters:

Returns:

delta_method_se¶

delta_method_se(X_pred: np.ndarray, vcov: np.ndarray) -> np.ndarray

Compute standard errors for predictions via delta method.

For predictions Xβ, the variance is: Var(Xβ) = X Var(β) X^T

Parameters:

Returns:

Note: Uses efficient computation: SE_i = sqrt(sum((X_pred @ vcov) * X_pred, axis=1)) which avoids forming the full quadratic form.

Examples:

>>> X_new = np.array([[1, 2], [1, 3]])
>>> vcov = np.array([[0.1, 0.01], [0.01, 0.05]])
>>> se = delta_method_se(X_new, vcov)

extract_ci_bound¶

extract_ci_bound(spline: UnivariateSpline, zeta_target: float, lower_bound: float, upper_guess: float) -> float

Extract CI bound by finding where spline equals target zeta.

Parameters:

Returns:

extract_factors_from_formula¶

extract_factors_from_formula(formula: str, data: pl.DataFrame) -> list[str]

Extract factor (categorical) column names from a model formula.

Parses the RHS of the formula and identifies which columns are categorical (String, Categorical, or Enum type in polars).

Parameters:

Returns:

Examples:

import polars as pl
df = pl.DataFrame({
    "y": [1.0, 2.0, 3.0],
    "group": ["A", "B", "A"],
    "x": [0.1, 0.2, 0.3],
})
extract_factors_from_formula("y ~ group + x", df)
# ['group']

# Multiple factors
df2 = df.with_columns(pl.col("group").alias("treatment"))
extract_factors_from_formula("y ~ group * treatment", df2)
# ['group', 'treatment']

Notes:

• Handles interaction terms (*, :) by extracting component columns

• Handles transforms like center(x) by extracting the inner column

• Skips intercept terms (1, 0, -1)

• Only returns columns that exist in the data

format_pvalue_with_stars¶

format_pvalue_with_stars(p_val: float) -> str

Format p-value with R-style significance codes.

Parameters:

Returns:

Examples:

>>> format_pvalue_with_stars(0.0001)
'< 0.001 *** '
>>> format_pvalue_with_stars(0.023)
'  0.023 *   '

parse_conf_int¶

parse_conf_int(conf_int: float | int | str) -> float

Parse flexible confidence interval input to float.

Parameters:

Returns:

Examples:

>>> parse_conf_int(0.95)
0.95
>>> parse_conf_int(95)
0.95
>>> parse_conf_int("95%")
0.95

profile_likelihood¶

profile_likelihood(model: MixedModelProtocol, conf_level: float = 0.95, threshold: float = 4.0, n_steps: int = 15, verbose: bool = False) -> dict

Compute profile likelihood confidence intervals for variance components.

Profiles theta parameters (relative Cholesky factors) to compute profile likelihood CIs that are more accurate than Wald CIs, especially when variance components are near the boundary (zero).

Sigma (residual SD) is derived analytically at each profile step. Sigma CI is computed via the Wald approximation.

Parameters:

Returns:

Notes:

• Profile CIs respect parameter bounds (variances >= 0)

• Computational cost: ~2 * n_steps * n_theta deviance evaluations

• For large models, consider using Wald CIs (cheaper) or bootstrap

Examples:

>>> # Legacy API
>>> m = lmer("y ~ x + (1|group)", data).fit()
>>> prof = profile_likelihood(m)
>>> prof.ci_sd  # CIs on SD scale
{'group:Intercept': (15.21, 40.23), 'Residual': (22.83, 28.69)}

>>> # New unified API
>>> m = model("y ~ x + (1|group)", data).fit()
>>> prof = profile_likelihood(m)

profile_theta_parameter¶

profile_theta_parameter(param_idx: int, param_name: str, param_opt: float, theta_opt: np.ndarray, dev_opt: float, deviance_fn: Callable[[np.ndarray], float], lower_bounds: np.ndarray, n_steps: int, threshold: float, verbose: bool) -> dict

Profile a single theta parameter bidirectionally from its MLE.

At each step, fixes the focal theta and re-optimizes all other theta elements to find the conditional MLE deviance.

Parameters:

Returns:

satterthwaite_df_for_contrasts¶

satterthwaite_df_for_contrasts(L: np.ndarray, vcov_beta: np.ndarray, vcov_varpar: np.ndarray, jac_list: list[np.ndarray], min_df: float = 1.0, max_df: float = 1000000.0) -> np.ndarray

Compute Satterthwaite degrees of freedom for arbitrary contrasts.

This generalizes compute_satterthwaite_df() to handle arbitrary linear contrasts L @ β, not just individual coefficients. This is needed for emmeans where each EMM is a linear combination of coefficients.

The Satterthwaite formula for contrast L is:

df = 2 * var_contrast² / Var(var_contrast)

var_contrast = L @ Vcov(β) @ L.T var_contrast = L @ Vcov(β) @ L.T grad[i] = L @ Jac_i @ L.T (gradient w.r.t. variance parameter i) Var(var_contrast) = grad @ Vcov(varpar) @ grad

Parameters:

Returns:

Examples:

# Compute df for EMMs where X_ref is the prediction matrix
df = satterthwaite_df_for_contrasts(
    X_ref, vcov_beta, vcov_varpar, jac_list
)
# df[i] is the denominator df for EMM i

Modules¶

diagnostics¶

Residuals, influence measures, and leverage statistics.

Functions:

Attributes¶

Classes¶

Functions¶

compute_cooks_distance¶

compute_cooks_distance(residuals: np.ndarray, hat: np.ndarray, sigma: float, p: int) -> np.ndarray

Compute Cook’s distance for influence.

Cook’s distance measures the influence of each observation on the fitted coefficients.

Parameters:

Returns:

Note: Cook’s distance formula (matching R’s cooks.distance()): D_i = (e_i² / (p * σ²)) * (h_i / (1 - h_i)²) where e_i is the raw residual, σ is sigma, h_i is leverage.

Values > 1 are traditionally considered influential. More conservatively, D_i > 4/(n-p) is used.

Examples:

>>> D = compute_cooks_distance(residuals, hat, sigma=2.5, p=3)
>>> influential = np.where(D > 1)[0]

compute_leverage¶

compute_leverage(X: np.ndarray, weights: np.ndarray | None = None, XtWX_inv: np.ndarray | None = None) -> np.ndarray

Compute diagonal of hat matrix (leverage values).

For the unweighted case (OLS), uses efficient QR decomposition. For the weighted case (GLM), can use pre-computed (X’WX)^{-1} to avoid redundant QR decomposition when this matrix is already available from IRLS.

Parameters:

Returns:

Note: Unweighted case: H = X(X’X)^{-1} X’, computed via QR as diag(H) = ||Q[i, :]||². Weighted case: H = W^{1/2} X (X’WX)^{-1} X’ W^{1/2}. The sum of leverage values equals p (number of predictors).

Examples:

>>> X = np.array([[1, 2], [1, 3], [1, 4]])
>>> h = compute_leverage(X)
>>> np.sum(h)  # Should equal p (number of predictors)
2.0

>>> # GLM with pre-computed XtWX_inv (avoids redundant QR)
>>> weights = np.array([1.0, 0.5, 0.8])
>>> XtWX = X.T @ np.diag(weights) @ X
>>> XtWX_inv = np.linalg.inv(XtWX)
>>> h = compute_leverage(X, weights=weights, XtWX_inv=XtWX_inv)

compute_studentized_residuals¶

compute_studentized_residuals(residuals: np.ndarray, hat: np.ndarray, sigma: float) -> np.ndarray

Compute internally studentized (standardized) residuals.

Parameters:

Returns:

Note: Formula: r_i / (σ * sqrt(1 - h_i)). Also called “internally studentized” or “standardized” residuals. They have approximately unit variance under model assumptions.

Examples:

>>> residuals = np.array([0.5, -0.3, 0.8])
>>> hat = np.array([0.2, 0.3, 0.5])
>>> sigma = 1.2
>>> r_std = compute_studentized_residuals(residuals, hat, sigma)

compute_vif¶

compute_vif(X: np.ndarray, X_names: list[str]) -> pl.DataFrame

Compute variance inflation factors.

Parameters:

Returns:

Note: VIF measures how much the variance of a coefficient is inflated due to multicollinearity. VIF_j = 1 / (1 - R²_j) where R²_j is from regressing predictor j on all other predictors.

Excludes intercept column. Rules of thumb: VIF > 10: severe multicollinearity, VIF > 5: moderate multicollinearity, VIF < 5: acceptable.

Examples:

>>> X = np.array([[1, 2, 3], [1, 3, 4], [1, 4, 5]])
>>> names = ['Intercept', 'x1', 'x2']
>>> vif_df = compute_vif(X, names)

estimation¶

Standard errors, confidence intervals, and coefficient inference.

Classes:

Functions:

Classes¶

InferenceResult¶

InferenceResult(se: np.ndarray, statistic: np.ndarray, p_values: np.ndarray, ci_lower: np.ndarray, ci_upper: np.ndarray, df_values: list[float]) -> None

Results from coefficient inference computation.

Attributes:

Attributes¶

ci_lower¶

ci_lower: np.ndarray

ci_upper¶

ci_upper: np.ndarray

df_values¶

df_values: list[float]

p_values¶

p_values: np.ndarray

se¶

se: np.ndarray

statistic¶

statistic: np.ndarray

Functions¶

compute_ci¶

compute_ci(estimate: np.ndarray, se: np.ndarray, critical: float | np.ndarray, alternative: str = 'two-sided') -> tuple[np.ndarray, np.ndarray]

Compute confidence interval bounds.

Parameters:

Returns:

Examples:

>>> estimate = np.array([1.5, 2.0])
>>> se = np.array([0.2, 0.3])
>>> critical = 1.96
>>> ci_lower, ci_upper = compute_ci(estimate, se, critical)

compute_coefficient_inference¶

compute_coefficient_inference(coef: np.ndarray, vcov: np.ndarray, df: float | np.ndarray | None, conf_level: float = 0.95, null: float = 0.0, alternative: str = 'two-sided') -> InferenceResult

Compute inference statistics for regression coefficients.

Orchestrates computation of standard errors, test statistics, p-values, and confidence intervals for coefficient estimates.

Parameters:

Returns:

Examples:

>>> coef = np.array([1.5, 2.0])
>>> vcov = np.array([[0.04, 0.01], [0.01, 0.09]])
>>> result = compute_coefficient_inference(coef, vcov, df=20)
>>> result.se.shape
(2,)

compute_se_from_vcov¶

compute_se_from_vcov(vcov: np.ndarray) -> np.ndarray

Compute standard errors from variance-covariance matrix.

Parameters:

Returns:

Examples:

>>> vcov = np.array([[0.04, 0.01], [0.01, 0.09]])
>>> se = compute_se_from_vcov(vcov)
>>> np.allclose(se, [0.2, 0.3])
True

compute_t_critical¶

compute_t_critical(conf_int: float, df: float | np.ndarray, alternative: str = 'two-sided') -> float | np.ndarray

Compute t-distribution critical value for confidence interval.

Parameters:

Returns:

Examples:

>>> crit = compute_t_critical(0.95, df=10)
>>> abs(crit - 2.228) < 0.01
True

compute_z_critical¶

compute_z_critical(conf_int: float, alternative: str = 'two-sided') -> float

Compute z-distribution critical value for confidence interval.

Parameters:

Returns:

Examples:

>>> crit = compute_z_critical(0.95)
>>> abs(crit - 1.96) < 0.01
True

delta_method_se¶

delta_method_se(X_pred: np.ndarray, vcov: np.ndarray) -> np.ndarray

Compute standard errors for predictions via delta method.

For predictions Xβ, the variance is: Var(Xβ) = X Var(β) X^T

Parameters:

Returns:

Note: Uses efficient computation: SE_i = sqrt(sum((X_pred @ vcov) * X_pred, axis=1)) which avoids forming the full quadratic form.

Examples:

>>> X_new = np.array([[1, 2], [1, 3]])
>>> vcov = np.array([[0.1, 0.01], [0.01, 0.05]])
>>> se = delta_method_se(X_new, vcov)

parse_conf_int¶

parse_conf_int(conf_int: float | int | str) -> float

Parse flexible confidence interval input to float.

Parameters:

Returns:

Examples:

>>> parse_conf_int(0.95)
0.95
>>> parse_conf_int(95)
0.95
>>> parse_conf_int("95%")
0.95

hypothesis¶

Statistical hypothesis tests (F, chi-square, t) for linear contrasts.

Classes:

Functions:

Classes¶

Chi2TestResult¶

Chi2TestResult(num_df: int, chi2: float, p_value: float) -> None

Result container for chi-square test.

Attributes:

Attributes¶

chi2¶

chi2: float

num_df¶

num_df: int

p_value¶

p_value: float

FTestResult¶

FTestResult(num_df: int, den_df: float, F_value: float, p_value: float) -> None

Result container for F-test.

Attributes:

Attributes¶

F_value¶

F_value: float

den_df¶

den_df: float

num_df¶

num_df: int

p_value¶

p_value: float

TTestResult¶

TTestResult(estimate: float, se: float, t_value: float, df: float, p_value: float) -> None

Result container for t-test.

Attributes:

Attributes¶

df¶

df: float

estimate¶

estimate: float

p_value¶

p_value: float

se¶

se: float

t_value¶

t_value: float

Functions¶

adjust_pvalues¶

adjust_pvalues(pvalues: np.ndarray, method: str = 'none') -> np.ndarray

Adjust p-values for multiple comparisons.

Implements standard p-value adjustment methods matching R’s p.adjust().

Parameters:

Returns:

Examples:

>>> p = np.array([0.01, 0.04, 0.03, 0.005])
>>> adjust_pvalues(p, method="bonferroni")
array([0.04, 0.16, 0.12, 0.02])
>>> adjust_pvalues(p, method="holm")
array([0.03, 0.08, 0.06, 0.02])

compute_chi2_test¶

compute_chi2_test(L: np.ndarray, coef: np.ndarray, vcov: np.ndarray) -> Chi2TestResult

Compute Wald chi-square test for L @ β = 0.

Used for GLMs and GLMMs where we don’t have residual df. The Wald statistic is asymptotically chi-squared distributed.

Parameters:

Returns:

Examples:

>>> L = np.array([[-1, 1, 0], [-1, 0, 1]])
>>> coef = np.array([0.5, 1.2, 0.8])
>>> vcov = np.diag([0.1, 0.15, 0.12])
>>> result = compute_chi2_test(L, coef, vcov)
>>> result.num_df
2

compute_contrast_variance¶

compute_contrast_variance(L: np.ndarray, vcov: np.ndarray) -> np.ndarray

Compute variance-covariance of linear contrasts L @ β.

For a contrast matrix L and coefficient vcov V, computes: Var(L @ β) = L @ V @ L’

Parameters:

Returns:

Examples:

>>> L = np.array([[-1, 1, 0], [-1, 0, 1]])
>>> vcov = np.diag([1.0, 1.5, 2.0])
>>> var_L = compute_contrast_variance(L, vcov)
>>> var_L.shape
(2, 2)

compute_f_pvalue¶

compute_f_pvalue(f_stat: float, df1: int, df2: float) -> float

Compute p-value from F-statistic.

Uses the upper tail of the F-distribution: P(F > f_stat | df1, df2).

Parameters:

Returns:

Examples:

>>> p = compute_f_pvalue(4.0, df1=2, df2=97)
>>> 0 < p < 0.05
True
>>> compute_f_pvalue(0.0, df1=1, df2=100)
1.0

compute_f_test¶

compute_f_test(L: np.ndarray, coef: np.ndarray, vcov: np.ndarray, df_resid: float) -> FTestResult

Compute F-test for linear hypothesis L @ β = 0.

Uses the Wald statistic divided by the number of constraints to produce an F-statistic with appropriate degrees of freedom.

Parameters:

Returns:

Note: The Wald statistic is: W = (L @ β)’ @ (L @ V @ L’)⁻¹ @ (L @ β)

The F-statistic is F = W / q, distributed as F(q, df_resid).

Examples:

>>> L = np.array([[-1, 1, 0], [-1, 0, 1]])  # 2 contrasts
>>> coef = np.array([10.0, 12.0, 15.0])
>>> vcov = np.diag([1.0, 1.5, 2.0])
>>> result = compute_f_test(L, coef, vcov, df_resid=97)
>>> result.num_df
2

compute_pvalue¶

compute_pvalue(statistic: np.ndarray, df: float | np.ndarray | None = None, alternative: str = 'two-sided') -> np.ndarray

Compute p-values from test statistics.

Uses t-distribution if df is provided, otherwise uses normal distribution. Supports two-sided, greater, and less alternatives.

Parameters:

Returns:

Examples:

>>> statistic = np.array([2.5, -1.8])
>>> p = compute_pvalue(statistic, df=20)
>>> all((p >= 0) & (p <= 1))
True

compute_t_test¶

compute_t_test(L: np.ndarray, coef: np.ndarray, vcov: np.ndarray, df: float) -> TTestResult

Compute t-test for a single contrast L @ β = 0.

For a single contrast (q=1), returns t-statistic instead of F.

Parameters:

Returns:

Examples:

>>> L = np.array([-1, 1, 0])  # Compare coef[1] to coef[0]
>>> coef = np.array([10.0, 12.0, 15.0])
>>> vcov = np.diag([1.0, 1.5, 2.0])
>>> result = compute_t_test(L, coef, vcov, df=97)
>>> result.estimate
2.0

compute_wald_statistic¶

compute_wald_statistic(contrast_values: np.ndarray, contrast_vcov: np.ndarray) -> float

Compute Wald statistic for testing L @ β = 0.

Computes: W = (L @ β)’ @ (L @ V @ L’)⁻¹ @ (L @ β)

Uses solve for efficiency, with pinv fallback for singular matrices.

Parameters:

Returns:

Examples:

>>> contrast_values = np.array([2.0, 5.0])
>>> contrast_vcov = np.array([[2.5, 0.5], [0.5, 3.0]])
>>> W = compute_wald_statistic(contrast_values, contrast_vcov)
>>> W > 0
True

format_pvalue_with_stars¶

format_pvalue_with_stars(p_val: float) -> str

Format p-value with R-style significance codes.

Parameters:

Returns:

Examples:

>>> format_pvalue_with_stars(0.0001)
'< 0.001 *** '
>>> format_pvalue_with_stars(0.023)
'  0.023 *   '

information_criteria¶

Deviance, AIC, BIC, and Akaike weights.

Functions:

Functions¶

compute_aic¶

compute_aic(loglik: float, k: int) -> float

Compute Akaike Information Criterion.

AIC = -2 * loglik + 2 * k. Trades off model fit against complexity, penalizing by 2 per parameter.

Parameters:

Returns:

Examples:

>>> compute_aic(-100.0, k=3)
206.0
>>> compute_aic(-50.0, k=5)
110.0

compute_akaike_weights¶

compute_akaike_weights(ic_values: np.ndarray) -> np.ndarray

Compute Akaike weights from information criterion values.

Weights quantify relative probability that each model is the best model in the set. Works with AIC or BIC values.

w_i = exp(-0.5 * delta_i) / sum(exp(-0.5 * delta_j))

where delta_i = IC_i - min(IC).

Parameters:

Returns:

Examples:

>>> weights = compute_akaike_weights(np.array([100.0, 102.0, 110.0]))
>>> np.isclose(weights.sum(), 1.0)
True
>>> weights[0] > weights[1] > weights[2]
True

compute_bic¶

compute_bic(loglik: float, k: int, n: int) -> float

Compute Bayesian Information Criterion.

BIC = -2 * loglik + k * log(n). Penalizes complexity more heavily than AIC for n >= 8, preferring simpler models.

Parameters:

Returns:

Examples:

>>> compute_bic(-100.0, k=3, n=100)
213.81...
>>> compute_bic(-100.0, k=3, n=50)
211.73...

compute_deviance¶

compute_deviance(loglik: float) -> float

Compute deviance from log-likelihood.

Deviance = -2 * loglik. Measures lack of fit relative to a saturated model.

Parameters:

Returns:

Examples:

>>> compute_deviance(-100.0)
200.0
>>> compute_deviance(0.0)
0.0

multiplicity¶

Multiplicity adjustment for contrast confidence intervals.

Provides critical value computation for multiplicity-adjusted CIs, matching R’s emmeans defaults:

• Tukey HSD for pairwise contrasts (studentized range).

• Multivariate-t (MVT) for sequential, treatment, sum, helmert contrasts (Genz-Bretz CDF + root-finding).

• No adjustment for polynomial contrasts.

compute_tukey_critical: Tukey HSD critical value via studentized range. compute_tukey_critical: Tukey HSD critical value via studentized range. compute_mvt_critical: MVT critical value via CDF inversion.

Functions:

Functions¶

compute_mvt_critical¶

compute_mvt_critical(conf_level: float, corr: np.ndarray, df: float, *, tol: float = 1e-06) -> float

Compute multivariate-t critical value for simultaneous inference.

Finds c such that P(-c < T_i < c for all i) = conf_level where T ~ MVT(df, corr). Uses scipy’s Genz-Bretz CDF integration with Brent root-finding, matching R’s mvtnorm::qmvt().

For df = inf (asymptotic), uses multivariate normal instead.

Parameters:

Returns:

Examples:

>>> corr = np.array([[1.0, 0.5], [0.5, 1.0]])
>>> crit = compute_mvt_critical(0.95, corr, df=50)
>>> crit > 1.96  # wider than unadjusted
True

compute_tukey_critical¶

compute_tukey_critical(conf_level: float, k: int, df: float) -> float

Compute Tukey HSD critical value for pairwise comparisons.

Uses the studentized range distribution to compute the critical value for simultaneous pairwise comparisons of k group means. Matches R’s qtukey(conf_level, k, df) / sqrt(2).

Parameters:

Returns:

Examples:

>>> crit = compute_tukey_critical(0.95, k=3, df=28)
>>> abs(crit - 2.474) < 0.001
True

profile¶

Profile likelihood for variance components in mixed-effects models.

This module implements the profile likelihood algorithm for computing confidence intervals for variance components (theta parameters and sigma) in linear mixed-effects models.

The algorithm follows MixedModels.jl and lme4:

1. For each theta parameter, step bidirectionally from the MLE

2. At each step, fix that theta and re-optimize all other thetas

3. Sigma is derived analytically at each step (profile concentrated likelihood)

4. Compute the signed likelihood root: zeta = sign * sqrt(dev - dev_opt)

5. Fit B-splines for interpolation

6. Extract CIs where |zeta| = sqrt(chi2(conf_level, df=1))

• Bates et al. (2015). Fitting Linear Mixed-Effects Models Using lme4.

• Bates et al. (2015). Fitting Linear Mixed-Effects Models Using lme4.

• MixedModels.jl: https://github.com/JuliaStats/MixedModels.jl

Classes:

Functions:

Classes¶

MixedModelProtocol¶

Bases: Protocol

Protocol for mixed-effects models compatible with profile likelihood.

This protocol defines the minimal interface required for profile_likelihood(), supporting both legacy lmer/glmer classes and the new unified model() class.

The protocol is runtime checkable, allowing isinstance() checks if needed, but the function primarily uses duck typing for flexibility.

Functions:

Functions¶

get_theta_lower_bounds¶

get_theta_lower_bounds(n_theta: int) -> list[float]

Get lower bounds for theta parameters.

Functions¶

compute_sigma_se_wald¶

compute_sigma_se_wald(model: MixedModelProtocol) -> float

Compute Wald standard error for sigma.

Uses the delta method approximation from the Hessian.

Parameters:

Returns:

convert_theta_ci_to_sd¶

convert_theta_ci_to_sd(ci_theta: dict[str, tuple[float, float]], theta_opt: np.ndarray, sigma_opt: float, group_names: list[str], random_names: list[str] | dict[str, list[str]], re_structure: str | list[str] | dict[str, str]) -> tuple[dict[str, tuple[float, float]], np.ndarray, np.ndarray]

Convert theta-scale CIs to SD-scale CIs.

For variance components, SD = sigma * |theta| for simple cases. For slope models with Cholesky factors, we use the diagonal elements of sigma^2 * L @ L.T.

Parameters:

Returns:

extract_ci_bound¶

extract_ci_bound(spline: UnivariateSpline, zeta_target: float, lower_bound: float, upper_guess: float) -> float

Extract CI bound by finding where spline equals target zeta.

Parameters:

Returns:

profile_likelihood¶

profile_likelihood(model: MixedModelProtocol, conf_level: float = 0.95, threshold: float = 4.0, n_steps: int = 15, verbose: bool = False) -> dict

Compute profile likelihood confidence intervals for variance components.

Profiles theta parameters (relative Cholesky factors) to compute profile likelihood CIs that are more accurate than Wald CIs, especially when variance components are near the boundary (zero).

Sigma (residual SD) is derived analytically at each profile step. Sigma CI is computed via the Wald approximation.

Parameters:

Returns:

Notes:

• Profile CIs respect parameter bounds (variances >= 0)

• Computational cost: ~2 * n_steps * n_theta deviance evaluations

• For large models, consider using Wald CIs (cheaper) or bootstrap

Examples:

>>> # Legacy API
>>> m = lmer("y ~ x + (1|group)", data).fit()
>>> prof = profile_likelihood(m)
>>> prof.ci_sd  # CIs on SD scale
{'group:Intercept': (15.21, 40.23), 'Residual': (22.83, 28.69)}

>>> # New unified API
>>> m = model("y ~ x + (1|group)", data).fit()
>>> prof = profile_likelihood(m)

profile_theta_parameter¶

profile_theta_parameter(param_idx: int, param_name: str, param_opt: float, theta_opt: np.ndarray, dev_opt: float, deviance_fn: Callable[[np.ndarray], float], lower_bounds: np.ndarray, n_steps: int, threshold: float, verbose: bool) -> dict

Profile a single theta parameter bidirectionally from its MLE.

At each step, fixes the focal theta and re-optimizes all other theta elements to find the conditional MLE deviance.

Parameters:

Returns:

sandwich¶

Sandwich covariance matrix estimators for robust standard errors.

Provides:

• HC0-HC3: Heteroscedasticity-consistent estimators for OLS/GLM.

• CR0-CR3: Cluster-robust estimators for mixed models, clustering by the random-effects grouping variable.

Functions:

Functions¶

compute_cr_vcov¶

compute_cr_vcov(X: NDArray[np.float64], residuals: NDArray[np.float64], cluster_ids: NDArray[np.intp], XtX_inv: NDArray[np.float64], cr_type: str = 'CR0') -> NDArray[np.float64]

Compute cluster-robust covariance matrix for Gaussian mixed models.

Implements the sandwich estimator clustered by random-effects grouping variable: (X'X)^{-1} M (X'X)^{-1} where M sums outer products of per-cluster score vectors.

Matches sandwich::vcovCL() in R.

Parameters:

Returns:

compute_glm_cr_vcov¶

compute_glm_cr_vcov(X: NDArray[np.float64], residuals: NDArray[np.float64], irls_weights: NDArray[np.float64], cluster_ids: NDArray[np.intp], XtWX_inv: NDArray[np.float64], cr_type: str = 'CR0') -> NDArray[np.float64]

Compute cluster-robust covariance matrix for non-Gaussian mixed models.

Implements the weighted sandwich estimator clustered by random-effects grouping variable. Only CR0 and CR1 are supported (CR2/CR3 require leverage adjustments that are not well-defined for GLM).

Parameters:

Returns:

compute_glm_hc_vcov¶

compute_glm_hc_vcov(X: NDArray[np.float64], residuals: NDArray[np.float64], irls_weights: NDArray[np.float64], XtWX_inv: NDArray[np.float64], hc_type: str = 'HC0') -> NDArray[np.float64]

Compute heteroscedasticity-consistent covariance matrix for GLM.

(X’WX)^{-1} X’W Omega W X (X’WX)^{-1} (X’WX)^{-1} X’W Omega W X (X’WX)^{-1}

where W are the IRLS weights and Omega accounts for variance misspecification.

Parameters:

Returns:

compute_hc_vcov¶

compute_hc_vcov(X: NDArray[np.float64], residuals: NDArray[np.float64], XtX_inv: NDArray[np.float64], hc_type: str = 'HC3') -> NDArray[np.float64]

Compute heteroscedasticity-consistent covariance matrix.

Implements the sandwich estimator: (X’X)^{-1} X’ Omega X (X’X)^{-1} where Omega is a diagonal matrix of squared (possibly adjusted) residuals.

Parameters:

Returns:

satterthwaite¶

Satterthwaite degrees of freedom approximation for linear mixed models.

Functions:

Functions¶

compute_satterthwaite_df¶

compute_satterthwaite_df(vcov_beta: np.ndarray, jacobian_vcov: np.ndarray, hessian_deviance: np.ndarray, min_df: float = 1.0, max_df: float = 1000000.0, tol: float = 1e-08) -> np.ndarray

Compute Satterthwaite degrees of freedom for each fixed effect.

Implements the Satterthwaite approximation for denominator degrees of freedom:

df_j = 2 * (SE(β̂_j)²)² / Var(SE(β̂_j)²)

Where Var(SE(β̂_j)²) is computed via the delta method using the Jacobian of Vcov(beta) w.r.t. variance parameters and the inverse Hessian of the deviance function.

Parameters:

Returns:

Notes:

• Uses Moore-Penrose pseudo-inverse for the Hessian to handle negative or near-zero eigenvalues (convergence issues)

• Warns when negative eigenvalues are detected

• Clips df to [min_df, max_df] range for numerical stability

• Returns max_df when variance of variance is near zero

Examples:

# After fitting lmer model and computing Hessian/Jacobian
df = compute_satterthwaite_df(vcov_beta, jac, hess)
# df[0] is the denominator df for testing beta[0]

compute_satterthwaite_summary_table¶

compute_satterthwaite_summary_table(beta: np.ndarray, beta_names: list[str], vcov_beta: np.ndarray, vcov_varpar: np.ndarray, jac_list: list[np.ndarray]) -> dict[str, list]

Compute full coefficient table with Satterthwaite df and p-values.

Computes t-statistics, degrees of freedom, and p-values for all fixed effects using Satterthwaite’s approximation.

This is a convenience function that combines compute_satterthwaite_df() and compute_satterthwaite_t_test() into a single summary table.

Parameters:

Returns:

Notes:

• Computes one t-test per fixed effect

• Uses standard t-distribution for p-values

• All computations vectorized for efficiency

Examples:

# After fitting lmer and computing Hessian/Jacobian
table = compute_satterthwaite_summary_table(
    beta, beta_names, vcov_beta, vcov_varpar, jac_list
)
import pandas as pd
df = pd.DataFrame(table)
print(df)

compute_satterthwaite_t_test¶

compute_satterthwaite_t_test(beta: np.ndarray, se: np.ndarray, df: np.ndarray, conf_level: float = 0.95) -> dict[str, np.ndarray]

Compute t-statistics, p-values, and confidence intervals.

Uses Satterthwaite degrees of freedom for each coefficient to compute t-test statistics and p-values.

Parameters:

Returns:

Examples:

# After computing df from compute_satterthwaite_df()
se = np.sqrt(np.diag(vcov_beta))
results = compute_satterthwaite_t_test(beta, se, df)
print(results['p_value'])  # Two-tailed p-values

satterthwaite_df_for_contrasts¶

satterthwaite_df_for_contrasts(L: np.ndarray, vcov_beta: np.ndarray, vcov_varpar: np.ndarray, jac_list: list[np.ndarray], min_df: float = 1.0, max_df: float = 1000000.0) -> np.ndarray

Compute Satterthwaite degrees of freedom for arbitrary contrasts.

This generalizes compute_satterthwaite_df() to handle arbitrary linear contrasts L @ β, not just individual coefficients. This is needed for emmeans where each EMM is a linear combination of coefficients.

The Satterthwaite formula for contrast L is:

df = 2 * var_contrast² / Var(var_contrast)

var_contrast = L @ Vcov(β) @ L.T var_contrast = L @ Vcov(β) @ L.T grad[i] = L @ Jac_i @ L.T (gradient w.r.t. variance parameter i) Var(var_contrast) = grad @ Vcov(varpar) @ grad

Parameters:

Returns:

Examples:

# Compute df for EMMs where X_ref is the prediction matrix
df = satterthwaite_df_for_contrasts(
    X_ref, vcov_beta, vcov_varpar, jac_list
)
# df[i] is the denominator df for EMM i

wald_variance¶

Wald confidence intervals for variance components using delta method.

This module computes asymptotic (Wald) confidence intervals for variance components in mixed-effects models. It uses the delta method to transform the asymptotic covariance of variance parameters (theta, sigma) to the standard deviation scale.

The key insight is that the Hessian of the deviance w.r.t. variance parameters is already computed for Satterthwaite inference, and we can reuse it here: vcov_varpar = 2 * H^{-1}

For transformation from theta scale to SD scale, we use: SD[i] = sigma * sqrt(diag(L @ L.T)[i])

where L is the Cholesky factor reconstructed from theta.

Functions:

Functions¶

build_cholesky_with_derivs¶

build_cholesky_with_derivs(theta_group: np.ndarray, n_re: int, structure: str) -> tuple[np.ndarray, list[np.ndarray]]

Build Cholesky factor L and its derivatives w.r.t. theta.

Parameters:

Returns:

compute_sd_jacobian¶

compute_sd_jacobian(theta: np.ndarray, sigma: float, group_names: list[str], random_names: list[str] | dict[str, list[str]], re_structure: str | list[str] | dict[str, str]) -> tuple[np.ndarray, np.ndarray]

Compute SDs and Jacobian of SDs w.r.t. varpar = [theta, sigma].

Sigma = sigma^2 * L @ L.T Sigma = sigma^2 * L @ L.T SD[i] = sqrt(Sigma[i,i]) = sigma * sqrt((L @ L.T)[i,i])

The Jacobian has shape (n_sds, n_theta + 1) where the last column corresponds to sigma.

Parameters:

Returns:

compute_wald_ci_varying¶

compute_wald_ci_varying(theta: np.ndarray, sigma: float, vcov_varpar: np.ndarray, group_names: list[str], random_names: list[str] | dict[str, list[str]], re_structure: str | list[str] | dict[str, str], conf_level: float = 0.95) -> tuple[np.ndarray, np.ndarray]

Compute Wald CIs for variance components on SD scale.

Uses the delta method to transform asymptotic variance from varpar (theta, sigma) scale to the standard deviation scale.

Parameters:

Returns:

Notes: The delta method approximation assumes:

1. Asymptotic normality of theta and sigma estimates

2. Sample size large enough for the approximation to hold

For variance components near zero (boundary), profile likelihood CIs are more accurate as they respect the parameter bounds.

welch¶

Welch-Satterthwaite degrees of freedom for unequal variance inference.

Classes:

Functions: