Dear Mixed Ones,
While waiting for Doug Bates's book, I'm trying trying to understand
the differences between the recommendations of Pinheiro & Bates using
the nlme packages, and the use of lmer in the lme4 package. I've
gotten myself confused.
#because lme4 and nlme sometimes clash:
detach('package:lme4')
require(nlme)
data(Machines)
mach <- Machines
# Machines and Workers are balanced
with(mach, table(Machine, Worker))
# Here is how it's analyzed in Pinheiro & Bates
m1 <- lme(score ~ Machine, data = mach, random = ~ 1 | Worker)
m2 <- update(m1, random = ~ 1 | Worker / Machine)
anova(m1, m2)
# Back to lme4
detach('package:nlme')
require(lme4)
mr1 <- lmer(score ~ Machine + (1 | Worker), data = mach)
mr2 <- update(mr1, . ~ Machine + (1 | Worker / Machine))
mr3 <- update(mr1, . ~ Machine + (1 + Machine | Worker))
anova(mr1, mr2, mr3)
require(gmodels)
options(digits = 3)
set.seed(20080427)
zapsmall(ci(m2))
zapsmall(ci(mr2, sim.lmer = T, n.sim = 10000))
# The following produces SEs an order of magnitude larger ci(mr2)
zapsmall(ci(mr3, sim.lmer = T, n.sim = 10000))
# From this it would seem that mr2 is the proper model.
# But I have never seen a crossed design written that way in lme4.
# Another example:
# Baayen (in the page proofs of Analyzing Linguistic Data:
# A practical introduction to statistics, p. 283) takes an example from
# Raaijmakers, J. G. W.; Schrijnemakers, J. M. C. & Gremmen, F.
# How to Deal with ``The Language-as-Fixed-Effect Fallacy'':
# Common Misconceptions and Alternative Solutions
# Journal of Memory and Language, 1999, 41, 416-426
subj <- factor(rep(paste('S', 1:8, sep = ''), each = 8))
item <- factor(rep(paste('W', 1:8, sep = ''), 8))
soa <- factor(rep(c("Short", "Long"), each = 4, 8))
rt <- c(
546, 567, 547, 566, 554, 545, 594, 522,
566, 556, 538, 566, 512, 523, 569, 524,
567, 598, 568, 584, 536, 539, 589, 521,
556, 565, 536, 550, 516, 522, 560, 486,
595, 609, 585, 588, 578, 540, 615, 546,
569, 578, 560, 583, 501, 535, 568, 514,
527, 554, 535, 527, 480, 467, 540, 473,
551, 575, 558, 556, 588, 563, 631, 558)
sp <- data.frame(subj = subj, item = item, soa = soa, rt = rt)
# soa and subj are crossed.
with(sp, table(subj, item))
# item and subj also
with(sp, table(subj, soa))
# items nested under soa
with(sp, table(item, soa))
# pr2 is analogous to m2 and mr2
pr2 <- lmer(rt ~ soa + (1 | subj / soa) + (1 | item), sp)
# pr3 is analogous to m3 and mr3 This is how Baayen analyzes it
# (the results aren't identical to his; I don't know why):
pr3 <- lmer(rt ~ soa + (1 + soa | subj) + (1 | item), sp)
ci(pr2, sim.lmer = TRUE, n.sim = 10000)
ci(pr3, sim.lmer = TRUE, n.sim = 10000)
cm <- unique(model.matrix(~ soa, sp))
estimable(pr2, cm = cm, sim.lmer = TRUE, n.sim = 10000, conf.int = 0.95)
estimable(pr3, cm = cm, sim.lmer = TRUE, n.sim = 10000, conf.int = 0.95)
_____________________________
Professor Michael Kubovy
University of Virginia
Department of Psychology
USPS: P.O.Box 400400 Charlottesville, VA 22904-4400
Parcels: Room 102 Gilmer Hall
McCormick Road Charlottesville, VA 22903
Office: B011 +1-434-982-4729
Lab: B019 +1-434-982-4751
Fax: +1-434-982-4766
WWW: http://www.people.virginia.edu/~mk9y/