[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

Fixed effects:
            Estimate Std. Error    df t value Pr(>|t|)
(Intercept)    84.00       2.48 11.07   33.94  1.5e-12
treatB          1.00       2.75 12.00    0.36    0.722
treatC          5.00       2.75 12.00    1.82    0.094
treatD          2.00       2.75 12.00    0.73    0.480

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 of type 3  with  Satterthwaite 
approximation for degrees of freedom
      Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
treat     70    23.3     3    12    1.24   0.34

These results match those previously obtained.

Can also test the random effects term:

rand(mmod)
Analysis of Random effects Table:
      Chi.sq Chi.DF p.value
blend   2.76      1     0.1

The produces the LRT statistic and uses the chi-squared approximation. This result is not reliable and the parametric bootstrap result from RLRsim should be preferred.

We can also obtain pairwise difference between the fixed effects:

difflsmeans(mmod)
Differences of LSMEANS:
            Estimate Standard Error DF t-value Lower CI Upper CI p-value
treat A - B       -1           2.74 12   -0.36    -6.98     4.98    0.72
treat A - C       -5           2.74 12   -1.82   -10.98     0.98    0.09
treat A - D       -2           2.74 12   -0.73    -7.98     3.98    0.48
treat B - C       -4           2.74 12   -1.46    -9.98     1.98    0.17
treat B - D       -1           2.74 12   -0.36    -6.98     4.98    0.72
treat C - D        3           2.74 12    1.09    -2.98     8.98    0.30

MCMCglmm package

See the One way anova analysis for an introduction to MCMCglmm. We proceed directly with the final prior (using parameter expansion):

library(MCMCglmm)
eprior <- list(R=list(V=1,nu=0.02),G=list(G1=list(V=1,nu=0.02,alpha.V=1000)))
bmod <- MCMCglmm(yield ~ treat, ~blend, data=penicillin,verbose=FALSE,prior=eprior,pr=TRUE)
xyplot(bmod$Sol)

plot of chunk unnamed-chunk-15

xyplot(log(bmod$VCV))

plot of chunk unnamed-chunk-15

The diagnostics are satisfactory

summary(bmod)

 Iterations = 3001:12991
 Thinning interval  = 10
 Sample size  = 1000 

 DIC: 125.3 

 G-structure:  ~blend

      post.mean l-95% CI u-95% CI eff.samp
blend        26 0.000219     80.6     1000

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units      23.5     8.45     43.4     1000

 Location effects: yield ~ treat 

            post.mean l-95% CI u-95% CI eff.samp  pMCMC
(Intercept)    84.008   78.130   89.886     1000 <0.001
treatB          0.997   -4.945    7.312      813  0.728
treatC          4.933   -1.039   10.696     1000  0.094
treatD          2.105   -4.179    7.726     1000  0.430

The output shows pMCMC which is twice the posterior probability observed on the other side of zero from the posterior mean. This is the same way in which a bootstrap p-value would be calculated and might be viewed as a “Bayesian p-value” (a dubious analogy). These suggest no strong evidence of a difference in the treatments.

We can plot the posterior distributions.

library(reshape2)
ref <- data.frame(bmod$Sol[,2:4])
rdf <- melt(ref)
ggplot(rdf, aes(x=value,color=variable))+geom_density()

plot of chunk unnamed-chunk-17

These are differences with the A treatment levels. As we see all three overlap zero substantially hence the large p-value above. Furthermore, the three distributions all overlap considerably so there is no evidence of any difference between the four treatment levels.

We can also look at the blend posterior distributions:

ref <- data.frame(bmod$Sol[,5:9])
colnames(ref) <- as.character(1:5)
rdf <- melt(ref)
colnames(rdf) <- c("blend","yield")
ggplot(rdf, aes(x=yield,color=blend))+geom_density()

plot of chunk unnamed-chunk-18

We see that the blends overlap considerably and it would be difficult to tell them apart reliably. We can look at the posterior distribution for the blend SD:

hist(sqrt(bmod$VCV[,1]), 50, xlab="yield", main="Blend SD")

plot of chunk unnamed-chunk-19

and also for the error SD:

hist(sqrt(bmod$VCV[,2]), 50, xlab="yield", main="Error SD")

plot of chunk unnamed-chunk-20

We can also look at the joint distribution:

rdf <- data.frame(bmod$VCV)
ggplot(rdf, aes(x=sqrt(blend),y=sqrt(units)))+geom_density2d()+geom_abline(int=0,slope=1)

plot of chunk unnamed-chunk-21

This suggests that the error SD might be bigger than the blend SD although it is not clearcut.

Check on versions of R packages used:

sessionInfo()
R version 3.1.1 (2014-07-10)
Platform: x86_64-apple-darwin13.1.0 (64-bit)

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] reshape2_1.2.2   MCMCglmm_2.21    ape_3.1-4        coda_0.16-1     
 [5] lattice_0.20-29  lmerTest_2.0-20  RLRsim_3.0       pbkrtest_0.4-1  
 [9] MASS_7.3-33      ggplot2_0.9.3.1  lme4_1.1-7       Rcpp_0.11.3     
[13] Matrix_1.1-4     knitr_1.6        rmarkdown_0.2.68

loaded via a namespace (and not attached):
 [1] bitops_1.0-6        caTools_1.16        cluster_1.15.2     
 [4] colorspace_1.2-4    corpcor_1.6.7       dichromat_2.0-0    
 [7] digest_0.6.4        evaluate_0.5.3      formatR_0.10       
[10] Formula_1.1-1       gdata_2.13.3        gplots_2.13.0      
[13] grid_3.1.1          gtable_0.1.2        gtools_3.3.1       
[16] Hmisc_3.14-3        htmltools_0.2.4     KernSmooth_2.23-12 
[19] labeling_0.2        latticeExtra_0.6-26 minqa_1.2.3        
[22] munsell_0.4.2       nlme_3.1-117        nloptr_1.0.4       
[25] numDeriv_2012.9-1   plyr_1.8.1          proto_0.3-10       
[28] RColorBrewer_1.0-5  scales_0.2.3        splines_3.1.1      
[31] stringr_0.6.2       survival_2.37-7     tensorA_0.36       
[34] tools_3.1.1         yaml_2.1.11