数据准备
df <- read.table(file = "D:/Documents/R wd/df.csv", header = T, sep = ",") # 数据导入。
df # 查看数据。
df$nitrogen <- as.factor(df$nitrogen) # 将nitrogen转为因子。
## 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
9.4 单因素协方差分析
单因素协方差分析（ANCOVA）扩展了单因素方差分析（ANOVA），包含一个或多个定量的
协变量。
attach(df) # 添加df为路径存储索引。
table(nitrogen) # 查看分组变量信息。
## nitrogen
## N1 N2
## 12 12
aggregate(v1, by = list(nitrogen), FUN = mean) # 分组统计v1，按nitrogen分组。
## Group.1 x
## 1 N1 1.2000
## 2 N2 1.3125
fit14 <- aov(v1 ~ block + nitrogen, data = df) # 协方差分析，协变量为block。
summary(fit14) # 返回分析结果。
## Df Sum Sq Mean Sq F value Pr(>F)
## block 1 0.00562 0.00562 0.753 0.39523
## nitrogen 1 0.07594 0.07594 10.170 0.00441 **
## Residuals 21 0.15680 0.00747
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
结果解读：block对应的p值大于0.05，说明block与v1不相关，v1的值显著受nitrogen的影响。
由于使用了协变量，你可能想要获取调整的组均值——即去除协变量效应后的组均值。可使
用effects包中的effects()函数来计算调整的均值：
library(effects) # 调用effects包。
effect("nitrogen", fit14) # 计算调整后的均值。
##
## nitrogen effect
## nitrogen
## N1 N2
## 1.2000 1.3125
library(multcomp) # 调用multcomp包。
contrast <- rbind("N1 vs N2" = c(1,-1)) # N1和N2组的比较。
summary(glht(fit14, linfct=mcp(nitrogen=contrast))) # 返回比较结果。
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: aov(formula = v1 ~ block + nitrogen, data = df)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## N1 vs N2 == 0 -0.11250 0.03528 -3.189 0.00441 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
9.4.1 评估检验的假设条件
ANCOVA与ANOVA相同，都需要正态性和同方差性假设，ANCOVA还假定回归斜率相同。ANCOVA模型包含block*nitrogen的交互项时，可对回归斜率的同质性进行检验。交互效应若显著，则意味着block和v1间的关系依赖于nitrogen的水平。若不显著，支持了斜率相等的假设。
library(multcomp) # 调用multcomp包。
fit15 <- aov(v1 ~ block*nitrogen, data = df) # 检验斜率相等的假设。
summary(fit15) # 返回结果。
## Df Sum Sq Mean Sq F value Pr(>F)
## block 1 0.00562 0.00562 0.812 0.37830
## nitrogen 1 0.07594 0.07594 10.960 0.00349 **
## block:nitrogen 1 0.01822 0.01822 2.630 0.12050
## Residuals 20 0.13858 0.00693
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
9.4.2 结果可视化
library(HH) # 调用HH包。
ancova(v1 ~ block + nitrogen, data = df) # 可视化结果。
## Analysis of Variance Table
##
## Response: v1
## Df Sum Sq Mean Sq F value Pr(>F)
## block 1 0.005625 0.005625 0.7533 0.395229
## nitrogen 1 0.075937 0.075937 10.1702 0.004414 **
## Residuals 21 0.156800 0.007467
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
image.png
图形解读：用block预测v1的回归线相互平行，N2下截距较N1大。
参考资料：
- 《R语言实战》（中文版），人民邮电出版社，2013.
