R语言学习系列27方差分析Word文档下载推荐.docx

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

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

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

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

5、着是双因素的交互项，再接着是三因素的交互项。若研究 不是正交的，一定要谨慎设置效应的顺序。有三种类型的方法可以分解 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:R默认调用类型 I 方法，其他软件（比如 SAS和 SPSS）默

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

7、x1.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#各组标准差D 2.572 aggregate（wear,by=list（

8、brand）,sd）1A 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（F） group 3 0.6987 0.566416fitF）brand 3 0.

9、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, 拒绝 原假设，即不同品牌的磨损（均值）有显著差别。library（gplots）plotmeans（wearbrand,xlab= 品牌,

10、 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 -0.006 -0.18735345 0.1753535 0.9996826C-A

11、 -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（） 函数提供了更为全面的多重均值比较方法。 library（mult

12、comp） attach（datas）tuk |t|）B -A = 0-0.006000.06339-0.0950.9997- A =-0.36400-5.7420.001 *D0.162002.5560.0886 .- B =-0.35800-5.6480.168002.6500.0743 .- C =0.526008.298 0.001 * 0.01 * 0.05（Adjusted p values reported - single-step method）plot（cld（tuk, level = 0.05）, col = lightgrey）说明：标记相同字母（标记 b）的品牌 AB

13、D认为是无显著差异，在 同一亚组，而品牌 C（标记 a）与另外三个品牌有显著差异。另外，也可以进行多重 t 检验，使用函数：pairwise.t.test（x, g, p.adjust.method=,.）其中， x 为因变量， g 为因子型的分组变量；p.adjust.method 设置 p 值的修正方法，由于多次重复 t 检验会 大大增加犯第一类错误的概率，为此要进行 p 值的修正，使用 bonferroni 法修正效果较好。pairwise.t.test（wear,brand,p.adjust.method=bonferroniPairwise comparisons using t t

14、ests 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无差异， C与 ABD有显著 差异。最后，方差分析对离群点非常敏感，检验是否有离群点：outlierTest（fit）No Studentized residuals with Bonferonni p Largest |rstudent|:rstudent unadjusted p

15、-value Bonferonni p9 2.528103 0.023182 0.46364经检验无离群点。三、两因素方差分析1 个因变量， 2 个影响因素：总差异 Yijk = 平均差异 + 因素 1 差异 i + 因素 2 差异 i+ 因素 1,2 交互作用差异 ij + 随机差异 ijk例 2 研究 60 只豚鼠的牙齿生长数据，按 2 种喂食方法：橙汁、维生 素 C，各喂食方法中抗坏血酸含量都有 3 个水平： 0.5mg、1mg、2mg，分配为 6 组，每组各 10 只，牙齿长度为因变量。 做两因素方差分析。attach（ToothGrowth） head（ToothGrowth）len

16、 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 10VC 10 10 10 aggregate（len, by=list（supp, dose）, mean）Group.1 Group.2 x#计算各组均值1OJ0.5 13.232VC0.5 7.9831.0 22.701.0 16.7752.0 26.0662.0 26.14#计算各组标准差0.5 4.459

17、7090.5 2.7466341.0 3.9109531.0 2.5153092.0 2.6550582.0 4.797731aggregate（len, by=list（supp, dose）, sd）bartlett.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=

18、ToothGrowth） #关于因素 dose 的方差齐性检验 len by doses K-squared = 0.66547, df = 2, p-value = 0.717 fit-aov（lensupp*dose,data=ToothGrowth） #做两因素方差分析，考虑全部效应supp 1 205.4 205.4 12.317 0.000894 *dose 1 2224.3 2224.3 133.415 2e-16 *Signif.codes:可以看出，主效应 supp和 dose都非常显著（ p 值都远小于 0.05 ），交互效应也显著（p 值=0.02460.05 ）。若交互作

19、用不显著，可以可以只做去掉交互效应的方差分析。图形化展示两因素方差分析的交互效应：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 Suppinteraction.plot（supp, dose, len, type=, col=c（Interaction between Supp and Dose有一个图的线有交叉，说明有交互作用。可以看出随着橙 汁和维生素 C中的抗坏血酸剂量的增加，牙

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

21、e（组间因子）为植物类型，自变量 conc（组内因子）为七种水平的 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） #先只考虑寒带植物-aov（uptake（conc*Type）+Error（Plant/conc）, data=w1b1） summary（fit）Error: PlantType 1 2667.2 2667.2 60.41 0.00148 *Residuals 4 176.6 44.1 Plant:concconc 1 888.6 888.6 215.46 0.000125 *conc:Type 1Residuals 4239.216.54.158.01 0.001595 *0.001 * 0.01* 0.