22.2 驱虫喷雾的效果

InsectSprays 数据集 (Beall 1942) 来源于农业实验,记录了不同杀虫剂的效果,即杀虫剂过后,单位实验区域内虫子的数量,如图22.1所示,横轴表示杀虫剂种类,纵轴表示虫子数量。

ggplot(data = InsectSprays, aes(x = spray, y = count, color = spray)) +
  geom_boxplot() +
  geom_jitter() +
  theme_minimal()
不同杀虫剂的效果

图 22.1: 不同杀虫剂的效果

先创建一个 aov 对象,把它命名为 mod1,见下方

mod1 <- aov(count ~ spray, data = InsectSprays)

第一个参数告诉 R count 是响应变量,spray 是协变量,第二个参数告诉 R 去对象 InsectSprays 中寻找这些变量。下面把分析结果以一种漂亮的格式打印出来

summary(mod1)
##             Df Sum Sq Mean Sq F value Pr(>F)    
## spray        5   2669   533.8    34.7 <2e-16 ***
## Residuals   66   1015    15.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

表格中的条目是很容易理解的,比如最右边的列表示 P 值。如果我们想做固定显著性水平下的检验,比如 \(\alpha = 0.075\) 时的 F 统计量的值,

qf(0.075, 5, 66, lower.tail = F)
## [1] 2.110783

上面的命令是说 \(F(5, 66)\) 分布的 0.075 分位点,最后一个参数很关键,因为默认情况下 R 计算下分位点,详情见 ?qf

方差分析做了三个假设

  1. 残差 \(\epsilon_{ij}\) 是相互独立的随机变量;
  2. 残差 \(\epsilon_{ij}\) 服从正态分布;
  3. 残差 \(\epsilon_{ij}\) 均值为 0,方差是固定的常数。

假设 1 和 3 通过图来检验,假设 2 通过 QQ 图来检验。值得一提的是 mod1 对象除了打印出来,还有很多方法

names(mod1)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"

比如获取残差,考虑到篇幅,这里仅显示前 10 个

head(mod1$residuals, 10)
##    1    2    3    4    5    6    7    8    9   10 
## -4.5 -7.5  5.5 -0.5 -0.5 -2.5 -4.5  8.5  2.5  5.5
par(mar = c(4, 4, 2, 2), mfrow = c(2,2))
plot(mod1)

plot(mod1$fitted.values, mod1$residuals, main = "Residuals vs. Fitted", pch = 20)
abline(h = 0, lty = 2)

plot(mod1$model$spray, mod1$residuals, main = "Residuals vs. Levels" )

plot(1:72, mod1$residuals, main = "Residuals vs. Time Order")
abline(h = 0, lty = 2)

qqnorm(mod1$residuals, pch = 20)
qqline(mod1$residuals)

如果上面的假设显著失效,我们要采用非参数检验

mod2 <- kruskal.test(count ~ spray, data = InsectSprays)
mod2
## 
##  Kruskal-Wallis rank sum test
## 
## data:  count by spray
## Kruskal-Wallis chi-squared = 54.691, df = 5, p-value = 1.511e-10

计算给定水平下的置信区间,构造置信水平为 95% 的区间

\[ \bar{X} \pm t_{1-\alpha/2}(s/\sqrt{n}) \] 以 A 号杀虫剂为例,

xbar = mean(InsectSprays[InsectSprays$spray == "A", "count"])
t_crit <- qt(0.025, mod1$df.residual, lower.tail = F)
s <- sqrt(sum((mod1$residuals)^2) / mod1$df.residual)
n <- sum(InsectSprays$spray == "A")
# 最后置信区间的上下限
c(xbar - t_crit * (s/ sqrt(n)), xbar + t_crit * (s/ sqrt(n)))
## [1] 12.23958 16.76042

比较 A 号和 C 号杀虫剂的效果,计算两个均值差的置信区间

\[ \bar{X}_1 - \bar{X}_2 \pm t_{1-\alpha/2}(s/\sqrt{1/n_1 + 1/n_2}) \]

n1 <- sum(InsectSprays$spray == "A")
n2 <- sum(InsectSprays$spray == "C")

x1bar = mean(InsectSprays[InsectSprays$spray == "A", "count"])
x2bar = mean(InsectSprays[InsectSprays$spray == "C", "count"])

代入公式即可计算得到置信区间

(x1bar - x2bar) - t_crit * s * sqrt( 1/ n1 + 1/n2)
## [1] 9.219948
(x1bar - x2bar) + t_crit * s * sqrt( 1/ n1 + 1/n2)
## [1] 15.61339

Fisher’s 最小显著性检验 (Fisher’s Least Significant Difference Test) 即

t_crit * s * sqrt( 1/ n1 + 1/n2)
## [1] 3.196719

Tukey’s Honestly Significant Difference Test 主要测量成对实验的误差比率,假定每个水平下的实验次数是相等的,只需将上面的 aov 对象传递给函数 TukeyHSD()

mod3 <- TukeyHSD(mod1, ordered = TRUE)
mod3
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
##     factor levels have been ordered
## 
## Fit: aov(formula = count ~ spray, data = InsectSprays)
## 
## $spray
##           diff       lwr       upr     p adj
## E-C  1.4166667 -3.282742  6.116075 0.9488669
## D-C  2.8333333 -1.866075  7.532742 0.4920707
## A-C 12.4166667  7.717258 17.116075 0.0000000
## B-C 13.2500000  8.550591 17.949409 0.0000000
## F-C 14.5833333  9.883925 19.282742 0.0000000
## D-E  1.4166667 -3.282742  6.116075 0.9488669
## A-E 11.0000000  6.300591 15.699409 0.0000000
## B-E 11.8333333  7.133925 16.532742 0.0000000
## F-E 13.1666667  8.467258 17.866075 0.0000000
## A-D  9.5833333  4.883925 14.282742 0.0000014
## B-D 10.4166667  5.717258 15.116075 0.0000002
## F-D 11.7500000  7.050591 16.449409 0.0000000
## B-A  0.8333333 -3.866075  5.532742 0.9951810
## F-A  2.1666667 -2.532742  6.866075 0.7542147
## F-B  1.3333333 -3.366075  6.032742 0.9603075

其中,diff 表示均值之差,lwrupr 表示置信区间的上下限,\(p_{adj}\) 是对应的。检查一下,看看哪些置信区间包含 0 ,包含 0 的表示不显著,从第三行来看, A 和 C 之间差别显著。之前计算过 A、C 均值,均值之差即

(x1bar - x2bar)
## [1] 12.41667

在误差比率 \(\alpha = 0.05\) 的情况下,如果你想手动计算 HSD 值

q_crit <- qtukey(p = 0.05, nmeans = length(mod1$xlevels[[1]]), df = mod1$df.residual, lower.tail = F)
# mod1$df.residual 是 6 
hsd <- q_crit * s / sqrt(6)
hsd
## [1] 6.645967

将模型结果 mod3 用图画出来,见下图

plot(mod3)
成对显著性水平

图 22.2: 成对显著性水平

关于多重比较请见 Frank Bretz, Torsten Hothorn, Peter Westfall 的书《Multiple Comparisons Using R》及配套 R 包 multcomp,该 R 包现由 Torsten Hothorn 维护,他还维护了一个由数据集构成的 R 包 TH.data,我们后续章节也会用到。

参考文献

Beall, Geoffrey. 1942. “The Transformation of Data from Entomological Field Experiments so That the Analysis of Variance Becomes Applicable.” Biometrika 32 (3/4): 243–62. https://doi.org/10.2307/2332128.