One-Way Random Effects Model

Apex Enterprises Example (pp 191-197)

> # Function gl is used to generate factor levels.

> # The first argument gives the number of levels.

> # Argument 2 gives the number of replications.

> # Argument 3 gives the total sample size (balanced).

> # Argument labels gives the labels of the levels.

> gl(5,1,20,labels=LETTERS[1:5])

 [1] A B C D E A B C D E A B C D E A B C D E
Levels: A B C D E

> apex <- data.frame(rating=scan(),

+     officer=gl(5,1,20,labels=LETTERS[1:5]))

1:  76 59 49 74 66

6:  65 75 63 71 84

11: 85 81 61 85 80

16: 74 67 46 89 79

21:

Read 20 items

> # The package lme4 (which in turn requires R (≥2.4.1))

> # needs to be installed. This package can handle

> # linear mixed effect (representation) models and beyond.

> # In fact, it can handle generalized mixed effect models

> # (i.e. non-Gaussian variance components).

> # I have installed it, so the library can be loaded now.

> require(lme4)

Loading required package: lme4
Loading required package: lattice
Loading required package: Matrix

> # The rhs (right hand side) of the model specification

> # consists of two parts, all fixed effects are specified

> # without parentheses and random effects are enclosed

> # within parentheses. The variable following a '|'

> # indicates a grouping variable (factor).

> apex.lmer <- lmer(rating~1 + (1|officer), data=apex)

> summary(apex.lmer)

Linear mixed model fit by REML ['lmerMod']
Formula: rating ~ 1 + (1 | officer) 
   Data: apex 

REML criterion at convergence: 145.2452 

Random effects:
 Groups   Name        Variance Std.Dev.
 officer  (Intercept) 80.41    8.967   
 Residual             73.28    8.561   
Number of obs: 20, groups: officer, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)   71.450      4.444   16.08

> # Note that the variance components were reported.

> # To get a 90% c.i. for the grand mean μ, use

> # the information reported in the Fixed effects

> # section (see page 1039 of textbook):

> 71.450 + c(-1,1) * print(qt(.05,4,lower=F)*4.444)

[1] 9.473927    ⇐ half interval width
[1] 61.97607 80.92393    ⇐ 90% c.i.

> #

> # The session below gives random effects and their

> # graphical representation (known as caterpillar

> # plot). These random effects resemble a sample

> # from N(0,σu2) where σu2 is estimated by

> # the variance component corresponding to the

> # random effect factor officer: 80.42.

> print(re <- ranef(apex.lmer, condVar=TRUE))

An object of class "ranef.lmer"
[[1]]
  (Intercept)
A   2.8913263
B  -0.7737352
C -13.6014503
D   6.7600022
E   4.7238570

attr(,"class")
[1] "ranef.mer"

> windows(width=5,height=3,pointsize=10)

> plot(dotplot(re, main="Cleveland Dotplot",

+              scale=list(x=list(alternating=3)))$office,

+      position=c(0,0,0.5,1), more=T)

> plot(qqmath(re, main="Caterpillar Plot")$office,

+      position=c(0.45,0,1,1), more=F)

> dev.off()

null device 
          1

(apex01.png here)