Chapter 39 Nested and Mixed Models
39.1 Particle count data
We have triplicate particle counts on each of three filters from each of two manufacturers for a total of 18 counts (data from Kuehl). Filter is nested in manufacturer.
> library(cfcdae)
> library(lme4)
> library(car)
> counts <- read.table("kuehl3.dat.txt",header=TRUE)
> counts <- within(counts, {manu <- factor(manu);filter<-factor(filter)})
> counts
manu filter partcount
1 1 1 1.12
2 1 1 1.10
3 1 1 1.12
4 1 2 0.16
5 1 2 0.11
6 1 2 0.26
7 1 3 0.15
8 1 3 0.12
9 1 3 0.12
10 2 1 0.91
11 2 1 0.83
12 2 1 0.95
13 2 2 0.66
14 2 2 0.83
15 2 2 0.61
16 2 3 2.17
17 2 3 1.52
18 2 3 1.58
Begin by assuming that we randomly selected the two manufacturers, so manufacturer is a random effect. we have manufacturer (random), filter nested in manufacturer (random), and residual.
There are
a couple of ways to set up the nested model. The
first is with slash notation: (1|manu/filter)
means do the random effects (here just the
intercept or constant) for manufacturer and for filter nested in manufacturer. In effect,
the (1|manu/filter)
expands to (1|manu) + (1|manu:filter)
.
The alternative is to just use (1|manu) + (1|manu:filter)
directly. These fit exactly the same model.
> partcount.lmer <- lmer(partcount ~ 1 + (1|manu/filter),data=counts)
> partcount.lmer.alt <- lmer(partcount ~ 1 + (1|manu) + (1|manu:filter),data=counts)
The estimated variances indicate that count to count variability on a filter is considerably smaller than variability between manufacturers, but filter to filter variability (within manufacturer) is larger still.
> summary(partcount.lmer)
Linear mixed model fit by REML ['lmerMod']
Formula: partcount ~ 1 + (1 | manu/filter)
Data: counts
REML criterion at convergence: 7.4
Scaled residuals:
Min 1Q Median 3Q Max
-1.34899 -0.39513 -0.07542 0.12340 2.73036
Random effects:
Groups Name Variance Std.Dev.
filter:manu (Intercept) 0.30331 0.5507
manu (Intercept) 0.10373 0.3221
Residual 0.02539 0.1593
Number of obs: 18, groups: filter:manu, 6; manu, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.7956 0.3222 2.469
Before going far into inference, we should check on assumptions. It looks like residual variability increases with mean.
> plot(partcount.lmer)

Figure 39.1: Residual plot for particle counts
And it looks like we have long tails.
> qqnorm(resid(partcount.lmer))

Figure 39.2: NPP of residuals for particle counts
There is no standard Box-Cox for linear mixed models, so we have a few alternatives:
- Use Box-Cox on a fixed effects analog model and hope that it is informative.
- When possible, use first principles to suggest a transformation.
- Just try some transformations and see if any of them help.
First principles suggests a square root (variance stabilizing transformation for the Poisson count distribution) or log (helpful if different counts have been rescaled by different factors). If we were just going to try something, we might try log first. And Box-Cox on the fixed effects analog suggests square root is probably a bit better than log, which is just on the edge of reasonable transformations.
> boxCox(lm(partcount~manu:filter,data=counts))

Figure 39.3: Box-Cox plot for particle counts
Refit with square roots and recheck residuals. Residual variance looks better (and if you try the log it goes too far).
> rpc.lmer <- lmer(sqrt(partcount)~1 + (1|manu/filter),data=counts)
> plot(rpc.lmer)

Figure 39.4: Residual plot for square root particle counts
Normality of residuals is better, but still not great.
> qqnorm(residuals(rpc.lmer))

Figure 39.5: NPP of residuals for square root particle counts
There are only six filters total, so it’s a bit of a lost cause to assess their normality.
> qqnorm(ranef(rpc.lmer)$"filter:manu"[[1]],main="Filter in Manufacturer effects")

