I admit defeat and ask for help.
I am trying to analyze the results of whole lake experiments on the effect
of fish farming in lakes on nutrients and carbon budgets. Two lakes (Lake)
have been monitored for 4 years (Year). In one of the lakes, a fish farm
was operated in the last two years. This is a standard
Before-After-Control-Impact (BACI) design. Period (before or after)
correspond to either the first 2 years or the last 2 years of the
experiment. The dependent variable (CRRC) is monitored weekly in each lake
during the summer. Over the summer, there is a strong seasonal variation and
we expect these longitudinal data to exhibit some temporal autocorrelation.
Following the best I can Pinheiro and Bates,I fitted models with and without
correlation structures:
mod1<-lme(CRRC~Lake*Period, random=~1|Year)
mod2<-update(mod1,corr=corAR1())
mod2 has a significantly better fit than mod1.
But the model term of interest (the Lake:Period interaction) is not
significant in mod2. It is very significant in mod1, and indeed, the effect
is so obvious on graphs that it is not even funny. During the first two
years, the two lakes have an almost identical temporal pattern, and there is
a complete divergence in the last two years where one like become
essentially constant while the other has a pattern similar to the first 2
years.
The size of the interaction term in mod2 appears incorrect: in mod1, the
interaction effect is twice that of the lake effect, and 4 times that of the
period, whereas in the mod2 model, the interaction is more than order of
magnitude smaller than either of the main effects. I have the nasty feeling
that I have incorrectly specified the model and/or grouped data.
Are mod1 and mod2 correctly specified?
Antoine Morin
Biologie, Université d'Ottawa
|