> # 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.
 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(),
1: 76 59 49 74 66
6: 65 75 63 71 84
11: 85 81 61 85 80
16: 74 67 46 89 79
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.
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)
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)
 9.473927 ⇐ half interval width  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" [] (Intercept) A 2.8913263 B -0.7737352 C -13.6014503 D 6.7600022 E 4.7238570 attr(,"class")  "ranef.mer"
> plot(dotplot(re, main="Cleveland Dotplot",
+ 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)
null device 1