[R-sig-ME] random effects specification

[Andrew Robinson]

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/