Generalized Linear Models

Generate $n = 80$ predictor values $x_i \sim \mathcal{N}(0, 1)$ and binary responses $y_i \sim \text{Bernoulli}(\sigma(0.5 + 1.8 x_i))$ where $\sigma$ is the logistic sigmoid. Fit OLS: $\hat y_i = \hat\beta_0 + \hat\beta_1 x_i$. The fitted line is approximately $\hat y \approx 0.5 + 0.30 x$, which crosses zero at $x \approx -1.7$ and exceeds one at $x \approx 1.7$. Roughly 12 of the 80 fitted values are negative; roughly 14 exceed unity. The model is internally inconsistent: it predicts probabilities outside $[0, 1]$, and any decision rule built on $\hat y$ either has to clip those values or admit the model is wrong about its own outputs.

The same data fit by logistic regression with the logit link $g(p) = \log(p/(1-p))$ produces $\hat\beta_0 \approx 0.51, \hat\beta_1 \approx 1.94$ on the latent scale, with $\hat p_i = \sigma(\hat\beta_0 + \hat\beta_1 x_i) \in (0, 1)$ for every $i$. The probabilities are bounded by construction; coefficient interpretation shifts from “additive on $y$” to “additive on $\log\text{odds}(y=1)$.”

The three failures above are not isolated — they are the three structural assumptions of the Normal linear model failing one at a time:

1. Bounded support. Binary, multinomial, and ordinal outcomes have support that does not match $\mathbb{R}$. The link function maps the bounded mean (e.g. probability) back to $\mathbb{R}$.
2. Variance–mean dependence. Counts have $\text{Var}(Y) = E[Y]$ (Poisson); proportions have $\text{Var}(Y) = E[Y](1 - E[Y])/n$ (binomial); positive-skewed amounts have $\text{Var}(Y) \propto E[Y]^2$ (Gamma). Homoscedasticity is the exception, not the rule. The variance function $V(\mu)$ encodes this dependence.
3. Non-Normal tails. Skewed positives, heavy-tailed counts, zero-inflated mixtures — none have Normal sampling distributions even at moderate $n$. The exponential-family response replaces “Normal errors” with “the right family for this support.”

GLMs address all three at once by composing a linear predictor with a link and pairing it with an exponential-family response. The linear-predictor structure is what survives; everything else generalizes.

Before 1972, statisticians had logistic regression (Bliss 1935, Berkson 1944), probit regression (Bliss 1934), Poisson regression (Cochran 1940, Birch 1963), log-linear models for contingency tables (Goodman 1968), and Gamma regression for survival times (Cox 1972) — all developed independently with bespoke fitting algorithms. Nelder & Wedderburn (1972) observed that all of these are the same model: an exponential-family response with mean $\mu$, a link function $g$ mapping $\mu$ to a linear predictor $\eta = \mathbf{x}^\top\boldsymbol\beta$, and a single fitting algorithm — iteratively reweighted least squares (IRLS) — that handles every family. The unification was both conceptual (one framework, not five) and practical (one algorithm to implement, not five). McCullagh & Nelder (1989) is the canonical book-length development; Agresti (2015) is the current graduate standard.

Topic 22 reproduces the unification mathematically. §22.2 lays out the framework; §22.3 proves the canonical-link score identity that drives the unified algorithm; §22.4 / §22.5 / §22.6 fit logistic / Poisson / Gamma regression as worked examples; §22.7 introduces the deviance as the universal goodness-of-fit statistic; §22.8 / §22.9 handle inference and misspecification.

Given a design matrix $\mathbf{X} \in \mathbb{R}^{n \times (p+1)}$ (with a leading column of ones for the intercept) and a coefficient vector $\boldsymbol\beta \in \mathbb{R}^{p+1}$, the linear predictor for observation $i$ is

In stacked form, $\boldsymbol\eta = \mathbf{X}\boldsymbol\beta \in \mathbb{R}^n$. The linear predictor lives on the real line — no constraints — even when the response $Y_i$ has bounded or discrete support.

