[R-sig-ME] help with repeated measures on a split-plot experiment

[Reinhold Kliegl]

On Tue, Jun 10, 2008 at 3:56 PM, James Hudson <jmghudson at gmail.com> wrote:
> Thank you Thierry, Reinhold and Hadley for the comments. The graphic
> plots are extremely helpful.
>
> Here's more about the experiment:
> The data comes from a simulated climate change experiment where we
> measured arctic tundra vegetation in 1m2 plots over time. Warm "P"=
> non-warmed, Warm "T" = warmed. Snow "A" = "decrease growing season
> length, Snow "R" = increase growing season length, Snow "C" = no
> manipulation. The dataset provided (for deciduous shrubs) is just one
> of about 15 we have for different vegetation classes. I will use the
> same code for each class.

So using C-level of snow as reference in a treatment contrast is
justified a-priori; also the Helmert contrast. But why does both
decreasing and increasing the growing season increase your response
linearly over the years? And why would you get a stronger increase
over the years for non-warmed  (P) than warmed (T) plots?

>
> Two design questions:
>
> A. Like I mentioned above, I want to re-use this code repeatedly. Are
> there methodological issues with repeating these univariate analyses?

Why would you want to do this? You can easily include the data from
the 15 different vegation classes in one analysis and use the full
statistical power of your design. This will provide you with important
additional information. For example,  in such an analyses you can test
differences between clusters of vegetation classes that you would
cluster together (i.e., by using contrasts on the factor veg. class).
If you are not interested in differences between vegetation classes,
you could specify it as another random factor. Then "plot" can be
specified as nested within vegetation class or as crossed with it and
you can estimate the variances of these factors as well as the
variances of effects associated with them. I would recommend to do a
bit of reading on this matter. (For presentation of results, it may be
useful to run the 15 separate analyses, too.)

> Afterwards, I will also conduct other community-level analyses such as
> adonis.
Not sure what this means.

>
> B. There was an error in the dataset I posted previously. When the
> experiment was set-up, the 36 plots were spatially paired. Each pair
> consisted of a warmed and non-warmed plot (with one Snow level). I
> should have labelled the plots 1-18 x2 rather than 1-36. Correct?
>
Correct.
In this case, the repeated-measures anova statement  must be changed
to take this into consideration

summary(model0 <- aov(response ~  year*snow*warm +
Error(plot/(year*warm)), data=shrubs))

lmer figures it out by itself. The mean plots will change in the
standard errors because I computed them after aggregating to the plot
level.

Now you can also test wether there is reliable variance in the
warm-effect between plots. It turns out, there is:

> print(model7.REML<-lmer(response~(year.L+snow.C+warm)^2  + (warm|plot), method="REML", data=shrubs), cor=FALSE)
Linear mixed model fit by REML
Formula: response ~ (year.L + snow.C + warm)^2 + (warm | plot)
   Data: shrubs
   AIC   BIC logLik deviance REMLdev
 638.7 676.2 -305.3      618   610.7
Random effects:
 Groups   Name        Variance Std.Dev. Corr
 plot     (Intercept) 10.998   3.3163
          warm T      16.402   4.0499   -0.429
 Residual             12.402   3.5216
Number of obs: 108, groups: plot, 18

Fixed effects:
                Estimate Std. Error t value
(Intercept)       9.2260     1.5907   5.800
year.L            0.1609     0.1377   1.169
snow.C1          -0.5507     2.2487  -0.245
snow.C2          -0.6121     2.2487  -0.272
warm T            0.1748     2.0298   0.086
year.L:snow.C1    0.4094     0.1687   2.427
year.L:snow.C2    0.4182     0.1687   2.480
year.L:warm T    -0.3456     0.1377  -2.509
snow.C1:warm T   -0.2778     2.8676  -0.097
snow.C2:warm T   -1.3889     2.8676  -0.484
> anova(model5.REML, model7.REML)  # not statistically reliable
Data: shrubs
Models:
model5.REML: response ~ (year.L + snow.C + warm)^2 + (1 | plot)
model7.REML: response ~ (year.L + snow.C + warm)^2 + (warm | plot)
            Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)
model5.REML 12  650.24  682.42 -313.12
model7.REML 14  645.95  683.50 -308.98 8.2811      2    0.01591 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1