《R语言学习系列27-方差分析讲课教案》由会员分享，可在线阅读，更多相关《R语言学习系列27-方差分析讲课教案（21页珍藏版）》请在金锄头文库上搜索。
1、R语言学习系列27-方差分析精品文档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）。当交互影响项呈现统计不显著时，表明各个因素独立，当呈现统计显著时，就需要列出这个交互影响项的效应，以助于作出正确的统计推断。举例解释上述概念：要考察焦虑症的治疗疗效，一个因素是治疗方案，有2种治疗方案，即该因素有2个水平
6、；（治疗方案称为组间因子，因为每个患者只能被分配到一个组别中，没有患者同时接受两种治疗）；再考虑一个因素治疗时间，也有两个水平：治疗5周和治疗6个月，同一患者在5周和6个月不止一次地被测量（两次），称为重复测量（治疗时间称为组内因子，因为每个患者在所有水平下都进行了测量）。建立方差分析模型时，既要考虑两个因素治疗方案和治疗时间（主效应），又要考虑治疗方案和时间的交互影响（交互效应），此时即两因素混合模型方差分析。当某个因素的各个水平下的因变量的均值呈现统计显著性差异时，必要时可作两两水平间的比较，称为均值间的两两比较。 二、R语言实现方差分析对数据的要求：满足正态性（来自同一正态总体）和方差齐
7、性（各组方差相等），在这两个条件下，若各组有差异，则只可能是来自影响因素的不同水平。用aov()函数进行方差分析，基本格式为：aov(formula, data=NULL, projections=FALSE, qr=TRUE,contrasts=NULL, .)其中，formula为方差分析公式；data为数据框；projection设置是否返回预测结果；qr设置是否返回QR分解结果；contrasts为公式中一些因子的列表。formula公式的表示：（y为因变量，ABC为分组因子）符号用法分隔符号，左边为响应变量，右边为解释变量eg：yA+B+C+分隔解释变量：表示变量的交互项eg：yA+
8、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是对被试者独有的标识变量）设计表达式单因素ANOVAyA含单个协变量的单因素ANCOVAyx+A双因素ANOVAyA*B含两个协变量的双因素ANCOVAyx1+x2+A*B随机化区组yB+A, B为区组因子单因素组内ANOVA
9、yA+Error(Subject/A)含单个组内因子(W)和单个组间因子(B)的重复测量ANOVAyB*W+Error(Subject/W)注意：非均衡设计时或存在协变量时，效应项的顺序对结果影响较大，越基础的效应越需要放在表达式前面，首先是协变量、然后是主效应、接着是双因素的交互项，再接着是三因素的交互项。若研究不是正交的，一定要谨慎设置效应的顺序。有三种类型的方法可以分解yA+B+A:B右边各效应对y所解释的方差：类型I（序贯型）效应根据表达式中先出现的效应做调整。A不做调整，B根据A调整，A:B交互项根据A和B调整。类型II（分层型）效应根据同水平或低水平的效应做调整。A根据B调整，B依
10、据A调整，A:B交互项同时根据A和B调整。类型III（边界型）每个效应根据模型其他各效应做相应调整。A根据B和A:B做调整，A:B交互项根据A和B调整。R默认调用类型I方法，其他软件（比如SAS和SPSS）默认调用类型III方法。car包中的Anova()函数（不要与标准anova()函数混淆）提供了使用类型II或类型III方法的选项，而aov()函数使用的是类型I方法。若想使结果与其他软件（如SAS和SPSS）提供的结果保持一致，可以使用Anova()函数。三、单因素方差分析1个因变量，1个影响因素：总差异Yij = 平均差异 + 因素差异i + 随机差异ij例1 比较4种品牌的胶合板的耐磨
11、性，各抽取5个样品，相同转速磨损相同时间测得磨损深度（mm），比较4个品牌胶合板的耐磨性有无差异？部分数据如下（ex27_ex1.Rdata）：setwd(E:/办公资料/R语言/R语言学习系列/codes)load(ex27_ex1.Rdata)head(datas) 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
12、 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.03271085library(car)qqPlot(lm(wearbrand,data=datas),simulate=TRUE) #用Q-Q图检验数据的正态性leveneTest(wearas.factor(brand),data=datas) #方差齐性检验Levenes Test for Homogeneity of Var
13、iance (center = median) Df F value Pr(F)group 3 0.6987 0.5664 16 fitF) brand 3 0.7398 0.24660 24.55 3.15e-06 *Residuals 16 0.1607 0.01005 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1 说明：方差齐性检验，原假设H0：方差齐，p值=0.56640.05, 故接受原假设，即方差齐。单因素方差分析结果，brand是因素，Residuals是残差，各列依次为自由度、平方和、均方和、F统计量，p值=3.15e-060.
14、05, 拒绝原假设，即不同品牌的磨损（均值）有显著差别。library(gplots)plotmeans(wearbrand,xlab=品牌, ylab=磨损) #图形展示带95%置信区间的各组均值通过前面的分析知道，不同品牌的磨损（均值）有显著差别，但并不知道哪个品牌与其它品牌有显著差别。TukeyHSD()函数提供了对各组均值差异的成对检验。TukeyHSD(fit) Tukey multiple comparisons of means 95% family-wise confidence levelFit: aov(formula = wear brand, data = datas)$brand diff lwr upr p adjB-A -0.006 -0.18735345 0.1753535 0.9996826C-A -0.364 -0.54535345 -0.1826465 0.0001610D-A 0.162
