方差分析笔记

方差分析的基本思想及应用条件

方差分析的基本思想

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

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

方差分析的目的就是在$H_{0}:\mu_{1}=\mu_{2}\dots=\mu_$ 成立的条件下，通过分析各处理组均数之间的差别大小，推荐$g$个总体均数$\bar{X}_{i}$间有无差别，从面说明处理因素的效应是否存在。

记总均数为$\bar{X}=\sum_{i=1}^g\sum_{j=1}^{n_{j}}X_{ij}/N$，各处理组均数为$\bar{X}=\sum_{j=1}^{n_{i}}X_{ij}/n_{i}$,总例数为$N=n_{1}+n_{2}+\dots+n_{g}$，g为处理组数。

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

2. 组间变异。各处理组由于接受处理的水平不同，各组的样本均数$\bar{X}_{i}（i=1,2,\dots,g）$也大小不等，这种变异称为组间变异，其大小用各组均数$\bar{X_{i}}$与总均数$\bar{X}$的离均差平方和表示，记为$SS_{组间}$，计算公式略。

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

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

总离均差平方和分解为组间离均差平方和与组内离均差平方和，就有了以下公式：
$SS_{总}=SS_{组间}+SS_{组内}$

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

$\nu_{总}=N-1，\nu_{组间}=g-1,\nu_{组内}=N-g$

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

$MS_{组间}=\frac{SS_{组间}}{\nu_{组间}}$
$MS_{组内}=\frac{SS_{组内}}{\nu_{组内}}$

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

$F=\frac{MS_{组间}}{MS_{组内}}$

如果F值接近于1，就没有理由拒绝$H_{0}$;反之，F值越大，拒绝$H_{0}$的理由越充分，数理统计的理论证明，当$H_{0}$成立时，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为自变量：

符号	用法
～	分隔符，左边为因变量，右边为自变量。例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～.

单因素方差分析

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

ANOVA-案例A

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

计算过程如下所示：

1. 导入数据

anova1 <- read.csv("https://raw.githubusercontent.com/20170505a/raw_data/master/data_szq_402.csv",sep=",")
library(reshape2)
anova1 <- melt(anova1)

2.正态检验

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

group <- c("control","high","middle","control")

for(i in group){
  print(shapiro.test(anova1[anova1$variable==i,]$value))
}

结果如下所示：

	Shapiro-Wilk normality test

data:  anova1[anova1$variable == i, ]$value
W = 0.95229, p-value = 0.1947


	Shapiro-Wilk normality test

data:  anova1[anova1$variable == i, ]$value
W = 0.95833, p-value = 0.2806


	Shapiro-Wilk normality test

data:  anova1[anova1$variable == i, ]$value
W = 0.9445, p-value = 0.1202


	Shapiro-Wilk normality test

data:  anova1[anova1$variable == i, ]$value
W = 0.95229, p-value = 0.1947

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

3. 方差齐性检验:

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

方差齐性通常用bartlett检验

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

结果如下所示：

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

或者用levene检验：

Levene检验既可以用于正态分布的数据，也可以用于非正态分布的数据或分布不明的数据，具有比较稳健的特点，检验效果也比较理理想。

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

结果如下所示：

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

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

方差例如技术的目的就是要比较因素A的r水平上，试验结果是否有显著差异，以样本作为检验的标准，写出检验假设，如下所示：

$H_{0}:\alpha_{1}=\alpha_{2}=\cdots=\alpha_{r},$，$H_{1}:\alpha_{1},\alpha_{2},\cdots,\alpha_{r}$

如果拒绝原假设，则说明样本来自不同的正态总体，则由因素A的各个水平所造成的均值差异有统计意义；若不能拒绝原假设，说明样本来自相同的正态总体，因素的 不同水平之间无差异。在这个案例中，我们研究的就是不同剂量的药物对低密度脂蛋白的水平是否产生影响，我们可以使用3个函数来进行计算。

aov()

aov函数的使用如下所示：

aov(formula,data=NULL,projections=FALSE,qr=TRUE,contrasts=NULL,...)

其中formula表示主方差分析的公式，在单因素方差分析中，即为x~A；data表示方差分析的数据框；projections为逻辑值，表示是否返回预测结果；qr是逻辑值，表示是否返回QR分解结果，默认为TRUE；contrasts是公式中的一些因子的对比列表，计算过程如下所示：

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

结果如下所示：

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

oneway.test()

也可以采用oneway.test()函数进行检验，

oneway.test()函数的用法如下所示：

oneway.test(formual,data,subset,na.action,var.equal=FALSE)

若各水平下的总体方差相等，且参数var.equal=TRUE，那么等于同函数aov()所作的方差分析结果。在默认情况下，各水平的总体方差不相等(var.equal=FALSE)，则使用Welch的近似方法，在任何多个样本的情况下作两样本Welch检验，在满足方差齐发的条件下，使用oneway.test()的计算过程如下所示：

