Inference

Hypothesis Tests

aurora.inference.hypothesis.wald.wald_test(result, contrast, *, value=0.0, include_intercept=True)[source]

Perform a Wald test for (possibly multi-constraint) linear hypotheses.

Parameters:

  • result (GLMResult) – Fitted GLM result providing coefficient estimates and covariance.

  • contrast (Iterable[float] | Sequence[Sequence[float]] | ndarray) – Coefficient weights defining the hypothesis C beta = value. Accepts one-dimensional vectors (single constraint) or matrices with one row per constraint.

  • value (float | Iterable[float] | ndarray) – Hypothesised values for the linear combination(s). Can be a scalar or an iterable matching the number of constraints.

  • include_intercept (bool) – Whether the contrast coefficients include the model intercept as the first element when present.

Returns:

Mapping with keys statistic, p_value and df.

Return type:

dict

Confidence Intervals

aurora.inference.intervals.confidence.confidence_intervals(result, *, level=0.95, method='wald', include_intercept=True)[source]

Compute confidence intervals for the coefficients of a fitted GLM.

Parameters:
Return type:

ConfidenceIntervalResult

Diagnostics

aurora.inference.diagnostics.glm.glm_diagnostics(result)[source]

Compute residual and influence diagnostics for a fitted GLM.

Parameters:

result (GLMResult)

Return type:

GLMDiagnosticResult

Robust Covariance

aurora.inference.robust.robust_covariance(result, *, hc_type='HC3')[source]

Compute heteroscedasticity-consistent (sandwich) covariance matrix.

Computes robust standard errors that are valid under heteroscedasticity and potential model misspecification. The sandwich estimator has the form:

Var(β) = (X’WX)⁻¹ M (X’WX)⁻¹

where W is the diagonal of final IRLS weights and M is the meat, built from Pearson residuals with an HC-type correction.

Parameters:

  • result (GLMResult) – Fitted GLM result object

  • hc_type ({'HC0', 'HC1', 'HC2', 'HC3', 'HC4'}, default='HC3') –

    Type of heteroscedasticity-consistent standard errors:

    • ’HC0’: White (1980) - basic sandwich estimator, no correction

    • ’HC1’: Degrees of freedom correction (n/(n-p))

    • ’HC2’: Leverage correction (1/(1-h_i))

    • ’HC3’: MacKinnon & White (1985) - 1/(1-h_i)², better for small samples

    • ’HC4’: Cribari-Neto (2004) - best for influential points

    HC3 is recommended for general use.

Returns:

Container with robust standard errors and covariance matrix

Return type:

RobustInferenceResult

Raises:

  • RuntimeError – If design matrix or response are not available in result

  • ValueError – If hc_type is not recognized

Notes

The GLM sandwich (Zeileis 2006) uses, unlike the plain OLS version:

  • bread = (XᵀWX)⁻¹ with the final IRLS weights W = diag(prior_wᵢ / [g’(μᵢ)² V(μᵢ)]). For Gaussian/identity models this reduces to (XᵀX)⁻¹. The dispersion estimate φ̂ cancels between bread and meat, so it is not applied explicitly (this matches statsmodels, whose GLM robust covariances are scale-free).

  • meat = Xᵀ diag(Wᵢ · rᵢ² · cᵢ) X with Pearson residuals rᵢ = (yᵢ − μᵢ)/√V(μᵢ) and the HC correction cᵢ.

  • leverage hᵢ from the weighted hat matrix H = W¹ᐟ²X(XᵀWX)⁻¹XᵀW¹ᐟ², used by HC2/HC3/HC4.

Verified against statsmodels.GLM.fit(cov_type='HC0') (which, for GLMs, returns identical values for HC0–HC3: statsmodels does not apply leverage corrections to GLM sandwiches) and, for Gaussian/identity models, against statsmodels.OLS.fit(cov_type='HC0'...'HC3') where the corrections do apply.

Examples

