R语言学习系列27方差分析.docx

1、R语言学习系列27方差分析22. 方差分析一、方差分析原理1.方差分析概述 方差分析可用来研究多个分组的均值有无差异， 其中分组是按影 响因素的不同水平值组合进行划分的。方差分析是对总变异进行分析。看总变异是由哪些部分组成的， 这些部分间的关系如何。方差分析， 是用来检验两个或两个以上均值间差别显著性 （影响 观察结果的因素： 原因变量（列变量）的个数大于 2，或分组变量 （行 变量）的个数大于 1）。一元时常用 F 检验（也称一元方差分析） ，多 元时用多元方差分析（最常用 Wilks 检验）。方差分析可用于：（ 1）完全随机设计（单因素） 、随机区组设计（双因素） 、析因 设计、拉丁方设计

2、和正交设计等资料；（2）可对两因素间交互作用差异进行显著性检验；（ 3）进行方差齐性检验。要比较几组均值时， 理论上抽得的几个样本， 都假定来自正态总 体，且有一个相同的方差， 仅仅均值可以不相同。还需假定每一个观 察值都由若干部分累加而成， 也即总的效果可分成若干部分， 而每一 部分都有一个特定的含义， 称之谓效应的可加性。 所谓的方差是离均 差平方和除以自由度，在方差分析中常简称为均方（ Mean Square）2.基本思想基本思想是，将所有测量值上的总变异按照其变异的来源分解为 多个部份，然后进行比较， 评价由某种因素所引起的变异是否具有统 计学意义。根据效应的可加性， 将总的 离均差平

3、方和 分解成若干部分， 每一 部分都与某一种效应相对应， 总自由度也被分成相应的各个部分， 各 部分的离均差平方除以各自的自由度得出各部分的均方， 然后列出方 差分析表算出 F 检验值，作出统计推断。方差分析的关键是 总离均差平方和 的分解，分解越细致， 各部分 的含义就越明确，对各种效应的作用就越了解，统计推断就越准确。效应项与试验设计或统计分析的目的有关， 一般有：主效应（包 括各种因素），交互影响项（因素间的多级交互影响） ，协变量（来自 回归的变异项），等等。当分析和确定了各个效应项 S 后，根据原始观察资料可计算出各 个离均差平方和 SS，再根据相应的自由度 df ，由公式 MS=S

4、Sd/ f ，求出 均方 MS，最后由相应的均方，求出各个变异项的 F 值，F 值实际上是 两个均方之比值，通常情况下，分母的均方是误差项的均方。根据 F值的分子、分母均方的自由度 f1和 f 2，在确定显著性水平 为 情况下，由 F（f 1, f2） 临界值表查得单侧 F界限值。当 F，不拒绝原假设 H0，说明不拒绝这个效应项的效应为 0 的原 假设，也即这个效应项是可能对总变异没有实质影响的；若 FF 则 P值，拒绝原假设 H0，也即这个效应项是很可能对总变异有实质影 响的。3.方差分析的实验设计 为了确定方差分析表中各个有关效应项，需要在试验设计阶段就 作出安排，再根据设计要求进行试验，

5、得出原始观察值，按原来设计 方案算出方差分析表中的各项。在试验设计阶段通常需要考虑如下 4 个方面： （1）研究的因变量 即试验所要观察的主要指标，一次试验时可以有多个观察指标， 方差分析时也可以同时对多个因变量进行分析；（2）因素和水平试验的因素（ factor ）可以是品种、人员、方法、时间、地区等 等，因素所处的状态叫水平（ level ）。在每一个因素下面可以分成若 干水平。（3）因素间的交互影响 多因素的试验设计，有时需要分析因素间的交互影响 （interaction ），2 个因素间的交互影响称为一级交互影响 （AB）；3 个因素间的交互影响称为二级交互影响（ ABC）。当交互影响

6、项呈现统计不显著时，表明各个因素独立，当呈现统 计显著时，就需要列出这个交互影响项的效应，以助于作出正确的统 计推断。举例解释上述概念 ：要考察焦虑症的治疗疗效，一个因素是治疗 方案，有 2 种治疗方案，即该因素有 2 个水平；（治疗方案称为 组间因 子，因为每个患者只能被分配到一个组别中，没有患者同时接受两种 治疗）；再考虑一个因素治疗时间，也有两个水平：治疗 5 周和治疗 6 个月，同一患者在 5 周和 6 个月不止一次地被测量（两次） ，称为重复 测量（治疗时间称为 组内因子 ，因为每个患者在所有水平下都进行了 测量）。建立方差分析模型时， 既要考虑两个因素治疗方案和治疗时间 （主 效应

