Linear Mixed Models

Chapter 1: Random Effects

Craig Anderson
 Course info

Lectures, tutorials, labs

20 lectures: 10 am Mon (Maths & Stats 116) & 11 am Tue

(5 The Square 330)
4 tutorials: Every two weeks starting week of 8/10
2 labs: 12-2 pm Tues 16/10 & 13/11 (Boyd Orr 418)

Lecturer

Craig Anderson, Room 343, Maths & Stats Building

Email: craig.anderson@glasgow.ac.uk
Twitter: @CAndersonStats

2/86
 Recommended textbooks

Ch. 4 of Ruppert, Wand & Carroll. Semiparametric

Regression. Cambridge University Press, 2003.
Ch. 10 of Venables and Ripley. Modern Applied Statistics
with S. 4th edition, Springer, 2002.
Faraway. Extending the Linear Model with R. Chapman &
Hall/CRC, 2006.
Pinheiro & Bates. Mixed-effects models in S and S-PLUS,
Springer, 2000.

All are available electronically via the library website.

3/86
 Aims and intended learning outcomes

Aims
To introduce students to Gaussian linear mixed effects
modelling in general and to the active use of modern mixed
effects modelling software.

Intended learning outcomes

By the end of the course students will be able to:
explain the notion of a random effect, why and when it is
useful and, in particular, how it differs from a fixed effect;
demonstrate detailed understanding of the theory
underpinning simple mixed effects models for balanced
designs;

4/86
Intended learning outcomes

demonstrate general understanding of the theory

underpinning general (e.g. unbalanced) mixed effects
models;
apply their knowledge of mixed effects modelling to
practical situations through the use of general mixed
effects modelling software;
check the assumptions of mixed effects models, both
graphically and by hypothesis testing based model
comparison;
explain when to use restricted maximum likelihood
(REML) when fitting mixed-effects models;

5/86
Intended learning outcomes

extend the treatment to cover scenarios involving general

covariance structures; interpret the output of R procedures
for linear mixed-effects models.
describe the basic ideas of multilevel models and of
generalised estimating equations as applied especially to
models for longitudinal data;
translate fluently between verbal, mathematical and
computational descriptions of models;
critically interpret analyses based on mixed effect models;
interpret the output of R procedures for linear
mixed-effects models.

6/86
 Course outline

Topics

Review of linear models

Random effects
Mixed models
Estimation
Prediction
Inference, diagnostics
Repeated measures and longitudinal data
Generalised estimating equations

7/86
 Revision
One-way ANOVA
One-way analysis of variance describes how the mean of a
continuous response variable depends on a categorical
explanatory variable.
Goal: to estimate the population group means µ1 , . . . , µI
and to test the hypothesis that all means are equal.
H0 : µ1 = · · · = µI
H1 : not all the µi are equal.
Note: some of the µi could be equal under H1 .
Main questions of interest:
1. Are there any differences in the I population means?
2. If yes, which of the means differ?

8/86
 Revision

Data structure in one-way ANOVA

Group 1 Group 2 . . . Group I

y11 y21 ... yI1
y12 y22 ... yI2
.. .. ..
. . .
y1n1 y2n2 ... yInI

Data
yij is the jth observation of the ith sample. (i = 1, . . . , I and
j = 1, . . . , ni ).

9/86
 Terminology

Factor
Any categorical variable such as treatment (trt1 , . . . , trtI ).

Levels
The categories of a factor, e.g. trt2 is one particular level of the
treatment factor.

Effects
The differences in the mean of the response variable among
different levels of a factor.

10/86
 Analysis of variance

The goal of one-way ANOVA is to assess whether a factor

has a significant effect on the response, Y.

Two ways of looking at this:

looking for differences in the mean of Y among the levels
of a factor or
assessing whether such differences explain some of the
variability in Y.

11/86
 Analysis of variance

Model assumptions

Yij : independent random variables for all i and j

Yij ∼ N(µi , σ 2 ) (common variance for all populations).
This can be written as a linear model as

yij = µ + αi + eij

where µi = µ + αi and the eij are independent, identically

distributed (i.i.d.) N(0, σ 2 ).

12/86
 Identifiability

There are I + 1 parameters in the model, but only I

population means to estimate.

We need some identifiability constraint for the αi .

Reference group constraint: set αI = 0, so that µ is the