Figure 39.6: NPP of filter effects for square root particle counts
On the other hand, I can guarantee that the NPP of the two manufacturer effects will look linear. :-)
Filter variance is the largest (at about 20 times residual), and manufacturer variance is in between (at about 10 times residual). You would think those would both be very significant, but we only have two manufacturers (one df between), so we have very little information about manufacturer.
> summary(rpc.lmer)
Linear mixed model fit by REML ['lmerMod']
Formula: sqrt(partcount) ~ 1 + (1 | manu/filter)
Data: counts
REML criterion at convergence: -16.5
Scaled residuals:
Min 1Q Median 3Q Max
-1.1919 -0.4497 -0.1293 0.2553 2.1725
Random effects:
Groups Name Variance Std.Dev.
filter:manu (Intercept) 0.105426 0.32469
manu (Intercept) 0.054346 0.23312
Residual 0.005305 0.07283
Number of obs: 18, groups: filter:manu, 6; manu, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.8219 0.2122 3.873
To
test terms using exactRLRT()
, we need the full model, the model with only
the term of interest, and the full model less the term of interest. Since there are
only two terms other than error in the model, the “only the term of interest” models
are also the “without the term of interest” models for the other term.
> library(RLRsim)
> rpc.manuonly.lmer <- lmer(sqrt(partcount) ~ 1 + (1|manu),data=counts)
> rpc.filteronly.lmer <- lmer(sqrt(partcount) ~ 1 + (1|manu:filter),data=counts)
Filter effects are highly significant.
> exactRLRT(rpc.filteronly.lmer,rpc.lmer,rpc.manuonly.lmer)
simulated finite sample distribution of RLRT.
(p-value based on 10000 simulated values)
data:
RLRT = 27.848, p-value < 2.2e-16
Manufacturer effects are not significant. Again, this is in large part because we only have two manufacturers (1 df between).
> exactRLRT(rpc.manuonly.lmer,rpc.lmer,rpc.filteronly.lmer)
simulated finite sample distribution of RLRT.
(p-value based on 10000 simulated values)
data:
RLRT = 0.40327, p-value = 0.1558
For completeness, here are the (often dubious) confidence intervals. Note that even for the highly significant filter effect, the upper endpoint of the CI is almost 3.8 times as large as the lower endpoint (almost 15 as variance instead of standard deviation). Variances are hard to estimate.
> confint(rpc.lmer,oldNames=FALSE)
Computing profile confidence intervals ...
2.5 % 97.5 %
sd_(Intercept)|filter:manu 0.18141409 0.6803578
sd_(Intercept)|manu 0.00000000 0.8682317
sigma 0.05116342 0.1155260
(Intercept) 0.30971341 1.3340794
39.1.1 Only two manufacturers?
What if there are only two manufacturers? In that case we are hardly taking a random sample, we are looking at all of them. We should now fit manufacturer as a fixed effect, but filters would still be random nested in manufacturer.
To include a fixed term in an lmer()
model, simply
don’t enclose it in parentheses.
> rpc.manfixed.lmer <- lmer(sqrt(partcount) ~ manu + (1|manu:filter),data=counts)
Looking at the summary, we have a fixed coefficient for the first manufacturer, and random effects only for filter and residuals. In this nicely balanced case the error and filter estimated variances are the same as in the fully random model.
> summary(rpc.manfixed.lmer)
Linear mixed model fit by REML ['lmerMod']
Formula: sqrt(partcount) ~ manu + (1 | manu:filter)
Data: counts
REML criterion at convergence: -16.2
Scaled residuals:
Min 1Q Median 3Q Max
-1.1729 -0.4688 -0.1103 0.2744 2.1534
Random effects:
Groups Name Variance Std.Dev.
manu:filter (Intercept) 0.105426 0.32469
Residual 0.005305 0.07283
Number of obs: 18, groups: manu:filter, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.8219 0.1337 6.149
manu1 -0.2122 0.1337 -1.588
Correlation of Fixed Effects:
(Intr)
manu1 0.000
Because filter is the only random effect other
than residual, the call to exactRLRT()
is simple.
As before, filter is highly significant.
> exactRLRT(rpc.manfixed.lmer)
simulated finite sample distribution of RLRT.
(p-value based on 10000 simulated values)
data:
RLRT = 27.848, p-value < 2.2e-16
You can do anova()
for an lmer
model with fixed
effects, and it
will give you an F test statistic for each fixed term
along with numerator df. Unfortunately, it doesn’t give
a denominator df and thus gives no p-value.
The issue is that lmer()
can fit crossed random
terms, and crossed random terms can make finding
denominator df very tricky. There are ways to
overcome that, but the authors of lme4
have
philosophical objections and have chosen not to
use them. Furthermore, lme4
won’t try to get
denominator df even in the absence of crossed
random terms.
Note: lme()
only fits
nested random terms,
so anova()
after lme()
gives a denominator df.
> anova(rpc.manfixed.lmer)
Analysis of Variance Table
npar Sum Sq Mean Sq F value
manu 1 0.013373 0.013373 2.5209
The car
package includes a function Anova()
,
which, when used with the test="F"
option, will produce an anova for fixed effects via
the Kenward-Roger approximation. This approximation is exact in many cases.
> Anova(rpc.manfixed.lmer,test="F")
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: sqrt(partcount)
F Df Df.res Pr(>F)
manu 2.5209 1 4 0.1875
KRmodcomp()
from package pbkrtest
compares two models that differ in their fixed effects via a Kenward-Roger test.
Here it
gives the same result as Anova did, but KRmodcomp()
can also compare models that differ
by more than one fixed term.
> library(pbkrtest)
> KRmodcomp(rpc.manfixed.lmer,rpc.filteronly.lmer)
large : sqrt(partcount) ~ manu + (1 | manu:filter)
small : sqrt(partcount) ~ 1 + (1 | manu:filter)
stat ndf ddf F.scaling p.value
Ftest 2.5209 1.0000 4.0000 1 0.1875
There is also a parametric bootstrap version that gives similar results, albeit much more slowly.
> library(pbkrtest)
> # not run, too slow
> # PBmodcomp(rpc.manfixed.lmer,rpc.filteronly.lmer)
39.1.2 Using lme()
We can fit the same models using lme()
instead of lmer()
.
> library(nlme)
> rpc.lme <- lme(sqrt(partcount)~1,random=~1|manu/filter,data=counts)
lme()
will nest random effects but not cross them. Thus
if we want just filter nested in manufacturer without
manufacturer itself, we must create a new factor that
enumerates all of the filters across manufacturers.
We can use the join()
function to do that.
> filterinmanu <- with(counts,conf.design::join(manu,filter))
Now fit with manufacturer fixed.
> rpc.manfixed.lme <- lme(sqrt(partcount)~manu,random=~1|filterinmanu,data=counts)
And presto, anova()
does what we want.
> anova(rpc.manfixed.lme)
numDF denDF F-value p-value
(Intercept) 1 12 37.81081 <.0001
manu 1 4 2.52093 0.1875
39.2 Soil porosity data
Data from problem 5-7 of Kuehl (1994 Duxbury). Fifteen fields are chosen at random. Two subsections are chosen at random from each field. Soil porosity is measured at a random location of each subsection; some subsections are measured at two locations.
Field and subsection are random, with subsection nested in field. The data set is not balanced.
Note that the variable section
enumerates all the sections,
but the variable sect is 1 or 2 within each field.
> soils <- read.table("kuehl5.dat.txt",header=TRUE)
> soils <- within(soils, {field <- factor(field);
+ section <- factor(section);sect <- factor(sect)})
> soils
field section sect porosity
1 1 1 1 3.846
2 1 1 1 3.712
3 1 2 2 5.629
4 1 2 2 2.021
5 2 3 1 5.087
6 2 4 2 4.621
7 3 5 1 4.411
8 3 6 2 3.357
9 4 7 1 3.991
10 4 8 2 5.766
11 5 9 1 5.677
12 5 10 2 3.333
13 6 11 1 4.355
14 6 11 1 6.292
15 6 12 2 4.940
16 6 12 2 4.810
17 7 13 1 2.983
18 7 14 2 4.396
19 8 15 1 5.603
20 8 16 2 3.683
21 9 17 1 5.942
22 9 18 2 5.014
23 10 19 1 5.143
24 10 20 2 4.061
25 11 21 1 3.835
26 11 21 1 2.964
27 11 22 2 4.584
28 11 22 2 4.398
29 12 23 1 4.193
30 12 24 2 4.125
31 13 25 1 3.074
32 13 26 2 3.483
33 14 27 1 3.867
34 14 28 2 4.212
35 15 29 1 6.247
36 15 30 2 4.730
Fit the model with sect nested within field. Fit the model with field and section. Since section enumerates all the sections individually, we don’t need to do “nesting” in the model when we use section. We’ll get the same results either way.
> porosity1 <- lmer(porosity~1+(1|field/sect),data=soils)
boundary (singular) fit: see help('isSingular')
> porosity2 <- lmer(porosity~1+(1|field) + (1|section),data=soils)
boundary (singular) fit: see help('isSingular')
There is some evidence of field to field variability, but the section within field is estimated at practically zero.
> summary(porosity1)
Linear mixed model fit by REML ['lmerMod']
Formula: porosity ~ 1 + (1 | field/sect)
Data: soils
REML criterion at convergence: 102.6
Scaled residuals:
Min 1Q Median 3Q Max
-2.33677 -0.53145 -0.04379 0.54268 1.80611
Random effects:
Groups Name Variance Std.Dev.
sect:field (Intercept) 3.967e-09 6.299e-05
field (Intercept) 5.947e-02 2.439e-01
Residual 9.360e-01 9.675e-01
Number of obs: 36, groups: sect:field, 30; field, 15
Fixed effects:
Estimate Std. Error t value
(Intercept) 4.4037 0.1741 25.29
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
We should check residuals before doing inference.
Normality looks OK.
> qqnorm(residuals(porosity1))

