[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