ST346 Lecture Notes

Contents

Preface 4

Introduction 5
Assumptions of the normal linear model . . . . . . . . . . . . . . . . . . . . . . . . . 5
Relaxing the assumptions for GLMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
Why learn about GLMs? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1 Weighted linear regression 6

1.1 The Gauss-Markov theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Weighted least squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3 The hat matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.4 Weighted residuals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2 Binomial models 15
2.1 Bernoulli distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 Latent linear model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3 Binomial outcomes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4 Scaled binomial distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.5 Implementation in R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.6 Logistic regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.7 Tolerance distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3 Poisson regression models 37

3.1 Moments and cumulants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2 The Poisson distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.3 Poisson regression for counts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.4 Poisson regression for rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

4 Exponential dispersion models 50

4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.2 Definition of an EDM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.3 Weighted EDMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.4 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.5 Cumulants for EDMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.6 The canonical link . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.7 Examples of canonical link functions . . . . . . . . . . . . . . . . . . . . . . . . 57
4.8 Deviance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5 Generalized linear models 64

5.1 The variance function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.2 Quasi-likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

2
CONTENTS 3

5.3 The gamma distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

5.4 The inverse Gaussian distribution . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.5 The Poisson limit of the binomial distribution . . . . . . . . . . . . . . . . . . . 71
5.6 Choosing between the Poisson and Binomial GLM . . . . . . . . . . . . . . . . . 72

6 Maximum likelihood estimation for GLMs 75

6.1 Review of maximum likelihood estimation . . . . . . . . . . . . . . . . . . . . . 75
6.2 Maximum likelihood for GLMs . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
6.3 Numerical maximum likelihood estimation for GLMs . . . . . . . . . . . . . . . 85
6.4 Estimating the dispersion parameter . . . . . . . . . . . . . . . . . . . . . . . . 89

7 Diagnostics for GLMs 92

7.1 Residuals for GLMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
7.2 Leverage and standardized residuals . . . . . . . . . . . . . . . . . . . . . . . . . 95
7.3 Influence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
7.4 Diagnostic plots in R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

8 Model choice and Hypothesis tests 104

8.1 Confidence intervals based on asymptotic normality . . . . . . . . . . . . . . . . 104
8.2 Confidence intervals based on profile likelihood . . . . . . . . . . . . . . . . . . . 105
8.3 Hypothesis tests for GLMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
8.4 The summary function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
8.5 Model choice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
8.6 Variable selection - Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

9 Overdispersion, prediction and classification 120

9.1 Overdispersion for count data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
9.2 Introduction to prediction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
9.3 Scoring rules for probabilistic predictions . . . . . . . . . . . . . . . . . . . . . . 127
9.4 Prediction versus Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
9.5 ROC curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135

10 Further topics 138

10.1 The saddlepoint approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
10.2 Applications of the saddlepoint approximation . . . . . . . . . . . . . . . . . . . 141
10.3 Quasi-likelihood for GLMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142

Some Revision Questions 144

Preface

These notes are an adaptation from the original ST346 slides developed by Prof Martyn Plummer
for the module.
If you find any typos, please inform the module leader.
These materials are solely for your own use and you must not distribute these in any format.
Do not upload these materials to the internet or any filesharing sites nor provide
them to any third party or forum.
All rights are reserved.

4
Introduction

Assumptions of the normal linear model

Recall the assumptions for a normal linear model as introduced in ST231 Linear Statistical
Modelling in R.

Y ∼ Nn (µ, σ 2 I n )
where for all i = 1, . . . , n:
• E(Yi | xi ) = µi with µi = xTi β;
• Var(Yi | xi ) = σ 2 ;
• Cov(Yi , Yj | xi , xj ) = 0 for i ̸= j;
• Yi | xi ∼ N (µi , σ 2 ) and so Yi is continuous and unbounded.

Relaxing the assumptions for GLMs

Generalized linear models (GLMs) extend from the linear models by relaxing the assumptions
as follows:
1. The response distribution is no longer necessarily normal, but a member of the exponential
family of distributions (which includes bounded and/or discrete distributions).
2. A link function g allows for a non-linear relationship between the parameter vector β
and the mean µi :
g(µi ) = xTi β.

3. A variance function V (µ) allows the variance of the response to depend on the mean:
Var(Yi ) = ϕV (µi ).

Note:
a. We recover the normal linear model by using the normal distribution as the response
distribution and setting g(µ) = µ, V (µ) = 1 and ϕ = σ 2 .
b. Point 3. is a consequence of Point 1., that is by choosing a specific distribution from the
exponential family we determine the relationship between the mean and the variance of
the response.

Why learn about GLMs?

• GLMs greatly expand the range of problems accessible to regression modelling.
• Many methods from normal linear models and its extensions can be generalised to GLMs.

5
Chapter 1 Weighted linear regression

Before we discuss GLMs, we introduce weighted linear models. These come into play when we
discuss how to produce maximum likelihood estimates for the parameters of a GLM.
First we consider the Gauss-Markov theorem which shows that the (ordinary) least squares
estimator for the parameters of a linear model is, in some sense, optimal.
We then consider a linear model where the variance is not constant but known up to a constant of
proportionality. We show that in this case, we can still estimate the parameter vector. However
to satisfy the same optimality criterion as in the standard case we need to adapt our estimator
to the generalised least squares estimator.

1.1 The Gauss-Markov theorem

Consider the linear model

Yi = xTi β + ϵi for i = 1, . . . , n

where the errors ϵ1 , . . . , ϵn are uncorrelated, have mean zero and variance σ 2 .
The (ordinary) least squares estimator is
 −1
β
b = XT X XT Y

and satisfies

E(β)
b = β;
 −1
Var(β)
b = σ2 X T X .

Suppose now that the scalar parameter γ = aT β is the focus of interest.

For example, if δik denotes Kronecker’s delta, that is

1 if i = k,
(
δik =
0 if i =
̸ k,

then we may set ai = δik for given k so that γ = aT β = βk selects a single element of the
parameter vector.

Theorem 1.1 (Gauss-Markov theorem). Let β b be the least squares estimator for β. Then
γb = aT β
b is the unique unbiased linear estimator of γ = aT β with minimum variance.

6
1.1. THE GAUSS-MARKOV THEOREM 7

Proof of the Gauss-Markov theorem (not examinable)

Our proof of the Gauss-Markov theorem is divided into four steps.
1. We show that γb = aT βb is a linear unbiased estimator for γ = aT β.
2. We then consider a rival estimator γ̃ for γ and investigate what conditions γ̃ needs to
satisfy to be unbiased.
3. Next, we prove that the variance of γ̃ is greater than or equal to the variance of γb .
4. Finally, we show that if the variance of γ̃ is the same as the variance of γb , then γ̃ = γb .
This proves the uniqueness of γb as the linear unbiased estimator of minimum variance.
Step 1
We start by proving that γ̂ belongs to the class of linear unbiased estimators.
We have
 −1  −1
γb = aT β
b = aT X T X XT Y = cT Y where c = X X T X a.

Because we can write γb as cT Y we have that γb is a linear combination of the response variables
Y1 , . . . , Yn and thus a linear estimator.
Next we show that γb is unbiased. As β
b is unbiased
   
E (γb ) = E aT β
b = aT E β
b = aT β = γ.

Hence E (γb ) = γ and so γb is unbiased.

In the last step we have shown that γ̂ is both a linear and an unbiased estimator. Our next
step is to introduce an alternative linear estimator γ̃, and show that such an estimator cannot
simultaneously be unbiased and have a smaller variance than γ̂. We start by identifying
conditions that the linear estimator γ̃, needs to satisfy in order to be unbiased.
Step 2
Let γ̃ be another linear estimator for γ. As it is a linear estimator, it must take the form
γ̃ = bT Y . Hence,
 
E (γ̃) = E bT Y = bT E (Y ) = bT Xβ

If γ̃ is unbiased, then this puts some constraints on the possible values of b, as we will show
next. For γ̃ to be unbiased we need E (γ̃) = γ = aT β and so bT Xβ = aT β for all values of
β. Hence
bT X = aT (1.1)

Therefore, if γ̃ is unbiased then b must satisfy (1.1).

8 CHAPTER 1. WEIGHTED LINEAR REGRESSION

Step 3
We now calculate the variances of γb and γ̃.
   
Var (γ̃) = Var bT Y = bT Var (Y ) b = bT σ 2 I n b = σ 2 bT b.

Similarly we can compute the variance of γb as

     −1  −1
Var (γb ) = Var aT β
b = aT Var β
b a = aT σ 2 X T X a = σ 2 aT X T X a.

Now, using property (1.1) we can replace aT with bT X and so obtain

 −1  −1
Var (γb ) = σ 2 aT X T X a = σ 2 bT X X T X XT b = σ 2 bT Hb
 −1
where H = X X T X X T is the so-called hat matrix.
Next we take the difference between the two variances and use the properties of the hat matrix
to simplify the formula
Var (γ̃) − Var (γb ) = σ 2 bT b − σ 2 bT Hb
= σ 2 bT (I n − H) b
= σ 2 bT (I n − H)T (I n − H) b (1.2)
where the last step uses the fact that (I n − H) is both symmetric and idempotent. Now, let
d = (I n − H) b then (1.2) becomes
Var (γ̃) − Var (γb ) = σ 2 dT d ⩾ 0
as dT d is a sum of squares. Therefore γb has minimum variance in the class of linear unbiased
estimators of γ.
The one thing still to prove is that there is no other linear unbiased estimator that has as small
a variance as γ̂. We therefore prove the uniqueness of γ̂ as regards to its variance.
Step 4
Suppose both estimators γ̃ and γb have the same variance. Then the difference in the variances is
Var (γ̃) − Var (γb ) = σ 2 dT d = 0.
As dT d is a sum of squares, this can only be true if di = 0 for all i = 1 . . . n. Therefore
dT = bT (I n − H) = 0
and so
bT = bT H (1.3)
We now use (1.1) and (1.3) to prove that γ̃ = γb .
γ̃ = bT Y = bT HY by (1.3)
= bT X β
b by the definition of H
= a βbT
by (1.1)
= γb
Hence if Var(γ̃) = Var(γb ) then γ̃ = γb .
1.2. WEIGHTED LEAST SQUARES 9

1.2 Weighted least squares

Let w1 , . . . , wn be known, positive weights and define the diagonal matrix W as

W = diag(w1 , . . . , wn ).

Now consider the linear model

Y = Xβ + ϵ, where ϵ ∼ Nn (0, ϕW −1 ).

Then

E(Yi | xi ) = xTi β, and

ϕ
Var(Yi | xi ) = ,
wi
where w1 , . . . , wn are known while the dispersion parameter ϕ may be unknown.
To determine the optimal estimator for β, we rescale our predictors xi and outcomes Yi . Let
√
Yi∗ = wi Yi ,
√
x∗i = wi xi .

Then

E(Yi∗ | x∗i ) = (x∗i )T β,

Var(Yi∗ | x∗i ) = ϕ,

and the problem is reduced to a homoscedastic linear model.

From the Gauss-Markov theorem, on the transformed scale (X ∗ , Y ∗ ), the optimal estimator for
β is OLS and thus minimizes the deviance (sum of squares)
X 2
D(β) = Yi∗ − (x∗i )T β
i

This can be written on the original scale (X, Y ) as

 2
= xTi β
X
D(β) wi Yi −
i

Hence, on the original scale, we minimize a weighted sum of squares.

The estimator β
b that minimizes the deviance D(β) can be written in closed form as

h i−1
β
b = (X ∗ )T X ∗ (X ∗ )T Y ∗

On the original scale (X, Y ) this can be rewritten as

h i−1
β
b = XT W X XT W Y

where W = diag(w1 , . . . wn ) is a diagonal matrix of weights.

10 CHAPTER 1. WEIGHTED LINEAR REGRESSION

Note
1. For a weighted linear model the weighted (or generalised) least squares estimator is the
unbiased linear estimator for β that has minimum variance. This generalisation of the
Gauss-Markov theorem is due to Aitken (1935).
2. Observations with a higher weight make a greater contribution to D(β). This is sensible as
an observation with higher weight has a smaller variance and so contains more information
about β. This explains why we divide the variance by the weight.

1.2.1 Example 1: weighted mean

Suppose Y1 , . . . , Yn are independent with common mean but different variances, that is

E(Yi ) = µ,
Var(Yi ) = ϕ/wi .

Then the optimal linear unbiased estimator is the weighted mean

Pn
wi Yi
µb = Pn
i=1

i=1 wi

where the weights are inversely proportional to the variance.

Proof:
This is a weighted linear model with an intercept term and no predictors.
Here X = 1n×1 and W = diag(w1 , . . . , wn ) and so
n
X T W X = 11×n W 1n×1 =
X
wi
i=1
n
T
= 11×n W Y =
X
X WY wi Yi
i=1

Hence the weighted least squares estimator

Pn
 −1 w i Yi
µb = T
X WX T
X WY = Pn
i=1

i=1 wi

has minimum variance as stated in the Gauss-Markov theorem.

1.2.2 Example 2: Galton’s peas

Consider Galton’s sweet pea data (galtonpeas in alr4 package). The recorded variables are:
• Parent: mean diameter of the parent (in 0.01 inches),
• Progeny: mean diameter of offspring (in 0.01 inches),
• SD: offspring diameter standard deviation (in 0.01 inches).
1.2. WEIGHTED LEAST SQUARES 11

The diameter standard deviation tends to be larger for offspring of parents with large diameter.
We compare an ordinary least squares regression of Progeny on Parent against a weighted least
squares regression using SD2 as weights.
lm1 <- lm(Progeny ~ Parent, data = galtonpeas)
lm2 <- lm(Progeny ~ Parent, weights=1/SDˆ2, data=galtonpeas)

OLS
WLS
Progeny mean diameter

17.0
16.5
16.0

15 16 17 18 19 20 21

Parent mean diameter

As we move from the ordinary least squares regression line to the weighted least squares
regression line we note a small shift due to observations with a small parent diameter being
given a larger weight.

Exercise 1
Consider an experiment in which it is known that the variance of the errors of the first two
observations is four times as large as the variance of the errors for the next three observations.
The errors are assumed to be independent. You are given the following model:

µ + ϵi i = 1 or i = 4,



Yi = µ + λ + ϵi i = 2 or i = 5,

µ + 2λ + ϵ3 i = 3.



Find the design matrix for this model and then calculate the generalised least squares estimator,
given that y T = (125, 62.5, 12.5, 100, 50).
Solution:
4 0 0 0 0 0 0 0 0
  1 
4
0 4 0 0 0 0 1
0 0 0
  
4
Var(ϵ) = σ 0
2
0 1 0 0 and so W = 0 0 1 0 0
   
  .
0 0 0 1 0 0 0 0 1 0
 
 
0 0 0 0 1 0 0 0 0 1
1 0
 

1 1

14 13 209.375
" # " #
Furthermore, X = 2 1
 and so X W X = 4 13 21 and X W y = 90.625 .
 T 1 T

0 1

1 1
12 CHAPTER 1. WEIGHTED LINEAR REGRESSION

21 −13
" #
With (X W X)
T −1
= 4
it follows that
125 −13 14

b = (X T W X)−1 X T W y = (103, −46.5)T .

β

1.3 The hat matrix

The fitted values of a weighted linear model are given by
 −1
Y
c = Xβ
b = X XT W X XT W Y = HY

where H is the hat matrix

 −1
H = X XT W X XT W

which projects the observations Y onto the fitted values Y

c

Y
c = HY .

The hat matrix for the weighted linear model is

• idempotent, that is HH = H, but
• not symmetric H T ̸= H.

Some sources use a different definition of the hat matrix

1
 −1 1
H∗ = W 2 X XT W X XT W 2

which is both symmetric and idempotent. This is the hat matrix on the scale of (X ∗ , Y ∗ ), that
is the matrix such that
c∗
Y = H ∗Y ∗.

This disagreement is not important.

• In practice we only ever calculate the diagonal elements of the hat matrix.
• The diagonal elements are the same under both definitions.

The hat-value or leverage hi is the ith diagonal element of the hat matrix

 −1
hi = wi xTi X T W X xi .
As we have seen in ST231, the sum of the hat values is equal to the number of parameters, that
is n
=
X
hi p.
i=1
1.4. WEIGHTED RESIDUALS 13

1.4 Weighted residuals

We have two different representations of the model which give us two different sets of residuals.
On the original scale (X, Y ) we have the response residuals defined as the difference between
the observed and fitted values:
ri = Yi − xTi β.b

For the transformed scale (X ∗ , Y ∗ ) we have the weighted residuals:

ri∗ = Yi∗ − (x∗i )T β

b
√

= wi Yi − xTi β)b
√
= w i ri

Both residuals have zero expectation

E(Ri | X) = 0,
E(Ri∗ | X) = 0,

but different variances

(1 − hi )ϕ
Var(Ri | X) = ,
wi

Var(Ri∗ | X) = (1 − hi )ϕ.

The weighted residuals are an example of what we will later call Pearson residuals when we
look at GLMs.
Unlike for weighted residuals, the variance of the response residual depends on the weight.
Therefore, we do not know if a response residual is large or small unless we also know its weight.
In contrast for weighted residuals, if there are no observations with high leverage, then we have
approximately

Var(Ri∗ | X) ≈ ϕ.

This makes the weighted residual more useful for model diagnostics.
The weighted residuals retain all the useful properties we associate with residuals in linear
models.
√
• If the model has an intercept term, then ni=1 wi ri∗ = ni=1 wi ri = 0.
P P

• An unbiased estimate of the dispersion parameter ϕ is the sum of squares of the weighted
residuals, divided by the residual degrees of freedom n − p

1 X n
ϕb = (r∗ )2 .
n − p i=1 i
14 CHAPTER 1. WEIGHTED LINEAR REGRESSION

The standardized residual is the residual divided by an estimate of its standard deviation.
Both residuals lead to the same definition:
√
(s) wi (Yi − xTi β)
ri = q
(1 − hi )ϕb

The standardized residual has asymptotic variance 1 as n → ∞.

We use standardized residuals to compare observations controlling for
• different weights w1 , . . . , wn , and
• different leverages h1 , . . . , hn .

Exercise 2
Show that the hat matrix on the scale of (X ∗ , Y ∗ ), that is the matrix such that
∗
Y
c = H ∗Y ∗

is given by −1
1 1

H∗ = W 2 X XT W X XT W 2 .
Solution
1 √ √ 1 1 1
Let W 2 = diag( w1 , . . . , wn ) and note that (W 2 )T = W 2 as W 2 is a diagonal matrix. Also
1
recall that we defined X ∗ = W 2 X where X is the design matrix on the original scale. Then,
using the expression for the hat matrix of an ordinary linear model, we obtain
 −1 
∗ T
 T
H ∗
= X ∗
X X ∗
X∗
 T −1 
1 1 1 1
= W 2X W 2X W 2X W 2 X)T
 −1
1 1
= W X X WX
2
T
XT W 2 .
Chapter 2 Binomial models

2.1 Bernoulli distribution

The simplest probability distribution is named after the Swiss mathematician Jacob Bernoulli
(1655–1705). Suppose Y ∈ {0, 1} and let µ = E(Y ), then

P(Y = 1) = µ, P(Y = 0) = 1−µ

and Y has the probability mass function

p(y) = µy (1 − µ)1−y .

The variance of the Bernoulli distribution depends on the mean µ as Var(Y ) = µ(1 − µ). We
write Y ∼ Bernoulli(µ) and Var(Y ) = V (µ). The variance function is depicted in the figure
below.
0.25
0.20
0.15
V(µ)

0.10
0.05
0.00

0.0 0.2 0.4 0.6 0.8 1.0

Note: Heteroscedasticity is built into the Bernoulli distribution as the variance is a quadratic
function in µ.
In linear models, we may try to resolve heteroscedasticity by transforming the outcome variable.
This will not work with the Bernoulli distribution.
• Any transformation of Y to a new variable (say Y ∗ ) will map all Y = 0 values to the same
value (say Y ∗ = a) and all Y = 1 values to the same value (say Y ∗ = b).
• The mean of the transformed variable is

µ∗ = E(Y ∗ ) = µ b + (1 − µ)a.

• The variance of the transformed variable is

Var(Y ∗ ) = (µ∗ − a) (b − µ∗ ).

Whatever we do, the variance still depends on the mean.

15
16 CHAPTER 2. BINOMIAL MODELS

2.2 Latent linear model

2.2.1 Definition
We can describe a binary response using a normal linear model by introducing latent variables.
Latent variables are random variables that are not observed. Let Zi ∼ N (xTi β, 1) where
Z1 , . . . , Zn are independent. Then, set

