r - 按组拟合线性模型/方差分析

Question

这个问题已经有答案了:

Linear Regression and group by in R (10 个回答)

已关闭 7 年前。

我正在尝试运行 anova()在 R 中遇到一些困难。这就是我到目前为止所做的，以帮助阐明我的问题。

这是 str()到目前为止我的数据。

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

r 列是一个数值，指示单个图位于字段中的哪一行c 列是一个数值，指示单个图位于哪一列
Column Quad 对应于每个地 block 所在田地的地理位置

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)

我已经安装了lm()如下

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

这是一个anova()对于整个田地，正在测试数据集中每个地 block 的 Cereal 产量与秸秆产量。

我的麻烦是我想运行一个个人anova()对于我的数据的四列来测试每个象限的粮食产量和秸秆产量。

也许是with()可能会解决这个问题。我以前从未使用过它，目前正在学习 R。任何帮助将不胜感激。

Answer 1

我认为您正在寻找 R 中的 by 工具。

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

由于 Quad 中有 4 个级别，因此 fit 中最终会得到 4 个线性模型，即 fit 是一个“by”长度为 4 的类对象(“列表”类型)。

要获取每个模型的系数，您可以使用

sapply(fit, coef)

要生成模型摘要，请使用

lapply(fit, summary)

要导出方差分析表，请使用

lapply(fit, anova)

---

作为一个可重现的示例，我从 ?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

---

我知道这个选项，但并不熟悉。感谢@Roland为上述可重现示例提供代码:

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

对于你的数据，我认为是

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

您不需要安装nlme；它附带 R 作为推荐包之一。