mean of the Ith population (referred to as baseline
population) and αi is the difference from the baseline
population.

13/86
 Hypothesis testing

Hypotheses
H 0 : µ1 = · · · = µI
H1 : not all the µi are equal

Main idea: compare the variability between groups with

that within groups (also known as error variability).

If the between-group variability is much greater than the

within-group variability, conclude that the µi vary → reject
H0 in favour of H1 .

14/86
 Analysis of variance

Notation
I
P
Total number of observations: N = ni
i=1
ni
1
P
ith sample mean: ȳi· = ni
yij
j=1
ni
I P
1
P
Overall mean: ȳ = N
yij
i=1 j=1

15/86
 Analysis of variance

Sums of squares
ni
I P
(yij − ȳ)2 measures the total variability in Y.
P
SST =
i=1 j=1
This can be rewritten as
X ni
I X ni
I X
X
2
SST = (ȳi· − ȳ) + (yij − ȳi· )2
i=1 j=1 i=1 j=1
I
X ni
I X
X
= ni (ȳi· − ȳ)2 + (yij − ȳi· )2 = SSB + SSE
i=1 i=1 j=1

SSB is between group and SSE is within group error.

16/86
 Analysis of variance

ANOVA table

Source SS df MS F

SSB MSB
Between-group SSB I−1 MSB =
I−1 MSE
SSE
Within-group (Error) SSE N−I MSE =
N−I

Total SST N−1

17/86
 Sources of variability

Between group variability

I
1
ni (ȳi· − ȳ)2
P
MSB = I−1
i=1

Within group variability

I P nj I
1 1
(yij − ȳi· )2 = (ni − 1)s2i ,
P P
MSE = N−I N−I
i=i j=1 i=1

where s2i is the sample variance of the ith group.

MSB is analogous to sample variance of group means.

MSE is a weighted average of sample variances.

18/86
 F-test

We use an F-test to test our hypothesis H0 : µ1 = · · · = µI .

This is based on a ratio of two variances.

Test Statistic
MSB
F= MSE

Numerator
I
E(MSB) = σ 2 + ni (µi − µ̄)2 /(I − 1),
P
i=1

Denominator
E(MSE) = σ 2

19/86
 F-test

Under H0 , E(MSB) = σ 2 and so we would expect F to take

values close to 1 if H0 is true.

Otherwise E(MSB) > σ 2 and the F-statistic will tend to be

large.

Under H0 the F-statistic follows the F(I − 1, N − I)

distribution. We reject H0 if F > F1−α (I − 1, N − I) for a
test significance level of 100α%.

Note that rejecting H0 tells us that not all the µi are equal,
but it does not specify which of the µi values differ.

20/86
Medicine Dosage for Ill Rabbits

21/86
 Example

Example adapted from Applied Longitudinal Analysis by

G. Fitzmaurice, N. Laird and J. Ware.

Evaluating the differences in drug dosage between four

different breeds of rabbits.

The dosage which is required to treat each breed is

recorded.

10 rabbits from each breed.

Goal: assess any differences in dosages across the four

breeds.

22/86
 Example

The sample means for the four breeds were:

ȳ1 = 25.9, ȳ2 = 22.2, ȳ3 = 20.0, ȳ4 = 19.6.

ANOVA table
Source SS df MS F
Between-group (Breed) 249.9 3 83.3 8.5
Within-group (Error) 350.9 36 9.7
Total 600.8 39

F = 83.3/9.7 = 8.5 (p < 0.001)

Estimated (common) standard
√ deviation of the distribution
of scores in each group = 9.7 = 3.11.

23/86
 Estimation

If we are interested in estimating a single mean, µi then we

can compute a point estimate and CI as follows.

Point Estimate
µ̂i = ȳi

Variance
σ2
Var (ȳi ) = ni
, where σ̂ 2 = MSE.

Confidence Interval
ȳi ± t0.975 (N − I) √σ̂ni .

24/86
 Estimating linear combinations

We can compute estimates for linear combinations of the

PI
form ci µi in a similar manner.
i=1

Point Estimate
I
P I
P
ci µ̂i = ci ȳi
i=1 i=1

Variance
 I
  I

c2i /ni σ2.
P P
Var ci ȳi =
i=1 i=1

25/86
 Estimating linear combinations - Example