> result2 <- oneway.test(value~variable,data=anova1,var.equal=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

lm()

单因素方差分析的数学模型可以看作是一种特殊的线性模型，因此方差分析还可以通过线性模型函数lm()计算得到，再利用anova()提取其中的方差分析表，因此aov(formula)等价于anova(lm(formula))，如下所示：

> anova(lm(value~variable,data=anova1))
Analysis of Variance Table

Response: value
           Df Sum Sq Mean Sq F value    Pr(>F)    
variable    3 32.156 10.7187  24.884 1.674e-12 ***
Residuals 116 49.967  0.4308                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

5. 均数间的两两比较

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

5.1 常见的两两比较方法

LSD：相当于t检验，只不过它需要在方差分析一定要有统计学差异的情况下才用。所以LSD法并没有控制假阳性错误。一般情况下，如果你在设计初期就有很明确的目的，可以考虑这种方法，因为每一对比较都是有特定意义的，不用非得控制假阳性错误。

SNK法是先按多组均值大小排序，然后按一个有点类似于t检验的公式分别比较（不过误差计算不同）。比如a、b、c三组均值分别是a最小，c最大，b居中。那么比较时，很显然a和c差别最大，所以在最后的界值上做一些调整。例如，原来可能t检验的界值是2.83，只要t值大于2.83就算有差异，现在把界值调大，比如调到3.40，a和c的比较计算的t值，只有大于3.4才算有差异。这种调整的界值叫做q界值。这种方法对假阳性错误的控制依然不到位。所以如果你想控制假阳性，不要选这种方法。

Bonferroni法大概是最流行的一种方法，因为简单。它的思想是调整检验水准，根据比较的次数重新设定检验水准，然后根据P值做出结论。比如常规的检验水准是0.05，只要P小于0.05就认为有统计学差异。但是如果用Bonferroni法调整，则需要0.05除以比较次数，如比较6次，这时调整后的检验水准是0.05/6=0.0083，也就是说，P值小于0.0083才算有差异。这种方法其实更像是一种检验水准调整方法，简单易用，不仅可以用于均值的两两比较，也可用于率的两两比较。比如你要做三组的卡方检验，想进一步两两比较，就可以考虑用这种方式，即分别做两组的卡方，但是一共比较3次，所以P值小于0.0167才算有差异。这种方法的缺点是，当比较次数多时，结果过于保守，比如比较10次，那就需要P值小于0.005才算有差异，有时很难达到。

再说点其它国内课本上不大介绍但是国外介绍比较多的几种方法。

Tukey法，有时也叫Tukey HSD法（Honestly Significant Difference test），这种方法大概是统计学家最认可一致的方法了，绝大多数统计学家都相信，Tukey法的检验效率可能是最高的。Tukey法也是基于q检验，大概意思是先确定一个最大差异的临界值，然后分别对其中两组比较，看看哪两组差值大于这个界值，就算有差异。Tukey法是大多数统计学家首先推荐的两两比较方法，不过这种方法只适用于组间例数相等的情况。对于组间例数不等的时候，可用修正的Tukey法，也叫作Tukey-Kramer法。

现在计算各组之间的均数与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 )

图片如下所示：

单因素方差分析是从总体的角度上说明各效应的均值之间存在着显著差异，但具体在哪些水平下的均值存在较大差异，还不清楚，因此我们要对每一对样本均值进行一一比较，既要进行组间均值的两两比较。最常用，也是最直接的方法就是多重t检验方法，其对因子A每个水平下的数据两两比较均值进行t检验，但估计方差时用的是全部数据的误差均方和$MS_{E}$，检验的假设为：

$H_{0}:\mu_{i}=\nu_{j},i\neq j,i,j=1,2,\cdots,r$

构建t检验统计量，如下所示：

在原假设成立的情况下，该统计量服从自由度为n-r的t分布，所以检验拒绝为$|T_{ij}>t_{\alpha/2}(n-r)$。

多重t检验的思路简单，但是这样的方法不能直接使用，这是因为在样本之间多次重复使用t检验会大大增加犯第一类错误的概率。例如因子A有三个水平，需要比较3次(n=k(k-1)/2)，显著性水平为0.05，回此每次比较犯第一类错误的概率就是0.05，那么截取次比较同时进行，犯第一次错误的总概率为$1-(1-\alpha)^n=1-0.95^3=0.1426$，这个数字就比较大了，因此直接使用多重t检验的结果不一定是可靠的，统计学家们提供了很多方法对P值进行修正，R中的P值修正函数是p.ajdust()，其调用格式如下所示：

p.adjust(p,method=p.ajdust.methods,n=length(p))

在R中，输入如下指令可以看到调整方法的列表，如下所示：

> p.adjust.methods
[1] "holm"       "hochberg"   "hommel"     "bonferroni" "BH"         "BY"         "fdr"       
[8] "none"

