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

[James Hudson]

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.

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?
Afterwards, I will also conduct other community-level analyses such as
adonis.

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?

Here's the updated dataset:
date , warm , snow , plot , response
1995 , P , A , 1 , 11
1995 , P , C , 2 , 11
1995 , P , A , 3 , 4
1995 , P , C , 4 , 5
1995 , P , R , 5 , 6
1995 , P , A , 6 , 3
1995 , P , R , 7 , 5
1995 , P , A , 8 , 3
1995 , P , C , 9 , 7
1995 , P , R , 10 , 8
1995 , P , A , 11 , 7
1995 , P , C , 12 , 13
1995 , P , C , 13 , 4
1995 , P , R , 14 , 7
1995 , P , R , 15 , 1
1995 , P , C , 16 , 6
1995 , P , R , 17 , 5
1995 , P , A , 18 , 1
1995 , T , A , 1 , 9
1995 , T , C , 2 , 13
1995 , T , A , 3 , 5
1995 , T , C , 4 , 5
1995 , T , R , 5 , 6
1995 , T , A , 6 , 9
1995 , T , R , 7 , 5
1995 , T , A , 8 , 15
1995 , T , C , 9 , 12
1995 , T , R , 10 , 7
1995 , T , A , 11 , 8
1995 , T , C , 12 , 12
1995 , T , C , 13 , 19
1995 , T , R , 14 , 5
1995 , T , R , 15 , 10
1995 , T , C , 16 , 3
1995 , T , R , 17 , 5
1995 , T , A , 18 , 2
2000 , P , A , 1 , 22
2000 , P , C , 2 , 19
2000 , P , A , 3 , 9
2000 , P , C , 4 , 3
2000 , P , R , 5 , 5
2000 , P , A , 6 , 5
2000 , P , R , 7 , 16
2000 , P , A , 8 , 7
2000 , P , C , 9 , 12
2000 , P , R , 10 , 7
2000 , P , A , 11 , 12
2000 , P , C , 12 , 10
2000 , P , C , 13 , 11
2000 , P , R , 14 , 12
2000 , P , R , 15 , 3
2000 , P , C , 16 , 11
2000 , P , R , 17 , 12
2000 , P , A , 18 , 5
2000 , T , A , 1 , 10
2000 , T , C , 2 , 13
2000 , T , A , 3 , 3
2000 , T , C , 4 , 4
2000 , T , R , 5 , 7
2000 , T , A , 6 , 5
2000 , T , R , 7 , 4
2000 , T , A , 8 , 23
2000 , T , C , 9 , 11
2000 , T , R , 10 , 10
2000 , T , A , 11 , 5
2000 , T , C , 12 , 8
2000 , T , C , 13 , 14
2000 , T , R , 14 , 7
2000 , T , R , 15 , 9
2000 , T , C , 16 , 0
2000 , T , R , 17 , 7
2000 , T , A , 18 , 2
2007 , P , A , 1 , 15
2007 , P , C , 2 , 19
2007 , P , A , 3 , 5
2007 , P , C , 4 , 1
2007 , P , R , 5 , 18
2007 , P , A , 6 , 14
2007 , P , R , 7 , 15
2007 , P , A , 8 , 9
2007 , P , C , 9 , 11
2007 , P , R , 10 , 9
2007 , P , A , 11 , 22
2007 , P , C , 12 , 3
2007 , P , C , 13 , 13
2007 , P , R , 14 , 6
2007 , P , R , 15 , 11
2007 , P , C , 16 , 9
2007 , P , R , 17 , 16
2007 , P , A , 18 , 9
2007 , T , A , 1 , 9
2007 , T , C , 2 , 20
2007 , T , A , 3 , 9
2007 , T , C , 4 , 5
2007 , T , R , 5 , 5
2007 , T , A , 6 , 7
2007 , T , R , 7 , 7
2007 , T , A , 8 , 18
2007 , T , C , 9 , 14
2007 , T , R , 10 , 12
2007 , T , A , 11 , 16
2007 , T , C , 12 , 7
2007 , T , C , 13 , 3
2007 , T , R , 14 , 15
2007 , T , R , 15 , 9
2007 , T , C , 16 , 4
2007 , T , R , 17 , 6
2007 , T , A , 18 , 2



