[R-sig-ME] random effects specification
Andrew Robinson
A.Robinson at ms.unimelb.edu.au
Thu May 1 23:30:27 CEST 2008
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/
More information about the R-sig-mixed-models
mailing list