[R-sig-ME] random effects specification

[Sebastian P. Luque]

On Thu, 3 Apr 2008 17:32:40 -0500,
"Douglas Bates" <bates at stat.wisc.edu> wrote:

[...]

> It seems that the observations are indexed by subject and treatment so
> the number of levels in the factor treatment:id equals the number of
> observations.  You can't estimate a variance for such a term and also
> estimate a residual variance.

> I would start with

> n ~ treatment * community +(1|id)

Yes, the observations are indexed by subject and treatment in the sense
that id levels are the same within treatments of the same community, but
are different among communities.  This is a subset of the data:

---<---------------cut here---------------start-------------->---
id  community treatment     n
 1          A         1 13.93
 2          A         1 14.42
 3          A         1 13.56
 1          A         2 14.61
 2          A         2 14.74
 3          A         2 15.59
 1          A         3 13.95
 2          A         3 15.21
 3          A         3 14.51
 1          A         4 13.61
 2          A         4 14.99
 3          A         4 15.13
 4          B         1 14.79
 5          B         1 13.41
 6          B         1 14.71
 4          B         2 14.69
 5          B         2 13.46
 6          B         2 14.28
 4          B         3 14.30
 5          B         3 13.18
 6          B         3 13.58
 4          B         4 14.54
 5          B         4 13.25
 6          B         4 14.09
---<---------------cut here---------------end---------------->---

Of course, there are many more individuals, but the levels of id differ
among communities, and are the same among treatments.  lmer did converge
rapidly with your suggested formula though:

---<---------------cut here---------------start-------------->---
Linear mixed-effects model fit by REML 
Formula: n ~ treatment * community + (1 | id) 
   Data: isotope.m.ph 
 AIC BIC logLik MLdeviance REMLdeviance
 450 481   -216        410          432
Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.193    0.439   
 Residual             0.232    0.481   
number of obs: 240, groups: id, 61

Fixed effects:
                               Estimate Std. Error t value
(Intercept)                     14.9748     0.1170   128.0
treatment2                       0.0884     0.1222     0.7
treatment3                      -0.2829     0.1222    -2.3
treatment4                       0.3568     0.1222     2.9
communitysanikiluaq             -0.5749     0.1678    -3.4
treatment2:communitysanikiluaq  -0.2471     0.1763    -1.4
treatment3:communitysanikiluaq  -0.7479     0.1763    -4.2
treatment4:communitysanikiluaq  -0.6169     0.1763    -3.5

Correlation of Fixed Effects:
            (Intr) trtmn2 trtmn3 trtmn4 cmmnty trtm2: trtm3:
treatment2  -0.522                                          
treatment3  -0.522  0.500                                   
treatment4  -0.522  0.500  0.500                            
cmmntysnklq -0.697  0.364  0.364  0.364                     
trtmnt2:cmm  0.362 -0.693 -0.347 -0.347 -0.524              
trtmnt3:cmm  0.362 -0.347 -0.693 -0.347 -0.524  0.503       
trtmnt4:cmm  0.362 -0.347 -0.347 -0.693 -0.524  0.503  0.503
---<---------------cut here---------------end---------------->---


However, I don't understand how (1 | id) accounts for treatment being
nested within community.  Maybe it's time for me to re-read some more
examples from "Mixed-effects models in S and S-plus".  Thanks.


Cheers,

-- 
Seb