当多重检验次数较多时，bonferroni修正方法的效果比较好，思路也很简单：如果在同一个数据集上同时进行n个独立的假设检验，那么用于每一个假设的统计显著水平，应为仅检验一个假设时显著水平的1/n，即$\alpha’=\alpha/n$，其中n为多重t检验的次数。

在R中计算均值的多重t检验时，用函数pairwise.t.test()返回多重比较的修正P值，其格式如下所示：

pairwise.t.test(x, g, p.adjust.method = p.adjust.methods,
                pool.sd = !paired, paired = FALSE,
                alternative = c("two.sided", "less", "greater"),
                ...)

其中参数x是响应向量，g是因子向量，p.adjust.methods给出了P值的修正方法，默认按Holm方法修正，若不作任何调整则设为none；paired为逻辑值，指示是否要配对t检验；alternative指明检验的方向问题。现在对前面的结果进行多重t检验，使用bonferroni修正方法，如下所示：

> pairwise.t.test(anova1$value,anova1$variable,p.adjust.method = "bonferroni")

	Pairwise comparisons using t tests with pooled SD 

data:  anova1$value and anova1$variable 

       control low     middle 
low    0.00029 -       -      
middle 0.00020 1.00000 -      
high   2.1e-13 0.00013 0.00020

P value adjustment method: bonferroni

经过修正后的p值会比原理增大很多，这在一定程度上克服了多重t检验增加犯第一类错误概率的缺点，从上面结果来看，除了middle与low组外，其余组之间的p值都很小，说明这几个样本之间有差异。

再看其它的两两比较方法，这些方法有LSD，TukeyHSD，Scheffe检验，

LSD检验

