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｛缺省值)表示用樣本的協(xié)方差陣S作主成分分析.covmat是協(xié)方差陣，如果數據不用x提供，可由協(xié)方差陣提供.其他參數的意義見在線幫助.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()函數是顯示主成分分析或因子分析中l(wèi)oadings載荷，見因子分析)的內容.在主成分分析中，該內容實際上是主成分對應的各列，即前面分析的正交矩陣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主成分分析的應用這一小節(jié)講兩個問題作為主成分分析的應用，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出發(fā)進疔主成分分析，對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),發(fā)現它并不合理.冋到問題本身，丫是進口肚，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.