1 if Zi ≥ 0,
(
Yi =
0 if Zi < 0.

2.2.2 Illustration

The figure below shows an iid random sample of Z ∼ N (xT β, 1) of size n = 30.

(Here we have chosen xT β = −1.)

Index (i)

−2 xTβ 0

Now recall:

• for values of zi ≥ 0, the corresponding yi = 1;

• for values of zi < 0, the corresponding yi = 0.

To illustrate this derivation we have taken the illustration of our sample and coloured observations
where yi = 1 as black circles, while the observations where yi = 0 are still depicted as white
circles.
2.2. LATENT LINEAR MODEL 17

Index (i)

y=1
y=0

−2 xTβ 0

However, in practice we do not observe z1 , . . . , zn , only y1 , . . . , yn , represented here by the colour

of the points as illustrated in the figure below.
Index (i)

y=1
y=0

−1.0 −0.5 0.0 0.5 1.0

Hence we need to use a probability model for Yi based on the latent value of Zi .

Consider the density of Zi ∼ N (xTi β, 1) shown in the figure below. The shaded tail area shows

µi = P(Yi = 1) = P(Zi ≥ 0).

18 CHAPTER 2. BINOMIAL MODELS

0.35
density of Z

0.25
0.15
0.05 xTβ 0

We have

µi = P(Yi = 1) = P(Zi ≥ 0)
= P(N (0, 1) ≥ −xTi β)
= P(N (0, 1) ≤ xTi β)
= Φ(xTi β).

where Φ (“Phi”) is the cumulative distribution function of the standard normal distribution.
To avoid writing xTi β we introduce a new notation

ηi = xTi β

and new terminology. We call this quantity, represented by the Greek letter “eta”, the linear
predictor.
2.3. BINOMIAL OUTCOMES 19

The expectation µi of the outcome variable Yi is not the same as the linear predictor ηi = xTi β.
The two quantities are related by a smooth increasing function, namely Φ, as
µi = Φ(ηi ).

The mean function maps the linear predictor η onto the mean µ. The link function is the
inverse of the mean function. It maps the mean onto the linear predictor.
For the model considered here the link function is the inverse of Φ, that is
ηi = Φ−1 (µi ).
This link function is called the probit function. The name comes from an abbreviation of
“probability unit”.
The figure below on the left shows the mean function which, as discussed above, is given by the
cdf of a standard normal distribution. The figure below on the right shows the corresponding
link function, the so-called probit link function.

Mean function of the probit model Link function of the probit model
1.0

3
2
0.8

1
0.6

η
µ

0
0.4

−1
0.2

−2
0.0

−3

−3 −2 −1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0

η µ

Exercise 3
Work through Sections 1 - 3 (inclusive) from Computer Practical 1.
Solution: See the solutions for Computer Practical 1 available on moodle.

2.3 Binomial outcomes

The name “binomial” comes from the Latin for “two names”. A binomial outcome is used
whenever the response is a choice between two categories, for example
• Yes versus No,
• Success versus Failure,
• Claim versus No Claim.
For a group of observations with the same predictor variables, we can add together the responses
into a single outcome value for the whole group. As an example, this was done for the turbines
dataset.
20 CHAPTER 2. BINOMIAL MODELS

The binomial distribution is defined via a sequence of m independent Bernoulli random

variables with a common mean.
iid
Z1 , . . . , Zm ∼ Bernoulli(µ).

The likelihood function of the success probability µ given the observations z1 , . . . , zm is

m
L(µ | z1 , . . . , zm ) =
Y
p(zi | µ)
i=1
m
= µzi (1 − µ)1−zi
Y

i=1
Pm Pm
= µ z
i=1 i (1 − µ)m− z
i=1 i

= µy (1 − µ)m−y
where y =
Pm
i=1 zi .

Note: As the likelihood function depends on the data z1 , . . . , zm only through y, we have that
y is sufficient for µ.
The sum of m Bernoulli trials with common mean µ has a binomial distribution:
Y ∼ Binomial(m, µ)
where Y has probability mass function:
!
m
p(y) = µy (1 − µ)m−y for y ∈ {0, 1, . . . , m}.
y

Then
E(Y ) = mµ, and
Var(Y ) = mµ(1 − µ).

2.4 Scaled binomial distribution

In GLMs we use a different version of the binomial distribution. Let
1 X m
Y∗ = Zi .
m i=1
Then Y ∗ has a scaled binomial distribution with probability mass function
1 2
!
m ∗ ∗)
p(y )
∗
= µmy (1 − µ)m(1−y for y ∗ ∈ {0, , , . . . , 1}.
my ∗ m m

Then,
E(Y ∗ ) = µ, and
1 1
Var(Y ∗ ) = µ(1 − µ) = V (µ).
m m
for variance function V (µ) = µ(1 − µ).
2.5. IMPLEMENTATION IN R 21

The corresponding log-likelihood function is

 
l(µ) = log p(y | µ)
 
= m y log(µ) + (1 − y ) log(1 − µ) + . . .
∗ ∗

where we discard terms depending only on data (y ∗ , m).

This formulation allows us to use
• the proportion Y ∗ = Y /m as the response variable;
• the sample size m as a weight.

2.5 Implementation in R
Let’s consider the turbines data from the GLMsData package.
We can fit the probit binomial GLM in two ways.
The first approach uses the proportion y/m as response and weights m:
glm.out1 <- glm(Fissures/Turbines ~ Hours, weights=Turbines,
family=binomial(link="probit"), data=turbines)

The second approach uses the vector (y, m − y) as the response without any weights:
turbines <- transform(turbines, NoFissures = Turbines - Fissures)
glm.out2 <- glm(cbind(Fissures, NoFissures) ~ Hours,
family=binomial(link="probit"), data=turbines)

We use the sumary() function from the faraway package to produce an abbreviated summary
of the two models.
library(faraway)
sumary(glm.out1)

Estimate Std. Error z value Pr(>|z|)

(Intercept) -2.2758e+00 1.9742e-01 -11.5278 < 2.2e-16
Hours 5.7832e-04 6.2597e-05 9.2388 < 2.2e-16

n = 11 p = 2
Deviance = 9.81484 Null Deviance = 112.67005 (Difference = 102.85521)
sumary(glm.out2)

Estimate Std. Error z value Pr(>|z|)

(Intercept) -2.2758e+00 1.9742e-01 -11.5278 < 2.2e-16
Hours 5.7832e-04 6.2597e-05 9.2388 < 2.2e-16

n = 11 p = 2
Deviance = 9.81484 Null Deviance = 112.67005 (Difference = 102.85521)
22 CHAPTER 2. BINOMIAL MODELS

The two commands result in the same fitted model, and thus the summary output reports the
same estimates, standard errors, deviance, etc.
Let’s explore the fitted model graphically. The fitted line is non-linear so we need to use
the predict() function with type = "response" to indicate that we would like to get the
predictions in the scale of the response.
plot(Fissures/Turbines ~ Hours, data=turbines, pch=16)
hvec <- seq(from=min(turbines$Hours), to=max(turbines$Hours), length=101)
ypred <- predict(glm.out1, newdata=list(Hours=hvec), type="response")
lines(hvec, ypred, lwd=2)

If we examine the fitted curve for very small and very large number of hours we observe
horizontal asymptotes at 0 and at 1 (beware though that in the example here this is a substantial
extrapolation!)

1.0
0.0 0.1 0.2 0.3 0.4 0.5 0.6

0.8
Fissures/Turbines

Fissures/Turbines

0.6
0.4
0.2
0.0

1000 2000 3000 4000 0 2000 4000 6000 8000 10000

Hours Hours

2.6 Logistic regression

To interpret the parameters of a binary regression model we need to be able to translate between
the linear predictor and the mean. We do this via the link function ηi = g(µi ) and mean function
µi = h(ηi ). But the probit link has no closed form.
The logit link function is an alternative that has a simple enough closed form.

2.6.1 The logit link function

The logistic function is defined as
exp(η) 1
h(η) = = .
1 + exp(η) 1 + exp(−η)
It is monotonically increasing with horizontal asymptotes at 0 and at 1 and thus is suited to
modelling probabilities.
The corresponding link function is the logit link function:
µ
 
g(µ) = log = η.
1−µ
2.6. LOGISTIC REGRESSION 23

The figure below shows the logistic function on the left and the logit link function on the right.
1.0
µ
logit(µ) = log( )

4
1−µ
0.8

2
0.6

η
µ

0
0.4

exp(η)
logistic(η) =

−2
1 + exp(η)
0.2

−4
0.0

−4 −2 0 2 4 0.0 0.2 0.4 0.6 0.8 1.0

η µ

A GLM for binary response data using a logistic mean function and logit link is known as
a logistic regression model. Before we discuss how to interpret the parameters of a logistic
regression model in terms of log odds ratios, we first introduce the concept of odds.

2.6.2 Odds of an event occurring

Odds are another way to express probabilities and are commonly used in gambling.
Suppose an event E has known probability µ. You think it will happen; I do not. I bet 1 pound
that event E will not occur.
• If E occurs, then you get to keep my stake of 1 pound.
• If E does not occur, then you return my stake and, in addition, give me L pounds.
Suppose we both agree on the probability µ but disagree only in our predictions of what will
happen. If the game is fair, then your expected gain and mine are both equal to 0.

µ × 1 − (1 − µ) × L = 0.

Rearranging gives
µ
L = .
1−µ
These are the odds of E occurring.
Next let us consider how we make use of odds when interpreting the parameters of a logistic
regression model.
24 CHAPTER 2. BINOMIAL MODELS

2.6.3 Interpreting parameters for binary predictors

In the following we will use terminology from epidemiology:
• If Yi = 1, then the individual is diseased.
• If Yi = 0, then the individual is healthy.
• If Xi = 1, then the individual is exposed.
• If Xi = 0, then the individual is unexposed.
Let
µ(0) = P(Yi = 1 | Xi = 0),
µ(1) = P(Yi = 1 | Xi = 1).

µ(0)
The odds of disease (Yi = 1) in the unexposed (Xi = 0) is 1−µ(0)
.
µ(1)
The odds of disease (Yi = 1) in the exposed (Xi = 1) is 1−µ(1)
.
The odds ratio is given by
µ(1) 1 − µ(0)
Ω = ×
1 − µ(1) µ(0)

Odds of disease for the exposed

= .
Odds of disease for the unexposed

If µi = E(Yi ), then using the logit link we have

µi
 
g(µi ) = log = α + βxi .
1 − µi

It follows that
µi
 
= exp α + βxi .
1 − µi

The odds of disease (Yi = 1) in the unexposed (Xi = 0) is

µ(0)
 
= exp α + β × 0 = exp(α).
1 − µ(0)

So the intercept α is the log odds of disease when Xi = 0.

The odds of disease (Yi = 1) in the exposed (Xi = 1) is
µ(1)
 
= exp α + β × 1 = exp(α) × exp(β).
1 − µ(1)

So the odds ratio is

exp(α) × exp(β)
Ω = = exp(β).
exp(α)

Hence β = log(Ω) is the log odds ratio for Xi = 1 relative to Xi = 0.

2.6. LOGISTIC REGRESSION 25

2.6.4 Interpreting parameters for continuous predictors

Suppose Xi is continuous instead of binary. We still assume a logistic regression model with
µi
 
log = α + βxi .
1 − µi
Then β is the log odds ratio for a unit increase in Xi .
• If we compare the odds for Xi = x + 1 relative to the odds for Xi = x, then the odds ratio
is exp(β).
• The magnitude of β depends on the units (e.g. recall the example of cigarettes versus
packets discussed in the practical class).

2.6.5 Interpreting parameters for categorical predictors

If we have a categorical predictor variable (factor) with K levels {1, . . . , K} and we use treatment
coding, then the logistic regression model is

K
µi
 
log = α+
X
βk 1[Xi =k]
1 − µi k=2

where the original categorical variable has been replaced by K − 1 indicator variables.
• The reference category (k=1) is “absorbed” into the intercept α.
• For k ∈ {2, . . . , K}, the parameter βk is the log odds ratio for level k relative to level 1.

2.6.6 Interpreting the intercept

For all cases, the intercept parameter α is the log odds of the outcome Yi = 1 when all predictors
are at their reference value, that is
• zero for binary variables,
• zero for continuous variables,
• reference level for factors when using using treatment coding.

2.6.7 Interpreting multiple logistic regression

We may have multiple predictors in a logistic regression model, for example
µi
 
log = α + βxi + γzi .
1 − µi
The interpretation of β is the same as before but assumes holding all other predictor variables
fixed (that is zi in the above example). The terms “adjusting for Z” or “controlling for Z”
indicate that we consider the effect of X while keeping Z fixed.
26 CHAPTER 2. BINOMIAL MODELS

2.6.8 Epidemiological example

Suppose we have a single binary predictor Xi ∈ {0, 1} and a binary response Yi ∈ {0, 1} for N
independent observations.
We can summarise the data in a 2 × 2 table.

Xi = 1 Xi = 0 Total
Yi = 1 D(1) D(0) D
Yi = 0 H (1) H (0) H
Total N (1) N (0) N

Then

D(0) ∼ Binomial(N (0) , µ(0) ),

D(1) ∼ Binomial(N (1) , µ(1) ),

where

µ(0) = P(Yi = 1 | Xi = 0),

µ(1) = P(Yi = 1 | Xi = 1).

Recall, the odds ratio is given by

µ(1) 1 − µ(0)
Ω = ×
1 − µ(1) µ(0)

Odds of disease for the exposed

= .
Odds of disease for the unexposed

Ω = 1 if and only if µ(1) = µ(0) , so the risk is the same in both exposed and unexposed.
The maximum likelihood estimate for Ω is the cross-product ratio

D(1) H (0)
Ω
b = .
H (1) D(0)

We note the balance between the numerator and denominator

P(Yi = 1 | Xi = 1) P(Yi = 0 | Xi = 0)
Ω = .
P(Yi = 1 | Xi = 0) P(Yi = 0 | Xi = 1)

We can use Bayes’ theorem to invert the conditional probabilities. For a, b ∈ {0, 1}

P(Yi = a)
P(Yi = a | Xi = b) = P(Xi = b | Yi = a) × .
P(Xi = b)
2.6. LOGISTIC REGRESSION 27

Due to the balance between numerator and denominator in the expression for Ω, the second
factor cancels out and so
P(Xi = 1 | Yi = 1) P(Xi = 0 | Yi = 0)
Ω = .
P(Xi = 1 | Yi = 0) P(Xi = 0 | Yi = 1)

So we can express Ω as the odds ratio of Xi = 1 for Yi = 1 relative to Yi = 0.

Inversion of the odds ratio is useful in epidemiology.
We can estimate the odds ratio by recruiting participants using their disease status.
• Identify D individuals with a given disease (Yi = 1).
• Find a group of H healthy individuals (Yi = 0).
• Measure the exposure X and classify participants as exposed (Xi = 1) or unexposed
(Xi = 0).
This is called a case control study, which is very useful to study rare diseases.

2.6.9 The rare outcome interpretation

Suppose the outcome Yi = 1 is rare that is, P(Yi = 1) is very small for both groups (exposed
and unexposed):

µ(0) ≪ 1,
µ(1) ≪ 1.

Then
µ(1) (1 − µ(0) ) µ(1)
Ω = ≈
µ(0) (1 − µ(1) ) µ(0)

and so the odds ratio is approximately equal to the relative risk, or risk ratio.

Exercise 4
Work through Sections 4 - 5 from Computer Practical 1.
Solution: See the solutions for Computer Practical 1 on moodle.
28 CHAPTER 2. BINOMIAL MODELS

2.7 Tolerance distributions

2.7.1 Generalizing the latent linear model
The earlier construction of the latent linear model for a probit regression can be generalized
beyond the normal distribution. We can write the latent variable Zi as

Zi = ηi + ϵi

where ϵi has a continuous, real valued distribution. The derived random variable Yi remains a
deterministic function of Zi given by

1 if Zi ≥ 0,
(
Yi =
0 if Zi < 0.

We have

µi = P(Yi = 1) = P(Zi ≥ 0)
= P(ηi + ϵi ≥ 0)
= P(ϵi ≥ −ηi )
= 1 − F (−ηi ).

where F is the cumulative distribution function of the error term ϵi .

The distribution of the latent variable Z (or equivalently the error term ϵ) is called the tolerance
distribution. If the tolerance distribution is symmetric, then we can simplify

µ = 1 − F (−η) = F (η)

and so the mean/response function is the cumulative distribution function F .

2.7. TOLERANCE DISTRIBUTIONS 29

2.7.2 Example tolerance distributions

[Link] The logistic distribution
The logistic distribution has a symmetric bell-shaped density, like the normal distribution, but
with longer tails.

Comparison of the standard normal and the logistic density

normal
logistic

−4 −2 0 2 4

The standard logistic probability density function is given by

exp(x)
f (x) =
(1 + exp(x))2

which gives the cumulative distribution function

exp(x)
F (x) = .
1 + exp(x)

Unlike the normal cdf, the logistic cdf can be written in closed form. This allows us to interpret
the parameters more easily.
A binary regression model using the logistic distribution as the tolerance distribution is called a
logistic regression model. It has mean function h : R → (0, 1) with

exp(η)
µ = h(η) = .
1 + exp(η)

The corresponding link function g : (0, 1) → R is given by

µ
 
η = g(µ) = h−1 (µ) = log .
1−µ
This link function is called the logit link.
30 CHAPTER 2. BINOMIAL MODELS

[Link] The Gumbel distribution

The Gumbel distribution has a skewed density and is commonly used in applications of extreme
value theory.

Standard Gumbel density

0.3
0.2
0.1
0.0

−4 −2 0 2 4

The probability density function of the standard Gumbel distribution is given by

f (x) = exp(−x − exp(−x))

which leads to the cumulative distribution function

 
F (x) = exp − exp(−x) .

Recall, the Gumbel distribution is asymmetric, so h(η) ̸= F (η).

The formula for the mean function is
 
µ = h(η) = 1 − F (−η) = 1 − exp − exp(η) .

The link function is given by

η = g(µ) = h−1 (µ) = log(− log(1 − µ)).

This link is called the complementary log-log link function.

2.7. TOLERANCE DISTRIBUTIONS 31

[Link] The Cauchy distribution

The Cauchy distribution is the special name given to a t-distribution on 1 degree of freedom.

Comparison of the Cauchy and the standard normal density.

normal
Cauchy

−4 −2 0 2 4

The probability density function of the Cauchy distribution

1
f (x) =
π (1 + x2 )

which corresponds to the cumulative distribution function

1 1
F (x) = arctan(x) + .
π 2
The mean and variance of the Cauchy distribution are not defined.
As the Cauchy distribution is symmetric, the mean function is equal to the cdf.
1 1
h(η) = arctan(η) + .
π 2
The link function is
g(µ) = tan(π (µ − 0.5)).
and is called the cauchit link.

[Link] Summary of tolerance distributions

Tolerance distribution Link function R link name

Normal Probit “probit”
Logistic Logit “logit”
Gumbel Complementary log log “cloglog”
Cauchy Cauchit “cauchit”
32 CHAPTER 2. BINOMIAL MODELS

[Link] R illustration of the different link functions

In the following we define the various link functions for a binary GLM and fit a model for each
link function to the turbines data.
library(GLMsData)
library(kableExtra)
data(turbines)

links <- c(Probit="probit", Logit="logit", CLogLog="cloglog")

[Link] <- function(L) {

glm(Fissures/Turbines ~ Hours, weights=Turbines,
family=binomial(link=L), data=turbines)}

[Link] <- lapply(links, [Link])

kable(round(sapply([Link], coef), 6)) %>%
kable_styling(position = "center", latex_options = "hold_position")

Probit Logit CLogLog

(Intercept) -2.275807 -3.923597 -3.60328
Hours 0.000578 0.000999 0.00081

We note differences between the estimated coefficients, for example between the model using the
probit link and the one using the logit link. However, while the coefficients differ substantially,
the models with different link functions may still produce very similar predictions as can be
seen in the figure below. It shows a scatterplot of the data together with the predictions from
the various models.

probit
0.0 0.1 0.2 0.3 0.4 0.5 0.6

logit
cloglog
Fissures/Turbines

1000 2000 3000 4000

Hours

Exercise 5
1. For each of the tolerance distributions presented verify
a) F ′ (x) = f (x) (the derivative of its cdf is its pdf);
2.7. TOLERANCE DISTRIBUTIONS 33

b) g(h(η)) = η and h(g(µ)) = µ (the link function is the inverse of the mean function).
2. Repeat the experiment from Section [Link] but this time include the cauchit link function.

Solutions for Part 1:

Logistic distribution:
a) We have F (x) = exp(x)
1+exp(x)
and f (x) = exp(x)
(1+exp(x))2
.

Note that if a(x) = x

x+1
, then a′ (x) = 1
(1+x)2
.
 
As F (x) = a exp(x) we have

  exp(x)
F ′ (x) = a′ exp(x) exp(x) =  2 = f (x).
1 + exp(x)

 
b) Recall that h(η) = exp(η)
1+exp(η)
and g(µ) = log µ
1−µ
.
Then, exp(g(µ)) = µ
1−µ
and so
 
  exp g(µ) µ
1−µ
µ
1−µ
h g(µ) = = = = µ.
1 + 1−µ
µ 1
 
1 + exp g(µ) 1−µ

Also  
h(η) exp(η)/ 1 + exp(η)
= = exp(η)
1 − h(η)
 
1/ 1 + exp(η)
and so
h(η)
     
g h(η) = log = log exp(η) = η.
1 − h(η)
34 CHAPTER 2. BINOMIAL MODELS

Standard Gumbel distribution:

   
a) We have F (x) = exp − exp(−x) and f (x) = exp − x − exp(−x) .

Note that if a(x) = − exp(−x), then a′ (x) = exp(−x).

As F (x) = exp(a(x)) we have
     
F (x) = exp a(x) a (x) = exp −exp(−x) exp(−x) = exp −x−exp(−x)
′ ′
= f (x).

   
b) Recall that h(η) = 1 − exp − exp(η) and g(µ) = log − log(1 − µ) .
 
Now note that − exp g(µ) = log(1 − µ) and so
      
h g(µ) = 1 − exp − exp g(µ) = 1 − exp log(1 − µ) = µ.
 
Also log 1 − h(η) = − exp(η) and so
      
g h(η) = log − log 1 − h(η) = log exp(η) = η.

Cauchy distribution:
a) For the Cauchy distribution F (x) = 1
π
arctan(x) + 1
2
and f (x) = 1
π (1+x2 )
.

The derivative of arctan(x) is 1/(1 + x2 ) and so F ′ (x) = 1 1

π (1+x2 )
= f (x).
 
b) Recall that h(η) = 1
π
arctan(η) + 1
2
and g(µ) = tan π (µ − 0.5) .
 
We have arctan g(µ) = π(µ − 1/2) and so
  1   1
h g(µ) = arctan g(µ) + = µ.
π 2
 
Also π h(η) − 1/2 = arctan(η) and so
      
g h(η) = tan π h(η) − 0.5 = tan arctan(η) = η.
2.7. TOLERANCE DISTRIBUTIONS 35

Solutions for Part 2:

Adapting the code from the lectures:
library(GLMsData); library(kableExtra); data(turbines)

## Fitting a binomial GLM for each link function

links <- c(Probit="probit", Logit="logit", CLogLog="cloglog", Cauchit="cauchit")
[Link] <- function(L) {
glm(Fissures/Turbines ~ Hours, weights=Turbines,
family=binomial(link=L), data=turbines)
}

## Print out of the estimated coefficients

[Link] <- lapply(links, [Link])
kable(round(sapply([Link], coef), 6)) %>%
kable_styling(position = "center", latex_options = "hold_position")

Probit Logit CLogLog Cauchit

(Intercept) -2.275807 -3.923597 -3.60328 -4.467548
Hours 0.000578 0.000999 0.00081 0.001139

## Plotting the predicted probability of fissures as a function of hours

newHours <- seq(0, 5000, length=100)
[Link] <- function(model) {
predict(model, type="response", newdata=[Link](Hours = newHours))}
predicted <- sapply([Link], [Link])
plot(Fissures/Turbines ~ Hours, data=turbines, pch=16)
library(RColorBrewer)
palette <- [Link](4, "Dark2")
matlines(newHours, predicted, lty=1:4, lwd=2, col=palette)
legend("topleft", lty=1:4, lwd=2, col=palette, legend=links)

probit
0.0 0.1 0.2 0.3 0.4 0.5 0.6

logit
cloglog
Fissures/Turbines

cauchit

1000 2000 3000 4000

Hours

As previously seen, the predictions of the fitted models for the probit, logit and cloglog link
are all very similar. The fitted Binomial GLM using the cauchit link, however, gives different
36 CHAPTER 2. BINOMIAL MODELS

predictions. We note that it consistently over-predicts for runtimes of less than 2000 hours. For
example, the group of turbines that were run for 400 hours had no fissures, but the model gives
a predicted probability of fissures of roughly 8%.
Chapter 3 Poisson regression models

3.1 Moments and cumulants

3.1.1 Moments
The moments of a random variable Y are the expected values of positive integer powers:

m1 = E(Y ),
m2 = E(Y 2 ),
m3 = E(Y 3 ),

and so on.
The moments of a distribution do not always exist. For example, for the t-distribution on d
degrees of freedom, mk does not exist for k ≥ d.
When all moments exist, they uniquely characterize the distribution, if Carleman’s condition is
satisfied: ∞
1
(m2k )− 2k = +∞.
X

k=1

Informally this condition ensures that the absolute value of moment mk does not grow too
quickly with k. For example, the lognormal distribution has moments

mk = exp(k 2 /2)

which grow rapidly with k. Thus the lognormal distribution does not satisfy Carleman’s
condition and is not uniquely determined by its moments.
If the moment generating function M (t) can be defined for t in some neighbourhood of zero,
then it uniquely characterizes the distribution. Recall the definition of the moment generating
function:

∞
tr Y r ∞
tr E (Y r ) ∞
tr mr
!
M (t) = E (exp(tY )) = = =
X X X
E
r=0 r! r=0 r! r=0 r!

By definition  
M (0) = E exp(0) = 1.

The moments can be recovered from the derivatives of M (t) at zero:

dr M (0)
mr = .
dtr
37
38 CHAPTER 3. POISSON REGRESSION MODELS

3.1.2 Cumulant generating function

The cumulant generating function is the log of the moment generating function:
h i h  i
K(t) = log M (t) = log E exp(tY )

The rth cumulant κr is defined as the rth derivative of the cumulant generating function
evaluated at t = 0:
dr K(0)
κr = .
dtr
The cumulants can be expressed in terms of moments up to the same order (you don’t need to
memorise these!):
κ1 = m1
κ2 = m2 − m21
κ3 = m3 − 2m1 m22 + 2m31
...

Note:
• The first 2 cumulants are well-known quantities: the mean and the variance.
• The 3rd cumulant is a measure of asymmetry, related to skewness:
κ3
Skewness = 3/2
.
κ2
• Higher order cumulants are harder to interpret.
Unlike moments, cumulants are additive for independent random variables.
If S = Yi for independent Y1 , . . . , Yn , then
Pn
i=1
n
KS (t) = KYi (t).
X

i=1

Proof:
 h i  h n i n
 hY i
KS (t) = log E exp(tS) = log E exp(t Yi ) = log E exp(tYi )
X

i=1 i=1
n
Y h i n  h i n
= log E exp(tYi ) = log E exp(tYi ) = KYi (t).
X X

i=1 i=1 i=1

It follows that the rth cumulant of S is the sum of the rth cumulants of Y1 , . . . , Yn :
dr KS (0) n
dr KYi (0)
=
X
.
dtr i=1 dtr

We will be dealing with a family of distributions called exponential dispersion models

(EDMs).
• All moments and cumulants are finite for EDMs. M (t) and K(t) are always defined.
• K(t) takes a particularly simple form for EDMs.
• We will use K(t) to derive the mean and variance of EDMs.
• In some proofs we will also use the fact that M (t) or K(t) uniquely defines a distribution.
3.1. MOMENTS AND CUMULANTS 39

3.1.3 Summary on moments and cumulants

• For the distributions considered in this module, the moment generation function M (t)
and the cumulant function K(t) are always defined and uniquely characterise the
distribution.
• The moment generating function M (t) of Y is given by
∞
tr mr
M (t) = E (exp(tY )) =
X
.
r=0 r!

• mr = E(Y r ) is the rth moment and satisfies

dr M (0)
mr = .
dtr

• The cumulant generating function K(t) of Y is given by

h i
K(t) = log M (t)

where the rth cumulant κr is defined as

dr K(0)
κr = .
dtr

• κ1 = E(Y ) and κ2 = Var(Y ).

• If S = Yi for independent Y1 , . . . , Yn , then
Pn
i=1

n
KS (t) = KYi (t).
X

i=1

• The rth cumulant of S is the sum of the rth cumulants of Y1 , . . . , Yn :

dr KS (0) n
dr KYi (0)
=
X
.
dtr i=1 dtr

Exercise 6
Show that

dr M (0)
mr = ,
dtr
κ1 = m1 ,
κ2 = m2 − m21 .

Solution to Exercise 6:
40 CHAPTER 3. POISSON REGRESSION MODELS

∞
tk mk
M (t) = E (exp(tY )) =
X
.
k=0 k!

Therefore
dr M (t) ∞
tk−r
=
X
mk
dtr k=r (k − r)!

and so
dr M (0)
= mr .
dtr
Furthermore as K(t) = log (M (t)) we have

dK(t) 1 dM (t)
=
dt M (t) dt

and
d2 K(t) 1 d2 M (t) 1 dM (t) 2
 
= − .
dt2 M (t) dt2 M (t)2 dt
Hence, setting t = 0 and noting that M (0) = 1 we have
1
κ1 = K ′ (0) = M ′ (0) = m1
M (0)

and
1 1
 2
κ2 = K ′′ (0) = M ′′ (0) − M ′ (0) = m2 − m21 .
M (0) M (0)2

Exercise 7
Work through Sections 1 - 3 of Computer Practical 2.
3.2. THE POISSON DISTRIBUTION 41

3.2 The Poisson distribution

The Poisson distribution is a discrete distribution with support on N0 :
µr
P(Y = r) = exp(−µ)
r!
for r = 0, 1, 2, . . . where
µ = E(Y ).

The Poisson distribution is commonly used to model count data.

Below are some example histograms for an iid sample of size n = 10, 000 from a Poisson(µ)
distribution. For small values of µ, the Poisson distribution is skewed with a high probability of
Y = 0 and a right tail that decays quickly.
For moderate values of µ the Poisson distribution is more bell-shaped, but still positively skewed.

µ=1 µ=5
1500
3000
Frequency

Frequency

1000
2000
1000

500
0

0 2 4 6 0 5 10 15

y y

For large values of µ the Poisson distribution converges to a normal distribution.

µ = 100 Normal Q−Q Plot

140
400

120
Sample Quantiles
300
Frequency

100
200
100

80
0

80 100 120 140 −4 −2 0 2 4

y Theoretical Quantiles
42 CHAPTER 3. POISSON REGRESSION MODELS

3.2.1 Characteristics of the Poisson distribution

The moment generating function of a Poisson(µ) distribution is
   h i
M (t) = E exp(tY ) = exp µ exp(t) − 1