Figure 39.7: NPP of residuals for porosity data
> plot(porosity1)

Figure 39.8: Residual plot for porosity data
Remember how I told you that you could get some strange looking residual plots when you have random/mixed effects, especially if the residual variance is large compared to effects variance? Well, here is an example.
Random effects are “shrunk” back towards zero. This happens more when the variance of the random effect is small compared to the error variance. In this case, residual variance is much larger than field variance, and field variance and its random effects are being shrunk towards zero. This shows up in the residual plot as trend: negative random effects aren’t negative enough (at least by what we’re used to in fixed effects) and positive random effects aren’t positive enough.
Let’s test the field effect. It’s not significant.
> porosity.nofield <- lmer(porosity ~ 1 + (1|field:sect),data=soils)
boundary (singular) fit: see help('isSingular')
> porosity.onlyfield <- lmer(porosity ~ 1 + (1|field),data=soils)
> exactRLRT(porosity.onlyfield,porosity1,porosity.nofield)
simulated finite sample distribution of RLRT.
(p-value based on 10000 simulated values)
data:
RLRT = 0.11968, p-value = 0.3423
OK, now
suppose that we are only worried about these 15 fields, and we’re not
thinking of them as sampled from some larger population. In this case,
we would treat field as fixed. Note, we could also use (1|field:sect)
for
the random term.
Field is not significant as a fixed effect either.
> porosity.fieldfixed <- lmer(porosity~field+(1|section),data=soils)
boundary (singular) fit: see help('isSingular')
> Anova(porosity.fieldfixed,test="F")
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: porosity
F Df Df.res Pr(>F)
field 1.069 14 12.749 0.4558
> # or
> KRmodcomp(porosity.fieldfixed,porosity.nofield)
large : porosity ~ field + (1 | section)
small : porosity ~ 1 + (1 | field:sect)
stat ndf ddf F.scaling p.value
Ftest 1.069 14.000 12.749 1.0041 0.4558
How do the fixed field effects compare to the random field effects? The random effects are a lot smaller. We can see that either numerically or graphically, where the random effects are shrunk to practically 0.
> model.effects(porosity.fieldfixed,"field")
1 2 3 4 5 6 7 8
-0.62106667 0.43093333 -0.53906667 0.45543333 0.08193333 0.67618333 -0.73356667 0.21993333
9 10 11 12 13 14 15
1.05493333 0.17893333 -0.47781667 -0.26406667 -1.14456667 -0.38356667 1.06543333
> ranef(porosity1)$field[[1]]
[1] -0.12192887 0.05077098 -0.05859276 0.05353326 0.01142258 0.14095212 -0.08052188 0.02698154
[9] 0.12112455 0.02235895 -0.09290000 -0.02758758 -0.12686053 -0.04106074 0.12230839
> par(pty="s")
> plot(model.effects(porosity.fieldfixed,"field"),ranef(porosity1)$field[[1]],xlim=c(-1.1,1.1),ylim=c(-1.1,1.1),xlab="Fixed field effects",ylab="Random field effects")
> abline(0,1)