>>> from aurora.models import fit_glm
>>> from aurora.inference import robust_covariance
>>> import numpy as np
>>>
>>> # Fit a model with heteroscedastic errors
>>> x = np.random.randn(100)
>>> y = 2 * x + np.random.randn(100) * (1 + 0.5 * np.abs(x))
>>> result = fit_glm(x.reshape(-1, 1), y, family='gaussian')
>>>
>>> # Get robust standard errors
>>> robust_result = robust_covariance(result, hc_type='HC3')
>>> print(robust_result.std_errors)

References

White, H. (1980). A heteroskedasticity-consistent covariance matrix estimator

and a direct test for heteroskedasticity. Econometrica, 48(4), 817-838.

MacKinnon, J.G., & White, H. (1985). Some heteroskedasticity-consistent

covariance matrix estimators with improved finite sample properties. Journal of Econometrics, 29(3), 305-325.

ANOVA

Analysis of variance (analysis of deviance) for GLM and GAM models.

This module provides ANOVA functionality for analyzing model effects and comparing nested models.

Single-model table

For a single fitted GLM, anova_glm reports partial (Type III-style) Wald tests per coefficient:

W_j = β̂_j² / Var(β̂_j)

where Var(β̂_j) is the diagonal of the model covariance matrix (already scaled by the dispersion estimate φ̂). For Gaussian models with test='F' the reference distribution is F(1, df_resid) (equivalent to the squared t-statistic, matching statsmodels.stats.anova_lm with typ=3); otherwise W_j ~ χ²₁ asymptotically.

Model comparison

For several nested models, anova_glm reports sequential deviance differences (analysis of deviance, R anova.glm convention):

  • test='Chisq' / 'LRT': ΔD ~ χ²_Δdf — valid for nested GLMs of the same family with known dispersion (Poisson, Binomial), where the deviance difference equals the likelihood-ratio statistic.

  • test='F': F = (ΔD/Δdf) / φ̂ with φ̂ taken from the largest model — the convention for Gaussian models with estimated dispersion. For non-Gaussian families test='F' falls back to the χ² test with a warning.

Examples

>>> from aurora.inference.anova import anova_glm, likelihood_ratio_test
>>>
>>> # Single-model Wald table (Type III-style)
>>> anova_glm(result, type=3)
>>>
>>> # Compare nested models
>>> anova_glm(reduced_model, full_model, test="Chisq")
>>> likelihood_ratio_test(reduced_model, full_model)

References

class aurora.inference.anova.ANOVAResult(df, ss, ms, f_statistic, p_value, source, anova_type, residual_df, residual_ss, test='F')[source]

Bases: object

Result of an ANOVA test.

Variables:

  • df (ndarray) – Degrees of freedom for each source.

  • ss (ndarray) – Sum of squares for each source. For single-model Wald tables this holds the per-coefficient Wald statistics (χ² scale); for model comparisons it holds the (non-negative) deviance differences.

  • ms (ndarray) – Mean squares (SS / df).

  • f_statistic (ndarray) – Test statistics: F values when test='F', χ² values otherwise.

  • p_value (ndarray) – P-values for each test.

  • source (list of str) – Names of sources of variation.

  • anova_type (int) – Type of ANOVA performed (1, 2, or 3).

  • residual_df (int) – Residual degrees of freedom of the (largest) model.

  • residual_ss (float) – Residual deviance (Gaussian: residual sum of squares) of the (largest) model.

  • test (str) – Reference distribution used for the p-values (‘F’ or ‘Chisq’).

Parameters:
df: ndarray
ss: ndarray
ms: ndarray
f_statistic: ndarray
p_value: ndarray
source: list[str]
anova_type: int
residual_df: int
residual_ss: float
test: str = 'F'
summary()[source]

Return formatted ANOVA table.

Return type:

str

to_dict()[source]

Convert to dictionary.

Return type:

dict[str, Any]

class aurora.inference.anova.LRTResult(statistic, df, p_value, model_names, ll_reduced, ll_full, boundary_conditions=None, boundary_correction_applied=False)[source]

Bases: object

Result of a likelihood ratio test.

