方差分析(analysis of variation,简写为ANOVA)又称变异数分析或F检验,用于两个及两个以上样本均值差别的显著性检验,从函数的形式看,方差分析和回归都是广义线性模型的特例,回归分析lm()也能作方差分析。其目的是推断两组或多组数据的总体均值是否相同,检验两个或多个样本均值的差异是否有统计学意义。方差分析的基本思路为:将试验数据的总变异分解为来源于不同因素的相应变异,并作出数量估计,从而明确各个变异因素在总变异中所占的重要程度;也就是将试验数据的总变异方差分解成各变因方差,并以其中的误差方差作为和其他变因方差比较的标准,以推断其它变因所引起的变异量是否真实的一种统计分析方法。把对试验结果发生影响和起作用的自变量称为因素(factor),即我们所要检验的对象。如果方差分析研究的是一个因素对于试验结果的影响和作用,就称为单因素方差分析。因素的不同选择方案称之为因素的水平(level of factor)或处理(treatment)。因素的水平实际上就是因素的取值或者是因素的分组。样本数据之间差异如果是由于抽样的随机性造成的,称之为随机误差;如果是由于因素水平本身不同引起的差异,称之为系统误差。

方差分析的基本前提各组的观察数据,要能够被看作是从服从正态分布的总体中随机抽得的样本。各组的观察数据,是从具有相同方差的总体中抽取得到的。观察值是相互独立的。 方差分析的原假设:H0 θ1=θ2=...=θk 即因素的不同水平对实验结果没有显著差异或影响。备择假设:不是所有的θi都相等(i=1,2,...,k),即因素的不同水平对实验结果有显著差异或影响。

aov()函数的语法为aov(formula,data=dataframe),formula可使用的特殊符号如下,其中y为因变量,A、B、C为自变量。

符号 用法
分隔符,左边为因变量,右边为自变量。例y~A+B+C
+ 分隔自变量
表示交互项,如y~A+B+A:B
* 表示所有可能的交互项,如y~A * B *C等价于y~A+B+C+A:B+A:C+B:C+A:B:C
^ 表示交互项达到的某个次数,如y~(A+B+C)^2等价于y~A+B+C+A:B+A:C+B:C
. 表示包含除因变量以外的所有变量。如y~.

常用的方差设计表达式如下,其中小写字母表示定量变量,大写字母表示组别因子,Subject是被试着的标识变量。

设计 表达式
单因素ANOVA Y~A
含但个协变量的单因素ANCOVA Y~x+A
双因素ANOVA Y~A*B
含两个协变量的双因素ANCOVA Y~x1+x2+A*B
随机化区组 y~B+A(B是区组因子)
单因素组内ANOVA y~A+Error(Subject/A)
含单个组内因子(w)和单个组间因子(b)的重复测量ANOVA Y~B*W+Error(Subject/W)

组别间观测数相等的设计均衡设计(balanced design),观测数不等的设计为非均衡设计(unbalanced design)。如果因子不止一个,且别是非平衡设计,或者存在协变量,表达式中的顺序会对结果造成影响。样本大小越不平衡,效应项的顺序对结果影响越大。通常,越基础的效应需要风在表达式的前面,如,先协变量,然后主效应,接着双因素的交互项,再接着是三因素的交互项。标准的anova()默认类型为序贯型,car包中的Anova()函数提供使用分层型和边界型(SAS和SPSS默认类型)的选项。

单因素方差分析(one-way ANOVA)

单因素方差分析是指对单因素试验结果进行分析,检验因素对试验结果有无显著性影响的方法。单因素方差分析是用来检验多个平均数之间的差异,从而确定因素对试验结果有无显著性影响的一种统计方法。对于完全随机设计试验且处理数大于2时可以用单因素方差分析(等于2 时用t检验)。离差平方和的分解公式为:SST(总和)=SSR(组间)+SSE(组内),F统计量为MSR/MSE,MSR=SSR/k-1,MSE=SSE/n-k。其中SST为总离差、SSR为组间平方和、SSE为组内平方和或残差平方和、MSR为组间均方差、MSE为组内均方差。

例 某医院欲研究A、B、C三种降血脂药物对家兔血清肾素血管紧张素转化酶(ACE)的影响,将家兔随机分为三组,均喂以高脂饮食,分别给予不同的降血脂药物。一定时间后测定家兔血清ACE浓度(u/ml),A组(45 44 43 47 48 44 46 44 40 45 42 40 43 46 47 45 46 45 43 44),B组(45 48 47 43 46 47 48 46 43 49 46 43 47 46 47 46 45 46 44 45 46 44 43 42 45),c组(47 48 45 46 46 44 45 48 49 50 49 48 47 44 45 46 45 43 44 45 46 43 42),问三组家兔血清ACE浓度是否相同?

