R语言实战-topic7方差分析

1、Topic7 方差分析一、相关术语以焦虑症治疗为例，现有两种治疗方案：认知行为疗法（CBT）和眼动脱敏再加工法（EMDR）。我们招募10位焦虑症患者作为志愿者，随机分配一半的人接受为期五周的CBT，另外一半接受为期五周的EMDR，设计方案如表9-1所示。在治疗结束时，要求每位患者都填写状态特质焦虑问卷（STAI），也就是一份焦虑度测量的自我评测报告。在这个实验设计中，治疗方案是两水平（CBT、 EMDR）的组间因子。之所以称其为组间因子，是因为每位患者都仅被分配到一个组别中，没有患者同时接受CBT和EMDR。表中字母s代表受试者（患者）。 STAI是因变量，治疗方案是自变量。由于在每种治疗方案

2、下观测数相等，因此这种设计也称为均衡设计（balanced design）；若观测数不同，则称作非均衡设计（unbalanceddesign）。因为仅有一个类别型变量，表9-1的统计设计又称为单因素方差分析（one-way ANOVA），或进一步称为单因素组间方差分析。方差分析主要通过F检验来进行效果评测，若治疗方案的F检验显著，则说明五周后两种疗法的STAI得分均值不同。假设你只对CBT的效果感兴趣，则需将10个患者都放在CBT组中，然后在治疗五周和六个月后分别评价疗效，设计方案如表9-2所示。疗法（therapy）和时间（time）都作为因子时，我们既可分析疗法的影响（时间跨度上的平均）和

3、时间的影响（疗法类型跨度上的平均），又可分析疗法和时间的交互影响。前两个称作主效应，交互部分称作交互效应。当设计包含两个甚至更多的因子时，便是多因素方差分析设计，比如两因子时称作双因素方差分析，三因子时称作三因素方差分析，以此类推。若因子设计包括组内和组间因子，又称作混合模型方差分析，当前的例子就是典型的双因素混合模型方差分析。本例中，你将做三次F检验：疗法因素一次，时间因素一次，两者交互因素一次。若疗法结果显著，说明CBT和EMDR对焦虑症的治疗效果不同；若时间结果显著，说明焦虑度从五周到六个月发生了变化；若两者交互效应显著，说明两种疗法随着时间变化对焦虑症治疗影响不同（也就是说，焦虑度从五

4、周到六个月的改变程度在两种疗法间是不同的）。现在，我们对上面的实验设计稍微做些扩展。众所周知，抑郁症对病症治疗有影响，而且抑郁症和焦虑症常常同时出现。即使受试者被随机分配到不同的治疗方案中，在研究开始时，两组疗法中的患者抑郁水平就可能不同，任何治疗后的差异都有可能是最初的抑郁水平不同导致的，而不是由于实验的操作问题。抑郁症也可以解释因变量的组间差异，因此它常称为混淆因素（confounding factor）。由于你对抑郁症不感兴趣，它也被称作干扰变数（nuisance variable）。假设招募患者时使用抑郁症的自我评测报告，比如白氏抑郁症量表（BDI），记录了他们的抑郁水平，那么你可以在

5、评测疗法类型的影响前，对任何抑郁水平的组间差异进行统计性调整。本案例中， BDI为协变量，该设计为协方差分析（ANCOVA）。以上设计只记录了单个因变量情况（STAI），为增强研究的有效性，可以对焦虑症进行其他的测量（比如家庭评分、医师评分，以及焦虑症对日常行为的影响评价）。当因变量不止一个时，设 计 被 称 作 多 元 方 差 分 析 （MANOVA ）， 若 协 变 量 也 存 在 ， 那 么 就 叫 多 元 协 方 差 分析（MANCOVA）。二、单因素方差分析（1）单因素方差分析以multcomp包中的cholesterol数据集为例（取自Westfall、 Tobia、 Rom、 H

6、ochberg， 1999）， 50个患者均接受降低胆固醇药物治疗（trt）五种疗法中的一种疗法。其中三种治疗条件使用药物相同，分别是20mg一天一次（1time）、 10mg一天两次（2times）和5mg一天四次（4times）。剩下的两种方式（drugD和drugE）代表候选药物。哪种药物疗法降低胆固醇（响应变量）最多呢？#单因素方差分析（感兴趣地是比较分类因子定义的两个或多个组别中的因变量均值） install.packages("multcomp") library(multcomp) attach(cholesterol) str(cholesterol) ch

7、olesterol table(trt) aggregate(response,by=list(trt),FUN=mean) aggregate(response,by=list(trt),FUN=sd) fit<-aov(responsetrt) summary(fit) library(gplots) plotmeans(responsetrt,xlab="Treatment",ylab="Response",main="Mean Plotn with 95%CI")（2）多重比较虽然ANOVA对各疗法的F检验表明五种药物疗