Suppose in the rabbit example we are interested in

estimating the difference between breeds 1 and 2.

We can compute the difference between the two breed

sample means.

Difference in Means
µ̂1 − µ̂2 = ȳ1 − ȳ2 = 25.9 − 22.2 = 3.7

Variance
Var (µ̂1 − µ̂2 ) = ( 101 + 1
10
)σ̂ 2 = 0.2 × 9.7 = 1.9.

26/86
 Multiple comparisons

When the F-test is significant, as is the case in the rabbit

example, we wish to identify the pairs that differ.

We could attempt to do this by checking if the 95%

confidence interval for each contrast between means
contains zero.

However, this method would result in a high rate of false

positive results because of multiple comparisons.

27/86
 Multiple comparison procedures

Bonferroni adjustment

For Q contrasts tested, divide the significance level by Q.

For example, there are Q = 6 pairwise comparisons for

four groups, and each can be tested at level 0.05/6.

This method is good if we are interested in just a few

contrasts, but it becomes very conservative when Q is
large.

28/86
 Multiple comparison procedures

Scheffé’s method

Simultaneous confidence intervals for all possible pairs:

v
I u I
X u X c2i
ci ȳi ± δ tσ̂ 2
i=1 i=1
ni

where δ 2 = (I − 1)F1−α (I − 1, N − I) is a larger critical value.

For pairwise comparisons this reduces to

s  
2
1 1
ȳi − ȳj ± δ σ̂ + .
ni nj

29/86
 Multiple comparison procedures

Example: rabbit data

√
δ = 3 × 2.886 = 2.94 and
s   s  
2
1 1 1 1
δ σ̂ + = 2.94 9.7 + = 4.1.
ni nj 10 10

Hence, for instance, the 95% confidence interval for

µ1 − µ4 is

25.9 − 19.6 ± 4.1 = (2.2, 10.4).

In general, differences in means of more than 4.1 are

significant.

30/86
 Multiple comparison procedures

Differences in Means

For Scheffé’s method the F-test is significant if and only if

at least one contrast is significantly different from zero.
Q: Which means differ significantly from each other?

Graphical representation

19.6 20.0 22.2 25.9

µ̂4 µ̂3 µ̂2 µ̂1

31/86
 Two-way ANOVA

One-way ANOVA

One factor with I different levels

Compare mean response between factor levels

Two-way ANOVA

Two factors
Compare mean response across the levels of factor 1 and
factor 2.

32/86
 Two-way ANOVA

Questions of interest

1 Does the mean outcome differ between levels of factor 1?

2 Does the mean outcome differ between levels of factor 2?

3 Is there an interaction between factors 1 and 2? Do

differences between the levels of factor 2 depend on the
levels of factor 1 (or vice versa)?

33/86
 Linear models
ANOVA as a linear model

ANOVA is a normal linear model and can be represented

as a regression model with dummy variables.
Consider a factor with I levels:
Define X1 = 1 if the subject belongs to level 1, and 0
otherwise;
Define X2 , . . . , XI−1 similarly.

In this parameterisation the last level of the factor is

selected as a “reference”.
For example, for a subject at level 2:
X1 X2 X3 . . . XI−1
0 1 0 ... 0

34/86
 Linear models

Model formulation
Thus the ANOVA model Yij = µi + eij can be expressed as

Yij = β0 + β1 X1ij + β2 X2ij + · · · + βI−1 XI−1;ij + eij

with

µ1 =β0 + β1
µ2 =β0 + β2
..
.
µI =β0

35/86
 Linear models

The normal linear model representation of ANOVA is more

attractive because:
it can handle balanced (i.e. equal observations in each
group) and unbalanced data in a seamless fashion;
in addition to the usual ANOVA table summaries, it
provides other useful and interpretable results such as
estimates of effects and standard errors;
generalizations of ANOVA to include continuous
explanatory variables (and interactions among categorical
and continuous explanatory variables) are straightforward.

36/86
 Fixed or random?

Fixed effects

Only finitely many levels of interest

Data collected from every level
Conclusions from modelling fixed effects apply only to
those levels

Random effects

Large population of levels

Data collected from a randomly chosen set of levels
Conclusions apply to the whole population

37/86
 Examples
Factory machine output

We measure the output of a machine that makes iPad

components. The machine is operated by a number of
workers.