with the corresponding cumulant generating function

   
K(t) = log M (t) = µ exp(t) − 1 .

Proof:
  ∞
M (t) = E exp(tY ) = exp(tr) P(Y = r)
X

r=0

∞
µr exp(−µ) ∞
exp(−µ)
= exp(tr) = exp(t)r µr
X X

r=0 r! r=0 r!
h ir
∞ exp(t)µ  
= exp(−µ) = exp(−µ) exp exp(t)µ
X

r=0 r!
   h i
= exp exp(t)µ − µ = exp µ exp(t) − 1 .

It follows that K(t) = µ [exp(t) − 1] and thus

K ′ (t) = µ exp(t),
K ′′ (t) = µ exp(t).

Therefore

E(Y ) = K ′ (0) = µ(exp(0)) = µ,

Var(Y ) = K ′′ (0) = µ(exp(0)) = µ.

Thus the variance of a Poisson distribution is equal to its mean.

3.2. THE POISSON DISTRIBUTION 43

3.2.2 Example - Boston bicyle hires

The data shows the count of bicycles hired in each hour.

Bicycle hire in Boston in 2011

500
number of hires

300
100
0

0 2 4 6 8 10 12 14 16 18 20 22

hour of the day (24 hour clock)

There is more variation at times of high demand than at times of low demand.

3.2.3 Sums of independent Poisson random variables.

Consider independent random variables Y1 , . . . , Yn where

Yi ∼ Poisson(µi ).

Let S = Yi , then
Pn
i=1
S ∼ Poisson(µS )
where n
=
X
µS µi .
i=1

Proof:
The cumulant generating function of S is given by
n
KS (t) = KYi (t)
X

i=1
n h i
= µi exp(t) − 1
X

i=1

h n
iX
= exp(t) − 1 µi
i=1

h i
= exp(t) − 1 µS

which is the cumulant generating function of a Poisson distribution with mean µS .

44 CHAPTER 3. POISSON REGRESSION MODELS

3.3 Poisson regression for counts

In a regression model, we explain the expectation µi in terms of a p-vector of predictor variables

xi and parameters β.

For Poisson outcome data, we normally use a log link to construct a regression model:

log(µi ) = xTi β.

The log link ensures that the expectation µi is always positive as

 
µi = exp xTi β > 0

for all β ∈ Rp .

Other link functions are possible, but they constrain the possible values of β.

3.4 Poisson regression for rates

3.4.1 Poisson process

A Poisson process is a continuous-time stochastic process, often used to represent queues.

Events occur in time with constant rate λ, that is

P(event ∈ (t, t + δ])

lim = λ,
δ→0 δ

independent of the history of the process and of t.

Event times T1 , T2 , . . . can be represented as points on the time line.

3.4. POISSON REGRESSION FOR RATES 45

T1 T2 T3 T4 T5
0
Time

The inter-event times have independent exponential distributions with common mean λ1 .

Suppose we observe the Poisson process up to a censoring time C.

46 CHAPTER 3. POISSON REGRESSION MODELS

T1 T2 T3 C T4 T5
0
Time

Events occurring after the censoring time are not counted. Let Y be the number of events
occurring before time C then

Y ∼ Poisson(µ) with µ = λC.

Observations may have different censoring times. With different observation windows, we have

Yi ∼ Poisson(µi ) with µi = λi C i .

The expected event count will vary between observations, even if they have the same rate λi = λ.
3.4. POISSON REGRESSION FOR RATES 47

C1
0 Y1 = 1

C2
0 Y2 = 3

How do we account for this in our Poisson regression model?

Using the log link we have

log(µi ) = log(λi Ci )
= log(λi ) + log(Ci )

If we assume
λi = exp(xTi β)
then
log(µi ) = xTi β + log(Ci ).
The term log(Ci ) is called an offset.
An offset is not the same as a predictor variable. If we put log(Ci ) in our model as a predictor,
then we get
log(µi ) = xTi β + γ log(Ci )
with a parameter γ to be estimated from the data. With an offset, we fix γ = 1.
The Poisson model with log link and an offset is used whenever we have a denominator Di such
that
E(Yi ) = Di exp(xTi β)

The denominator does not need to be time. It could be, for example
48 CHAPTER 3. POISSON REGRESSION MODELS

• spatial area,
• population size,
• volume.

3.4.2 Aggregation of Poisson data

Suppose we have Poisson observations Y1 , Y2 , . . . , Ym that have the same observed predictor
values – hence the same rate λ – but different observation times C1 , C2 , . . . , Cm .
Then the log-likelihood function is given by

m
= (Yi log(λ) − λCi ) + . . .
X
l(λ)
i=1
m
X  m
X 
= log(λ) Yi − λ Ci + . . .
i=1 i=1

where terms not depending on λ have been discarded. So the total event count Yi and the
Pm
i=1
sum of censoring times mi=1 Ci are sufficient for λ as
P

m m
 X 
Poisson λ
X
Yi ∼ Ci .
i=1 i=1

(See also the result on sums of independent Poisson random variables in Section 3.2.3!)
If all our predictor variables are categorical, we can aggregate the data. For each combination
of predictor variables, all we need is
• the total number of events Yi and
Pm
i=1

• the sum of observation times Ci .

Pm
i=1

We apply Poisson regression with a log link and use log( Ci ) as an offset. The next section
Pm
i=1
discusses an example.

3.4.3 Example - The ships data

The dataset ships in the MASS package gives the number of damage incidents and aggregate
months of service by ship type, year of construction, and period of operation.
library(MASS)
data(ships)
head(ships,7)
type year period service incidents
1 A 60 60 127 0
2 A 60 75 63 0
3 A 65 60 1095 3
4 A 65 75 1095 4
5 A 70 60 1512 6
6 A 70 75 3353 18
7 A 75 60 0 0
3.4. POISSON REGRESSION FOR RATES 49

This is aggregated data. As we have 5 ship types, 4 different year of construction labels and 2
periods of operation, the data was aggregated over 5 × 4 × 2 = 40 groups. The number of ships
for each group was not recorded. However, the total services times (sum of censoring times)
and the total number of damage incidents (sum of events) are provided for each group and are
sufficient to fit a Poisson model for rates.
The variable service gives the aggregate months of service for each group in the data. However
as aggregating the data leads to groups where the period of operation occurs before the year of
construction, the ships in those groups have a service time of zero months. These groups are
excluded from further analysis.
The remaining groups have different total service times and thus different exposure to the risk
of a damage incident and so we are interested in the rate of damage incidents per month of
service time.
ships <- subset(ships, service > 0)

rate <- sum(ships$incidents)/sum(ships$service)

round(rate, 5)
[1] 0.00218

We use log(service) as an offset and fit a Poisson GLM with log-link. Here we are fitting a
null model, that is a model with no predictor variables.
[Link] <- glm(incidents ~ offset(log(service)),
family=poisson(link="log"), data=ships)
round(coef([Link]), 2)
(Intercept)
-6.13
round(exp(coef([Link])), 5)
(Intercept)
0.00218

Suppose we would like to examine the incident risks for different types of ships, then we may fit
the following Poisson GLM:
[Link] <- glm(incidents ~ 0 + type + offset(log(service)),
family=poisson(link="log"), data=ships)
round(exp(coef([Link])), 4)
typeA typeB typeC typeD typeE
0.0044 0.0018 0.0019 0.0038 0.0062

Exercise 8
Work through Sections 3 and 4 of Computer Practical 2.
Chapter 4 Exponential dispersion models

4.1 Motivation
Before we define the class of exponential dispersion models, we take a moment to understand
why they are needed. Consider a simple problem of estimating a common mean.
• Let Y1 , . . . , Yn be i.i.d. random variables with mean µ.

• We propose a family of probability density functions p(y | µ) parameterized by µ.

We may estimate the mean µ using two alternative approaches.
1. Sample mean:
1X n
µ = yi .
n i=1
2. Maximum likelihood estimate:
n  
= arg max log p(yi | µ) .
X
µb
µ i=1

Question: Is the maximum likelihood estimate µb identical to the sample mean µ?

Answer: The maximum likelihood estimate µb is identical to the sample mean µ for some
probability models (e.g. normal, gamma). But this is not true for all probability models
(e.g. t-distribution).
• Exponential dispersion models (EDMs) are probability models for which µb = µ for i.i.d.
observations with common mean µ.

• This property uniquely characterizes EDMs.

• For non-EDMs, the sample mean may still be a consistent and efficient estimator of µ.

• Hence µ may be “close to” µ,

b but is not necessarily identical.

50
4.2. DEFINITION OF AN EDM 51

4.2 Definition of an EDM

An EDM is a distribution from the exponential family of distributions. The probability density
function (or probability mass function) of an EDM can be put in the canonical form:

θy − b(θ)
 
p(y | θ, ϕ) = a(y, ϕ) exp
ϕ

where
• θ ∈ Θ is the canonical parameter and Θ = {θ ∈ R : |b(θ)| < ∞},
• b(θ) is the cumulant function,
• ϕ ∈ R+ is the dispersion parameter. The dispersion parameter may be free, in which
case it is an additional parameter to be estimated, or it may be fixed to a known value
(usually ϕ = 1).
• a(y, ϕ) is the normalizing function. It ensures that
Z
p(y | θ, ϕ)dy = 1
y∈S

where S is the support (the permitted values of y). The normalizing function does not
depend on θ and plays no role in inference on θ.
Examples of EDMs include the following standard distributions:
• Normal (Gaussian),
• Binomial,
• Poisson,
• Gamma,
• Inverse Gaussian,
• Negative Binomial.
The support S of an EDM is determined by its normalizing function a(y, ϕ). Different EDMs
have different support, see the table below.

Distribution Support S
Normal R
Poisson N0
Scaled Binomial {0, m1 , m2 , . . . , 1}
Negative Binomial N0
Gamma R+
Inverse Gaussian R+
52 CHAPTER 4. EXPONENTIAL DISPERSION MODELS

4.3 Weighted EDMs

Suppose we have independent Y1 , . . . , Yn from an EDM with the same canonical parameter θ
but different dispersion parameters ϕ1 , . . . , ϕn . We can extend our definition of EDMs to include
this case if we assume
ϕ
ϕi =
wi
for known weights w1 , . . . , wn and common dispersion parameter ϕ.
The density function is then
h i
w
i θyi − b(θ) 
p(yi | θ, ϕ) = a(yi , ϕ/wi ) exp .
ϕ

For fixed ϕ, the log likelihood of θ is

1 Xn
  h i
log L(θ | ϕ, y) = wi θyi − b(θ) + . . .
ϕ i=1

where terms depending on the normalizing function a(y, ϕ/wi ) have been omitted.
An observation with weight wi ∈ N makes the same contribution to the log likelihood as wi
identical observations with weight 1.

4.4 Examples
4.4.1 Normal distribution
Recall the canonical form of the pdf/pmf of an EDM:

θy − b(θ)
 
p(y | θ, ϕ) = a(y, ϕ) exp .
ϕ

The density of a N (µ, σ 2 ) can be written in canonical form as

1 (y − µ)2
 
p(y | µ, σ 2 ) = √ exp −
2πσ 2 2σ 2

1 −y 2 + 2µy − µ2
 
= √ exp
2πσ 2 2σ 2

1 −y 2 µy − 12 µ2
   
= √ exp exp .
2πσ 2 2σ 2 σ2

Comparing this to the general canonical form

4.4. EXAMPLES 53

θy − b(θ)
 
p(y | θ, ϕ) = a(y, ϕ) exp
ϕ
we deduce that
1 −y 2 θy − 12 θ2
   
p(y | µ, σ 2 ) = p(y | θ, ϕ) = √ exp exp
2πϕ 2ϕ ϕ
is a pdf in canonical form with

θ = µ
ϕ = σ2,
θ2
b(θ) = , and
2
1 −y 2
 
a(y, ϕ) = √ exp .
2πϕ 2ϕ

4.4.2 Scaled binomial distribution

We have
!
m
p(y | µ, m) = µmy (1 − µ)m(1−y)
my
! my
m µ
= (1 − µ)m
my 1−µ
!
m µ
 
= exp my log( ) + m log(1 − µ)
my 1−µ
!
m µ
  
= exp m y log( ) + log(1 − µ) .
my 1−µ

To obtain the canonical form we set

! !
µ m
θ = log , w = m, ϕ=1 and a(y, ϕ) =
1−µ my
and so then
  
p(y | µ, m) = a(y, ϕ) exp w yθ − b(θ) ,

where b(θ) = − log(1 − µ).

To determine b(θ) in terms of θ we need to express µ as a function of θ. We have
µ exp(θ)
 
θ = log if and only if µ = .
1−µ 1 + exp(θ)
54 CHAPTER 4. EXPONENTIAL DISPERSION MODELS

Therefore
exp(θ)
 
b(θ) = − log(1 − µ) = − log 1 −
1 + exp(θ)
1
   
= − log = log 1 + exp(θ) .
1 + exp(θ)

The scaled Binomial distribution is an example of an EDM where the dispersion parameter ϕ is
fixed at 1. Another such example is the Poisson distribution as the following exercise shows.

Exercise 9: canonical form

Derive the canonical form of the probability mass function for the Poisson distribution.
Solution: We have
µy
p(y | µ) = exp(−µ)
y!
1
= exp(y log(µ)) exp(−µ)
y!
1
 
= exp y log(µ) − µ .
y!

Now set θ = log(µ) and ϕ = 1, then

1
 
p(y | µ) = p(y | θ, ϕ = 1) = exp yθ − exp(θ)
y!
which is a density in canonical form with

b(θ) = exp(θ), and

1
a(y, ϕ) = a(y, ϕ = 1) = .
y!

4.5 Cumulants for EDMs

Recall the canonical form of an EDM
θy − b(θ)
 
p(y | θ, ϕ) = a(y, ϕ) exp .
ϕ

The function b(θ) is called the cumulant function.

The cumulant generating function of an EDM is
b(θ + tϕ) − b(θ)
K(t) = .
ϕ

Note that the recommended textbook by Dunn and Smyth uses κ(θ) (kappa) for the cumulant
function, whereas we use b(θ).
4.5. CUMULANTS FOR EDMS 55

We can obtain the mean and variance of an EDM from the derivatives of its cumulant generating
function:

b(θ + tϕ) − b(θ)

K(t) = ,
ϕ
ϕ b′ (θ + tϕ)
K ′ (t) = = b′ (θ + tϕ),
ϕ
K ′′ (t) = ϕb′′ (θ + tϕ).

Hence

E(Y | ϕ, θ) = K ′ (0) = b′ (θ)

Var(Y | ϕ, θ) = K ′′ (0) = ϕb′′ (θ)

Therefore, the mean is independent of ϕ and the variance is proportional to ϕ (hence the name
“dispersion parameter”.)

Next we prove that for EDMs the cumulant generating function is given by

b(θ + tϕ) − b(θ)

K(t) = .
ϕ

Proof Recall the canonical form of an EDM

θy − b(θ)
 
p(y | θ, ϕ) = a(y, ϕ) exp .
ϕ

Then,

 
M (t) = E exp(tY )

Z
= exp(ty) p(y |, θ, ϕ)dy
y∈S

θy − b(θ)
Z  
= a(y, ϕ) exp ty + dy
y∈S ϕ

(tϕ + θ)y − b(θ)

Z  
= a(y, ϕ) exp dy
y∈S ϕ
56 CHAPTER 4. EXPONENTIAL DISPERSION MODELS

Let θ∗ = θ + tϕ, then

 ∗
θ y − b(θ)
Z 
M (t) = a(y, ϕ) exp dy
y∈S ϕ

− b(θ∗ ) + b(θ∗ ) − b(θ)

 ∗
θ y
Z 
= a(y, ϕ) exp dy
y∈S ϕ

− b(θ∗ ) b(θ∗ ) − b(θ)

 ∗
θ y
Z   
= a(y, ϕ) exp exp dy
y∈S ϕ ϕ

b(θ∗ ) − b(θ) Z
 
= exp p(y | θ∗ , ϕ)dy
ϕ y∈S

b(θ∗ ) − b(θ)
 
= exp .
ϕ

Hence the cumulant generating function is given by

  b(θ∗ ) − b(θ) b(θ + tϕ) − b(θ)
K(t) = log M (t) = = .
ϕ ϕ

4.6 The canonical link

4.6.1 Definition
Recall that the mean and variance of an EDM can be derived from the cumulant function:
E(Y | θ) = b′ (θ)
Var(Y | ϕ, θ) = ϕb′′ (θ)
Let µ = E(Y | θ) = b′ (θ), then
dµ Var(Y | ϕ, θ)
= b′′ (θ) = > 0.
dθ ϕ
So µ is a strictly increasing function of θ (and vice versa) and thus, there is a one-to-one
correspondence between the canonical parameter θ and the mean µ.
For every EDM there is a function g that maps µ onto the canonical parameter θ
g(µ) = θ.
This is the canonical link function.
We can derive the canonical link from the cumulant function. The canonical mean function h(θ)
is the inverse of the canonical link function:
h(θ) = µ = b′ (θ).
So we invert h() to get
θ = g(µ).
4.7. EXAMPLES OF CANONICAL LINK FUNCTIONS 57

4.7 Examples of canonical link functions

4.7.1 The normal distribution
The cumulant function is
θ2
b(θ) = .
2
Hence
µ = b′ (θ) = θ.
Therefore, the canonical link g() for the normal distribution is the identity link, that is
g(µ) = µ.

4.7.2 The Poisson distribution

As shown in Exercise 9, the cumulant function for the Poisson distribution is

b(θ) = exp(θ).

Hence
µ = b′ (θ) = exp(θ).
Solving for θ gives
θ = log(µ).
Therefore, the canonical link g() for the Poisson distribution is the log link, that is g(µ) = log(µ).

4.7.3 The scaled Binomial distribution

The cumulant function is  
b(θ) = log 1 + exp(θ) .
Hence
exp(θ)
µ = b′ (θ) = .
1 + exp(θ)
Solving for θ gives
µ
 
θ = log .
1−µ
Therefore, the
 canonical
 link g() for the (scaled) Binomial distribution is the logit link, that is
g(µ) = log 1−µ .
µ
58 CHAPTER 4. EXPONENTIAL DISPERSION MODELS

Exercise 10: canonical link

Find the canonical link given the cumulant function for the following distributions:
• The gamma distribution

b(θ) = − log(−θ) for θ < 0.

• The negative binomial distribution

 
b(θ) = −k log 1 − exp(θ) for θ < 0.

• The inverse Gaussian distribution

√
b(θ) = − −2θ for θ < 0.

Solution:
Recall that θ ∈ Θ such that Θ = {θ ∈ R : |b(θ)| < ∞}. Note that in the examples above θ
has been restricted to lie on the negative real line. This is to ensure that |b(θ)| < ∞ for all
permissible θ.
Gamma distribution b(θ) = − log(−θ) and so b′ (θ) = − 1θ . Setting µ = b′ (θ) = − 1θ and solving
for θ we find θ = − µ1 and so the canonical link function of the gamma distribution is given by

1
g(µ) = − .
µ
 
Negative binomial distribution b(θ) = −k log 1 − exp(θ) and so b′ (θ) = k exp(θ)
1−exp(θ)
. Setting
µ = b′ (θ) and solving for θ we find
µ
 
θ = log
k+µ
 
and so the canonical link for the negative binomial distribution is g(µ) = log k+µ
µ
.
√
Inverse Gaussian distribution b(θ) = − −2θ and so b′ (θ) = √−2θ 1
. Setting µ = b′ (θ) and
solving for θ we find θ = − 2µ1 2 and so the canonical link function for the inverse Gaussian
distribution is g(µ) = − 2µ1 2 .

4.8 Deviance
In ST231 we derived our parameter estimates by minimizing the residual sum of squares function,
or deviance. n
= (yi − µi )2 .
X
D(y, µ)
i=1

We saw that, because the normal density is given by

(y − µ)2
 
p(y | µ, σ ) 2
∝ exp − ,
2σ 2
4.8. DEVIANCE 59

the likelihood
i=1 (yi − µi )2
 Pn 
p(y | µ, σ )
2
∝ exp −
2σ 2
is maximized when the deviance is minimized.
In week 1 of ST346 we extended this to weighted models:
n
= wi (yi − µi )2 .
X
D(y, µ)
i=1

The concept of deviance arises naturally from EDMs, but different EDMs will have different
formulae for the deviance. The one-to-one correspondence between θ and µ implies that we can
re-write the density function of an EDM in terms of µ, ϕ instead of θ, ϕ. Let

t(y, µ) = θy − b(θ),

then the canonical form of the EDM is given by

θy − b(θ)
 
p(y | θ, ϕ) = a(y, ϕ) exp
ϕ

t(y, µ)
 
= a(y, ϕ) exp
ϕ

t(y, y) t(y, µ) − t(y, y)

   
= a(y, ϕ) exp exp
ϕ ϕ
h i
 2 t(y, y) − t(y, µ) 
= a∗ (y, ϕ) exp −
2ϕ

 
where a∗ (y, ϕ) = a(y, ϕ) exp t(y, y)/ϕ .
Setting  
d(y, µ) = 2 t(y, y) − t(y, µ) .

gives the dispersion model form of the EDM:

d(y, µ)
 
p(y | µ, ϕ) = a (y, ϕ) exp
∗
−
2ϕ

where d(y, µ) is the unit deviance.

Proposition 4.1 (Unit deviance). The unit deviance d(y, µ) is non-negative and exactly zero if
and only if µ = y. It is thus is a measure of the discrepancy between the expected value µ and
the observed value y.
60 CHAPTER 4. EXPONENTIAL DISPERSION MODELS

Proof
If t(y, µ) has a unique maximum at µ = y, then the unit deviance
 
d(y, µ) = 2 t(y, y) − t(y, µ)

is zero at µ = y and positive otherwise.

Consider t as a function of θ. Then

dt
= y − b′ (θ) = y − µ
dθ
and so
dt
= 0 ⇐⇒ µ = y.
dθ
Moreover
d2 t
= −b′′ (θ) = −Var(Y )/ϕ < 0.
dθ2

Hence t is a strictly concave function of θ with a unique maximum at y = µ.

It follows that the unit deviance is non-negative and is exactly zero if and only if µ = y.

4.8.1 Example: normal unit deviance

Let’s derive the unit deviance for the normal distribution from the canonical form.
Recall that b(θ) = θ2 /2 and the canonical link is µ = θ. Then

t(y, µ) = yθ − b(θ)
= yθ − θ2 /2
= yµ − µ2 /2.

This is maximised at
t(y, y) = y 2 − y 2 /2 = y 2 /2.

The unit deviance for the normal distribution is

 
d(y, µ) = 2 t(y, y) − t(y, µ)
 
= 2 y /2 − yµ + µ /2
2 2

= y 2 − 2yµ + µ2
= (y − µ)2 .
4.8. DEVIANCE 61

4.8.2 Unit deviance on the boundary

If the parameter space for µ is bounded we need to take extra care.
If y lies on the boundary of the possible values of µ, it is possible that t(y, y) is not defined.
So we modify the definition of the unit deviance.
Bounded below:
   
d(y, µ) = lim 2 t(y, y + ϵ) − t(y, µ) = 2 lim t(y, y + ϵ) − t(y, µ) where ϵ > 0.
ϵ→0+ ϵ→0+

Bounded above:
   
d(y, µ) = lim 2 t(y, y − ϵ) − t(y, µ) = 2 lim t(y, y − ϵ) − t(y, µ) where ϵ > 0.
ϵ→0+ ϵ→0+

4.8.3 Example: Poisson unit deviance

Recall that the cumulant function is b(θ) = exp(θ) and the canonical link is log(µ) = θ.
Hence
t(y, µ) = yθ − b(θ) = y log(µ) − µ.

Note that t(y, y) is not defined for y = 0.

First case: if y > 0, then
 
d(y, µ) = 2 t(y, y) − t(y, µ)
 
= 2 y log(y) − y − y log(µ) + µ
 
= 2 y log(y/µ) − (y − µ)

So the unit deviance depends partly on the ratio y/µ and partly on the difference y − µ.
Second case: if y = 0, then
 
d(0, µ) = lim 2 t(0, ϵ) − t(0, µ)
ϵ→0+
 
= lim 2 0 × log(ϵ) − ϵ − 0 × log(µ) + µ)
ϵ→0+

= 2µ.
62 CHAPTER 4. EXPONENTIAL DISPERSION MODELS

4.8.4 Total deviance

Suppose we have independent Yi ∼ EDM(µi , ϕ/wi ) for i = 1, . . . , n, where w1 , . . . , wn are fixed
weights.
The total deviance is
n
= wi d(yi , µi ).
X
D(y, µ)
i=1

The scaled deviance is

D(y, µ)
D∗ (y, µ) = .
ϕ
• The total deviance D(y, µ) and the scaled deviance D∗ (y, µ) measure the discrepancy
between the observed values y1 , . . . , yn and the corresponding mean values predicted by
the model µ1 , . . . , µn .
• The smaller the deviance the better the fit. Hence the total deviance measures relative
goodness-of-fit of the model.
• Later we will see that we can compare the deviance of nested models and generalize the
Analysis of Variance (ANOVA) to an Analysis of Deviance.

Exercise 11 - Binomial unit deviance

Derive the unit deviance for the scaled binomial distribution. There are three cases to consider
1. y = r/m for r = 1, . . . , m − 1,
2. y = 0,
3. y = 1.
Solution: The canonical link function for the scaled binomial distribution is the logit link.
Therefore we have θ = log 1−µ
µ
and so exp(θ) = 1−µ
µ
.
Then the cumulant function is
µ 
  
b(θ) = log 1 + exp(θ) = log 1 + = − log(1 − µ).
1−µ

Therefore

t(y, µ) = yθ − b(θ)
!
µ
= y log + log(1 − µ)
1−µ
= y log(µ) + (1 − y) log(1 − µ).

Case 1: y = r/m for r = 1, . . . , m − 1 and so 0 < y < 1.

We have
t(y, y) = y log(y) + (1 − y) log(1 − y).
4.8. DEVIANCE 63