Figure 39.9: Random vs fixed field effects for porosity data
> par(pty="m")
With fields as fixed, the residual plot looks more like what you would expect.
> plot(porosity.fieldfixed,which=1)

Figure 39.10: Residual plot with fields fixed for porosity data
39.3 Serum glucose
These are the data from Table 7.10 of Kuehl. We are investigating how well an instrument measures serum glucose. We measure at three fixed levels of glucose. The machine may work differently on different days, so we choose three days at random. Also, we will do two runs or batches on the machine each day. That is, do everything once, and then do it all over again. There may be a run effect nested in day. Finally, each concentration is measured twice during each run.
Day and run random; run nested in day; day and run crossing concentration.
> glucose <- read.table("kuehl4.dat.txt",header=TRUE)
> glucose <- within(glucose,{conc <- factor(conc);
+ day <- factor(day);run <- factor(rn)})
> glucose
conc day rn y run
1 1 1 1 41.2 1
2 1 1 1 42.6 1
3 1 1 2 41.2 2
4 1 1 2 41.4 2
5 1 2 1 39.8 1
6 1 2 1 40.3 1
7 1 2 2 41.5 2
8 1 2 2 43.0 2
9 1 3 1 41.9 1
10 1 3 1 42.7 1
11 1 3 2 45.5 2
12 1 3 2 44.7 2
13 2 1 1 135.7 1
14 2 1 1 136.8 1
15 2 1 2 143.0 2
16 2 1 2 143.3 2
17 2 2 1 132.4 1
18 2 2 1 130.3 1
19 2 2 2 134.4 2
20 2 2 2 130.0 2
21 2 3 1 137.4 1
22 2 3 1 135.2 1
23 2 3 2 141.1 2
24 2 3 2 139.1 2
25 3 1 1 163.2 1
26 3 1 1 163.3 1
27 3 1 2 181.4 2
28 3 1 2 180.3 2
29 3 2 1 173.6 1
30 3 2 1 173.9 1
31 3 2 2 174.9 2
32 3 2 2 175.6 2
33 3 3 1 166.6 1
34 3 3 1 165.5 1
35 3 3 2 175.0 2
36 3 3 2 172.0 2
You always wondered where my Hasse diagrams came from, right?
> mixed.hasse(~Day+Conc+Day:Run+Day:Conc+Day:Conc:Run,c(3,3,2,2),random=c("Day","Run"),lwd=1,cex=1)

