Generalized Linear Models

Intuition Beginner

Ordinary straight-line fitting answers one question: how does the average of some measured number move as you turn a few dials? You pick weights on the dials, add them up, and that sum is your prediction of the average. It works beautifully when the thing you measure is a quantity that can wander up or down freely — a height, a price, a temperature — and when the scatter around the average is the same size everywhere.

But much of what you want to predict is not like that. A count of accidents at an intersection cannot go negative, and a busy intersection scatters more than a quiet one. A yes-or-no outcome, like whether a loan defaults, is a probability trapped between $0$ and $1$. A waiting time is positive and skewed. For these, the plain straight-line recipe predicts impossible values and misjudges the scatter.

A generalized linear model keeps the part that works — a weighted sum of the dials — and fixes the two parts that break. First, it runs the weighted sum through a link, a simple transform that bends the unbounded sum onto the right range: through a logarithm so a count stays positive, through an S-curve so a probability stays inside $0$ to $1$. Second, it lets the scatter grow or shrink with the average in the way each kind of data actually behaves. One framework, many response types, all fitted by the same procedure.

Visual Beginner

The same skeleton serves every case: build a weighted sum, bend it onto the right range with a link, and read off the average. Only the link and the scatter-pattern change from row to row.

Response type	Range of the average	Link (the bend)	How scatter grows
Free quantity (height, price)	any number	none — straight line	same everywhere
Yes/no outcome (default)	between $0$ and $1$	S-curve (logit)	largest near one-half
Count (accidents)	$0$ and up	logarithm	grows with the average
Waiting time (positive, skewed)	above $0$	logarithm or reciprocal	grows with the square of the average

The picture to keep: one engine, the weighted sum, feeding a chosen bend, feeding a needle whose wobble matches the kind of measurement. Plain regression is the top row with the bend removed; everything else is the same engine wearing a different link and a different wobble.

Worked example Beginner

A shop counts customers per hour and wants to predict the count from a single dial, the temperature outside. Counts cannot be negative, so we bend the weighted sum through a logarithm: the log of the average count is the straight-line part. Suppose fitting gives the rule

$$\log (\text{average count})=1.0+0.05\times ({\text{temperature in }}^{\circ }\text{C}).$$

Step 1. Predict at $20^{\circ }$C. The straight-line part is $1.0+0.05\times 20=1.0+1.0=2.0$.

Step 2. Undo the logarithm to get the average count. The opposite of a logarithm is the exponential, so the average count is $e^{2.0}\approx 7.39$. About $7$ customers an hour.

Step 3. Predict at $30^{\circ }$C. The straight-line part is $1.0+0.05\times 30=2.5$, and the average count is $e^{2.5}\approx 12.18$. About $12$ customers an hour.

Step 4. Read the effect of the dial. A $10^{\circ }$ rise added $0.5$ to the log, and adding $0.5$ to a logarithm multiplies the count by $e^{0.5}\approx 1.65$. Check: $7.39\times 1.65\approx 12.18$, matching Step 3.

What this tells us: because the link is a logarithm, a fixed change in the dial does not add a fixed number of customers — it multiplies the count by a fixed factor. The bend turns the additive straight line into a multiplicative effect on the count, which is exactly what keeps the prediction positive.

Check your understanding Beginner

Formal definition Intermediate+

A generalized linear model (GLM) for a response $Y$ given covariates $x\in {\mathbb{R}}^p$ is specified by three components, following McCullagh & Nelder ^{[McCullagh 1989]} Ch. 2.

The random component (the exponential-dispersion family). The response has a density or mass function of exponential-dispersion form

$$f(y;\theta ,\varphi )=\exp \{\frac{y\theta -b(\theta )}{a(\varphi )}+c(y,\varphi )\},$$

with $\theta$ the canonical parameter, $\varphi >0$ the dispersion parameter, $b$ the (smooth, strictly convex) cumulant function, and $a,c$ known functions; commonly $a(\varphi )=\varphi /w$ for a known prior weight $w$. This is the one-parameter exponential family of 45.01.02 indexed in its mean, carrying the extra scale $\varphi$. Differentiating the identity $\int fdy=1$ under the integral sign gives the mean-variance relationship

$$\mathbb{E}[Y]=\mu =b^{\prime }(\theta ),\mathrm{Var}(Y)=a(\varphi )b^{\prime \prime }(\theta )=a(\varphi )V(\mu ),$$

where the variance function $V(\mu ):=b^{\prime \prime }(\theta (\mu ))$ records how the variance depends on the mean through the inverse of the strictly increasing map $\theta \mapsto \mu =b^{\prime }(\theta )$. The variance function, together with $a(\varphi )$, determines the second-moment structure entirely.

The systematic component. A linear predictor $\eta =x^{\top }\beta$ assembles the covariates through unknown coefficients $\beta \in {\mathbb{R}}^p$. Stacked over $n$ observations this is $\eta =X\beta$ with design matrix $X\in {\mathbb{R}}^{n\times p}$.