It follows that the unit deviance in this case is given by

 
d(y, µ) = 2 t(y, y) − t(y, µ)
 
= 2 y log(y) + (1 − y) log(1 − y) − y log(µ) − (1 − y) log(1 − µ)
1−y
! !
y

= 2 y log + (1 − y) log .
µ 1−µ

Case 2: y = 0.
For a scaled binomial distribution we have 0 < µ < 1. Hence, for y = 0, we need to replace
t(y, y) with limϵ→0+ t(y, y + ϵ).
We have
 
lim t(0, ϵ) = lim 0 × log(ϵ) + 1 × log(1 − ϵ) = log(1) = 0.
ϵ→0+ ϵ→0+

It follows that
 
d(0, µ) = 2 lim t(0, ϵ) − t(0, µ)
ϵ→0+
 
= 2 0 − 0 × log(µ) − 1 × log(1 − µ) = −2 log(1 − µ)

Case 3: y = 1.
For y = 1, we need to replace t(y, y) with limϵ→0+ t(y, y − ϵ).
We have
 
lim t(1, 1 − ϵ) = lim 1 × log(1 − ϵ) + 0 × log(ϵ) = log(1) = 0.
ϵ→0+ ϵ→0

It follows that
 
d(1, µ) = 2 lim t(1, 1 − ϵ) − t(1, µ)
ϵ→0
 
= 2 0 − 1 × log(µ) − 0 × log(1 − µ) = −2 log(µ).

In summary

−2 log(1 − µ) when y = 0,







  
=
  
d(y, µ)
 2 y log y
µ
+ (1 − y) log 1−y
1−µ
when y ∈ { m1 , . . . , m−1
m
}



−2 log(µ) when y = 1.



Chapter 5 Generalized linear models

5.1 The variance function

5.1.1 Introduction
Recall that the mean and variance of an EDM can be derived from the cumulant function b(θ).

E(Y | ϕ, θ) = b′ (θ) and Var(Y | ϕ, θ) = ϕb′′ (θ).

Furthermore, there is a one-to-one correspondence between the canonical parameter θ and the
mean parameter µ = E(Y ) via the canonical link. Hence we can rewrite the variance in terms
of µ and ϕ as
Var(Y | µ, ϕ) = ϕV (µ)
where V (µ) is the variance function.

5.1.2 Properties of the variance function

Proposition 5.1. An EDM is uniquely defined by the relation

Var(Y | µ, ϕ) = ϕV (µ).

To prove the proposition we need the following lemmas.

Lemma 5.1. Z
1
θ(µ) = dµ.
V (µ)

Proof. (Lemma 5.1) Recall that µ = b′ (θ) and so

dµ
= b′′ (θ) = V (µ)
dθ

dθ 1
=⇒ =
dµ V (µ)
Z
1
=⇒ θ = dµ + C
V (µ)

NB: Integrals are indefinite, so we can add an arbitrary constant C to θ.

64
5.1. THE VARIANCE FUNCTION 65

Lemma 5.2.
µ
  Z
b θ(µ) = dµ.
V (µ)

Proof. (Lemma 5.2) We have b′ (θ) = µ and so

Z Z
dθ Z
µ
b(θ) = µ(θ)dθ = µ(θ) dµ = dµ + C
dµ V (µ)

NB: Integrals are indefinite, so we can add an arbitrary constant C to b(θ).

Next, the proof of Proposition 5.1.

Proof. (Proposition 5.1) Lemmas 5.1 and 5.2 show that we can construct the canonical parameter
θ and the cumulant function b(θ) from the variance function. As we have seen earlier, the
cumulant generating function is then defined as

b(θ + tϕ) − b(θ))

K(t) = .
ϕ
For EDMs, if the cumulant generating function exists then it uniquely determines the distribution.
NB: Any constant terms added to b(θ) cancel out.

Proposition 5.1 means that when modelling real-world data, we need to check the mean-variance
relationship. This will determine which EDM to use.

5.1.3 Example - normal distribution

For the normal distribution V (µ) = 1. Hence
Z
1 Z
θ(µ) = dµ = 1dµ = µ
V (µ)

and

Z
µ Z
µ2
b(θ(µ)) = dµ = µdµ = .
V (µ) 2

Substituting µ = θ gives
θ2
b(θ) = .
2

Exercise 12 - variance function

Derive the cumulant function for the scaled binomial and the Poisson distributions from their
variance functions. Recall
a. the scaled binomial variance function is V (µ) = µ(1 − µ);
b. the Poisson variance function is V (µ) = µ.
66 CHAPTER 5. GENERALIZED LINEAR MODELS

Solution
a. For the scaled binomial distribution we have
Z
1 Z
1
θ = dµ = dµ.
V (µ) µ(1 − µ)

Next we apply a change of variable, setting z = log(µ/(1 − µ)). Then

exp(z) dµ exp(z)
µ = and = h i2 = µ(1 − µ).
1 + exp(z) dz 1 + exp(z)

Therefore
Z
1 Z
1 dµ Z
1
θ = dµ = dz = µ(1 − µ)dz = z.
µ(1 − µ) µ(1 − µ) dz µ(1 − µ)

As z = log(µ/(1 − µ)) we conclude that θ = log(µ/(1 − µ)). Next,

Z
µ Z
µ Z
1
b(θ) = dµ = dµ = dµ = − log(1 − µ).
V (µ) µ(1 − µ) (1 − µ)

To express the cumulant function as a function of θ we use that θ = log(µ/(1 − µ)). Then
µ = 1+exp(θ)
exp(θ)
and so

1
!
b(θ) = − log(1 − µ) = − log = log (1 + exp(θ)) .
1 + exp(θ)

b. For the Poisson distribution we have

Z
1 Z
1
θ = dµ = dµ = log(µ).
V (µ) µ

Next, Z
µ Z
µ Z
b(θ) = dµ = dµ = 1dµ = µ.
V (µ) µ
Thus, with θ = log(µ), we have

b(θ) = µ = exp(θ).

5.2 Quasi-likelihood
5.2.1 Introduction
It is possible that an EDM does not exist for a given combination of dispersion parameter ϕ
and variance function V (µ). For example, the Poisson and binomial distributions assume ϕ = 1.
But if we allow ϕ to be free, then there is no EDM for ϕ ̸= 1. Surprisingly, this does not stop us.
5.3. THE GAMMA DISTRIBUTION 67

We can define a quasi-likelihood:

• Quasi-Poisson: ϕ free, V (µ) = µ.
• Quasi-binomial: ϕ free, V (µ) = µ(1 − µ).
Later we will see that the maximum likelihood estimator and its large sample properties depend
only on:
• the dispersion parameter ϕ,
• the variance function V (µ), and
• the link function g(µ).
So we can still conduct inference for quasi-likelihoods even when there is no corresponding EDM.

5.2.2 Example - Noisy miner data

A study1 examined the effect of the density of eucalyptus trees on the invasion of woodland
habitats by noisy miners, a small but aggressive Australian bird. The data from the study is
available as the dataframe nminer in the GLMsData package and discussed in Example 5.9 in
the textbook by Dunn and Smyth.2
The noisy miner data is an example of over-dispersed count data, where the Poisson variance
function V (µ) = µ fits the data, but with a dispersion parameter ϕ > 1.

5.3 The gamma distribution

5.3.1 Properties
Let Y have a gamma distribution with shape s > 0 and rate r > 0, that is,

Y ∼ Γ(s, r),

then Y has probability density function

exp(−ry)rs y s−1
p(y | r, s) = for y > 0.
Γ(s)

The graph below shows the gamma density for different values of the shape parameter but with
a fixed mean.
• For s = 1, the gamma distribution is the exponential distribution.
• For s > 1, the density is unimodal and becomes increasingly symmetric as s → ∞.
• For s < 1, the density tends to ∞ as y → 0.
1
M. Maron (2007): Threshold effect of eucalypt density on an aggressive avian competitor. Biological
Conservation, 136, 100–107.
2
Dunn, P. K. and Smyth, G.K (2018): Generalized linear models with examples in R Vol. 53. New York:
Springer.
68 CHAPTER 5. GENERALIZED LINEAR MODELS

2.0
s = 1/2
s=1
s=2

1.5
s=4

density

1.0
0.5
0.0

0.0 0.5 1.0 1.5 2.0 2.5 3.0

The gamma distribution is an EDM with

s 1
µ = and ϕ = .
r s
The density can be rewritten as

1 1
1
y y
  
ϕ
p(y | µ, ϕ) = exp −
ϕµ y ϕµ Γ( ϕ1 )

[You are not expected to memorize this density.]

The cumulant function of the gamma distribution is

b(θ) = − log(−θ).

This yields mean µ and variance function V (µ) as

1
µ = b′ (θ) = − ,
θ
1
V (µ) = b′′ (θ) = = µ2 .
θ2

The canonical link function is

1
g(µ) = − .
µ
The default link for the gamma distribution in R is slightly different as it omits the minus sign:
1
g(µ) = .
µ

5.3.2 Implementation in R
Use the argument family=Gamma in the command glm()! Note the capital letter G.
5.4. THE INVERSE GAUSSIAN DISTRIBUTION 69

In practice we do not use the canonical link. The canonical link function only permits negative
values for the linear predictor η = xT β and hence imposes a boundary on the parameter space
for the coefficients β. In contrast, the log link does not impose such a boundary, is more
numerically stable and makes it easier to interpret parameters.
For a log link function use family=Gamma(link="log").

5.3.3 Example: the lime dataset

The lime data set in the package GLMsData concerns measurements of leaf biomass from small-
leaved lime trees (Tilia cordata). See Sections 11.2 and 11.3 in the textbook by Dunn and
Smyth.3

5.4 The inverse Gaussian distribution

5.4.1 Properties
Let Y have an inverse Gaussian distribution with parameters µ and ϕ, that is

Y ∼ IG(µ, ϕ),

then Y has density function

1 (y − µ)2
s  
p(y | µ, ϕ) = exp −
2πϕy 3 2ϕµ2 y

for y > 0, µ > 0 and ϕ > 0.

The inverse Gaussian is an EDM with
1 √
θ = − and b(θ) = − −2θ.
2µ2

Its variance function is

V (µ) = µ3 .

The inverse Gaussian distribution is useful for continuous distributions on the positive real line
where the variance increases very rapidly with the mean.

5.4.2 Relationship to Brownian Motion

Consider a stochastic process with

X0 = 0,
t
Xt = + σWt .
µ
3
Dunn, P. K. and Smyth, G.K (2018): Generalized linear models with examples in R Vol. 53. New York:
Springer.
70 CHAPTER 5. GENERALIZED LINEAR MODELS

where Wt is a Brownian motion (or Wiener process). This is a Brownian motion with drift µ1 . A
random walk is superimposed on top of a systematic tendency for the process {Xt } to increase
with time t. The hitting time for the process to hit fixed level Xt = 1 has an inverse Gaussian
distribution:
T = inf{t > 0 | Xt = 1} ∼ IG(µ, σ 2 ).

The plot below shows three realisations of the Brownian motion with drift, using the same
parameters for µ and ϕ. The dots show the hitting times for the three realisations.
Brownian Motion

Time

The canonical link is

1
g(µ) = − .
2µ2
The default link for the inverse Gaussian in R is slightly different as it omits the minus sign

1
g(µ) = .
2µ2

In practice it is rarely used. The log link is numerically more stable and makes it easier to
interpret parameters.

5.4.3 Example: the perm dataset

The perm dataset in the R package GLMsData concerns measurements of permeability of building
materials.
• Permeability is a measure of the ease of passage of liquids or gases through the material.
• The Brownian motion with drift is a useful physical model for permeability and hence the
inverse Gaussian is a useful statistical model for the results of the permeability tests.
See Sections 11.4 and 11.7.1 in the recommended textbook by Dunn and Smyth.4
4
Dunn, P. K. and Smyth, G.K (2018): Generalized linear models with examples in R Vol. 53. New York:
Springer.
5.5. THE POISSON LIMIT OF THE BINOMIAL DISTRIBUTION 71

5.5 The Poisson limit of the binomial distribution

5.5.1 Binomial distribution for rare outcomes

Suppose that Y has an unscaled binomial distribution

Y ∼ Binomial(m, µ).

Then Y has expectation E(Y ) = mµ.

Suppose we increase the size parameter m but keep E(Y ) constant. What happens to the
distribution of Y ?

Example: Consider E(Y ) = 5.

• The light bars show the probability mass function (pmf) for the binomial distribution as
the size parameter m increases.

• The dark bars show the pmf of the Poisson distribution.

Size m = 10 Size m = 30 Size m = 200

0.20
0.25

0.15
0.20
0.15

0.15
0.10
pmf

pmf

pmf

0.10

0.10

0.05 0.05
0.05

0.00 0.00 0.00

0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
x x x

More formally, suppose

Y ∼ Binomial(m, µ)

with µ = λ/m so that

E(Y ) = mµ = λ.

Then the limiting distribution of Y as m → ∞ with λ fixed is

Y ∼ Poisson(λ).

Proof Consider the moment generating function (mgf) of the (unscaled) binomial distribution
72 CHAPTER 5. GENERALIZED LINEAR MODELS

 
M (t) = E exp(tY )
m
!
m r
= exp(rt) µ (1 − µ)m−r
X

r=0 r
m
! r
m
= exp(t)µ (1 − µ)m−r
X

r=0 r
 m
= exp(t)µ + 1 − µ
 h im
= 1 + µ exp(t) − 1 .

Substituting µ = λ/m we obtain

h i
 λ exp(t) − 1 m
M (t) = 1+
m
m
a

= 1+ for a = λ[exp(t) − 1].
m

Taking the limit m → ∞ with t, λ fixed we obtain

 
M (t) −→ exp(a) = exp λ[exp(t) − 1]

which is the mgf of the Poisson distribution with mean λ.

5.6 Choosing between the Poisson and Binomial GLM

We have seen that we can use a Poisson GLM to model rates and a Binomial GLM to model
proportions. As the Binomial distribution converges to the Poisson distribution, for rare events
we have a choice between the two models. Prof Martyn Plummer designed the following flowchart
to help you decide. In the flow chart Y is the count response variable and D the denominator
which turns the count into a rate or a proportion.
5.6. CHOOSING BETWEEN THE POISSON AND BINOMIAL GLM 73

Exercise 13 - GLM scenarios

Consider the following scenarios. Determine which model might be appropriate to use in each
case. Give the exponential dispersion model, the link function, the outcome variable, and the
predictor variables.
a. A cohort of 18 year-old school leavers were surveyed regarding their future plans, with
specific interest in higher education. Several variables were measured, including age,
gender, smoking status, as well as whether or not they had a place at university or other
higher education institute. We wish to investigate the variables which are associated with
higher education attendance.
b. Sixty rats of the same age were divided into 2 groups of 30. At the start of the experiment,
all animals were weighed. Then one group was fed a control diet, whilst the other was fed
a diet supplemented with vitamin D. After 6 weeks, all rats were weighed again. We wish
to understand how vitamin D supplementation affects weight gain.
c. A car manufacturer conducted a study to investigate the reliability of their cars. They
measured the number of times that each of 1000 vehicles had broken down in the 10 years
since it was made. They also measured the number of miles that each car had been driven
as well as the model of each vehicle and the number of times that it had been serviced.
Solution:
a. The outcome variable is the binary indicator of a university place, therefore the family
must be binomial. Predictor variables are age, gender, and smoking status. A number of
link functions are compatible with the binomial family (logit, probit, cloglog, cauchit, and
log).
b. There are several possible correct answers. The fundamental choice is the outcome variable,
which could be:
• final weight (possibly with a log transformation);
• difference between final weight and initial weight;
74 CHAPTER 5. GENERALIZED LINEAR MODELS

• ratio between the final weight and initial weight (possibly with a log transformation,
in which case this is also the difference of the logs)
The family could be gamma, since the outcome is positive and continuous, with a log or
identity link. The family could also be gaussian with an identity link, particularly if the
response variable is log-transformed.
If the outcome is final weight, then the initial weight needs to be a predictor. The treatment
group (vitamin D supplementation versus control diet) should always be a predictor.
c. The outcome variable is the number of times a given car has broken down in 10 years. As
this is count data we may consider a Poisson model. The predictor variables are: mileage,
model, and service history. An offset is not required as the cars all have the same follow-up
time (10 years). The Poisson family is compatible with the log link, square root link, and
identity link.
Chapter 6 Maximum likelihood estimation for
GLMs

6.1 Review of maximum likelihood estimation

6.1.1 Likelihood
In the following we assume suitable regularity conditions to be satisfied by the models
considered in this module.
Suppose we observe independent random variables Y1 , . . . , Yn , where the pdf of Yi is

pi (y | β) for β ∈ Rp ,

which is assumed to be a member of the exponential family of distributions.

• Y1 , . . . , Yn are independent but not identically distributed.
• The distribution of Yi is parameterized by β.
The likelihood L(β) is the joint pdf considered as a function of the parameters:
n
L : Rp → R with pi (yi | β)
Y
β 7→
i=1

We normally work in terms of the log likelihood, as the log likelihood is the sum of individual
contributions from independent observations
  n  
= log L(β | y) = log pi (yi | β)
X
l(β | y)
i=1

and this is usually more convenient than taking products.

Note The above also applies to discrete distributions where pi (yi | β) is a pmf rather than a
pdf.
The likelihood is a relative measure of consistency between the parameters β and the data y.
• l(β | y) is defined up to an additive constant.
• Differences in log likelihood are always well defined.
• We may omit terms that are constant when deriving the log likelihood from the pdf/pmf.
(NB: what is constant may depend on the context!)
Suppose β (1) and β (2) are two candidate values for the unknown parameter β. If

l(β (1) | y) − l(β (2) | y) > 0,

then β (1) has more support from the data y than β (2) .

75
76 CHAPTER 6. MAXIMUM LIKELIHOOD ESTIMATION FOR GLMS

(Under suitable regularity conditions) the maximum likelihood estimate β

b satisfies