Figure 39.11: Hasse diagram for glucose experiment
Fit the model, and look at residuals. There appears to be increasing variance.
> glu.lmer <- lmer(y ~ conc + (1|day/run) + (1|conc:day) + (1|conc:day:run),data=glucose)
boundary (singular) fit: see help('isSingular')
> plot(glu.lmer)

Figure 39.12: Residual plot for glucose data
Try again on the square root scale. This is better.
> rglu.lmer <- lmer(sqrt(y) ~ conc + (1|day/run) + (1|conc:day) + (1|conc:day:run),data=glucose)
> plot(rglu.lmer)

Figure 39.13: Residual plot for sqrt glucose data
Looking at the fit summary we see:
- The concentration levels have large t-statistics. We still need to do the KR test, but these certainly look like they will be significant (and since we manipulated the concentration, we certainly hope that our instrument can find differences).
- Day effect variance is essentially 0.
- Run within day and concentration by day variances are about the same size as error variance. Maybe significant, maybe not.
- Concentration by run within day looks like it might be important.
> summary(rglu.lmer)
Linear mixed model fit by REML ['lmerMod']
Formula: sqrt(y) ~ conc + (1 | day/run) + (1 | conc:day) + (1 | conc:day:run)
Data: glucose
REML criterion at convergence: -39.4
Scaled residuals:
Min 1Q Median 3Q Max
-1.84022 -0.47269 -0.01539 0.58780 1.54751
Random effects:
Groups Name Variance Std.Dev.
conc:day:run (Intercept) 2.383e-02 1.544e-01
conc:day (Intercept) 5.048e-03 7.105e-02
run:day (Intercept) 9.300e-03 9.644e-02
day (Intercept) 3.309e-10 1.819e-05
Residual 3.190e-03 5.648e-02
Number of obs: 36, groups: conc:day:run, 18; conc:day, 9; run:day, 6; day, 3
Fixed effects:
Estimate Std. Error t value
(Intercept) 10.43086 0.05936 175.72
conc1 -3.93972 0.06283 -62.71
conc2 1.25351 0.06283 19.95
Correlation of Fixed Effects:
(Intr) conc1
conc1 0.000
conc2 0.000 -0.500
Test concentration; yes, it’s highly significant.
> Anova(rglu.lmer,test="F")
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: sqrt(y)
F Df Df.res Pr(>F)
conc 2052.7 2 4 9.475e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Test run; not significant.
> rglu.norun <-lmer(sqrt(y) ~ conc + (1|day) + (1|conc:day) + (1|conc:day:run),data=glucose)
boundary (singular) fit: see help('isSingular')
> rglu.runonly <- lmer(sqrt(y)~conc+(1|day:run),data=glucose)
> exactRLRT(rglu.runonly,rglu.lmer,rglu.norun)
simulated finite sample distribution of RLRT.
(p-value based on 10000 simulated values)
data:
RLRT = 0.60624, p-value = 0.18
Test concentration by run; highly significant.
> rglu.noconcrun <- lmer(sqrt(y) ~ conc + (1|day/run) + (1|conc:day),data=glucose)
boundary (singular) fit: see help('isSingular')
> rglu.concrunonly <- lmer(sqrt(y) ~ conc + (1|conc:day:run),data=glucose)
> exactRLRT(rglu.concrunonly, rglu.lmer,rglu.noconcrun)
simulated finite sample distribution of RLRT.
(p-value based on 10000 simulated values)
data:
RLRT = 19.644, p-value < 2.2e-16
Test concentration by day; not significant.
> rglu.noconcday <- lmer(sqrt(y)~conc+(1|day/run)+(1|conc:day:run),data=glucose)
boundary (singular) fit: see help('isSingular')
> rglu.concdayonly <- lmer(sqrt(y)~conc+(1|conc:day),data=glucose)
> exactRLRT(rglu.concdayonly,rglu.lmer,rglu.noconcday)
simulated finite sample distribution of RLRT.
(p-value based on 10000 simulated values)
data:
RLRT = 0.12869, p-value = 0.3148
What does it all mean? There are concentration differences, but we still need to compare what we are estimating with spiked concentrations that are in the samples to see if they are matching up.
There is no evidence that some days read high (or low), or that some runs read high (or low), or that some days read high at some concentrations and low at other concentrations. However, there is strong evidence that some runs read high at some concentrations and low at other concentrations. The size of this effect is an order of magnitude smaller than the concentration effects, but it is still somewhat worrisome.