The link. A monotone, differentiable link function $g$ connects the mean to the linear predictor,

$$g(\mu )=\eta =x^{\top }\beta ,\mu =g^{-1}(\eta ).$$

The link is canonical when $g=(b^{\prime }{)}^{-1}$, equivalently when $\theta =\eta$: the canonical parameter equals the linear predictor. For the canonical link the natural statistic $X^{\top }y$ is sufficient for $\beta$, which is the structural reason the canonical-link fit is the cleanest.

Standard members. Each pair (family, canonical link) is a line of the framework:

• Gaussian, identity link. $b(\theta )={\theta }^2/2$, $V(\mu )=1$, $g(\mu )=\mu$. This is ordinary linear regression 45.06.01.
• Bernoulli/binomial, logit link. $b(\theta )=\log (1+e^{\theta })$, $V(\mu )=\mu (1-\mu )$, $g(\mu )=\log \{\mu /(1-\mu )\}$. This is logistic regression 45.08.03.
• Poisson, log link. $b(\theta )=e^{\theta }$, $V(\mu )=\mu$, $g(\mu )=\log \mu$. This is log-linear regression for counts.
• Gamma, reciprocal (canonical) or log link. $b(\theta )=-\log (-\theta )$, $V(\mu )={\mu }^2$, $g(\mu )=1/\mu$ (or $\log \mu$ in practice).

Counterexamples to common slips

• The link transforms the mean, not the data. Modelling $g(Y)$ by least squares is not a GLM: the GLM applies $g$ to $\mathbb{E}[Y]$, leaving the response on its own scale and its own variance function. Log-transforming $Y$ and the log-link Poisson GLM give different estimators with different mean targets ($\mathbb{E}[\log Y]$ versus $\log \mathbb{E}[Y]$).
• Canonical is a choice, not a necessity. The probit and complementary-log-log links are valid for binary data though neither is canonical; only the logit is. Using a non-canonical link costs the sufficiency of $X^{\top }y$ and makes observed and expected information differ, but the model is well-posed.
• The variance function is not a free parameter. Once the family is fixed, $V(\mu )$ is determined by $b$; one cannot pair Poisson's $V(\mu )=\mu$ with a Gaussian likelihood. Postulating $V$ without a full family is the separate quasi-likelihood device.
• Constant dispersion is an assumption. The factor $a(\varphi )$ is taken known and common ($\varphi =1$ for Poisson and binomial). When the data scatter more than $\varphi =1$ allows — overdispersion — the nominal standard errors are too small and $\varphi$ must be estimated.

Key theorem with proof Intermediate+

The unifying computational fact is that the maximum-likelihood fit of any GLM is a sequence of weighted least-squares fits: Fisher scoring on the GLM log-likelihood is iteratively reweighted least squares (IRLS), with weights and a working response that the current fitted means supply. This generalises the logistic case of 45.08.03 to every exponential-dispersion family and any monotone link. The result follows McCullagh & Nelder ^{[McCullagh 1989]} §2.5 and Nelder & Wedderburn ^{[Nelder 1972]}.

Theorem (GLM fitting by IRLS). Let $y_1,\dots ,y_n$ be independent with $y_i\sim$ exponential-dispersion density of mean ${\mu }_i$, $a_i(\varphi )=\varphi /w_i$, linear predictor ${\eta }_i=x_i^{\top }\beta$, and link $g({\mu }_i)={\eta }_i$. Write $\mu =({\mu }_i)$, $\eta =({\eta }_i)$, and define the diagonal Fisher weights and working response

$$W=\mathrm{diag}\left(\frac{w_i}{V({\mu }_i)g^{\prime }({\mu }_i{)}^2}\right),z=\eta +(y-\mu )\odot g^{\prime }(\mu ),$$

where $g^{\prime }(\mu )=(g^{\prime }({\mu }_i))$ acts entrywise. Then the score has the form $\nabla \ell (\beta )=\frac{1}{\varphi }{X}^{\top }W(z-\eta )$, the expected (Fisher) information is $\mathcal{I}(\beta )=\frac{1}{\varphi }{X}^{\top }WX$, and the Fisher-scoring update equals the weighted least-squares solution

$${\beta }^{\mathrm{new}}=(X^{\top }WX{)}^{-1}{X}^{\top }Wz.$$

Each step regresses the working response $z$ on $X$ with weights $W$; for the canonical link the observed and expected information coincide, so Fisher scoring equals Newton-Raphson and the log-likelihood is concave.

Proof. The log-likelihood contribution of observation $i$ is ${\ell }_i=\{y_i{\theta }_i-b({\theta }_i)\}/a_i(\varphi )+c(y_i,\varphi )$. By the chain rule through ${\theta }_i\to {\mu }_i\to {\eta }_i\to \beta$,