b | y) − l(β | y)
l(β > 0 ∀β ̸= β,
b

so has the highest support from the data among all possible parameter values.
The score function U : Rp → Rp is the first derivative of the log likelihood:
∂l(β | y)
U (β | y) = .
∂β

The maximum likelihood estimate β

b satisfies the score equations, that is

U (β
b | y) = 0.

Score equations generalize the normal equations for linear models.

The score function may be viewed as a random vector by replacing the observed data y1 , . . . , yn
with the corresponding random variables Y1 , . . . , Yn . This random vector has expectation zero:
 
E U (β | Y ) = 0.

We say that the score function is an unbiased estimating function.

The variance of the score function is called the Fisher (expected) information matrix:
   
I(β) = Var U (β | Y ) = E U (β | Y ) U (β | Y ) T
.

The Fisher information (or expected information) matrix I(β) is positive semi-definite, that is

aT I(β) a ≥ 0 for any a ∈ Rp .

It can be shown that

∂ 2 l(β | Y )
   
I(β) = E − = E J(β | Y ) ,
∂β∂β T

where J(β | y) is the observed information matrix defined as

∂ 2 l(β | y)
J(β | y) = − .
∂β∂β T

Exercise 14 - Poisson maximum likelihood estimation

Suppose y1 , . . . , yn is an iid sample from a Poisson distribution with mean µ.
a. Derive the likelihood function, score function and expectation of the score function.
b. Determine an expression for the observed information and for the Fisher information.
Solution
6.1. REVIEW OF MAXIMUM LIKELIHOOD ESTIMATION 77

a. Let L denote the likelihood function and U the score function, then for a Poisson distribu-
tion with mean µ we have
n
µyi −µ
=
Y
L(µ | y) e ,
i=1 yi !
n
= log(L(µ|y)) = yi log(µ) − µ − log(yi !)
X
l(µ | y)
i=1
n
hX i
= log(µ) yi − nµ + constant,
i=1

n P
dl(µ|y) i=1 yi
U (µ | y) = = − n,
dµ µ
1 n
nµ
  X 
E U (µ | Y ) = E Yi − n = −n = 0.
µ i=1 µ

b. Let J denote the observed information and I the Fisher information, then
dU (µ | y)
Pn
d2 l(µ | y) yi
J(µ | y) = − = − = i=1
,
dµ2 dµ µ2
1 n
n
  X 
I(µ) = E J(µ | y) = E Yi = .
µ2 i=1 µ

6.1.2 Properties of the maximum likelihood estimator

Suppose the number of parameters p is fixed as n → ∞. Assuming suitable regularity conditions,
the maximum likelihood estimator βb has the following properties:

1. Consistency:
p
As n → ∞ we have β
b → β.
2. Invariance under reparameterisation:
Suppose γ is an alternative parameterization to β. Then for some invertible function s
we have
γ = s(β) and β = s−1 (γ).
The maximum likelihood estimates then satisfy
γb = s(β)
b and b = s−1 (γ).
β b

3. Asymptotic unbiasedness:
As n → ∞,
√  b 
n E(β) − β → 0.

4. Asymptotic Efficiency (Generalization of Gauss-Markov theorem):

b is the unique asymptotically unbiased estimator with minimum variance.
β
78 CHAPTER 6. MAXIMUM LIKELIHOOD ESTIMATION FOR GLMS

5. Asymptotic normality:
For sufficiently large n we can use the approximation:

β
b ∼ N (β, I(β)−1 ).

For further details see Sections 4.4 - 4.9 in the recommended textbook by Dunn and Smyth.1

6.2 Maximum likelihood for GLMs

6.2.1 Recap
Note: to simplify notation we omit the explicit conditioning on y and Y but this is still assumed.
The maximum likelihood estimate β
b solves the score equations

U (β)
b = 0.

We need an expression for the score function U (β) for GLMs. We also need an expression for
the Fisher information matrix for GLMs, that is
∂ 2 l(β)
 
I(β) = E −
∂β∂β T
so that we can calculate standard errors using the large sample approximation
 
β
b ∼ N β, I(β)−1 .

6.2.2 Log likelihood for GLMs

Definition 6.1 (Generalized linear model). A generalized linear model for outcomes
Y1 , . . . , Yn and predictor variables x1 , . . . , xn is defined by a combination of an exponential
dispersion model and a link function

Yi ∼ EDM(µi , ϕ/wi )
g(µi ) = xTi β

where E(Yi ) = µi = g −1 (xTi β).

There is a common dispersion parameter ϕ which is modified by individual prior weights
(w1 , . . . , wn ).

Recall from Section 4.3 that Yi as defined above has pdf/pmf given by
h i
w
i yi θi − b(θi ) 
p(yi |θi , ϕ) = a(yi , ϕ/wi ) exp .
ϕ
1
Dunn, P. K. and Smyth, G.K (2018): Generalized linear models with examples in R Vol. 53. New York:
Springer.
6.2. MAXIMUM LIKELIHOOD FOR GLMS 79

Also recall that the mean parameter µi for observation i can be mapped onto the canonical
parameter θi using the canonical link function. This allows us to write the log likelihood in
canonical form:
n
l(β, ϕ) = li (β, ϕ)
X

i=1
h i
n w
i θi yi − b(θi )  
= + log a(yi , ϕ/wi )
X

i=1 ϕ

where implicitly θi = θi (β) is a function of β. We can use this to derive the score function.

Lemma 6.1. h i
n wi yi − µi ∂θi
U (β) =
X
.
i=1 ϕ ∂β

Proof of Lemma 6.1

n n
∂li (θi )
U (β) = Ui (β) =
X X

i=1 i=1 ∂β

n
∂li (θi ) ∂θi
=
X

i=1 ∂θi ∂β
h i
wi yi − b′ (θi ) ∂θi
n
=
X

i=1 ϕ ∂β
h i
wi yi − µi ∂θi
n
=
X
.
i=1 ϕ ∂β

We consider two distinct cases:

• Canonical link. The score function and information matrix take a particularly simple
form.
• General link. The score function is more complex but can be expressed in terms of µi
and ϕ.

[Link] Overview on key results

If g() is the canonical link function, then

∂θi
= xi ,
∂β
h i
n w i y i − µi
U (β) =
X
xi ,
i=1 ϕ
n
wi V (µi )
= xi xTi .
X
I(β)
i=1 ϕ
80 CHAPTER 6. MAXIMUM LIKELIHOOD ESTIMATION FOR GLMS

If g() is a general link function, then

∂θi xi
= ,
∂β g ′ (µ i )V (µi )

h i
n wi yi − µi xi n 
Wi g ′ (µi ) 
 
U (β) = =
X X
yi − µi xi ,
i=1 ϕ g ′ (µi )V (µi ) i=1 ϕ

n n 
wi Wi

= x xT = xi xTi ,
X X
I(β) 2 V (µ ) i i
i=1 ϕ[g ′ (µ i )] i i=1 ϕ

wi
where Wi = is the ith working weight.
[g ′ (µ i )] V (µi )
2

[Link] Canonical link

Proposition 6.1 (Score function for GLM with canonical link). If g() is the canonical link
function, then

∂θi
= xi , and
∂β
h i h i
n wi yi − µi ∂θi n w i y i − µi x i
U (β) = =
X X
.
i=1 ϕ ∂β i=1 ϕ

Proof of Proposition 6.1

Recall that the canonical link maps the mean parameter onto the canonical parameter, that is
g(µ) = θ. Hence

θi = g(µi ) = xTi β, and so

∂θi
= xi .
∂β

Therefore the score function is

h i h i
n wi yi − µi ∂θi n wi yi − µi
U (β) = =
X X
xi .
i=1 ϕ ∂β i=1 ϕ

Note:
• The maximum likelihood estimate β
b solves the score equation U (β)
b = 0 which we can
write as
n h i
wi yi − µi (β) =
X
b x
i 0
i=1

independent of ϕ.
6.2. MAXIMUM LIKELIHOOD FOR GLMS 81

• Suppose the model has an intercept term, then the first column of the design matrix X is
a column of ones (xi1 = 1 ∀i). Writing µb i = µi (β),
b the score equation for column 1 solves

n h i
= 0.
X
wi yi − µb i
i=1

This generalizes the result for linear models with an intercept term that the (weighted)
sum of residuals is zero.

Proposition 6.2 (Fisher information for GLM with canonical link).

n
∂ 2 li (β) n
wi V (µi )xi xTi
 
= = (6.1)
X X
I(β) E − T .
i=1 ∂β∂β i=1 ϕ

Proof of Proposition 6.2

The Fisher information I(β) = E(J(β)), where J(β) is the observed information. Using the
fact that µi = b′ (θi ) we have
h i
n wi yi − b′ (θi )
U (β) =
X
xi .
i=1 ϕ

Recall that ∂θi

∂β
= xi . Then,

n
∂ 2 li (β) n
∂Ui
= =
X X
J(β) − T − T
i=1 ∂β∂β i=1 ∂β

n
∂Ui ∂θi
=
X
− T
i=1 ∂θi ∂β

h i
n
∂ wi yi − µi xi ∂θi
=
X
−
i=1 ∂θi ϕ ∂β T
h i
∂ wi yi − b (θi ) xi T
n ′
=
X
− xi
i=1 ∂θi ϕ

n
wi b′′ (θi )xi xTi
=
X

i=1 ϕ

n
wi V (µi )xi xTi
=
X

i=1 ϕ

This does not depend on Y1 , . . . , Yn , hence I(β) = J(β).

82 CHAPTER 6. MAXIMUM LIKELIHOOD ESTIMATION FOR GLMS

Exercise 15 - Score function for weighted normal linear model

Using the canonical form of the normal density, derive an expression for the score function and
the Fisher information for a weighted normal linear model.
Solution: The canonical form for a N (µi , wϕi ) density is given by
h i
ϕ
w
i θi yi − 12 θi2 
p(yi | θ, ) = a(yi , ϕ/wi ) exp .
wi ϕ

Recall that the canonical link function of the normal distribution is the identity link. Therefore

θi = θi (β) = µi (β) = xTi β.

It follows that ∂θi

∂β
= xi .
From the above density we derive the log-likelihood function as
n
wi h 1 i
= θi yi − θi2 .
X
l(θ | y)
i=1 ϕ 2

Differentiating with respect to β we obtain

n n
∂l(θ | y) wi h i ∂θ
i wi h i
U (θ|y) = = =
X X
yi − θi yi − θi xi .
∂β i=1 ϕ ∂β i=1 ϕ

Compare this to the expression in Proposition 6.1. Note that for a normal linear model with
wi = 1 for all i = 1, . . . , n:

1X n h i 1X n h i
U (θ|y) = yi − θi xi = yi − xTi β xi .
ϕ i=1 ϕ i=1

Then setting the score function to zero we obtain the normal equations we derived in
ST221/ST231:
n h i
yi − xTi β xi =
X
0
i=1

or in matrix form
(y − Xβ)T X = 0
which is solved by setting β = (X T X)−1 X T y.
Next, to obtain the Fisher information, we take the expectation of minus the Hessian of the
score function, which gives

∂U (θ|y) n
wi ∂θi n
wi
  X 
= = = = xi xTi .
X
I(β) E(J(β)) −E −E xi T
∂β T i=1 ϕ ∂β i=1 ϕ

Compare this to the expression in Proposition 6.2 using the fact that V (µ) = 1 for a (weighted)
normal linear model.
6.2. MAXIMUM LIKELIHOOD FOR GLMS 83

[Link] General link function

Define the working weights W1 , . . . , Wn as
wi
Wi = .
[g ′ (µ i )] V (µi )
2

The working weights should not be confused with the prior weights w1 , . . . , wn defined by us
when we fit the model.

Proposition 6.3 (Score function for GLM with a general link). For a general link function g()
we have
∂θi xi
= and
∂β g ′ (µi )V (µi )
n 
Wi g ′ (µi ) 
 
U (β) =
X
y i − µi x i
i=1 ϕ

Proof of Proposition 6.3 We have

∂θi ∂θi ∂µi

= .
∂β ∂µi ∂β

As µi = b′ (θi ) we have
∂µi
= b′′ (θi ) = V (µi )
∂θi
and so
∂θi 1
= .
∂µi V (µi )

Furthermore g(µi ) = β T xi . Hence with the chain rule

∂µi
g ′ (µi ) = xi
∂β

and so
∂µi xi
= .
∂β g (µ
′
i)

Therefore it follows that

∂θi ∂θi ∂µi 1 xi
= = × ′ .
∂β ∂µi ∂β V (µi ) g (µi )

Next we derive the expression for the score function.

84 CHAPTER 6. MAXIMUM LIKELIHOOD ESTIMATION FOR GLMS

h i
n wi yi − µi ∂θi
U (β) =
X

i=1 ϕ ∂β
h i
n wi yi − µi 1 xi
=
X

i=1 ϕ V (µi ) g ′ (µi )

n 
Wi g ′ (µi ) 
 
=
X
yi − µi xi .
i=1 ϕ
where we used the working weights defined earlier as
wi
Wi = .
[g (µi )]2 V (µi )
′

Proposition 6.4 (Fisher information for a GLM with a general link). For a general link function
g() the Fisher information matrix is
n 
Wi

= xi xTi .
X
I(β)
i=1 ϕ

Proof of Proposition 6.4 (non-examinable): Again, we will use Lemma 6.1, that is
h i
n wi yi − b′ (θi ) ∂θi
U (β) =
X
.
i=1 ϕ ∂β
We first consider the observed information
h i
∂U (β) ∂ n
X wi yi − b′ (θi ) ∂θi 
J(β) = − = −
∂β T ∂β T i=1 ϕ ∂β
h i
n w y − b′ (θ )
n h
wi b′′ (θi ) ∂θi i ∂θi i i i ∂ 2 θi
=
X X
−
i=1 ϕ ∂β ∂β T i=1 ϕ ∂β∂β T
where the last step follows by the product rule of differentiation. Now, in the above result the
first term does not depend on y. The second term, when we replace yi with Yi , has expectation
zero as E(Yi ) = µi = b′ (θi ). Hence
n
wi b′′ (θi ) ∂θi ∂θi
=
X
I(β)
i=1 ϕ ∂β ∂β T
Now b′′ (θi ) = V (µi ) and by Proposition 6.3
∂θi xi
= .
∂β g ′ (µ i )V (µi )

This allows us to write the Fisher (or expected) information matrix in terms of µi :
n
wi
= xi xTi .
X
I(β)
i=1 ϕ[g ′ (µ i )] V (µi )
2
6.3. NUMERICAL MAXIMUM LIKELIHOOD ESTIMATION FOR GLMS 85

6.3 Numerical maximum likelihood estimation for GLMs

6.3.1 The Newton-Raphson algorithm
Under suitable regularity conditions we find the maximum likelihood estimate γb of a scalar
parameter γ by solving U (γ) = 0. With an initial guess γe , a first order Taylor expansion of the
score function U gives
h i
U (γb ) ≈ U (γe ) + U ′ (γe ) γb − γe
h i
= U (γe ) − J(γe ) γb − γe

as U ′ (γe ) = −J(γe ) where J() is the observed information.

Now, using the fact that U (γb ) = 0, we have
h i−1
γb ≈ γe + J(γe ) U (γe ).

To determine γb we can now compute an iterative sequence of approximations

h i−1
γ (k+1) ≈ γ (k) + J(γ (k) ) U (γ (k) ), k = 0, 1, 2, . . . ,

until convergence is reached. This is the so-called Newton-Raphson algorithm.

If we approximate the observed information by the Fisher information and so compute the
sequence h i−1
γ (k+1) ≈ γ (k) + I(γ (k) ) U (γ (k) ), k = 0, 1, 2, . . . ,
then this algorithm is referred to as Fisher scoring.
For GLMs, it turns out that the iterations of the Fisher scoring algorithm can be performed by
solving the score equation of a weighted normal linear model. The algorithm is therefore known
as the IWLS (iterated weighted least squares) algorithm.

6.3.2 The IWLS algorithm for GLMs

Maximum likelihood (or quasi-likelihood) estimates for GLMs are obtained from the IWLS
algorithm. We take a local linear approximation to reduce a GLM to a linear model:
(0)
1. Start with an initial estimate β̃ .
Then, for k = 0, . . . ,
(k)
2. Take a linear approximation for β “close” to β̃ :
2.1 Approximate the likelihood using a weighted linear model.
(k+1)
2.2 Obtain new estimate β̃ of β from this linear model.
3. Continue with step 2 until convergence.
86 CHAPTER 6. MAXIMUM LIKELIHOOD ESTIMATION FOR GLMS

(k)
Let β̃ be our current estimate of β. Using first order Taylor approximation we can approximate
(k)
the score function in the neighbourhood of β̃ :
(k)
(k) ∂U (β̃ ) h (k)
i
U (β) ≈ U (β̃ )+ β − β̃
∂β T
(k) (k) (k)
= U (β̃ ) − J(β̃ )(β − β̃ )

(k) (k) (k)

h i
≈ U (β̃ ) − I(β̃ ) β − β̃ .

We previously derived that

1X n  
U (β) = Wi g ′ (µi ) yi − µi xi
ϕ i=1
1X n
I(β) = Wi xi xTi ,
ϕ i=1
wi
where Wi = .
[g (µi )]2 V (µi )
′

Hence
(k) (k) (k)
h i
U (β) ≈ U (β̃ ) − I(β̃ ) β − β̃

n g(k) n g(k)
Wi   Wi h (k) i
= g (µ
′
fi ) yi − µ
(k) (k)
xi xTi β − β
X X
fi xi − e
i=1 ϕ i=1 ϕ

n g(k) 
Wi

e (k)
 
= g (µ
′
fi ) yi − µ
(k) (k)
+ xTi β xTi β
X
fi − xi
i=1 ϕ

n g(k) 
Wi

= zei (k) xTi β
X
− xi
i=1 ϕ
where
(k) wi
W
g
i =
[g ′ (µ
fi )]2 V (µ
(k)
fi (k) )
is the current working weight and
  (k)
zei (k) = g ′ (µ
fi (k) ) yi − µ
fi (k) + xTi β
e

is the current working observation for observation i.

We solve the approximate score equations
n g(k) 
Wi

U (β) zei (k) − xTi β xi =
X
≈ 0.
i=1 ϕ

This is equivalent to estimating β for a linear model on the working observations, that is
ϕ
 
T
Z
f
i ∼ N x i β, .
g(k)
W i
6.3. NUMERICAL MAXIMUM LIKELIHOOD ESTIMATION FOR GLMS 87

The maximum likelihood estimate for this weighted linear model solves the approximate score
equation and is given by
 −1
(k+1) (k) (k) (k) (k) (k)
β
e = XT W
g X XT W
g zei (k) where W
g = diag(W
g
1 ,...,W
g
n ).

(k+1)
This new estimate β
e becomes our new value of βe for the next iteration. We continue until
e (k+1) is the same as the previous one, that is β
the new estimate β e (k) .

e (k+1) solves the approximate score equations and so

Then β

n g(k) 
Wi

(k+1)
zei (k) − xTi β =
X
e xi 0.
i=1 ϕ

But then
(k+1)   (k) (k+1)
zei (k) − xTi β
e = g ′ (µ
fi (k) ) yi − µ
fi (k) + xTi β
e − xTi β
e

e (k) = β
e (k+1) .
 
= fi (k) ) yi − µ
g ′ (µ fi (k) as β

e (k+1) also solves the exact score equations:

Thus β

e (k+1) ) 1X n
g(k) g ′ (µ
 
U (β = Wi fi (k) ) yi − µ
fi (k) xi = 0.
ϕ i=1

e (k+1) is equal to the maximum likelihood estimate.

Therefore β
The information matrix from our approximate linear model is equal to

1 Xn
g(k) x xT
Wi i i
ϕ i=1

and thus exactly equal to I(βe(k) ), the Fisher information matrix for our GLM.
Hence we can use the asymptotic result
 n 
g(k) x
X
i xi )
T −1
βb ∼ N β, ϕ W i .
i=1

Exercise 16 - IWLS algorithm

Show that for a normal linear model, the IWLS algorithm converges to the maximum likelihood
estimate in one iteration, whatever starting value we use for β.
88 CHAPTER 6. MAXIMUM LIKELIHOOD ESTIMATION FOR GLMS

Solution:
Let β
e be our starting value for β. We then solve the approximate score equation

n g
Wi

zei − xTi β xi =
X
0
i=1 ϕ

for β. Note that this is equivalent to solving

n  
g ze − xT β x =
X
W i i i i 0
i=1

where

µ
fi = E(Yi | β)
e
 
zei = g ′ (µ
fi ) yi − µfi + xTi β
e and
wi
W
g
i =
[g ′ (µ
fi )]2 V (µ
fi )

For a normal linear model, we have

g(µi ) = µ,
g ′ (µi ) = 1.
V (µi ) = 1.
wi = 1.

Hence
wi
W
g
i = = 1.
[g (µ
′ fi )]2 V (µ
fi )

Furthermore
   
zei = g ′ (µ
fi ) yi − µ
fi + xTi β
e = e + xT β
yi − xTi β i
e = yi .

The approximate score equation now becomes

n   n  
g ze − xT β x = yi − xTi β xi =
X X
Wi i i i 0,
i=1 i=1

which are the usual normal equations for a normal linear model. Note that these equations no
longer depend on the starting value β
e for β. We can write the equations in matrix notation as

 T
y − Xβ X = 0

which can be solved in closed form as

β
b = (X T X)−1 X T y.
6.4. ESTIMATING THE DISPERSION PARAMETER 89

6.3.3 Convergence in practice

In practice we stop the IWLS algorithm after a finite number of iterations.
• The convergence criterion is based on relative changes in the deviance.
• By maximizing the log likelihood over β for fixed ϕ we are minimizing the deviance.
• When relative changes to the deviance are sufficiently small, then we have converged.
The glm function in R derives starting values from the data.
(0)
We only need an initial estimate µe i for µi , for example
(0)
• normal, gamma, inverse Gaussian: µe i = yi
(0)
• Poisson: µe i = yi + 0.1
(0)
• scaled binomial: µe i = (y(m
i mi +0.5)
i +1)

and then our initial working observations can be derived from

 
(0) (0) (0) (0)
zei = g ′ (µe i ) yi − µe i + g(µe i ).

This only works for certain link functions.

Exercise 17 - Computer Practical 3

Solutions to Computer Practical 3 are provided on moodle.

6.4 Estimating the dispersion parameter

6.4.1 Overview
• The maximum likelihood estimator β b does not depend on ϕ.
• We estimate β first and then estimate ϕ in a second step.
• The IWLS algorithm gives us an estimator for ϕ based on the linear model.
• This estimator reduces to an intuitively clear form based on the sum of squares of the
residuals.

6.4.2 Derivation
The last iteration of the IWLS algorithm is based on a linear approximation

ϕ
 
Z
c
i ∼ N xTi β,
Wi
d

where  
zbi = g ′ (µb i ) yi − µb i + g(µb i )
with
wi
g(µb i ) = xTi β
b and W
d
i = h i2 .
V (µb i ) g ′ (µb i )
90 CHAPTER 6. MAXIMUM LIKELIHOOD ESTIMATION FOR GLMS

Under the normal weighted linear model approximation the estimator of ϕ is given as
1 X n  2
ϕb = Wi i
d z
b − x Tb
i β .
n − p i=1

Note the denominator is n − p as we lose p degrees of freedom from estimating the p parameters.
Substituting expressions for W
d and zb :
i i

1 X n  2
ϕb = Tb
Wi zbi − xi β
d
n − p i=1
1 X n 2
wi
  
= g ′
(µ i ) y i − µ i +x Tb
i β − xTb
i β
n − p i=1 V (µb i ) g ′ (µb i ) 2
h i b b

1 X n
wi h i2  2
= i2 g (µ
′
bi yi − µb i
n − p i=1 V (µb i ) g ′ (µb i )
h

1 X n
wi  2
= yi − µb i
n − p i=1 V (µb i )

The Pearson residual for observation i is

s
(p) wi
ri = (yi − µb i ).
V (µb i )

We can thus write the estimator of ϕ as

1 X n  2
(p)
ϕb = ri .
n − p i=1

In R we can get the Pearson residuals from a fitted GLM with the command
residuals([Link], type="pearson")
We will see later that there are many possible types of residual for a GLM.
The Pearson estimator for ϕ is the one used by the summary() function in R.

6.4.3 Other estimators for the dispersion parameter

Other possible estimators discussed by Dunn & Smyth in Section 6.8 of the recommended
textbook.2
• Modified profile likelihood. Optimal estimator but requires stronger assumptions
about the distribution of Y and usually requires numerical maximization.
• Mean deviance. Not suitable for Poisson or binomial models with small counts (y < 3
or, for binomial, m − y < 3.)
For normal linear models, all three estimators of ϕ are the same.
2
Dunn, P. K. and Smyth, G.K (2018): Generalized linear models with examples in R Vol. 53. New York:
Springer.
6.4. ESTIMATING THE DISPERSION PARAMETER 91

6.4.4 Example: cherry tree data

We illustrate the estimation of the dispersion parameter with the trees example from the
datasets package. (You will recall the dataset from ST221/ST231!)
• We want to predict the volume (V) of wood from the height of the tree (H) and its diameter
(G).
• Suppose the tree is a cylinder, then
 2
V = πH G/2
log(V ) = log(π/4) + 2 log(G) + log(H)

This suggests the following model for µ = E(V ):

log(µ) = β0 + β1 log(G) + β2 log(H)

As V is positive and real-valued we use the gamma EDM:

V ∼ Γ(µ, ϕ).

The R summary of the model is given below.

Call:
glm(formula = Volume ~ log(Height) + log(Girth), family = Gamma(link = "log"),
data = trees)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.69111 0.78784 -8.493 3.11e-09 ***
log(Height) 1.13288 0.20138 5.625 5.04e-06 ***
log(Girth) 1.98041 0.07389 26.802 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Gamma family taken to be 0.006427286)

Null deviance: 8.31720 on 30 degrees of freedom

Residual deviance: 0.18352 on 28 degrees of freedom
AIC: 139.9

Number of Fisher Scoring iterations: 3

Chapter 7 Diagnostics for GLMs

7.1 Residuals for GLMs

For linear models we used the notation
Yi = xTi β + ϵi
which clearly separates the structural part of the model xTi β from the random error part ϵi .
The residual
ri = yi − xTi β b

can be thought of as an estimate of ϵi . In contrast, for GLMs there is no clear separation of

structural and error components. Therefore we need a more general concept of residuals.
There are three basic requirements for residuals:
1. A residual measures the discrepancy between the observed value yi and the fitted value
µb i = µb i (β).
b
2. A residual is zero if yi = µb i .
3. A residual is signed:
• positive if yi > µb i ;
• negative if yi < µb i .
For GLMs there are (at least) four definitions of residuals that satisfy these requirements.

7.1.1 Response residuals

The response residual is the difference between the observed and fitted values
(r)
ri = yi − µb i .
For GLMs we need additional information about prior weight wi and variance function V (µi ) to
interpret the response residuals.

7.1.2 Pearson residuals

The Pearson residual adjusts for the variance function and prior weight.
(p) yi − µb i
ri = q .
V (µb i )/wi
As sample size n → ∞,

 
(p)
E ri −→ 0,
 
(p)
Var ri −→ ϕ(1 − hi )
where hi is the leverage of the ith observation (see later).

92
7.1. RESIDUALS FOR GLMS 93

7.1.3 Deviance residuals

The deviance residual is
q
(d)
ri = sign(yi − µb i ) wi d(yi , µb i )

where d(y, µ) is the unit deviance and

1 if T > 0



sign(T ) = 0 if T = 0
−1 if T < 0.



Recall from Chapter 4 that the total deviance from a fitted model is
n
= wi d(yi , µb i )
X
D(y, µ)
b
i=1

This can be rewritten in terms of the deviance residuals

n  2
(d)
=
X
D(y, µ)
b ri .
i=1

This generalizes the result for (unweighted) linear models that the deviance is the sum of squares
of the residuals.
In general, the deviance residual is a non-linear function of the outcome variable Yi and is not
guaranteed to have expectation zero:
(d)
E(ri ) ̸= 0.

This limits the usefulness of deviance residuals for diagnostic plots.

7.1.4 Working residuals

The working residuals are the residuals from the last iteration of the IWLS algorithm
(w)
ri = zbi − xTi β
b

where zbi is the working observation

zbi = (yi − µb i )g ′ (µb i ) + xTi β

b

The working residual can be simplified to

(w)
ri = (yi − µb i )g ′ (µb i )

Like the response residuals, the working residuals do not account for difference in variance
between observations. You need the working weights Wc to interpret them.
i
94 CHAPTER 7. DIAGNOSTICS FOR GLMS

Exercise 18 - Residuals in a linear model

Show that for an unweighted normal linear model

Yi ∼ N (xTi β, ϕ)

the different types of residuals are all the same, that is

(r) (p) (d) (w)
ri = ri = ri = ri .

Hint: for a unweighted normal linear model we have

V (µ) = 1,
g(µ) = µ,
d(y, µ) = (y − µ)2 .

Solution: We start with the response residuals. As g(µ) = µ, we have µb i = xTi β.

b Then,

(r)
ri = yi − µb i = yi − xTi β.
b

Next consider the Pearson residuals. Recall that V (µb i ) = 1 and wi = 1. Hence

(p) yi − µb i
ri = q = yi − µb i = yi − xTi β.
b
V (µb i )/wi

For the deviance residuals recall that the unit deviance of a normal linear model is given by
d(y, µ) = (y − µ)2 and so
q
(d)
ri = sign(yi − µi ) wi d(yi , µb i )
b = sign(yi − µb i ) |yi − µb i | = yi − µb i = yi − xTi β.
b

Finally, for the working residuals, recall that the working observations are defined as

zbi = (yi − µb i )g ′ (µb i ) + xTi β

b = (yi − µb i ) + µb i = yi ,

where we have used the fact that g ′ (µ) = 1. Then

(w)
ri = zbi − xTi β
b = yi − xTi β.
b

Note that if we consider a weighted linear model, that is wi ̸= 1 for some i, then we have two
sets of residuals
(p) (d) √
ri = ri = wi (yi − xTi β)
b
(r) (w)
ri = ri = yi − xTi β.
b

Note that these correspond to the weighted and unweighted residuals discussed in Chapter 1.
7.2. LEVERAGE AND STANDARDIZED RESIDUALS 95

7.1.5 Residuals and weights in R

The residuals() function will extract the various residuals from a fitted model.
Use the type argument to specify which residuals you want.
• residuals([Link], type="response")
• residuals([Link], type="pearson")
• residuals([Link], type="deviance")
• residuals([Link], type="working")
If you omit the type argument, then you get the deviance residuals by default. Therefore it is
strongly recommend to always specify type.
The weights() function will extract the weights from a fitted model.
Use the type argument to specify which weights you want.
• weights([Link], type="prior")
• weights([Link], type="working")
If you omit the type argument, then you get the prior weights by default.

7.1.6 Which residuals to use?

• Use Pearson residuals primarily for diagnostics.
• Deviance residuals may be more useful than Pearson residuals for QQ plots (see later).
• Pearson and deviance residuals are often the same when the amount of information in one
observation is large. This is not the case for Poisson and binomial data with small counts.
• Response and working residuals exist but are not useful for model criticism. However, we
use them to calculate other quantities.

7.2 Leverage and standardized residuals

7.2.1 Leverage in normal linear models - Revision
Consider the weighted linear model

Yi ∼ N (xTi β, ϕ/wi ).

The residual Ri = Yi − µb i has mean and variance

E(Ri ) = 0,
(1 − hi )ϕ
Var(Ri ) = ,
wi
where hi = H ii is the ith diagonal element of the hat matrix
 −1
H = X XT W X XT W
96 CHAPTER 7. DIAGNOSTICS FOR GLMS

and W = diag(w1 , . . . , wn ). See Chapter 1.

The hat value hi ∈ [0, 1] is a measure of leverage for observation i. Observations with high
leverage force the fitted value from the model to be close to the observed value. If observation i
has maximum leverage hi = 1, then
Var(Ri ) = 0.
Hence ri = E(Ri ) = 0 and µb i = yi . So the observed value and the fitted value are identical.
Example: Anscombe’s Quartet
These are four bivariate data sets with identical basic statistics (means, variances, correlations,
R2 etc.) but totally different shapes when the data are plotted.
In this example from Anscombe’s quartet, the observation on the right has maximum leverage
hi = 1 forcing the fitted line through the observed point.
12
10
Y

8
6

8 10 12 14 16 18

Example: Simple linear regression

Suppose µi = α + βxi . Then
1 (xi − x )2
hi = + Pn .
j=1 (xj − x )
n 2

