生物统计(4)-单因素方差分析

方差分析的基本思想

在进行科学研究时，有时要按实验设计将所研究的对象分为多个处理组进行不同的处理，其中处理因素(treatment)至少有两个水平（level）。这类科研资料的统计分析，是通过所获得的样本信息来推断各处理组均数间的差别是否有统计学意义，即处理是否有影响。常用采用的分析方法就是方差分析（ANOVA，analysis of variance），这是由英国统计学家R.A.Fisher首创，以F命名，故方差分析又称为F检验。

设处理因素有g(g>= 2)个不同水平，实验对象随机分为g组，分别接受不同水平的干预，第i(i=1,2,...,g)组的样本含量为n_{i}，第i处理组的第j(j=1,2,…ni个观测值用Xij来表示，其计算结果可能可以整理成以下面的形式，如下所示：

方差分析的目的就是在 成立的条件下，通过分析各处理组均数之间的差别大小，推断g个总体均数之间有无差别，从面说明处理因素的效应是否存在。

记总均数为

各处理组均数为

总例数为其中，g为处理组数。

实验数据有三个不同的变异：
1. 总变异。全部观测值大小不同，这种变异称为总变异。总变异的大小可能用离均差平方和(sum of squares of deviations from eman,SS)来表示，即各观测值与总均数X差值的平方和，记为。公式略。

2. 组间变异。各处理组由于接受处理的水平不同，各组的样本均数也大小不等，这种变异称为组间变异，其大小用各组均数与总均数的离均差平方和表示，记为SS组间，计算公式略。

各组均数之间相关越悬殊，它们与总均数的差值越在在，就越大，反之就越小。反应了各组均数的变异，存在这种变异的原因有：①随机误差；②处理的不同水平可能对实验结果的影响。

3. 组内变异。在同一处理组中，虽然每个实验对象接受的处理相同，但观测值仍各不相同，这种变异称为组内变异（误差）。组内变异用组内各观测值与其所在组的均数的差值的平方和表示，记为，表示随机误差的影响。公式略。

总离均差平方和分解为组间离均差平方和与组内离均差平方和，就有了以下公式：

mark

各离均差平方和的自由度为：

mark

变异程序除与离均差平方和的大小有关外，还与其自由度有关，由于各部分自由度不相等，因此积分离均差平方和不能直接比较，须将各部分离均差平方和除以相应的自由度，其比值称为均方差，简称为均方（mean square, MS)。公式为：

mark

如果各组样本的总体均数相等()，即各处理组的样本来自相同的总体，无处理因素的作用，则组间变异同组内变异一样，只反映随机误差作用的大小，组间均方与组内无方的比值称为F统计量，如下所示：

mark

如果F值接近于1，就没有理由拒绝H0;反之，F值越大，拒绝H0的理由越充分，数理统计的理论证明，当H0成立时，F统计量服从F分布，方差分析是单侧F检验。

变异是方差分析的基本思想

上面的话可能不太好理解。现在用大白话来理解一下，例如我们要研究某个化合物是否有改善肥胖的效果，我们最初肯定是要做动物实验，动物实验的话，例如采用C57的小鼠，分为5组，第1组，给生理盐水，第2组，给减肥药（相当于阳性对照），第3组，给高剂量的化合物，第4组，给中剂量的化合物，第5组，给低剂量的化合物。分别喂一段时间后，我们发现小鼠的体重有所变化，这个变化由两部分构成，第一个就是外界的刺激因素导致的，第二种就是小鼠自身导致的。这种变化我们可以称为变异（variance），方差分析研究的本质就是这种变异（体重的变化，不是基因变异那种变异），方差分析的英语就是Analysis of variance，如果外界的刺激的因素导致的变异程度远远大于小鼠自身因素导致的变异，那么我们就可以认为，导致小鼠体重这种变异是由外界刺激引起的。

不过这样还有一个问题，因为数据越多，变异程度就越大，为了解决这个问题，就需要用变异除以自由度（例数-1），这样比较的就是平均的变异，因此方差分析中就出现了均方（MS）和组内均方的概念。组间均方/组内均方就是通常所说的F值，实际上代表了这样一个含义：如果组间变异远远大于组内变异，那么组间均方除以组内均方的值肯定很大，反之，这一值就会很小。但是，到底大到什么程度才认为有统计学意义呢，那就得根据F分布来判断。

方差分析的应用条件

多个样本均数比较的方差分析其应用条件为：①各样本是相互独立的随机样本，均来自**正态分布总体；②相互比较的各样本的总体方差相等，即具有等方差齐性。

R中的方差分析函数

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

mark

单因素方差分析

单因素方差分析（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为组内均方差。

案例分析

单因素方差分析例4-2：某医生为了研究一种降血脂新药的临床疗效，按统一纳入标准选择120名高血脂患者，采用完全随机设计方法将患者等分为4组（具体分组方法见例4-1），进行双盲试验。6周后测得低密度脂蛋白作为试验结果，见表4-3。问4个处理组患者的低密度脂蛋白含量总体均数有无差别?
数据如下所示：

计算过程如下所示：

1. 导入数据

anova1 <>'https://raw./20170505a/raw_data/master/data_szq_402.csv',sep=',')
library(reshape2)
anova1 <>


2.正态检验

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

shapiro.test(anova1$value)


