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

整理文档很辛苦,赏杯茶钱您下走!

免费阅读已结束,点击下载阅读编辑剩下 ...

阅读已结束,您可以下载文档离线阅读编辑

资源描述

22.方差分析一、方差分析原理1.方差分析概述方差分析可用来研究多个分组的均值有无差异,其中分组是按影响因素的不同水平值组合进行划分的。方差分析是对总变异进行分析。看总变异是由哪些部分组成的,这些部分间的关系如何。方差分析,是用来检验两个或两个以上均值间差别显著性(影响观察结果的因素:原因变量(列变量)的个数大于2,或分组变量(行变量)的个数大于1)。一元时常用F检验(也称一元方差分析),多元时用多元方差分析(最常用Wilks’∧检验)。方差分析可用于:(1)完全随机设计(单因素)、随机区组设计(双因素)、析因设计、拉丁方设计和正交设计等资料;(2)可对两因素间交互作用差异进行显著性检验;(3)进行方差齐性检验。要比较几组均值时,理论上抽得的几个样本,都假定来自正态总体,且有一个相同的方差,仅仅均值可以不相同。还需假定每一个观察值都由若干部分累加而成,也即总的效果可分成若干部分,而每一部分都有一个特定的含义,称之谓效应的可加性。所谓的方差是离均差平方和除以自由度,在方差分析中常简称为均方(MeanSquare)。2.基本思想基本思想是,将所有测量值上的总变异按照其变异的来源分解为多个部份,然后进行比较,评价由某种因素所引起的变异是否具有统计学意义。根据效应的可加性,将总的离均差平方和分解成若干部分,每一部分都与某一种效应相对应,总自由度也被分成相应的各个部分,各部分的离均差平方除以各自的自由度得出各部分的均方,然后列出方差分析表算出F检验值,作出统计推断。方差分析的关键是总离均差平方和的分解,分解越细致,各部分的含义就越明确,对各种效应的作用就越了解,统计推断就越准确。效应项与试验设计或统计分析的目的有关,一般有:主效应(包括各种因素),交互影响项(因素间的多级交互影响),协变量(来自回归的变异项),等等。当分析和确定了各个效应项S后,根据原始观察资料可计算出各个离均差平方和SS,再根据相应的自由度df,由公式MS=SS/df,求出均方MS,最后由相应的均方,求出各个变异项的F值,F值实际上是两个均方之比值,通常情况下,分母的均方是误差项的均方。根据F值的分子、分母均方的自由度f1和f2,在确定显著性水平为α情况下,由F(f1,f2)临界值表查得单侧Fα界限值。当FFα时,则P值α,不拒绝原假设H0,说明不拒绝这个效应项的效应为0的原假设,也即这个效应项是可能对总变异没有实质影响的;若FFα则P值≤α,拒绝原假设H0,也即这个效应项是很可能对总变异有实质影响的。3.方差分析的实验设计为了确定方差分析表中各个有关效应项,需要在试验设计阶段就作出安排,再根据设计要求进行试验,得出原始观察值,按原来设计方案算出方差分析表中的各项。在试验设计阶段通常需要考虑如下4个方面:(1)研究的因变量即试验所要观察的主要指标,一次试验时可以有多个观察指标,方差分析时也可以同时对多个因变量进行分析;(2)因素和水平试验的因素(factor)可以是品种、人员、方法、时间、地区等等,因素所处的状态叫水平(level)。在每一个因素下面可以分成若干水平。(3)因素间的交互影响多因素的试验设计,有时需要分析因素间的交互影响(interaction),2个因素间的交互影响称为一级交互影响(A×B);3个因素间的交互影响称为二级交互影响(A×B×C)。当交互影响项呈现统计不显著时,表明各个因素独立,当呈现统计显著时,就需要列出这个交互影响项的效应,以助于作出正确的统计推断。举例解释上述概念:要考察焦虑症的治疗疗效,一个因素是治疗方案,有2种治疗方案,即该因素有2个水平;(治疗方案称为组间因子,因为每个患者只能被分配到一个组别中,没有患者同时接受两种治疗);再考虑一个因素治疗时间,也有两个水平:治疗5周和治疗6个月,同一患者在5周和6个月不止一次地被测量(两次),称为重复测量(治疗时间称为组内因子,因为每个患者在所有水平下都进行了测量)。建立方差分析模型时,既要考虑两个因素治疗方案和治疗时间(主效应),又要考虑治疗方案和时间的交互影响(交互效应),此时即两因素混合模型方差分析。当某个因素的各个水平下的因变量的均值呈现统计显著性差异时,必要时可作两两水平间的比较,称为均值间的两两比较。二、R语言实现方差分析对数据的要求:满足正态性(来自同一正态总体)和方差齐性(各组方差相等),在这两个条件下,若各组有差异,则只可能是来自影响因素的不同水平。用aov()函数进行方差分析,基本格式为:aov(formula,data=NULL,projections=FALSE,qr=TRUE,contrasts=NULL,...)其中,formula为方差分析公式;data为数据框;projection设置是否返回预测结果;qr设置是否返回QR分解结果;contrasts为公式中一些因子的列表。formula公式的表示:(y为因变量,ABC为分组因子)符号用法~分隔符号,左边为响应变量,右边为解释变量eg:y~A+B+C+分隔解释变量:表示变量的交互项eg:y~A+B+A:B*表示所有可能交互项eg:y~A*B*C可展开为:y~A+B+C+A:B+A:C+B:C+A:B:C^表示交互项达到次数eg:y~(A+B+C)^2展开为:y~A+B+C+A:B+A:C+B:C.表示包含除因变量外的所有变量eg:若一个数据框包括变量y,A、B和C,代码y~.可展开为y~A+B+C常见研究设计的表达式:(小写字母表示定量变量,大写字母表示组别因子,Subject是对被试者独有的标识变量)设计表达式单因素ANOVAy~A含单个协变量的单因素ANCOVAy~x+A双因素ANOVAy~A*B含两个协变量的双因素ANCOVAy~x1+x2+A*B随机化区组y~B+A,B为区组因子单因素组内ANOVAy~A+Error(Subject/A)含单个组内因子(W)和单个组间因子(B)的重复测量ANOVAy~B*W+Error(Subject/W)注意:非均衡设计时或存在协变量时,效应项的顺序对结果影响较大,越基础的效应越需要放在表达式前面,首先是协变量、然后是主效应、接着是双因素的交互项,再接着是三因素的交互项。若研究不是正交的,一定要谨慎设置效应的顺序。有三种类型的方法可以分解y~A+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调整。R默认调用类型I方法,其他软件(比如SAS和SPSS)默认调用类型III方法。car包中的Anova()函数(不要与标准anova()函数混淆)提供了使用类型II或类型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(datas)wearbrand12.30A22.32A32.40A42.45A52.58A62.35Battach(datas)table(brand)#各组的样本数brandABCD5555aggregate(wear,by=list(brand),mean)#各组均值Group.1x1A2.4102B2.4043C2.0464D2.572aggregate(wear,by=list(brand),sd)#各组标准差Group.1x1A0.112694282B0.117601023C0.112160604D0.03271085library(car)qqPlot(lm(wear~brand,data=datas),simulate=TRUE)#用Q-Q图检验数据的正态性leveneTest(wear~as.factor(brand),data=datas)#方差齐性检验Levene'sTestforHomogeneityofVariance(center=median)DfFvaluePr(F)group30.69870.566416fit-aov(wear~brand,data=datas)#单因素方差分析(检验组间差异)summary(fit)DfSumSqMeanSqFvaluePr(F)brand30.73980.2466024.553.15e-06***Residuals160.16070.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.05,拒绝原假设,即不同品牌的磨损(均值)有显著差别。library(gplots)plotmeans(wear~brand,xlab=品牌,ylab=磨损)#图形展示带95%置信区间的各组均值通过前面的分析知道,不同品牌的磨损(均值)有显著差别,但并不知道哪个品牌与其它品牌有显著差别。TukeyHSD()函数提供了对各组均值差异的成对检验。TukeyHSD(fit)Tukeymultiplecomparisonsofmeans95%family-wiseconfidencelevelFit:aov(formula=wear~brand,data=datas)$branddifflwruprpadjB-A-0.006-0.187353450.17535350.9996826C-A-0.364-0.54535345-0.18264650.0001610D-A0.162-0.019353450.34335350.0886142C-B-0.358-0.53935345-0.17664650.0001929D-B0.168-0.013353450.34935350.0744337D-C0.5260.344646550.70735350.0000019说明:可以看出(H0:无差异),B与A的差异非常不显著,C与A、C与B、D与C的差异非常显著。multcomp包中的glht()函数提供了更为全面的多重均值比较方法。library(multcomp)attach(datas)tuk-glht(fit,linfct=mcp(brand=Tukey))#注意datas$brand必须是因子型summary(tuk)SimultaneousTestsforGeneralLinearHypothesesMultipleComparisonsofMeans:TukeyContrastsFit:aov(formula=wear~brand,data=datas)LinearHypotheses:EstimateStd.ErrortvaluePr(|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

1 / 21
下载文档,编辑使用

©2015-2020 m.777doc.com 三七文档.

备案号:鲁ICP备2024069028号-1 客服联系 QQ:2149211541

×
保存成功