A link function $g: \mathcal{M} \to \mathbb{R}$ is a smooth, strictly monotone map from the response’s mean space $\mathcal{M}$ (e.g. $(0, 1)$ for Bernoulli, $(0, \infty)$ for Poisson and Gamma, $\mathbb{R}$ for Normal) to the linear-predictor scale $\mathbb{R}$. It satisfies

Strict monotonicity guarantees the inverse $g^{-1}$ is well-defined, so the model uniquely specifies $\mu_i$ given $\boldsymbol\beta$ and $\mathbf{x}_i$. Smoothness is needed for the IRLS derivation in §22.3.

The response $Y_i$ follows an exponential family with density (or mass function)

where $\theta_i$ is the natural parameter, $A(\theta_i)$ is the cumulant function, $\phi$ is the dispersion, and $h$ is a base-measure factor that does not depend on $\theta$. By Topic 7 Thm 3 and Cor 2, $E[Y_i] = A'(\theta_i) = \mu_i$ and $\text{Var}(Y_i) = \phi A''(\theta_i) = \phi V(\mu_i)$, where $V(\mu) = A''((A')^{-1}(\mu))$ is the variance function.

A link function is canonical when $g = (A')^{-1}$, i.e. when the natural parameter coincides with the linear predictor: $\theta_i = \eta_i = \mathbf{x}_i^\top\boldsymbol\beta$. Under the canonical link the score equation collapses to the design-times-residual identity that makes §22.3 Thm 1 the featured theorem of the topic.

Take $Y_i \sim \mathcal{N}(\mu_i, \sigma^2)$ with the identity link $g(\mu) = \mu$. The exponential-family form has $\theta = \mu$ (so the link is canonical), $A(\theta) = \theta^2/2$, $h(y, \phi) = (2\pi\phi)^{-1/2}\exp(-y^2/(2\phi))$, and dispersion $\phi = \sigma^2$. The cumulant identities give $A'(\theta) = \theta = \mu$ (mean) and $A''(\theta) = 1$, so $V(\mu) = 1$ — constant variance, recovering homoscedasticity.

The log-likelihood is

so maximizing $\ell$ over $\boldsymbol\beta$ is equivalent to minimizing $\sum_i (y_i - \mathbf{x}_i^\top\boldsymbol\beta)^2$ — the OLS objective from §21.2. The score equation $\nabla_{\boldsymbol\beta}\ell = (1/\sigma^2)\mathbf{X}^\top(\mathbf{y} - \mathbf{X}\boldsymbol\beta) = \mathbf{0}$ is $\mathbf{X}^\top\mathbf{X}\hat{\boldsymbol\beta} = \mathbf{X}^\top\mathbf{y}$, which is §21.4 Thm 3 (normal equations) verbatim. Topic 21’s entire estimation theory is the Normal+identity-link special case of the GLM framework. We will see this reduction again, more abstractly, in §22.3 Proof 1 Step 4.

For the four GLM families covered in this topic, the canonical link is determined by inverting $A'$:

The Bernoulli logit and Poisson log are the two canonical links you will see most often in practice; the Gamma inverse link is unusual and §22.6 Rem 14 discusses why the log link is often preferred there despite being non-canonical.

The canonical link is mathematically privileged — it is the unique link under which the score equation is exactly $\mathbf{X}^\top(\mathbf{y} - \boldsymbol\mu) = \mathbf{0}$ (§22.3 Thm 1) and observed Fisher information equals expected Fisher information (§22.3 Thm 2 Step 2). Non-canonical links — probit and cloglog for binary data, identity link for Poisson, log link for Gamma — break both properties: the score acquires an extra Jacobian factor and the Hessian’s expected and observed forms diverge.

So why use them? Three reasons. Latent-variable interpretation: the probit link $g(p) = \Phi^{-1}(p)$ corresponds to a latent Normal variable whose threshold determines the binary outcome — natural for psychological-test response models. Hazard interpretation: the cloglog link $g(p) = \log(-\log(1-p))$ corresponds to proportional hazards in discrete time — natural for survival analysis. Coefficient interpretability: the log link on Gamma gives $\hat\beta_j$ as a log-multiplicative effect on the mean (a one-unit change in $x_j$ multiplies $\mu$ by $e^{\hat\beta_j}$), which is more useful in practice than the inverse-link’s “additive on $1/\mu$” interpretation.

In the GLM framework, non-canonical links work — IRLS still converges, asymptotic normality still holds — they just do not enjoy the elegant identities that make the canonical case the featured theorem of §22.3.

The dispersion $\phi$ scales the variance: $\text{Var}(Y_i) = \phi V(\mu_i)$. For Bernoulli and Poisson, $\phi$ is fixed at 1 by the structure of the family — there is no separate variance parameter to estimate, because mean determines variance. For Gamma and Normal, $\phi$ is a free parameter that must be estimated alongside $\boldsymbol\beta$. The estimator is the Pearson statistic:

When $\phi$ is estimated rather than known, the deviance test (§22.7) uses an $F$-reference instead of $\chi^2$, exactly as in Topic 21 §21.8 where $\hat\sigma^2$ replaced $\sigma^2$ in the $F$-statistic. The §22.9 sandwich estimator generalizes this further: when the variance specification $\phi V(\mu_i)$ is wrong (overdispersed Poisson, misspecified Gamma scale, clustered Bernoulli), the QMLE is still consistent for $\boldsymbol\beta$ as long as the mean specification $\mu_i = g^{-1}(\mathbf{x}_i^\top\boldsymbol\beta)$ is correct; only the variance estimate needs adjustment.

Let $Y_i$ follow an exponential-family density $f(y_i \mid \theta_i, \phi) = h(y_i, \phi)\exp\!\left((\theta_i y_i - A(\theta_i))/\phi\right)$ with canonical link $\theta_i = \eta_i = \mathbf{x}_i^\top\boldsymbol\beta$. Then:

1. Score identity. $\displaystyle \nabla_{\boldsymbol\beta}\ell(\boldsymbol\beta) = \frac{1}{\phi}\mathbf{X}^\top(\mathbf{y} - \boldsymbol\mu(\boldsymbol\beta))$, where $\mu_i = A'(\theta_i) = A'(\mathbf{x}_i^\top\boldsymbol\beta)$.

2. Hessian. $\displaystyle \mathcal{H}(\boldsymbol\beta) = -\frac{1}{\phi}\mathbf{X}^\top\mathbf{W}(\boldsymbol\beta)\mathbf{X}$, where $\mathbf{W}(\boldsymbol\beta)$ is diagonal with $W_{ii} = V(\mu_i) > 0$.

3. Concavity. When $\mathbf{X}$ has full column rank $p+1$, the log-likelihood $\ell$ is strictly concave in $\boldsymbol\beta$. Consequently, any interior critical point of $\ell$ is the unique global maximum.

In particular, the MLE $\hat{\boldsymbol\beta}$, when it exists in the interior of the parameter space, is the unique solution to $\mathbf{X}^\top(\mathbf{y} - \boldsymbol\mu(\hat{\boldsymbol\beta})) = \mathbf{0}$.

For the logistic regression model $Y_i \sim \text{Bernoulli}(p_i)$ with logit link $\eta_i = \log[p_i/(1-p_i)] = \mathbf{x}_i^\top\boldsymbol\beta$, the canonical-link score identity gives

where $\mathbf{p} = (p_1, \ldots, p_n)^\top$ with $p_i = \sigma(\mathbf{x}_i^\top\boldsymbol\beta)$ and $\sigma(z) = 1/(1 + e^{-z})$ is the sigmoid. Setting the score to zero gives the system $\sum_i \mathbf{x}_i (y_i - p_i) = \mathbf{0}$ — the binary-classification analog of “residual orthogonal to the design columns.” Unlike OLS, the system is non-linear in $\boldsymbol\beta$ (because $p_i$ depends on $\boldsymbol\beta$ through the sigmoid), so there is no closed-form solution; we need the IRLS algorithm of Thm 2.

The Hessian is $\mathcal{H} = -\mathbf{X}^\top\mathbf{W}\mathbf{X}$ with $W_{ii} = p_i(1 - p_i)$ — the variance of a Bernoulli with parameter $p_i$. The weight is largest at $p_i = 0.5$ (variance $0.25$) and shrinks to zero at the boundary $p_i \in \{0, 1\}$, which is the source of the “near-separation” pathology Rem 8 will treat in §22.4.

Under the canonical link, Newton’s method on the GLM log-likelihood produces iterates

where $\mathbf{W}^{(t)} = \text{diag}(V(\mu_i^{(t)}))$ and the adjusted response is

The right-hand side is the closed-form weighted-least-squares estimator for the linear model $\mathbf{z}^{(t)} \approx \mathbf{X}\boldsymbol\beta$ with weights $\mathbf{W}^{(t)}$. Therefore Newton’s method on the canonical-link GLM log-likelihood is iteratively reweighted least squares. Furthermore, because $\mathbf{W}^{(t)}$ depends on $\boldsymbol\beta$ only through the deterministic mean $\boldsymbol\mu^{(t)}$ (no randomness in $\mathbf{W}$), the observed information equals the expected information, so Newton’s method is Fisher scoring.

The IRLS update has a deeper interpretation in information geometry: the exponential family is a Riemannian manifold whose metric is the Fisher information matrix $\mathcal{I}(\boldsymbol\beta) = \mathbf{X}^\top\mathbf{W}\mathbf{X}$, and the IRLS step is the natural-gradient step under that metric — gradient ascent in the geometry the family naturally lives in, rather than the flat Euclidean geometry of $\boldsymbol\beta \in \mathbb{R}^{p+1}$. The natural gradient is invariant under reparameterization (a property the ordinary gradient lacks), which is exactly the kind of invariance that makes IRLS work uniformly across families and links.

The deep-learning community rediscovered natural-gradient descent in Amari (1998) for neural networks, where it is approximated by methods like K-FAC and Shampoo. The full information-geometry framing, including connections to the dual coordinate systems of exponential families, is at formalML: Information geometry .

We can now make the “OLS is GLM” claim of Step 4 of Proof 1 more emphatic. For Normal response with identity link and $\phi = \sigma^2$: $A(\theta) = \theta^2/2$ gives $V(\mu) = A''(\theta) = 1$, so $\mathbf{W} = \mathbf{I}$. The IRLS update collapses to

The adjusted response collapses to $\mathbf{y}$, so the WLS solve becomes the OLS solve, and one iteration converges. Topic 21’s normal equations are not just contained in the canonical-link score identity — Topic 21’s fitting algorithm is contained in IRLS as the special case where $\mathbf{W} = \mathbf{I}$ and the iteration converges in one step. Gauss–Markov BLUE optimality (§21.6) is the homoscedastic-case statement; under heteroscedastic errors the right BLUE is WLS with the true weights, which is what IRLS with the true variance function delivers (Topic 21 §21.9 Rem 21 forward-promised this; we have now discharged it).

Under GLM regularity (canonical link, $\mathbf{X}$ full column rank, MLE $\hat{\boldsymbol\beta}$ in the interior of the parameter space, finite Fisher information), as $n \to \infty$,

where $\mathcal{I}(\boldsymbol\beta_0) = (1/n)\mathbf{X}^\top\mathbf{W}(\boldsymbol\beta_0)\mathbf{X}$ is the per-observation Fisher information matrix. Equivalently, in finite-sample-approximation form,

The estimated covariance $\widehat{\text{Cov}}(\hat{\boldsymbol\beta}) = \hat\phi (\mathbf{X}^\top\hat{\mathbf{W}}\mathbf{X})^{-1}$ replaces the unknown $\boldsymbol\beta_0$ by $\hat{\boldsymbol\beta}$ and the unknown $\phi$ by the Pearson estimator (§22.2 Rem 4), and is what every coefficient SE in §22.8 is built on.

On a synthetic logistic dataset with $n = 60$, $p = 1$ predictor, and true $\boldsymbol\beta = (-0.5, 1.8)^\top$, IRLS started from $\boldsymbol\beta^{(0)} = (0, 0)^\top$ produces:

Five iterations to convergence at $\|\Delta\boldsymbol\beta\| < 10^{-8}$. The error $\|\boldsymbol\beta^{(t)} - \hat{\boldsymbol\beta}\|$ shrinks roughly quadratically per step ($1.51 \to 0.39 \to 0.026 \to \cdots$, each step squaring the previous gap up to a constant) — the signature of Newton-method convergence in a neighborhood of the optimum (Topic 14 Thm 6). The right-panel readout of IRLSVisualizer reproduces this trajectory live.

Given binary responses $Y_i \in \{0, 1\}$ and predictors $\mathbf{x}_i \in \mathbb{R}^{p+1}$ (with leading intercept column), the logistic regression model is

Equivalently, $p_i = \sigma(\mathbf{x}_i^\top\boldsymbol\beta)$ where $\sigma(z) = 1/(1 + e^{-z})$ is the logistic sigmoid. The coefficient $\beta_j$ is the log-odds-ratio associated with a one-unit increase in $x_j$ (holding others fixed): $e^{\beta_j}$ is the multiplicative change in the odds $p/(1-p)$.

Under GLM regularity, the asymptotic distribution $\hat{\boldsymbol\beta} \overset{\text{approx}}{\sim} \mathcal{N}_{p+1}(\boldsymbol\beta_0, (\mathbf{X}^\top\hat{\mathbf{W}}\mathbf{X})^{-1})$ from Thm 3 gives the per-coefficient Wald-z confidence interval

with asymptotic coverage $1 - \alpha$. The Wald test $z_j = \hat\beta_j/\widehat{\text{SE}}(\hat\beta_j)$ for $H_0: \beta_j = 0$ is the dual test (Topic 19 Thm 1). For $\phi = 1$ Bernoulli, no F-correction is needed; the $z$ asymptotic applies directly.

A consumer-finance synthetic example: $n = 200$ loan applicants with two predictors: standardized age $(z_{age})$ and standardized log-income $(z_{inc})$. The default response is $Y_i = 1$ for default, $Y_i = 0$ for repayment, with true coefficients $\boldsymbol\beta_0 = (-0.5, 1.0, 2.0)$ on the logit scale (intercept, age, income). The notebook fits the model with sm.GLM(y, X, family=Binomial()):

Both predictors are highly significant. Coefficient interpretation: a one-standard-deviation increase in age raises log-odds-of-default by $1.32$ (multiplicative odds change $e^{1.32} \approx 3.74$); the analogous income coefficient is $e^{1.87} \approx 6.48$. The full model deviance is $D_{\text{full}} = 158.96$ on $n - p - 1 = 197$ residual df, vs the intercept-only deviance $D_{\text{null}} = 264.63$ on $199$ df. The deviance reduction $D_{\text{null}} - D_{\text{full}} = 105.67$ on 2 df has $p \ll 0.001$ under $\chi^2_2$ — overwhelming evidence the predictors carry information about default risk.

The logit is canonical, but two other links see practical use. Probit: $g(p) = \Phi^{-1}(p)$ where $\Phi$ is the standard Normal CDF; the model corresponds to a latent Normal variable $Y_i^* = \mathbf{x}_i^\top\boldsymbol\beta + \varepsilon_i$ with $\varepsilon_i \sim \mathcal{N}(0, 1)$ and $Y_i = \mathbb{1}[Y_i^* > 0]$. Common in psychometric / response-threshold modeling. Complementary log-log (cloglog): $g(p) = \log(-\log(1 - p))$ corresponds to a discrete-time hazard interpretation — natural for survival-analysis-style binary outcomes. Both are non-canonical: IRLS still fits them (with weights $W_{ii} = (g'(\mu_i))^{-2}/V(\mu_i)$ rather than the simpler $V(\mu_i)$ form), and asymptotic normality still applies, but the score identity is no longer $\mathbf{X}^\top(\mathbf{y} - \boldsymbol\mu)$ exactly. In practice, logit and probit give nearly identical fitted probabilities and indistinguishable LRT $p$-values; the choice is interpretive (log-odds vs latent-Normal) rather than statistical.