《R语言实验九.doc》由会员分享，可在线阅读，更多相关《R语言实验九.doc（15页珍藏版）》请在三一办公上搜索。
1、精选优质文档-倾情为你奉上实验九 单因素和双因素方差分析、协方差分析【实验类型】验证性【实验学时】2 学时【实验目的】1、掌握方差分析的基本原理、数学模型及检验；2、掌握方差分析与线性回归的联系；3、掌握协方差分析的基本原理及求解方法。【实验内容】第一部分、课件例题：#例9.10 的求解：x-scan(F:/文档/大学课程/R语言/ch09/competition.dat, n=28) # 读取共sum(1:7)=28 个数据names-scan(F:/文档/大学课程/R语言/ch09/competition.dat, what= , skip=7)# 读取 8 个变量R-matrix(1,
2、nrow=8, ncol=8, dimnames=list(names, names)for (i in 2:8) for (j in 1:(i-1) Ri,j-x(i-2)*(i-1)/2+j; Rj,i-Ri,j # 生成相关矩阵fa-factanal(factors=2, covmat=R) # 因子数为 2fa分析：Uniquenesses为特殊方差；Loadings为因子载荷矩阵；SS loadings为公共因子方对变量的总方差贡献；Proportion Var为方差贡献率；Cumulative Var为累积方差贡献率。在计算结果中，因子 f 1 后几个变量 (X 6 , X 7 ,
3、 X 8 ) 的载荷因子接近于 l ，这些变量涉及的是长跑，因此可称 f l 是耐力因子耐力因子；而因子 f 2 中前几个变量 (X l , X 2 ) 接近 1 ，涉及的是短跑，因此可称 f 2 是 速度因子。#例9.10 的求解：# 不选择旋转，先计算出载荷矩阵后，再使用varimax()函数作最大方差旋转函数作最大方差旋转fac-factanal(factors=2, covmat=R, rotation=none)varm-varimax(fac$loadings);varmt(varm$rotmat)%*%varm$rotmat # 旋转矩阵为正交阵 得到的结果与前面直接计算的结果完
4、全相同；此时的旋转矩阵 $rotmat 仍为正交矩阵。#例9.10 的求解：# 作斜交变换pro-promax(fac$loadings);prot(pro$rotmat)%*%pro$rotmat # 旋转矩阵不再是正交阵solve(t(pro$rotmat)%*%pro$rotmat) # 逆矩阵为旋转因子的相关矩阵 此时的旋转矩阵 $rotmat 不再是正交矩阵。#例9.11 的求解：rt-read.table(F:/文档/大学课程/R语言/ch09/employ.dat) # 读取数据fa-factanal(., factors=5, data=rt); fa # 选取 5 个因子#例
5、9.12 的求解：fa-factanal(., factors=5, data=rt, scores=regression) #使用回归方法计算因子得分fa$scores # 因子得分plot(fa$scores,1:2,type=n) # 画出应聘者在第一、第二公共因子下的散点图text(fa$scores,1,fa$scores,2)#例 9.13 的求解：test-read.table(F:/文档/大学课程/R语言/ch09/recovery club.dat)test-scale(test) #将数据标准化ca-cancor(test,1:3,test,4:6)ca# 例 9.13 的
6、求解：# 计算样本数据在典型变量下的得分U-as.matrix(test,1:3)%*%ca$xcoef;UV-as.matrix(test,4:6)%*%ca$ycoef;V# 画出相关变量 U1,V1 和 U3,V3 为坐标的数据散点图plot(U,1, V,1, xlab=U1, ylab=V1)plot(U,3, V,3, xlab=U3, ylab=V3)左图为第 1 典型变量的散点图，其相关系数为 0.796 ，接近于 1 ，所以在一直线附近；而右图是第 3 典型变量的散点图，其相关系数为 0.0726 ，接近于 0 ，所以很分散。因此，我们只利用第 1 典型变量分析问题，达到降维
7、的目的。#例 9.14 的求解：data-data.frame(很不满意=c(42,35,13,7,3),有些不满=c(82,62,28,18,7),比较满意 =c(67,165,92,54,32),很满意=c(55,118,81,75,54),row.names=c(10 万 ) )datachisq.test(data) # 卡方检验由于 P 值 0.05 ，拒绝原假设，即认为收入与满意度之间有密切联系，可以进一步作对应分析。#例 9.14 的求解：library(MASS)cra-corresp(data,nf=2) # 对应分析，因子个数为 2crabiplot(cra) # 绘制对应
8、图abline(v=0,h=0,lty=3)结果表明 : 收入在 1 万以下的人对自己的职业有些不满或很不满意 ; 收入在 15 万的人对自己的职业感到比较满意 ;而收入在 5 万以上的人大都对自己的职业感到很满意。第二部分、教材例题：#例 8.1.1X-c(25.6, 22.2, 28.0, 29.8, 24.4, 30.0, 29.0, 27.5, 25.0, 27.7, 23.0, 32.2, 28.8, 28.0, 31.5, 25.9, 20.6, 21.2, 22.0, 21.2)A-factor(rep(1:5, each=4)miscellany-data.frame(X, A
9、)aov.misF)表示检验的p值; A就是因素A;Residuals为残差. p=0.016180.05,说明有理由拒绝原假设, 即认为五种除杂方法有显著差异#通过函数plot( )绘图可直观描述5种不同除杂方法之间的差异plot(miscellany$Xmiscellany$A)#例8.1.3sales-data.frame( X=c(23, 19, 21, 13, 24, 25, 28, 27, 20, 18,19, 15, 22, 25, 26, 23, 24, 23, 26, 27), A=factor(rep(1:5, c(4, 4, 4, 4, 4)summary(aov(XA,
10、 sales) #其次进行方差分析TukeyHSD(aov(XA, sales) #求均值之差的同时置信区间可见不同的销售方式有差异可以看出, 共有10个两两比较的结果, A 3 ? A 1 、A 4 ? A 2 、A 5 ? A 2 和A 5 ? A 4 的差异是显著的, 其它两两比较的结果均是不显著的.第三部分、课后习题：8.1X-c(38.7,41.5 ,43.8 ,44.5, 45.5 ,46.0 ,47.7, 58.0,39.2, 39.3, 39.7 ,41.4 ,41.8 ,42.9, 43.3, 45.8,34.0 ,35.0, 39.0 ,40.0 ,43.0, 43.0,
11、44.0 ,45.0,34.0 ,34.8 ,34.8 ,35.4 ,37.2 ,37.8 ,41.2 ,42.8)miscellany-data.frame(X, A)aov.mis-aov(XA, data=miscellany)summary(aov.mis)plot(miscellany$Xmiscellany$A)p=0.002660.05,说明有理由拒绝原假设, 即认为各个实验室生产的纸张的光滑度有显著差异8.2A-c(20.4,30.2,210.4,365.0,56.8,37.8,265.3,175.0,169.8,356.4,254.0,262.3,170.5,360.0,78
12、.4,86.4,128.0,24.1,28.5,108.5,472.5,158.6,238.7,253.6,57.0,189.6,59.3,259.3,380.2,210.5,64.6,87.3)B-c(281.0,377.1,230.0,537.9,248.7,571.4,766.2,495.0,87.3,389.8,423.9,577.3,66.8,521.3,327.8,421.4,149.7,47.5,425.7,270.8,378.5,228.0,538.7,245.6,584.1,64.8,485.6,110.8,398.7,452.6,587.7,86.8,532.1,311.6
13、,442.2)C-c(480.0,488.9,350.7,652.8,1400.0,850.0,725.6,590.0,765.0,1200.0,231.2,485.3,600.0,1380.0,438.5,652.4,432.8,286.1,464.8,608.4,688.5,630.5,750.0,815.0,664.0,348.6,550.0,640.0)X - c(A,B,C)a=factor(rep(1:3,c(32,35,28)ai - data.frame(X,a)#显著性检验summary(aov(Xa,ai)TukeyHSD(aov(Xa, ai)#正态shapiro.tes
14、t(X)#方差齐性检验bartlett.test(Xa, data=ai)#bartlettlibrary(carData)library(car)leveneTest(ai$X, ai$a) 由于方差齐性检验的P值（= 4.792e-05）0.05，因此可以认为三个群体中CEA含量的方差具有显著差异。8.3A1 - c(1073,1058,1071,1037,1066,1026,1053,1049,1065,1051) #鱼粉A2 - c(1061,1058,1038,1042,1020,1045,1044,1061,1034,1049) #槐树粉A3 - c(1084,1069,1160,
15、1078,1075,1090,1079,1094,1111,1092) #苜蓿A - c(A1,A2,A3)a - factor(rep(1:3,each=10)miscellany-data.frame(A,a)aov.misF)表示检验的p值; a就是因素a;Residuals为残差. p= 1.16e-050.05,说明有理由拒绝原假设, 即认为三组小鸡体重有显著差异。8.4tanxin - data.frame( X - c(71.73,73.75,76.73,71.73,72.73,76.74,79.77,73.72,75.73,78.77,74.75,70.71,77.75,76.
16、74,74.73,69.69), A=gl(4,4), B=gl(4,1,16)tanxin.aov-aov(XA+B, data=tanxin) summary(tanxin.aov)#方差齐性检验bartlett.test(XA, data=tanxin) # 对因素A bartlett.test(XB, data=tanxin) #对因素B结论: p值说明因素B对纤维弹性有影响, 而没有充分理由说明因素A对纤维弹性有影响.结论:对因素A, p值(0.9243)远大于0.05, 接受原假设, 认为因素A的各水平下的数据是等方差的; 对因素B, p值(0.8905)大于0.05, 接受原假设
17、, 认为因素B的各水平下的数据是等方差的。8.5s - data.frame( X - c(19.3,19.2,24.0,27.3,26.0,28.5,27.8,28.5,21.7,22.6,27.5,30.3,29.0,28.7,30.2,29.8,20.0,20.1,24.2,27.3,24.5,27.1,28.1,27.7), A=gl(3,8), B=gl(4,2,24)op-par(mfrow=c(1, 2)plot(XA+B, data=s)上图显示两因素的各水平均有较大差异存在.#用函数interaction.plot( )作出交互效应图, 以考查因素之间交互作用是否存在with
18、(s,interaction.plot(A, B, X, trace.label=B) with(s,interaction.plot(B, A, X, trace.label=A)由于两图中的曲线有明显的相交情况出现, 因此我们初步认为两个因素有交互作用。#方差分析考虑交互作用s.aov-aov(XA*B, data=s) summary(s.aov) 根据p值知, 因素A和B对Y的影响是高度显著的, 而交互作用对Y的影响却不显著.8.6d - data.frame( X - c(52,43,39,41,47,53,49,38,42,48,37,39,50,41,30,36,48,47,34
19、,42,38,36,39,44,37,40,32,45,58,42,44,46,60,43,56,41), A=gl(4,9,36), B=gl(3,3,36)#箱线图op-par(mfrow=c(1, 2)plot(XA+B, data=d)#方差分析考虑交互作用d.aov-aov(XA*B, data=d) summary(d.aov)根据p值知, 因素A对Y的影响是高度显著的, 而因素B与交互作用对Y的影响却不显著.8.7s - data.frame(X-c(4.6,4.3,6.1,6.5,6.8,6.4,6.3,6.7,3.4,3.8,4.0,3.8,4.7,4.3,3.9,3.5,6
20、.5,7.0), A=gl(3,6,18), B=gl(3,2,18)#方差分析考虑交互作用s.aov-aov(XA*B, data=s) summary(s.aov)根据p值知, 因素A、B与其交互作用对产量的影响是高度显著的.#均值的多重比较a=factor(rep(1:9,each=2)pairwise.t.test(X, a, p.adjust.method=none)pairwise.t.test(X, a, p.adjust.method=holm)pairwise.t.test(X,a,p.adjust.method=bonferroni)9与其它8个差异明显，后者差异不明显作调
21、整后p值增大, 在一定程度上克服了多重t检验的缺点.8.8h - rep(c(wu,30,300),each=7)X - c(1.40,1.5,1.8,2.2,3.4,3.6,4.6,1.6,2.0,2.3,2.9,4.5,5.1,6.0,2.2,2.3,3.0,3.2,4.5,5.9,7.0)Y - c(0,3,5,10,2,25,30,0,3,5,10,20,25,25,0,3,10,5,20,30,30)m - data.frame(h,X,Y)library(HH)ancova(Y X+h , data=m)可见胸腺嘧啶核苷X的含量对DNA的含成量Y有显著差异8.9r - rep(c(
22、白人,黑人,西班牙人,亚洲人),each=6)X - c(260,275,278,280,282,288,260,263,270,278,281,285,262,264,270,275,280,284,260,271,274,279,281,283)Y - c(130,135,138,142,146,149,115,118,120,125,128,132,113,115,120,121,127,132,111,174,117,118,120,122)m - data.frame(r,X,Y)library(HH)ancova(Y X+r , data=m)ancova(Y X*r , data=m)说明不同种族对出生体重的影响是显著的。专心-专注-专业