$$\frac{\partial {\ell }_i}{\partial \beta }=\frac{\partial {\ell }_i}{\partial {\theta }_i}\frac{d{\theta }_i}{d{\mu }_i}\frac{d{\mu }_i}{d{\eta }_i}\frac{\partial {\eta }_i}{\partial \beta }.$$

Each factor is explicit: $\partial {\ell }_i/\partial {\theta }_i=(y_i-b^{\prime }({\theta }_i))/a_i(\varphi )=(y_i-{\mu }_i)/a_i(\varphi )$ since ${\mu }_i=b^{\prime }({\theta }_i)$; $d{\mu }_i/d{\theta }_i=b^{\prime \prime }({\theta }_i)=V({\mu }_i)$ so $d{\theta }_i/d{\mu }_i=1/V({\mu }_i)$; $d{\mu }_i/d{\eta }_i=1/g^{\prime }({\mu }_i)$ because $g({\mu }_i)={\eta }_i$; and $\partial {\eta }_i/\partial \beta =x_i$. Multiplying and substituting $a_i(\varphi )=\varphi /w_i$,

$$\frac{\partial {\ell }_i}{\partial \beta }=\frac{(y_i-{\mu }_i)}{\varphi /w_i}\cdot \frac{1}{V({\mu }_i)}\cdot \frac{1}{g^{\prime }({\mu }_i)}{x}_i=\frac{1}{\varphi }\frac{w_i(y_i-{\mu }_i)}{V({\mu }_i)g^{\prime }({\mu }_i)}{x}_i.$$

Define $W_i=w_i/\{V({\mu }_i)g^{\prime }({\mu }_i{)}^2\}$, the $i$-th diagonal of $W$. Then $\partial {\ell }_i/\partial \beta =\frac{1}{\varphi }{W}_i{g}^{\prime }({\mu }_i)(y_i-{\mu }_i)x_i$, and summing gives the score

$$\nabla \ell (\beta )=\frac{1}{\varphi }\sum _i{W}_i{g}^{\prime }({\mu }_i)(y_i-{\mu }_i)x_i=\frac{1}{\varphi }{X}^{\top }W[g^{\prime }(\mu )\odot (y-\mu )]=\frac{1}{\varphi }{X}^{\top }W(z-\eta ),$$

where $z-\eta =g^{\prime }(\mu )\odot (y-\mu )$ is the definition of the working response. Setting the score to zero gives the GLM score equation $X^{\top }W(z-\eta )=0$.

For the information, compute the negative expected second derivative. Differentiating the per-observation score $s_i(\beta )=\frac{1}{\varphi }{W}_i{g}^{\prime }({\mu }_i)(y_i-{\mu }_i)x_i$ in $\beta$ produces two groups of terms: one proportional to $(y_i-{\mu }_i)$ and one to $\partial {\mu }_i/\partial \beta$. The first group has mean zero because $\mathbb{E}[y_i-{\mu }_i]=0$, so taking expectations annihilates it. The surviving term is, using $\partial {\mu }_i/\partial \beta =(1/g^{\prime }({\mu }_i))x_i$,

$$\mathbb{E}\left[-\frac{\partial s_i}{\partial {\beta }^{\top }}\right]=\frac{1}{\varphi }{W}_i{g}^{\prime }({\mu }_i)\frac{\partial {\mu }_i}{\partial {\beta }^{\top }}{x}_i=\frac{1}{\varphi }{W}_i{x}_i{x}_i^{\top }.$$

Summing, $\mathcal{I}(\beta )=\frac{1}{\varphi }{X}^{\top }WX$. The Fisher-scoring step replaces the Hessian by $-\mathcal{I}$:

$${\beta }^{\mathrm{new}}=\beta +\mathcal{I}(\beta {)}^{-1}\nabla \ell (\beta )=\beta +(X^{\top }WX{)}^{-1}{X}^{\top }W(z-\eta ).$$

Because $\eta =X\beta$, the identity $\beta =(X^{\top }WX{)}^{-1}{X}^{\top }WX\beta =(X^{\top }WX{)}^{-1}{X}^{\top }W\eta$ holds, so

$${\beta }^{\mathrm{new}}=(X^{\top }WX{)}^{-1}[X^{\top }W\eta +X^{\top }W(z-\eta )]=(X^{\top }WX{)}^{-1}{X}^{\top }Wz,$$

the weighted normal equation for regressing $z$ on $X$ with weights $W$. The dispersion $\varphi$ cancels from the update, entering only the scale of the information. For the canonical link, $g=(b^{\prime }{)}^{-1}$ gives $g^{\prime }(\mu )=1/b^{\prime \prime }(\theta )=1/V(\mu )$, whence $\theta =\eta$ and the per-observation Hessian ${\partial }^2{\ell }_i/\partial \beta \partial {\beta }^{\top }=-\frac{w_i}{\varphi }{b}^{\prime \prime }({\theta }_i)x_i{x}_i^{\top }$ is already non-random — the $(y_i-{\mu }_i)$ terms vanish identically rather than only in expectation — so observed information equals expected information, Newton equals Fisher scoring, and the Hessian $-\frac{1}{\varphi }{X}^{\top }WX⪯0$ makes $\ell$ concave. $\square$

