文章详细介绍了R语言在执行方差分析（ANOVA）中的应用，包括单因素、双因素以及重复测量方差分析。通过具体的数据集和实例，展示了如何分析不同因素对结果的影响，如治疗方案（CBT和EMDR）对焦虑度的影响，以及药物疗法对胆固醇降低的效果。此外，还涉及了正态性、方差齐性和多重比较的检验，以及协方差分析和交互效应的探讨。
一 基本术语
例子：在治疗焦虑症的方案中，有两种，一种为认知行为疗法（CBT），另一种为眼动脱敏再加工法（EMDR），在治疗结束时要求每位患者填写焦虑度测量的自我评测报告（STAI）
因此在这个试验中，
组间因子：治疗方案，是两水平的（CBT、EMDR），也是自变量
因变量：STAI
均衡设计：自变量不同水平的观测数相同，反之则是非均衡设计
如果只对CBT的治疗效果感兴趣，则将患者都放在CBT中，然后在不同的治疗时间中评价疗效，则设计方案为如下所示：
此时只考虑CBT时，时间就是它的组内因子，因为每位患者在所有水平下都进行了测量，因此这种统计设计称为单因素方差分析，而每个受试者都不止一次被测量，也称作重复测量方差分析，如果F值显著，说明因变量得分均值在不同时间内发生了改变。
若既考虑治疗方案，又考虑时间的差异
此时疗法和时间都作为因子时，可以同时分析单独的主效应影响，还可分析两因子的交互影响。
二 ANOVA模型拟合
R语言中，可使用单独的aov（）函数拟合，也可使用回归分析中的lm（）函数分析，在这里主要讨论aov（）函数，语法为aov（formula，data=dataframe），如下为formula中的一些符号说明
在分析时，因子的顺序很重要，R语言默认序贯型方法计算ANOVA效应，如
三 单因素方差分析案例
单因素方差分析中，你感兴趣的是比较分类因子定义的两个或多个组别中的因变量均值。以multcomp包中的cholesterol数据集为例，50个患者均接受降低胆固醇药物治疗（trt）五种疗法中的一种疗法。其中三种治疗条件使用药物相同，分别是20 mg一天一次（1time）、10 mg一天两次（2times）和5 mg一天四次（4times）。剩下的两种方式（drugD和drugE）代表候选药物。进而探究哪种药物疗法降低胆固醇（响应变量）最多
因此这里的自变量只有一个，为治疗（trt），因变量则为response
#方差分析步骤
library(multcomp)
attach(cholesterol)
table(trt)
#trt
1time 2times 4times drugD drugE
10 10 10 10 10
aggregate(response, by=list(trt), FUN=mean) #检验各组均值
# Group.1 x
1 1time 5.78197
2 2times 9.22497
3 4times 12.37478
4 drugD 15.36117
5 drugE 20.94752
aggregate(response, by=list(trt), FUN=sd) #检验各组标准差
# Group.1 x
1 1time 2.878113
2 2times 3.483054
3 4times 2.923119
4 drugD 3.454636
5 drugE 3.345003
fit <- aov(response ~ trt) #检验组间差异
summary(fit)
# Df Sum Sq Mean Sq F value Pr(>F)
trt 4 1351.4 337.8 32.43 9.82e-13 ***
Residuals 45 468.8 10.4
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
library(gplots) #绘制各组均值和置信区间的图形
plotmeans(response ~ trt, xlab="Treatment", ylab="Response",
main="Mean Plot\nwith 95% CI")
detach(cholesterol)
从输出结果可以看到每10个患者接受其中一个药物疗法，均值显示drugE降低胆固醇最多，而1time降低的最少，各组的标准差相对恒定，在2.88到3.48间浮动，F值为显著，说明五种疗法的效果确实有不同的效果。
在以上的方差分析中，只是得出了五种疗法的效果不同，但没有给出哪种疗法与其他疗法不同，因此这里的多重比较可以解决这个问题
#多重比较
TukeyHSD(fit)
# Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = response ~ trt)
$trt
diff lwr upr p adj
2times-1time 3.44300 -0.6582817 7.544282 0.1380949
4times-1time 6.59281 2.4915283 10.694092 0.0003542
drugD-1time 9.57920 5.4779183 13.680482 0.0000003
drugE-1time 15.16555 11.0642683 19.266832 0.0000000
4times-2times 3.14981 -0.9514717 7.251092 0.2050382
drugD-2times 6.13620 2.0349183 10.237482 0.0009611
drugE-2times 11.72255 7.6212683 15.823832 0.0000000
drugD-4times 2.98639 -1.1148917 7.087672 0.2512446
drugE-4times 8.57274 4.4714583 12.674022 0.0000037
drugE-drugD 5.58635 1.4850683 9.687632 0.0030633
par(las=2) # 旋转轴标签
par(mar=c(5,8,4,2)) #调整图形的单位英分
plot(TukeyHSD(fit))
par(opar)
上图为Tukey HSD成对比较图，图形中的置信区间包括0的疗法说明差异不显著（p>0.5）,以下使用了glht（）函数，它提供了更为全面的比较方法
library(multcomp)
par(mar=c(5,4,6,2))
tuk <- glht(fit, linfct=mcp(trt="Tukey"))
plot(cld(tuk, level=.05),col="lightgrey")
par(opar)
在该图的上部分，有相同的字母说明均值差异不显著，例如1times和2times差异不显著，因为有相同字母a，该图也得出drugE比其他药物和疗法更优的结论。
对于结果来说，在做统计检验时，数据满足假设条件的程度可以使我们更信赖得到的结果。
#正态性假设验证
library(car)
qqPlot(lm(response ~ trt, data=cholesterol),
simulate=TRUE, main="Q-Q Plot", labels=FALSE)
qqplot（）要求用lm（）拟合，数据落在95%置信区间范围内，说明满足正态性假设。
# 方差齐性验证
bartlett.test(response ~ trt, data=cholesterol)
#Bartlett test of homogeneity of variances
data: response by trt
Bartlett's K-squared = 0.57975, df = 4, p-value = 0.9653
# 检验离群点
library(car)
outlierTest(fit)
# No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
rstudent unadjusted p-value Bonferroni p
19 2.251149 0.029422 NA在bartlett检验中表明五组的方差并没有显著不同（p=0.97），由于方差齐性分析对离群点非常敏感，因此可利用car包中的outlierTest()函数检测离群点，从输出结果中也可判断出没有证据说明数据中含有离群点。因此根据三个检验可说明ANOVA模型拟合效果较好的结论。
四 单因素协方差分析
单因素协方差分析包含一个或多个定量的协变量，下面以multcomp包中的litter数据集为例，设怀孕小鼠分为四个小组，每个小组接受不同剂量（0、5、50或500）的药物处理，产下幼崽的体重均值设为因变量，怀孕时间为协变量，代码如下所示
data(litter, package="multcomp")
attach(litter)
table(dose)
#dose
0 5 50 500
20 19 18 17
aggregate(weight, by=list(dose), FUN=mean) #查看不同剂量下的体重均值
# Group.1 x
1 0 32.30850
2 5 29.30842
3 50 29.86611
4 500 29.64647
fit <- aov(weight ~ gesttime + dose)
summary(fit)
# Df Sum Sq Mean Sq F value Pr(>F)
gesttime 1 134.3 134.30 8.049 0.00597 **
dose 3 137.1 45.71 2.739 0.04988 *
Residuals 69 1151.3 16.69
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1根据ANOVA的F检验，可知：（1）怀孕时间与幼崽的出生体重有关；（2）控制怀孕时间，发现药物剂量的不同确定导致出生体重的均值不同。
由于使用了协变量，若需要获取去除协变量效应后的组均值，可以使用effects包中的函数计算
library(effects)
effect("dose", fit)
# dose effect
dose
0 5 50 500
32.35367 28.87672 30.56614 29.33460 接下来利用多重比较，得出哪种剂量对体重均值的影响最大
#多重比较
library(multcomp)
contrast <- rbind("no drug vs. drug" = c(3, -1, -1, -1)) #设定第一组和其他三组作对比
summary(glht(fit, linfct=mcp(dose=contrast)))
# Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: User-defined Contrasts
Fit: aov(formula = weight ~ gesttime + dose)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
no drug vs. drug == 0 8.284 3.209 2.581 0.012 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)根据检验统计量p值<0.05，可得出未用药组比其他组的出生体重更高的结论（同样，也可利用其他组作为对照，比较和另外三组的统计量是否显著，从而得出结论）
最后，单因素协方差分析（AMCOVA）除了正态性和同方差假定，还需假定回归斜率相同，因此在本例中需要假定四个组通过怀孕时间来预测出生体重的回归斜率都相同，ANCOVA模型包含怀孕时间*剂量的交互项时，可对回归斜率的同质性进行检验。交互效应若显著，则意味着时间和幼崽出生体重间的关系依赖于药物剂量的水平。
#检验回归斜率同质性
library(multcomp)
fit2 <- aov(weight ~ gesttime*dose, data=litter)
summary(fit2)
#Analysis of Variance Table
Response: weight
Df Sum Sq Mean Sq F value Pr(>F)
gesttime 1 134.30 134.304 8.0493 0.005971 **
dose 3 137.12 45.708 2.7394 0.049883 *
Residuals 69 1151.27 16.685
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1可以看出，交互效应不显著，支持斜率相同的假设（若假设不成立则可尝试变换协变量或因变量），接下来是对结果的可视化
library(HH)
ancova(weight ~ gesttime + dose, data=litter)
该函数可以绘制因变量、协变量和因子之间的关系图。从图中可以看到，用怀孕时间来预测出生体重的回归线相互平行，只是截距项不同。随着怀孕时间增加，幼崽出生体重也会增加。另外，还可以看到0剂量组截距项最大，5剂量组截距项最小。由于上面的设置，直线会保持平行。
五 双因素方差分析
在该分析中，受试者被分配到两因子的交叉类别组中，以ToothGrowth数据集为例，随机分配60只豚鼠，分别采用两种喂食方法（橙汁或维生素C），各喂食方法中抗坏血酸含量有三种水平（0.5 mg、1 mg或2 mg），每种处理方式组合都被分配10只豚鼠。使牙齿长度为因变量，分析结果如下
attach(ToothGrowth)
table(supp,dose) #输出的结果表明该设计是均衡设计
# dose
supp 0.5 1 2
OJ 10 10 10
VC 10 10 10
aggregate(len, by=list(supp,dose), FUN=mean) #查看不同喂食方法和喂食水平下的牙齿长度均值
# Group.1 Group.2 x
1 OJ 0.5 13.23
2 VC 0.5 7.98
3 OJ 1.0 22.70
4 VC 1.0 16.77
5 OJ 2.0 26.06
6 VC 2.0 26.14
aggregate(len, by=list(supp,dose), FUN=sd) #查看不同喂食方法和喂食水平下的牙齿长度标准差
# Group.1 Group.2 x
1 OJ 0.5 4.459709
2 VC 0.5 2.746634
3 OJ 1.0 3.910953
4 VC 1.0 2.515309
5 OJ 2.0 2.655058
6 VC 2.0 4.797731
dose <- factor(dose) #dose转化为因子变量
fit <- aov(len ~ supp*dose)
summary(fit)
# Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.4 205.4 15.572 0.000231 ***
dose 2 2426.4 1213.2 92.000 < 2e-16 ***
supp:dose 2 108.3 54.2 4.107 0.021860 *
Residuals 54 712.1 13.2
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1从方差分析的结果中可以看出主效应和交互效应都是显著的，因此我们也可展示出交互效应的关系，如下所示（三种方法）：
#第一种
interaction.plot(dose, supp, len, type="b",
col=c("red","blue"), pch=c(16, 18),
main = "Interaction between Dose and Supplement Type")
#第二种
library(gplots)
plotmeans(len ~ interaction(supp, dose, sep=" "),
connect=list(c(1, 3, 5),c(2, 4, 6)),
col=c("red","darkgreen"),
main = "Interaction Plot with 95% CIs",
xlab="Treatment and Dose Combination")
library(HH)
interaction2wt(len~supp*dose)
以上三幅图形都表明随着橙汁和维生素C中的抗坏血酸剂量的增加，牙齿长度变长。对于0.5 mg
和1 mg剂量，橙汁比维生素C更能促进牙齿生长；对于2 mg剂量的抗坏血酸，两种喂食方法下牙
齿长度增长相同。在三种绘图方法中HH包中的interaction2wt()函数能展示任意复杂度设计（双因素方差分析、三因素方差分析等）的主效应（箱线图）和交互效应。
六 重复测量方差分析
所谓重复测量方差分析，即受试者被测量不止一次。这里以包含一个组内和一个组间因子的重复测量方差分析为例。示例来源于生理生态学领域，研究方向是生命系统的生理和生化过程如何响应环境因素的变异（此为应对全球变暖的一个非常重要的研究领域）。基础安装包中的CO2数据集包含了北方和南方牧草类植物Echinochloa crus-galli 的寒冷容忍度研究结果，在某浓度二氧化碳的环境中，对寒带植物与非寒带植物的光合作用率进行了比较。
首先，我们关注寒带植物。因变量是二氧化碳吸收量（uptake），单位为ml/L，自变量是植
物类型Type（魁北克VS密西西比州）和七种水平（95~1000 umol/m^2 sec）和二氧化碳浓度
（conc）。另外，Type是组间因子，conc是组内因子。
CO2$conc <- factor(CO2$conc) #转化为因子
w1b1 <- subset(CO2, Treatment=='chilled')
fit <- aov(uptake ~ (conc*Type) + Error(Plant/(conc)), w1b1)
summary(fit)
#Error: Plant
Df Sum Sq Mean Sq F value Pr(>F)
Type 1 2667.2 2667.2 60.41 0.00148 **
Residuals 4 176.6 44.1
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: Plant:conc
Df Sum Sq Mean Sq F value Pr(>F)
conc 6 1472.4 245.40 52.52 1.26e-12 ***
conc:Type 6 428.8 71.47 15.30 3.75e-07 ***
Residuals 24 112.1 4.67
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
par(las=2)
par(mar=c(10,4,4,2))
with(w1b1,
interaction.plot(conc,Type,uptake,
type="b", col=c("red","blue"), pch=c(16,18),
main="Interaction Plot for Plant Type and Concentration"))
boxplot(uptake ~ Type*conc, data=w1b1, col=(c("gold","green")),
main="Chilled Quebec and Mississippi Plants",
ylab="Carbon dioxide uptake rate (umol/m^2 sec)")方差分析表表明在0.01的水平下，主效应类型和浓度以及交叉效应类型x浓度都非常显著，下面展示交互效应
从以上任意一幅图都可以看出，魁北克省的植物比密西西比州的植物二氧化碳吸收率高，而
且随着CO2浓度的升高，差异越来越明显。
七 多元方差分析
当因变量不止一个时，可选择多元方差分析（MANOVA）对他们同时分析，使用MASS包里的UScereal数据集为例，研究美国谷物中的卡路里、脂肪和糖含量是否会因为储存架位置的不同而发生变化。其中1代表底层货架，2代表中层货架，3代表顶层货架。卡路里、脂肪和糖含量是因变量，货架是三水平（1、2、3）的自变量。分析过程如下
library(MASS)
attach(UScereal)
shelf <- factor(shelf) #转为因子变量
y <- cbind(calories, fat, sugars) #将三个因变量合成为矩阵
aggregate(y, by=list(shelf), FUN=mean) #查看不同货架下的因变量均值
#Group.1 calories fat sugars
1 1 119.4774 0.6621338 6.295493
2 2 129.8162 1.3413488 12.507670
3 3 180.1466 1.9449071 10.856821
cov(y) #矩阵的协方差值
# calories fat sugars
calories 3895.24210 60.674383 180.380317
fat 60.67438 2.713399 3.995474
sugars 180.38032 3.995474 34.050018
fit <- manova(y ~ shelf)
summary(fit)
# Df Pillai approx F num Df den Df Pr(>F)
shelf 2 0.4021 5.1167 6 122 0.0001015 ***
Residuals 62
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
summary.aov(fit） #对每一个因变量都做单因素方差分析
Response calories :
Df Sum Sq Mean Sq F value Pr(>F)
shelf 2 50435 25217.6 7.8623 0.0009054 ***
Residuals 62 198860 3207.4
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Response fat :
Df Sum Sq Mean Sq F value Pr(>F)
shelf 2 18.44 9.2199 3.6828 0.03081 *
Residuals 62 155.22 2.5035
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Response sugars :
Df Sum Sq Mean Sq F value Pr(>F)
shelf 2 381.33 190.667 6.5752 0.002572 **
Residuals 62 1797.87 28.998
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1由于上面的多元检验是显著的，因此可对每一个因变量进行单因素方差分析，从上述结果可以看到，三组中每种营养成分的测量值都是不同的。另外，还可以用均值比较步骤（比如TukeyHSD）来判断对于每个因变量，哪种货架与其他货架都是不同的。
下面进行统计检验
center <- colMeans(y)
n <- nrow(y)
p <- ncol(y)
cov <- cov(y)
d <- mahalanobis(y,center,cov)
coord <- qqplot(qchisq(ppoints(n),df=p),
d, main="QQ Plot Assessing Multivariate Normality",
ylab="Mahalanobis D2")
abline(a=0,b=1)
identify(coord$x, coord$y, labels=row.names(UScereal))
若数据服从多元正态分布，那么点将落在直线上。你能通过identify()函数交互性地对图中的点进行鉴别。从图形上看，有两个观测点是异常状态，数据集似乎违反了多元正态性。可以删除这两个点再重新分析。
扫一扫
4万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