7、），又要考虑治疗方案和时间的交互影响（交互效应） ，此时即两 因素混合模型方差分析。当某个因素的各个水平下的因变量的均值呈现统计显著性差异时， 必要时可作两两水平间的比较，称为均值间的两两比较。二、 R语言实现 方差分析对数据的要求：满足正态性（来自同一正态总体）和方 差齐性（各组方差相等） ，在这两个条件下，若各组有差异，则只可 能是来自影响因素的不同水平。用 aov（） 函数进行方差分析，基本格式为：aov（formula, data=NULL, projections=FALSE, qr=TRUE, contrasts=NULL, .） 其中， formula 为方差分析公式；data

8、为数据框；projection 设置是否返回预测结果； qr 设置是否返回 QR分解结果； contrasts 为公式中一些因子的列表formula 公式的表示：(y 为因变量， ABC为分组因子)符号用法分隔符号，左边为响应变量，右边为解释变量eg： yA+B+C+分隔解释变量：表示变量的交互项 eg： yA+B+A:B*表示所有可能交互项eg： yA*B*C 可展开为： yA+B+C+A:B+A:C+B:C+A:B:C表示交互项达到次数eg： y(A+B+C)2展开为： yA+B+C+A:B+A:C+B:C.表示包含除因变量外的所有变量 eg：若一个数据框包括变量 y,A、B和 C，代码

9、y. 可展开为 yA+B+C常见研究设计的表达式 ：(小写字母表示定量变量，大写字母表示组别因子， Subject 是对被试者独有的标识变量)设计表达式单因素 ANOVAyA含单个协变量的单因素 ANCOVAyx+A双因素 ANOVAyA*B含两个协变量的双因素 ANCOVAyx1+x2+A*B随机化区组yB+A, B 为区组因子单因素组内 ANOVAyA+Error(Subject/A)含单个组内因子 (W)和单个组间因子(B) 的重复测量 ANOVAyB*W+Error(Subject/W)注意 ：非均衡设计时或存在协变量时， 效应项的顺序对结果影响 较大，越基础的效应越需要放在表达式前面

10、，首先是协变量、然后是 主效应、接着是双因素的交互项，再接着是三因素的交互项。若研究 不是正交的，一定要谨慎设置效应的顺序。有三种类型的方法可以分解 yA+B+A:B右边各效应对 y 所解释的 方差：类型 I （序贯型）效应根据表达式中先出现的效应做调整。 A 不做调整， B 根据 A 调整， A:B 交互项根据 A和 B调整。类型 II （分层型）效应根据同水平或低水平的效应做调整。 A根据 B 调整， B 依据 A调整， A:B 交互项同时根据 A和 B调整。类型 III （边界型）每个效应根据模型其他各效应做相应调整。 A根据 B 和 A:B 做调 整，A:B 交互项根据 A和 B调整。

11、R默认调用类型 I 方法，其他软件（比如 SAS和 SPSS）默认调用 类型 III 方法。 car 包中的 Anova（） 函数（不要与标准 anova（） 函数 混淆）提供了使用类型 II 或类型 III 方法的选项，而 aov（） 函数使 用的是类型 I 方法。若想使结果与其他软件（如 SAS和 SPSS）提供 的结果保持一致，可以使用 Anova（） 函数。三、单因素方差分析1 个因变量， 1 个影响因素：总差异 Yij = 平均差异 + 因素差异 i + 随机差异 ij 例 1 比较 4 种品牌的胶合板的耐磨性，各抽取 5 个样品，相同转速 磨损相同时间测得磨损深度( mm)，比较

12、4 个品牌胶合板的耐磨性有 无差异？部分数据如下( ex27_ex1.Rdata )：setwd(E:/ 办公资料 /R 语言 /R 语言学习系列 /codes) load(ex27_ex1.Rdata) head(datas)wear brand12.30 A22.32 A32.40 A42.45 A52.58 A62.35 B attach(datas) table(brand) #各组的样本数brandA B C D5 5 5 5aggregate(wear,by=list(brand),mean) #各组均值Group.1 x1A 2.4102B 2.4043C 2.0464#各组标准