library(agricolae) #此包中有LSD检验
result <- LSD.test(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 <- aov(value~variable,data=anova1)
result.tukey <- TukeyHSD(result2)
result.tukey
plot(result.tukey)

结果如下所示：

图片结果：

glht()函数做Tukey检验

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

library(multcomp)
result4 <- aov(value ~ variable, data = anova1)
tukey4 <- glht(result4, linfct=mcp(variable="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%的置信区间

7.残差分析

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

residual <- rstudent(result4)
qqnorm(residual, pch=20, cex=2)
qqline(residual, col="gray60", lwd=2)

残差的QQ图如下所示：

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="各组的残差图")

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

> 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，可以认为残差满足方差齐性。

ANOVA-案例B

再看一个案例，进行一次完整的ANOVA分析，如下所示：

某单位研究不同药物对小白鼠的镇咳作用，抽取40只小白鼠随机分配到各个药物组中，实验时先用0.2mL NH4OH对小白鼠喷雾，测定其发生咳嗽时间。以给药前后发生咳嗽时间的差值衡量不同药物的镇咳作用，结果见表5.1.试比较3种药物的平均推迟咳嗽时间的差异有无差异。（王炳顺. 医学统计学及SAS应用[M]. 上海交通大学出版社, 2007.）

计算过程如下所示：

anova2 <- read.csv("https://raw.githubusercontent.com/20170505a/raw_data/master/data_anova_01.csv",sep=",")
# input raw data 
# 王炳顺. 医学统计学及SAS应用[M]. 上海交通大学出版社, 2007.例5.1

anova2$group <- as.factor(anova2$group)
library(car)
leveneTest(anova2$time~anova2$group)
result <- aov(anova2$time~anova2$group)
summary(result)
compareresult <- TukeyHSD(result)
compareresult
plot(compareresult)

结果如下所示：

> for(i in 1:3){
+   result <- shapiro.test(filter(anova2,group==i)$time)
+   print(result)
+ }

	Shapiro-Wilk normality test

data:  filter(anova2, group == i)$time
W = 0.9778, p-value = 0.9523


	Shapiro-Wilk normality test

data:  filter(anova2, group == i)$time
W = 0.98349, p-value = 0.9878


	Shapiro-Wilk normality test

data:  filter(anova2, group == i)$time
W = 0.94005, p-value = 0.5536

> 
> # test for homogeneity var
> leveneTest(anova2$time~anova2$group)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  2  1.4002 0.2593
      37               
> 
> ## test for mean of all samples
> result <- aov(anova2$time~anova2$group)
> summary(result)
             Df Sum Sq Mean Sq F value Pr(>F)  
anova2$group  2   5062  2531.2   3.485 0.0411 *
Residuals    37  26877   726.4                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> ## GLM for mean of all samples
> anova(lm(anova2$time~anova2$group))
Analysis of Variance Table

Response: anova2$time
             Df  Sum Sq Mean Sq F value  Pr(>F)  
anova2$group  2  5062.5 2531.23  3.4845 0.04107 *
Residuals    37 26877.4  726.42                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> ## test for each pair of samples
> ## TukeyHSD
> TukeyHSD(result)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = anova2$time ~ anova2$group)

$`anova2$group`
        diff        lwr      upr     p adj
2-1 12.33333 -11.694604 36.36127 0.4302024
3-1 29.03333   2.169283 55.89738 0.0317315
3-2 16.70000 -10.164050 43.56405 0.2944100

> pairwise.t.test(anova2$time,anova2$group,p.adjust.method = "bonferroni")

	Pairwise comparisons using t tests with pooled SD 

data:  anova2$time and anova2$group 

  1     2    
2 0.654 -    
3 0.036 0.413

P value adjustment method: bonferroni

重复测量方差分析

未设立平等对照的前后测量

重复测量资料(repeated measurement data)最常见的情况是同一个试验对象前后有2次测量结果的试验结果，也称作前后测量设计(premeasure-postmeasure design)，如下所示：

… …

虽然这个表和配对设计的t检验结果表达完全相同，但却是属于两种不同类型的设计，其区别在于

1. 配对设计中同一对子的两个实验单位可以随机分配处理，两个实验单位同期观察试验结果，哦可以比较处理组间差别。前后测量设计不能同期观察试验结果，虽然可以在前后测量之间安排处理，但本质上比较的是前后差别，推论处理是否有效是有条件的，即假定测量时间对观察结果没有影响。
2. 配对t检验要求同一对子的两个实验单位的观察结果分别与差值相互独立，差值服从正态分布。而胶后测量设计前后两次观察结果通常与差值不独立，大多数情况下第一次观察结果与差值存在负相关关系，例如上表，治疗前舒张压与差值的相关系数为-0.602。
3. 配对设计用平均差值推论处理的作用，前后测量测序除了分析平均差值外，还可以进行相关回归分析，由上表计算，治疗前后舒张压的相关系数为0.963，P<0.01，用治疗前舒张压(X)推论治疗后舒张压(Y)的回归方向为$\hat{Y}=-49.534+1.266X$，截矩检验P=0.014，回归系数检验P<0.01。

设立平等对照的前后测量

第1张表中的高血压患者治疗后的舒张压平均下降了16mmHg，虽然经配对t检验，t=16.18，P<0.01，这也未必说明治疗有效，因为住院休息，环境和情绪的改变同样可以使血压恢复平衡。因此，确定疗效的前后测量设计必须增加平行对照，例如将20位轻度高血压患者随机分配到处理组和对照组，试验结果见表12-2.如下所示：

比较这张表中处理组与对照组之间的差别，似乎可以用前后测量的差值(d)做两组均数比较的t检验，但上表中处理组与对照组的差值d方差不齐（F值为6.58，P<0.01），不符合两无数比较t检验的前提条件，更为重要的是，如果处理组与对照组基线值（治疗前的舒张压）不同，那么差值不能作为判别组间差别的依据。

重复测试设计

当前后测量设计的重复测量次数m>=3时，称为重复测试设计或重复测量数据，重复测试数据类似于随机区间设计数据，如下表所示，而且同样可以计算出随机区组设计的方差分析表：

案例分析B

下表是随机区组数据的方差分析表：

但重复测量设计与随机区级设计有区别，它们的区别在于：

（一）重复测量设计中的“处理”是在区组（受试者）间随机分配，区组内的各时间点是固定的，不能随机分配，如下表所示：

A、B两种处理随机分配给各个患者后，每个患者测量的时间是相同的，随机区组设计要求每个区组内实验单位彼此独立，处理只能在区组内随机分配，每个实验单位接受的处理是不相同的，如下所示（随机区间设计）：

（二）重复测量设计区组内实验单位彼此不独立，如在表12-3中，每个受试者血糖浓度的个体特征是用4个时间点的测量值刻画的，如下图所示：

同一受试者的血样重复测量结果是高度相关的，相关系数如下所示：

更准确地说，表12-3重复测量数据用随机区组方差分析比较处理组间差异，前提条件是要满足“球对称(sphericity)”假设，即重复测量误差的协方差矩阵经正交对比变换后，与单位矩阵$I_{4\times4}$成比较，表12-3数据Mauchly“球对称”检验结果见表12-7，其中$\chi^2=15.44,\nu=5,P=0.01$，拒绝“球对称”假设，即表12-3的数据不满足随机区组方差分析（表12-4）“球对称”的条件，如下所示：

此时，就需要对组内效应的F界值进行校正，校正的方法是用“球对称”系数$\epsilon$乘以组内效应F界值的自由度$\nu_{1}$和$\nu_{2}$，得到$\tau_{1}=\nu_{1}\epsilon,\tau_{2}=\nu_{2}\epsilon$，用$F_{\alpha,(\tau_{1},\tau_{2})}$作为检验界值，“球对称”系数的常用估计方法有Greenhouse-Geisser、Huynh-Feldt和Lower-bound三种方法，例如，表12-4组内效应的F界值为$F_{0.05,(3,21)}=3.07$，Greenhouse-Geisser的校正系数0.536校正后的F界值为$F_{0.05,(1.6,11.25)}\approx F_{0.05,(1,11)}=4.84$，大于未校正的界值3.07，也就是说当重复测量数据不满足”球对称“假设时，采用随机区组设计方差分析时，增大了I型错误（无差别判断为有差别）的概率，当样本含量较小时，Mauchly”球对称“检验交通较低，易发生II型错误（不满足”球对称“判断为满足），此时，对处理组内效应F界值的自由度$\nu_{1}$和$\nu_{2}$仍要进行校正，以得到可信的统计结论。

重复测量数据的两因素两水平分析

两因素离均差平方和分解

将重复测量测序的干预因素作为A因素，共两个水平，1水平为”对照“，2水平为”干预“；前后两次测量时间作为B因素，共两个水平，1水平为”第1次测量时间“，2水平为”第2次测量时间“，这样共有(a1b1)、(a1b2)、(a2b1)、(a2b2)4个”处理“组，各组观察值的合计分别为T1、T2、T3、T4表示，A因素两水平的小计分别用A1、A2表示，B因素两水平的小计分别用B1、B2表示，处理组和对照组的例数相等，即n=n1=n2，SS处理可分解为A因素的离均差平方和$SS_{A}$，B因素的离均差平方和$SS_{B}$，AB交互作用的离均差平方和$SS_{AB}$三部分，见表12-8，如下所示：

将表12-2的数据按照表12-8进行两因素离均差平方和进行分解，如下所示：

主效应与交互作用的离均差平方和分解：

组间、组内离均差平方和分解与方差分析

重复测量数据的差异由两部分组成。一是观察对象的个体差异，其离均差平方和记作SS组间，二是每个观察对象前后两次观察之间的差异，其离均差平方和记作SS组内，设两组观察对象都为n，每个观察对象前后两次观察的合计值为Mj，SS组间和SS组内的计算方法见表12-9，如下所示：

SS组间、SS组内分解完毕后，就可以用方差分析对A因素主效应、B因素主效应和AB因素交互作用进行假设检验了，由于A因素（干预分组）的作用于观察以刚，作用效应在SS组间，方差分析按表12-10进行，B因素主效应（测量前后）和AB交互作用的效应在SS组内,方差分析按表12-11进行，如下所示：

现在根据表12-2中的数据，对处理组与对照组、治疗前后舒张压的差别进行统计分析。

1.计算$SS_{组间}$、$SS_{组内}$

表12-2共有20名患者，每个患者舒张压的合计为M1=130+114=244，M2=124+110=234，…，M19=120+124=244，M20=134+128=262，由于前面已经计算出C=580328.1，按表12-9中的公式计算，如下所示：

2.测量前后比较与交互作用分析

由于已经计算出了$SS_{B}=1020.1,SS_{AB}=348.1$，现在SS组内=1702，代入表12-11，得到方差分析表，如下所示：

3.处理组与对照组比较

前面已经计算出$SS_{A}=202.5$，现在SS组间=2517.9，代入表12-10，得到方差分析表，如下所示：

4.结论

①由表12-12的方差分析表可知，治疗前后舒张压的交效应有差别(P<0.01)，由于表12-2计算可知，治疗前的平均舒张压为(126.2+124.8) 2="125.5mmHg，治疗后的平均舒张压为(110.2+120.6)/2=115.4mmHg，治疗前后舒张压主效应的估计值为115.4-125.5=-9.6mmHg。①由表12-13方差分析表，不考虑测量时间，处理组与对照组舒张压的主效应未见差别(P">0.05)，由原始数据过计算，处理组的平均舒张压为(126.2+110.2)/2=118.2mmHg，对照组的平均舒张压为(124.8+120.6)/2=122.7mmHg，处理组与对照组舒张压主效应的估计值为118.2-122.7=-4.5mmHg。③由表12-12方差分析表可知，测量前后与处理存在交互作用(P<0.01)，即处理组和对照组治疗前后的舒张压的变化幅度不同。由单独效应、主效应和交互作用的概念可知，若存在交互作用，单独分析主效应的意义不大，需要爱一分析各因素的单位效应，按表12-2计算处理组和对照组的单独效应，处理组和对照组治疗前的基线数据没有统计学差异，如下所示：