Is the worker effect a fixed or random effect?

Drug study

A researcher wants to compare the effect of three drugs (A,

B, and C). She is only interested in the comparison of these
three drugs.

Is the drug effect a fixed or random effect?

38/86
 Examples

Drug study

Let’s extend the previous drug study example.

Suppose that four clinics are randomly selected from a

population of clinics in a region.

The researcher wants to make inference for the drug effects

across all clinics, not only the ones included in the study.

Is the clinic effect:

A: fixed effect
B: random effect

39/86
Random, fixed, and mixed effects models

A model in which all effects are fixed (e.g. original drug

example) is called a fixed effects model.

A model with only random effects (e.g. factory machine

output example) is called a random effects model.

A model with both fixed and random effects (e.g. drug

trials at different clinics example) is called a mixed effects
model.

40/86
 One-factor random effects model

Consider the iPad factory example.

The workers in the study can be considered as being a

sample from a larger population of workers.

The worker effect is the only variable, and is modelled as

random.

This can therefore be described as a one-factor random

effects model.

41/86
 Example - Factory machine output
Data

Let yij be the output measurement for worker i on day j.

i = 1, . . . , I j = 1, . . . , ni .

Model

We can consider the random effects model:

yij = µ + ai + eij

Here µ is the overall mean, ai is the random effect for

worker i and eij is residual error.
We assume the ai are i.i.d. N(0, σA2 ) independently of eij
which are i.i.d. N(0, σE2 ).
42/86
 Example - Factory machine output

Question of interest

Is there variation in output in this population?

If there is no variation in output, then there will be no
worker-specific error.
Test HA : σA2 = 0.
Note: in a fixed effects model (one-way ANOVA),
I
1
σA2 = I−1 αi2 .
P
i=1

43/86
 Variance components

Under this model Var (Yij ) = σA2 + σE2 (the Yij share a
common variance).
However, the Yij are no longer independent.


Cov Yi1 j , Yi2 j = 0 (different workers are independent of
each other).
Within the same level the Yij are not independent of each
other:

Cov Yij1 , Yij2 = Cov ai + eij1 , ai + eij2 = σA2

(measurements from the same worker are correlated).

44/86
 Intraclass correlation

The correlation between observations at the same level is

known as the intraclass correlation coefficient (ICC):

Intraclass correlation coefficient

Cov Yij1 , Yij2 σ2
ρ= p = 2 A 2.
Var (Yij1 ) Var (Yij2 ) σA + σE

When the variation between levels (workers) is much

larger than the variation within the levels, ρ approaches 1.
When the variation between levels (workers) is much
smaller than the variation within the levels, ρ approaches 0.

45/86
 Repeatability and Reproducibility

Repeatability refers to the variation in measurements

made on the same subject (machine) under identical
circumstances (same worker).

This is equivalent to the ICC expressed as a percentage.

Reproducibility refers to the variation in measurements

made on the same subject (machine) under changing
circumstances (different workers).

This can be computed as 100% − repeatability.

46/86
 Estimation
Assume for simplicity that we have an equal number of
observations, J, for each level of factor A.
Decomposing the variation

SST = SSA + SSE

I X
X J
SST = (yij − ȳ)2 ,
i=1 j=1

I
X
SSA = J (ȳi· − ȳ)2 ,
i=1
I X
X J
SSE = (yij − ȳi· )2 .
i=1 j=1

47/86
 Estimation

ANOVA table

Source SS df MS F

SSA MSA
A SSA I−1 MSA =
I−1 MSE
SSE
Error SSE I(J − 1) MSE =
I(J − 1)

Total SST IJ − 1

48/86
 Distributional results

We can derive distributional results for SSE and SSA.

Recall that ȳi· = µ + ai + ei· and ȳ = µ + a· + e··

Distribution of SSE

I X
X J
SSE = (yij − ȳi· )2
i=1 j=1
I X
X J
= (eij − ei· )2 ∼ σE2 χ2 (I[J − 1])
i=1 j=1

49/86
 Distributional results

Distribution of SSA

I
X
SSA = J (ȳi· − ȳ)2
i=1
XI
=J (ai − a· + ei· − e·· )2
i=1
σE2
X  
2 2
=J (gi − ḡ) ∼ J σA + χ2 (I − 1)
J

where gi = ai + ei· with Var (gi ) = σA2 + σE2 /J.

