5CCM242A / 6CCM242B

Statistical Modelling
Vasiliki Koutra
Department of Mathematics, King’s College London

Contents
1 Introduction to Statistical Modelling 2
1.1 Linear Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2 General linear model 3

3 Estimation 4
3.1 Least squares estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
3.2 Maximum likelihood estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.3 Properties of least squares estimators . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.4 Gauss–Markov Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

4 Inference in linear models 8

4.1 Tests and confidence intervals for the regression parameter . . . . . . . . . . . . . 8
4.2 Estimating σ 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
4.3 Sampling distribution of s2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
4.4 Hypothesis test on a single parameter . . . . . . . . . . . . . . . . . . . . . . . . 11
4.5 Confidence interval for an individual parameter . . . . . . . . . . . . . . . . . . . 12
4.6 Example: Portrait Studios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
4.7 Inference for the expected response . . . . . . . . . . . . . . . . . . . . . . . . . . 13

5 Prediction 14

6 Analysis of Variance (ANOVA) 14

6.1 Sampling Distribution of M SR and the Global F -Test . . . . . . . . . . . . . . . 15

7 Model Checking 16
7.1 Coefficient of Determination . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
7.2 Residual plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
7.3 Normal probability plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
7.4 Other residual plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
7.5 Transformation of the response . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
7.6 Box–Cox transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

1
 7.7 Outliers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

8 Problems in fitting regression models 22

8.1 Leverage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
8.2 Influence (Cook’s distance) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
8.3 Multicollinearity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

9 Common mistakes when interpreting fitted models 23

10 Model selection: functional marginality 23

1 Introduction to Statistical Modelling

Statistical modelling provides a language for turning data into evidence. The central idea is
simple: we represent the mechanism that produced the observations using a probability model,
and then use that model to learn about unknown quantities and to quantify uncertainty.
In most applications we do not observe the full mechanism that generates the data. Instead, we
observe a finite sample, often affected by measurement error, natural variability, and sampling
variation. A statistical model provides a structured way to connect the data we see to the
questions we want to answer.
These notes develop the main ideas of statistical modelling through the lens of linear models.
We begin by motivating why models are needed, how variables are classified, and how random
variation is represented mathematically.
We use data to learn about the process from which the data were collected. A statistical model,
which includes a random component, can be useful. We have seen simple models, e.g.
Yi ∼ N (µ, σ 2 ),
where y1 , y2 , . . . , yn are observed data, assumed to be realisations of random variables (r.v.s)
Y1 , Y2 , . . . , Yn and we assume Y1 , Y2 , . . . , Yn are independent. Typically, we will estimate the
parameters µ and σ 2 from the data, produce confidence intervals, test hypotheses, predict new
responses, etc.
In practice, most data sets involve several variables, which might be:
1. qualitative (categorical);
2. quantitative, either:
(a) discrete; or
(b) continuous.
We will assume that the primary interest is in one particular variable Y , the response variable,
which we are interested in explaining the variation in, for understanding the process, inference,
prediction etc.
The other variables X1 , X2 , . . . , Xq are explanatory variables, which might be useful for explain-
ing the variation in Y .

1.1 Linear Models

We start by considering a response variable Y which is continuous. It is natural to assume Yi

has a normal distribution, but instead of constant expectation, we assume dependence on the

2
levels (or categories) of one or more explanatory variables, e.g., simple linear regression (SLR):
ind
Yi ∼ N (β0 + β1 xi , σ 2 ), i = 1, 2, . . . , n .

This can be written

ind
Yi ∼ N (µi , σ 2 ), µi = β0 + β1 xi ,
iid
Yi = µi + εi , µi = β0 + β1 xi , εi ∼ N (0, σ 2 ) ,
iid
Yi = β0 + β1 xi + εi , εi ∼ N (0, σ 2 ) .
The εi are known as errors or departures.
This generalises to the normal linear model:
ind
Yi ∼ N (µi , σ 2 ),

µi = β1 φ1 (xi ) + · · · + βp φp (xi ),
where φj (xi ) are known basis functions. The model is linear in the parameters β1 , . . . , βp .
Examples:
• µi = β0 + β1 xi . This is simple linear regression with basis functions φ0 (x) = 1, φ1 (x) = x.
• µi = β0 + β1 xi + β11 x2i (φ0 (x) = 1, φ1 (x) = x, φ11 (x) = x2 ). This is quadratic regression
with basis functions φ0 (x) = 1, φ1 (x) = x, φ11 (x) = x2 .
• µi = β0 + β1 x1i + β2 x2i . This is multiple linear regression (with two explanatory variables)
with basis functions φ0 (x) = 1, φ1 (x) = x1 , φ2 (x) = x2 .
• µi = β0 + β1 Z1i + β2 Z2i with Z1i = 1 if Ti = 1, 0 otherwise; Z2i = 1 if Ti = 2, 0 otherwise.
This corresponds to a one-way ANOVA model with three groups (levels 0, 1, 2) and basis
functions φ0 (T ) = 1, φ1 (T ) = Z1 , φ2 (T ) = Z2 .

2 General linear model

We can relax the normality assumption and write the general linear model as

E(Yi ) = µi = β1 φ1 (xi ) + β2 φ2 (xi ) + · · · + βp φp (xi ) ,

with
Var(Yi ) = σ 2 , Cov(Yi , Yi0 ) = 0 .
We make no assumptions about the distribution of Yi beyond the first two moments.
The general linear model can be written

