[R-sig-ME] random effects specification

[Douglas Bates]

On Thu, Apr 3, 2008 at 6:32 PM, Sebastian P. Luque <spluque at gmail.com> wrote:
> 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.

I'm not sure that I understand what you mean by "treatment being
nested within community".  Does this mean that there are really 8
different treatments because treatment 1 in community A is different
from treatment 1 in community B?  If so, then it would make sense to
me to simply create a new factor that is the interaction of treatment
and community.

Perhaps I am approaching the community factor incorrectly.  In your
data there are two communities so, even if it would be reasonable to
model community effects as random effects, that would be difficult.
With only two levels I think it is best modeled as a fixed effect,
which would mean that questions about treatment and community are
related to the fixed effects.