[R-sig-ME] random effects specification
Sebastian P. Luque
spluque at gmail.com
Tue Apr 8 18:22:30 CEST 2008
Hi,
I think I'm misunderstanding something because I expected a different
result below (building up a reproducible example).
---<---------------cut here---------------start-------------->---
## Simulating data
library(lattice)
library(lme4)
set.seed(1000)
rCom <- rnorm(2, mean=5, sd=0.5)
rTr <- rep(rCom / 1.1, 2)
nbase <- rnorm(240, 10, 0.1)
dta <- within(expand.grid(community=LETTERS[1:2], treatment=letters[1:4],
id=factor(1:30)), {
n <- rCom[as.numeric(community)] +
rTr[as.numeric(treatment)] + nbase
})
dta <- dta[order(dta$community, dta$treatment), ]
## Simulate an interaction
dta$n[dta$community == "A"] <- rev(dta$n[dta$community == "A"])
## Have a look
xyplot(n ~ treatment | community, data=dta, groups=id,
type="b", pch=19, cex=0.3)
## We fit LME with community and treatment fixed effects, their
## interactions, and use random effects for subject.
n.lmer1 <- lmer(n ~ community * treatment + (1 | id), dta)
## Remove the interaction, for comparison.
n.lmer2 <- lmer(n ~ community + treatment + (1 | id), dta)
## Compare these models.
anova(n.lmer1, n.lmer2) # interaction term needed
## So let's test the treatment effect at each community separately.
n.lmerA <- lmer(n ~ treatment + (1 | id), dta,
subset=community == "A")
## Here I expected some terms to be significantly different
## from zero, but that's not the case:
summary(n.lmerA)
---<---------------cut here---------------end---------------->---
Looking at the data, there seems to be an obvious treatment effect, yet
this doesn't show in the lmer summary. What am I misunderstanding here?
Thanks.
Cheers,
--
Seb
More information about the R-sig-mixed-models
mailing list