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

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

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

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

4、下，分母的均方是误差项的均方。根据 F 值的分子、分母均方的自由度 f1 和 f2，在确定显著性水平为 情况下，由 F(f1, f2)临界值表查得单侧 F 界限值。当 F，不拒绝原假设 H0，说明不拒绝这个效应项的效应为 0 的原假设，也即这个效应项是可能对总变异没有实质影响的；若 FF 则 P值 ，拒绝原假设 H0，也即这个效应项是很可能对总变异有实质影响的。3.方差分析的实验设计为了确定方差分析表中各个有关效应项，需要在试验设计阶段就作出安排，再根据设计要求进行试验，得出原始观察值，按原来设计方案算出方差分析表中的各项。在试验设计阶段通常需要考虑如下 4 个方面：（1）研究的因变量即试验所

5、要观察的主要指标，一次试验时可以有多个观察指标，方差分析时也可以同时对多个因变量进行分析；（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 为数据框；projection 设置是否返回预测结果；qr 设置是否返回 QR 分解结果；contrasts 为公式中一些因子的列表。formula 公式的表示： （ y 为因变量，AB

8、C 为分组因子）符号 用法 分隔符号，左边为响应变量，右边为解释变量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，代码 y.可展开为yA+B+C常见研究设计的表达式：（小写字母表示定量变量，大写字母表示组别因子，Subject 是对被试者独有的标识变量）设计 表达式单因素 ANOVA yA含单个协变

9、量的单因素 ANCOVA yx+A双因素 ANOVA yA*B含两个协变量的双因素 ANCOVA yx1+x2+A*B随机化区组 yB+A, B 为区组因子单因素组内 ANOVA yA+Error(Subject/A)含单个组内因子(W)和单个组间因子(B)的重复测量 ANOVA yB*W+Error(Subject/W)注意：非均衡设计时或存在协变量时，效应项的顺序对结果影响较大，越基础的效应越需要放在表达式前面，首先是协变量、然后是主效应、接着是双因素的交互项，再接着是三因素的交互项。若研究不是正交的，一定要谨慎设置效应的顺序。有三种类型的方法可以分解 yA+B+A:B 右边各效应对 y

10、所解释的方差：类型 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 调整。R 默认调用类型 I 方法，其他软件（比如 SAS 和 SPSS）默认调用类型 III 方法。car 包中的 Anova()函数（不要与标准 anova()函数混淆）提供了使用类型 II 或类

11、型 III 方法的选项，而 aov()函数使用的是类型 I 方法。若想使结果与其他软件（如 SAS 和 SPSS）提供的结果保持一致，可以使用 Anova()函数。三、单因素方差分析1 个因变量，1 个影响因素：总差异 Yij = 平均差异 + 因素差异 i + 随机差异 ij例 1 比较 4 种品牌的胶合板的耐磨性，各抽取 5 个样品，相同转速磨损相同时间测得磨损深度（mm ） ，比较 4 个品牌胶合板的耐磨性有无差异？部分数据如下（ex27_ex1.Rdata）：setwd(“E:/办公资料/R 语言/R 语言学习系列/codes“)load(“ex27_ex1.Rdata“)head(d

12、atas)wear brand1 2.30 A2 2.32 A3 2.40 A4 2.45 A5 2.58 A6 2.35 Battach(datas)table(brand) #各组的样本数brandA B C D 5 5 5 5 aggregate(wear,by=list(brand),mean) #各组均值Group.1 x1 A 2.4102 B 2.4043 C 2.0464 D 2.572aggregate(wear,by=list(brand),sd) #各组标准差Group.1 x1 A 0.112694282 B 0.117601023 C 0.112160604 D 0.

13、03271085library(car)qqPlot(lm(wearbrand,data=datas),simulate=TRUE) #用 Q-Q图检验数据的正态性leveneTest(wearas.factor(brand),data=datas) #方差齐性检验Levenes Test for Homogeneity of Variance (center = median)Df F value Pr(F)group 3 0.6987 0.566416 fitF) brand 3 0.7398 0.24660 24.55 3.15e-06 *Residuals 16 0.1607 0.01

14、005 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1说明：方差齐性检验，原假设 H0：方差齐，p 值=0.56640.05, 故接受原假设，即方差齐。单因素方差分析结果，brand是因素，Residuals 是残差，各列依次为自由度、平方和、均方和、F统计量， p值=3.15e-06|t|) B - A = 0 -0.00600 0.06339 -0.095 0.9997 C - A = 0 -0.36400 0.06339 -5.742 F) supp 1 205.4 205.4 12.317 0.000894 *dose 1 2224.3

15、2224.3 133.415 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 1Error: Plant:concDf Sum Sq Mean Sq F value Pr(F) conc 1 888.6 888.6 215.46 0.000125 *conc:Type 1 239.2 239.2 58.01 0.001595 * Residuals 4 16.5 4.1 -Signif.codes: 0 * 0.001 * 0.0

