1、R语言学习系列27方差分析R语言学习系列27-方差分析22. 方差分析一、方差分析原理1. 方差分析概述方差分析可用来研究多个分组的均值有无差异,其中分组是按影响因素的不同水平值组合进行划分的。方差分析是对总变异进行分析。看总变异是由哪些部分组成的,这些部分间的关系如何。方差分析,是用来检验两个或两个以上均值间差别显著性(影响观察结果的因素:原因变量(列变量)的个数大于2,或分组变量(行变量)的个数大于1)。一元时常用F检验(也称一元方差分析),多元时用多元方差分析(最常用Wilks检验)。方差分析可用于:(1)完全随机设计(单因素)、随机区组设计(双因素)、析因设计、拉丁方设计和正交设计等资
2、料;(2)可对两因素间交互作用差异进行显著性检验;(3)进行方差齐性检验。要比较几组均值时,理论上抽得的几个样本,都假定来自正态总体,且有一个相同的方差,仅仅均值可以不相同。还需假定每一个观察值都由若干部分累加而成,也即总的效果可分成若干部分,而每一部分都有一个特定的含义,称之谓效应的可加性。所谓的方差是离均差平方和除以自由度,在方差分析中常简称为均方(Mean Square)。计推断。举例解释上述概念:要考察焦虑症的治疗疗效,一个因素是治疗方案,有2种治疗方案,即该因素有2个水平;(治疗方案称为组间因子,因为每个患者只能被分配到一个组别中,没有患者同时接受两种治疗);再考虑一个因素治疗时间,
3、也有两个水平:治疗5周和治疗6个月,同一患者在5周和6个月不止一次地被测量(两次),称为重复测量(治疗时间称为组内因子,因为每个患者在所有水平下都进行了测量)。建立方差分析模型时,既要考虑两个因素治疗方案和治疗时间(主效应),又要考虑治疗方案和时间的交互影响(交互效应),此时即两因素混合模型方差分析。当某个因素的各个水平下的因变量的均值呈现统计显著性差异时,必要时可作两两水平间的比较,称为均值间的两两比较。 二、R语言实现方差分析对数据的要求:满足正态性(来自同一正态总体)和方差齐性(各组方差相等),在这两个条件下,若各组有差异,则只可能是来自影响因素的不同水平。用aov()函数进行方差分析,
4、基本格式为: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+B+A:B*表示所有可能交互项eg:yA*B*C可展开为:yA+B+C+A:B+A:C+B:C+A:B:C表示交互
5、项达到次数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为区组因子单因素组内ANOVAyA+Error(Subject/A)含单个组内因子(W)和单个组间因子(B)的重复测量ANOVAyB*W+Err
6、or(Subject/W)注意:非均衡设计时或存在协变量时,效应项的顺序对结果影响较大,越基础的效应越需要放在表达式前面,首先是协变量、然后是主效应、接着是双因素的交互项,再接着是三因素的交互项。若研究不是正交的,一定要谨慎设置效应的顺序。有三种类型的方法可以分解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:
7、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个品牌胶合板的耐磨性有无差异?部分数据如下(e
8、x27_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 A 2.4102 B 2.4043 C 2.0464 D 2.572aggregate(wear,by=list
9、(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 Variance (center = median) Df F value Pr(F)group 3 0.6987 0.
10、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.05, 拒绝原假设,即不同品牌的磨损(均值)有显著差别。library(gplots)plotmeans(wear
11、brand,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
12、 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的差异非常显著。multcomp包中的glht()函数提供了更为全面的多重均值比较方法
13、。library(multcomp)attach(datas)tuk |t|) B - A = 0 -0.00600 0.06339 -0.095 0.9997 C - A = 0 -0.36400 0.06339 -5.742 0.001 *D - A = 0 0.16200 0.06339 2.556 0.0886 . C - B = 0 -0.35800 0.06339 -5.648 0.001 *D - B = 0 0.16800 0.06339 2.650 0.0743 . D - C = 0 0.52600 0.06339 8.298 0.001 *-Signif. codes:
14、0 * 0.001 * 0.01 * 0.05 . 0.1 1(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检验会大大增
15、加犯第一类错误的概率,为此要进行p值的修正,使用bonferroni法修正效果较好。pairwise.t.test(wear,brand,p.adjust.method=bonferroni) Pairwise comparisons using t tests with pooled SD data: wear and brand A B C B 1.00000 - - C 0.00018 0.00022 - D 0.12695 0.10474 2.1e-06P value adjustment method: bonferroni 说明:原假设H0: 无差异,可见A与B无差异,C与ABD有
16、显著差异。最后,方差分析对离群点非常敏感,检验是否有离群点: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 + 随机差异ijk例2 研究60只豚鼠的牙齿生长数据,
17、按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 2 OJ 10 10 10 VC 10 10 10aggregate(len, by=list
18、(supp, dose), mean) #计算各组均值 Group.1 Group.2 x1 OJ 0.5 13.232 VC 0.5 7.983 OJ 1.0 22.704 VC 1.0 16.775 OJ 2.0 26.066 VC 2.0 26.14aggregate(len, by=list(supp, dose), sd) #计算各组标准差 Group.1 Group.2 x1 OJ 0.5 4.4597092 VC 0.5 2.7466343 OJ 1.0 3.9109534 VC 1.0 2.5153095 OJ 2.0 2.6550586 VC 2.0 4.797731bart
19、lett.test(lensupp,data=ToothGrowth) #关于因素supp的方差齐性检验 Bartlett test of homogeneity of variancesdata: len by suppBartletts K-squared = 1.4217, df = 1, p-value = 0.2331bartlett.test(lendose,data=ToothGrowth) #关于因素dose的方差齐性检验 Bartlett test of homogeneity of variancesdata: len by doseBartletts K-squared
20、= 0.66547, df = 2, p-value = 0.717fitF) supp 1 205.4 205.4 12.317 0.000894 *dose 1 2224.3 2224.3 133.415 2e-16 *supp:dose 1 88.9 88.9 5.333 0.024631 * Residuals 56 933.6 16.7 -Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1说明:可以看出,主效应supp和dose都非常显著(p值都远小于0.05),交互效应也显著(p值=0.02460.05)。若交互作用不显著,可以可以只做去掉交
21、互效应的方差分析。图形化展示两因素方差分析的交互效应:par(mfrow=c(1,2)interaction.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) 说明:有一个图的线有交叉,说明有交互作用。可
22、以看出随着橙汁和维生素C中的抗坏血酸剂量的增加,牙齿长度变长;。对于0.5 mg和1 mg剂量,橙汁比维生素C更能促进牙齿生长;对于2 mg剂量的抗坏血酸,两种喂食方法下牙齿长度增长相同。也可以用HH包中的interaction2wt()函数(也适合三因素方差分析)来展示更全面的可视化结果:library(HH)interaction2wt(lensupp*dose)三、重复测量方差分析 重复测量方差分析,即受试者被测量不止一次。例3(1个组内1个组间因子的重复测量)在某浓度CO2的环境中,对寒带植物(来自魁北克)和非寒带植物的(来自密西西比)光合作用率进行比较。因变量uptake为CO2吸收
23、量,自变量Type(组间因子)为植物类型,自变量conc(组内因子)为七种水平的CO2浓度。attach(CO2)head(CO2) #注意CO2是长格式的数据 Plant Type Treatment conc uptake1 Qn1 Quebec nonchilled 95 16.02 Qn1 Quebec nonchilled 175 30.43 Qn1 Quebec nonchilled 250 34.84 Qn1 Quebec nonchilled 350 37.25 Qn1 Quebec nonchilled 500 35.36 Qn1 Quebec nonchilled 675
24、39.2w1b1-subset(CO2, Treatment=chilled) #先只考虑寒带植物fitF) 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:conc Df 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 1
25、6.5 4.1 -Signif.codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Error: Within Df 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
26、=Interaction Plot for Plant 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:重复测量设计时,需要有长格式数据才能拟合模型,若是宽数据,需要用res
27、hape2包中的melt()函数转化为长数据;注2:这里是用的传统的重复测量方差分析,假设任意组内因子的协方差矩阵为球形,并且任意组内因子两水平间的方差之差都相等。实际中,该假设一般不满足,可以尝试:使用lme4包中的lmer()函数拟合线性混合模型(Bates,2005); 使用car包中的Anova()函数调整传统检验统计量以弥补球形假设的不满足(例如Geisser-Greenhouse校正); 使用nlme包中的gls()函数拟合给定方差-协方差结构的广义最小二乘模型(UCLA,2009); 用多元方差分析对重复测量数据进行建模(Hand,1987)。四、多元方差分析1. 因变量不只一个
28、的方差分析,称为多元方差分析。要求数据满足:多元正态性、方差协方差矩阵同质性(即指各组的协方差矩阵相同,通常可用Boxs M检验来评估该假设)。例4 使用MASS包中的UScereal数据集,研究美国谷物中的卡路里、脂肪、糖含量是否会因为储物架位置的不同而发生变化。自变量货架位置shelf有1, 2, 3三个水平。library(MASS)attach(UScereal)y-cbind(calories, fat, sugars)head(y) calories fat sugars1, 212.1212 3.030303 18.181822, 212.1212 3.030303 15.151513, 100.0000 0.000000 0.000004, 146.6667 2.666667 13.333335, 110.0000 0.000000 14.000006, 173.3333 2.666667 10.66667aggregate(y, by=list(shelf), mean) Group.
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1