a <- c(45, 44, 43, 47, 48, 44, 46, 44, 40, 45, 42, 40, 43, 46, 47, 45, 
    46, 45, 43, 44)
b <- c(45, 48, 47, 43, 46, 47, 48, 46, 43, 49, 46, 43, 47, 46, 47, 46, 
    45, 46, 44, 45, 46, 44, 43, 42, 45)
c <- c(47, 48, 45, 46, 46, 44, 45, 48, 49, 50, 49, 48, 47, 44, 45, 46, 
    45, 43, 44, 45, 46, 43, 42)
dfCRp <- data.frame(length = c(a, b, c), site = factor(c(rep("1", 20), 
    rep("2", 25), rep("3", 23))))
boxplot(length ~ site, data = dfCRp, xlab = "Sites", ylab = "Length")

plot.design(length ~ site, fun = mean, data = dfCRp, main = "Group means")

箱形图中可观察到不同的因素对于因变量的影响。

假设检验

方差分析需要一定的假设,即数据集应该符合正态和各组的方差相等,可以分别用shapiro.test和bartlett.test检验从P值观察到这两个假设是符合的。对于不符合假设的情况,我们就要用到非参数方法,例如Kruskal-Wallis秩和检验

shapiro.test(dfCRp$length)
## 
##  Shapiro-Wilk normality test
## 
## data:  dfCRp$length
## W = 0.97397, p-value = 0.1654
bartlett.test(length ~ site,data = dfCRp) #Fligner-Killeen(fligner.test()函数)和Brown-Forsythe检验(HH包中的hov()函数)也可以用来检验方差齐性
## 
##  Bartlett test of homogeneity of variances
## 
## data:  length by site
## Bartlett's K-squared = 0.76406, df = 2, p-value = 0.6825

正态性检验和方差齐性检验P值均大于0.05,可以认为数据满足正态性和方差齐性的要求。

oneway.test()和aov()函数进行方差分析

aovCRp =aov(length ~ site, data = dfCRp)
summary(aovCRp)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## site         2  26.29  13.146   3.244 0.0454 *
## Residuals   65 263.40   4.052                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#oneway.test(length ~ site, data=dfCRp, var.equal=TRUE),与aov()结果基本相同。
plotmeans(length ~ site,data =dfCRp ) #绘制有置信区间的组均值图

par(mfrow=c(2,2))
plot(aovCRp)

par(mfrow=c(1,1))

用aov函数建立单因子方差模型,从结果的P值可看到各组均值有显著不同。Sum Sq = deviance (within groups, and residual),总方差和(分别有groups和residual的),Mean Sq = variance (within groups, and residual),平均方差和(分别有groups和residual的)。单因子方差分析结果显示F value = 3.24 ,Pr(>F) = 0.045,因此拒绝原假设,即认为三组组家兔血清ACE浓度在统计学上有显著差异。

模型比较

(anovaCRp <- anova(lm(length ~ site, data=dfCRp)))
## Analysis of Variance Table
## 
## Response: length
##           Df  Sum Sq Mean Sq F value Pr(>F)  
## site       2  26.292 13.1462  3.2442 0.0454 *
## Residuals 65 263.399  4.0523                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(length ~ 1, data=dfCRp), lm(length ~ site, data=dfCRp))
## Analysis of Variance Table
## 
## Model 1: length ~ 1
## Model 2: length ~ site
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)  
## 1     67 289.69                             
## 2     65 263.40  2    26.293 3.2442 0.0454 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anovaCRp["Residuals", "Sum Sq"]
## [1] 263.3987

比较不含有自变量和含有一个自变量site模型,含有site变量的模型结果较好(残差的总方差和较小)。

效果大小(Effect size)

效果大小是指某个特定总体中的某种特殊的非零的数值,这个数值越大,就表明由研究者所处理的研究现象所造成的效果越大。效果大小本身可以被视为是一种参数:当原假设为真时,效果大小的值为零;当原假设为假时,效果大小为某种非零的值。因此,可以把效果大小视为某种与原假设分离程度的指标 。方差分析效果大小的含义也基本上与Z检验或t检验的效果大小的含义相同只不过它反映的是多组实验处理下不同组之间实验效果差异大小的指标。常用的指标如下η2=SSSStotal,f=η21η2,ω2=SSDFMSESStotal+MSE

dfSSb <- anovaCRp["site",        "Df"]
SSb   <- anovaCRp["site",        "Sum Sq"]
MSb   <- anovaCRp["site",        "Mean Sq"]
SSw   <- anovaCRp["Residuals", "Sum Sq"]
MSw   <- anovaCRp["Residuals", "Mean Sq"]