Bridge. The proof's engine is that the score factorises through the chain $\theta \to \mu \to \eta \to \beta$ and the two link-and-variance factors collapse into the single weight $W_i=w_i/\{V({\mu }_i)g^{\prime }({\mu }_i{)}^2\}$, so the Fisher information is the weighted Gram matrix $X^{\top }WX$ of the linear model 45.06.01 and each scoring step is the generalized-least-squares solve of that unit applied to the working response $z$. This is exactly the IRLS identity of 45.08.03, now freed of the Bernoulli specifics: the foundational reason the logistic update reduced to clean weighted least squares is that the logit is the canonical link of the binomial, and the same canonical-link mechanism — observed information equal to expected information — operates for every family, which is the central insight that one algorithm fits them all. The construction generalises the exponential-family moment identities of 45.01.02, where $b^{\prime }(\theta )=\mathbb{E}[Y]$ and $b^{\prime \prime }(\theta )=\mathrm{Var}$ are the cumulant derivatives that here become $\mu$ and $V(\mu )$; putting these together, the GLM is the exponential family of that unit driven by a linear predictor, and the score equation $X^{\top }W(z-\eta )=0$ is its likelihood equation. The result builds toward the deviance and the analysis of deviance below, where the converged fit is compared against the saturated model, and it appears again in the asymptotic theory of 45.04.03, where $(X^{\top }WX{)}^{-1}$ at the optimum is the asymptotic covariance of $\hat{\beta }$.

Exercises Intermediate+

Lean formalization Intermediate+

Mathlib does not host the GLM layer, so no module is wired in (lean_status: none). The intended statements — the mean-variance identities from the cumulant function and the IRLS update from Fisher scoring — would read roughly as follows.

-- Intended shape; not part of the current Babel Bible Lean build.
-- Requires: an exponential-dispersion family with cumulant function b,
-- a design matrix X, a link g, and the matrix/positive-definite API.
variable {n p : ℕ} (X : Matrix (Fin n) (Fin p) ℝ)

-- cumulant function b, canonical parameter θ, mean μ = b'(θ)
noncomputable def mean (b : ℝ → ℝ) (θ : ℝ) : ℝ := deriv b θ
noncomputable def varFun (b : ℝ → ℝ) (θ : ℝ) : ℝ := deriv (deriv b) θ

-- mean-variance identity: E[Y] = b'(θ), Var(Y) = a(φ) b''(θ)
theorem mean_eq_b' (b : ℝ → ℝ) (θ : ℝ) :
    mean b θ = deriv b θ := rfl

-- Fisher information of a GLM is the weighted Gram matrix XᵀWX
noncomputable def fisherInfo (W : Matrix (Fin n) (Fin n) ℝ) :
    Matrix (Fin p) (Fin p) ℝ := Xᵀ * W * X

-- canonical-link concavity: -XᵀWX ⪯ 0 with W ⪰ 0 diagonal
theorem glm_concave_canonical (W : Matrix (Fin n) (Fin n) ℝ)
    (hW : W.PosSemidef) : (fisherInfo X W).PosSemidef := by
  sorry

the Mathlib gap analysis records what is missing: the exponential-dispersion family object, the mean-variance identities $\mu =b^{\prime }(\theta )$ and $\mathrm{Var}=a(\varphi )b^{\prime \prime }(\theta )$, the variance function, the canonical link $g=(b^{\prime }{)}^{-1}$, the IRLS update ${\beta }^{\mathrm{new}}=(X^{\top }WX{)}^{-1}{X}^{\top }Wz$, the canonical-link concavity, and the deviance with its asymptotic distribution.

Advanced results Master

The framework's coherence rests on the exponential-dispersion structure, and the principal results follow from the cumulant function $b$ and its derivatives.

Theorem 1 (mean, variance, and the variance function from the cumulant). For the exponential-dispersion density $f(y;\theta ,\varphi )=\exp \{(y\theta -b(\theta ))/a(\varphi )+c(y,\varphi )\}$ with $a(\varphi )>0$ and $b$ strictly convex on the interior of the canonical parameter set, the mean is $\mu =b^{\prime }(\theta )$ and the variance is $\mathrm{Var}(Y)=a(\varphi )b^{\prime \prime }(\theta )=a(\varphi )V(\mu )$ with $V(\mu )=b^{\prime \prime }((b^{\prime }{)}^{-1}(\mu ))$. The map $\theta \mapsto \mu =b^{\prime }(\theta )$ is a strictly increasing bijection of the interior onto the mean space, so each family is equivalently parametrised by $\theta$ or by $\mu$. These are the cumulant identities of 45.01.02 specialised to the one-parameter mean-indexed family carrying a dispersion.