50/86
 Distributional results

Note that SSA is independent of SSE, since

eij is independent of ai for all j
eij − ei· is independent of ei· .

Because the distribution of eij is normal, the latter can by

proved by showing that the covariance of the two terms is
zero.

51/86
 Distributional results

Proof of independence

P P 
j eij j eij



Cov eij − ei· , ei· = Cov eij − ,
J J

Var (eij )
= − Var (ei· )
J

σE2 σ2
= − E = 0.
J J

52/86
 Hypothesis Test

Test the null hypothesis HA : σA2 = 0 using the F-statistic.

F-statistic
MSA
F=
MSE
SSA/(I − 1)
=
SSE/(N − I)
∼ F(I − 1, N − I)

where N = IJ is the total number of observations.

The computational procedure here is the same as that of a

fixed-effect model (one-way layout).

53/86
 Estimation of variance components

Method of moments

Since

E(MSE) = σE2
E(MSA) = σE2 + JσA2 ,

we can take

σ̂E2 = MSE and σ̂A2 = (MSA − MSE)/J

which are unbiased estimates of σE2 and σA2 .

These are sometimes also called ANOVA estimators.

54/86
 Method of moments estimators

Ease of calculation by hand was an advantage in the

pre-computing days.

Potential problem: these could return a negative value for

σ̂A2 if MSE > MSA.

Also, while for a balanced design (equal observations per

cell) the ANOVA decomposition into sums of squares is
unique, this is not true for unbalanced data.

55/86
 Example
Nitrogen concentrations in Mississippi river
We have data about 37 nitrogen concentrations at six randomly
selected influents of the Mississippi river. This is an unbalanced
dataset as the number of observations is not the same for all
levels of the factor influent.
This dataset comes from Littel, R. C., Milliken, G. A., Stroup, W. W., and Wolfinger, R. D.
(1996), SAS System for Mixed Models, SAS Institute (Data Set 4.2).

56/86
 Nitrogen data

40
Nitrogen concentration

30
20
10

1 2 3 4 5 6

Influent
57/86
 Example

Random effects model

Fit a model for nitrogen concentration with influent as a

random factor.
Test the hypothesis HA : σA2 = 0, i.e. whether there is
significant variation between influents.

58/86
 Example

ANOVA table
Source SS DF MS F p
Influent 1925.19 5 385.04 9.04 < 0.0001
Error 1319.78 31
Total 3244.97 36

F-test

The F-statistic is significant, indicating that there is

variation in nitrogen concentrations between influents.

59/86
Variance component estimates

Method of moments estimator of σE2 :

1319.779
σ̂E2 = = 42.57.
31
Unbalanced data → E(MSA) = σE2 + J̃σA2 where
PI
n− i=1 n2i /n
J̃ = = 6.0973.
I−1
Method of moments estimator of σA2 :

1925.194/5 − σ̂E2
σ̂A2 = = 56.17.
6.0793
60/86
 Two-factor random effects model

Consider a model with two random effects of the form

yijk = µ + ai + bj + (ab)ij + eijk

for i = 1, . . . , I, j = 1, . . . , J and k = 1, . . . , K.

Here, {ai }, {bj }, {(ab)ij } and {eijk } are independent

normal random variables with mean zero and variances σA2 ,
σB2 , σAB
2
and σE2 respectively.

61/86
 Sums of squares - Factor A

Factor A
I
X
SSA = JK (ȳi·· − ȳ)2 .
i=1

We have
ȳi·· = µ + ai + b· + (ab)i· + ei··
ȳ = µ + a· + b· + (ab)·· + e···
so the sum of squares can be expressed as
I
X
SSA = JK (gi − ḡ)2 where gi = ai + (ab)i· + ei·· .
i=1

62/86
 Sums of squares

Factor A

2
σAB σ2
SSA ∼ JKVar (gi ) χ2 (I − 1) with Var (gi ) = σA2 + + E.
J JK
This gives
2
σE2
   2 
2 σAB χ (I − 1)
E(MSA) = JK σA + + E
J JK I−1
2 2
 
σ σ
= JK σA2 + AB + E .
J JK

63/86
 Sums of squares

Factor B
Similarly

σ2 σ2
 
SSB ∼ IK σB2 + AB + E χ2 (J − 1)
I IK