So the hat value hi measures how far away the predictor variable xi is from the mean x .
Suppose x ∈ {0, 1} is a binary predictor with
• n0 observations with xi = 0 and
• n1 observations with xi = 1.
Then
if xi = 0,
( 1
hi = n0
1
n1
if xi = 1.

You should expect to see observations with high leverage if they have unusual values for the
predictor variables:
• for continuous variables if they have values far from mean;
• for categorical variables if they belong to a sparse category.
7.2. LEVERAGE AND STANDARDIZED RESIDUALS 97

7.2.2 Standardized residuals in normal linear models - Revision

The standardized residual for a normal linear model is the residual divided by an estimate of its
standard deviation. √
(s) wi (yi − µb i )
ri = q
(1 − hi )ϕb
where ϕb is the unbiased estimator of the dispersion parameter ϕ, that is
1 X n
ϕb = wi (yi − µb i )2 .
n − p i=1

The standardized residual then has asymptotic variance 1 as n → ∞.

7.2.3 Leverage for GLMs

The final iteration of the IWLS algorithm is based on the approximation
ϕ
 
T
zbi ∼ N xi β, c
Wi
for working observation zbi and working weight W
c evaluated at β.
i
b

Hence for the working residuals

(w) (w)
E(ri ) ≈ 0, Var(ri ) ≈ (1 − h
b )ϕ/W
i
c,
i

where h
b =H
i
c for the working hat matrix
ii
 −1
H
c = X XT W
dX XT W
d

and W
d = diag(W c ) is a diagonal matrix of working weights.
c ,...,W
1 n

If we standardize the working residuals, then we get

q
d(zb − xT β)
W b
(sw) i i
ri = q i
.
(1 − h
c)ϕ
i
b

Substituting the definitions of zbi and W d (after some cancellations) gives

i
√ (p)
(sw) wi (yi − µb i ) r (sp)
ri = q = q i = ri .
(1 − hi )V (µb i )ϕ
c b (1 − hi )ϕ
c b

(sp)
This is also the standardized Person residual ri .
We can also standardize the deviance residual using the working hat matrix
q
(sd) wi d(yi , µb i )
ri =
sign(yi − µb i ) q .
(1 − hc)ϕ
i
b

The rstandard() function extracts standardized residuals from a fitted model

• rstandard([Link], type="deviance")
• rstandard([Link], type="pearson")
98 CHAPTER 7. DIAGNOSTICS FOR GLMS

7.3 Influence
Model fitting should not depend on a handful of observations. We therefore need to consider
what happens if we drop an observation from the data set and refit the model. If the estimates
or predicted values change substantially, we say the observation is influential.

7.3.1 Influence in normal linear models

Consider the weighted normal linear model:

ϕ
 
Yi ∼ N µi , ,
wi
µi = xi β T .

We will use a superscript (i) to denote a quantity estimated after deleting observation i. Then,
the effect on the estimated parameters can be measured as
ri
 
(i)
b −β
β b = (XW X)−1 xi wi
1 − hi
b (Proof omitted).
where W = diag(w1 , . . . , wn ) and ri = yi − xTi β
This quantity is known as dfbeta. It can be calculated without refitting the model.
The effect on fitted values is given by
ri ri
   
b (i)
 
xTj (XW X)−1 xi wi
(i)
µb j − µb j = xTj β
b −β = = H ji .
1 − hi 1 − hi

Note that dropping observation i perturbs the fitted values for j ̸= i.

The Cook’s distance for a linear model summarizes the changes to the fitted values across all
observations.
1 X n  2
(i)
DiCook
= w j µj − µj
b b
p ϕb j=1

where ϕb is the unbiased estimator of ϕ.

(i)
We may simplify the expression for the Cook’s distance as follows. Substituting µb j − µb j gives

1 X n
ri 2 ri2 n
 
DiCook = 2
= wj H 2ji
X
wj H ji
pϕb j=1 1 − hi (1 − hi )2 pϕb j=1
ri2   ri2
= H T
W H = (W H)ii
(1 − hi )2 pϕb ii (1 − hi )2 pϕb
ri2 1 wi ri2 hi
= h i w i = b − h ) (1 − hi )
(1 − hi )2 pϕb p ϕ(1 i

The above uses the result H T W H = W H (proof omitted).

7.3. INFLUENCE 99

Recall the standardized residual for a linear model

√
(s) w i ri
ri = q .
b −h)
ϕ(1 i

and so the Cook’s distance is given by

1 (s)
2
hi

DiCook = r .
p i (1 − hi )

7.3.2 Cook’s distance for GLMs

We take the values from the last iteration of the IWLS algorithm
1 (sw)
2
ĥi

DiCook = r
p i (1 − ĥi )
(sw)
where ri is the standardized working residual.
(sw) (sp)
Recall that ri = ri , the standardized Pearson residual and so
1 (sp)
2
h
 b
i
DiCook = r .
p i (1 − h b )
i

Thus the Cook’s distance of observation i is large if

• the observation has high leverage (h
b close to 1) and/or
i
(sp)
• the absolute value of the standardized Pearson residual |ri | is large.
As sample size n → ∞
n
X 
E DiCook −→ 1.
i=1
(Proof omitted.)
Therefore we expect the influence of individual observations to diminish at the rate O(n−1 ).
There is no commonly accepted threshold for when an observation should be considered
influential. Instead look at the distribution of DiCook across observations and look for values
that are much larger than the others.
Question: If we do identify influential observations, should we delete them?
Answer: Not necessarily. We should only delete data if we have a very strong justification for
doing so! Investigate first, there may be good reasons why an observation is influential.

7.3.3 Summary
1. Model fitting should not depend on a handful of observations.
2. Look for observations with high influence:
• large absolute residuals and/or
• high leverage.
3. Investigate observations with much larger influence than the rest (if any).
4. Sensitivity analysis to check stability of results if you do delete an observation.
100 CHAPTER 7. DIAGNOSTICS FOR GLMS

7.4 Diagnostic plots in R

The plot() method for fitted GLM objects in R produces a sequence of diagnostic plots, using
• Pearson residuals,
• standardized Pearson residuals,
• standardized deviance residuals,
• fitted values,
• hat values,
• Cook’s distance.
Not all of these plots are appropriate all the time. You can select a single plot using the which
argument. In the following we illustrate the plot function using a data set that was produced
from a normal linear model and to which we fit a normal linear model. Thus we know that the
model is appropriate for the data.

7.4.1 Residuals vs fitted values

With which=1 you get a scatter plot of:
(p)
• Pearson residuals on the y-axis (ri );
• fitted values on the scale of the linear predictor on the x-axis (ηbi = xTi β).
b

A red line shows a smooth curve fitted to the residuals. If the structural part of the model is
correct (link function and predictor variables), the smooth curve should remain close to the
value zero on the y-axis.

Residuals vs Fitted
3

106 174
2
Pearson Residuals

1
0
−1
−2
−3

132

0 2 4 6 8

Predicted values
glm(y ~ x)

7.4.2 QQ plot
With which=2 you get a quantile-quantile (QQ) plot with
• sorted absolute values of the standardized deviance residuals on the y-axis;
• theoretical quantiles of the standard half-normal distribution on the x-axis.
Ideally, the points should line up in a straight line through the origin and with slope 1. The
(sd)
observations with the 3 largest values of ri are labelled with their row numbers.
7.4. DIAGNOSTIC PLOTS IN R 101

Q−Q Residuals

0.0 0.5 1.0 1.5 2.0 2.5 3.0

132
174
106

Std. Deviance resid.

0.0 0.5 1.0 1.5 2.0 2.5 3.0

Theoretical Quantiles
glm(y ~ x)

The absolute value of the standardized deviance residuals has a half-normal distribution when
the saddlepoint approximation applies:
• exact for normal and inverse Gaussian,
• approximate for
– Gamma for shape > 3 (ϕ < 1/3),
– Poisson for y ≥ 3,
– Binomial for 3 ≤ y ≤ m − 3.

7.4.3 Scale Location plot

With which=3 you getqa scatterplot very similar to the plot of residuals vs fitted values (which=1),
q
(sp) (p) (sp)
but on the y-axis is |ri | instead of ri . A red line shows a smooth curve fitted to |ri |.
If the variance function is correct, the red line should show no upward or downward trend.
Observations with the largest standardized Pearson residuals (in absolute value) are labelled
with the row numbers.

Scale−Location
132
174
106
1.5
Std. Pearson resid.

1.0
0.5
0.0

0 2 4 6 8

Predicted values
glm(y ~ x)

7.4.4 Cook’s distance

With which=4 you get a plot of Cook’s distance against observation [Link] plot is an
index plot with vertical lines rising from zero on the y-axis to the Cook’s distance value. The
102 CHAPTER 7. DIAGNOSTICS FOR GLMS

observations with the 3 largest values of Cook’s distance are labelled with the row number. This
plot is normally omitted if you do not specify the which argument.

Cook's distance

0.08
132

0.06
Cook's distance
47

0.04
106

0.02
0.00

0 50 100 150 200

Obs. number
glm(y ~ x)

7.4.5 Standardized residuals versus leverage

With which=5 you get a scatter plot of:
(sp)
• standardized Pearson residuals on the y-axis (ri );
• hat values on the x-axis (h
b ).
i

A red line shows a smooth curve fitted to the residuals. Note that Cook’s distance can be
(sp)
calculated from ri and h
b .
i

• The observations with the 3 largest values of DiCook are labelled with their row numbers.
• Dashed grey lines show contours of the Cook’s statistic for DCook = 0.5 and DCook = 1 to
highlight highly influential observations.
• Contours may not be shown if they lie outside the boundaries of the plot, but a legend is
still produced.

Residuals vs Leverage
3

106
47
2
Std. Pearson resid.

1
0
−1
−2
−3

132
Cook's distance

0.000 0.005 0.010 0.015 0.020

Leverage
glm(y ~ x)

7.4.6 When diagnostics break down

For discrete outcomes (Poisson and binomial) with small counts, the diagnostic plots do not
work well.
7.4. DIAGNOSTIC PLOTS IN R 103

For binary data, only the first diagnostic plot is of interest. Do not try to interpret the other
plots (except leverage plot which=4).
Binned residual plots may be more visually appealing for discrete data.

Exercise 19 - Diagnostic plots

Perform a residual analysis for the examples encountered in the computer practicals.
Solution: R code to perform the residual analysis is available on moodle. The diagnostics plot
were discussed in the lectures.
Chapter 8 Model choice and Hypothesis tests

8.1 Confidence intervals based on asymptotic normality

8.1.1 Asymptotic normality for GLMs
Under regularity conditions, as n → ∞, β
b has a multivariate normal distribution

β
b ∼ Np (β, I(β)−1 ),
where I(β) is the Fisher (or expected) information matrix. Here

I(β)−1 = ϕ (X T W
dX)−1 .

where Wd is a diagonal matrix of working weights weights from the final iteration of the IWLS
algorithm.

8.1.2 Wald confidence intervals

Wald confidence intervals are based on the asymptotic normal distribution.
If ϕ is fixed:
The 95% Wald confidence interval for βj is given by

βbj ± 1.96 × SE(βbj )

where h i1/2
SE(βbj ) = ϕ (X T W
dX)−1
jj

and 1.96 is the 97.5 percentile of the standard normal.

If ϕ is free:
h i1/2
SE(βbj ) = ϕb (X T W
dX)−1
jj

where ϕb is estimated from the Pearson residuals

1 n h
(P ) 2
i
=
X
ϕb ri
(n − p) i=1

As in the normal linear model, we use the t-distribution to account for additional uncertainty
due to estimation of ϕ, for example confidence limits are

βbj ± tn−p (0.975) × SE(βbj )

where td (π) is the π-quantile of the t-distribution on d degrees of freedom.

104
8.2. CONFIDENCE INTERVALS BASED ON PROFILE LIKELIHOOD 105

8.1.3 Quantiles in R
All distributions in R have a function with prefix q that gives the quantiles of the distribution,
for example
qnorm(0.975) #Standard normal distribution
[1] 1.959964
qt(0.975, df=20) # t-distribution on 20 df
[1] 2.085963
qt(0.975, df=2) # t-distribution on 2 df
[1] 4.302653

(note that the uncertainty becomes much larger when you have very few d.f.)
The general rule for a confidence interval with coverage 100 × (1 − α) is to call these quantile
functions with argument 1 − α/2.

8.2 Confidence intervals based on profile likelihood

8.2.1 Deviance and log-likelihood
The deviance measures the discrepancy between the data y and the expected value under the
model µ(β). We can consider the deviance as a function of the parameters β
D(β) = D(y, µ(β)).

Recall the dispersion model form of an EDM

wi d(yi , µi )
!
p(yi | µi , ϕ/wi ) = a (yi , ϕ/wi ) exp −
∗
2ϕ
Considering µi as a function of β, and holding ϕ fixed, we can write the log likelihood as
log L(β) = log p(yi | µi (β), ϕ/wi )
X

i
1 X
= wi d(yi , µi ) + log a∗ (yi , ϕ/wi )
X
−
2ϕ i i
1
= − D(y, µ(β)) + log a∗ (yi , ϕ/wi )
X
2ϕ i

Suppose β A and β B are two competing possible values for β. We can compare the relative
support for the two values by comparing the log likelihood
1
log L(β A ) − log L(β B ) = − [D(y, µ(β A )) − D(y, µ(β B ))]
2ϕ
Hence differences in deviance are scaled differences in log likelihood
D(y, µ(β A )) − D(y, µ(β B )) = −2ϕ [log L(β A ) − log L(β B )]
L(β A )
!
= −2ϕ log
L(β B )
106 CHAPTER 8. MODEL CHOICE AND HYPOTHESIS TESTS

The profile deviance is a function of a single element of β

Diprof (b) = inf {D(y, µ(β)) : βi = b}

We focus interest on the ith element of β, and eliminate nuisance parameters βj for j ̸= i by
minimizing the deviance over all values of β that have βi = b.
There is a different profile deviance function for each of the p elements of β. Hence the subscript
i in Diprof (b).
We can create a set of values for parameter βi most consistent with the data by choosing a set

n o
b : Diprof (b) ≤ K

for some cut-off value K. Note that such a confidence interval is not restricted to be symmetric!
K is chosen based on the repeated sampling principle:
• We want to use statistical procedures that give the right answer most of the time.
• We want to choose procedures that have good properties under hypothetical repetition,
assuming the model is correct.
• For confidence intervals the key property is coverage.
– Coverage is the proportion of times that a confidence interval generated using the
same procedure contains the true value under repeated sampling.
– By convention we use 95% coverage.

8.2.2 Calibrating the profile deviance

If ϕ is fixed:
For n → ∞
Diprof (βi ) − Diprof (βbi ) ∼ ϕχ21 .
By convention, for an EDM with fixed dispersion we set ϕ = 1. Hence we choose our cut-off as
the 95th percentile of χ21 , and our confidence interval for βi is
n o
b : Diprof (b) − Diprof (βbi ) ≤ 3.84

If ϕ is free:
We plug in the estimate ϕb and adjust for the additional uncertainty from estimating ϕ. Then,
for n → ∞,
Diprof (βi ) − Diprof (βbi )
∼ F1,n−p .
ϕb

Hence we choose our cut-off as the 95th percentile of the F1,n−p distribution.
8.2. CONFIDENCE INTERVALS BASED ON PROFILE LIKELIHOOD 107

In R:
qchisq(0.95, df=1)
[1] 3.841459
qf(0.95, df1=1, df2=20)
[1] 4.351244
qf(0.95, df1=1, df2=2)
[1] 18.51282

For a confidence interval with coverage 100 × (1 − α), call these quantile functions with argument
1 − α.
Recall the Wald statistics used to derive the CI based on asymptotic normality. The Wald
statistics are based on a quadratic approximation (for n → ∞) to the likelihood.

log L(β) ≈ b T I(β)(β

−(β − β) b − β)/2
b

This is exact for the normal linear model.

Profile deviance confidence intervals take on the simplified Wald form when the quadratic
approximation holds.

8.2.3 Confidence intervals for transformed parameters

Often we do not present the parameters β directly, but work with transformations of these
parameters:
• In a logistic model, parameter βx corresponding to predictor x is the log odds ratio.
• We present results in terms of the odds ratio exp(βx ).
We calculate the confidence limits for β first, then transform them to give the limits of the
transformed parameter.
• If βbx ± 1.96 × SE(βbx ) are the limits for βx then the limits for exp(βx ) are
h i
exp(βbx − 1.96 × SE(βbx )), exp(βbx + 1.96 × SE(βbx )) .

8.2.4 Example
Below the computation of confidence intervals in R is demonstrated. The relevant R code is
available on Moodle.
The Wald confidence intervals can be computed from the quantities provided by the model
summary output.
To compute profile confidence intervals we use the command confint(). This is a numerically
expensive operation, so R will give a warning message.
library(faraway)
data(wcgs)
wcgs$height = wcgs$height * 0.0254
wcgs$weight = wcgs$weight * 0.453592
108 CHAPTER 8. MODEL CHOICE AND HYPOTHESIS TESTS

wcgs$bmi = wcgs$weight/wcgs$heightˆ2

wcgs$bmicat = cut(wcgs$bmi, breaks=c(18.5, 25, 30, Inf), right=FALSE,

labels=c("normal", "overweight", "obese"))

[Link] <- glm(chd ~ bmicat + age + cigs + dbp + sdp + dibep,

family=binomial, data=wcgs)
faraway::sumary([Link])

Estimate Std. Error z value Pr(>|z|)

(Intercept) -8.3087702 0.7649847 -10.8614 < 2.2e-16
bmicatoverweight 0.1269456 0.1411712 0.8992 0.368529
bmicatobese 0.4673957 0.3438657 1.3592 0.174071
age 0.0628743 0.0118178 5.3203 1.036e-07
cigs 0.0224560 0.0041512 5.4095 6.322e-08
dbp 0.0029435 0.0103221 0.2852 0.775521
sdp 0.0196150 0.0062805 3.1232 0.001789
dibepB -0.6845462 0.1435946 -4.7672 1.868e-06

n = 3131 p = 8
Deviance = 1639.18499 Null Deviance = 1777.31898 (Difference = 138.13399)
To compute the Wald confidence intervals we will need to extract the appropriate summary
statistics, as below.
[Link] <- coef([Link])
SE <- sqrt(diag(vcov([Link])))
CI <- cbind(estimate=[Link],
lower=[Link] - 1.96 * SE,
upper=[Link] + 1.96 * SE)
CI
estimate lower upper
(Intercept) -8.308770209 -9.808140159 -6.80940026
bmicatoverweight 0.126945642 -0.149749854 0.40364114
bmicatobese 0.467395714 -0.206581090 1.14137252
age 0.062874304 0.039711377 0.08603723
cigs 0.022455987 0.014319544 0.03059243
dbp 0.002943465 -0.017287862 0.02317479
sdp 0.019614957 0.007305223 0.03192469
dibepB -0.684546170 -0.965991606 -0.40310073

The above gives the Wald confidence intervals on the scale of the parameter. We transform the
intervals to obtain the intervals for the relevant odds rations.
round(exp(CI[-1,]), 2)
estimate lower upper
bmicatoverweight 1.14 0.86 1.50
bmicatobese 1.60 0.81 3.13
8.3. HYPOTHESIS TESTS FOR GLMS 109

age 1.06 1.04 1.09

cigs 1.02 1.01 1.03
dbp 1.00 0.98 1.02
sdp 1.02 1.01 1.03
dibepB 0.50 0.38 0.67

The confint() function computes confidence intervals from the profile deviance.
CI2 <- confint([Link])

Waiting for profiling to be done...

round(exp(CI2[-1,]), 2)

2.5 % 97.5 %
bmicatoverweight 0.86 1.50
bmicatobese 0.78 3.03
age 1.04 1.09
cigs 1.01 1.03
dbp 0.98 1.02
sdp 1.01 1.03
dibepB 0.38 0.67
The profile likelihood based confidence intervals are similar but not all identical to the Wald
confidence intervals.
Please try for yourself the code provided in Moodle for the lime data where you also need to
estimate the dispersion parameter ϕ.

8.3 Hypothesis tests for GLMs

8.3.1 Overview
1. Wald tests for individual parameters.
• Available from the summary() function.
2. Analysis of deviance for nested models.
• Use the anova() function.
3. Goodness-of-fit.
• Sometimes we can check the overall goodness of fit of a model by comparing the
residual deviance with the number of degrees of freedom. However, circumstances
in which we can do this are limited and this topic requires understanding of the
saddlepoint approximation.

8.3.2 Wald tests

Asymptotically as sample size n → ∞
  −1 
Td
β
b ∼ Np β, ϕ X W X
110 CHAPTER 8. MODEL CHOICE AND HYPOTHESIS TESTS

where X is the design matrix, and W

d is a diagonal matrix of working weights from the last
iteration of the IWLS algorithm.
We consider two cases:
• ϕ fixed (Poisson, binomial, . . . )
• ϕ free (normal, gamma, inverse Gaussian, . . . )

[Link] Wald test for fixed ϕ

Suppose we are interested in βi , the ith element of β.
 −1
Var(βbi ) = ϕ XT W
dX
ii

and so the standard error is given by

  −1 1/2
SE(βbi ) = ϕ XT W
dX
ii

To test H0 : βi = 0 against H1 : βi ̸= 0 we use the Z-statistic

βbi
Zi = .
SE(βbi )

Under H0 : βi = 0,
Zi ∼ N (0, 1)
The p-value is the probability under hypothetical repeated sampling that Zi is more extreme
than the observed value ziobs .
Since we have a two-sided alternative hypothesis, we take the absolute value.
  h  i
P |Zi | ≥ ziobs = 2 1 − Φ ziobs

where Φ is the cumulative distribution function (cdf) of a standard normal distribution.

[Link] Wald test for free ϕ

As before, we have
 −1
Var(βbi ) = ϕ XT W
dX
ii

but now ϕ is unknown. So we plug in our estimate of ϕ from the Pearson residuals to get
standard error   −1 1/2
SE(βi ) =
b b Td
ϕ X WX .
ii

To test H0 : βi = 0 against H1 : βi ̸= 0 we use the T -statistic:

βbi
Ti = .
SE(βbi )
8.3. HYPOTHESIS TESTS FOR GLMS 111

Under H0 : βi = 0, we have Ti ∼ tn−p , that is Ti has a t-distribution on n − p degrees of

freedom.
The p-value is the probability under hypothetical repeated sampling that the Ti is more extreme
than the observed value tobs
i .

Since we have a two-sided alternative hypothesis, we take the absolute value.

  h  i
P |Ti | ≥ tobs
i = 2 1 − Ft tobs
i , n−p

where Ft (·, d) is the cumulative distribution function (cdf) of a td -distribution.

8.3.3 Analysis of deviance

Consider two GLMs A and B such that A is nested in B. This means that
• A and B have the same EDM,
• A and B have the same link function,
• A and B have a common parameter space:
– β ∈ Rp ,
– if βi is free in Model A, then it is free in Model B,
– ϕ is either free in both models or fixed to the same value in both models.
Consider two GLM models with the same EDM and link function and with their linear predictor
being specified by the following R formulae:
• Model A: y ~ age,
• Model B: y ~ age + education
We can consider both models to have the same set of parameters:
• β0 (intercept)
• βage (effect of age)
• βed (effect of years of education)
Hence β = (β0 , βage , βed )T ∈ R3 for both models. But βed = 0 for Model A.
Let
• pA be the number of free parameters in model A,
• b be the maximum likelihood estimate from model A,
β A
• b A = µ(β
µ b ),
A
• ϕA be the estimate of ϕ from model A (based on the Pearson residuals),
b

with pB , β
b ,µ
B b B , ϕB defined similarly.
b

[Link] Analysis of deviance for fixed ϕ

As n → ∞, under model A

b A ) − D(Y , µ
D(Y , µ b B)
∆ = ∼ χ2pB −pA .
ϕ
112 CHAPTER 8. MODEL CHOICE AND HYPOTHESIS TESTS

A large difference in deviance is evidence against model A. Hence the p-value is

P(∆ ≥ ∆obs ) = 1 − Fχ2 (∆obs , pB − pA )

and Fχ2 (·, d) is the cumulative distribution function of the χ2d distribution.

[Link] Analysis of deviance for free ϕ

As n → ∞, under model A
b A ) − D(Y , µ
D(Y , µ b B)
∆ = ∼ (pB − pA ) FpB −pA ,n−pB
ϕbB

Note we use the estimate ϕbB from the larger model, which is a valid estimator for both models.
A large difference in deviance is evidence against model A. Hence the p-value is

∆obs
!
 
P ∆≥∆ obs
= 1 − FF , p B − pA , n − pB
pB − pA

and FF (·, d1 , d2 ) is the cumulative distribution function of the Fd1 ,d2 distribution.
8.3. HYPOTHESIS TESTS FOR GLMS 113

[Link] Analysis of deviance in R

We use the anova() function in R to compare two models, for example
anova(modelA, modelB, test="Chisq").
Note
• the first argument is the smaller model,
• the second argument is the larger model,
• in older versions of R we need to specify the test argument:
– for ϕ fixed: test="Chisq" or test="LRT" (likelihood ratio test)
– for ϕ free: test="F"
The anova() function will NOT give an error or warning if the two models are not nested.
You must check this yourself!
Both models must use the same observations. This might be violated if the there are missing
values. The glm() function will drop observations with missing values and so the larger model
may be fitted to fewer observations. Thus make sure to remove missing values before fitting the
two models.
Consider following model:
[Link] <- glm(chd ~ bmicat + age,
family=binomial, data=wcgs)
sumary([Link])

Estimate Std. Error z value Pr(>|z|)

(Intercept) -6.038707 0.553979 -10.9006 < 2.2e-16
bmicatoverweight 0.215450 0.135127 1.5944 0.11084
bmicatobese 0.760373 0.329026 2.3110 0.02083
age 0.074285 0.011334 6.5543 5.59e-11

n = 3131 p = 4
Deviance = 1728.09442 Null Deviance = 1777.31898 (Difference = 49.22456)
Let’s consider what happens if we remove bmicat from the model.
glm.wcgs2 <- glm(chd ~ age,
family=binomial, data=wcgs)
sumary(glm.wcgs2)

Estimate Std. Error z value Pr(>|z|)

(Intercept) -5.939516 0.549318 -10.8125 < 2.2e-16
age 0.074423 0.011302 6.5847 4.558e-11

n = 3154 p = 2
Deviance = 1738.35613 Null Deviance = 1781.24374 (Difference = 42.88761)
Note the change in the sample size n that is reported at the bottom of the summary. The
variable bmicat has missing data as all participants with BMI less than 18.5 were classified as
114 CHAPTER 8. MODEL CHOICE AND HYPOTHESIS TESTS

NA.
We remove these participants from the data and then refit the model.
# Fitting the model on the reduced dataset
glm.wcgs2 <- glm(chd ~ age,
family=binomial, data=subset(wcgs, bmi>=18.5))
sumary(glm.wcgs2)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -5.938165 0.549709 -10.8024 < 2.2e-16
age 0.074565 0.011311 6.5923 4.332e-11

n = 3131 p = 2
Deviance = 1734.33688 Null Deviance = 1777.31898 (Difference = 42.98210)

We can now perform an Analysis of Deviance comparing model [Link] against [Link]:
anova(glm.wcgs2, [Link])

Analysis of Deviance Table

Model 1: chd ~ age

Model 2: chd ~ bmicat + age
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 3129 1734.3
2 3127 1728.1 2 6.2425 0.0441 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
There is evidence to suggest that, controlling for age, BMI as measured in the variable bmicat,
is associated with the risk of developing chronic heart disease.
If we use the anova() comment with just one model as input it produced a sequential Analysis
of Deviance which is interpreted in a manner analogous to the sequential ANOVA covered in
ST231. Recall that the result will depend on the order in which the predictor variables were
specified when fitting the model.
anova([Link])

Analysis of Deviance Table

Model: binomial, link: logit

Response: chd

Terms added sequentially (first to last)

Df Deviance Resid. Df Resid. Dev Pr(>Chi)

NULL 3130 1777.3
bmicat 2 6.732 3128 1770.6 0.03453 *
8.4. THE SUMMARY FUNCTION 115

age 1 42.493 3127 1728.1 7.095e-11 ***

---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8.4 The summary function

When the summary() function is applied to a fitted GLM, this gives similar output to when it
is applied to a fitted normal linear model. To compare the output we fit a normal linear model
and a Gaussian GLM to the same data.
The table below gives an overview on the differences between the two summary outputs.

normal linear model generalized linear model

Residual statistics —

Residual standard error Estimate of dispersion parameter

R2 , Radj
2
—

F-statistics (with p-value) Deviance statistics (no p-value)

— AIC

— Fisher’s scoring iterations

Call:
lm(formula = log(DMFT) ~ Sugar + Indus, data = dental)

Residuals:
Min 1Q Median 3Q Max
-1.57067 -0.31668 0.02093 0.33304 1.27374

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.373565 0.199409 1.873 0.0644 .
Sugar 0.018351 0.004076 4.502 2.08e-05 ***
IndusNonInd -0.203221 0.143950 -1.412 0.1616
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.557 on 87 degrees of freedom

Multiple R-squared: 0.3005, Adjusted R-squared: 0.2844
116 CHAPTER 8. MODEL CHOICE AND HYPOTHESIS TESTS

F-statistic: 18.69 on 2 and 87 DF, p-value: 1.768e-07

Call:
glm(formula = log(DMFT) ~ Sugar + Indus, family = gaussian(),
data = dental)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.373565 0.199409 1.873 0.0644 .
Sugar 0.018351 0.004076 4.502 2.08e-05 ***
IndusNonInd -0.203221 0.143950 -1.412 0.1616
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.3102132)

