Generalized Linear Models (GLM)¶

This chapter covers GLM fitting in detail: distribution families, link functions, diagnostics, prediction, and inference.

Overview¶

A Generalized Linear Model relates a response variable to predictors through a link function:

\[g(\mathbb{E}[Y]) = X \beta\]

Aurora-GLM fits GLMs using Iteratively Reweighted Least Squares (IRLS). The main entry point is fit_glm():

from aurora.models.glm import fit_glm

result = fit_glm(
    X,                     # Design matrix (n, p)
    y,                     # Response vector (n,)
    family="gaussian",     # Distribution family
    link=None,             # Link function (default: family canonical)
    weights=None,          # Sample weights
    offset=None,           # Offset term
    max_iter=25,           # Max IRLS iterations
    tol=1e-8,             # Convergence tolerance
    fit_intercept=True,    # Include intercept
)

The family parameter accepts either a string or a Family instance. The link parameter works the same way.

Distribution families¶

Aurora-GLM supports 9 exponential family distributions:

Family	Use case	Variance function	Canonical link
`GaussianFamily`	Continuous outcomes	V(mu) = 1	Identity
`PoissonFamily`	Count data	V(mu) = mu	Log
`BinomialFamily`	Binary/proportion data	V(mu) = mu(1-mu)	Logit
`GammaFamily`	Positive continuous	V(mu) = mu^2	Inverse
`BetaFamily`	Proportions in (0,1)	Inverse quadratic	Logit
`InverseGaussianFamily`	Positive durations	V(mu) = mu^3	InverseSquare
`NegativeBinomialFamily`	Overdispersed counts	V(mu) = mu + mu^2/k	Log
`StudentTFamily`	Robust regression	Heavy-tailed	Identity
`TweedieFamily`	Zero-inflated continuous	Variance power law	Log

Families can be specified by string or instance:

from aurora.distributions.families import GaussianFamily, PoissonFamily

# String form (uses default link)
result = fit_glm(X, y, family="poisson")

# Instance form (uses default link)
result = fit_glm(X, y, family=PoissonFamily())

Negative Binomial with estimated dispersion. For overdispersed count data where the dispersion theta is unknown, pass NegativeBinomialFamily(theta="estimate"). fit_glm then runs the two-step MASS::glm.nb flow: a provisional Poisson-style fit with theta = 1, maximum-likelihood estimation of theta from the fitted means (Lawless 1987), and a final refit with the estimated theta. The estimate is exposed as result.theta_:

from aurora.distributions.families import NegativeBinomialFamily
from aurora.models.glm import fit_glm

# theta estimated from the data (glm.nb flow)
result = fit_glm(X, y, family=NegativeBinomialFamily(theta="estimate"))
print(result.theta_)  # estimated dispersion

# Fixed theta behaves as a standard GLM family (result.theta_ is None)
result = fit_glm(X, y, family=NegativeBinomialFamily(theta=2.0))

Link functions¶

Nine link functions connect the mean to the linear predictor:

Link	Formula	Typical use	Inverse
`IdentityLink`	g(mu) = mu	Gaussian (OLS)	g^-1(eta) = eta
`LogLink`	g(mu) = log(mu)	Poisson, Gamma	g^-1(eta) = exp(eta)
`LogitLink`	g(mu) = log(mu/(1-mu))	Binomial	g^-1(eta) = 1/(1+exp(-eta))
`ProbitLink`	g(mu) = Phi^-1(mu)	Binomial (probit)	g^-1(eta) = Phi(eta)
`CLogLogLink`	g(mu) = log(-log(1-mu))	Binomial (extreme)	g^-1(eta) = 1-exp(-exp(eta))
`InverseLink`	g(mu) = 1/mu	Gamma (canonical)	g^-1(eta) = 1/eta
`InverseSquareLink`	g(mu) = 1/mu^2	Inverse Gaussian	g^-1(eta) = 1/sqrt(eta)
`SqrtLink`	g(mu) = sqrt(mu)	Count data alt	g^-1(eta) = eta^2
`PowerLink`	g(mu) = mu^p	General purpose	g^-1(eta) = eta^(1/p)

Override the default link:

from aurora.distributions.links import LogLink, LogitLink

# Poisson with identity link (instead of canonical log)
result = fit_glm(X, y, family="poisson", link=LogLink())

# Binomial with probit link (instead of canonical logit)
result = fit_glm(X, y, family="binomial", link="probit")

Inference¶

Standard errors, p-values, and confidence intervals are computed lazily on first access:

result = fit_glm(X, y)

# Standard errors (computed from Fisher information)
print(result.std_errors_)

# Wald p-values (two-sided)
print(result.p_values_)

# Coefficient covariance matrix
print(result.coef_cov_)

# Confidence intervals
mean, lower, upper = result.predict(X_new, interval="confidence", level=0.95)

Diagnostics¶

The glm_diagnostics() function computes comprehensive residual diagnostics:

from aurora.inference.diagnostics import glm_diagnostics

diag = glm_diagnostics(result)
print(diag.summary())

Available diagnostics include deviance residuals, Pearson residuals, leverage values, Cook’s distance, and DFBETAS.

Prediction¶

The .predict() method supports both the response and link scale:

# Predict on response scale (default)
y_pred = result.predict(X_new)

# Predict on link scale
eta_pred = result.predict(X_new, type="link")

# With 95% Wald confidence intervals
mean, lower, upper = result.predict(X_new, interval="confidence", level=0.95)

Comparison with statsmodels¶

Aurora-GLM produces results validated against both statsmodels and R:

import numpy as np
from aurora.models.glm import fit_glm

# Aurora-GLM
X = np.random.randn(200, 3)
y = 1.5 + X @ [0.8, -0.5, 0.3] + np.random.randn(200) * 0.5
result = fit_glm(X, y)

# Results match statsmodels within 1e-6 tolerance
# Coefficients, deviance, AIC all agree with R's glm()