R作主成分分析主成分分类和主成分回归

9.1.3相关的R函数以及实例卜•面介绍与主成分分析有关的函数.i.princomp函数作主成分分析最主要的函数是pnncompO函数，其使用格式为princomp(formula,data=NULL,subset,na.action,••・)其中formula是没有响应变飒的公式(类似冋归分析、方差分析，但无响应变僦)・data是数据框(类似于冋归分析、方差分析).或者princomp(x,cor=FALSE,scores=TRUE,co\nnat=NULL,subset=rep(TRUE,nrow(as.niaXrix(x))),•••)其中X是用于主成分分析的数据，以数值矩阵或数据框的形式给出.cor是逻辑变ht,当cor=TRUE表示用样本的相关矩阵R作主成分分析，当cor=FALSE｛缺省值)表示用样本的协方差阵S作主成分分析.covmat是协方差阵，如果数据不用x提供，可由协方差阵提供.其他参数的意义见在线帮助.prcompO函数的意义与使用方法与princomp0函数相同.2.summary函数summary()与冋归分析中的用法相同，其目的是提取主成分的信息，其作用格式为summary(object,loadings=FALSE,cutoff=0.1,...)其中object是由princomp()得至lj的对象.loadings是逻辑变ht,当loadings=TRUE表示显示loadings的内容具体含义在|、•面的loadings()函数］.当loadings=FALSE则不显示.3・loadings函数loadings()函数是显示主成分分析或因子分析中loadings载荷，见因子分析)的内容.在主成分分析中，该内容实际上是主成分对应的各列，即前面分析的正交矩阵Q.析的正交矩阵Q.4子分析中，其内容就是载荷因子矩阵loadings()函数的使用格式为loadings(x)其中x是由函数princomp()或factanal()(JSL因子分析得到的对象.predict函数P"di“()函数是预测主成分的ff[(类似于冋归分析中的使用方法).苴使用格式为predict(object,newdata,•••)其中object是由princompO得到的对象.newdata是曲预测值构成的数据框，当newdata缺省时，预测已有数据的主成分值.screeplot函数screeplorO函数是画出主成分的碎石图，其使用格式为screeplot(x,npcs=minClO,length(x$sdev)),type=c(l,barplotM,MlinesH),main=deparse(substitute(x))f…)比中X是由princompO得到的对象.npcs是画出的主成分的个数.type是描述画出的碎石图的类型，-barplof是直方图类型，-lines-是直线图类型.biplot函数biplot()是画出数据关于主成分的散点图和原坐标在主成分卜•的方向，其使用格式为biplot(x,choices=1:2,scale=1,pc.biplot=FALSE,...)其中x是由princompO得至lj的对象.choices是选择的主成分，缺省值是第1、第2主成分.pc.biplot是逻辑变肚:缺省值为FALSE),当pc.biplot=TRUE.用Gabriel(1971)提出的画图方法.实例I、•面用一个例子说明前面介绍的网数的使用方法.例9.1(中学生身体四项指标的主成分分析)在某中学随机抽取某年级％名学生，測量其身高(XJ、体重(XJ、胸围(X3)和坐高(X4),数据如表9.1所示.试对这30名中学生身体四项指标数据做主成分分析.解:用数据框的形式输入数据.用princompO作主成分分析，由前面的分析，选择相关矩阵作主成分分析更合理，因此这里选择的参数是cor=TRUE.最后用summary()列出主成分分析的值，这里选择loadings=TRUE.U|、•是相应的

表9・1:30名中学丰身体四呗指标数据序号&X2%3X4序号Xi X2X3X4114841727816152357379213934717617149478279316049■■((86181453570■■414936677919160477487■015945808620156447885614231667621151427382••15343768322147387378815043Ml■((7923157396880915142■■1(802414730657510139316871251574880881114029647426151367480121614778842714436687613158■19788328141306776141403367■■«291393268731513731667330148387078>student<-data•:frame>student<-data•:frame(Xl=c(148,139,160,149,140,161,158,140,151,147,157,147,X2=c(41,34,49,36,45,29,47,49,33,31,42,38,39,30,48,X3=c(72,71,77,67,80,64,78,78,67,66,73,73,68,65,80,X4=c(78,76,86,79,86,74,84,83,77,73,82,78,80,75,88,####用数据框形式输入数据)159,142,1153,150,151,139,137,152,1-49,145,160,156,157,151,1<44,141,139,148)31,43,43,42,31,35,47,35,47,44,36,36,30,32,38),66,76,77,77,68,73,82,70,74,78,74,68,67,68,70),76,83,79,80,74,79,79,77,87,85,80,76,76,73,78)####作主成分分析，并显示分析结果student.pr<・princomp(student,cor=TRUE)summary(student.pr,loadings=TRUE)Importanceofcomponents:Comp・1Comp•2Comp・3Comp•4Standarddeviation1.88178050.559806360.281795940.25711844ProportionofVariance0.88527450.078345790.019852240.01652747CumulativeProportion0.88527450.963620290.983472531.00000000Standarddeviation1.88178050.559806360.281795940.25711844ProportionofVariance0.88527450.078345790.019852240.01652747CumulativeProportion0.88527450.963620290.983472531.00000000Loadings:Comp.1Comp・Comp.1Comp・2Comp・3Comp.4XI-0.4970.543-0.4500.506X2-0.515-0.210-0・462-0.691X3-0.481-0.7250.1750.461X4-0.XI-0.4970.543-0.4500.506X2-0.515-0.210-0・462-0.691X3-0.481-0.7250.1750.461X4-0.5070.3680.744-0.232在上述程序中，语句studenr.pr<-princomp(srudenr,cor=TRUE)可以student.pr<-princomp(*X1*X2*X3*X4,data=student,cor二TRUE).两者是等价的.summary()函数列出了主成分分析的重要倍息，Standarddeviation行表示的是主成分的标准差，即主成分的方差的开方，也就是相就的持征值入1.A2.入3.A4的开方. ProportionofVariance行表示的是方差的贡献率.emulativeProportion行表示的是方差的累积贡献率.由TtEsummarypfq数的参数中选取了loadings=TRUE.因此列出了loadings:载荷;的内容，它实际上是主成分对应于原始变Ift 的系数，即前面介绍的矩阵Q.因此，得到Z：=-0.497Xf-0.515X；-0.481X"-0.507X；,=0.543X：-0.210X；-0.725X；+0.368X；.由干前两个主成分的累积贡献率已达到96%另外二个主成分可以舍去，达到降维的目的.第1主成分对应系数的符号都相同，其值在0.5左右，它反映了中学生身材魁梧程度：身体高人的学生，他的4个部分的尺寸都比较大，因此第1主成分的值就较小(因为系数均为负值)；而身材矮小的学生，他的1部分的尺寸都比校小因此，第1主成分绝对值就较人.我们称第1主成分为人小因子.第2主成

分是高度与围度的差，第2主成分值人的学生表明该学生叱［0高••而第2主成分值越小的学生表明该学生••矮胖S因此，祢笫2主成分为体形因子.我们看一下各样本的主成分的值（用predict0函数）.####作预测>predict(student.pr)Comp.1Comp.2Comp.3Comp.410.06990950-0.23813701-0.35509248-0.26612013921.59526340-0.718473990.32813232-0.1180566463-2.847931510.38956679-0.09731731-0.27948248740.759969880.80604335-0.04945722-0.1629492985-2.739667770.017180870.360126150.35865304462.105831680.322843930.18600422-0.0364560847-1.42105591-0.060531650.21093321-0.0442230928-0.82583977-0.78102576-0.275577980.0572885729-0.93464402-0.58469242-0.088141360.181037746102.36463820-0.365321990.088404760.045520127112.837419160.348758410.03310423-0.03114693012-2.608512240.21278728-0.333980370.21015757413-2.44253342-0.16769496-0.46918095-0.162987830141.866306690.050213840.37720280-0.358821916152.81347421-0.31790107-0.03291329-0.222035112160.063929830.207184480.043343400.70353362417-1.55561022-1.70439674-0.331264060.007551879181.07392251-0.067634180.022836480.04860668019-2.521742120.972743010.12164633-0.39066799120-2.140723770.022178810.374109720.12954896021-0.796244220.163078870.12781270-0.294140762220.28708321-0.35744666-0.039621160.08099198923-0.251510751.25555188-0.556173250.109068939242.057060320.78894494-0.265521090.38808864325-3.08596855-0.057753180.62110421-0.21893961226-0.163675550.043179320.244818500.560248997271.372650530.02220972・0・23378320・0.257399715

282.160977780.137332330.355897390.093123683292.40434827-0.48613137-0.16154441-0.007914021300.502874680.14734317-0.20590831-0.122078819从第1主成分来看，较小的几个值是25号样本.3号样本和5号样本,因此说明这几个学生身材魁怖而11号样本、15号样本和29号样本的的值较大，说明这几个学生身材瘦小.从第2主成分来看，较人的几个值是23号样本、19号样本和4号样本，因此说明这几个学生属于••细高”型；而17号样本、8号样本和2号样本的的值较小，说明这几个学生身材属于•矮聆型.画出主成分的的碎石图.>screeplot(student.pr>type=,,lines,1)参数选择的直线型，其图形如图9.1所示.student.pr图9.1:30名中学生身体指标数据主成分的碎石图还可以画出关于第1主成分和第2主成分样本的散点图，其图形如图9.2所示.从该散点图可以很容易看出：哪些学生属于髙人魁桥型，如25号学生，哪些学生属于身材瘦小型，如11号或15号；哪些学生属于“细討型，如23号哪些学生属于•矮聆型，如17号还有哪些学生属于正常体形，如26号，等表92表9216项身体指标数据的相关矩阵X2X3 X5X6XzX8X9X10XnX12X13XuX15%16-4 -4 -2 0 2 4Comp.1图9.2:30名中学生身体指标数据关于第1主成分和第2主成分的散点图9.1.4主成分分析的应用这一小节讲两个问题作为主成分分析的应用，f是变吊分类问题；另一个是主成分回归问题.1.主成分分类例9.2对胆8个成年男子的身材进疔测量，每人各測得"项指标：身高（XJ、坐高Q®、胸團（X3丿、头高（XJ、裨长（X5）、下档八手长（X*、领围（XJ、前胸（X』、后背（X】。丿、肩厚（X】J、肩宽（X\2）、袖长久3八肋ffl（Xj、腰围丿和腿肚（Xw丿.16项指标的相关矩阵R如表9.2所示（由于相关矩阵是对称的，只给出下三角部分丿.迖从相关矩阵R出发进疔主成分分析，对16项/脱进疔分类.

X1X2X3X4X5X6X?A8X9X1()X1lHxl3XUX!5X161.000.790.360.960.S90.790.760.260.211.000.311.000.740.380.580.310.580.300.550.350.190.580.070.281.000.900.780.750.250.200.260.160.330.221.000.790.740.250.180.230.070.210.380.08-0.020.520.410.350.770.470.410.250.170.640.510.350.580.210.160.511.000.73LOO0.180.241.000.180.29-0.040.230.250.000.101.000.490.440.530.790.270.570.260.480.790.270.510.230.380.440.690.670.140.160.260.380.000,120.30X1X2X3X4X5X6X?A8X9X1()X1lHxl3XUX!5X161.000.790.360.960.S90.790.760.260.211.000.311.000.740.380.580.310.580.300.550.350.190.580.070.281.000.900.780.750.250.200.260.160.330.221.000.790.740.250.180.230.070.210.380.08-0.020.520.410.350.770.470.410.250.170.640.510.350.580.210.160.511.000.73LOO0.180.241.000.180.29-0.040.230.250.000.101.000.490.440.530.790.270.570.260.480.790.270.510.230.380.440.690.670.140.160.260.380.000,120.300.320.510.510.38-0.341.00-0.160.231.00-0.050.500.241.000.230.210.150.180.310.100.621.000.150.310.170.261.000.290.280.410.500.631.000.140.310.180.240.500.651.00解：首先输入相关矩阵，再用princompO对相关矩阵作主成分分析，最后画岀各变笊在第一、第二主成分下的散点图（程序名：exan0902.R）####输入数据，按下三角输入，构成向量x<-c(l.00,0.79,1.00,0.36,0.31,1.00,0.96,0.74,0.38,1.00,0.89,0.58,0.31,0.90,1.00,0.79,0.58,0.30,0.78,0.79,1.00,0.76,0.55,0.35,0.75,0.74,0.73,1.30,0.26,0.19,0.58,0.25,0.25,0.18,0.24,1.00,0.21,0.07,0.28,0.20,0.18,0.18,0.29,--0.04,1.00,0.26,0.16,0.33,0.22,0.23,0.23,0.25,0.49,-•0.34,1.00,0.07,0.21,0.38,0.08,-■0.02,0.00,0.10,0.44,-0.16,0.23,1.00,0.52,0.41,0.35,0.53,0.48,0.38,0.44,0.30,-0.05,0.50,0.24,1.00,0.77,0.47,0.41,0.79,0.79,0.69,0.67,0.32,0.23,0.31,0.10,0.62,1.00,0.25,0.17,0.64,0.2乙0.27,0.14,0.16,0.51,0.21,0.15,0.31,0.17,0.26,1.00,0.51,0.35,0.58,0.57,0.51,0.26,0.38,0.51,0.15,0.29,0.28,0.41,0.50,0.63,1.00,0.21,0.16,0.51,0.26,0.23,0.00,0.12,0.38,0.18,0.14,0.31,0.18,0.24,0.50,0.65,1.00)####输入变量名称nani€sJc(”Xl”，"X2","X3",,,X4",nX5","X6","X7","X8","X9","X10","Xll","X12","X13","X14","X15","X16")胖料将矩阵生成相关矩阵R<-Eatrix(0,nrow=16,ncol=16,dimnames=list(names,names))for(iin1:16){for(jin1:i){R[i,j]<-x[(i-l)*i/2+j];R[j,i]<-R[i,j]}}料料作主成分分析pr<・princomp(covmat=R);load<-loadmgs(pr)###*»画散点图plotCload[fl:2]);text(load[,1]9load[.2],adj=c(-0.4f0.3))得到的图形由图9.3所示・图9.3中左上角的点看成一类，它们是欲类，即身高(Xi)、坐髙(X2)、头盒(%4)、裤长(疋)、下档(%6)、手长(%7)、袖长(X13).右下角的点看成一类，它们是W类，即身胸围(Xd、领围(XJ、肩厚(Xu)、肋围(X14)、腰围(%15)、腿B上(X®.中间的点看成一类，为体形持征指标,即前胸(XJ、后背(Xl0)、肩宽(%12).06CM©"00o7o2Lo13oo9ooo12LCMo10o15o3o8o14o16o11-0.35 -0.30 -0.25 -0.20 -0.15 -0.10Comp.1图9.3:16个变球在第一、第二主成分卜.的散点图2.主成分回归在冋归分析一章中，曾经讲过，当自变ht出现多重共线性时，经典冋归方法作冋归系数的最小二乘估计，V效果会较差，而采用主成分冋归能够克服直接冋归的不足.I、•面用一个例子来说明如何作主成分冋归，并J1是如何克服经典冋归的不足.例9.3（法国经济分析数据）考虑进口总额y与三个自变量：国内总产值Xi,存储量x2,总消彷量X3/单位为M亿法郎丿之间的关系•现收集了1949年至"59年共门年有數据，如表9.3所示.试对此数据作经典回归分析和主成分回归分析.解：输入数据（采用数据框形式人再用一般线性冋归方法作冋归分析（程序名：exam0903.R\####用数据框的形式输入数据>conomy<・dara・frame（xl-c（149.3f161.2,171.5,175.5,180.8,190.7,

202.1,212.4,226.1,231.9,239.0),x2=c(4.2t4.1,3.1,3.1t1.1,2.2,2.1f5.6,5.0,5.1,0.7),x3=c(108.1,114.8,123.2,126.9,132.1,137.7,146.0,154.1,162.3,164.3,167.6),y=c(15.9,16.4,19.0,19.1,18.8,20.4,22.7,26.5,28.1,27.6,26.3)衣"2；序号XiX2x31149.34.210&115.92161.24.1114.816.43171.53.1123.219.041r0.53.1126.919.15180.81.1132.118.86190.72.2137.720.4■202.12.1146.022.78212.45.6154.126.59226.15.0162.328.110231.95.1164.327.611239.00.7167.626.3####作线性回归Izi.sol<-lm(y*xl+x2+x3,data=cononiy)summary(Im.sol)Call:Im(formula=y"xl*x2*x3,data=conomy)Residuals:MmIQMedian3QMmIQMedian3QMax-0.52367-0.38953 0.05424 0•22644 0.78313Coefficients:EstimateStd.ErrortvaluePrOItl)(Intercept)-10.127991.21216-8.3556.9e-05X1-0.051400.07028-0.7310.488344x20.586950.094626.2030.000444♦♦♦x30.286850.102212.8070.026277*Signif.codes:0朴0.001"杆0.01 0.05 0.1J11Residualstandarderror:0.4889on7degreesoffreedomMultipleR-Squared:0.9919,AdjustedR-squared:0.9884F-statistic:285.6on3and7DF,p-value:1.112e-07从计算结果可以看出，按三个变秋得到冋归方程Y=-10.12799一0.05140X1-0.58695X2*0.28685X3. (9.14)仔细分析方程(9-14),发现它并不合理.冋到问题本身，丫是进口肚，X】是国内总产值，而对应系数的符号确为负，也就是说，国内的总产值越高，其进nht确越少，这与实际情况是不相符的.问其原因，三个变ht存在着多重共线性(后面我们将会看到最才喑征值接近于0).为克服多重共线性的影响，对变伙作作主成分冋归.先作主成分分析.#轴#作主成分分析conomy・pr<・princonip(~xl.・x2*x3,data=conomy>cor=T)summary(conomy・pr,loadmgs=TRUE)Importanceofcomponents:Comp•1Comp.2 Comp・3Standarddeviation1.4139150.99907670.0518737839ProportionofVariance0.6663850.33271810.0008969632CumulativeProportion0.6663850.99910301.0000000000

Loadings:Comp・1Comp・2Comp.3xl0.7060.707x2-0.999x30.707-0・707前两个主成分已达到99%的贡献率.第1主成分是关于国内总产值和总消费，因此称第1主成分为产销因子.第2主成分只与存储尿有关，称为存储因子.注意，入3=0.05187378392=0.(X)2690889%0,所以变ht存在着多重共线性.卜•面作主成分冋归.首先计算样本的主成分的预测值，并将第1主成分的预测值和第2主成分的预测值存放在数据框conomy中，然后再对主成分作冋归分析.其命令格式如F####预测测样本主成分，并作主成分分析pre<-predict(conomy.pr)conomySzl<-pre[,1];conomy$z2<・pre[,2]Im.sol<-lm(y*zl*z2>data=conomy)siimmary(lm.sol)Call:Im(formula=yzl+z2,data=conomy)Residuals:Min IQMedian 3QMax■0・89838・0・26050 0.08435 0.35677 0.