对于多组样本平均数差异的显著性检验, 这类问题通常是研究影响随机变量值的因素。例如,同一类作物有 5 种甚至更多的品系, 不同品系的产量可以相等或者存在差异, 可以对不同品系的产量平均数进行差异检验以选择高产品系; 又如, 吸烟对肺功能是否有影响? 吸烟因素可以分为不吸烟、被动吸烟、少量吸烟、中等程度吸烟、重度吸烟等不同水平, 通过随机选择在不同水平上的人群, 利用能反映肺功能的生理指标一一用力呼气量 (FEV), 检测不同人群 FEV 的平均数差异。

这两个例子是单因素多水平的平均数检验问题, 用双样本假设检验进行两两水平对的分析不是好方法, 存在计算量大、多重检验带来高假阳性率等问题,一般采用单因素方差分析的方法。

此外, 影响变量值的因素有多个, 例如影响肺功能的因素除了吸烟之外,还与年龄、性别、基础疾病、遗传背景等因素相关,这是多个因素、每个因素有多个水平的问题, 需要用到多因子方差分析的方法。方差分析可以减少两两样本检验带来的假阳性率过高的情况, 根据需要可进一步采用多重检验方法进一步 判断样本组之间的差异, 也有多种分析方法来替代双样本 $t$ 检验法。

单因素方差分析

双样本平均数比较可使用成组的独立样本 $t$ 检验。对于多组样本的平均数检验, 零假设是这些多组样本来自同一总体, 各组的平均数相等, 备择假设是平均数并不都相等, 存在两对或多对样本的平均数不相等, 对于这种情况, 可以采用单因素多水平的方差分析。

方差分析的原理是检验多组样本的平均数的方差是否足够大,如果多组样本来自同一个总体, 不同组平均数的差异是随机产生的, 那么组间平方和接近于随机误差产生的平方和; 反之,组间平方和过大可以说明存在处理效应,即这几组样本来自不同总体。

单因素方差分析是针对 (单个因素 $a$ 个水平 (或处理)下得到的样本数据)进行的分析, 其线性统计模型是: $$ x_{i j}=\mu+\alpha_i+\varepsilon_{i j}, i=1,2, \cdots, a ; j=1,2, \cdots, n $$ 这里, $\mu$ 是总体平均数; $\alpha_i$ 是第 $i$ 个水平的相对平均数或相对处理效应; $\varepsilon_{i j}$ 是随机误差, 独立且服从正态分布。 $a$ 是水平的数目, $n$ 是每个水平下的样本数目, 在这里每个水平有相等 数目的样本数, 也称平衡设计。

方差分析的零假设是 $\alpha_i=0$ 或者 $\alpha_i$ 的方差为 0 , 变量 $x$ 的 平方和可以分解为处理平方和 (也称组间平方和) 和误差平方和, 分别除以各自的自由度 到处理均方和误差均方, 它们的比值服从 $F$ 分布, 可以使用 $F$ 检验进行判断。

方差分析采用函数 aov(), 它调用 lm() 函数进行线性模型的拟合, 与直接调用 lm()的 区别在于拟合结果的输出方式不同。因此, 也可以使用 lm()和anova() 函数联合进行方差分析。

aov() 函数的基本语法格式为： aov(formula, data = NULL, projections =FALSE, qr=TRUE, contrasts = NULL,...)

这里的参数 formula 是以公式的形式来表示模型, 模型模拟需要的数据在 data 中, projections、qr 都是逻辑值, 与返回的结果有关,一般使用缺省值。contrasts 是公式中用于因子水平对照(比较)方式，传递给 lm()函数。

例 不同均值的 4 组随机数的平均数差异分析, 使用 rnorm()函数随机产生 4 组服从不同参数正态分布的数据, 比较它们均值的差异。

data  <- round(c(rnorm(5),rnorm(5,2),rnorm(5,3),rnorm(5,1)),2)
V1 <- data.frame(data,FA=factor(c(rep(1,5),rep(2,5),rep(3,5),rep(4,5))))
V1.aov <- aov(data~FA,data=V1)
summary(V1.aov)

结果如下

            Df Sum Sq Mean Sq F value   Pr(>F)    
FA           3  29.36   9.785   9.318 0.000846 ***
Residuals   16  16.80   1.050                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(V1.aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = data ~ FA, data = V1)

FA
      diff        lwr        upr     p adj