(etaSq <- SSb / (SSb + SSw))#   DescTools包中EtaSq(aovCRp, type=1)函数可以计算
## [1] 0.09076038
(omegaSq <- dfSSb * (MSb-MSw) / (SSb + SSw + MSw)) 
## [1] 0.06191765
(f <- sqrt(etaSq / (1-etaSq)))
## [1] 0.3159432

η2,ω2,f2值如上,如η2实验处理之后各组间平方和在总体平方和中所占的比重,η2越大反映实验效果大。一般η2大于0.14,就认为效果有大的效果。

多重比较

方差分析只告诉我们这三组之间是不同的,但没有告诉哪两组之间有明显差别,此时需要使用TukeyHSD函数进行均值的多重比较分析,从结果中观察到有一个两两比较是不显著的。 ####计划好的多重比较(Planned comparisons - a-priori) 在收集数据之前就已确定。它与实验目的有关,反映了实验者的意图。可以直接进行计划好的多重比较,不用考虑基本的“均值相等的 F-test”。

cntrMat <- rbind("a-c"          =c(1,0,-1),
                 "1/3*(a+b)-c"=c(1/3,1/3,-1),
                 "b-c"          =c(0,1,-1))
summary(glht(aovCRp, linfct=mcp(site=cntrMat), alternative="less"),
        test=adjusted("none"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = length ~ site, data = dfCRp)
## 
## Linear Hypotheses:
##                  Estimate Std. Error t value  Pr(<t)   
## a-c >= 0          -1.5196     0.6155  -2.469 0.00809 **
## 1/3*(a+b)-c >= 0  -1.1429     0.5331  -2.144 0.01790 * 
## b-c >= 0          -0.3896     0.5816  -0.670 0.25268   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
#pairwise.t.test(dfCRp$length, dfCRp$site, p.adjust.method="bonferroni") #结果与glht()函数类似。

依据事先实验的目的,进行多重比较,a组和c组,a、b组和c组的差异有显著意义。

非计划的多重比较(Planned comparisons - post-hoc)

在查看数据之后,并且“均值相等的 F-test” 结果显著情况下才进行。它用于探究研究者感兴趣但头脑中没有特定假设。

#ScheffeTest检验
ScheffeTest(aovCRp, which="site", contrasts=t(cntrMat))  #DescTools包
## 
##   Posthoc multiple comparisons of means : Scheffe Test 
##     95% family-wise confidence level
## 
## $site
##              diff     lwr.ci       upr.ci   pval    
## 1-3    -1.5195652  -3.061467   0.02233664 0.0543 .  
## 1,2-3 -15.9262319 -17.092478 -14.75998618 <2e-16 ***
## 2-3    -0.3895652  -1.846661   1.06753062 0.7997    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Tukey HSD检验
(tHSD <- TukeyHSD(aovCRp))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = length ~ site, data = dfCRp)
## 
## $site
##          diff         lwr      upr     p adj
## 2-1 1.1300000 -0.31850529 2.578505 0.1552673
## 3-1 1.5195652  0.04333482 2.995796 0.0422495
## 3-2 0.3895652 -1.00547115 1.784602 0.7817904
plot(tHSD)

置信区间包含0说明差异不显著。

multcomp包中glht()函数提供了多重均值更全面的方法,适用于线性模型和广义线性模型。下面的代码重现Tukey HSD检验。

tukey <- glht(aovCRp, linfct=mcp(site="Tukey"))
summary(tukey)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = length ~ site, data = dfCRp)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)  
## 2 - 1 == 0   1.1300     0.6039   1.871   0.1552  
## 3 - 1 == 0   1.5196     0.6155   2.469   0.0422 *
## 3 - 2 == 0   0.3896     0.5816   0.670   0.7817  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
plot(cld(tukey,level = .05),col="lightgrey") #cld()函数中level选项设置了使用显著水平0.05,即95%的置信区间

有相同字母的组(箱线图表示)说明均值差异不显著。

离群点检测

outlierTest(aovCRp)
## 
## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferonni p
## 9 -2.288156           0.025441           NA

离群点检测结果显示,数据中没有离群点(当p>1时产生NA)。

残差的相关检验

残差的正态性检验

Estud <- rstudent(aovCRp)
shapiro.test(Estud)
## 
##  Shapiro-Wilk normality test
## 
## data:  Estud
## W = 0.9876, p-value = 0.7381
qqnorm(Estud, pch=20, cex=2)
qqline(Estud, col="gray60", lwd=2)