2
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Glen_b -Reinstate Monica Apr 26 '13 at 2:07
  • $\begingroup$ 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 $\endgroup$ – Relena Ribbons Apr 26 '13 at 2:12
  • $\begingroup$ 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). $\endgroup$ – gung - Reinstate Monica Apr 26 '13 at 2:13
  • $\begingroup$ 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)? $\endgroup$ – Relena Ribbons Apr 26 '13 at 2:15
  • $\begingroup$ @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. $\endgroup$ – Glen_b -Reinstate Monica Apr 26 '13 at 2:18
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Relena Ribbons Apr 26 '13 at 2:21
  • $\begingroup$ @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? $\endgroup$ – Glen_b -Reinstate Monica Apr 26 '13 at 2:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.