We have now covered all the core elements of generalised linear models. This chapter consolidates those ideas by revisiting the three-component structure of a GLM, summarising the key results from earlier chapters, and working through a complete analysis from start to finish. We then look forward to Part II, where we will relax the linearity assumption on the systematic component.
Every generalised linear model has exactly three components:
Random component. The response variable \(Y_i\) follows a distribution from the exponential family: \[f(y_i;\theta_i,\phi) = \exp\left\{\frac{y_i\theta_i - b(\theta_i)}{\phi} + c(y_i,\phi)\right\}.\] The mean and variance are determined by \(b(\theta_i)\): \[\mathrm{E}[Y_i] = b'(\theta_i) = \mu_i, \qquad \mathrm{Var}[Y_i] = \phi\, b''(\theta_i).\]
Systematic component. A linear predictor in the explanatory variables: \[\eta_i = \mathbf{x}_i^T\boldsymbol{\beta} = \sum_{j=1}^p \beta_j x_{ij}.\]
Link function. A monotone differentiable function \(g\) connecting mean to linear predictor: \[g(\mu_i) = \eta_i.\] When \(g = (b')^{-1}\), the link is canonical and simplifies the score equations.
The following table summarises the choices for the three distributions we have studied in detail.
| Distribution | Canonical link | Variance function | Typical use |
|---|---|---|---|
| Normal | Identity: \(\eta = \mu\) | \(b''(\theta)=1\) | Continuous response |
| Binomial | Logit: \(\eta = \log\frac{p}{1-p}\) | \(b''(\theta)=mp(1-p)\) | Counts / proportions |
| Poisson | Log: \(\eta = \log\lambda\) | \(b''(\theta)=\lambda\) | Count data |
Estimation. Parameters are estimated by maximum likelihood using the iteratively reweighted least squares (IRLS) algorithm (see the supplementary notes proof_IRLS.pdf for the full derivation). The MLE satisfies asymptotically: \[\hat{\boldsymbol\beta} \sim N_p\!\left(\boldsymbol\beta,\, \mathcal{J}^{-1}(\boldsymbol\beta)\right).\]
Deviance. The deviance measures how much a fitted model \(M\) falls short of the saturated model: \[D = 2\phi\left\{l(\tilde{\boldsymbol\theta};\mathbf{y},\phi) - l(\hat{\boldsymbol\theta};\mathbf{y},\phi)\right\}.\] For the three standard families, \(D\) reduces to:
Model comparison. For two nested models \(M_1 \subset M_2\) with \(r_1 < r_2\) parameters:
Residuals. Deviance residuals \(e_i^\text{dev} = \text{sign}(y_i-\hat\mu_i)\sqrt{d_i}\) and Pearson residuals \(e_i^P = (y_i-\hat\mu_i)/\sqrt{b''(\hat\theta_i)}\) are the main diagnostic tools.
This example illustrates the complete GLM workflow: data exploration, model formulation, fitting, inference, and model criticism.
The data (Dobson and Barnett, 2008, p.127) relate to a study of coronary heart disease (CHD). A total of 1,329 men were cross-classified by cholesterol level (4 ordered groups: 1 = lowest, 4 = highest) and blood pressure (4 ordered groups: 1 = lowest, 4 = highest). For each of the \(4 \times 4 = 16\) cells, the number of men with CHD (\(y_{ij}\)) and the total number of men (\(m_{ij}\)) are recorded. The data are in https://www.richardpmann.com/MATH3701/Datasets/chd.txt.
chd = read.table("https://www.richardpmann.com/MATH3701/Datasets/chd.txt", header=TRUE)
chd y m chol bp
1 2 119 1 1
2 3 124 1 2
3 3 50 1 3
4 4 26 1 4
5 3 88 2 1
6 2 100 2 2
7 0 43 2 3
8 3 23 2 4
9 8 127 3 1
10 11 220 3 2
11 6 74 3 3
12 6 49 3 4
13 7 74 4 1
14 12 111 4 2
15 11 57 4 3
16 11 44 4 4The proportion with CHD in each group:
chd$prop = chd$y / chd$m
par(mar=c(4,4,1,1))
plot(chd$bp, chd$prop,
type="n",
xlab="Blood pressure group", ylab="Proportion with CHD",
xlim=c(0.8, 4.2), ylim=c(0, 0.35))
cols = c("blue","darkgreen","orange","red")
for (k in 1:4) {
sub = chd[chd$chol == k, ]
lines(sub$bp, sub$prop, col=cols[k], pch=16, type="b")
}
legend("topleft", legend=paste("Chol group", 1:4),
col=cols, lty=1, pch=16, bty="n")
Both increasing cholesterol and increasing blood pressure appear associated with higher CHD proportions.
Each \(Y_{ij}\) is a count of CHD cases out of \(m_{ij}\) men, so a Binomial GLM with logit link is appropriate: \[Y_{ij} \sim \text{B}(m_{ij}, p_{ij}), \qquad \text{logit}(p_{ij}) = \mu + \alpha_i + \beta_j,\] where \(\alpha_i\) is the effect of cholesterol group \(i\) and \(\beta_j\) is the effect of blood pressure group \(j\). We will also consider a saturated model (with an interaction term) and the null model.
Since both factors are ordered, we will first treat them as factors (to allow flexible effect estimates) and then consider whether a simpler model with numeric scores is adequate.
chd$cholF = as.factor(chd$chol)
chd$bpF = as.factor(chd$bp)
# Response matrix: successes and failures
ym = cbind(chd$y, chd$m - chd$y)
# Null model
glm0 = glm(ym ~ 1, family=binomial, data=chd)
# Main effects model (independence)
glm1 = glm(ym ~ cholF + bpF, family=binomial, data=chd)
# Saturated model (with interaction)
glm2 = glm(ym ~ cholF * bpF, family=binomial, data=chd)
# Summary of deviances
cat("Null deviance: ", round(glm0$deviance, 2),
"on", glm0$df.residual, "df\n")Null deviance: 58.73 on 15 dfcat("Main effects deviance:", round(glm1$deviance, 2),
"on", glm1$df.residual, "df\n")Main effects deviance: 8.08 on 9 dfcat("Saturated deviance: ", round(glm2$deviance, 2),
"on", glm2$df.residual, "df\n")Saturated deviance: 0 on 0 dfTest main effects vs null:
pchisq(glm0$deviance - glm1$deviance,
glm0$df.residual - glm1$df.residual,
lower.tail=FALSE)[1] 3.48186e-09Both cholesterol and blood pressure are highly significant.
Test interaction vs main effects:
pchisq(glm1$deviance - glm2$deviance,
glm1$df.residual - glm2$df.residual,
lower.tail=FALSE)[1] 0.5264905The interaction is not significant, so the main-effects model is adequate.
Goodness of fit for main-effects model:
pchisq(glm1$deviance, glm1$df.residual, lower.tail=FALSE)[1] 0.5264905The main-effects model fits the data well (\(p > 0.05\)).
summary(glm1)
Call:
glm(formula = ym ~ cholF + bpF, family = binomial, data = chd)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.48194 0.34865 -9.987 < 2e-16 ***
cholF2 -0.20798 0.46641 -0.446 0.65566
cholF3 0.56223 0.35080 1.603 0.10900
cholF4 1.34412 0.34297 3.919 8.89e-05 ***
bpF2 -0.04146 0.30365 -0.137 0.89139
bpF3 0.53236 0.33240 1.602 0.10925
bpF4 1.20042 0.32689 3.672 0.00024 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 58.7262 on 15 degrees of freedom
Residual deviance: 8.0762 on 9 degrees of freedom
AIC: 73.18
Number of Fisher Scoring iterations: 4The estimated log-odds ratios for cholesterol and blood pressure groups (relative to group 1 in each case) show a clear increasing trend: higher cholesterol and higher blood pressure each increase the log-odds of CHD.
par(mar=c(4,4,1,1))
plot(fitted(glm1), residuals(glm1, type="deviance"),
pch=16,
xlab="Fitted probability", ylab="Deviance residual")
abline(h=0, lty=2)
No strong pattern is apparent, consistent with a well-fitting model.
The GLM framework is powerful, but the linearity constraint \(\eta_i = \mathbf{x}_i^T\boldsymbol{\beta}\) can be restrictive.
Consider a situation where the true relationship between a continuous explanatory variable \(x\) and the response is smooth but nonlinear. We have two broad options:
Polynomial regression: add \(x^2, x^3, \dots\) to the linear predictor. This is easy but can produce unstable fits and the polynomial degree must be chosen in advance.
Non-parametric smoothing: replace \(\beta x\) with a smooth function \(f(x)\) estimated from the data. This forms the basis of Generalised Additive Models (GAMs), studied in Part II.
The GAM model is written: \[g(\mu_i) = \alpha + f_1(x_{i1}) + f_2(x_{i2}) + \cdots + f_p(x_{ip}),\] where each \(f_j\) is a smooth function (a spline) chosen to fit the data while penalising excessive wiggliness.
We will see that GAMs retain the interpretability of GLMs while gaining considerable flexibility, and that the estimation and inference methods are natural generalisations of what we have developed in Part I.
Key points from Part I