On Tue, Jun 10, 2008 at 3:59 AM, Reinhold Kliegl
<reinhold.kliegl at gmail.com> wrote:
> Hi,
>
> Nice example, also to contrast with standard repeated-measures ANOVA,
> plot of significant interactions with ggplot, and use of contrasts.
> What, by the way, are three levels of snow?
>
> library(lme4, lib.loc=.libUsr)
> library(ggplot2)
>
> shrubs<-as.data.frame(read.csv("dshrubs.csv", header=TRUE))
> shrubs$plot <- as.factor(shrubs$plot)
> shrubs$year <- as.factor(shrubs$date)
> contrasts(shrubs$year) <- contr.poly
>
> # standard repeated-measures ANOVA for this balanced design: shows
> significant interaction of year x warm
> summary(model0 <- aov(response ~  year*snow*warm + Error(plot/year),
> data=shrubs))
>
> # lmer equivalen, with polynomial contrasts for year and treatment
> contrast for snow (reference=A)
> print(model1<-lmer(response~ year*snow*warm + (1|plot), method="ML",
> data=shrubs), cor=FALSE)
>
> # Using only linear trend of year
> shrubs$yearL <- as.numeric(as.character(shrubs$year))-2000  # center
> at year 2000
> print(model2<-lmer(response~yearL*snow*warm  + (1|plot), method="ML",
> data=shrubs), cor=FALSE)
>
> # dropping 3-factor interactions, clearly suggests two significant
> interactions with year
> print(model3<-lmer(response~(yearL+snow+warm)^2  + (1|plot),
> method="ML", data=shrubs), cor=FALSE)
> anova(model1, model2, model3)
>
> # plot two significant interactions
> shrubs.rs <- melt(shrubs, id=c("plot", "yearL", "warm", "snow"),
> measure="response")
>
> # (1) year x temperature
> table<-cast(shrubs.rs, yearL+warm ~ .,
>           function(x) c(M=mean(x), SE=sd(x)/sqrt(length(x)), N=length(x) ))
> p <- qplot(x=yearL, y=M, data=table, shape=warm, colour=warm) +
> scale_x_continuous("Year", breaks=c(-5, 0, 7), labels=c("1995",
> "2000", "2007"))
> p <- p  + geom_point(size=3) + geom_line(aes(group=warm), size=1) +
> geom_errorbar(aes(max=M+SE, min=M-SE), width=0.1)
> p
>
> # (2) year x snow
> table<-cast(shrubs.rs, yearL+snow ~ .,
>           function(x) c(M=mean(x), SE=sd(x)/sqrt(length(x)), N=length(x) ))
> p <- qplot(x=yearL, y=M, data=table, shape=snow, colour=snow) +
> scale_x_continuous("Year", breaks=c(-5, 0, 7), labels=c("1995",
> "2000", "2007"))
> p <- p  + geom_point(size=3) + geom_line(aes(group=snow), size=1) +
> geom_errorbar(aes(max=M+SE, min=M-SE), width=0.1)
> p
>
> # Options on factor snow (tailored to pattern of means; should ideally
> be specified a priori)
> # ... refit lmer model with helmert contrasts for snow: (1) c vs. a+r;
> (2) a vs. r
> shrubs$snowH <- C(shrubs$snow, matrix(c(-1, +2, -1,    -1,  0, +1), 3, 2), 2)
> contrasts(shrubs$snowH)
> print(model4<-lmer(response~(yearL+snowH+warm)^2  + (1|plot),
> method="ML", data=shrubs), cor=FALSE)
>
> # ... refit lmer model with treatment contrasts for snow using C as
> reference: (1) c vs. a+r; (2) a vs. r
> shrubs$snowC <- C(shrubs$snow, matrix(c(1, 0, 0,     0, 0, 1), 3, 2), 2)
> contrasts(shrubs$snowC)
> print(model5<-lmer(response~(yearL+snowC+warm)^2  + (1|plot),
> method="ML", data=shrubs), cor=FALSE)
>
> # re-specification of snow factor does not change overall goodness of fit
> anova(model3, model4, model5)
>
> # check for reliable between-plot variance for linear effect of year;
> switch from full ML to restricted ML
> print(model5.REML<-lmer(response~(yearL+snowC+warm)^2  + (1|plot),
> method="REML", data=shrubs), cor=FALSE)
> print(model6.REML<-lmer(response~(yearL+snowC+warm)^2  + (yearL|plot),
> method="REML", data=shrubs), cor=FALSE)
> anova(model5.REML, model6.REML)  # not statistically reliable
>
> Reinhold Kliegl
>
> On Tue, Jun 10, 2008 at 6:41 AM, hadley wickham <h.wickham at gmail.com> wrote:
>>>> If you can make do with time as a linear function then I think that's OK. I've still got a a feeling that if you have time as a factor then you run out of dfs. If you can post a self-contained example with model and data then others can probably comment if this is the case (I received some great help from the list on this topic about a year ago.
>>>>
>>>
>>> Here is my script with data below (including time as a linear function
>>> rather than a factor):
>>> library(nlme)
>>> dataset<-as.data.frame(read.csv("dshrubs.csv", header=TRUE))
>>> attach(dataset)
>>> names(dataset)
>>> model1<-lme(response~date*snow*warm, random =~ 1|plot, data=dataset)
>>> anova(model1)
>>> plot(model1,resid(.)~plot,abline=0)
>>> qqnorm(model1)
>>
>> I'd suggest you start with some good explanatory plots:
>>
>> install.packages("ggplot2")
>> library(ggplot2)
>> qplot(date, response, data=shrubs, colour = warm, group=plot,
>> geom=c("line", "point"), facets = snow ~ .)
>>
>> This seems revealing to me: not much going on, with the possible
>> exception of snow = "R", which has a lower variance and a slight
>> upward trend, particularly for warm = P.  Group-wise linear models
>> support this interpretation:
>>
>> qplot(date, response, data=shrubs, colour = warm, geom="point", facets
>> = snow ~ .) +
>> geom_smooth(method=lm)
>>
>> Although if the group-level variances truly are equal, you will get
>> more power from the mixed effects model.
>>
>> Hadley
>>
>>
>> --
>> http://had.co.nz/
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>