[R-sig-ME] random effects specification
Sebastian P. Luque
spluque at gmail.com
Thu May 1 15:28:11 CEST 2008
Hi again,
I've made further explorations into lmer, toying with the example I
showed earlier:
---<---------------cut here---------------start-------------->---
set.seed(1000)
rCom <- rnorm(2, mean=5, sd=0.5)
rTr <- rep(rCom / 1.1, 2)
nbase <- rnorm(60, 10, 0.1)
## 20 individuals; 10 in community "A" and 10 in "B", each receiving 3
## different treatments once.
dta <- within(expand.grid(community=LETTERS[1:2], treatment=letters[1:3],
id=1:10), {
id[community == "B"] <- id[community == "B"] + 10
id <- as.factor(id)
n <- rCom[as.numeric(community)] +
rTr[as.numeric(treatment)] + nbase
})
dta <- dta[order(dta$id, 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 test for community (A, B) and treatment (a, b, c) fixed effects,
## their interactions, and use random effects for subject (1:20). Am I
## writing this correctly?
##
## y_{ijk} = B_j + B_k + B_{jk} + b_i + e_{ijk} i=1,...,20; j=A,B; k=a,b,c
n.lmer1 <- lmer(n ~ community * treatment + (1 | id), dta)
---<---------------cut here---------------end---------------->---
I'm a bit confused whether I'm describing the model being fit correctly
(not the lmer call, but the model description in the comment above), and
how it could be described in matrix form. I think this type of
exercises would help me grasp the syntax conventions better.
Another issue is that the lmer call results in a warning:
---<---------------cut here---------------start-------------->---
Warning message:
In .local(x, ..., value) :
Estimated variance for factor ‘id’ is effectively zero
---<---------------cut here---------------end---------------->---
which I presume is due to the fact that the data are unreplicated,
i.e. individuals get each treatment only once. Are there any gotchas in
the interpretation of the results after this warning?
Thanks once again!
--
Seb
More information about the R-sig-mixed-models
mailing list