Y = Xβ + ε, E(ε) = 0, Var(ε) = σ 2 I ,

where
       
Y1 β1 ε1 φ1 (x1 ) φ2 (x1 ) · · · φp (x1 )
 Y2  β2   ε2   φ1 (x2 ) φ2 (x2 ) · · · φp (x2 ) 
Y =  . , β =  . , ε =  . , X= . ..  .
       
..
 ..   ..   ..   .. . . 
Yn βp εn φ1 (xn ) φ2 (xn ) · · · φp (xn )

The normal linear model,

Y ∼ Nn (Xβ, σ 2 I) ,

3
adds the assumption of normality for the errors, where Nn (µ, Σ) denotes an n-dimensional
multivariate normal distribution with expectation µ and variance-covariance matrix Σ.
Aside: A k × 1 random vector A with joint p.d.f. of A1 , A2 , . . . , Ak given by
 
−k/2 −1/2 1 > −1
fA (a; µ, Σ) = (2π) |Σ| exp − (a − µ) Σ (a − µ) ,
2

is said to have a k-dimensional multivariate normal distribution, denoted A ∼ Nk (µ, Σ).

It has E(A) = µ and Var(A) = Σ.
     
A1 µ1 Σ11 Σ12
If A = ,µ= and Σ = , then A1 ∼ Nk1 (µ1 , Σ11 ).
A2 µ2 Σ21 Σ22
Aside: a random vector is a vector whose elements are random variables. The expectation of
random vector Y is given by  
E(Y1 )
 E(Y2 ) 
E(Y) =  .  ,
 
 .. 
E(Yn )
with variance-covariance matrix
 
Var(Y1 ) Cov(Y1 , Y2 ) · · · Cov(Y1 , Yn )
 .. 
 Cov(Y2 , Y1 ) Var(Y2 ) ··· . 
Var(Y) = 
 .. ..
.

 . . 
Cov(Yn , Y1 ) ··· Var(Yn )

Here I is the n × n identity matrix. Note that Cov(Yi , Yj ) = Cov(Yj , Yi ), so the variance-
covariance matrix is symmetric.

3 Estimation
3.1 Least squares estimation

Given a general linear model

Y = Xβ + ε , E(ε) = 0 , Var(ε) = σ 2 I ,

we will obtain estimates β̂ of the unknown parameter β.

These imply estimates
\ = µ̂ ,
E(Y)
commonly labelled ŷ, of the expected responses given the observed X. These are

ŷ = Xβ̂.

The idea of least squares is that we want ŷi to be close to yi . We choose β to minimise
n
X
S(β) = (yi − µi )2 = (y − µ)> (y − µ) .
i=1

Substituting µ = Xβ,

S(β) = (y − Xβ)> (y − Xβ) = y> y − 2β > X> y + β > X> Xβ.

4
Aside: if u is a scalar and v is a vector
 ∂u 
∂v
∂u  . 1 
=  ..  .
∂v ∂u
∂vt

Let a be a constant vector and B be a constant symmetric matrix. It follows that

∂(v> a) ∂(v> Bv) ∂(Bv)

= a, = 2Bv, = B.
∂v ∂v ∂v

To find a stationary point of S(β) we need

∂S
= −2X> y + 2X> Xβ = 0 .
∂β
This give us the normal equations:
X> Xβ̂ = X> y
If X> X is non-singular,
β̂ = (X> X)−1 X> y.
The second derivative is 2X> X, which is positive definite if X> X is non-singular, so the sta-
tionary point is a minimum.

Example

Let E(Yi ) = β0 + β1 xi , with    

12 0
y = 10 ,
 x = 1 .

6 2
Then,
 
1 0      
> 3 3 > > −1 1 5 −3 5/6 −1/2
X = 1 1 . X X = , |X X| = 6, (X X) = = .
3 5 6 −3 3 −1/2 1/2
1 2
   
>28 37/3
X y= , β̂ = .
22 −3

Note: We have defined the least squares estimate

β̂ = (X> X)−1 X> y

for a given data vector y. The least squares estimator (LSE) is the random vector

β̂ = (X> X)−1 X> Y ,

where Y is the random vector of responses.

5
3.2 Maximum likelihood estimation

If we assume the normal linear model,

Y ∼ Nn (Xβ, σ 2 I) ,

the likelihood function is

 
1
L(β, σ 2 ; y) = (2πσ 2 )−n/2 exp − 2 (y − Xβ)> (y − Xβ) .
2σ

Maximising L(β, σ 2 ; y) is equivalent to minimising

(y − Xβ)> (y − Xβ) ,

which is the same as the least squares criterion. Hence, the maximum likelihood estimate (MLE)
of β is the same as the LSE, i.e.,

β̂ = (X> X)−1 X> y .

3.3 Properties of least squares estimators

Aside: For a constant n × n matrix A and random vector Y,

E(AY) = AE(Y), E(Y> A) = E(Y> )A, Var(AY) = AVar(Y)A> .

Lemma 3.1. β̂ is an unbiased estimator of β, i.e., E(β̂) = β.

Proof.

E(β̂) = E[(X> X)−1 X> Y] = (X> X)−1 X> E(Y) = (X> X)−1 X> Xβ = β.

Lemma 3.2. Var(β̂) = σ 2 (X> X)−1 .

