数据准备
df <- read.table(file = "D:/Documents/R wd/df.csv", header = T, sep = ",") # 数据导入。
df # 查看数据。
## year nitrogen variety block v1 v2 v3 v4 v5
## 1 2020 N1 a 1 1.26 2.14 0.4 5.0 3.25
## 2 2020 N1 a 2 1.20 2.90 0.1 5.3 1.27
## 3 2020 N1 a 3 1.30 3.00 0.3 5.6 2.24
## 4 2020 N1 b 1 1.08 1.72 1.8 2.8 1.00
## 5 2020 N1 b 2 1.05 1.65 1.7 2.5 3.12
## 6 2020 N1 b 3 1.15 1.35 1.5 3.1 4.57
## 7 2020 N2 a 1 1.32 3.78 1.6 6.0 5.85
## 8 2020 N2 a 2 1.28 4.32 1.4 6.1 6.48
## 9 2020 N2 a 3 1.35 3.95 1.3 6.2 7.21
## 10 2020 N2 b 1 1.33 3.47 2.8 4.1 6.56
## 11 2020 N2 b 2 1.28 2.72 2.4 4.3 8.43
## 12 2020 N2 b 3 1.30 3.90 2.2 4.5 7.55
## 13 2021 N1 a 1 1.19 3.61 0.8 6.0 3.11
## 14 2021 N1 a 2 1.21 3.29 0.5 5.7 2.54
## 15 2021 N1 a 3 1.24 3.26 0.7 5.6 1.28
## 16 2021 N1 b 1 1.09 2.71 1.8 4.0 3.24
## 17 2021 N1 b 2 1.28 2.32 1.6 4.2 1.27
## 18 2021 N1 b 3 1.35 1.95 1.3 4.3 1.15
## 19 2021 N2 a 1 1.45 4.35 1.8 7.2 5.74
## 20 2021 N2 a 2 1.40 3.80 1.2 7.0 6.85
## 21 2021 N2 a 3 1.37 4.23 1.6 6.8 7.42
## 22 2021 N2 b 1 1.28 2.72 2.4 5.1 8.20
## 23 2021 N2 b 2 1.15 3.35 2.5 5.5 5.70
## 24 2021 N2 b 3 1.24 3.46 2.7 4.9 6.00
df$nitrogen <- as.factor(df$nitrogen) # 将nitrogen转为因子。
9.5 双因素方差分析
attach(df) # 将df加入搜索路径。
df$variety <- as.factor(df$variety) # 将variety转为因子。
table(nitrogen,variety) # 查看因子基本信息。
## variety
## nitrogen a b
## N1 6 6
## N2 6 6
aggregate(v1, by = list(nitrogen, variety), FUN = mean) # 分组统计平均值。
## Group.1 Group.2 x
## 1 N1 a 1.233333
## 2 N2 a 1.361667
## 3 N1 b 1.166667
## 4 N2 b 1.263333
fit15 <- aov(v1 ~ nitrogen*variety) # 双因素方差分析。
summary(fit15) # 返回fit15结果。
## Df Sum Sq Mean Sq F value Pr(>F)
## nitrogen 1 0.07594 0.07594 12.647 0.00198 **
## variety 1 0.04084 0.04084 6.802 0.01683 *
## nitrogen:variety 1 0.00150 0.00150 0.251 0.62217
## Residuals 20 0.12008 0.00600
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
结果可视化
方法1：interaction.plot()函数来展示双因素方差分析的交互效应。
interaction.plot(nitrogen,variety,v1) # 结果可视化。
image.png
图形解读：无论哪个品种，v1值均是N2显著高于N1，品种来看，a品种显著高于b品种。
方法2：gplots包中的plotmeans()函数来展示交互效应。
library(gplots) # 调用gplots包。
plotmeans(v1 ~ interaction(nitrogen, variety, sep = "")) # 绘图。
image.png
图形解读：每个样点都显示了置信区间，对应横坐标上方显示了样本量。本例中，同一品种下，均是N2显著高于N1
方法3：HH包中的interaction2wt()函数来可视化结果。
library(HH) # 调用HH包。
interaction2wt(v1 ~ nitrogen*variety) # 绘图。
image.png
图形解读：两个箱图（左下角和右上角）是主效应图，左下角图，单就nitrogen看，N2显著高于N1；右上角图，单就variety看，a品种显著高于b品种；两个线图（左上角和右下角）为交互效应图，左上图，nitrogen应用在不同variety下的影响，即不同品种下，都是N2显著高于N1，同时a品种显著高于b品种；右下图，不同品种在不同nitrogen下的影响，即不同nitrogen下，都是a品种显著高于b品种，同时N2显著高于N1。
参考资料：
- 《R语言实战》（中文版），人民邮电出版社，2013.