13、差D 2.572 aggregate(wear,by=list(brand),sd)Group.1 x1A 0.112694282B 0.117601023C 0.112160604D 0.03271085library(car)#用 Q-Q图检验数据qqPlot(lm(wearbrand,data=datas),simulate=TRUE) 的正态性leveneTest(wearas.factor(brand),data=datas) #方差齐性检验Levenes Test for Homogeneity of Variance (center = median)Df F value Pr(

14、F) group 3 0.6987 0.566416fitF)brand 3 0.7398 0.24660 24.55 3.15e-06 *Residuals 16 0.1607 0.01005Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1说明 ：方差齐性检验， 原假设 H0：方差齐，p 值=0.56640.05, 故 接受原假设，即方差齐。单因素方差分析结果， brand是因素， Residuals 是残差，各列依 次为自由度、平方和、均方和、 F统计量， p值=3.15e-060.05, 拒绝 原假设，即不同品牌的磨损(均值)有显著差别。lib

15、rary(gplots)plotmeans（wearbrand,xlab= 品牌, ylab= 磨损） #图形展示带 95%置信 区间的各组均值通过前面的分析知道，不同品牌的磨损（均值）有显著差别，但并不知道哪个品牌与其它品牌有显著差别。 TukeyHSD（）函数提供了对各组均值差异的成对检验TukeyHSD(fit)Tukey multiple comparisons of means95% family-wise confidence levelFit: aov(formula = wear brand, data = datas) $branddiff lwr upr p adjB-A

16、-0.006 -0.18735345 0.1753535 0.9996826C-A -0.364 -0.54535345 -0.1826465 0.0001610D-A 0.162 -0.01935345 0.3433535 0.0886142C-B -0.358 -0.53935345 -0.1766465 0.0001929D-B 0.168 -0.01335345 0.3493535 0.0744337D-C 0.526 0.34464655 0.7073535 0.0000019说明 ：可以看出（H0：无差异），B与A的差异非常不显著， C与A、C与B、D与C的差异非常显著multco

17、mp包中的 glht() 函数提供了更为全面的多重均值比较方法。 library(multcomp) attach(datas)tuk |t|)B -A = 0-0.006000.06339-0.0950.9997C- A = 0-0.364000.06339-5.7420.001 *D- A = 00.162000.063392.5560.0886 .C- B = 0-0.358000.06339-5.6480.001 *D- B = 00.168000.063392.6500.0743 .D- C = 00.526000.063398.2980.001 *Signif. codes: 0

18、* 0.001 * 0.01 * 0.05(Adjusted p values reported - single-step method)plot(cld(tuk, level = 0.05), col = lightgrey)说明：标记相同字母（标记 b）的品牌 ABD认为是无显著差异，在 同一亚组，而品牌 C（标记 a）与另外三个品牌有显著差异。另外，也可以进行多重 t 检验，使用函数：pairwise.t.test（x, g, p.adjust.method=,.）其中， x 为因变量， g 为因子型的分组变量；p.adjust.method 设置 p 值的修正方法，由于多次重复 t

19、检验会 大大增加犯第一类错误的概率，为此要进行 p 值的修正，使用 bonferroni 法修正效果较好。pairwise.t.test(wear,brand,p.adjust.method=bonferroni)Pairwise comparisons using t tests with pooled SD data: wear and brandA B CB 1.00000 - -C 0.00018 0.00022 -D 0.12695 0.10474 2.1e-06P value adjustment method: bonferroni说明：原假设 H0: 无差异，可见 A与 B无差

20、异， C与 ABD有显著 差异。最后，方差分析对离群点非常敏感，检验是否有离群点：library(car)outlierTest(fit)No Studentized residuals with Bonferonni p 0.05Largest |rstudent|:rstudent unadjusted p-value Bonferonni p9 2.528103 0.023182 0.46364说明：经检验无离群点。三、两因素方差分析1 个因变量， 2 个影响因素：总差异 Yijk = 平均差异 + 因素 1 差异 i + 因素 2 差异 i+ 因素 1,2 交互作用差异 ij + 随机