Theorem 2 (Fisher scoring is IRLS; canonical-link concavity). As established in the Key theorem, the GLM score is $\nabla \ell =\frac{1}{\varphi }{X}^{\top }W(z-\eta )$ with $W=\mathrm{diag}(w_i/\{V({\mu }_i)g^{\prime }({\mu }_i{)}^2\})$ and working response $z=\eta +(y-\mu )\odot g^{\prime }(\mu )$, the Fisher information is $\frac{1}{\varphi }{X}^{\top }WX$, and the scoring step is the weighted least-squares fit $(X^{\top }WX{)}^{-1}{X}^{\top }Wz$. For a canonical link the score is $X^{\top }(y-\mu )/\varphi$, the observed and expected information agree, and the Hessian $-\frac{1}{\varphi }{X}^{\top }WX$ is negative semidefinite, so the log-likelihood is concave and (under full column rank with $W\succ 0$) strictly concave, giving a unique maximiser when one exists. Non-canonical links retain consistency but lose the observed-equals-expected coincidence, so Fisher scoring and Newton-Raphson then differ.

Theorem 3 (the deviance and the analysis of deviance). Define the deviance $D=2\{\ell ({\hat{\mu }}^{\text{sat}};y)-\ell (\hat{\mu };y)\}a(\varphi )={\sum }_i{d}_i$, where the saturated model fits ${\hat{\mu }}_i^{\text{sat}}=y_i$ and $d_i$ is the unit deviance; the scaled deviance is $D^{\ast }=D/\varphi$. For two nested GLMs $M_0\subset M_1$ with $p_0<p_1$ parameters, the scaled-deviance difference $D_0^{\ast }-D_1^{\ast }=2\{\ell ({\hat{M}}_1)-\ell ({\hat{M}}_0)\}$ equals the likelihood-ratio statistic for the nested comparison (the saturated terms cancel), so under $M_0$ and regular conditions it is asymptotically ${\chi }_{p_1-p_0}^2$ by Wilks' theorem 45.04.09. Sequencing such comparisons over a chain of nested models is the analysis of deviance, the GLM analogue of the analysis of variance, which it recovers exactly in the Gaussian-identity case where $d_i=(y_i-{\hat{\mu }}_i{)}^2$ and $D$ is the residual sum of squares.

Theorem 4 (quasi-likelihood and overdispersion). Suppose only the mean-variance relationship $\mathbb{E}[Y_i]={\mu }_i(\beta )$, $\mathrm{Var}(Y_i)=\varphi V({\mu }_i)$ is specified, with no full distribution. The quasi-score $U(\beta )={\sum }_i\{w_i(y_i-{\mu }_i)/(\varphi V({\mu }_i)g^{\prime }({\mu }_i))\}{x}_i=0$ coincides with the genuine GLM score ^{[Wedderburn 1974]}, so the quasi-likelihood estimator ${\hat{\beta }}_Q$ solves the same estimating equation and is consistent and asymptotically normal with covariance $\varphi (X^{\top }WX{)}^{-1}$ whenever the mean and variance models are correct, regardless of higher moments. The dispersion is estimated by the Pearson statistic $\hat{\varphi }=(n-p{)}^{-1}{\sum }_i(y_i-{\hat{\mu }}_i{)}^2/V({\hat{\mu }}_i)$; when the data are overdispersed ($\hat{\varphi }>1$ for nominal Poisson or binomial) this inflates the naive standard errors, restoring honest inference that the nominal $\varphi =1$ would understate.

Theorem 5 (asymptotic normality of $\hat{\beta }$). Under the regularity conditions of maximum-likelihood asymptotics 45.04.03 — interior true parameter, full column rank with bounded covariates, finite and non-singular limiting information $\bar{\mathcal{I}}=\lim n^{-1}\mathcal{I}({\beta }_0)$ — the GLM maximum-likelihood estimator satisfies $\sqrt{n}(\hat{\beta }-{\beta }_0)\to \mathcal{N}(0,\varphi {\bar{\mathcal{I}}}^{-1})$, with the plug-in estimate $\widehat{\mathrm{Var}}(\hat{\beta })=\hat{\varphi }(X^{\top }\hat{W}X{)}^{-1}$. Wald, score, and likelihood-ratio tests for $\beta$ then follow the standard trinity of 45.04.09, the score test using the information $X^{\top }WX$ evaluated under the null and the likelihood-ratio test reading off the analysis of deviance.