但处理组和舒张压平均下降16.0mmmHg，对照组平均仅下降4.2mmHg，AB交互作用F=18.8，P<0.01，说明处理组的降压效应优于对照组。

随机区组设计

随机区组设计(randomized block design)又称为配伍组设计，是配对设计的扩展。具体做法是：先按影响实验结果的非处理因素（如性别、体重、年龄、职业、病情、病程等）将实验对象配成区组(block)，再分别将各区组内的实验对象随机分配到各处理组或对照组。与完全随机设计相比，随机区组设计的特点是随机同分配的次数要重复多次，每次随机分配都对同一个区组内的实验对象进行，且各个处理组的实验对象数量相同，区组内均衡。在进行统计分析时，将区组变异离均差平方和从完全随机设计的组内离均差平方和中分离出来，从而减小组内离均差平方和（误差平方和），提高统计检验效率。若将区组作为另一处理因素的不同水平，随机区组设计等同于无重复观察的两因素设计。

现在看一下随机区组设计的案例。

如何按随机区组设计，分配5个区组的15只小白鼠接受甲、乙、丙三种抗癌药物？

方法为先将小白鼠的体重从轻到重编号，体重相近的3只小白鼠配成一个区组（下表中第一行和第二行），然后对这些小鼠分配随机数，在每个区组内将随机数按大小排序（下表第四行），各区组中序号为1的接受甲药、序号为2的接受乙药、序号为3的接受丙药，如下所示：

