ANOVA: 自由度几乎全部相等为1

Question

我有一个数据集，起始如下：

> d.weight
    R   N   P  C D.weight
1   1   0   0 GO     45.3
2   2   0   0 GO     34.0
3   3   0   0 GO     19.1
4   4   0   0 GO     26.6
5   5   0   0 GO     23.5
6   1  45   0 GO     22.1
7   2  45   0 GO     15.5
8   3  45   0 GO     23.4
9   4  45   0 GO     15.8
10  5  45   0 GO     42.9
...

• R 代表重复，共有5个（1-5）。
• N 代表氮水平，也有5个（0、45、90、180、360）。
• P 代表磷水平，同样有5个（0、35、70、140、280）。
• C 代表植物组合，共有4种（GO、GB、LO、LB）。
• D.weight 是干重，以克为单位。

然而，当我进行方差分析时，自由度会出现错误。通常情况下，我会在完整数据集的子集上运行ANOVA，但是让我们做一次我不会实际执行的分析，这样您就可以看到几乎所有的Df（自由度）都是错误的。

> example.aov=aov(D.weight ~ R+N+P+C, data=d.weight)
> summary(example.aov)
         Df Sum Sq Mean Sq F value  Pr(>F)    
R             1   1158    1158   9.484 0.00226 ** 
N             1    202     202   1.657 0.19900    
P             1  11040   11040  90.408 < 2e-16 ***
C             3  41032   13677 112.010 < 2e-16 ***
Residuals   313  38220     122

基本上，唯一正确的是C因子。这是因为它具有字母而不是数字吗？

我在某个地方发现，如果我在每个术语中写入interaction()，我会得到正确的Df，但我不知道这是否总体上是正确的做法。例如：

> example.aov2=aov(D.weight ~ interaction(R)+interaction(N)+interaction(P)+interaction(C), data=d.weight)
> summary(example.aov2)
                Df Sum Sq Mean Sq F value   Pr(>F)    
interaction(R)   4   7423    1856  19.544 2.51e-14 ***
interaction(N)   4    543     136   1.429    0.224    
interaction(P)   4  13788    3447  36.301  < 2e-16 ***
interaction(C)   3  41032   13677 144.042  < 2e-16 ***
Residuals      304  28866      95

我只尝试了使用C因子来检查是否会出现问题：

> example.aov3=aov(D.weight ~ C, data=d.weight)
> summary(example.aov3)
             Df Sum Sq Mean Sq F value Pr(>F)    
C             3  41032   13677   85.38 <2e-16 ***
Residuals   316  50620     160                   
> 
> example.aov4=aov(D.weight ~ interaction(C), data=d.weight)
> summary(example.aov4)
                Df Sum Sq Mean Sq F value Pr(>F)    
interaction(C)   3  41032   13677   85.38 <2e-16 ***
Residuals      316  50620     160

看起来都一样。我是否应该在各个地方添加interaction()？

- XGF

---

将您的数值变量转换为因子变量...例如：facs <- c("R","N","P"); d_weight[facs] <- lapply(d.weight[facs],factor) - Ben Bolker

这要看情况；你可以将原始数据保留为数字，并在统计模型语句中使用显式的 factor() 语句（例如 D.weight ~ factor(R)+factor(N)+factor(P)），或者将新的因子变量作为辅助变量（例如 fR、fN、fP）而不是覆盖旧变量...如果你提供一个完全可重现的示例 http://tinyurl.com/reproducible-000，说明你想要的图形类型，你可能会得到更明确的帮助。 - Ben Bolker

真的吗？你和你的导师都不知道如何编写数学算法？转化为代码只是过程中的一个小点。 - Carl Witthoft

@BenBolker：谢谢，这两个选项都可以解决问题！我会将它们作为额外的选择添加到您的答案中。 - XGF

@CarlWitthoft：我是一名农业学生，只上了一个学期的统计课。我的导师已经70岁，几乎失明，所以...怎么说呢，我尽力而为。 - XGF

Answer 1

R通过检查变量是numeric还是factor变量来确定是否应将变量视为分类变量（ANOVA类型分析）或连续变量（回归类型分析）。最简单的方法是，您可以通过以下方式将您的预测（独立）变量转换为因子：

facs <- c("R","N","P")
d.weight[facs] <- lapply(d.weight[facs],factor) 

如果您希望创建辅助变量而不是覆盖原始变量，可以这样做：

for (varname in facs) {
   d.weight[[paste0("f",varname)]] <- factor(d.weight[[varname]])
}

可能有更紧凑的方法来做到这一点，但这应该足够了...

- Ben Bolker

---

跟进问题：当你写下“d_weight[facs]”时，你是想创建一个新的变量，还是想写成“d.weight[facs]”？ - daOnlyBG