# Blocked repeated measures ANOVA in R

I am trying to test if richness varies by treatment (Severity) over time (Block) using a repeated measures ANOVA in R. Any suggestions on how to do this correctly? I have tried, but get an error message stating I have an extra = in the formula:

TWP.aov <- aov(Richness ~ Severity * Year + Error(Plot/Severity), data = TWP)
summary(TWP.aov)


Here is a subset of the data (called TWP which I read in as a csv file):

Block       Severity    Plot    Richness
2003-2004   High        A       18
2003-2004   High        B       24
2003-2004   High        C       21
2005-2006   High        A       28
2005-2006   High        B       24
2005-2006   High        C       20
2007-2009   High        A       14
2007-2009   High        B       27
2007-2009   High        C       29
2003-2004   Low         A       12
2003-2004   Low         B       10
2003-2004   Low         C       14
2005-2006   Low         A       18
2005-2006   Low         B       16
2005-2006   Low         C       14
2007-2009   Low         A       8
2007-2009   Low         B       19
2007-2009   Low         C       20


When I input the above subset I get the correct summary table:

Error: Plot
Df Sum Sq Mean Sq F value Pr(>F)
Residuals  2  49.33   24.67

Error: Plot:Severity
Df Sum Sq Mean Sq F value Pr(>F)
Severity   1 304.22  304.22   85.56 0.0115 *
Residuals  2   7.11    3.56
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Block           2  43.00  21.500   0.745  0.505
Severity:Block  2   1.44   0.722   0.025  0.975
Residuals       8 230.89  28.861


My actual data's summary table looks like this though:

Error: Plot
Df Sum Sq Mean Sq F value Pr(>F)
Severity        1    248  247.66   4.013 0.0479 *
Block           2    493  246.67   3.997 0.0214 *
Severity:Block  2    244  122.17   1.980 0.1436
Residuals      99   6110   61.71
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


As of June 12, 2013: As an update, I still have not resolved the "Error() model is singular" issue with the aov function, so I have tried to run this in the nlme package using the following code:

TWP.lme <- lme(Richness ~Severity * Block, random = ~1|(Plot/Severity), data = TWP)
summary(TWP.lme)
anova(TWP.lme)


I do not get any error message, and I believe I have structured the error term correctly. I have checked the results visually against the graphs of the data and all appears correct. Is there a way to check the results are correct? Why is there no error message regarding the model is singular?

• Hi Relena and welcome to CV. I have taken the liberty of using the formatting tools to make your post more readable. Can you paste in the exact message you get? Please note that the data you show there doesn't have a Year variable at all, it has one called Block instead. Apr 26, 2013 at 2:07
• yes, sorry about that- I was trying to make it reproduce-able for CV. Thank you for the formatting help. Here is my error message: Warning message: In aov(Richness ~ Severity * Block + Error(Plot/Severity), data = TWP) : Error() model is singular Apr 26, 2013 at 2:12
• This appears to be only about how to get something done in R, & not about any related statistical content, as such, it would be off-topic for CV (see our FAQ), but on-topic at Stack Overflow. If you have a substantive question about statistics, please edit to clarify, if not, flag your Q & we'll migrate it for you (please don't cross-post, though). Apr 26, 2013 at 2:13
• I am concerned as to whether I am using the appropriate type of ANOVA for this data, but I can post this on Stack Overflow instead (how do I migrate the question)? Apr 26, 2013 at 2:15
• @gung the 'model is singular' message the OP mentioned a few minutes ago (doubtless while you were typing there) may arguably relate to a statistical issue. I'd argue to leave it here for the moment. Apr 26, 2013 at 2:18

## 1 Answer

Your data there works fine for me as soon as I put 'Block' in the model instead of 'Year'

Check the data reads in okay:

> summary(TWP)
Block   Severity Plot     Richness
2003-2004:6   High:9   A:6   Min.   : 8.00
2005-2006:6   Low :9   B:6   1st Qu.:14.00
2007-2009:6            C:6   Median :18.50
Mean   :18.67
3rd Qu.:23.25
Max.   :29.00


Fit the model, then fix it so it has the right names:

> TWP.aov <- aov(Richness ~ Severity * Year + Error(Plot/Severity), data = TWP)
Error in eval(expr, envir, enclos) : object 'Year' not found
> TWP.aov <- aov(Richness ~ Severity * Block + Error(Plot/Severity), data = TWP)


After that, it seems to work:

> summary(TWP.aov)

Error: Plot
Df Sum Sq Mean Sq F value Pr(>F)
Residuals  2  49.33   24.67

Error: Plot:Severity
Df Sum Sq Mean Sq F value Pr(>F)
Severity   1 304.22  304.22   85.56 0.0115 *
Residuals  2   7.11    3.56
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Block           2  43.00  21.500   0.745  0.505
Severity:Block  2   1.44   0.722   0.025  0.975
Residuals       8 230.89  28.861


What message did you get?

---

Warning message: In aov(Richness ~ Severity * Block + Error(Plot/Severity), data = TWP) : Error() model is singular


Please note that "model is singular" has nothing to do with there being extra "=" in the formula. (How did you get that interpretation of the message?).

Looking at the code for aov, it looks like it relates to the rank of the design matrix being lower than the number of coefficients on completing the QR deomposition. This would suggest you have multicollinearity.

• I tried with block instead of year, and the code worked but then I believe my p-values are incorrect because when the full dataset is graphed out with 95% confidence intervals, the error bars do not overlap for the final block (so a statistically sig. difference between high and low severity). Shouldn't the p-value for the interaction Severity:Block also be significant in that case? Apr 26, 2013 at 2:21
• @RelenaRibbons Could you edit the summary table into the bottom of your question, please, and then use the {} tool to turn it into 'code' (which will make it monospaced), and thereby render it readable? Apr 26, 2013 at 2:31