> # 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,σ_{u}^{2}) where
σ_{u}^{2} 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