Null deviance: 38.584 on 89 degrees of freedom

Residual deviance: 26.989 on 87 degrees of freedom
AIC: 155.01

Number of Fisher Scoring iterations: 2

We note:
• The GLM summary reports the estimated dispersion parameter as ϕb = 0.3102. The
summary of the linear model gives the estimated residual standard error as σb = 0.557 which
corresponds to an unbiased estimate of the error variance given by σb 2 = (0.557)2 = 0.3102
(modulo rounding errors.)
• The summary of the linear model gives an F -statistic equal to Fobs = 18.69 on 2 and
87 degrees of freedom. Recall that this F -statistic compares the null model M0 against
the fitted model M1 . As M1 has an intercept term (and uses the same EDM and link
function), M0 is nested in M1 . The GLM summary gives the deviance of the null model
(null deviance DM0 = 38.584) and the deviance of the fitted model (residual deviance
DM1 = 26.989). As we discussed in Section [Link],

∆ DM0 − DM1
= ∼ Fp−1,n−p ,
p−1 ϕbB (p − 1)

where in our example n = 90 and p = 3 and so the relevant F -distribution has 2 and 87
degrees of freedom. The observed value of p−1
∆
is given by

38.584 − 26.989
= 18.69
2 × 0.3102

Therefore the information provided on the null and residual deviance for the GLM summary
is equivalent to the information on the F -statistics in the summary of the normal linear
model.
8.5. MODEL CHOICE 117

Furthermore, the Akaike Information Criterion (AIC) introduced in ST231 for the normal
linear model can be extended to GLMs and is routinely reported in the glm summary output.
The AIC is defined as
AIC = −2 log L(β, b ϕb
MLE ) + 2p ,
′

where
• L is the likelihood of the chosen model,
• ϕbMLE is the maximum likelihood estimate of ϕ and not the estimate introduced earlier
based on the Pearson residuals,
• p′ is the number of parameters which includes ϕ if the chosen EDM has a free dispersion
parameter.
We note that
• only differences in AIC are relevant,
• the AIC can be used to compare non-nested models,
• the AIC is not consistent.
An alternative that is consistent is the Bayesian information criterion (BIC) which is
defined as
BIC = −2 log L(β, b ϕ
b
MLE ) + p log(n).
′

Thus, provided log(n) > 2, the BIC imposes a heavier penalty for model complexity and will
therefore tend to choose a simpler model than the AIC.

8.5 Model choice

8.5.1 Creating the design matrix
• We start with a data frame Z containing q potential predictor variables (numeric variables
or factors), but then we choose a design matrix X with p columns. How do we get from
one to the other? We may use any of the following:
– transformations (e.g. xi = log(zi ));
– polynomial terms (zi , zi2 , zi3 , . . .);
– interactions;
– variable selection.
Recall from ST231 some of the guiding principles in choosing the design matrix:
• principle of parsimony,
• bias - variance trade-off.

8.5.2 Approaches to model choice

The approach we take to choosing a model is influenced by its purpose.
• Description:
– reduce a large amount of data to simple summaries;
118 CHAPTER 8. MODEL CHOICE AND HYPOTHESIS TESTS

– maintain consistency of methodology for comparability.

• Prediction:

– make good out-of-sample predictions;

– we need to define the cost of wrong predictions.

• Causal reasoning:

– model represents the “true” data generating mechanism;

– we want to interpret the parameters.

8.5.3 All models are wrong, but some of useful.

Image Source: DavidMCEddy at [Link]. CC BY-SA 3.0

George Box (1919-2013), one of the greatest applied statisticians of the 20th century, wrote:

It has been said that “all models are wrong but some models are useful.” In
other words, any model is at best a useful fiction. There never was, or ever will
be, an exactly normal distribution or an exact linear relationship. Nevertheless,
enormous progress has been made by entertaining such fictions and using them as
approximations.

8.6 Variable selection - Review

8.6.1 Graphical representation of nesting
Note that for GLMs, nesting not only considers the linear predictor, but also assumes that the
same EDM and same link function is used.

• Nesting defines a partial ordering between models, which can be represented graphically
by a directed acyclic graph.
• Mi ⪯ Mj if there is a path from Mi to Mj following the arrows.
• If there is no path either way, then the two models are not nested.
8.6. VARIABLE SELECTION - REVIEW 119

8.6.2 Stepwise variable selection methods

The stepwise variable selection methods introduced in ST221/231 can also be applied to GLMs
using the AIC or BIC is a model selection criterion.
Recall there are three possible search strategies:
• forward: start with minimal model, add one predictor at a time.
• backward: start with maximal model, remove one predictor at a time.
• hybrid: start with an initial model, add or remove a predictor term one at a time.
• step(initial_model, direction, scope)
Notes on automatic variable selection:
• It should not replace sound judgement and a good understanding of the data.
• The regression coefficients and standard errors obtained after variable selection are biased.
• Automatic variable selection does not resolve the problem of multi-collinearity.
More modern methods exist, namely shrinkage methods such as ridge regression and LASSO.
This is normally covered in ST404.

Exercise 20 - Computer Practical 5

Explore variable selection in Sections 1, 2 and 5 of Computer Practical 5.
Solutions: See computer practical solutions on Moodle.
Chapter 9 Overdispersion, prediction and clas-
sification

9.1 Overdispersion for count data

Poisson regression makes the assumption that the conditional mean and variance of the outcome
variable, given the predictor variables, are equal.

Var(Yi | xi ) = E(Yi | xi ) = µi (9.1)

Or, in terms of the variance function and dispersion parameter,

V (µ) = µ and ϕ = 1.

Count data that exhibit the behaviour

Var(Yi | xi ) > E(Yi | xi )

are called overdispersed and cannot be represented with a Poisson model. There are several
alternatives:
• the quasi-Poisson distribution,
• robust "sandwich" standard errors,
• the negative binomial distribution.

9.1.1 Robust variance estimation

If the variance function V (µ) is misspecified, then the estimate β
b is inefficient. (We can always
do better, in terms of asymptotic variance, by correctly specifying the variance function.)
In this case
Var(β)
b ̸= b −1
I(β)
and the standard errors reported by the summary() function in R are not correct.
We can approximate the estimating function in the neighbourhood of β
b with a Taylor series.

∂U (β)
b
U (β) ≈ U (β)
b + (β − β)
b
∂β T

= J(β)(
b β b − β)

≈ I(β)(
b β b − β)

120
9.1. OVERDISPERSION FOR COUNT DATA 121

Hence
Var [U (β)] ≈ I(β)
b Var(β)
b I(β)
b

Next we obtain an empirical estimate of the variance of U (β) by taking a scaled sum of squares
of its individual components:
n
n cT
Ω
X
b ≡ U
cU
i i
(n − p) i=1

where
d (y − µ
W b i )g ′ (µ
b i )xi
i i
U
c
i =
ϕ
and we make the necessary adjustment for degrees of freedom lost from estimating β.
Hence
Var [U (β)] ≈ Ω
b

and so
Var(β)
b ≈ b −1 Ω
I(β) b −1
b I(β)

This is called the “information sandwich” estimate of the variance.

In R, this can be computed by the vcovHC() function from the sandwich package.

9.1.2 The negative binomial distribution

The negative binomial distribution is an alternative distribution for count data.
• It includes the Poisson distribution as a limiting case.
• It is typically used when we have over-dispersion relative to the Poisson distribution:

Var(Y ) > µ.

• It is an alternative to the quasi-Poisson family.

Consider a series of Bernoulli trials Zi ∈ {0, 1} with common mean π = P (Zi = 1).
• F is the number of failures F = i I{Zi =0} .
P

• S is the number of successes S = i I{Zi =1} .

P

• M = F + S is the total number of trials.

The binomial distribution arises when the number of trials is fixed (M = m) and we count
the number of successes S.
!
m
P(S = s) = π s (1 − π)m−s
s

for s ∈ {0, 1, . . . m}.

122 CHAPTER 9. OVERDISPERSION, PREDICTION AND CLASSIFICATION

The negative binomial distribution arises when we fix the target number of successes S = s
and we count the number of failures F until we reach the target.

f +s−1
!
P(F = f ) = π s (1 − π)f
s−1

for f ∈ N0 .
The expectation of the negative binomial distribution is given by

(1 − π) s
µ = E(F ) =
π
We can rewrite the probability mass function in terms of µ as

f +s−1 s s µf
!
P(F = f ) =
s−1 (µ + s)f +s

9.1.3 The negative binomial as a mixture distribution

Suppose Yi has a Poisson distribution, but its mean is also a random variable drawn from a
gamma distribution

Yi ∼ Poisson(λi ),
λi ∼ Γ(µ, 1/s),

where we have parameterized the gamma distribution in mean value form with mean µ and
dispersion 1/s.
Then (skipping the derivation)

Γ(y + s) s s µy
P(Yi = y) = .
Γ(s)y! (µ + s)s+y

This generalizes the previous formula to s ∈ R+ .

9.1.4 The negative binomial as an EDM

If s is fixed, then the negative binomial is an exponential dispersion model with canonical
parameter !
µ
θ = log ,
µ+s
cumulant function h i
b(θ) = −s log 1 − exp(θ)
mean and variance function
s exp(θ) µ2
µ = , V (µ) = µ+
1 − exp(θ) s
9.2. INTRODUCTION TO PREDICTION 123

and dispersion parameter ϕ = 1.

As s → ∞, V (µ) → µ and the negative binomial converges to the Poisson distribution.
Recall that an EDM is uniquely determined by its variance function.
For finite values of s, the negative binomial distribution is overdispersed relative to the Poisson
distribution:
µ2
Var(Y ) = ϕV (µ) = µ + > µ.
s

9.1.5 The negative binomial GLM in R

The MASS package is included in every R distribution. MASS provides the [Link]
family for use with glm().
You must specify the overdispersion parameter which is named theta, for example
family=[Link](theta=2).
The default link is log.
MASS also provides a function [Link] used for fitting negative binomial models when the
overdispersion parameter is unknown.
The [Link]() function alternates between two estimation steps:
• For fixed µ, it maximises the likelihood to estimate θ (= overdispersion parameter s).
• For fixed θ, it uses glm() with family="[Link]" to estimate µ.
Initial estimates of µ are obtained by fitting a Poisson model.

Exercise 21 - Overdispersion
Explore overdispersion in Computer Practical 4.
Solution: Solutions to Practical can be found on moodle.

9.2 Introduction to prediction

9.2.1 Predictive models
Often the aim of a statistical model is to make predictions about the future.
For example, when you shop online, a recommender algorithm will often draw your attention to
products that may interest you. This may be based on
• your own personal purchase history, or
• associations between products frequently brought together.
These algorithms are statistical models and their aim is to maximize sales.
• An optimal model is one that makes the best predictions about what you may purchase.
• There is no attempt to understand why you make your choice.
124 CHAPTER 9. OVERDISPERSION, PREDICTION AND CLASSIFICATION

9.2.2 Sequential prediction

Some prediction problems have a natural time ordering. For example,
• weather forecasting,
• prediction in sports.
In such cases we make a prediction for each new event based on all currently available data.

9.2.3 Data replication

Sometimes we have a genuine replicate of the data that can be used to test the model.
In the Washington bicycle hire data, the goal was to predict the demand for rental bicycles. We
have data for years 2011 and 2012. We can
• fit the model to the 2011 data, then
• test predictions on the 2012 data.

9.2.4 Data splitting

If there is no sequential ordering and no replicate data set, then we can randomly split an
existing data set into two parts:
• fit model to the training set, then
• test predictions on the validation set.
We lose information by holding out the validation set, but we gain the ability to assess predictions
on new data. Random sampling ensures that an observation in the training set and an observation
in the validation set are exchangeable.

9.2.5 Double dipping

“Double dipping” is the informal name given to the poor practice of over-fitting a model by
building and evaluating it on the same data.
We have seen three strategies that avoid double-dipping:
• sequential prediction,
• replication, and
• random data splitting.
Example of double dipping:
We generate an artificial data set with 100 units of observation for each variable. We sample 100
potential explanatory variables X1 , . . . , X100 , each iid N (0, 1) and independent of one another,
using the command mydata <- matrix(rnorm(100 * 100), 100, 100).
Next we sample Y from a normal distribution with variance 1 and mean given by

E(Y ) = 1 + X100
using the command y <- rnorm(100, mean=1+mydata[,100]).
9.2. INTRODUCTION TO PREDICTION 125

To select from the 100 potential explanatory variables, we fit 100 simple linear regressions,
using each explanatory variable in turn and choose the predictor variables that are significantly
associated with the response (at a 5% significance level).
We then fit a normal linear model with the chosen predictor variables which produces the
following results:

Call:
lm(formula = y ~ ., data = [Link])

Residuals:
Min 1Q Median 3Q Max
-2.13577 -0.55479 0.05966 0.57150 1.93034

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.08301 0.08889 12.184 < 2e-16 ***
V6 -0.22977 0.09327 -2.464 0.01557 *
V35 -0.27887 0.09728 -2.867 0.00512 **
V54 0.26847 0.08019 3.348 0.00117 **
V85 0.18997 0.08259 2.300 0.02365 *
V100 0.76339 0.08704 8.770 7.56e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8593 on 94 degrees of freedom

Multiple R-squared: 0.5251, Adjusted R-squared: 0.4999
F-statistic: 20.79 on 5 and 94 DF, p-value: 6.204e-14
Now let’s simulate another data set. We keep the same realisations of the 5 predictor variables
but re-sample Y from a normal distribution with variance 1 and mean given by

E(Y ) = 1 + X100 ,

where X100 is our fifth predictor variable. When we fit the a normal linear model to this new
data set, we obtain the following results:
126 CHAPTER 9. OVERDISPERSION, PREDICTION AND CLASSIFICATION

Call:
lm(formula = new.y ~ ., data = [Link])

Residuals:
Min 1Q Median 3Q Max
-3.2440 -0.6340 -0.0391 0.8171 1.9605

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.88582 0.10332 8.573 1.98e-13 ***
V6 0.24212 0.10842 2.233 0.0279 *
V35 -0.16283 0.11308 -1.440 0.1532
V54 -0.06810 0.09321 -0.731 0.4669
V85 0.01671 0.09601 0.174 0.8622
V100 0.83288 0.10118 8.231 1.04e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9989 on 94 degrees of freedom

Multiple R-squared: 0.4827, Adjusted R-squared: 0.4551
F-statistic: 17.54 on 5 and 94 DF, p-value: 3.084e-12
Although our initial model suggested that all five selected predictor variables were informative
about the response, applying the model to new response data tells a different story.

9.2.6 Cross-validation
If the sample size is too small to be able to divide the data into a training and a test set, we
may use cross-validation.
K-fold cross-validation works as follows:
• Divide the data into K subsets.
• One of the K subsets acts as the test set while the others are used as a training set.
• Repeat such that each of the K subsets are used once as a test set.
• Combine the evaluation results.
Leave-one-out cross-validation:
• One observation acts as the “test set” while the others are used for training.
• Repeat such that each observation acts once as the “test set”.
9.3. SCORING RULES FOR PROBABILISTIC PREDICTIONS 127

9.3 Scoring rules for probabilistic predictions

9.3.1 Assessing probabilistic predictions

Suppose our outcome is binary Yi ∈ {0, 1} and we consider predictions for observations in the
validation set. Let µb i be the predicted value of P(Yi = 1) from our [Link] do we assess the
accuracy of our predictions? We need some way to compare the prediction µb i with the observed
outcome yi .

9.3.2 Example: weather forecasting

Weather forecasters never give a deterministic forecast of future weather. They always use
probabilities. Suppose we aggregate all predictions together that have the same forecast
probability of rain.

Forecast (in %) 10 20 30 40 50 60 70 80 90
Days 10 40 105 192 282 200 116 44 11
Rainy Days 3 6 30 76 120 132 70 38 8
% Rainy Days 30 15 29 40 43 66 60 86 73

To simplify calculations we assume forecasts are given to the nearest 10%.

9.3.3 Calibration plots

A calibration plot shows how good the probabilistic predictions are. Each point represents
a group of days with the same forecast probability of rain. For each group we calculate the
observed frequency of rain and compare it with the probabilistic prediction from the model.
Well calibrated predictions should fall on the diagonal line.
128 CHAPTER 9. OVERDISPERSION, PREDICTION AND CLASSIFICATION

100
Frequency of rain (%)

60
20
0

0 20 40 60 80 100

Probability forecast (%)

9.3.4 Scoring rules

A scoring rule S(y, µ)
b is a function that associates any given combination of observation y
and prediction µ
b with a real-valued score that measures how good the prediction is.

Example: The weather forecast says the probability of rain tomorrow is 1%. But it does rain
tomorrow. The purpose of a scoring rule is to assess the quality of the probabilistic predictions.
A scoring rule penalizes the mismatch between prediction µb i = 0.01 and observation y1 = 1.
In this module we consistently use the notation µi to represent the expectation for response of
the ith unit of observation under the model.
• We need to consider what happens when the model is not true.
• Hence we need a new notation for P(Yi = 1) distinct from µi .
• For binary Yi we use πi to represent the true-but-unknown probability that Yi = 1.
• We call π1 , . . . , πn the out-of-model probabilities.
A scoring rule is
• proper if Eπ [S(Y , µ)] is minimized over µ when µi = πi for i = 1 . . . n.
• strictly proper if no other values of µ attain the minimum.
In this context, the expectation Eπ is the out-of-model expectation with P(Yi = 1) = πi .

[Link] The Brier score

The Brier score is a proper scoring rule given by

N
= (yi − µi )2 .
X
S(y, µ)
i=1
9.3. SCORING RULES FOR PROBABILISTIC PREDICTIONS 129

It is the sum of squares of the response residuals, and is equal to the deviance for a Gaussian
outcome variable.
Example The Brier score in the weather example is equal to 227.88.
Proof that Brier score is a proper scoring rule:
We have
N h i
Eπ [S(y, µ)] = Eπ (Yi −πi + πi − µi )2
X

i=1
N h i
= Eπ (Yi − πi )2 + 2(Yi − πi )(πi − µi ) + (πi − µi )2
X

i=1
N N
= πi (1 − πi ) + (µi − πi )2
X X

i=1 i=1

This is minimized over µ if and only if µi = πi .

[Link] The log score

The log score is
N
= [yi log(µi ) + (1 − yi ) log(1 − µi )] .
X
S(y, µ) −
i=1

The log score is equal to minus the log likelihood for a binary outcome:

S(y, µ) = − log L(µ | y).

Example The log score for the weather example is 648.1.

The log score is also a proper scoring rule. Consider

Eπ [S(Y , µ)] = [πi log(µi ) + (1 − πi ) log(1 − µi )] .

X
−
i=1

Differentiating with respect to µi gives:

∂Eπ (S) πi 1 − πi
= − +
∂µi µi 1 − µi
µi − π i
= .
µi (1 − µi )

Thus the expected score is minimized when

∂Eπ (S)
= 0,
∂µi

that is, when µi = πi .

130 CHAPTER 9. OVERDISPERSION, PREDICTION AND CLASSIFICATION

Exercise 21: Brier vs Log score

Large differences between the Brier and the log score appear when the predictions are very
wrong. Suppose we have two models that predict low values for P(Yi = 1).
• Model A: µb A
i = 0.01,
• Model B: µb B
i = 0.001.

Now suppose that we observe Yi = 1.

1. Calculate the difference in log score between model A and B.
2. Calculate the difference in Brier score between model A and B.
Solution:
1. The difference in log score between model A and model B is

Slog (1, 0.01) − Slog (1, 0.001) = − log(0.01) + log(0.001) = −2.3

2. The difference in Brier score between model A and model B is

SBrier (1, 0.01) − SBrier (1, 0.001) = (1 − 0.01)2 − (1 − 0.001)2 = −0.0179

9.4 Prediction versus Classification

Probabilistic predictions are useful, but sometimes we need to make a binary decision or choose
one of two actions based on our probability model.
• How do we choose an optimal decision rule?
• How do we evaluate the decision rule?
The way that classifications are evaluated is different from the way probabilistic predictions are
evaluated.
For consistency of notation we label the outcome of our binary variable:
• positive if yi = 1;
• negative if yi = 0.
We will also label our predictions based on the fitted model
• positive if µb i ≥ α,
• negative if µb i < α.
How do we choose the threshold α?

9.4.1 Naive Bayes classifier

For a naive Bayes classifier, α = 0.5. Each individual is classified in the group more likely to
include it, according to the model.
9.4. PREDICTION VERSUS CLASSIFICATION 131

However, sometimes we are not able to find individuals with a high probability of being positive.
The figure below comes from a model to predict credit defaults and shows the distribution of
fitted probability of default (µb i ) over the individuals in the data set.
• Very few individuals have high predicted risk of default.
• We can see three clusters: low, medium, and high risk.

^
Histogram of µ
1000 1500 2000
Frequency

500
0

0.0 0.2 0.4 0.6 0.8 1.0

^
µ

9.4.2 Minimizing total misclassification

We could try to choose α to minimize the total number of misclassified individuals.
• This depends on the values of πi = P (Yi = 1) in the population.
• The true values of πi are unknown but we can get an empirical estimate of the proportion
of misclassified individuals using a validation sample.
Problems with this approach:
• For a rare outcome, it is easy to get impressive-sounding values for total misclassification.
• Consequences of the two different kinds of misclassification are not the same.
Examples with asymmetric consequences:
1. Lateral flow-test for COVID-19:
• False positive means you need to self-isolate and seek PCR confirmation.
• False negative means you may be spreading the infection without knowing.
2. Biometric access (e.g. fingerprint, face recognition) to your mobile phone:
• False positive means someone else obtains access to your personal data.
• False negative means you may need to try again, or use PIN code access.

9.4.3 Conditional probabilities for classifiers

Asymmetric consequences of false positive and false negative results mean that a single probability
is not sufficient to characterize the predictive value of the classifier. Instead the properties of
132 CHAPTER 9. OVERDISPERSION, PREDICTION AND CLASSIFICATION

the classifier can be summarized by two conditional probabilities. The statistics we use depend
on what we condition on. We should be extremely skeptical if you hear reports that a classifier
is “90% accurate”.
We will use this simple example of a 2 × 2 table to illustrate the calculations. A table like this
is referred to as a confusion matrix.

predicted status true status total

positive negative
positive 90 30 120
negative 10 70 80
total 100 100 200

Choosing the denominator to be the true status we obtain the following two conditional
probabilities:
• Sensitivity:
P(Predict + | True +) = 90/100 = 0.9.
• Specificity:
P(Predict − | True −) = 70/100 = 0.7.
9.4. PREDICTION VERSUS CLASSIFICATION 133

Choosing the denominator to be the predicted status we obtain another set of conditional
probabilities:
• Positive predictive value (PPV)

P(True + | Predict +) = 90/120 = 0.75.

• Negative predictive value (NPV)

P(True − | Predict −) = 70/80 = 0.875.

Suppose we have a test for a rare disease:

• Sensitivity = Specificity = 99%.
• 0.1% of the population have the disease.
Consider the expected outcomes for 100,000 individuals and calculate the PPV and the NPV.

predicted status true status total

positive negative
positive 99 999 1,098
negative 1 98,901 98,902
total 100 99,900 100,000

Note that the PPV depends on sensitivity, specificity and the proportion of positive individuals
in the population.
A hybrid approach is used in machine learning using precision, which is identical to the PPV,
and recall, which is identical to sensitivity. As an example consider the confusion matrix:

predicted status true status total

positive negative
positive 90 30 120
negative 10 70 80
total 100 100 200

Then,
• Precision (= PPV)

