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).