Proof.

Var(β̂) = Var[(X> X)−1 X> Y] = (X> X)−1 X> Var(Y)X(X> X)−1 = σ 2 (X> X)−1 .

3.4 Gauss–Markov Theorem

We want to show that V ar(β̂) is “small”. We can use the Gauss–Markov theorem, which applies
to any estimable linear function of β,

`> β = `1 β1 + `2 β2 + · · · + `p βp ,

where ` = (`1 , `2 , . . . , `p )> is a known vector of constants. Note that `> β̂ is an unbiased
estimator of `> β and that
`> β̂ = `> (X> X)−1 X> Y
is a linear function of Y.

6
Theorem 3.1 (Gauss–Markov Theorem). Let `> β be an estimable function. Then `> β̂ has
minimum variance among all linear (in Y) unbiased estimators of `> β.

Note: ŷ = Xβ̂ = X(X> X)−1 X> Y = HY.

Definition 3.1. H = X(X> X)−1 X> is the hat matrix. The hat matrix is symmetric and
idempotent.
Lemma 3.3. I − H is symmetric and idempotent.

Proof. Symmetry:
(I − H)> = I> − H> = I − H .
Idempotency:
(I − H)2 = I − 2H + H2 = I − 2H + H = I − H .

Proof of Gauss–Markov Theorem

Proof. Let c> Y be another linear unbiased estimator of `> β. Then

E(c> Y) = c> E(Y) = c> Xβ = `> β ,

for all β so c> X = `> .

Now,
Var(c> Y) = c> Var(Y)c = σ 2 c> c .
Also,
Var(`> β̂) = σ 2 `> (X> X)−1 ` = σ 2 c> X(X> X)−1 X> c = σ 2 c> Hc .

Then,
Var(c> Y) − Var(`> β̂) = σ 2 c> (I − H)c .
Since I − H is symmetric and idempotent,

c> (I − H)c ≥ 0 ,