16、1 * 0.05 . 0.1 1Error: WithinDf Sum Sq Mean Sq F value Pr(F)Residuals 30 869 28.97 说明：在 0.01 的显著水平下，主效应“类型” （p 值=0.00148）和“浓度” （p 值=0.000125）以及交叉效应“类型*浓度” （p 值=0.001595）都非常显著。attach(w1b1)interaction.plot(conc, Type, uptake, type=“b“, col=c(“red“, “blue“), pch=c(16, 18), main=“Interaction Plot for Pl

17、ant Type and Concentration“)boxplot(uptakeType*conc, data=w1b1, col=(c(“gold“, “green“), main=“Chilled Quebec and Mississippi Plants“, ylab = “Carbon dioxide uptake rate (umol/m2 sec)“)detach(w1b1)可以看出，魁北克省的植物比密西西比州的植物二氧化碳吸收率高，而且随着 CO2 浓度的升高，差异越来越明显。注 1：重复测量设计时，需要有长格式数据才能拟合模型，若是宽数据，需要用 reshape2 包中的

18、melt()函数转化为长数据；注 2：这里是用的传统的重复测量方差分析，假设任意组内因子的协方差矩阵为球形，并且任意组内因子两水平间的方差之差都相等。实际中，该假设一般不满足，可以尝试：使用 lme4 包中的 lmer()函数拟合线性混合模型（Bates，2005） ； 使用 car 包中的 Anova()函数调整传统检验统计量以弥补球形假设的不满足（例如 Geisser-Greenhouse 校正） ； 使用 nlme 包中的 gls()函数拟合给定方差- 协方差结构的广义最小二乘模型（UCLA，2009） ； 用多元方差分析对重复测量数据进行建模（Hand ，1987） 。四、多元方差分析

19、1. 因变量不只一个的方差分析，称为多元方差分析。要求数据满足：多元正态性、方差协方差矩阵同质性（即指各组的协方差矩阵相同，通常可用 Boxs M 检验来评估该假设） 。例 4 使用 MASS 包中的 UScereal 数据集，研究美国谷物中的卡路里、脂肪、糖含量是否会因为储物架位置的不同而发生变化。自变量货架位置 shelf 有 1, 2, 3 三个水平。library(MASS)attach(UScereal)yF) shelf 1 0.19594 4.955 3 61 0.00383 *Residuals 63 -Signif.codes: 0 * 0.001 * 0.01 * 0.05

20、 . 0.1 1summary.aov(fit) #输出每个单变量的方差分析结果Response calories :Df Sum Sq Mean Sq F value Pr(F) shelf 1 45313 45313 13.995 0.0003983 *Residuals 63 203982 3238 -Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Response fat :Df Sum Sq Mean Sq F value Pr(F) shelf 1 18.421 18.4214 7.476 0.008108 *Residuals 63 15

21、5.236 2.4641 -Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Response sugars :Df Sum Sq Mean Sq F value Pr(F) shelf 1 183.34 183.34 5.787 0.01909 *Residuals 63 1995.87 31.68 -Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1说明：多元方差分析的 p 值=0.00383F) brand 3 0.7398 0.24660 24.55 3.15e-06 *Residuals 16 0.160

22、7 0.01005 -Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1fit.lm|t|) (Intercept) 2.41000 0.04482 53.768 2e-16 *brandB -0.00600 0.06339 -0.095 0.9258 brandC -0.36400 0.06339 -5.742 3.03e-05 *brandD 0.16200 0.06339 2.556 0.0212 * -Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 0.1

23、002 on 16 degrees of freedomMultiple R-squared: 0.8215, Adjusted R-squared: 0.7881 F-statistic: 24.55 on 3 and 16 DF, p-value: 3.152e-06说明：可以看到 p 值是一样的。线性模型要求预测变量是数值型，当 lm()函数遇到因子型自变量时，会用一系列与因子水平相对应的数值型对照变量来代替因子。若因子有 k 个水平，将会创建 k-1 个对照变量。R 语言中，对照变量的创建方法：contr.helmert第二个水平对照第一个水平，第三个水平对照前两个的均值，第四个水平对

24、照前三个的均值，以此类推；contr.poly基于正交多项式的对照，用于趋势分析（线性、二次、三次等）和等距水平的有序因子；contr.sum对照变量之和限制为 0，也称为偏差找对，对各水平的均值与所有水平的均值进行比较；contr.treatment各水平对照基线水平（默认第一个水平） ，也称作虚拟编码（虚拟变量、哑变量） ；contr.SAS类似于 contr.treatment，只是基线水平变成了最后一个水平，生成的系数类似于大部分 SAS 过程中使用的对照变量。注 1：默认将对照处理（contr.treatment ）用于无序因子，正交多项式（contr.poly）用于有序因子；注 2：可以通过 lm()函数的参数 contrasts=“修改默认的对照方法。contrasts(datas$brand) #查看对照编码，默认第 1 水平（品牌 A）作为参考基准B C DA 0 0 0B 1 0 0C 0 1 0D 0 0 1若样本是品牌 C，则变量 C=1，其它变量 B、D 都=0，无需列出品牌 A 的变量值，因为它的 4 个变量都为 0，就说明是品牌 A.回归模型的参数估计中的变量，brandB 就表示品牌 A 与品牌 B的一个对照，其余也是类似。