One fixed effect and one random effect
Julian Faraway
[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()
ggplot(penicillin,aes(x=treat,y=yield,group=blend,linetype=blend))+geom_line()
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.500anova(mmod)Analysis of Variance Table
Df Sum Sq Mean Sq F value
treat 3 70 23.3 1.24The 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.1367Inference 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.34The 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.26Several 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.0429Again 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.500anova(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.34These 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.1The 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.30MCMCglmm 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)
xyplot(log(bmod$VCV))
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.430The 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()
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()
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")
and also for the error SD:
hist(sqrt(bmod$VCV[,2]), 50, xlab="yield", main="Error SD")
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)
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