and
2
σE2
   2 
σAB χ (J − 1)
E(MSB) = IK +σB2 + E
I IK J−1
2 2
 
σ σ
= IK σB2 + AB + E .
I IK

64/86
 Sums of squares

Error

I X
X J X
K
SSE = (yijk − ȳij· )2
i=1 j=1 k=1

I X
X J X
K
= (eijk − eij· )2
i=1 j=1 k=1

∼ σE2 χ2 (IJ(K − 1)).

65/86
 Mean sums of squares

The mean squares MSA and MSE are independent because

ai + (ab)i· + ei·· are independent of eijk − eij· .
However, MSA/MSE does not follow an F-distribution
under HA , because the F-test requires E(MSA) = E(MSE)
under the null hypothesis.
Here this is not true because of the extra term for the
interaction variance.
The statistic MSA/MSE doesn’t have an F-distribution
unless there is no interaction.

66/86
 Interaction term

The interaction sum of squares SSAB is equal to

I X
X J X
K
((ab)ij − (ab)i· − (ab)·j + (ab)·· + eijk − ei·· − e·j· + e··· )2
i=1 j=1 k=1

From this, we can derive that:

2
E(SSAB) = KσAB (I − 1)(J − 1) + σE2 (I − 1)(J − 1)

2
E(MSAB) = KσAB + σE2 .

67/86
 ANOVA table

ANOVA table
Source MS E(MS) F

SSA MSA
A MSA = σE2 + KσAB
2
+ JKσA2 FA =
I−1 MSAB
SSB MSB
B MSB = σE2 + KσAB
2
+ IKσB2 FB =
J−1 MSAB
SSAB MSAB
AB MSAB = σE2 + KσAB
2
FAB =
(I − 1)(J − 1) MSE

SSE
Error MSE = σE2
IJ(K − 1)

68/86
 Test statistics for hypothesis tests
Under HA : σA2 = 0

MSA
FA = ∼ F(I − 1, (I − 1)(J − 1))
MSAB

Under HB : σB2 = 0

MSB
FB = ∼ F(J − 1, (I − 1)(J − 1))
MSAB

2
Under HAB : σAB =0

MSAB
FAB = ∼ F((I − 1)(J − 1), IJ(K − 1))
MSE
69/86
 Example
Assessment of heart X-rays: observer error
Two randomly selected radiologists at the Western
Infirmary in Glasgow calculated the diagnostic ratio (ratio
of width of the heart to width of the thorax).
Each radiologist assessed the heart X-rays of 15 patients
on two separate occasions.
X-rays were presented to the observers in random order.
The purpose of the study was to investigate the variability
in the assessments made by the different observers.

70/86
 Example

Data Source
The data are taken from J. Aitchison, J.W. Kay and I.J. Lauder’s
text Statistical Concepts and Applications in Clinical Medicine,
2004.

Question of interest
Would it be safe to make the assessments using different
observers, or should a single observer be used to make all
assessments of this type?

71/86
 Interaction plot of heart X-ray data

1
0.55

2 2 Observer
2
2 2
2
2
1 1
1
0.50

2
2 1 2 2 1
2
1 1
mean of DR

1 2
1 2 1
1 1
2
2
0.45

1 1
1
2
1
2 1
1
2 2
0.40

2 1 1
2
1 2
1
1 2
1
1

2
1
1
0.35

1 2 3 4 5 6 7 8 9 10 12 14

Case
72/86
 Random effects model

Fit a model for the diagnostic ratio with case and observer
included as random factors.

Each level of the factor case appears with every level of the
observer factor. Such factors are known as crossed effects.

Goal: to assess the statistical significance of the variability

in measurements of the diagnostic ratio due to cases,
observers and the interaction between cases and observers.

73/86
 ANOVA table

ANOVA table
Source DF SS MS F p
Case 14 0.1569397 0.011210 33.19 < 0.0001

Observer 1 0.0024448 0.002448 7.24 < 0.0001

Case:Observer 14 0.0047289 0.000338 2.97 0.006

Error 30 0.0034075 0.000114

74/86
 Variance component estimation
Method of moments (unbiased estimates)

2
MSA − K σ̂AB − σ̂E2 MSA − MSAB
σ̂A2 = =
JK JK

MSB − MSAB
σ̂B2 =
IK