Synthesis. The foundational reason the whole framework coheres is the exponential-dispersion family: the cumulant function $b$ supplies the mean $b^{\prime }(\theta )$, the variance function $b^{\prime \prime }$, and — through the canonical link $g=(b^{\prime }{)}^{-1}$ — the identification $\theta =\eta$ that makes the natural statistic $X^{\top }y$ sufficient for $\beta$, so the GLM is exactly the exponential family of 45.01.02 driven by a linear predictor. This is exactly why one algorithm fits them all: the chain rule collapses the link and the variance function into the single weight $W=\mathrm{diag}(w_i/\{V({\mu }_i)g^{\prime }({\mu }_i{)}^2\})$, and Fisher scoring becomes the weighted least squares of the linear model 45.06.01, generalising the logistic IRLS of 45.08.03 from the Bernoulli variance $\mu (1-\mu )$ to every variance function at once. The deviance is dual to the likelihood ratio — the scaled-deviance difference is the nested likelihood-ratio statistic, so the analysis of deviance and Wilks' theorem 45.04.09 are one object seen twice, recovering the analysis of variance in the Gaussian corner. Putting these together, the central insight is that quasi-likelihood needs even less than a family: because the score depends on the response only through $\mu$ and $V(\mu )$, postulating just the mean-variance relationship suffices for consistent estimation, which is the bridge from full-likelihood GLMs to the estimating-equation viewpoint and the practical remedy for overdispersion, while the asymptotic covariance $\varphi (X^{\top }WX{)}^{-1}$ ties the converged fit back to the MLE asymptotics of 45.04.03.

Full proof set Master

The Fisher-scoring/IRLS identity (Theorem 2) and the deviance/likelihood-ratio equality (Theorem 3, Exercise 7) are proved above. The remaining Master claims are recorded here.

Proposition 1 (mean and variance from the cumulant function). For $f(y;\theta ,\varphi )=\exp \{(y\theta -b(\theta ))/a(\varphi )+c(y,\varphi )\}$ with $b$ twice differentiable on the interior of the canonical set, $\mathbb{E}[Y]=b^{\prime }(\theta )$ and $\mathrm{Var}(Y)=a(\varphi )b^{\prime \prime }(\theta )$.

Proof. The score in $\theta$ of a single observation is ${\partial }_{\theta }\log f=(y-b^{\prime }(\theta ))/a(\varphi )$. A regular likelihood has zero-mean score, so $\mathbb{E}[(Y-b^{\prime }(\theta ))/a(\varphi )]=0$, giving $\mathbb{E}[Y]=b^{\prime }(\theta )=\mu$. The second log-derivative is ${\partial }_{\theta }^2\log f=-b^{\prime \prime }(\theta )/a(\varphi )$, a constant in $y$. The information identity $\mathrm{Var}({\partial }_{\theta }\log f)=-\mathbb{E}[{\partial }_{\theta }^2\log f]$ reads $\mathrm{Var}(Y)/a(\varphi {)}^2=b^{\prime \prime }(\theta )/a(\varphi )$, so $\mathrm{Var}(Y)=a(\varphi )b^{\prime \prime }(\theta )$. Writing $V(\mu )=b^{\prime \prime }(\theta (\mu ))$ for the inverse map $\theta (\mu )=(b^{\prime }{)}^{-1}(\mu )$, which is well-defined and increasing because $b^{\prime \prime }>0$ by the strict convexity of $b$, gives $\mathrm{Var}(Y)=a(\varphi )V(\mu )$. $\square$

Proposition 2 (the canonical link makes the likelihood concave). With a canonical link, $g=(b^{\prime }{)}^{-1}$, the GLM log-likelihood $\ell (\beta )=\frac{1}{\varphi }{\sum }_i{w}_i\{y_i{\eta }_i-b({\eta }_i)\}+\text{const}$ has gradient $\frac{1}{\varphi }{X}^{\top }{W}_0(y-\mu )$ and Hessian $-\frac{1}{\varphi }{X}^{\top }WX⪯0$, with $W=\mathrm{diag}(w_i{b}^{\prime \prime }({\eta }_i))=\mathrm{diag}(w_{i}V({\mu }_i))$; under full column rank of $X$ and $W\succ 0$ the Hessian is negative definite and $\ell$ is strictly concave.

Proof. For the canonical link ${\theta }_i={\eta }_i=x_i^{\top }\beta$, so ${\ell }_i=\frac{w_i}{\varphi }\{y_i{\eta }_i-b({\eta }_i)\}+c$. Then ${\partial }_{\beta }{\ell }_i=\frac{w_i}{\varphi }\{y_i-b^{\prime }({\eta }_i)\}{x}_i=\frac{w_i}{\varphi }(y_i-{\mu }_i)x_i$ since ${\mu }_i=b^{\prime }({\eta }_i)$, and summing gives $\nabla \ell =\frac{1}{\varphi }{X}^{\top }\mathrm{diag}(w_i)(y-\mu )$. Differentiating again, ${\partial }_{\beta {\beta }^{\top }}^2{\ell }_i=-\frac{w_i}{\varphi }{b}^{\prime \prime }({\eta }_i)x_i{x}_i^{\top }$, which is non-random because $b^{\prime \prime }({\eta }_i)$ does not involve $y_i$; summing gives ${\nabla }^2\ell =-\frac{1}{\varphi }{X}^{\top }WX$ with $W=\mathrm{diag}(w_i{b}^{\prime \prime }({\eta }_i))$ and $b^{\prime \prime }({\eta }_i)=V({\mu }_i)>0$. For any $v$, $v^{\top }{X}^{\top }WXv={\sum }_i{w}_{i}V({\mu }_i)(x_i^{\top }v{)}^2\ge 0$, so $\ell$ is concave; equality forces $x_i^{\top }v=0$ for all $i$, i.e. $Xv=0$, which under full column rank forces $v=0$. Hence $X^{\top }WX\succ 0$ and $-\ell$ is strictly convex. $\square$

