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?

Improve this question
$\endgroup$
6
  • 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
    Commented Apr 26, 2013 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
    Commented Apr 26, 2013 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
    Commented Apr 26, 2013 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
    Commented Apr 26, 2013 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
    Commented Apr 26, 2013 at 2:18

1 Answer 1

Reset to default
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.

Improve this answer
$\endgroup$
2
  • $\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
    Commented Apr 26, 2013 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
    Commented Apr 26, 2013 at 2:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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