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.


<!--chapter:end:507-NestingMixed.Rmd-->

# Randomized Complete Blocks




## Immer data

One early and
famous example of a Randomized Complete Block analyzed by Fisher
involves five varieties of barley grown at six locations (which
included Crookston, Waseca, ..., data from Minnesota!).  The original
experiment had two years of data; we will only look at the second year.
This data set is included in the MASS package.

<div class="borderboxc">

```{.r .RSource}
> library(cfcdae)
> library(MASS)
> library(lme4)
> library(car)
> immer
   Loc Var    Y1    Y2
1   UF   M  81.0  80.7
2   UF   S 105.4  82.3
3   UF   V 119.7  80.4
4   UF   T 109.7  87.2
5   UF   P  98.3  84.2
6    W   M 146.6 100.4
7    W   S 142.0 115.5
8    W   V 150.7 112.2
9    W   T 191.5 147.7
10   W   P 145.7 108.1
11   M   M  82.3 103.1
12   M   S  77.3 105.1
13   M   V  78.4 116.5
14   M   T 131.3 139.9
15   M   P  89.6 129.6
16   C   M 119.8  98.9
17   C   S 121.4  61.9
18   C   V 124.0  96.2
19   C   T 140.8 125.5
20   C   P 124.8  75.7
21  GR   M  98.9  66.4
22  GR   S  89.0  49.9
23  GR   V  69.1  96.7
24  GR   T  89.3  61.9
25  GR   P 104.1  80.3
26   D   M  86.9  67.7
27   D   S  77.1  66.7
28   D   V  78.9  67.4
29   D   T 101.8  91.8
30   D   P  96.0  94.1

Here are the results separately by block. There are large block to block differences (and all of the treatment to treatment differences are within each boxplot), so blocking should be a real help.

> boxplot(Y2~Loc,data=immer,main="Barley by location")
Boxplot of barley by location

Figure 39.14: Boxplot of barley by location

Here are the results separately by treatment. We can maybe see some treatment differences, but we also see a lot of variation within each treatment. Much of this is block to block variability that would show up in our error if we hadn’t blocked.

> boxplot(Y2~Var,data=immer,main="Barley by variety")
Boxplot of barley by variety

Figure 39.15: Boxplot of barley by variety

39.3.1 As Fixed Effects

Fit as fixed blocks first. Get into the habit of treatments after blocks (which is the order you’ll need if you have missing data).

> barley.lm <- lm(Y2~Loc+Var,data=immer)

Residuals aren’t too bad, although there is just a hint of a dip in the middle.

> plot(barley.lm,which=1)
Residual plot of barley data

Figure 39.16: Residual plot of barley data

Normality looks good.

> plot(barley.lm,which=2)
NPP of barley data residuals

Figure 39.17: NPP of barley data residuals

Here is the ANOVA, we just enter treatments after blocks. Variety effects are moderately significant. We do not test the block effects (even though R prints an F, it doesn’t know about blocking factors).

> anova(barley.lm)
Analysis of Variance Table

Response: Y2
          Df  Sum Sq Mean Sq F value    Pr(>F)    
Loc        5 10285.0 2056.99 10.3901 5.049e-05 ***
Var        4  2845.2  711.29  3.5928   0.02306 *  
Residuals 20  3959.5  197.98                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

You can do all the usual things: summaries, pairwise comparisons, marginal means, etc.

Varieties M and S have low yields, variety T has a high yield, and the others are in the middle. Again, ignore ``testing” related to block effects.

> summary(barley.lm)

Call:
lm(formula = Y2 ~ Loc + Var, data = immer)

Residuals:
    Min      1Q  Median      3Q     Max 
-25.007  -8.665  -0.900   8.185  23.893 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   93.133      2.569  36.254  < 2e-16 ***
Loc1          -1.493      5.744  -0.260 0.797543    
Loc2         -15.593      5.744  -2.715 0.013344 *  
Loc3         -22.093      5.744  -3.846 0.001008 ** 
Loc4          25.707      5.744   4.475 0.000232 ***
Loc5         -10.173      5.744  -1.771 0.091790 .  
Var1          -6.933      5.138  -1.349 0.192261    
Var2           2.200      5.138   0.428 0.673081    
Var3         -12.900      5.138  -2.511 0.020748 *  
Var4          15.867      5.138   3.088 0.005797 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.07 on 20 degrees of freedom
Multiple R-squared:  0.7683,    Adjusted R-squared:  0.664 
F-statistic: 7.369 on 9 and 20 DF,  p-value: 0.0001069
> model.effects(barley.lm,"Var")
         M          P          S          T          V 
 -6.933333   2.200000 -12.900000  15.866667   1.766667 

Using HSD, only S and T seem to differ, the others cannot be distinguished.

> pairwise(barley.lm,Var)
                                    
Pairwise comparisons ( hsd ) of Var 
           estimate signif diff      lower     upper
  M - P  -9.1333333    24.30866 -33.441989 15.175322
  M - S   5.9666667    24.30866 -18.341989 30.275322
  M - T -22.8000000    24.30866 -47.108656  1.508656
  M - V  -8.7000000    24.30866 -33.008656 15.608656
  P - S  15.1000000    24.30866  -9.208656 39.408656
  P - T -13.6666667    24.30866 -37.975322 10.641989
  P - V   0.4333333    24.30866 -23.875322 24.741989
* S - T -28.7666667    24.30866 -53.075322 -4.458011
  S - V -14.6666667    24.30866 -38.975322  9.641989
  T - V  14.1000000    24.30866 -10.208656 38.408656

“Old school” estimate of block variation based on barley.lm (will agree with REML).

> (2056.99-197.98)/5
[1] 371.802

39.3.2 As Random Effects

Maybe you think the blocks are supposed to be representative of all of Minnesota; then we would treat them as random.

> barley.lmer <- lmer(Y2~(1|Loc)+Var,data=immer)

For this simple model and balanced data, we get all the same results:

> summary(barley.lmer)
Linear mixed model fit by REML ['lmerMod']
Formula: Y2 ~ (1 | Loc) + Var
   Data: immer

REML criterion at convergence: 227

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.92838 -0.49474 -0.05208  0.72439  1.54701 

Random effects:
 Groups   Name        Variance Std.Dev.
 Loc      (Intercept) 371.8    19.28   
 Residual             198.0    14.07   
Number of obs: 30, groups:  Loc, 6

Fixed effects:
            Estimate Std. Error t value
(Intercept)   93.133      8.280  11.247
Var1          -6.933      5.138  -1.349
Var2           2.200      5.138   0.428
Var3         -12.900      5.138  -2.511
Var4          15.867      5.138   3.088

Correlation of Fixed Effects:
     (Intr) Var1   Var2   Var3  
Var1  0.000                     
Var2  0.000 -0.250              
Var3  0.000 -0.250 -0.250       
Var4  0.000 -0.250 -0.250 -0.250
> Anova(barley.lmer,test="F")
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: Y2
         F Df Df.res  Pr(>F)  
Var 3.5928  4     20 0.02306 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> pairwise(barley.lmer,Var)
                                    
Pairwise comparisons ( hsd ) of Var 
           estimate signif diff      lower     upper
  M - P  -9.1333333    24.30866 -33.441989 15.175322
  M - S   5.9666667    24.30866 -18.341989 30.275322
  M - T -22.8000000    24.30866 -47.108656  1.508656
  M - V  -8.7000000    24.30866 -33.008656 15.608656
  P - S  15.1000000    24.30866  -9.208656 39.408656
  P - T -13.6666667    24.30866 -37.975322 10.641989
  P - V   0.4333333    24.30866 -23.875322 24.741989
* S - T -28.7666667    24.30866 -53.075322 -4.458011
  S - V -14.6666667    24.30866 -38.975322  9.641989
  T - V  14.1000000    24.30866 -10.208656 38.408656

39.3.3 Relative Efficiency

Recall the fixed effects ANOVA:

> anova(barley.lm)
Analysis of Variance Table

Response: Y2
          Df  Sum Sq Mean Sq F value    Pr(>F)    
Loc        5 10285.0 2056.99 10.3901 5.049e-05 ***
Var        4  2845.2  711.29  3.5928   0.02306 *  
Residuals 20  3959.5  197.98                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This is our estimate of what the error variance would be in a CRD from the same kind of data we used in our RCB.

Noting that 518 is a lot more than 198, our blocking has been successful.

> (5*2056.99+24*197.98)/29
[1] 518.499

Here is the estimate of relative efficiency. Mostly it is a ratio of the MSE in the CRD (estimated) to the MSE in the RCB. In addition, there is a degrees of freedom correction term that adjusts for the loss of degrees of freedom in the RCB. This term is usually pretty close to 1. Here it is .98.

> (21*28/23/26) * 518.499/197.98
[1] 2.575151

The RE says that we would have needed a sample size 2.6 times as large to have the same power in a CRD as we achieved with blocking. That is very worthwhile.