2-1  1.666 -0.1883338  3.5203338 0.0862571
3-1  3.264  1.4096662  5.1183338 0.0006348
4-1  0.792 -1.0623338  2.6463338 0.6224982
3-2  1.598 -0.2563338  3.4523338 0.1043223
4-2 -0.874 -2.7283338  0.9803338 0.5473131
4-3 -2.472 -4.3263338 -0.6176662 0.0074814

在例中变量 V1 的数据是随机产生的, 有 4 组, 均服从不同均值的正态分布, 均值 别为 $0,2 、 3 、 1$, 它们的方差均为 1 。不同均值对应于因素 $\mathrm{FA}$ 的 4 个水平, 函数 aov() 可以判断这 4 个水平的变量均值是否相等, 分析得到 $\mathrm{F}$ 值为 $9.25, \mathrm{p}$ 值为 $0.00088$, 拒绝零假设, 除了F值和p值之外,还分别给出了处理(FA) 和误差(Residuals)项的自由度、平方和和均方值。

函数 TukeyHSD() 给出了 FA 因素 4 个水平之间的两两平均数差异分析的置信区间和 $\mathrm{p}$ 值, 可以得到 2-1,3-1,4-3 这三对水平之间存在平均数差异, 对照数据来源, 可以看到水 平 1 来自平均数为 0 的总体分布, 而水平 2 的总体平均数为 2 , 水平 3 的总体平均数为 3 , 水 平 4 的总体平均数为 $1 、 3-1$ 对之间的平均数差异相比较 $2-1,4-3$ 的更为显著, 与直观认识相符合。因素 $\mathrm{FA}$ 的水平数为 4 , 共有 6 对比较,这是个多重检验的问题,需要对 $\mathrm{p}$ 值进行 校正。图 2-1 给出了两两差异分析的结果,给出了平均数差异的区间, 2-1,3-1,4-3 的95\%置信区间不包含 0, 说明两者差异具有统计学显著性。

两因素和多因素方差分析

在考虑单因素的情况下, 再增加一个或多个因素, 与单独考虑单因素不同的是因素之间可能存在交互作用。在做方差分析时, 组间平方和(组间变异)包括了各个单因素主效应, 还应包括各因素各水平之间的相互作用效应, 对于计算来说就变得很复杂, 但原理与单因素分析是一样的, 计算各因素的组间变异 (处理间均方)和组内变异 (误差均方), 然后开展各因素和交互作用的处理间均方相对于误差均方的 $F$ 检验, 进而判断各因素的效应或交 互作用效应的统计学显著性。

例 判断新生儿体重与母亲的年龄、种族、吸烟状态、高血压史等因素是否有关。

本例选用 MASS 包中的 birthwt 数据集来示例多因素方差分析, 该数据集包括 189 个样本, 每个样本包括 10 个参数。

library(MASS)
data(birthwt)
birthwt.anova <- aov(bwt~race*smoke+ht+ptl*ui,data=birthwt)
summary(birthwt.anova)

结果

             Df   Sum Sq Mean Sq F value   Pr(>F)    
race          1  3790184 3790184   8.983 0.003107 ** 
smoke         1  7429202 7429202  17.608 4.25e-05 ***
ht            1  1857663 1857663   4.403 0.037266 *  
ptl           1  1057238 1057238   2.506 0.115175    
ui            1  6362550 6362550  15.080 0.000144 ***
race:smoke    1  1369773 1369773   3.247 0.073239 .  
ptl:ui        1  1735726 1735726   4.114 0.043998 *  
Residuals   181 76367320  421919                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

其中, $d$ 表示等级差异, $n$ 表示变量样本数。
R语言中计算相关系数的函数是 cor(), 其基本调用格式如下: cor(x, y=NULL, use = " everything ", method =c( "pearson", " kendall ", "spearman")) 其中, x表示矩阵或数据框; use 指定缺失数据的处理方式, 可选方式为 all,obs(不存在缺失 数据正常运行, 存在则报错)、everything(遇到缺失数据时, 相关系数的计算结果将被设为 missing)、complet.obs(行删 除)、pairwise, complete.obs (成对删除), 默认参数为 everything。若已过滤的数据可忽略此参数; method指定相关系数的类型, 默认为 pearson 相关系数。

函数 cor.test()对于两个测量变量之间的相关系数给出更多信息, 给出不相关的零假设并给出相关系数的置信区间。

例 计算随机变量之间的相关系数并检验, 随机产生符合正态分布的两个变量, 它们之间有相关性, 分析它们之间的相关关系。

x <- rnorm(20,4,1)
y <- 2*x+rnorm(20)
cor(x,y)

结果

0.949113956033621

摘自：