`[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()`

```
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
```

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
```

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.

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.

```
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
```