8、法效果不同，但是并没有告诉你哪种疗法与其他疗法不同。多重比较可以解决这个问题。TukeyHSD(fit) par(las=2) par(mar=c(5,4,6,2) plot(TukeyHSD(fit) （3）评估检验的假设条件当因变量服从正态分布，各组方差相等时，可用Q-Q图来检验正态性假设qqPlot（）要求用lm（）拟合，若数据落在95%的置信区间范围内，说明满足正态性假设。library(car)qqPlot(lm(response trt, data=cholesterol), simulate=TRUE, main="Q-Q Plot", labels=FALS

9、E)R提供了一些可用来做方差齐性检验的函数bartlett.test(response trt, data=cholesterol)离群点检验library(car)outlierTest(fit)三、单因素协方差分析单因素协方差分析（ANCOVA）扩展了单因素方差分析（ANOVA），包含一个或多个定量的协变量。下面的例子来自于multcomp包中的litter数据集（见Westfall et al.， 1999）。怀孕小鼠被分为四个小组，每个小组接受不同剂量（0、 5、 50或500）的药物处理。产下幼崽的体重均值为因变量，怀孕时间为协变量。data(litter, package=&quo

10、t;multcomp")attach(litter)table(dose) aggregate(weight, by=list(dose), FUN=mean)fit <- aov(weight gesttime + dose) summary(fit)（1）评估检验的假设条件ANCOVA还假定回归斜率相同。本例中，假定四个处理组通过怀孕时间来预测出生体重的回归斜率都相同。library(multcomp) fit2<-aov(weightgesttime*dose,data=litter) summary(fit2)（2）结果可视化HH包中的ancova()函数可以绘制

11、因变量、协变量和因子之间的关系图。library(HH)ancova(weight gesttime + dose, data=litter)四、双因素方差分析在 双 因 素 方 差 分 析 中 ， 受 试 者 被 分 配 到 两 因 子 的 交 叉 类 别 组 中 。 以 基 础 安 装 中 的ToothGrowth数据集为例，随机分配60只豚鼠，分别采用两种喂食方法（橙汁或维生素C），各喂食方法中抗坏血酸含量有三种水平（0.5mg、 1mg或2mg），每种处理方式组合都被分配10只豚鼠。牙齿长度为因变量attach(ToothGrowth)table(supp,dose)aggregate

12、(len, by=list(supp,dose), FUN=mean)aggregate(len, by=list(supp,dose), FUN=sd)dose <- factor(dose)fit <- aov(len supp*dose)summary(fit)五、重复测量方差分析所谓重复测量方差分析，即受试者被测量不止一次。本节重点关注含一个组内和一个组间因子的重复测量方差分析（这是一个常见的设计）。示例来源于生理生态学领域，研究方向是生命系统的生理和生化过程如何响应环境因素的变异（此为应对全球变暖的一个非常重要的研究领域）。基础安装包中的CO2数据集包含了北方和南方牧草类

13、植物Echinochloa crus-galli （Potvin、Lechowicz、 Tardif， 1990）的寒冷容忍度研究结果，在某浓度二氧化碳的环境中，对寒带植物与非寒带植物的光合作用率进行了比较。研究所用植物一半来自于加拿大的魁北克省（Quebec），另一半来自美国的密西西比州（Mississippi）。首先，我们关注寒带植物。因变量是二氧化碳吸收量（uptake），单位为ml/L，自变量是植物类型Type（魁北克VS.密西西比）和七种水平（951000 umol/m2 sec）的二氧化碳浓度（conc）。另外， Type是组间因子， conc是组内因子。 Type已经被存储为一个

14、因子变量，但你还需要先将conc转换为因子变量含一个组间因子和一个组内因子的重复测量方差分析CO2$conc <- factor(CO2$conc)w1b1 <- subset(CO2, Treatment='chilled')fit <- aov(uptake (conc*Type) + Error(Plant/(conc), w1b1)summary(fit)par(las=2)par(mar=c(10,4,4,2)with(w1b1, interaction.plot(conc,Type,uptake, type="b", col=c

15、("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/m2

16、 sec)")par(opar)六、多元方差分析当因变量（结果变量）不止一个时，可用多元方差分析（MANOVA）对它们同时进行分析。以MASS包中的UScereal数据集为例（Venables， Ripley（1999），我们将研究美国谷物中的卡路里、脂肪和糖含量是否会因为储存架位置的不同而发生变化；其中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)cov(y)fit <- manova(y shelf)summary(fit)summary.aov(fit)（1）评估假设检验单因素多元方差分析有两个前提假设，一个是多元正态性，一个是方差-协方差矩阵同质性。第一个假设即指因变量组合成的向量服从一个多元正态分布。可以用Q-Q图来检验该假设条件center <- colMeans(y)n <- nrow(y)p <