2 MSAB − MSE
σ̂AB =
K

σ̂E2 = MSE

75/86
 Variance component estimation

Variance component estimates

Variance component Estimate
Case 0.002718

Observer 0.000070

Case:Observer 0.000112

Error 0.000114

76/86
 Diagnostic plots

Residuals vs Fitted Normal Q−Q

0.03

4
13 ● 13 ●

Standardized residuals
● ●

2
0.01

● ● ●●
●●
Residuals

● ● ● ●
●●● ● ● ●●●●●
● ●● ● ●
● ● ●
● ●
●●●●
●● ●●●
● ●●● ● ● ●●
●
●
●●
●●
●
●
● ●● ●
●●
●●
●●
●

0
● ● ●
● ●
●●
●
●●● ● ●
●
●●
●●
●●
●
●●
●
−0.01

● ● ● ●● ●●
●●
● ●●●●
●
●●
● 30

−2
● 30
−0.03

14 ●
● 14

0.40 0.45 0.50 0.55 −2 −1 0 1 2

Fitted values Theoretical Quantiles

77/86
 Nested designs

A nested design is a design in which levels of factor B

vary only within the levels of another factor A.

Such a design is only possible for models which include

random effects.

Non-nested effects are called crossed effects.

An example of crossed effects would be the case and the

observer in in the heart X-ray example.

78/86
 Example

Silicon wafers

An experiment is conducted at a production line in which

y, the thickness of the oxide layer of silicon wafers
produced, is measured by first selecting 8 lots at random,
then from each lot randomly selecting 3 wafers, and for
each wafer randomly selecting 3 sides which are to be
measured.

Questions

Which type of design is appropriate here?

How many factors do we have?

79/86
 Example

Silicon wafers

We have a three-factor nested design.

Factor A is the lot, with 8 lots selected from the population.
Factor B is the wafer, with 3 wafers selected within each
lot.
Factor C is the side, with 3 sides selected for each wafer.
B is nested in A and C nested in A and B.

80/86
 Example

Model

The model can be written as

yijk = µ + ai + bj(i) + ck(ij)

for i = 1, . . . , 8, j = 1, 2, 3 and k = 1, 2, 3.
In the above notation, brackets denote nesting.
There is no random noise term in the model because no
replicate measurements were taken.
Consider the general model with i = 1, . . . , I, j = 1, . . . , J
and k = 1, . . . , K.

81/86
 Sums of squares
Factor A

I
X
SSA = JK (yi·· − ȳ)2
i=1

where
yi·· = µ + ai + b·(i) + c·(i·)
and
ȳ = µ + a· + b·(·) + c·(··) ,
so
I
X
SSA = JK (ai· − a·· + b·(i) − b·(·) + c·(i·) − c·(··) )2 .
i=1

82/86
 Sums of squares

Expected mean squares

Assuming independent normal distributions for the random
effects with variances σA2 , σB(A)
2 2
and σC(AB) = σE2 respectively,
we obtain
" 2 2
#
σB(A) σC(AB)
E(SSA) = JK (I − 1)σA2 + (I − 1) + (I − 1)
J JK

and !
2 2
σB(A) σC(AB)
E(MSA) = JK σA2 + +
J JK
with degrees of freedom equal to I − 1.

83/86
 Sums of squares

Factor B

The sum of squares for factor B (which is nested in A) is

I X
X J
SSB(A) = K (yij· − ȳi·· )2
i=1 j=1

E[MSB(A)] = σE2 + KσB(A)

2

Degrees of freedom = I(J − 1) (for different i the j’s are

different).

84/86
 Sums of squares

Factor C

The sum of squares for factor C (which is nested in A and

B) is
I X
X J
SSC(B) = SSE = K (yijk − ȳij· )2
i=1 j=1

E(MSE) = σE2
Degrees of freedom = IJ(K − 1)
We can think of this term as a measure of the variation
among k’s given i and j.

85/86
 Hypothesis testing

To test for A, divide MSA by MSB(A).

To test for B(A), divide MSB(A) by MSE.
For the silicon wafer data, these tests are shown in the
following ANOVA table:

ANOVA table
Source DF SS MS F p
A 7 9025.319 1289.331 10.73 < 0.0001
B(A) 16 1922.667 120.167 9.56 < 0.0001
C(B) 48 603.333 12.569

86/86