random effects specification
Hi Sebastian,
On Thu, May 01, 2008 at 08:28:11AM -0500, Sebastian P. Luque wrote:
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
I think so.
how it could be described in matrix form.
I'm not sure what you're looking for here. Matrix form as I understand it tends to be pretty generic. But, have a go and we can critique it :)
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.
No, I think that it's because you didn't assign any subject to subject variation to the observations. n comprises a contribution from community, from treatment, and from the base error, but not from id. Try adding id-based variation.
Are there any gotchas in the interpretation of the results after this warning?
I would ordinarily interpret that message as meaning that the model is overparameterized. That seems to be true, here. Andrew
Andrew Robinson Department of Mathematics and Statistics Tel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/