随机区组设计资料在进行统计分析时，需要根据数据特征选择方法，对于正态分布且方差齐同的资料，采用双向分类的方差分析(two-way classification ANOVA)或配对t检验(g=2)，当不满足方差分析或t检验条件时，可以进行变量变换后采用双向分析的方差分析或采用Fredman M检验。

重复测量数据的两因素多水平分析

大多数医学实验都有重复测量记录，如果统计分析时只分析最后一次的测量结果，那就会丧失很多“过程“信息，如果通过记录不同时间点的重复观察数据了解测量指标的时间变化趋势、药物在血中浓度变化的时间分析特征、不同时间的治疗效果等。在临床上，各个单个时间的测量值或许都在正常参考值范围以内（或以外），但如果发现该测量值持续上升（或下降），也说明患者的状态处理变化中，在统计分析时，保留”处理“前的基线(baseline)，如前面所述的患者治疗前收缩压，可以评价随机分级的均衡性，也能够提高统计分析的效率。

实验设计方法

一般来说，任何实验设计都可以采用重复测量计算，即在实验过得中定期记录观察结果。所谓重复测量数据的两因素多水平设计，两因素指干预（A因素）和测量时间（B因素），多水平指干预（A因素）有g（大于等于2）个水平，测量时间（B因素）有m（大于等于2）个水平（时间点），即每个观察对象有m个重复测量数据，随机化分组是将作用于观察对象的g个干预随机分配，前面提到的重复测量数据的两因素两水平是两因素多水平的设计的特征，即g=m=2.

试验结果的方差分析

重复测量数据的统计分析有许多统计方法供选用，可以用单因素方差分析，也可以用多变量方差分析(MANOVA)，其中单因素方差分析最简单。设有观察对象随机等分成g个干预组，每组倒数为n，重复测量次数为m，每个观察对象测量值合计为$M_{j}$，g个干预组，每组的测量值合计为$A_{i}$。多个干预组比较的方差分析用表12-14，设$B_{j}$表示第j个时间的测量值合计，$T_{ij}$表示第i个干预，第j个时间的点测量值合计，多个时间点测量前后与交互作用的方差分析见表12-15，如下所示：

案例分析B

将手术要求基本相同的15名患者随机分为3组，在手术过程中分别采用A、B、C三种麻醉诱导方法，在T0（诱导前）、T1、T2、T3、T4五个时间点测量患者的收缩压，数据记录见表12-16，试进行方差分析，如下所示：

1.分解$SS_{组间}，SS_{组内}$

g=3,m=5,n=5,15名患者的合计分别为：

2.分解$SS_{A}，SS_{B},SS_{AB}$

分组计算不同麻醉诱导、不同时间点患者的收缩压的合计值($T_{ij}$)，见表12-17，如下所示：

3.方差分析表

按表12-14、12-15列出方差分析表。

4.F值的界值校正

按表12-15的Lower-Bound自由度校正方法，$F_{B}$的校正自由度$\tau_{1}=1,\tau_{2}=g(n-1)=3 \times(5-1)=12$，查F界值表，$F_{B}$的校正界值为$F_{0.01,(1,12)=9.33}$（校正前为$F_{0.01,(4,48)=3.74}$）。$F_{AB}$的校正自由度$\tau_{1}=g-1=2,\tau_{2}=g(n-1)=3 \times(5-1)=12$，查F界值表，$F_{AB}$的校正界值为$F_{0.01,(2,12)=6.93}$（校正前为$F_{0.01,(8,48)=2.90}$），校正后的自由度和P值见下表：

5.结论

不同麻醉诱导方法存在组间差别，见表12-18：

患者的收缩压在不同的诱导方法下，不同诱导时相变化的趋势不同，见表12-19:

其中A组不同诱导时间收缩压相对稳定，见表12-20，如下所示：

重复测量数据的多重比较

重复测量数据的多重比较是一个比较复杂的问题，对不同研究有不同的要求，例如前面提到的，A、B、C三种麻醉诱导方法存在差别，每组患者不同诱导时相收缩压的趋势不同，如下所示：

