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

[Reinhold Kliegl]

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
>