Proposition 3 (Fisher scoring = Newton-Raphson for the canonical link). For a canonical link the observed information $-{\nabla }^2\ell$ equals the expected information $\mathcal{I}(\beta )$, so the Fisher-scoring and Newton-Raphson updates coincide.

Proof. Newton-Raphson uses the observed Hessian ${\nabla }^2\ell$; Fisher scoring uses $-\mathcal{I}(\beta )=\mathbb{E}[{\nabla }^2\ell ]$. By Proposition 2 the canonical-link Hessian $-\frac{1}{\varphi }{X}^{\top }WX$ with $W=\mathrm{diag}(w_i{b}^{\prime \prime }({\eta }_i))$ contains no $y_i$, so it equals its own expectation; thus ${\nabla }^2\ell =\mathbb{E}[{\nabla }^2\ell ]$ identically and the two updates are the same map $\beta \mapsto \beta +(X^{\top }WX{)}^{-1}{X}^{\top }\mathrm{diag}(w_i)(y-\mu )$. For a non-canonical link the omitted $(y_i-{\mu }_i)$ terms in the Hessian have nonzero realisation though zero mean, so observed and expected information differ and the two iterations diverge in their steps (though they share the same fixed point). $\square$

Proposition 4 (consistency of the quasi-likelihood estimator). If $\mathbb{E}[Y_i]={\mu }_i({\beta }_0)$ and the variance model $\mathrm{Var}(Y_i)=\varphi V({\mu }_i)$ holds, the quasi-score $U(\beta )={\sum }_i{w}_i(y_i-{\mu }_i)x_i/\{V({\mu }_i)g^{\prime }({\mu }_i)\}$ is an unbiased estimating function, and its root ${\hat{\beta }}_Q$ is consistent for ${\beta }_0$ with asymptotic covariance $\varphi (X^{\top }WX{)}^{-1}$.

Proof. Unbiasedness: ${\mathbb{E}}_{{\beta }_0}[U({\beta }_0)]={\sum }_i{w}_i\mathbb{E}[y_i-{\mu }_i({\beta }_0)]x_i/\{V{g}^{\prime }\}=0$ since $\mathbb{E}[y_i]={\mu }_i({\beta }_0)$, using nothing beyond the mean model. An unbiased estimating function with a non-singular sensitivity matrix $-\mathbb{E}[\partial U/\partial {\beta }^{\top }]=\frac{1}{\varphi }{X}^{\top }WX$ yields a consistent root by the standard estimating-equation argument 45.04.03: a Taylor expansion of $U({\hat{\beta }}_Q)=0$ about ${\beta }_0$ gives $\sqrt{n}({\hat{\beta }}_Q-{\beta }_0)=(\frac{1}{n}{X}^{\top }WX{)}^{-1}\frac{1}{\sqrt{n}}U({\beta }_0)+o_P(1)$. The variance of the quasi-score is $\mathrm{Var}(U({\beta }_0))={\sum }_i{w}_i^2\mathrm{Var}(y_i)x_i{x}_i^{\top }/\{V{g}^{\prime }{\}}^2=\varphi X^{\top }WX$ using $\mathrm{Var}(y_i)=\varphi V({\mu }_i)$, so by the central limit theorem $\frac{1}{\sqrt{n}}U({\beta }_0)\to \mathcal{N}(0,\varphi \bar{W})$ and $\sqrt{n}({\hat{\beta }}_Q-{\beta }_0)\to \mathcal{N}(0,\varphi {\bar{W}}^{-1})$ with $\bar{W}=\lim n^{-1}{X}^{\top }WX$. Higher moments of $Y$ never enter. $\square$

Connections Master

• Logistic regression and IRLS 45.08.03 is the canonical example this unit abstracts: the Bernoulli GLM with the logit link has variance function $V(\mu )=\mu (1-\mu )$ and weight $W_i={\mu }_i(1-{\mu }_i)$, and its update $(X^{\top }WX{)}^{-1}{X}^{\top }Wz$ is the general GLM IRLS step with that single variance function substituted, so everything proved there is the binomial line of the framework here.

