2

I'm trying to run anova() in R and running into some difficulty. This is what I've done up to now to help shed some light on my question.

Here is the str() of my data to this point.

 str(mhw)
'data.frame':   500 obs. of  5 variables:
$ r    : int  1 2 3 4 5 6 7 8 9 10 ...
$ c    : int  1 1 1 1 1 1 1 1 1 1 ...
$ grain: num  3.63 4.07 4.51 3.9 3.63 3.16 3.18 3.42 3.97 3.4 ...
$ straw: num  6.37 6.24 7.05 6.91 5.93 5.59 5.32 5.52 6.03 5.66 ...
$ Quad : Factor w/ 4 levels "NE","NW","SE",..: 2 2 2 2 2 2 2 2 2 2 ...

Column r is a numerical value indicating which row in the field an individual plot resides Column c is a numerical value indicating which column an individual plot resides
Column Quad corresponds to the geographical location in the field to which each plot resides

Quad <- ifelse(mhw$c > 13 & mhw$r < 11, "NE",ifelse(mhw$c < 13 & mhw$r < 11,"NW", ifelse(mhw$c < 13 & mhw$r >= 11, "SW","SE")))
mhw <- cbind(mhw, Quad)

I have fit a lm() as follows

nov.model <-lm(mhw$grain ~ mhw$straw)
anova(nov.model)

This is an anova() for the entire field, which is testing grain yield against straw yield for each plot in the dataset.

My trouble is that I want to run an individual anova() for the Quad column of my data to test grain yield and straw yield in each quadrant.

perhaps a with() might fix that. I have never used it before and I am in the process of learning R currently. Any help would be greatly appreciated.

Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
pc8807
  • 47
  • 1
  • 6

1 Answers1

5

I think you are looking for by facility in R.

fit <- with(mhw, by(mhw, Quad, function (dat) lm(grain ~ straw, data = dat)))

Since you have 4 levels in Quad, you end up with 4 linear models in fit, i.e., fit is a "by" class object (a type of "list") of length 4.

To get coefficient for each model, you can use

sapply(fit, coef)

To produce model summary, use

lapply(fit, summary)

To export ANOVA table, use

lapply(fit, anova)

As a reproducible example, I am taking the example from ?by:

tmp <- with(warpbreaks,
            by(warpbreaks, tension,
               function(x) lm(breaks ~ wool, data = x)))

class(tmp)
# [1] "by"

mode(tmp)
# [1] "list"

sapply(tmp, coef)

#                    L         M         H
#(Intercept)  44.55556 24.000000 24.555556
#woolB       -16.33333  4.777778 -5.777778

lapply(tmp, anova)

#$L
#Analysis of Variance Table
#
#Response: breaks
#          Df Sum Sq Mean Sq F value  Pr(>F)  
#wool       1 1200.5 1200.50  5.6531 0.03023 *
#Residuals 16 3397.8  212.36                  
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#$M
#Analysis of Variance Table
#
#Response: breaks
#          Df  Sum Sq Mean Sq F value Pr(>F)
#wool       1  102.72 102.722  1.2531 0.2795
#Residuals 16 1311.56  81.972               
#
#$H
#Analysis of Variance Table
#
#Response: breaks
#          Df  Sum Sq Mean Sq F value Pr(>F)
#wool       1  150.22 150.222  2.3205 0.1472
#Residuals 16 1035.78  64.736

I was aware of this option, but not familiar with it. Thanks to @Roland for providing code for the above reproducible example:

library(nlme)
lapply(lmList(breaks ~ wool | tension, data = warpbreaks), anova)

For your data I think it would be

fit <- lmList(grain ~ straw | Quad, data = mhw)
lapply(fit, anova)

You don't need to install nlme; it comes with R as one of recommended packages.

Community
  • 1
  • 1
Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
  • Thank you. This provides the `coef` of the data. If I wanted to produce the full anova summary to test the hypothesis, "Does grain yield differ in East Quadrants to West Quadrants". How can I produce that with `fit` – pc8807 Sep 29 '16 at 12:04
  • Worked Perfect. Thank you very much for the help. I apologize for the simple nature of the questions. As I said I'm in the process of learning R and get a little turned around on ,my syntax. – pc8807 Sep 29 '16 at 12:12
  • That link is quite helpful. That cleared up many of my side questions that came up as I was working with R. – pc8807 Sep 29 '16 at 12:37
  • Or use `lmList`: `library(nlme); lapply(lmList(breaks ~ wool | tension, data = warpbreaks), anova)` – Roland Sep 29 '16 at 13:59