so
Var(c> Y) ≥ Var(`> β̂) .

Hence `> β̂ has minimum variance among all linear unbiased estimators of `> β.

We say that β̂ is a minimum variance linear unbiased estimator (MVLUE) of β, or a best linear
unbiased estimator (BLUE) of β.
Note:
• A special case of the Gauss–Markov theorem is that each β̂j has minimum variance among
all linear unbiased estimators of βj .
• µ̂(x0 = x> >
0 β̂ is a BLUE of µ(x0 ) = φ(x0 ) β, where
 
φ1 (x0 )
φ2 (x0 )
φ(x0 ) =  .  .
 
 .. 
φp (x0 )

7
4 Inference in linear models
Recall that the LSE is
β̂ = (X> X)−1 X> Y ,
with E(β̂) = β and Var(β̂) = σ 2 (X> X)−1 . Also note that β̂ is a linear function of Y.
Hence, if the normal linear model holds, i.e., Y ∼ Nn (Xβ, σ 2 I), then β̂ is a linear function of
a multivariate normal random vector, so β̂ is also multivariate normal:
 
β̂ ∼ Np β, σ 2 (X> X)−1 .

That is, we know sampling distribution of β̂ (up to the unknown constant σ 2 ).

For the general linear model, we do not know the sampling distribution of β̂. However, by the
multivariate Central Limit Theorem, for large n, β̂ follows approximately the same multivariate
normal distribution.

Practical implications

1. We usually do inference assuming normality.

2. We should check the validity of this assumption.
3. CLT ensures inference is robust to departures from normality in large samples.

4.1 Tests and confidence intervals for the regression parameter

In the normal linear and, asymptotically as n → ∞, in the general linear model,

 
β̂ ∼ Np β, σ 2 (X> X)−1 .

Let V = (X> X)−1 . Hence

β̂j − βj
∼ N (0, 1), j = 1, 2, . . . , p
σ 2 vjj

This could be used to construct hypothesis tests and confidence intervals (CIs) for βj , if σ 2 were
known. However, in practice σ 2 is unknown.

4.2 Estimating σ 2

It seems reasonable to base the estimate of σ 2 on the attained minimum value of S(β), namely

S(β̂) = (y − Xβ̂)> (y − Xβ̂) .

S(β̂) is called the residual sum of squares (RSS or SSRes ) or error sum of squares (SSE ).
As usual, we also define a random variable,

SSE = (Y − Xβ̂)> (Y − Xβ̂) = Y> (I − H)Y .

Theorem 4.1. E(SSE ) = (n − p)σ 2 .

8
Proof.
E(SSE ) = E[Y> (I − H)Y] = E[tr(Y> (I − H)Y)] .

Aside: Recall for matrices A and B of appropriate size, tr(AB) = tr(BA). Also, Var(Y) =
E(YY> ) − E(Y)E(Y)> .
Hence,
E(SSE ) = tr[(I − H)E(YY> )] = tr[(I − H)(Var(Y) + E(Y)E(Y)> )] .
Since Var(Y) = σ 2 I and E(Y) = Xβ,

E(SSE ) = tr[(I − H)(σ 2 I + Xββ > X> )] .

Using properties of the trace operator,

E(SSE ) = σ 2 tr(I − H) + tr[(I − H)Xββ > X> ] .

But HX = X, so
tr[(I − H)Xββ > X> ] = 0 .
Also, tr(I − H) = n − tr(H) = n − p. Thus,

E(SSE ) = (n − p)σ 2 .

SSE
Corollary 4.1. n−p is an unbiased estimator of σ 2 .

Definition 4.1.
SSE
s2 = MSE =
n−p
is called the mean square error or residual mean square.

Note: s2 is not the sample variance of the yi , unless the model is yi = β0 + εi .

Note: It can be shown that the MLE of σ 2 is SSE /n, but this is a biased estimator of σ 2 . In
practice, we use the unbiased estimator s2 = SSE /(n − p).
Now we can estimate Var(β̂j ) by
d β̂j ) = s2 vjj ,
Var(
so an obvious test statistic to test H0 : βj = βj0 is

β̂j − βj0
p = Tj .
s2 vjj

To carry out hypothesis tests, we need the sampling distribution of s2 .

4.3 Sampling distribution of s2

Definition 4.2. ei = yi − ŷi is the ith residual, and e = y − ŷ is the vector of residuals.

Note:
1. e = y − ŷ = y − Hy = (I − H)y.
2. We use e to denote both the random vector (I − H)Y and its realised value.

9
It follows that
SSE = e> e .

Lemma 4.1. e = (I − H)ε.

Proof.
e = Y − Ŷ = Xβ + ε − H(Xβ + ε) = (I − H)ε .

Note: ei can be thought of as an estimator of εi .

Corollary 4.2. e ∼ Nn (0, σ 2 (I − H)).

Proof.
1. Each ei is a linear function of the normal random vector ε, so e is multivariate normal.
2. E(e) = E[(I − H)ε] = (I − H)E(ε) = 0.
3. Var(e) = (I − H)Var(ε)(I − H)> = σ 2 (I − H).

Two results from linear algebra

Lemma 4.2. For a symmetric idempotent matrix A of rank r, there exists an orthogonal matrix
C such that A = CDC> , where  
Ir 0
D= .
0 0

Lemma 4.3. For an idempotent matrix A, tr(A) = rank(A).

Corollary 4.3. rank(I − H) = n − p.

Theorem 4.2 (Cochran’s Theorem).

(n − p)s2
∼ χ2n−p .
σ2
Proof.
(n − p)s2 SSE e> e ε> (I − H)ε
= = = .
σ2 σ2 σ2 σ2
Hence, we can write
(n − p)s2 ε> CDC> ε (C> ε)> D(C> ε)
= = .
σ2 σ2 σ2
Let Z = C> ε/σ. Then,
(n − p)s2
= Z> DZ .
σ2
As ε ∼ Nn (0, σ 2 I) and C is orthogonal, Z ∼ Nn (0, I). Thus,
n−p
(n − p)s2 X 2
= Zi ∼ χ2n−p .
σ2
i=1

10
Aside: A chi-squared random variable with r degrees of freedom is the distribution of the sum of
squares of r independent standard normal random variables. If Z1 , Z2 , . . . , Zr ∼ i.i.d. N (0, 1),
then
Xk
Zi2 ∼ χ2r .
i=1

Lemma 4.4. If z ∼ N (0, 1) and w ∼ χ2ν are independent, then

z
p ∼ tν .
w/ν

We have
β̂j − βj
p ∼ N (0, 1) ,
σ 2 vjj
and, independently,
(n − p)s2
∼ χ2n−p .
σ2
Hence,
β̂j −βj
√
σ 2 vjj β̂j − βj
q =p ∼ tn−p .
(n−p)s2
2
s2 vjj
σ (n−p)

This result is used to obtain confidence intervals for βj and to test hypotheses such as H0j :
p
βj = 0 against H1j : βj 6= 0. Recall that s2 vjj is the estimated standard error (s.e.) of β̂j .

4.4 Hypothesis test on a single parameter

We can test H0j : βj = βj0 , usually βj0 = 0, using the test statistic

β̂j − βj0
Tj = p .
s2 vjj

If the null hypothesis is true, then Tj ∼ tn−p .

Note: R (and other statistical software) gives this test with βj0 = 0 for each j.
For a one-sided alternative H1j : βj > βj0 , reject H0j at the 100α% level of significance if tj ,
the obbserved value of Tj , has tj ≥ tn−p,1−α . Or for H0j : βj < βj0 , reject H0j if tj ≤ −tn−p,1−α .
For a two-sided alternative H1j : βj 6= βj0 , reject H0j at the 100α% level of significance if
|tj | ≥ tn−p,1−α/2 .
Alternatively, we can calculate the p-value of the test, for observed tj as:
• 1-sided >: P (Tj ≥ tj | H0j true)
• 1-sided <: P (Tj ≤ tj | H0j true)
• 2-sided: 2P (Tj ≥ |tj | | H0j true)
Recall: The p-value is the probability of observing a value of the test statistic at least as extreme
as that calculated from our data, given the null hypothesis is true.
Also, for a p-value < α, reject H0j at the 100α% level of significance.

11
4.5 Confidence interval for an individual parameter

Since
β̂j − βj
p ∼ tn−p ,
s2 vjj
we have: !
β̂j − βj
P −tn−p, 1−α/2 ≤p ≤ tn−p, 1−α/2 = 1 − α.
s2 vjj
Rearranging:
 q q 
P β̂j − tn−p, 1−α/2 s2 vjj ≤ βj ≤ β̂j + tn−p, 1−α/2 s2 vjj . = 1 − α .

p
Since s.e.(β̂j ) = s2 vjj , the 100(1 − α)% CI is:
h i
β̂j − tn−p, 1−α/2 s.e.(β̂j ), β̂j + tn−p, 1−α/2 s.e.(β̂j ) .

4.6 Example: Portrait Studios

A company operates portrait studios, specialising in children’s portraits, in 21 cities. They want
to investigate whether sales (Y , $’000) depend on the number of people under 16 in the city
(x1 , ’000) and per capita income (x2 , $’000). They have data for each city.
The multiple regression model:

Yi = β0 + β1 x1i + β2 x2i + εi

was fitted, producing:

Parameter Estimate s.e. (Estimate)
β0 −68.8571 60.0170
β1 1.4546 0.2118
β2 9.3655 4.0640

To test H01 : β1 = 0 vs HA1 : β1 6= 0, we use

β̂1 1.4546
t1 = = = 6.8678
s.e.(β̂1 ) 0.2118

If H01 is true, T1 ∼ tn−p with n = 21, p = 3 (i.e. t18 distribution).

The p-value is:

P (|T1 | > 6.8678) = 2P (T1 > 6.8678) = 2 × 0.0000 = 0.0000 (3 d.p.)

There is overwhelming evidence against H01 , so sales are related to the number of children in
the city.
To test H02 : β2 = 0 vs HA2 : β2 6= 0, we use:

β̂2 9.3655
t2 = = = 2.3045
s.e.(β̂2 ) 4.0640

Critical values:

t0.975,18 = 2.101 (α = 5%), t0.995,18 = 2.878 (α = 1%).

12
Since 2.101 < 2.3045 < 2.878, we reject H02 at α = 5% but not at α = 1%.
The p-value is:
p = 2 [1 − P (T2 ≤ 2.3045)] = 2 [1 − 0.9833] = 0.0334.
Thus, there is moderate evidence that income is related to sales.
Note: These tests only make sense for this model. We are asking:

(i) Does the number of children have any relationship with sales, given the relationship of
income with sales?
That is, we are comparing the models:

E(Yi ) = β0 + β1 x1i + β2 x2i with E(Yi ) = β0 + β2 x2i

not
E(Yi ) = β0 + β1 x1i with E(Yi ) = β0

(ii) Does income have any relationship with sales, beyond what is explained by the number
of children?

A 95% CI for β1 is
β̂1 ± tn−p, 1−α/2 s.e.(β̂1 )
1.4546 ± t18, 0.975 × 0.2118
1.4546 ± 2.101 × 0.2118
[1.0096, 1.8996]
Each increase of 1 child (in thousands) in a city is associated with increased sales of between
$1.01k and $1.90k.

4.7 Inference for the expected response

As µ0 = µ(x0 ) = x>0 β. inferences for the expected response at some particular values, x0 , of
the explanatory variables are easily obtained. The point estimator of µ0 is

µ̂0 = x>
0 β̂ .

with sampling distribution being normal with

E(µ̂0 ) = x>
0 β = µ0 ,
Var(µ̂0 ) = Var(x> 2 > > −1
0 β) = σ x0 (X X) x0 .

Hence, µ̂0 ∼ N µ0 , σ 2 x> > −1

0 (X X) x0 .
Replacing σ 2 by s2 = M SE yields the test statistic for H0 : µ0 = µ? ,
µ̂ − µ?
T0 = q 0 ∼ tn−p (if H0 is true).
s x>
0 (X > X)−1 x
0

A 100(1 − α)% confidence interval for µ0 is

q
µ̂0 ± tn−p; 1−α/2 s x> > −1
0 (X X) x0 .

13
5 Prediction
To predict an actual future observation Y0 at x0 , we use Ŷ0 = x> 0 β̂ = µ̂0 . Because we must
account for the variation in the observation, a confidence interval for µ0 is not adequate; instead
we use a prediction interval.
We have
Y0 ∼ N (µ0 , σ 2 ) ,
but µ0 is unknown. We also have

µ̂0 ∼ N (µ0 , σ 2 x> > −1

0 (X X) x0 ) .

Hence  
Ŷ0 − Y0 = µ̂0 − Y0 ∼ N 0, σ 2 x> > −1

0 (X X) x0 + 1 ,

as Y0 is independent of the sample Y1 , . . . , Yn (and hence independent of µ̂0 ).

Therefore,
Ŷ0 − Y0
q ∼ N (0, 1) .
σ x> > −1
0 (X X) x0 + 1

Then,
Ŷ0 − Y0
q ∼ tn−p ,
> > −1
s x0 (X X) x0 + 1
and a 100(1 − α)% prediction interval for Y0 is
q
Ŷ0 ± tn−p; 1−α/2 s x> > −1
0 (X X) x0 + 1 .

6 Analysis of Variance (ANOVA)

Recall that we have defined

SSE = (y − Xβ̂)> (y − Xβ̂) = y> (I − H)y ,

where H = X(X> X)−1 X> .

If wePnignore the explanatory variables, the best predictor of yi is the sample mean ȳ =
n −1
i=1 yi . The sample variance of the yi is
n
1 X
s2y = (yi − ȳ)2 .
n−1
i=1

Definition 6.1. The total sum of squares is

n
X
SST = (yi − ȳ)2 = y> y − nȳ 2 .
i=1

Definition 6.2. The regression sum of squares is

n
X
SSR = (ŷi − ȳ)2 = ŷ> ŷ − nȳ 2 = y> Hy − nȳ 2 .
i=1

These results give the analysis of variance identity.

14
Theorem 6.1 (ANOVA identity).

SST = SSR + SSE .

Total variation = Explained variation + Residual variation.

Proof. Using ŷ = Hy and I = H + (I − H),

SSR + SSE = y> Hy − nȳ 2 + y> (I − H)y

= y> H + I − H y − nȳ 2 = y> y − nȳ 2 = SST .

We know that SST is related to a χ2n−1 distribution, and SSE is related to a χ2n−p distribution.
We can interpret the corresponding degrees of freedom (d.f.) as “pieces of information” contained
in the sums of squares. The n observations contribute n − 1 pieces of information about the
total variability (since we estimate ȳ from the data). Fitting the regression model uses up p − 1
pieces of information (since we estimate p parameters, including the intercept), leaving n − p
pieces of information to estimate the residual variability.
We present all of this in the analysis of variance table:

Source d.f. Sum of Squares Mean Square

Regression (Model) p−1 SSR M SR = SSR /(p − 1)
Error (Residual) n−p SSE M SE = SSE /(n − p)
Total n−1 SST

The variance ratio (F-ratio) is F = M SR /M SE .

Notes.
1. The total mean square would equal the sample variance of the yi but is not usually shown
in the ANOVA table.
2. The variance ratio (F -ratio) measures model variance explained relative to the error (ran-
dom) variance between observations.

6.1 Sampling Distribution of M SR and the Global F -Test

Let β = (β0 , β1 , . . . , βp−1 )> and denote by β1 = (β1 , . . . , βp−1 )> the coefficients of the explana-
tory variables (excluding the intercept).

Theorem 6.2. If β1 = 0, then SSR /σ 2 ∼ χ2p−1 .

(p − 1)M SR SSR
Corollary 6.1. If β1 = 0, then 2
= 2 ∼ χ2p−1 .
σ σ

(W1 /ν1 )
Aside (distribution theory). If W1 ∼ χ2ν1 and W2 ∼ χ2ν2 are independent, then ∼ Fν1 ,ν2 .
(W2 /ν2 )
Since
(p − 1)M SR M SR
2
/(p − 1) =
σ σ2
and
(n − p)M SE M SE
2
/(n − p) = ,
σ σ2
the following result follows.

15
Theorem 6.3. If β1 = 0, then M SR and M SE are independent, and
M SR
F = ∼ Fp−1, n−p .
M SE

We use this result to test the null hypothesis

H0 : β1 = 0 (i.e., no relationship between Y and any Xj ) ,
against the alternative
H1 : β1 6= 0 (i.e., at least one Xj is related to Y ) .

This is called the (global) F -test (named after R. A. Fisher, who formalised it in 1921). The
test statistic F = M SR /M SE is the variance ratio in the ANOVA table.
Under H0 , F ∼ Fp−1, n−p . We reject H0 at level α if the observed value of F is sufficiently large.
Specifically, reject H0 at level α if
F ≥ Fp−1, n−p; 1−α .
In practice the global F -test is performed before testing individual parameters.
1. If H0 is not rejected, no further analysis is usually needed.
2. More formally, this guards the Type-I error rate against multiple testing (a single test
with P (Type I error) = α).

7 Model Checking
We have seen how to fit models and do inference on parameters. However, this estimation
and inference is conditional on an assumed model. An obvious question is whether the model
is appropriate for the data. That is, does it capture the main features of the process which
generated the data. In particular, are the assumptions of the form of the mean relationship,
constant variance, uncorrelated observations and normality at least a reasonable approximation?

7.1 Coefficient of Determination

The coefficient of determination, often called R2 , is

SSR SSE
R2 = =1− ,
SST SST
and is often expressed as a percentage 100 SS
SST . As SSR , SST ≥ 0 and SSR ≤ SST , we have
R

2
0 ≤ R ≤ 1.
For simple linear regression (only), R2 equals the squared sample correlation between x and
y. In general, R2 is interpreted as the proportion of variability in y1 , . . . , yn explained by the
model fitted.
An R2 ≥ 90% is considered good and R2 ≥ 70% satisfactory, if all other aspects of the model
are acceptable.
If R2 ≤ 70%, the model will have low predictive power but may still indicate relationships
between E(Y ) and the explanatory variables. Note that R2 can always be inceased by adding
more terms to the model and any model with n parameters has R2 = 100%.
A large R2 does not mean the model is good and a low R2 does not mean the model is poor. It
only tells us about one aspect of model fit.

16
7.2 Residual plots

We have defined the (raw) residuals as ei = yi − ŷi , i = 1, . . . , n, and showed that E(ei ) =
0 and Var(ei ) = (1 − hii )σ 2 , where hii is the ith diagonal element of the hat matrix H =
X(X> X)−1 X> .

Definition 7.1. Standardised (studenised) residuals are given by

ei
ri = p . i = 1, . . . , n .
(1 − hii ) s2

Aproximately, ri ∼ N (0, 1). The ri are correlated but for large samples, hij (i 6= j) is small so
that approximately they behave like i.i.d. draws for large.
We can check the assumption of constant variance, against the alternative that it varies with
the expected response, by plotting ri against ŷi . We should see a random scatter of points on
either side of ri = 0; see Figure 1.

ri

ri = 0
ŷi

Figure 1: Residuals vs fitted values: ri against ŷi : random scatter about 0 indicates constant
variance.

If the variance increases with expectation, we will see more scatter for larger ŷi (a funnel shape);
see Figure 2.

ri

ri = 0
ŷi

Figure 2: Residuals vs fitted values: ri against ŷi showing increasing variance (funnel shape).

17
7.3 Normal probability plot

To check if a collection of observation, or quantities calculated from observations, come from a

normal distribution, we could use a histogram. However, it is difficult to see small but important
depatures from normality using such a plot. The normal probability (or q-q) plot is an
alternative.
Assume we have observed, or calculated, w1 , . . . , wn and wish to check the assumption wi ∼
2 ). We do not know µ or σ 2 , so we estimate them by the sample mean and standard
N (µw , σw w w
deviation, and start by standardising, i.e., define
wi − w̄
zi = , i = 1, . . . , n .
sw
We now need to check if Zi ∼ N (0, 1) with all Zi independent.
Let z[ k] denote the kth order statistic of the zi , i.e., the kth smallest value. So z[ 1], . . . , z[ n] are
the zi arranged in increasing order.
What is the expected value of Z[k] ? A good approximation is given by noting that if Z ∼ N (0, 1),
then P (Z ≤ z) = Φ(z), where Φ is the standard normal c.d.f. If we set P (Z ≤ z[ k]) = (k−0.5)/n,
i.e., the proportion of observations less than or equal to z[ k] (with a continuity correction). This
gives  
−1 k − 0.5
z[ k] = Φ
n
Order statistics tend to have quite concentrated distributions, so the expected value of Z[k] is
close to this value. Thus, a good approximation to E{Z[k] } is
 
−1 k − 0.5
E{Z[k] } ≈ Φ , k = 1, . . . , n
n

We plot z(k) against Φ−1 ((k − 12 )/n). If the zi are approximately standard normal, the points
should lie approximately on the 45◦ line through the origin. Departures from this line indicate
departures from normality.
A normal probability plot of r[ k], the ordered standardised residuals, allows us to check the
assumption of normalty; see Figure 3.

7.4 Other residual plots

We can plot the standardised residuals against the explanatory variables to check for non-
linearity in the mean relationship or to see if the variance changes with any of the explanatory
variables; see Figure 4.
If the data are collected over time, we can plot the standardised residuals against time to check
for evidence of time trends or autocorrelation; see Figure 5.
Note:
1. We are looking for large departures from assumptions.
2. Variation on the horizontal axis tells us nothing about model adequacy.
3. We should produce these plots routinely when fitting linear models.

7.5 Transformation of the response

If we detect non-constant variance (heteroskedasticity), a response transformation might help.

In particular, if the variance increases with the mean of the response, for strictly positive

18
 Z(k)
Normal Positive skewness

Φ−1 ((k − 12 )/n)

Heavy tails / high kurtosis

Figure 3: Normal probability plot sketches: (top-left) approximately normal; (top-right) positive
skew; (bottom) heavy tails. The red line is the 45◦ reference.

responses Yi > 0 we can use:

√
1. Yi if Var(Yi ) ∝ E(Yi )
2. log(Yi ) if Var(Yi ) ∝ [E(Yi )]2
√
3. 1/ Yi if Var(Yi ) ∝ [E(Yi )]3
4. 1/Yi if Var(Yi ) ∝ [E(Yi )]4
These transformations are ordered from less to more extreme.

7.6 Box–Cox transformation

These ad-hoc choices can be embedded in the Box–Cox family of transformations:

 λ
 Yi − 1 , λ 6= 0,

(λ)
Yi = λ

log Y ,
i λ = 0,
where λ is an additional parameter to be estimated from the data.
The easiest way to estimate λ is to calculate the log-likelihood log L(β, σ 2 ; λ, y) for a grid of λ
values, and pick the maximiser λ̂. The log-likelihood is calculated by fitting the linear model
Y (λ) = Xβ + ε, ε ∼ N (0, σ 2 I) .

We often round λ̂ to a simple rational number for [Link] way to judge how to round
is to construct a 95% confidence interval for λ using the profile log-likelihood; any convenient
value inside this interval is acceptable.

19
(i) No evidence of nonlinearity or nonconstant variance (ii) Quadratic nonlinearity
ri∗ ri∗

r∗ = 0 r∗ = 0
xi xi

(iii) Non-constant variance (funnelling)

ri∗

r∗ = 0
xi

Figure 4: Three diagnostic plots of standardised residuals ri∗ against an explanatory variable
x: (i) no evidence of nonlinearity or nonconstant variance; (ii) quadratic nonlinearity; (iii)
non-constant variance (funnelling).

To see the effect of the transformaton, consider λ = 0, i.e., the log transformation:

log Yi = β0 + β1 x1i + · · · + βq xqi + εi .

Exponentiating both sides gives

 q
X 
Yi = exp β0 + βk xki exp(εi ).
k=1

Thus the errors are multiplicative on the original scale; the variance increases with the expec-
tation; and if εi ∼ N (0, σ 2 ) then Yi is log-normal (positively skewed).
√
Occasionally, other transformations are used, e.g., sin−1 ( Yi ) for proportions.

7.7 Outliers

An outlier is an observation whose response does not fit the pattern of the rest of the data. In
linear models an outlier has a large positive or negative residual. Some software flag potential

20
 (i) ok (ii) time trend
ri∗ ri∗

r∗ = 0 r∗ = 0
ti ti

(iii) autocorrelation
ri∗

r∗ = 0
ti

Figure 5: Three diagnostic plots of standardised residuals ri∗ against time t: (i) no evidence of
time trend or autocorrelation; (ii) time trend; (iii) autocorrelation.

outliers with |ri | > 2 or |ri | > 3 (with standardised residuals ri ); in large samples such rules
flag too many points. Outliers will show up in any of the residual plots.
If we find an outlier, we should:
1. check whether the observation has been recorded or copied incorrectly; if so, correct or
delete it;
2. check if the sampled unit is somehow different in kind from the rest; if so, remove the
observation but note that the analysis then applies only to the main population;
3. if neither applies, repeat the analysis excluding the observation and report both analyses
if conclusions differ.
Be wary of labelling outliers: some observation will always have the largest residual due to
random variation; the ones that matter are those that actually change conclusions. If the model
is inadequate, spurious outliers may appear.

21
8 Problems in fitting regression models
8.1 Leverage

Observations which have unusual x values do not indicate anything wrong with the model or
the data (as long as it has been recorded correctly) but can affect the fitted model.
Let H = X(X> X)−1 X> be the hat matrix. Then ŷ = Hy, so for the ith fitted value
n
X X
ŷi = hij yj = hii yi + hij yj .
j=1 j6=i

Usually hij , i 6= j, are small relative to hii , hii indicates how much yi contributes to ŷi and is
called the leverage of (xi , yi ). It measures the remoteness of xi in the covariate space.
Note that
n
X
hii = tr(H) = tr X(X> X)−1 X> = tr(Ip ) = p.
 

i=1
Hence the mean leverage is p/n. If hii > 2p/n the leverage is large enough to warrant attention;
hii > 3p/n is very large.
A point with high leverage is potentially influential (it may exert a disproportionate effect on
the fitted line). Actual influence is measured using Cook’s distance.

8.2 Influence (Cook’s distance)

Cook’s distance for the ith observation is defined as

n
1 X (i)
2
Di = ŷj − ŷj ,
ps2
j=1

(i)
where ŷj is the fitted value for yj when the ith observation is removed, s2 is the mean square
error, and p the number of parameters (including the intercept). A large Di indicates that the
ith observation is influential.
A useful identity for the leave-one-out estimator relates β̂ (i) , the estimate of β with the ith
observation removed, to β̂, the estimate using all data:
ei (X> X)−1 xi
β̂ (i) = β̂ − .
1 − hii
From this expression, we can derive an alternative formula for Cook’s distance:
ri2 hii
Di = ,
p (1 − hii )
p
where ri = ei / (1 − hii ) s2 is the standardised residual.
The influence of the ith observation depends on both the residual size and the leverage. An
influential observation signals a problem caused by the data and cannot be “fixed” by further
analysis.

8.3 Multicollinearity

If X> X is singular we cannot obtain unique least squares estimates. Singularity arises when
one column of X is a linear combination of other columns (column rank < p). If X> X is nearly
singular (one column is nearly a linear combination of others), the following problems arise:

22
 1. numerical procedures for inverting X> X can fail or be inaccurate/unstable;
2. estimates are sensitive to small changes in the data (e.g. rounding);
3. some estimators are highly correlated so different models may fit equally well;
4. parameter estimators have large variances.
One way to detect this multicollinearity is to compute the variance inflation factor (VIF) for
each covariate. For j ∈ {1, . . . , q}, regress xj on the remaining covariates x1 , . . . , xj−1 , xj+1 , . . . , xq
and let Rj2 be the resulting coefficient of multiple determination. Then
1
VIFj = .
1 − Rj2
If Rj2 = 0 (i.e. xj is orthogonal to the other covariates), VIFj = 1; if Rj2 = 1 (i.e. xj is an exact
linear combination of the others), VIFj = ∞. A common rule of thumb is that multicollinearity
is serious if VIFj > 10 for any j ∈ {1, . . . , q}. This is a property of the data, not a flaw of the
model itself.
One suggested solution for multicolinearity is ridge regression, where we use the estimator
β̂R = (X> X + λI)−1 X> y, λ > 0,
with λ chosen as small as possible while avoiding large variance inflation. A more formal ap-
proach is to use Bayesian linear models, where prior or external information on β is incorporated
through a prior distribution.

9 Common mistakes when interpreting fitted models

1. Drawing causal conclusions when the study only supports descriptive relationships. Causal
inference requires either an experiment or a carefully agreed model for which plausible
causal pathways are encoded and justified.
2. Extrapolation beyond the range of the data, including hidden extrapolation. The model
should only be used to draw conclusions for covariate values x that are close to those
observed.
3. Extrapolation to different conditions from those under which the data were collected.

10 Model selection: functional marginality

If we detect problems with a fitted model, or are unsure which model to fit, we may use the
data for model selection. One aspect is choosing the linear predictor: which variables to include
(variable selection) and which terms involving these variables (linear, quadratic, interactions,
etc.).
We should insist on functional marginality: if any term is in the model, all terms marginal to
it must also be in the model. Hence
1. the model must include an intercept;
2. if the model includes x2j , it must include xj ;
3. if the model includes an interaction xj xk , it must include both xj and xk .
A simple reason is that rescaling a covariate should not change the model specification.
Model selection will continue in the second part of the module.

23