多重比较可能包括3种情况，一是A、B、C三组收缩压均数差别的两两比较，需要重复检验3次（A与B，A与C，B与C）；二是各组T0，T1，T2，T3和T4这5个不同时相收缩压均数变化曲线是否有固定趋势（$H_{0}:\mu_{i}=\mu,i=1,2,3,4,5$)其中，$\mu$是收缩压均数，如果拒绝H0，对于5次重复测量结果，曲线经下次多项式(polynomial)变换，需要检验下次多项式1次(linear)、2次(quadratic)、3次(cubic)、4次(order 4)是否为0，上图中的三条曲线共需要重复检验12次；三是各组T0，T1，T2，T3，T4这5个不同时相收缩压均数的两两比较，每组需要比较10次，3组共要做30次重复检验 。因此，前面提到的例题的多重检比较需要做2+12+30=44次。下面我们列出了SPSS重复测量数据的多重比较的计算结果。

组间差别多重比较

在表12-18基础上进行两两比较，误差均方为$\sqrt{946.48/4}=15.38$，由表12-20可知，A、B、C三组收缩压均数分别为119.68、124.48、128.20，两两比较可以选择LSD、SNK、Bonferroni等多重比较方法，下图为LSD两两比较结果，其中A与C的差别有统计学意义，如下所示：

时间趋势比较

由表12-20可知，A组5个时相的收缩压分别为121.00，112.40，118.40，125.80，120.80，在SPSS中选择Polynomial正交多项式。表12-23就列出了A组收缩压均数5个时相的时间趋势检验结果，正交多项式1次项、3次项有统计学意义(P<0.05)，拒绝$H_{0}:\mu_{i}\neq \mu$。B组和C组的时间趋势检验结果见表12-24、表12-25，正交多项式1次项、2次项、4次项均有统计学意义，拒绝$H_{0}:\mu_{i}\neq \mu$。如下所示：

时间点的多重比较

一些医学研究关心同一个处理重复测量数据时间点的两两差别，如例12-3中A组5个时相收缩压均数之间的差别。由于当时间点较多时，重复检验次数将难以承受，如果假定5个处理组，10个时间点，重复检验次数就是5x45=225次。因此，时间点多重比较通常采用事前检验(prior tests)的策略。事前检验与事后检验(ost hoc tests)的不同之处在于，事前检验的检验次数是在试验开始前就确定了，无论方差分析结论如何都执行重复检验，比较的时间点按专业意义设定，如试验后各时间点均数与基线比较，重复检验次数比较少。事后检验则是方差分析结果拒绝$H_{0}$才进行两两比较，通常所有的时间上炽都要比较、重复检验次数较多。

按照事前检验策略，例12-3中A组5个时间收缩压均数差别的两两比较，只比较麻醉诱导后各时相T1，T2，T3，T4与诱导前T0的差别，只需要检验4次，采用Bonferroni多重比较方法，设定α=0.05，检验的名义水准为$a’=0.05/4=0.0125$，下表为SPSS的检验结果，可以发现，A组麻醉诱导后的4个时相与诱导前比较的重复检验，仅T1收缩压降低(P<0.05)。B组诱导后T2收缩压降低，T3、T4收缩压增高(P<0.05)，C组诱导后T2收缩压降低，T3收缩压升高(P<0.05)，因此A组麻醉诱导后收缩压的波动小于B组和C组，如下所示：

例11-1：2×2析因设计

将20只家兔随机等分4组，每组5只，进行神经损伤后的缝合试验。处理由两个因素组合而成，A因素为缝合方法，有两水平，一水平为外膜缝合，另一水平为束膜缝合，B因素为缝合后的时间，有两水平，一水平为缝合后1月，另一水平为缝合后2月。试验结果为家兔神经缝合后的轴突通过率(%)（注：测量指标，视为计量资料），见表11-1。欲比较不同缝合方法及缝合后时间对轴突通过率的影响，试做析因方差分析。

data1101 <- c(10,10,40,50,10,30,30,70,60,30,10,20,30,50,30,50,50,70,60,30)
factor_a <- gl(2,10) # a取两个水平，每个水平数量是10个，排列顺序是1...2...
factor_b <- rep(gl(2,5),2) # b取两个水平，每个水平是10个，排列顺序是1..2..1..2
data1101.aov <- aov(data1101~factor_a+factor_b+factor_a*factor_b) # factor_a*factor_b表示a与b的交互作用
summary(data1101.aov)

从上述分析结果可以看出，只有B因素的主效应达到了p<0.05，而a因素的主效应p>0.05，A与B的交互作用的p>0.05，结论为：沿不能认为两咱缝合方法对神经轴突通过率有影响，但缝合后2个月与1个月相比，神经轴突通过率提高了。

例11-2：I×J两因素析因分析

观察A,B两种镇痛药物联合运用在产妇分娩时的镇痛效果。A药取3个剂量：1.0mg，2.5mg，5.0mg；B药也取3个剂量：5 ，15 ，30 。共9个处理组。将27名产妇随机等分为9组，每组3名产妇，记录每名产妇分娩时的镇痛时间，结果见表11-7。试分析A,B两药联合运用的镇痛效果。