• Sufficiency and exponential families 45.01.02 is the random component: the exponential-dispersion density is the one-parameter exponential family of that unit indexed in its mean and carrying a dispersion, the mean-variance identities $\mu =b^{\prime }(\theta )$ and $V(\mu )=b^{\prime \prime }(\theta )$ are its cumulant identities, and the canonical link $g=(b^{\prime }{)}^{-1}$ is exactly the natural-parameter map that makes $X^{\top }y$ sufficient for $\beta$.

• The linear model and least squares 45.06.01 is reused inside every scoring step: the GLM Fisher information $X^{\top }WX$ is the weighted Gram matrix of that unit, and IRLS is a loop of its generalized-least-squares solves applied to the working response, recovering ordinary least squares exactly in the Gaussian-identity corner where $W=I$ and $z=y$.

• The likelihood-ratio, Wald, and score asymptotics 45.04.09 supply the inference: the scaled-deviance difference between nested GLMs is the likelihood-ratio statistic, so the analysis of deviance is Wilks' theorem applied along a chain of models, and the score test uses the information $X^{\top }WX$ under the null.

• MLE asymptotic normality and efficiency 45.04.03 governs the converged estimator: under interior-parameter, full-rank, and finite-information regularity the GLM MLE is $\sqrt{n}$-consistent and asymptotically normal with covariance $\varphi (X^{\top }WX{)}^{-1}$, the inverse Fisher information of that unit, and the quasi-likelihood estimator inherits the same sandwich from the estimating-equation theory.

Historical & philosophical context Master

The generalized linear model was named and unified by John Nelder and Robert Wedderburn in their 1972 paper Generalized linear models (J. R. Statist. Soc. A 135) ^{[Nelder 1972]}, which recognised that ordinary regression, probit and logit analysis of binary data, log-linear models for contingency tables and counts, and gamma models for positive continuous responses are instances of a single class: an exponential-family response, a linear predictor, and a link function joining the mean to that predictor. Their decisive contribution was algorithmic as well as conceptual — a common maximum-likelihood fitting procedure, Fisher scoring expressed as iteratively reweighted least squares, that fits every member with one routine, implemented in the GLIM software that carried the framework into practice.

The statistical antecedents were older and scattered. Fisher's scoring method and the exponential-family theory of sufficiency 45.01.02 supplied the likelihood machinery; probit analysis (Bliss, 1934) and the logit (Berkson, 1944) were the binary-response special cases; log-linear count models grew from contingency-table analysis. Robert Wedderburn's 1974 paper Quasi-likelihood functions, generalized linear models, and the Gauss-Newton method (Biometrika 61) ^{[Wedderburn 1974]} then loosened the requirement of a full distribution, observing that the scoring equations depend on the response only through its mean and variance function, so that postulating the mean-variance relationship alone yields consistent estimators and a principled treatment of overdispersion. The definitive synthesis is McCullagh and Nelder's monograph Generalized Linear Models (1989), whose treatment of the exponential-dispersion family, the deviance, and the analysis of deviance remains the standard reference; the quasi-likelihood line was later subsumed into the broader theory of generalized estimating equations by Liang and Zeger (1986) for correlated data.

Bibliography Master

@book{mccullagh1989generalized,
  author    = {McCullagh, Peter and Nelder, John A.},
  title     = {Generalized Linear Models},
  edition   = {2},
  publisher = {Chapman and Hall},
  year      = {1989}
}

@article{nelder1972generalized,
  author  = {Nelder, John A. and Wedderburn, Robert W. M.},
  title   = {Generalized linear models},
  journal = {Journal of the Royal Statistical Society, Series A},
  volume  = {135},
  number  = {3},
  pages   = {370--384},
  year    = {1972}
}

@article{wedderburn1974quasi,
  author  = {Wedderburn, Robert W. M.},
  title   = {Quasi-likelihood functions, generalized linear models, and the Gauss-Newton method},
  journal = {Biometrika},
  volume  = {61},
  number  = {3},
  pages   = {439--447},
  year    = {1974}
}

@book{hastie2009elements,
  author    = {Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome},
  title     = {The Elements of Statistical Learning: Data Mining, Inference, and Prediction},
  edition   = {2},
  publisher = {Springer},
  year      = {2009},
  note      = {Chapter 4, linear methods for classification}
}

@article{berkson1944application,
  author  = {Berkson, Joseph},
  title   = {Application of the logistic function to bio-assay},
  journal = {Journal of the American Statistical Association},
  volume  = {39},
  number  = {227},
  pages   = {357--365},
  year    = {1944}
}

@article{liang1986longitudinal,
  author  = {Liang, Kung-Yee and Zeger, Scott L.},
  title   = {Longitudinal data analysis using generalized linear models},
  journal = {Biometrika},
  volume  = {73},
  number  = {1},
  pages   = {13--22},
  year    = {1986}
}