R语言实验九

实验九单因素和双因素方差分析、协方差分析【实验类型】验证性【实验学时】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,nrow=8,ncol=8,dimnames=list(names,names))for(iin2:8){for(jin1:(i-1)){R[i,j]-x[(i-2)*(i-1)/2+j];R[j,i]-R[i,j]}}#生成相关矩阵fa-factanal(factors=2,covmat=R)#因子数为2fa分析：Uniquenesses为特殊方差；Loadings为因子载荷矩阵；SSloadings为公共因子方对变量的总方差贡献；ProportionVar为方差贡献率；CumulativeVar为累积方差贡献率。在计算结果中，因子f1后几个变量(X6,X7,X8)的载荷因子接近于l，这些变量涉及的是长跑，因此可称fl是耐力因子耐力因子；而因子f2中前几个变量(Xl,X2)接近1，涉及的是短跑，因此可称f2是速度因子。#例9.10的求解：#不选择旋转，先计算出载荷矩阵后，再使用varimax()函数作最大方差旋转函数作最大方差旋转fac-factanal(factors=2,covmat=R,rotation=none)varm-varimax(fac$loadings);varmt(varm$rotmat)%*%varm$rotmat#旋转矩阵为正交阵得到的结果与前面直接计算的结果完全相同；此时的旋转矩阵$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个因子#例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/recoveryclub.dat)test-scale(test)#将数据标准化ca-cancor(test[,1:3],test[,4:6])ca#例9.13的求解：#计算样本数据在典型变量下的得分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典型变量分析问题，达到降维的目的。#例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(1万,1~3万,3~5万,5~10万,10万))datachisq.test(data)#卡方检验由于P值0.05，拒绝原假设，即认为收入与满意度之间有密切联系，可以进一步作对应分析。#例9.14的求解：library(MASS)cra-corresp(data,nf=2)#对应分析，因子个数为2crabiplot(cra)#绘制对应图abline(v=0,h=0,lty=3)结果表明:收入在1万以下的人对自己的职业有些不满或很不满意;收入在1~5万的人对自己的职业感到比较满意;而收入在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)aov.mis-aov(X~A,data=miscellany)summary(aov.mis)说明:上述结果中,Df表示自由度;sumSq表示平方和;MeanSq表示均方和;Fvalue表示F检验统计量的值,即F比;Pr(F)表示检验的p值;A就是因素A;Residuals为残差.p=0.016180.05,说明有理由拒绝原假设,即认为五种除杂方法有显著差异#通过函数plot()绘图可直观描述5种不同除杂方法之间的差异plot(miscellany$X~miscellany$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(X~A,sales))#其次进行方差分析TukeyHSD(aov(X~A,sales))#求均值之差的同时置信区间可见不同的销售方式有差异可以看出,共有10个两两比较的结果,A3?A1、A4?A2、A5?A2和A5?A4的差异是显著的,其它两两比较的结果均是不显著的.第三部分、课后习题：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,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(X~A,data=miscellany)summary(aov.mis)plot(miscellany$X~miscellany$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.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,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(X~a,ai))TukeyHSD(aov(X~a,ai))##正态shapiro.test(X)#方差齐性检验bartlett.test(X~a,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,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.mis-aov(A~a,data=miscellany)summary(aov.mis)plot(miscellany$A~miscellany$a)说明:上述结果中,Df表示自由度;sumSq表示平方和;MeanSq表示均方和;Fvalue表示F检验统计量的值,即F比;Pr(F)表示检验的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.74,74.73,69.69),A=gl(4,4),B=gl(4,1,16))tanxin.aov-aov(X~A+B,data=tanxin)summary(tanxin.aov)#方差齐性检验bartlett.test(X~A,data=tanxin)#对因素Abartlett.test(X~B,data=tanxin)#对因素B结论:p值说明因素B对纤维弹性有影响,而没有充分理由说明因素A对纤维弹性有影响.结论:对因素A,p值(0.9243)远大于0.05,接受原假设,认为因素A的各水平下的数据是等方差的;对因素B,p值(0.8905)大于0.05,接受原假设,认为因素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