[1] "Created: Wed Apr  1 16:46:56 2015"

See the introduction for an overview.

library(lme4)
library(ggplot2)
options(digits=5,show.signif.stars=FALSE)

Load in and plot the data:

data(penicillin, package="faraway")
summary(penicillin)
 treat    blend       yield   
 A:5   Blend1:4   Min.   :77  
 B:5   Blend2:4   1st Qu.:81  
 C:5   Blend3:4   Median :87  
 D:5   Blend4:4   Mean   :86  
       Blend5:4   3rd Qu.:89  
                  Max.   :97  
ggplot(penicillin,aes(x=blend,y=yield,group=treat,linetype=treat))+geom_line()

plot of chunk unnamed-chunk-3

ggplot(penicillin,aes(x=treat,y=yield,group=blend,linetype=blend))+geom_line()

plot of chunk unnamed-chunk-3

Fit the model

mmod <- lmer(yield ~ treat + (1|blend), penicillin)
summary(mmod)
Linear mixed model fit by REML ['lmerMod']
Formula: yield ~ treat + (1 | blend)
   Data: penicillin

REML criterion at convergence: 103.8

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.415 -0.502 -0.164  0.683  1.284 

Random effects:
 Groups   Name        Variance Std.Dev.
 blend    (Intercept) 11.8     3.43    
 Residual             18.8     4.34    
Number of obs: 20, groups:  blend, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)    84.00       2.48    33.9
treatB          1.00       2.75     0.4
treatC          5.00       2.75     1.8
treatD          2.00       2.75     0.7

Correlation of Fixed Effects:
       (Intr) treatB treatC
treatB -0.555              
treatC -0.555  0.500       
treatD -0.555  0.500  0.500
anova(mmod)
Analysis of Variance Table
      Df Sum Sq Mean Sq F value
treat  3     70    23.3    1.24

The standard output provides no inference although we can get confidence intervals:

confint(mmod)
               2.5 %  97.5 %
.sig01       0.00000  7.6767
.sigma       2.81917  5.8261
(Intercept) 79.28841 88.7116
treatB      -4.13671  6.1367
treatC      -0.13671 10.1367
treatD      -3.13671  7.1367

Inference with AOV

The aov() function can be used to fit simple models with a single random effects component. The results are reliable only for balanced data such as this.

lmod <- aov(yield ~ treat + Error(blend), penicillin)
summary(lmod)

Error: blend
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  4    264      66               

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
treat      3     70    23.3    1.24   0.34
Residuals 12    226    18.8               

We see that the test of the significance for the fixed effects which is effectively the same as the original F-test presented in ELM. Note that the p-values are provided only for the fixed effects terms. The fixed effect coefficients may be obtained as

coef(lmod)
(Intercept) :
(Intercept) 
         86 

blend :
numeric(0)

Within :
treatB treatC treatD 
     1      5      2 

pbkrtest package

We test the treatment effect using:

library(pbkrtest)
amod <- lmer(yield ~ treat + (1|blend), penicillin, REML=FALSE)
nmod <- lmer(yield ~ 1 + (1|blend), penicillin, REML=FALSE)
KRmodcomp(amod, nmod)
F-test with Kenward-Roger approximation; computing time: 0.12 sec.
large : yield ~ treat + (1 | blend)
small : yield ~ 1 + (1 | blend)
       stat   ndf   ddf F.scaling p.value
Ftest  1.24  3.00 12.00         1    0.34

The results are identical to the original output in the text and to the aov output because the data are balanced in this example.

The pbkrtest package also implements the parametric bootstrap. The idea is the same as explained in the text but the implementation in this package saves us the trouble of explicitly coding the procedure. Furthermore, the package makes it easy to take advantage of the multiple processing cores available on many computers. Given that the parametric bootstrap is expensive to compute, this is well worthwhile for very little additional effort for the user. We set up the computing clusters:

library(parallel)
nc <- detectCores()
clus <- makeCluster(rep("localhost", nc))

We execute the parametric boostrap:

pmod <- PBmodcomp(amod, nmod, cl=clus)
summary(pmod)
Parametric bootstrap test; time: 10.45 sec; samples: 1000 extremes: 349;
large : yield ~ treat + (1 | blend)
small : yield ~ 1 + (1 | blend)
         stat   df  ddf p.value
PBtest   4.05              0.35
Gamma    4.05              0.34
Bartlett 3.33 3.00         0.34
F        1.35 3.00 11.3    0.31
LRT      4.05 3.00         0.26

Several different ways of computing the p-value are shown. First, the proportion of bootstrap samples in which the bootstrapped test statistic exceeds the observed value is given. This results in a p-value of 0.36. Next three different ways of approximating the null distribution are presented, all giving similar but not identical results. The default number of boostrap samples is 1000. On today’s computers this will take only a few seconds. In this example, the estimated p-value is far from 0.05 so there is no doubt about the statistical significance. We certainly did not need more than 1000 samples and could have survived with less. However, for more complex models involving larger datasets, computational time may become a more important consideration. The approximations, if appropriate, will require fewer samples to work. This would motivate us to consider this approach. But in this example, we can afford to be profligate.

Finally, the outcome of LRT is reported which matches our earlier calculations.

RLRsim package

We can execute the test of the random effect blend in the penicillin example:

library(RLRsim)
pmod <- lmer(yield ~ treat + (1|blend), penicillin, REML=FALSE)
lmod <- lm(yield ~ treat, penicillin)
exactLRT(pmod, lmod)

    simulated finite sample distribution of LRT. (p-value based on
    10000 simulated values)

data:  
LRT = 3.4536, p-value = 0.0429

Again the results are much the same as in the text. There is a slight difference in that here we have used the ML estimate whereas we used REML version in the text.

lmerTest package

library(lmerTest)
mmod <- lmer(yield ~ treat + (1|blend), penicillin)
summary(mmod)
Linear mixed model fit by REML t-tests use Satterthwaite approximations
  to degrees of freedom [merModLmerTest]
Formula: yield ~ treat + (1 | blend)
   Data: penicillin

REML criterion at convergence: 103.8

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.415 -0.502 -0.164  0.683  1.284 

Random effects:
 Groups   Name        Variance Std.Dev.
 blend    (Intercept) 11.8     3.43    
 Residual             18.8     4.34    
Number of obs: 20, groups:  blend, 5