Repeated measures ANOVA in R and Bonferroni adjusted intervals

If somebody can tell me what R commands I need to use for a repeated measures ANOVA, I'd really appreciate it. I have trouble with the random term. I've seen random=id, random=id/(treatment*group) and others.

Also can you please indicate to me what the formula for the Bonferroni adjusted intervals is?

-
What do you mean by Bonferroni adjusted intervals? Or rather, what do you mean by intervals? Why kind of intervals? Intervals around what? –  John Jul 24 '11 at 23:18
Thank you for your answer. I mean when the repeated measures ANOVA detects a difference between the means, what are the Bonferroni intervals to know which time point and between which groups contributed to this difference. –  sara Jul 25 '11 at 4:28
see my edited response and accept an answer –  John Jul 25 '11 at 11:49
@sara Please, don't forget to register your account, vote on response(s) that were helpful (or not), and accept the one that fit with your original question. –  chl Jul 25 '11 at 20:08

I agree with suncoolsu that it is difficult to tell exactly what you were looking for. And in general it's not recommended to do standard repeated measures ANOVAs anymore since there are generally better alternatives.

Nevertheless, perhaps you want to generate a simple stratified ANOVA. By stratified I mean that your effects are measured within another grouping variable, in your case the subject and thus a within subjects design. If your data frame is df and your response variable is y then you might have a within subjects predictor x1 and that crossed with a within subjects predictor x2, and perhaps a between subjects predictor z. To get the full model with all interactions you would use.

myModel <- aov( y ~ x1 * x2 * z + Error(id/(x1*x2)), data = df)
summary(myModel)


You'll note that within the Error term we are grouping x1, x2, and their interaction under id. Note that z is not in there because it is not a within subjects variable.

Keep in mind further that this data is laid out in long format and you probably need to aggregate it first to run this correctly since a repeated measures design often suggests more samples / subject than conditions in order to get good estimates of each subject's response value. Therefore, df above might be replaced with the following dfa.

dfa <- aggregate ( y ~ x1 + x2 + z + id, data = df, mean)


(BTW, suncoolsu gave a much more modern answer based on multi-level modelling. It's suggested you learn about that if you continue to do repeated measures designs because it is much more powerful, flexible, and allows one to ignore certain kinds of within subjects assumptions (notably sphericity). What I've described is how to do repeated measure ANOVA. You also might want to look at the car, or higher level ez packages in order to do it as well.)

As for your Bonferroni query... it should probably have been a separate question. Nevertheless, that's a bit of a hard one to answer with repeated measures. You could try ?pairwise.t.test. If you give the interactions of all your within variables as the group factor and set paired to true and the correction to bonf you're set. However, straight corrections like that probably are far too conservative. You state at the outset you're only going to use it if there is a significant effect, you probably also have a theoretical reason for making some comparisons, therefore it's not strictly the fishing expedition that Bonferroni (over) corrects for. So, something like...

with( df, pairwise.t.test(y, x1:x2, paired = TRUE, p.adj = 'bonf') )


will do what you want but that's probably not really what you want.

-
Thank you so much. I was actually trying to use the aov function and your response was of big help. –  sara Jul 25 '11 at 4:25
Thank you. That definitely helps. –  sara Jul 25 '11 at 16:19

You should provide more details about your data. From the limited details provided by you, assuming you have a data frame df which has response, trt, time, and subject information, then these are many ways to fit a LME model in R using lme4 package. However, I will illustrate three methods that I think will be useful for you.

library(lme4)

# Random intercepts for different subjects but time and trt effects are fixed
mmod1 <- lmer(response ~ time*trt + (1 | subject), df)

# Random intercepts and trt effects for different subjects, but the time effect is still fixed
mmod2 <- lmer(response ~ time*trt + (1 + trt | subject), df)

# Random intercepts, trt, and time effects for different subjects
mmod3 <- lmer(response ~ time*trt + (1 + trt + time | subject), df)


Once you have the p-values from the model fit above, you can use:

HPDinterval(mmod1, prob = 0.95, ...)
HPDinterval(mmod2, prob = 0.95, ...)
HPDinterval(mmod3, prob = 0.95, ...)


to obtain the 95% CI determined from MCMC sample. Since this CI is obtained from MCMC sampling, it takes into account of the random errors and you won't need to correct for multiple comparisons (I think so, please correct me if I am wrong).

-
Thank you so much. I have a data frame that has a response measured at different time points and different treatments. So response, time and treatment. Could you please help me when to choose a random intercept model and when a random intercept and treatment model? –  sara Jul 25 '11 at 4:23
@sara. I think going with a random intercept and slope is a good choice, unless you have a reason that the slopes won't vary with the time. –  suncoolsu Jul 25 '11 at 5:51
Thanks for the help :) –  sara Jul 25 '11 at 16:11
Hello again. How was the time taken into account here? Thanks –  sara Jul 25 '11 at 18:31
I edited my answer and made the time dependence explicit. HTH. –  suncoolsu Jul 25 '11 at 20:23