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)
Residual plot for particle counts

Figure 39.1: Residual plot for particle counts

And it looks like we have long tails.

> qqnorm(resid(partcount.lmer))
NPP of residuals for particle counts

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:

  1. Use Box-Cox on a fixed effects analog model and hope that it is informative.
  2. When possible, use first principles to suggest a transformation.
  3. 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))
Box-Cox plot for particle 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)
Residual plot for square root particle counts

Figure 39.4: Residual plot for square root particle counts

Normality of residuals is better, but still not great.

> qqnorm(residuals(rpc.lmer))
NPP of residuals for square root particle counts

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")
NPP of filter effects for square root particle counts

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))
NPP of residuals for porosity data

Figure 39.7: NPP of residuals for porosity data

But oh my gosh! Look at the residuals versus fitted plot!
> plot(porosity1)
Residual plot for porosity data

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)
Random vs fixed field effects for porosity data

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)
Residual plot with fields fixed for porosity data

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)
Hasse diagram for glucose experiment

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)
Residual plot for glucose data

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)
Residual plot for sqrt glucose data

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.