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.
<!--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")

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")

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)

Figure 39.16: Residual plot of barley data
Normality looks good.
> plot(barley.lm,which=2)

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.