Solved – Blocked repeated measures ANOVA in R

anovarrepeated measures

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?

Best 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.

Related Question