结果如下所示：

P值大于0.05说明数据正态P值大于0.05说明数据正态

3. 方差齐性检验:

方差齐性检验就是检验各组样本所代表的总体方差是否一致的检验，两样本方差齐性检验使用Bartlett法，同样，它也适用于多样本的方差齐性检验，它是它要求所检验的样本总体符合正态分页，当不符合正态分布的时候，就不能使用，则要用Levene检验。Levene检验不受数据颁贩限制，是一种稳健的检验，因而被广泛地认为是一种标准的检验方差齐性的检验。

方差齐性通常用bartlett检验

bartlett.test(anova1$value~anova1$variable)


结果如下所示：

结果显示p值大于0.05，可认为方差齐性。

或者用levene检验：

library(car)
leveneTest(anova1$value~anova1$variable)


结果如下所示：

结果发现sig值大于0.05,表明符合方差齐性假设,可以进行进一步的参数检验。

4. 检验整体均值是否有差异

result <>value~variable,data=anova1)
summary(result)


结果如下所示：

其中p值小于0.001，因此各组之间存在显著性差异。另外，R给出的F值是24.88，而书中的例子是24.93，书中的值是由查F表得来的，是个范围，R中的是具体值。

另外也可以采用oneway.test()函数进行检验，如下所示：

> result2 <>TRUE)
> result2

    One-way analysis of means

data:  value and variable
F = 24.884, num df = 3, denom df = 116, p-value = 1.674e-12


5. 均数间的多重比较

方差分析得出总体之间有差异，要进一步知道哪两组之间有差异，就要使用均数间的多重比较，常用的比较方法有SNK检验（q检验），LSD检验，Bonferroni检验，Dunnett检验，TurkeyHSD检验。

现在计算各组之间的均数与SD

aggregate(anova1$value,by=list(anova1$variable),FUN=mean)
aggregate(anova1$value,by=list(anova1$variable),FUN=sd)


结果如下所示：

绘制箱形图可能观察到不同因素对于因变量的影响

plot.design(value ~ variable,data =anova1, main = 'Group means') 
# plot.design was inclued in packages 'graphics'


绘制有置信区间的组均值图

library(gplots)
plotmeans(value ~ variable,data =anova1 )


图片如下所示：

6. 组间均值的两两比较

通过方差分析后，如果整体有差异，则进一步进行两两比较，常用的方法有LSD，TukeyHSD，Scheffe检验，如下所示：

LSD检验

library(agricolae) #此包中有LSD检验
result <>'variable',p.adj='bonferroni')
result


结果如下所示：

注：R给出的F值是24.88，而书中的例子是24.93，书中的值是由查F表得来的，是个范围，R中的是具体值。
Bonferroni校正('bonferroni')，用于多重比较的p值校正，次数不多时适用。如果在同一数据集上同时检验n个独立的假设，那么用于每一假设的统计显著水平，应为仅检验一个假设时的显著水平的1/n。举个例子：如要在同一数据集上检验两个独立的假设，显著水平设为常见的0.05。此时用于检验该两个假设应使用更严格的0.025。即0.05* (1/2)。该方法是由Carlo Emilio Bonferroni发展的，因此称Bonferroni校正。在药理实验中，动物数通常不会太多，用Bonferroni校正的居多。

TukeyHSD检验

result2 <>value~variable,data=anova1)
result.tukey <>
result.tukey
plot(result.tukey)


结果如下所示：

图片结果：

glht()函数做Tukey检验

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

library(multcomp)
result4 <>
tukey4 <>'Tukey'))
summary(tukey4)


计算结果如下所示：

> summary(tukey4)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = value ~ variable, data = anova1)

Linear Hypotheses:
                      Estimate Std. Error t value Pr(>|t|)    
low - control == 0    -0.71500    0.16946  -4.219 0.000291 ***
middle - control == 0 -0.73233    0.16946  -4.322 0.000171 ***
high - control == 0   -1.46400    0.16946  -8.639  <>1e-04 ***
middle - low == 0     -0.01733    0.16946  -0.102 0.999615    
high - low == 0       -0.74900    0.16946  -4.420 0.000108 ***
high - middle == 0    -0.73167    0.16946  -4.318 0.000170 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)


绘制不同组的箱线图，如下所示：

plot(cld(tukey4,level = .05),col='lightgrey')
#cld()函数中level选项设置了使用显著水平0.05，即95%的置信区间


mark

7.残差分析

这一步做残差分析就是为了再次验证原始数据是否服从正态分布，如下所示：

residual <>
qqnorm(residual, pch=20, cex=2)
qqline(residual, col='gray60', lwd=2)


残差的QQ图如下所示（如果不了解QQ图，可以参考这篇文章《分位数及其应用》：

mark

shapiro检验，如下所示：

> shapiro.test(residual)

    Shapiro-Wilk normality test

data:  residual
W = 0.9884, p-value = 0.4026


p值的为0.4026，可以认为残差满足正态性。

绘制残差图，如下所示：

plot(residual ~ anova1$variable, main='各组的残差图')


mark

对残差进行方差齐性检验，如下所示：

> library(car)
> leveneTest(result4)
Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   3   1.493 0.2201
      116


P值大于0.05，可以认为残差满足方差齐性。

参考资料

1. 白话统计.冯国双

2. 医学统计学.第四版.孙振球