# Repeated-measures ANOVA looking at interactions for a BACI design in R

Here are my data:

Time    Treatment    Ephemeroptera    Trichoptera    Plecoptera
Before  Control                 98             76            32
After   Control                 91             54            26
Before  Impact                  24             30            12
After   Impact                  29             21            45


Input (example):

Exp = read.csv("fast.csv", header = T)
attach(Exp)
names(Exp)
m1 = aov(Ephemeroptera ~ Time)
summary(m1)
#It gives me output with p-value, df etc. (which is expected)
m2 = aov(Ephemeroptera ~ Treatment)
summary(m2)
#It gives me output with p-value, df etc. (which is expected)
m3 = aov(Ephemeroptera ~ Time*Treatment)
summary(m3)
#It gives me the interaction output (Time:Treatment), but it doesn't provide a p-value.


I am not on my home computer, so I can't include what I got for the output right now. However, I will update it in a little bit (when I get home). Until then, is there any reason why it does not give me any p-values for the output when I try the repeated measures, but it provides it for the regular anovas? Should I try lme instead?

• See here. – gung Apr 8 '13 at 18:14
• Are you actually testing this with only four rows of data? Otherwise, I can't see what's going on. – Gala Jun 11 '13 at 12:23
• @gung But aov() is just regular lm(), not Bates' lmer(). – Gala Jun 11 '13 at 12:24
• @Billy Also note that the code you provided does not define a repeated measures ANOVA model, just a regular 2-way ANOVA. – Gala Jun 11 '13 at 12:25

@GaelLaurans makes a good point that you are thinking of this as a repeated-measures analysis, but you are actually fitting a regular ANOVA. Thus, this isn't the problem of how to correctly determine the denominator's degrees of freedom (which is what I had linked to in my comment above).

I think the issue here is simpler. You have four data points, hence four degrees of freedom. (Note that you lose one degree of freedom automatically when you fit a model for estimating the mean of Ephemeroptera). Your Time factor has two levels, which consumes one degree of freedom. Your Treatment factor has two levels, which also consumes one degree of freedom. If you entered both factors into your model, you would have only one degree of freedom left. If you also entered the interaction term (Time*Treatment), that would also consume one degree of freedom (i.e., your last), and your model would become completely saturated. That is, it is impossible to assess the error variance and / or conduct any statistical tests (for more information about this, see: What is a saturated model?). In short, you simply don't have enough data (as shown) to fit a model with both factors and the interaction.

Here, try this.

Make multiple columns (Velocity, Treatment, Time, Spot, Ephemeroptera, Plecoptera, Trichoptera) in your _.csv file. Then give your Time a value of '1' for Before and '2' for after. I would also suggest inputting another value such as velocity or something else for that value (if you are doing a stream experiment). Then run an lme model for each of your dependent values as your aquatic families. If you have different sections in the area that you sampled...keep them under the Spot variable for each Before and After (1, 2, 3, etc.).

m1 = lme(Ephemeroptera ~ Velocity*Treatment*Time, random = ~ Time|Spot)
summary(m1)


This should get your values. If you want to see your residuals and fitted values, just plug in: plot(m1).

• This could be a valuable approach but doesn't explain the specific problem faced by the OP. – Gala Jun 11 '13 at 12:24