Variables:

  • statistic (float) – Chi-squared test statistic (2 * (LL_full - LL_reduced)).

  • df (int) – Degrees of freedom (difference in parameters).

  • p_value (float) – P-value from chi-squared distribution.

  • model_names (tuple of str) – Names of the compared models.

  • boundary_conditions (list of str) – Variance parameters detected at or near boundary.

  • boundary_correction_applied (bool) – Whether Self & Liang (1987) mixture correction was applied.

Parameters:
statistic: float
df: int
p_value: float
model_names: tuple[str, str]
ll_reduced: float
ll_full: float
boundary_conditions: list[str] | None = None
boundary_correction_applied: bool = False
summary()[source]

Return formatted LRT summary.

Return type:

str

aurora.inference.anova.anova_glm(*models, type=3, test='F')[source]

Perform ANOVA on one or more GLM models.

Parameters:

  • *models (ModelResult) – One or more fitted model results. With one model, produces a table of partial Wald tests per coefficient. With multiple models, performs a sequential analysis of deviance (models are sorted by number of parameters).

  • type ({1, 2, 3}, default=3) – Type of sum of squares. Only type=3 (partial tests) is supported for a single model: it requires no refits and is order-independent. For multiple models the comparison is always sequential (Type I style).

  • test ({'F', 'Chisq', 'LRT'}, default='F') –

    Test statistic / reference distribution:

    • 'F': F test. For a single Gaussian model this is the squared t-statistic with F(1, df_resid) reference. For model comparison of Gaussian models, F = (ΔD/Δdf)/φ̂ with φ̂ from the largest model (R anova.glm convention). For non-Gaussian families it falls back to the χ² test with a warning.

    • 'Chisq' / 'LRT': χ² test on the Wald statistics (single model) or on the deviance differences (model comparison). The deviance difference equals the likelihood-ratio statistic for families with unit dispersion (Poisson, Binomial).

Returns:

result – ANOVA result with table and statistics.

Return type:

ANOVAResult

Examples

>>> from aurora import fit_glm, Gaussian
>>> from aurora.inference.anova import anova_glm
>>>
>>> result = fit_glm(X, y, family=Gaussian())
>>> anova_glm(result, type=3)

>>> # Compare nested models
>>> m1 = fit_glm(X[:, :2], y)
>>> m2 = fit_glm(X, y)
>>> anova_glm(m1, m2, test="Chisq")

Notes

Type III (partial) tests: each coefficient is tested adjusted for all other terms via its Wald statistic. Sequential (Type I) tests for a single model would require refitting sub-models and are not currently implemented.

aurora.inference.anova.anova(*models, type=3, test='F')

Perform ANOVA on one or more GLM models.

Parameters:

  • *models (ModelResult) – One or more fitted model results. With one model, produces a table of partial Wald tests per coefficient. With multiple models, performs a sequential analysis of deviance (models are sorted by number of parameters).

  • type ({1, 2, 3}, default=3) – Type of sum of squares. Only type=3 (partial tests) is supported for a single model: it requires no refits and is order-independent. For multiple models the comparison is always sequential (Type I style).

  • test ({'F', 'Chisq', 'LRT'}, default='F') –

    Test statistic / reference distribution:

    • 'F': F test. For a single Gaussian model this is the squared t-statistic with F(1, df_resid) reference. For model comparison of Gaussian models, F = (ΔD/Δdf)/φ̂ with φ̂ from the largest model (R anova.glm convention). For non-Gaussian families it falls back to the χ² test with a warning.

    • 'Chisq' / 'LRT': χ² test on the Wald statistics (single model) or on the deviance differences (model comparison). The deviance difference equals the likelihood-ratio statistic for families with unit dispersion (Poisson, Binomial).

Returns:

result – ANOVA result with table and statistics.

Return type:

ANOVAResult

Examples

>>> from aurora import fit_glm, Gaussian
>>> from aurora.inference.anova import anova_glm
>>>
>>> result = fit_glm(X, y, family=Gaussian())
>>> anova_glm(result, type=3)

