[R-sig-ME] lmer model specification

[Ben Bolker]

jude girard <jude.girard at ...> writes:

 [snip]
 
> In my study, we compared invertebrate biomass in paired organic and
> conventional fields.  I had 9 conventional and 9 organic fields, for a
> total of 18. Each field was sampled in May and in July.  At each
> sampling time, two transects, each of 4 traps, were placed in each
> field.
> 
> My fixed effects are management type and month (I expect invertebrate
> abundance to be higher in July than in May), and the interaction
> between them.  At first I specified my random effects as
> pair/month/transect.  However, when I run this model
> 
> model.edge<-lmer(biomass~management*month+(1|pair/month/transect), data=edge)
> summary(model.edge)
> 
> Linear mixed model fit by REML
> Formula: biomass ~ management + month + (1 | pair/month/transect)
>    Data: edge
>    AIC   BIC logLik deviance REMLdev
>  535.8 559.8 -260.9    513.4   521.8
> Random effects:
>  Groups                Name        Variance Std.Dev.
>  transect:(month:pair) (Intercept) 0.071400 0.26721
>  month:pair            (Intercept) 0.000000 0.00000
>  pair                  (Intercept) 0.018794 0.13709
>  Residual                          0.491406 0.70100
> Number of obs: 230, groups: transect:(month:pair), 36; 
> month:pair, 18; pair, 9
> 
> Fixed effects:
>               Estimate Std. Error t value
> (Intercept)   -1.84429    0.11147 -16.545
> managementorg  0.29042    0.09474   3.065
> monthR2        0.45889    0.12924   3.551

 [snip]

> the variance for month:pair is 0, which I think might be because the
> month is already specified in the fixed effects?  

  No: month is specified in the fixed effects only as its
main effect, whereas the variance specifies month:pair interaction
only.  You're getting an estimate of zero variance because there
is essentially not any detectable among-(month:pair) variation
over and above that expected from transect:month:pair variation.
This is essentially as expected when trying to fit several levels
of nesting to a moderate-sized data set.

> So, now I am
> wondering if a better way to run the model might be
> model2.edge<-lmer(biomass~management*month+(1|pair/transect.unique),
> data=edge)  where is transect is treated as a unique level,
> irrespective of which month it was run in?

   I believe that this would then leave out the effect
of "month within pair" (still trying to figure out why 0.1
of the variance comes out of the residual variance and is
counted within the transect:unique.pair variance instead here).

> When I run the second model I get
> 
> summary(model2.edge)
> Linear mixed model fit by REML
> Formula: biomass ~ management * month + (1 | pair/transect.unique)
>    Data: edge
>    AIC   BIC logLik deviance REMLdev
>  524.2 548.3 -255.1    501.7   510.2
> Random effects:
>  Groups               Name        Variance Std.Dev.
>  transect.unique:pair (Intercept) 0.17567  0.41913
>  pair                 (Intercept) 0.01887  0.13737
>  Residual                         0.39778  0.63070
> Number of obs: 230, groups: transect.unique:pair, 70; pair, 9
> 
> Fixed effects:
>                       Estimate Std. Error t value
> (Intercept)            -1.8028     0.1374 -13.122
> managementorg           0.1665     0.1866   0.892
> monthR2                 0.3745     0.1842   2.033
> managementorg:monthR2   0.1866     0.2645   0.705
> 

  The bigger difference in this model is that by putting the
interaction (month:management) in you are allowing for
month-to-month variation in the effect of management, which
wasn't in the previous model.  These two effects (whether
you include a pair:month variance term, and whether you
include a month:management fixed effect) are more or less
orthogonal, I believe.  (Note that adding the interaction
term screws up your power to detect the management effect ...)

  I would also consider running this in lme -- it should
give you the same answers, and it will give you denominator
df and p-values (which should be reasonably well justified 
in this classical, balanced linear mixed model).