data1102 <- c(105,80,65,75,115,80,85,120,125,115,105,80,125,130,90,65,120,100,75,95,85,135,120,150,180,190,160)
drug_a <- factor(rep(gl(3,3),3))
drug_b <- factor(gl(3,9))
data1102.aov <- aov(data1102~drug_a+drug_b+drug_a*drug_b)
summary(data1102.aov)

例11-3：I×J×K三因素析因分析

用5×2×2析因设计研究5种类型的军装在两种环境、两种活动状态下的散热效果，将100名受试者随机等分20组，观察指标是受试者的主观热感觉(从“冷”到“热”按等级评分)，结果见表11-11。试进行方差分析。

data1103 <- c(0.25,0.3,0.75,0.2,-0.1,-0.25,0.1,-0.5,-1,0,1.25,0.5,0.6,0.85,2.5,-0.75,-0.35,0.4,-0.5,0.1,0.4,0.05,-0.2,0.9,-0.1,4.75,4.6,4.55,4.25,4.72,3.45,4.8,3.5,3.1,4.3,4,4,4.25,4,4.1,4.85,5.2,4.1,5,4.8,4.55,4.3,4.4,4.2,3.6,0.5,1.5,0.75,-0.75,1.75,2.1,1.5,2.65,0.9,2.4,2.75,1.25,3,0.95,1.75,1,1.37,0.05,0.62,3.05,2.35,2.55,1.17,1.05,2.75,3.75,4,4.1,3.27,4.8,4,4.05,5,4.25,4.02,4,4.15,4.2,4,4.15,4.25,4.1,4.15,4.25,4.75,4.6,4.25,4.17,4.25,4.8)

a <- factor(rep(1:5,20))
b <- factor(rep(1:2,each=50))
c <- factor(rep(c(1,2,1,2),each=25))

data1103.aov <- aov(data1103~a+b+c+a*b+a*c+b*c+a*b*c)
summary(data1103.aov)

例11-4

研究雌螺产卵的最优条件，在20cm2的泥盒里饲养同龄雌螺10只，试验条件有4个因素（见表11-15），每个因素2个水平。试在考虑温度与含氧量对雌螺产卵有交互作用的情况下安排正交试验。

data1104 <- c(86, 95,91,94,91,96,83,88)
A <- factor(rep(1:2,each=4))
B <- factor(gl(2,2,8))
AB <- factor(c(1,1,2,2,2,2,1,1))
C <- factor(rep(1:2,4))
D <- factor(c(1,2,2,1,2,1,1,2))
data1104.aov <- aov(data1104~A+B+AB+C+D)
summary(data1104.aov)

例11-5

研究高频呼吸机A，B，C，D，E五个参数对通气量的影响，每个参数有高、低两个水平，试按正交设计安排试验。

data1105 <- c(16.26, 19.38,23.6,28.43,20.28,34.88,49.1,47.44,18.32,24.85,39.45,32.08,45.5,50.3,55.26,66.64)
A <- factor(rep(1:2,each=8)) # 产生2个水平的因子，每个因子8个
B <- factor(gl(2,4,16)) # 产生2个水平的因子，每个因子4个，并且把向量B补充到16个元素
C <- factor(gl(2,2,16))
D <- factor(rep(1:2,8)) # 产生2个水平的因子，2个水平的因子作为整体，复制8遍
E <- factor(c(1,2,2,1,2,1,1,2,2,1, 1,2,1,2,2,1))

data1105.aov <- aov(data1105~A+B+C+D+E)
summary(data1103.aov)

例11-6

试验甲、乙、丙三种催化剂在不同温度下对某化合物的转化作用。由于各催化剂所要求的温度范围不同，将催化剂作为一级实验因素（I=3），温度作为二级实验因素（J=3），采用嵌套设计，每个处理重复2次（n=2），试验结果见表11-25。试做方差分析。

catalyst <- factor(gl(3,6,18))
temperature <- factor(rep(c(70,80,90,55,65,75,90,95,100),each=2))
data1106 <- c(82, 84,91,88,85,83,65,61,62,59,56,60,71,67,75,78,85,89)

data1106.aov <- aov(data1106~catalyst+temperature)
summary(data1106.aov)

例11-7

试验一种全身注射抗毒素对皮肤损伤的保护作用试验，将10只家兔随机等分两组，一组注射抗毒素，一组注射生理盐水作对照。分组后，每只家兔取甲、乙两部位，分别随机分配注射低浓度毒素和高浓度毒素，观察指标为皮肤受损直径(mm)，试验结果见表11-31。试做方差分析。

A <- factor(gl(2,5,20))
B <- factor(rep(1:2,each=10))
data1107 <- c(15.75,15.5,15.5,17,16.5,18.25,18.5,19.75,21.5,20.75,19,20.75,18.5,20.5,20,22.25,21.5,23.5,24.75,23.75)
data1107.aov <- aov(data1107~A+B)
summary(data1107.aov)

参考资料

1. 白话统计.冯国双
2. 医学统计学.孙振球.第四版
3. 数据分析：R语言实战.李诗羽.电子工业出版社