>>> # Compare nested models
>>> m1 = fit_glm(X[:, :2], y)
>>> m2 = fit_glm(X, y)
>>> anova_glm(m1, m2, test="Chisq")

Notes

Type III (partial) tests: each coefficient is tested adjusted for all other terms via its Wald statistic. Sequential (Type I) tests for a single model would require refitting sub-models and are not currently implemented.

aurora.inference.anova.likelihood_ratio_test(model_reduced, model_full, *, names=None)[source]

Perform likelihood ratio test comparing two nested models.

Parameters:

  • model_reduced (ModelResult) – The simpler (nested) model.

  • model_full (ModelResult) – The more complex model.

  • names (tuple of str, optional) – Names for the models in output.

Returns:

result – Likelihood ratio test result.

Return type:

LRTResult

Examples

>>> from aurora.inference.anova import likelihood_ratio_test
>>>
>>> lrt = likelihood_ratio_test(reduced_model, full_model)
>>> print(lrt)

Notes

The test statistic is:

\[\chi^2 = 2 (\ell_{full} - \ell_{reduced})\]

This follows a chi-squared distribution with df equal to the difference in the number of parameters. The statistic is truncated at zero; a warning is issued when the raw difference is negative beyond numerical noise (the models are then likely not nested).

Boundary correction. When a variance component is detected at (or numerically near) the boundary of the parameter space (e.g. σ² ≈ 0), the null distribution is a 50:50 mixture of χ²_df and χ²_{df-1} (Self & Liang 1987, case 5; Stram & Lee 1994). For df = 1 the second component is χ²₀, a point mass at zero, whose survival function is 0 for any positive statistic, giving p = 0.5·P(χ²₁ > t). The correction is applied whenever any variance parameter of either model is at the boundary; this is conservative when the boundary parameter is not the one being tested.

GAMM caveats. Two warnings apply to GAMMResult inputs:

  • Non-Gaussian GAMMs are fitted by PQL, whose log-likelihood is a conditional (quasi-)likelihood given the variance components, not a marginal likelihood. Likelihood-ratio tests on it are not valid for comparing fixed effects (Pinheiro & Bates 2000, §2.4).

  • Gaussian GAMMs report the REML log-likelihood, which is not comparable between models with different fixed-effects structures.

aurora.inference.anova.lrt(model_reduced, model_full, *, names=None)

Perform likelihood ratio test comparing two nested models.

Parameters:

  • model_reduced (ModelResult) – The simpler (nested) model.

  • model_full (ModelResult) – The more complex model.

  • names (tuple of str, optional) – Names for the models in output.

Returns:

result – Likelihood ratio test result.

Return type:

LRTResult

Examples

>>> from aurora.inference.anova import likelihood_ratio_test
>>>
>>> lrt = likelihood_ratio_test(reduced_model, full_model)
>>> print(lrt)

Notes

The test statistic is:

\[\chi^2 = 2 (\ell_{full} - \ell_{reduced})\]

This follows a chi-squared distribution with df equal to the difference in the number of parameters. The statistic is truncated at zero; a warning is issued when the raw difference is negative beyond numerical noise (the models are then likely not nested).

Boundary correction. When a variance component is detected at (or numerically near) the boundary of the parameter space (e.g. σ² ≈ 0), the null distribution is a 50:50 mixture of χ²_df and χ²_{df-1} (Self & Liang 1987, case 5; Stram & Lee 1994). For df = 1 the second component is χ²₀, a point mass at zero, whose survival function is 0 for any positive statistic, giving p = 0.5·P(χ²₁ > t). The correction is applied whenever any variance parameter of either model is at the boundary; this is conservative when the boundary parameter is not the one being tested.

GAMM caveats. Two warnings apply to GAMMResult inputs:

  • Non-Gaussian GAMMs are fitted by PQL, whose log-likelihood is a conditional (quasi-)likelihood given the variance components, not a marginal likelihood. Likelihood-ratio tests on it are not valid for comparing fixed effects (Pinheiro & Bates 2000, §2.4).

  • Gaussian GAMMs report the REML log-likelihood, which is not comparable between models with different fixed-effects structures.