21、差异 ijk例 2 研究 60 只豚鼠的牙齿生长数据，按 2 种喂食方法：橙汁、维生 素 C，各喂食方法中抗坏血酸含量都有 3 个水平： 0.5mg、1mg、2mg，分配为 6 组，每组各 10 只，牙齿长度为因变量。 做两因素方差分析。attach(ToothGrowth) head(ToothGrowth)len supp dose1 4.2 VC 0.52 11.5 VC 0.53 7.3 VC 0.54 5.8 VC 0.55 6.4 VC 0.56 10.0 VC 0.5table(supp, dose) #各组样本数相同，即为均衡设计dosesupp 0.5 1 2OJ 10 10

22、 10VC 10 10 10 aggregate(len, by=list(supp, dose), mean)Group.1 Group.2 x#计算各组均值1OJ0.5 13.232VC0.5 7.983OJ1.0 22.704VC1.0 16.775OJ2.0 26.066VC2.0 26.14#计算各组标准差1OJ0.5 4.4597092VC0.5 2.7466343OJ1.0 3.9109534VC1.0 2.5153095OJ2.0 2.6550586VC2.0 4.797731aggregate(len, by=list(supp, dose), sd)Group.1 Grou

23、p.2 xbartlett.test(lensupp,data=ToothGrowth) 检验#关于因素 supp 的方差齐性Bartlett test of homogeneity of variancesdata: len by suppBartletts K-squared = 1.4217, df = 1, p-value = 0.2331 bartlett.test(lendose,data=ToothGrowth) #关于因素 dose 的方差齐性检验Bartlett test of homogeneity of variancesdata: len by doseBartlett

24、s K-squared = 0.66547, df = 2, p-value = 0.717 fitF)supp 1 205.4 205.4 12.317 0.000894 *dose 1 2224.3 2224.3 133.415 2e-16 *Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1说明：可以看出，主效应 supp和 dose都非常显著( p 值都远小于 0.05 ），交互效应也显著（p 值=0.02460.05 ）。若交互作用不显著，可以可以只做去掉交互效应的方差分析。图形化展示两因素方差分析的交互效应：par(mfrow=c(1,2)i

25、nteraction.plot(dose, supp, len, type=b, col = c(red, blue), pch = c(16, 18), main=Interaction between Dose and Supp)interaction.plot(supp, dose, len, type=b, col=c(red, blue), pch = c(16, 18), main=Interaction between Supp and Dose)说明：有一个图的线有交叉，说明有交互作用。可以看出随着橙 汁和维生素 C中的抗坏血酸剂量的增加，牙齿长度变长；。对于 0.5 mg 和

26、 1 mg剂量，橙汁比维生素 C更能促进牙齿生长；对于 2 mg剂量的抗坏血酸，两种喂食方法下牙齿长度增长相同也可以用 HH包中的 interaction2wt（） 函数（也适合三因素方差分析）来展示更全面的可视化结果：library(HH)interaction2wt(lensupp*dose)三、重复测量方差分析重复测量方差分析，即受试者被测量不止一次例 3（1 个组内 1 个组间因子的重复测量） 在某浓度 CO2 的环境中，对寒带植物（来自魁北克）和非寒带植物的（来自密西西比）光合作用率进行比较。因变量 uptake 为 CO2 吸收量，自变量 Type（组间因子）为植物类型，自变量 c

27、onc（组内因子）为七种水平的 CO2 浓度。attach(CO2)head(CO2) #注意 CO2是长格式的数据Plant Type Treatment conc uptake1Qn1 Quebec nonchilled 95 16.02Qn1 Quebec nonchilled 175 30.43Qn1 Quebec nonchilled 250 34.84Qn1 Quebec nonchilled 350 37.25Qn1 Quebec nonchilled 500 35.36Qn1 Quebec nonchilled 675 39.2 w1b1-subset(CO2, Treatment=chilled) #先只考虑寒带植物fitF)Type 1 2667.2 2667.2 60.41 0.00148 *Residuals 4 176.6 44.1Error: Plant:concDf Sum Sq Mean Sq F value Pr(F)conc 1 888.6 888.6 215.46 0.000125 *conc:Type 1Residuals 4239.216.5239.24.158.01 0.001595 *1Signif.codes:00.001 * 0.01* 0.05 . 0.