[1] "Created: Wed Apr  1 16:35:44 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(abrasion, package="faraway")
abrasion
   run position material wear
1    1        1        C  235
2    1        2        D  236
3    1        3        B  218
4    1        4        A  268
5    2        1        A  251
6    2        2        B  241
7    2        3        D  227
8    2        4        C  229
9    3        1        D  234
10   3        2        C  273
11   3        3        A  274
12   3        4        B  226
13   4        1        B  195
14   4        2        A  270
15   4        3        C  230
16   4        4        D  225
ggplot(abrasion,aes(x=material, y=wear, shape=run, color=position))+geom_point(position = position_jitter(width=0.1, height=0.0))

plot of chunk unnamed-chunk-3

Fit the model

mmod <- lmer(wear ~ material + (1|run) + (1|position), abrasion)
print(summary(mmod),correlation=FALSE)
Linear mixed model fit by REML ['lmerMod']
Formula: wear ~ material + (1 | run) + (1 | position)
   Data: abrasion

REML criterion at convergence: 100.3

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.090 -0.302  0.027  0.422  1.210 

Random effects:
 Groups   Name        Variance Std.Dev.
 run      (Intercept)  66.9     8.18   
 position (Intercept) 107.1    10.35   
 Residual              61.3     7.83   
Number of obs: 16, groups:  run, 4; position, 4

Fixed effects:
            Estimate Std. Error t value
(Intercept)   265.75       7.67    34.7
materialB     -45.75       5.53    -8.3
materialC     -24.00       5.53    -4.3
materialD     -35.25       5.53    -6.4
anova(mmod)
Analysis of Variance Table
         Df Sum Sq Mean Sq F value
material  3   4621    1540    25.1

No inferential output is provided (as expected)

Construct confidence intervals using the parametric bootstrap:

confint(mmod, method="boot")
                           2.5 %  97.5 %
sd_(Intercept)|run        0.0000  16.418
sd_(Intercept)|position   0.0000  19.511
sigma                     3.6584  11.830
(Intercept)             251.0119 281.808
materialB               -57.7385 -35.538
materialC               -34.5255 -12.980
materialD               -45.4872 -24.729

pbkrtest package

We can apply the Kenward-Roger adjustment to test the fixed effect like this:

library(pbkrtest)
mmod <- lmer(wear ~ material + (1|run) + (1|position), abrasion,REML=FALSE)
nmod <- lmer(wear ~ 1+ (1|run) + (1|position), abrasion,REML=FALSE)
KRmodcomp(mmod, nmod)
F-test with Kenward-Roger approximation; computing time: 0.16 sec.
large : wear ~ material + (1 | run) + (1 | position)
small : wear ~ 1 + (1 | run) + (1 | position)
      stat  ndf  ddf F.scaling p.value
Ftest 25.1  3.0  6.0         1 0.00085

We see strong evidence of a difference between materials.

RLRsim package

The RLRsim package cannot be used here as there is more than one random effect. We can implement the parametric bootstrap directly to test whether there is any variance in the runs (for example)

First we compute the observed LRT statistic

rmod <- lmer(wear ~ material + (1|position), abrasion,REML=FALSE)
(olrt <- 2*(logLik(mmod)-logLik(rmod)))
'log Lik.' 5.4164 (df=7)

Now we resample 1000 times:

lrstat <- numeric(1000)
for(i in 1:1000){
rwear <- unlist(simulate(rmod))
bmod <- refit(mmod,rwear)
smod <- refit(rmod,rwear)
lrstat[i] <- 2*(logLik(bmod)-logLik(smod))
}
mean(lrstat > olrt)
[1] 0.036

This gives a p-value close to 0.05 so we might run this a lot more bootstrap samples if we cared deeply about this hypothesis test.

lmerTest package

library(lmerTest)
mmod <- lmer(wear ~ material + (1|run) + (1|position), abrasion,REML=FALSE)
print(summary(mmod),correlation=FALSE)
Linear mixed model fit by maximum likelihood t-tests use Satterthwaite
  approximations to degrees of freedom [merModLmerTest]
Formula: wear ~ material + (1 | run) + (1 | position)
   Data: abrasion

     AIC      BIC   logLik deviance df.resid 
   134.3    139.7    -60.2    120.3        9 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3040 -0.3777  0.0011  0.5490  1.4368 

Random effects:
 Groups   Name        Variance Std.Dev.
 run      (Intercept) 61.4     7.84    
 position (Intercept) 91.1     9.54    
 Residual             41.1     6.41    
Number of obs: 16, groups:  run, 4; position, 4

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)
(Intercept)   265.75       6.96   8.52   38.20  7.9e-11
materialB     -45.75       4.53   8.87  -10.09  3.7e-06
materialC     -24.00       4.53   8.87   -5.29  0.00052
materialD     -35.25       4.53   8.87   -7.77  3.0e-05
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)
material   4621    1540     3  8.87    37.5 2.3e-05

We prefer the Kenward-Rogers (seen above) to the Satterthwaite adjustment shown here.

Can test the random effects:

rand(mmod)
Analysis of Random effects Table:
         Chi.sq Chi.DF p.value
run        3.05      1    0.08
position   4.59      1    0.03

but we prefer the parametric bootstrap result because of problems with asymptotic distribution of the LRT in such cases.

May also be useful to have the pairwise difference of the fixed effects:

difflsmeans(mmod)
Differences of LSMEANS:
               Estimate Standard Error    DF t-value Lower CI Upper CI
material A - B     45.8           4.53   8.9   10.09   35.469   56.031
material A - C     24.0           4.53   8.9    5.29   13.719   34.281
material A - D     35.2           4.53   8.9    7.77   24.969   45.531
material B - C    -21.8           4.53   8.9   -4.80  -32.031  -11.469
material B - D    -10.5           4.53   8.9   -2.32  -20.781   -0.219
material C - D     11.2           4.53   8.9    2.48    0.969   21.531
               p-value
material A - B  <2e-16
material A - C   5e-04
material A - D  <2e-16
material B - C   0.001
material B - D   0.046
material C - D   0.035

All pairwise difference appear significant.

MCMCglmm

We use the expanded prior for the reasons established in our one random effect example. In contrast to the previous nested example the labeling of the factors already indicates a crossed design for the random effects.

library(MCMCglmm)
eprior <- list(R=list(V=1,nu=0.02),G=list(G1=list(V=1,nu=0.02,alpha.V=1000),G2=list(V=1,nu=0.02,alpha.V=1000)))
bmod <- MCMCglmm(wear ~ material,random=~run+position,data=abrasion,verbose=FALSE,prior=eprior,pr=TRUE)
xyplot(log(bmod$VCV))