P(True + | Predict +) = 90/120 = 0.75.

• Recall (= Sensitivity)

P(Predict + | True +) = 90/100 = 0.90.

134 CHAPTER 9. OVERDISPERSION, PREDICTION AND CLASSIFICATION

9.4.4 Trade off

Whichever pair of statistics we use, changing the threshold α in our decision rule has the same
effect:
• one statistic will increase;
• the other statistic will decrease.
How do we find a good trade off?

[Link] F1 score
In machine learning, the F1 score sometimes used to evaluate predictive accuracy:
1 1 1 1
 
= +
F1 2 Precision Recall

F1 is the harmonic mean of the precision and recall.

Choosing threshold α to maximize F1 is a purely empirical rule designed to penalize small values
of precision or recall.

[Link] Fβ score
A generalization of the F1 score is

1 1 1
!
β2
= +
Fβ 1 + β2 Precision Recall

• β > 1 gives more weight to the recall, for example F2 .

• β < 1 gives more weight to the precision, for example F0.5 .
9.5. ROC CURVES 135

9.5 ROC curves

9.5.1 Definition and examples

An ROC (Receiver Operating Characteristic) curve shows the relationship between
sensitivity and specificity as we change the threshold parameter α from zero (top right) to one
(bottom left).

1.0
0.8
0.6
Sensitivity
0.4
0.2
0.0

1.0 0.5 0.0

Specificity

Note that the x-axis (specificity) is reversed.

Suppose the outcome is independent of all of our predictor variables. For the plot below a
predictor variable was simulated independent of the outcome and then a model was fitted to
the data.
1.0
0.8
0.6
Sensitivity
0.4
0.2
0.0

1.0 0.5 0.0

Specificity

Such a classifier appears as a diagonal line on an ROC plot. This is because what we predict is
independent of the true value.

Sensitivity = P(Predict + | True +) = P(Predict +).

Specificity = P(Predict − | True −) = P(Predict −).
136 CHAPTER 9. OVERDISPERSION, PREDICTION AND CLASSIFICATION

Hence

Sensitivity + Specificity = P(Predict +) + P(Predict −)

= 1.

Sometimes we can see that one classifier dominates another, if it always has higher specificity
for given sensitivity (and vice versa).
1.0
0.8
0.6
Sensitivity
0.4
0.2

with pay
without pay
0.0

1.0 0.5 0.0

Specificity

The above plot shows the ROC curves for two models fitted to the credit default data. One model
uses all available predictors. For the second model we remove the most important predictor
variable pay from the model. Predictions from this reduced model (in blue) are always worse
than the full model (in red).

9.5.2 The AUC statistic

The AUC statistic (Area Under the Curve) measures the area under an ROC curve.
The AUC varies between 0.5 (poor classifier) to 1 (perfect classifier).
Note: In theory, we can have AUC < 0.5, but then we can immediately improve your classifier
by swapping the two classifications.
The AUC measures the ability to rank individuals in order of probability that they are positive.
Suppose we draw one true positive individual and one true negative individual at random from
the population. We know we have one positive and one negative individual. We want to identify
the positive individual. The AUC of the classifier is the probability that we correctly classify
the pair.
The AUC is a clever way of evaluating a classifier without having to make a choice about
threshold α and hence fix the sensitivity and specificity of the classifier. When we misclassify the
pair then we simultaneously make one false positive and one false negative classification. Thus,
the relative cost of false positive and false negative is irrelevant. However, the pair scenario is
artificial. In practice we want to classify individuals not pairs.
9.5. ROC CURVES 137

Exercise 22 - Classification and Prediction

Explore the evaluation of predictive models in Sections 3 and 4 of Computer Practical 5.
Chapter 10 Further topics

10.1 The saddlepoint approximation

The saddlepoint approximation is an approximate density function for an EDM.
1
!
d(y, µ)
p(y | µ, ϕ) ≈ exp −
2ϕ
q
2πϕV (y)

which resembles the dispersion model form (Chapter 4)

!
d(y, µ)
p(y | µ, ϕ) = a (y, ϕ) exp −
∗
2ϕ
and so is equivalent to the approximation
1
a∗ (y, ϕ) ≈ q .
2πϕV (y)

Note that the expression uses V (y) and not V (µ).

10.1.1 The normal distribution

The normal density in dispersion model form is
1
!
d(y, µ)
p(y | µ, ϕ) = √ exp −
2πϕ 2ϕ
for unit deviance
d(y, µ) = (y − µ)2 .

This is exactly equal to the saddlepoint approximation with V (y) = 1.

10.1.2 The inverse Gaussian distribution

The inverse Gaussian density in dispersion model form is
1
!
d(y, µ)
p(y | µ, ϕ) = √ exp −
2πϕy 3 2ϕ
for unit deviance
(y − µ)2
d(y, µ) = .
µ2 y

This is exactly equal to the saddlepoint approximation with

V (y) = y3.

138
10.1. THE SADDLEPOINT APPROXIMATION 139

10.1.3 The Poisson distribution

The Poisson density in dispersion model form is
y y exp(−y)
!
d(y, µ)
p(y | µ) = exp −
y! 2
for unit deviance ! !
y
d(y, µ) = 2 y log − (y − µ)
µ
The saddlepoint approximation is
1
!
d(y, µ)
p(y | µ) ≈ √ exp −
2πy 2

The saddlepoint approximation is not exact for the Poisson distribution. However, we can
obtain the saddlepoint approximation applying Stirling’s formula:
√
y! ≈ 2πy exp(−y)y y

10.1.4 Accuracy of the saddlepoint approximation

For a general EDM, if pe is the approximate density from the saddlepoint approximation then
e | µ, ϕ)
p(y
= 1 + O(ϕ)
p(y | µ, ϕ)
• The saddlepoint approximation is accurate when ϕ is “small”.
• The error is relative. This allows us to carry out integration using the approximate density
and preserve the relative error.

10.1.5 Distribution of the unit deviance

When the saddlepoint approximation holds, the scaled unit deviance has a χ21 distribution
d(Y, µ)
∼ χ21
ϕ
or, equivalently
1 1
!
d(Y, µ) ∼ Γ ,
2 2ϕ
using the result that
1 1
 
χ21 ≡ Γ , .
2 2
Proof:
The mgf of the target Γ( 12 , 2ϕ
1
) distribution is

M (t) = (1 − 2ϕt)−1/2
140 CHAPTER 10. FURTHER TOPICS

We need to prove that this is the mgf of the unit deviance.

Under the saddlepoint approximation:

   Z
M (t) = E exp td = e | µ, ϕ) exp(td)dy
p(y

1
!
Z
d  
= exp − exp td dy
2ϕ
q
2πϕV (y)

1
!
Z
d
= exp − ∗ dy
2ϕ
q
2πϕV (y)

where
ϕ
ϕ∗ = .
1 − 2tϕ

Hence

1
s !
ϕ∗ Z d
M (t) = exp − ∗ dy
2ϕ
q
ϕ 2πϕ∗ V (y)
s
ϕ∗ Z
= e | µ, ϕ∗ )dy
p(y
ϕ
s
ϕ∗
=
ϕ

1
s
=
(1 − 2tϕ)

This completes the proof.

10.1.6 When can we use the saddlepoint approximation?

See Sections 5.4.4–5.4.5 in the textbook by Dunn and Smyth.1

• Normal, inverse Gaussian: always.

• Poisson: y ≥ 3.
• Binomial: y ≥ 3 and m − y ≥ 3.
• Scaled binomial: my ≥ 3 and m(1 − y) ≥ 3.
• Gamma: ϕ ≤ 1/3 (shape ≥ 3).
1
Dunn, P. K. and Smyth, G.K (2018): Generalized linear models with examples in R Vol. 53. New York:
Springer.
10.2. APPLICATIONS OF THE SADDLEPOINT APPROXIMATION 141

10.2 Applications of the saddlepoint approximation

10.2.1 Distribution of the deviance residuals
Recall that
d(Y, µ)
∼ χ21 .
ϕ
This result holds exactly for the normal and inverse Gaussian distributions, and approximately
for
• Poisson and Binomial distributions with large counts (y ≥ 3, m − y ≥ 3),
• Gamma distribution with large shape parameter (shape [= ϕ1 ] ≥ 3).
When applying the saddlepoint approximation to the unit deviance after fitting a model, we
must adjust for the leverage:
d(Yi , µb i )
∼ χ21 .
(1 − hi )ϕ
b

This implies that the absolute value of the standardized deviance residual has a half-normal
distribution. v
u d(Yi , µ bi)
u
(sd)
Ri = t ∼ N+ (0, 1).
(1 − h b )ϕ
i
b

When the saddlepoint approximation holds, we can compare the absolute value of the standard-
(sd)
ized deviance residuals ri with the quantiles of the half-normal distribution.
(sd)
Note that Ri does not, in general, have a normal distribution. We know that the deviance
residuals do not have expectation zero.

10.2.2 Distribution of the total deviance

For Poisson and binomial models with large counts, the total deviance has an asymptotic
chi-squared distribution:
n
= d(Yi , µb i ) χ2n−p .
X
D(Y , µ)
b ∼
i=1

Hence for large n we should have

D(y, µ)
b ≈ n − p.
If the residual deviance D(y, µ)
b is much larger than the residual degrees of freedom n − p then
this indicates poor fit.

10.2.3 Limitations
We cannot use the residual deviance to measure goodness of fit when the saddlepoint approxi-
mation does not apply, e.g. for binary outcomes.
We cannot use it for EDMs with a free dispersion parameter that must be estimated, e.g. normal,
gamma, inverse Gaussian.
142 CHAPTER 10. FURTHER TOPICS

10.3 Quasi-likelihood for GLMs

10.3.1 Estimating equations
So far we have thought of estimating β by maximizing a log likelihood.
An alternative perspective is that we are estimating β by solving the estimating equations
U (β)
b = 0.
where
1X n
U (β) = Wi (yi − µi ) g ′ (µi )xi
ϕ i=1
and µi , Wi are implicit functions of β.
It is sufficient to show that
E [U (β)] = 0
!
∂U
Var [U (β)] = E − T = I(β)
∂β
to have the asymptotic result  
β
b ∼ N β, I(β)−1
as n → ∞ (Wedderburn 1974).
We will show that these conditions hold if we correctly specify the mean and variance.
Proof:
Suppose that
E (Yi | xi ) = µi ,
ϕ
Var (Yi | xi ) = V (µi ).
wi

Then n
Wi g ′ (µi ) xi
!
E [U (β)] = E(Yi − µi ) =
X
0.
i=1 ϕ

Now,
∂U (β) n
Wi g ′ (µi )xi
! !
∂µi
=
X
E −
∂β i=1 ϕ ∂β

n
Wi g ′ (µi )xi xTi
!
=
X

i=1 ϕ g ′ (µi )

n
Wi xi xTi
=
X

i=1 ϕ

= I(β).
10.3. QUASI-LIKELIHOOD FOR GLMS 143

Furthermore,
!2
n
Wi g ′ (µi )
Var [U (β)] = xi xTi Var(Yi )
X

i=1 ϕ

n
Wi2 [g ′ (µi )]2 ϕV (µi )
! !
= xi xTi
X

i=1 ϕ2 wi

n
!
Wi
= xi xTi
X

i=1 ϕ

= I(β)

since
wi
Wi = .
V (µi ) [g ′ (µi )]2

10.3.2 Fitting a quasi-likelihood model

The large sample properties of β
b depend only on correct specification of the mean µ and the
i
variance function V (µi ).
This allows us to fit quasi-likelihood models where there is no EDM such that
ϕ
Var(Yi | xi ) = V (µi ).
wi
It also means that we only need to check that the assumptions about the mean and variance
are correct. We can do this by analyzing the residuals.

10.3.3 Unbiased estimating functions

The estimate β
b derived from the estimating equation

U (β)
b = 0

is consistent if the estimating function U (β) is unbiased, that is

E [U (β)] = 0.

This holds if E(Yi | xi ) = µi , even if we misspecify the variance function.

Some Revision Questions

1. Budworms are moth caterpillars that can destroy garden flowers. The output below shows
the results of fitting a model in R to a toxicology experiment where different doses of
insecticide was used to kill budworms. In addition to the dose ldose, the experiment also
recorded the sex of the budworm (M=male, F=female).

Call:
glm(formula = cbind(numalive, numdead) ~ -1 + sex + ldose, family = binomial,
data = budworm)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
sexF 3.4732 0.4685 7.413 1.23e-13 ***
sexM 2.3724 0.3855 6.154 7.56e-10 ***
ldose -1.0642 0.1311 -8.119 4.70e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 126.2269 on 12 degrees of freedom

Residual deviance: 6.7571 on 9 degrees of freedom
AIC: 42.867

Number of Fisher Scoring iterations: 4

a. Calculate value of ldose that will kill 50% of female budworms (This is called the
“LD50” dose, which is short for “lethal dose 50%”).
b. Calculate the probability that the same dose will kill a male budworm.
Solution This question is about a toxicology experiment which is a typical application of
binomial regression models. An unusual feature of this model is that it lacks a global intercept
term (This is what the -1 in the model formula means). In such cases, when there is a factor in
the model, there will be a separate intercept for each factor level. In this model we get separate
intercepts for female budworms (sexF) and male budworms (sexM).
a. The linear predictor for this model is

ηi = αF 1{sexi =F } + αM 1{sexi =M } + β × ldosei

The default link function for the binomial family is logit. So the mean µi is related to the
linear predictor by !
µi
log = ηi
1 − µi

144
10.3. QUASI-LIKELIHOOD FOR GLMS 145

To find the LD50 dose, which kills half the budworms, we want µi = 0.5. Hence
0.5
 
ηi = log = log (1) = 0.
0.5
Now we substitute the fitted values

αF = 3.4732,
β = −1.0642,

and solve the equation for the linear predictor

0 = 3.4732 − 1.0642 × ldosei

giving
ldosei = 3.4732/1.0642 = 3.26.
Thus a LD50 dose is 3.26 (units were not provided).
b. To calculate the probability of death we first calculate the linear predictor for male
budworms. (Note the use of the different intercept term αM ).

ηi = αM + β × ldosei
= 2.3724 − 1.0642 × 3.26
= −1.097.

Then apply the inverse link to find the mean

1
µi =
1 + exp(−ηi )
1
=
1 + exp(1.097)
= 0.25.

Finally, we must remember that the expectation µi is the probability that individual
i is alive (When we use the two-column outcome variable for the binomial model, the
expectation is the probability of the first event). Hence, the predicted probability that
the LD50 dose for a female budworm will kill a make budworm is p = 1 − 0.25 = 0.75.
If you ever get confused about which probability you are modelling (death or survival) in
a toxicology model then remember that more poison kills more things. Look at the
sign of the dose term. If it is negative then you are modelling the probability of survival
(decreasing probability of survival with dose). If it is positive then you are modelling the
probability of death (increasing probability of death with dose).
146 CHAPTER 10. FURTHER TOPICS

2. Consider the following box plots of outcome Y against a numeric variable X that takes
integer values from 4 to 10.
a. Suggest an appropriate variance function for modelling these data.
b. Suggest an appropriate link function for modelling these data.
Justify your answer in both cases.
12
10
8
y

6
4
2
0

4 5 6 7 8 9 10

Solution: This kind of question tests your ability to formulate a GLM based on a summary of
the data, such as a graph or table. Your answer does not have to be definitive. Some flexibility
is allowed in your answer as long as you justify your reasoning.
a. The mean of Y appears to be increasing with X, and so does the variance. This suggests
that the variance function should be a power function, i.e. either V (µ) = µ (Poisson) or
V (µ) = µ2 (gamma) are appropriate. We can rule out V (µ) = µ3 (inverse Gaussian) as
this implies a very rapid increase in variance with the mean, which is not visible in the
graph.
You may have spotted that the variable Y is not discrete. You can see this from the
individual outlier points plotted for X = 4. This is not sufficient reason to rule out the
Poisson variance function V (µ) = µ because this is shared by the quasi-Poisson family
Var(Y ) = ϕµ for ϕ ̸= 1, and the quasi-Poisson can be applied to non-integer data.
b. We can see from the graph that Y ≥ 0 for all observations. All the whiskers from the box
plot are non-negative and there are no negative outliers. This suggests that we need a
link function appropriate for non-negative data such as the log link.
10.3. QUASI-LIKELIHOOD FOR GLMS 147

3. The data set esoph concerns a study of esophageal cancer conducted in the Ille-et-Vilaine
department of France in the 1970s. Participants have been classified by age group (agegp),
alcohol consumption in g/day (alcgp) and tobacco consumption in g/day (tobgp). The
categories are:
• agegp: 25-34, 35-44, 45-54, 55-64, 65-74, 75+ years
• alcgp: 0-39, 40-79, 80-119, 120+ g/day
• tobgp: 0-9, 10-19, 20-29, 30+ g/day
In each class defined by the cross-classification of these variables, ncases gives the number
of participants in the study with esophageal cancer and ncontrols gives the number of
healthy controls without cancer.
The output on the below shows the results of fitting a model to these data in R.

Call:
glm(formula = cbind(ncases, ncontrols) ~ agegp + alcgp + tobgp,
family = binomial(), data = esoph)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.8954 1.0859 -6.350 2.16e-10
agegp35-44 1.9809 1.1041 1.794 0.072786
agegp45-54 3.7763 1.0680 3.536 0.000407
agegp55-64 4.3352 1.0650 4.070 4.69e-05
agegp65-74 4.8964 1.0764 4.549 5.39e-06
agegp75+ 4.8265 1.1213 4.304 1.67e-05
alcgp40-79 1.4346 0.2501 5.737 9.63e-09
alcgp80-119 1.9807 0.2848 6.956 3.51e-12
alcgp120+ 3.6029 0.3850 9.357 < 2e-16
tobgp10-19 0.4381 0.2283 1.919 0.055039
tobgp20-29 0.5126 0.2730 1.878 0.060398
tobgp30+ 1.6410 0.3441 4.769 1.85e-06

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 367.953 on 87 degrees of freedom

Residual deviance: 82.337 on 76 degrees of freedom
AIC: 221.39

Number of Fisher Scoring iterations: 6

a. Based on the output above, can you say how many rows are in the data set esoph?
b. Calculate a test statistic for testing the existence of regression, i.e. for testing this model
([Link]) against the null model with an intercept but no predictors. Give the distribution
of this test statistic under the null model.
c. Calculate the odds ratio for group tobgp30+ versus the reference category for tobacco
(which is omitted from the summary output but is labelled 0-9g/day in the data set
148 CHAPTER 10. FURTHER TOPICS

esoph). Calculate a 90% confidence interval for the odds ratio.

Solution:
a. The key to this question is the degrees of freedom for the null model, which has only an
intercept and no predictors. The null degrees of freedom is n − 1 = 87. Hence n = 88.
b. Testing the existence of regression compares the null model against the model that has
been fitted above. Here the null model is a binomial GLM which like the fitted model
uses the logit link function but its linear predictor consists only of an intercept term. It
is nested in the fitted model as the fitted model has an intercept term. For binomial
models the dispersion parameter is fixed at 1. Hence to test two nested models we take
the difference in deviance. Our test statistic is

∆ = Deviance1 − Deviance2
= 367.953 − 82.337
= 285.616

where model 2 is the larger model (with smaller deviance) and model 1 (here the null
model) is nested in model 2. Note that a GLM is nested in another, if both have the same
EDM, This has a χ2d distribution with

d = p2 − p1 = 12 − 1 = 11

degrees of freedom.
The question does not ask you to calculate the p-value, but it is clearly very small. A
difference of 285.6 on 11 degrees of freedom is a very large and highly significant difference.
c. This is a logistic regression model (binomial model with logit link). Hence for the predictor
variables, the coefficients are log odds ratios. To calculate the odds ratio we take the
exponential transformation

log OR = 1.6410
OR = exp(1.6410) = 5.16

The 100 × (1 − α)% confidence interval is calculated on the scale of the log odds ratio

Estimate ± Φ−1 (1 − α/2) × Std. Error = 1.6410 ± 1.65 × 0.3441

using the model output and Φ−1 (0.95) = 1.65.

This gives limits (1.073, 2.208). We then transform the confidence limits with the exponential
function to get the confidence interval (2.92, 9.10) for the odds ratio.
10.3. QUASI-LIKELIHOOD FOR GLMS 149

4. Continuing with the same data set we look at analysis of deviance. The single-argument
form of the anova() function tests the effect of adding each predictor in sequential order,
starting with the null model and ending with the full model with predictors agegp, alcgp,
and tobgp.
Since there are 3 predictors this gives us 3 p-values.
Analysis of Deviance Table

Model: binomial, link: logit

Response: cbind(ncases, ncontrols)

Terms added sequentially (first to last)

Df Deviance Resid. Df Resid. Dev Pr(>Chi)

NULL 87 367.95
agegp 5 121.045 82 246.91 < 2.2e-16
alcgp 3 141.028 79 105.88 < 2.2e-16
tobgp 3 23.544 76 82.34 3.11e-05
The two-argument form of anova() tests two nested models. Below we test the model
with agegp versus the model with agegp and alcgp as predictors.
Analysis of Deviance Table

Model 1: cbind(ncases, ncontrols) ~ agegp

Model 2: cbind(ncases, ncontrols) ~ agegp + alcgp
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 82 246.91
2 79 105.88 3 141.03 < 2.2e-16
Using the correspondence between the two tables above, produce a Analysis of Deviance
table that compares a model with agegp and alcgp versus the model with agegp, alcgp
and tobgp as predictors.
Model 1: cbind(ncases, ncontrols) ~ agegp + alcgp
Model 2: cbind(ncases, ncontrols) ~ agegp + alcgp + tobgp
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1
2

Solution:
This is less of a typical exam question and more a revision question designed to reinforce your
understanding of the two different forms of anova() in R.
• In the first form, there is a single argument which is a fitted model. The anova() function
tests a nested set of models starting with the null model (no predictors, but using the
150 CHAPTER 10. FURTHER TOPICS

same EDM and link function) and adding terms sequentially as they are given in the
model formula until the full model is reached.
• All models use the same EDM and link function. The full model contains an intercept
term and so has the null model nested in it.
• In the second form of anova(), there are two arguments which are both fitted models.
The first model is the smaller model nested in the second larger model. This is because
the first model uses the same EDM, the same link function and its predictor is a subset of
the predictors of the second model.
The key to understanding this question is that rows 2 and 3 of the first anova() output
correspond to the models
cbind(ncases, ncontrols) ~ agegp
and
cbind(ncases, ncontrols) ~ agegp + alcgp
which are the two models compared in the second anova() output. This second anova() output
contains exactly the same quantities found in rows 2 and 3 of the first output, but in a
different order.
Only one calculation is required. This is for column Df, which is the difference in the residual
degrees of freedom (Resid. Df) between the two models: Df = 82 − 79 = 3.
With this information, all you need to do is work out the correspondence between the two tables
and then apply the same principles to the third, empty table, which corresponds to rows 3 and
4 of the first anova() output. The completed table is
Analysis of Deviance Table

Model 1: cbind(ncases, ncontrols) ~ agegp + alcgp

Model 2: cbind(ncases, ncontrols) ~ agegp + alcgp + tobgp
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 79 105.881
2 76 82.337 3 23.544 3.11e-05
10.3. QUASI-LIKELIHOOD FOR GLMS 151

5. The density of the scaled binomial distribution with mean µ ∈ (0, 1) and sample size
m > 0 is !
m
p(y | µ, m) = µmy (1 − µ)m(1−y)
my

Recall that for the scaled binomial, considered as an EDM, the dispersion parameter ϕ = 1
and the prior weight w = m.
Write down the log likelihood of the the scaled binomial and use this to give an expression
for the unit deviance in terms of y and µ. Pay attention to the special cases y = 0 and
y = 1.
Solution: This is a typical bookwork question. Starting with the canonical form
!
w(θy − b(θ))
p(y | θ, ϕ) = a(y, ϕ/w) exp
ϕ

Let
t(y, µ) = θy − b(θ).
Then the unit deviance is defined as

d(y, µ) = 2 (t(y, y) − t(y, µ)) .

The first step is to put the log density into canonical form

log p(y | µ) = my log(µ) + m(1 − y) log(1 − µ) + C(m, y)

" ! #
µ
= m y log + log(1 − µ) + C(m, y)
1−µ

where C(m, y) is a normalizing term not depending on µ which we do not need to specify. This
is recognizable as an EDM in canonical form with
!
µ
θ = log ,
1−µ
b(θ) = − log(1 − µ),
ϕ = 1,
w = m.

It is important to recognize that m is the weight and is not included in the expression for the
unit deviance, which depends only on y and µ.
The expression for t(y, µ) is
!
µ
t(y, µ) = y log + log(1 − µ)
1−µ
= y log(µ) + (1 − y) log(1 − µ).
152 CHAPTER 10. FURTHER TOPICS

Hence

d(y, µ) = 2 (t(y, y) − t(y, µ))

1−y
! !
y
= 2y log + 2(1 − y) log .
µ 1−µ

The subtlety of unit deviance calculations is that the standard formula may not work on the
boundaries of the support of the distribution. For the binomial model, the above expression for
d(y, µ) contains terms in log(y) and log(1 − y) which are invalid for y = 0 and y = 1 respectively.
For the boundary cases we use the limiting formula
 
d(y, µ) = 2 lim
′
t(y, µ ) − t(y, µ) .
′
µ →y

For y = 0 we have

lim t(0, µ′ ) = lim log(1 − µ′ ) = log(1) = 0.

µ′ →0 µ′ →0

Hence  
d(0, µ) = 2 0 − t(0, µ) = −2 log(1 − µ).
Similarly, for y = 1 we have

lim t(1, µ′ ) = lim log(µ′ ) = log(1) = 0.

µ′ →1 µ′ →1

Hence  
d(1, µ) = 2 0 − t(1, µ) = −2 log(µ).

An important consistency check when you derive an expression for the unit deviance is that it
should be zero when y = µ. This can be confirmed for all three expressions given above.