R軟件應用多元分析(II)

9.應用多元分析(II)安徽師范大學數(shù)學計算機科學學院丁新濤現(xiàn)在是1頁\一共有75頁\編輯于星期一9.1主成分分析9.1.1總體主成分主成分的定義與導出假定你是一個公司的財務經(jīng)理，掌握了公司的所有數(shù)據(jù)，比如固定資產(chǎn)、流動資金、每一筆借貸的數(shù)額和期限、各種稅費、工資支出、原料消耗、產(chǎn)值、利潤、折舊、職工人數(shù)、職工的分工和教育程度等等。如果讓你向上面介紹公司狀況，你能夠把這些指標和數(shù)字都原封不動地擺出去嗎？當然不能。你必須要把各個方面作出高度概括，用一兩個指標簡單明了地把情況說清楚。本章介紹兩種把變量維數(shù)降低以便于描述、理解和分析的方法：主成分分析（principalcomponentanalysis）和因子分析（factoranalysis）。實際上主成分分析可以說是因子分析的一個特例。現(xiàn)在是2頁\一共有75頁\編輯于星期一例子:成績數(shù)據(jù)100個學生的數(shù)學、物理、化學、語文、歷史、英語的成績?nèi)缦卤恚ú糠郑?。目前的問題是，能不能把這個數(shù)據(jù)的6個變量用一兩個綜合變量來表示呢？現(xiàn)在是3頁\一共有75頁\編輯于星期一主成分分析例中的的數(shù)據(jù)點是六維的；也就是說，每個觀測值是6維空間中的一個點。我們希望把6維空間用低維空間表示。以二維為例,如果這些數(shù)據(jù)形成一個橢圓形狀的點陣(比如二維正態(tài)分布)在短軸Z2方向上，數(shù)據(jù)變化很少(方差小);進一步,短軸如果退化成一點，則長軸Z1即可解釋這些點的變化;這樣,由二維到一維的降維就自然完成了.這相當于在平面上做一個坐標變換.Z1,Z2是X1,X2的特殊線性組合.X1X2Z1Z2現(xiàn)在是4頁\一共有75頁\編輯于星期一推廣到p維情況:對于多維變量的情況和二維類似，也有高維的橢球，只不過無法直觀地看見罷了。首先把高維橢球的主軸找出來，再用代表大多數(shù)數(shù)據(jù)信息的最長(方差大)的幾個軸作為新變量；這樣，主成分分析就基本完成了。注意，和二維情況類似，高維橢球的主軸也是互相垂直的。這些互相正交的新變量是原先變量的線性組合，叫做主成分(principalcomponent)。正如二維橢圓有兩個主軸，三維橢球有三個主軸一樣，有幾個變量，就有幾個主成分。選擇越少的主成分，降維就越好。什么是標準呢？那就是這些被選的主成分所代表的主軸的長度之和占了主軸長度總和的大部分。有些文獻建議，所選的主軸總長度占所有主軸長度之和的大約85%即可，其實，這只是一個大體的說法；具體選幾個，要看實際情況而定?，F(xiàn)在是5頁\一共有75頁\編輯于星期一定義:設X是p維隨機變量,μ=E(X),Σ=var(X)考慮線性變換(p維坐標主軸→p維橢球主軸):我們希望Z1的方差達到最大(Z1是橢球主軸中最長的一個主軸),則a1滿足:max=λmax，a1是Σ最大特征值(λ1)的特征向量.Z1=a1TX稱為第一主成分.現(xiàn)在是6頁\一共有75頁\編輯于星期一特征方程，λ是Σ的特征值a是對應的特征向量現(xiàn)在是7頁\一共有75頁\編輯于星期一類似地,希望Z2的方差達到最大,為保證橢球的主軸也是互相垂直的,a2與a1正交,即:cov(Z1,Z2)=a1T

Σa2=0,類似地,a2是Σ的第二大特征值(設為λ2),Z2=a2TX稱為第二主成分.一般地,即為所求線性變換矩陣A.Zi=aiTX即為第i個主成分.,Q是正交陣.現(xiàn)在是8頁\一共有75頁\編輯于星期一2.主成分的性質(zhì)主成分的均值和協(xié)方差陣.主成分的總方差.主成分分析是把p個原始變量X1,X2,…,Xp的總方差分解成p個不相關變量(cov(Z1,Z2)=0,)Z1,Z2,…,Zp的方差之和.總方差中第i主成分Zi的比例稱為主成分Zi的貢獻率.總方差中前m個主成分Zi的貢獻率之和稱Z1,Z2,…,

Zm的累積貢獻率.Xj與Zi之間的相關系數(shù).現(xiàn)在是9頁\一共有75頁\編輯于星期一2.主成分的性質(zhì)m個主成分對原始變量的貢獻率.總方差中前m個主成分Zi的貢獻率之和稱Z1,Z2,…,

Zm對X1,X2,…,Xp的累積貢獻率.Z1,Z2,…,

Zm對Xj的累積貢獻率:原始變量對主成分的影響.qji稱為第i個主成分在第j個原始變量Xj上的載荷,它度量了Xj對Zi的重要程度.現(xiàn)在是10頁\一共有75頁\編輯于星期一3.從相關矩陣出發(fā)求主成分的方差矩陣是X的相關矩陣R.p個主成分為:相關矩陣R的主成分性質(zhì):E(Z*)=0;var(Z*)=Λ*,其中Λ*=diag(λ1*,λ2*,…,λp*).變量Xj*與主成分Zi*之間的相關系數(shù)為:主成分Z1*,Z2*,…,Zm*對Xj*的貢獻率為:.設是相關矩陣R的p個特征值，是相應的單位特征向量。現(xiàn)在是11頁\一共有75頁\編輯于星期一9.1.2樣本主成分總體的參數(shù)在實際問題中，通常是未知的，則通過樣本來估計。樣本變量樣本方差矩陣:樣本相關矩陣R:相關記號:現(xiàn)在是12頁\一共有75頁\編輯于星期一1.從S出發(fā)求主成分S的特征值:λ1≧λ2≧…≧λp對應的特征向量(標準化):a1,a2,…,ap,ai稱為主軸向量(簡稱主軸);則第i個主成分zi=aiTX=aiT(x1,x2,…,xp)T令:z=(z1,z2,…,zp)T=(a1,a2,…,ap)Tx=QTx則樣本主成分:z(k)=QTX(k)每一個樣本的主成分樣本投影到新坐標系下的坐標X的第2個主成分實際問題中,經(jīng)常將樣本進行中心化.協(xié)方差矩陣不變現(xiàn)在是13頁\一共有75頁\編輯于星期一樣本主成分性質(zhì):樣本主成分的總方差等于原變量樣本的總方差Xj與Zi的樣本相關系數(shù).現(xiàn)在是14頁\一共有75頁\編輯于星期一2.從R出發(fā)求主成分樣本相關矩陣R的特征值:λ1*≧λ2*≧…≧λp*對應的特征向量(標準化):a1*,a2*,…,ap*,z(k)*=QTX(k)*性質(zhì):略現(xiàn)在是15頁\一共有75頁\編輯于星期一9.1.3相關R函數(shù)以及實例princomp{stats}Description:princompperformsaprincipalcomponentsanalysisonthegivennumericdatamatrixandreturnstheresultsasanobjectofclassprincomp.Usageprincomp(x,...)##S3methodforclass'formula':princomp(formula,data=NULL,subset,na.action,...)Formula:aformulawithnoresponsevariable,referringonlytonumericvariables.##DefaultS3method:princomp(x,cor=FALSE,scores=TRUE,covmat=NULL,subset=rep(TRUE,nrow(as.matrix(x))),...)X:anumericmatrixordataframewhichprovidesthedatafortheprincipalcomponentsanalysis.Cor:alogicalvalueindicatingwhetherthecalculationshouldusethecorrelationmatrixorthecovariancematrix.(Thecorrelationmatrixcanonlybeusediftherearenoconstantvariables.)現(xiàn)在是16頁\一共有75頁\編輯于星期一3-63.Loadings():Extractorprintloadingsinfactoranalysis(orprincipalcomponentsanalysis).Usage:loadings(x)x:anobjectofclass"factanal"or"princomp"ortheloadingscomponentofsuchanobject.predict():predictisagenericfunctionforpredictionsfromtheresultsofvariousmodelfittingfunctions.Thefunctioninvokesparticularmethodswhichdependontheclassofthefirstargument.Usage:predict(object,...)screeplot():screeplot.defaultplotsthevariancesagainstthenumberoftheprincipalcomponent.Thisisalsotheplotmethodforclasses"princomp"and"prcomp".Usage:screeplot(x,npcs=min(10,length(x$sdev)),type=c("barplot","lines"),main=deparse(substitute(x)),...)npcs:thenumberofcomponentstobeplotted.biplot():畫出數(shù)據(jù)關于主成分的散點圖和原坐標在主成分下的方向.Biplot(x,choices=1:2,scale=1,pc.biplot=FALSE,…)Choices:是選擇的主成分,默認是第1，2主成分.現(xiàn)在是17頁\一共有75頁\編輯于星期一7.實例序號X1X2X3X4114841727821393471763160497786414936677951594580866142316676715343768381504377799151427780101393168741114029647412161477884131584978831414033677715137316673161523573791714947827918145357077191604774872015644788521151427382221473873782315739688024147306575251574880882615136748027144366876281413067762913932687330148387078現(xiàn)在是18頁\一共有75頁\編輯于星期一R實現(xiàn):student=read.table('dataexample901.txt')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.00000000Loadings:Comp.1Comp.2Comp.3Comp.4x1-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#采用從R出發(fā)求主成分主成分貢獻率:Q由于前2個主軸的貢獻率達到96.4%,這個例子的主成分可以認為是comp.1(Z1*),Comp.2(Z2*)現(xiàn)在是19頁\一共有75頁\編輯于星期一預測:求Z*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.257399715282.160977780.137332330.355897390.093123683292.40434827-0.48613137-0.16154441-0.007914021300.502874680.14734317-0.20590831-0.122078819Comp1(Z1*)對應的系數(shù)符號都是負號,反映了學生身材的魁梧程度,身材高大的學生,他的4個部分值都比較大,則comp1的值就比較小,反之,身材’矮小的學生comp1的值比較大;Comp2(Z2*)是X1,X4與X2,X3的差,即:縱向高度與橫向圍度的差,所以,”細高”的同學comp2值越大,反之,該值越小說明樣本越”矮胖”.藍色樣本身材魁梧;粉色樣本身材瘦小;紅色樣本身材細高;綠色樣本身材矮胖;現(xiàn)在是20頁\一共有75頁\編輯于星期一結(jié)果圖示化:screeplot(student.pr,type='lines')碎石圖:縱坐標:(1.88178050.559806360.281795940.25711844)2通過碎石圖可以直觀的觀察出樣本在變換主軸方向的波動程度biplot(student.pr,choice=1:2)橫坐標是Comp1縱坐標是comp2矮胖細長魁梧瘦小現(xiàn)在是21頁\一共有75頁\編輯于星期一9.1.4主成分分析的應用1.主成分分類X1X2X3X4X5X6X7X8X9X10X11X12X13X14X15X16X11X20.791X30.360.311X40.960.740.381X50.890.580.310.91X60.790.580.30.780.791X70.760.550.350.750.740.731X80.260.190.580.250.250.180.241X90.210.070.280.20.180.180.29-0.041X100.260.160.330.220.230.230.250.49-0.341X110.070.210.380.08-0.0200.10.44-0.160.231X120.520.410.350.530.480.380.440.3-0.050.50.241X130.770.470.410.790.790.690.670.320.230.310.10.621X140.250.170.640.270.270.140.160.510.210.150.310.170.261X150.510.350.580.570.510.260.380.510.150.290.280.410.50.631X160.210.160.510.260.2300.120.380.180.140.310.180.240.50.651現(xiàn)在是22頁\一共有75頁\編輯于星期一R實現(xiàn):x=scan('dataexample902.txt')names=c('x1','x2','x3','x4','x5','x6','x7','x8','x9','x10','x11','x12','x13','x14','x15','x16')r=matrix(0,nrow=16,ncol=16,dimnames=list(names,names))for(iin1:16){for(jin1:i){r[i,j]=x[(i-1)*i/2+j]r[j,i]=r[i,j]}}pr=princomp(covmat=r)load=loadings(pr)plot(load[,1:2])text(load[,1],load[,2],adj=c(-0.4,0.3))#rij=x[1+2+…+(i-1)+j]長類圍類體型指標>summary(pr)Importanceofcomponents:Comp.1Comp.2Comp.3Standarddeviation2.6521.61679711.2775386ProportionofVariance0.43977980.16337700.1020066CumulativeProportion0.43977980.60315690.7051634前7個主成分占貢獻率87%現(xiàn)在是23頁\一共有75頁\編輯于星期一主成分聚類:Comp.1Comp.2x1-0.341770990.20040027x2-0.264991630.14320222x3-0.23415224-0.32862480x4-0.344232670.18112449x5-0.326117710.19965028x6-0.285913570.26980664x7-0.295261390.19214958x8-0.18927312-0.37026699x9-0.084792910.06747164x10-0.15429508-0.17424610x11-0.09835526-0.34784952x12-0.24254582-0.01766472x13-0.317158240.11191444x14-0.18011330-0.37135294x15-0.26635929-0.27122509x16-0.15833266-0.3628239316個樣本,每個樣本是2維向量;可以用樣本的距離對其實施聚類;write(load[,1:2],'dataexample902load12.txt')x=scan('dataexample902load12.txt')x=as.matrix(x);dim(x)=c(16,2)d=dist(scale(x));hc=hclust(d)plot(hc,hang=-1);re3=rect.hclust(hc,k=3)現(xiàn)在是24頁\一共有75頁\編輯于星期一用相關矩陣直接聚類d=as.dist(1-r);hc=hclust(d)plot(hc,hang=-1);re3=rect.hclust(hc,k=3)長類圍類體型指標長類主成分聚類圖現(xiàn)在是25頁\一共有75頁\編輯于星期一2.主成分回歸例9.3考慮進口總額Y與3個自變量:國內(nèi)總產(chǎn)值X1,存儲量X2和總消費量X3(單位為10億法郎)之間的關系.現(xiàn)收集了1949-1959年共11年數(shù)據(jù),如表,試對此數(shù)據(jù)做經(jīng)典回歸分析和主成分回歸分析.序號X1X2X3Y1149.34.2108.115.92161.24.1114.816.43171.53.1123.2194175.53.1126.919.15180.81.1132.118.86190.72.2137.720.47202.12.114622.78212.45.6154.126.59226.15162.328.110231.95.1164.327.6112390.7167.626.3現(xiàn)在是26頁\一共有75頁\編輯于星期一R實現(xiàn):conomy=read.table('dataexample903.txt')lm.sol=lm(Y~X1+X2+X3,data=conomy)#普通線性回歸summary(lm.sol)Coefficients:EstimatePr(>|t|)(Intercept)-10.127996.9e-05***X1-0.051400.488344X20.586950.000444***X30.286850.026277*Residualstandarderror:0.4889MultipleR-squared:0.9919AdjustedR-squared:0.9884F-statistic:285.6p-value:1.112e-07Y=-10.12799-0.0514X1+0.58695X2+0.28685X3lm.sol=lm(Y~X1+X2,data=conomy)Coefficients:EstimatePr(>|t|)(Intercept)-8.440140.000370***X10.145313.14e-08***X20.622480.001243**Residualstandarderror:0.6667MultipleR-squared:0.9828,AdjustedR-squared:0.9785F-statistic:228.3p-value:8.796e-08Y=-8.44014+0.14531X1+0.62248X2現(xiàn)在是27頁\一共有75頁\編輯于星期一殘差圖:Y=-8.44014+0.14531X1+0.62248X2H0:同方差,不是小概率事件接受HO:同方差;異方差檢驗:H0:殘差不相關,不是小概率事件接受HO:殘差不相關;不相關檢驗:現(xiàn)在是28頁\一共有75頁\編輯于星期一分析:λ3≈0X1,X2共線性Y=-8.44014+0.14531X1+0.62248X2后果很嚴重:回歸系數(shù)不穩(wěn)定總之:回歸結(jié)果仍不可信(多重共線性)Y=-10.12799-0.0514X1+0.58695X2+0.28685X3回歸系數(shù)不能通過顯著性t檢驗,所以該回歸結(jié)果不可信.Y=-8.44014+0.14531X1+0.62248X2回歸診斷：誤差項是否滿足獨立性、等方差性、正態(tài)性；選擇線性模型是否合適；是否存在異常樣本；回歸分析的結(jié)果是否對某些樣本的依賴過重，即回歸模型是否具備穩(wěn)定性；變量之間是否存在高度相關，即是否有多重共線性問題存在；(不通過)回歸殘差通過白噪聲檢驗怎么辦?現(xiàn)在是29頁\一共有75頁\編輯于星期一多重共線性的處理方法:刪去模型中次要的或可替代的解釋變量Y~x2不能通過系數(shù)顯著性t檢驗;Y~x1:y=-6.558101+0.146199X1,

R-squared 0.931761模型的檢驗都能通過.是一個”能用的”模型.實際上,一個更好的模型是y~x2+x3:y=-9.74+0.596X1+0.212X3R-squared 0.991277R2比y~x1的要高一些.是一個”比較好的”模型.現(xiàn)在是30頁\一共有75頁\編輯于星期一主成分法:conomy.pr=princomp(~X1+X2+X3,data=conomy,cor=T)summary(conomy.pr,loadings=TRUE)Importanceofcomponents:Comp.1Comp.2Comp.3Standarddeviation1.4139150.99907670.0518737839ProportionofVariance0.6663850.33271810.0008969632CumulativeProportion0.6663850.99910301.0000000000Loadings:Comp.1Comp.2Comp.3X10.7060.707X2-0.999

X30.707-0.707λ3=0.05187378392≈0X1,X2共線性#Z1,Z2的貢獻率已達99.9%小于0.05的數(shù)沒有顯示出來可以通過conomy.pr$loadings[][[2]]等查看0.706330.0356890.7069820.043501-0.999030.0069710.7065440.02583-0.7072現(xiàn)在是31頁\一共有75頁\編輯于星期一X1starx2starx3starZ1star-1.509720.545705-1.53319-2.12589-1.113050.485071-1.20848-1.61893-0.76971-0.12127-0.8014-1.11517-0.63637-0.12127-0.62209-0.8943-0.4597-1.33395-0.37008-0.64421-0.1297-0.66697-0.09869-0.190350.250307-0.727610.3035530.3596220.5936461.394580.6961010.9718021.050321.0307761.0934961.5593161.2436571.091411.1904221.7669951.480327-1.576481.3503491.931103pre=predict(conomy.pr)conomy$Z1=pre[,1]conomy$Z2=pre[,2]X*=-0.70720.025830.7065440.006971-0.999030.0435010.7069820.0356890.70633Q=(q1,q2,q3)=-2.2296493-1.6979452-1.1695976-0.9379462-0.6756511-0.19964230.37717461.01923441.63542431.85324012.0253583Z1*=X*q1lm.sol=lm(Y~Z1+Z2,data=conomy)summary(lm.sol)Coefficients:EstimatePr(>|t|)(Intercept)21.89091.21e-14***Z12.98926.02e-09***Z2-0.82880.00106**

Residualstandarderror:0.55MultipleR-squared:0.9883AdjustedR-squared:0.9853F-statistic:337.2p-value:

1.888e-08

#回歸結(jié)果比較理想.Y=21.89+2.99Z1*-0.83Z2*現(xiàn)在是32頁\一共有75頁\編輯于星期一將Z*表達成X:beta=coef(lm.sol)a=loadings(conomy.pr)x.bar=conomy.pr$centerx.sd=conomy.pr$scalecoef=(beta[2]*a[,1]+beta[3]*a[,2])/x.sdbeta0=beta[1]-sum(x.bar*coef)c(beta0,coef)(Intercept)Z1Z221.89090912.9891518-0.8287678

(Intercept)X1X2X3-9.130107820.072779810.609220120.10625939主成分法回歸方程為:Y=-9.13010782+0.07277981X1+0.60922012X2+0.10625939X3現(xiàn)在是33頁\一共有75頁\編輯于星期一9.2.2因子分析模型一般地，設X=(x1,x2,…,xp)T為可觀測的隨機變量，且有f=(f1,f2,…,fm)’為公共（共性）因子(commonfactor)，簡稱因子(factor);e=(e1,e2,…,ep)’為特殊因子（specificfactor）f和e均為不可直接觀測的隨機變量;μ=(μ1,μ2,…,μp)’為隨機變量x的總體均值;A=(aij)p*m為因子負荷（載荷）（factorloading）矩陣因子分析的模型為:假設通常先對x作標準化處理，使標準化得到的新變量均值為零，方差為１．現(xiàn)在是34頁\一共有75頁\編輯于星期一如果再滿足（４）fi與fj相互獨立（i≠j），則稱該因子模型為正交因子模型。正交因子模型具有如下特性：x的方差可表示為設(1)hi2是m個公共因子對第i個變量的貢獻，稱為第i個共同度（communality）或共性方差，公因子方差（commonvariance）(2)δi稱為特殊方差（specificvariance），是不能由公共因子解釋的部分(3)因子載荷（負荷）aij是隨機變量xi與公共因子fj的相關系數(shù)，系數(shù)aij是用來度量Xi可由f1,f2,…,fm線性組合表示的程度。現(xiàn)在是35頁\一共有75頁\編輯于星期一設稱gj2為公共因子fj對x的“貢獻”，是衡量公共因子fj重要性的一個指標。現(xiàn)在是36頁\一共有75頁\編輯于星期一三、因子分析的步驟輸入原始數(shù)據(jù)xn*p，計算樣本均值和方差，進行標準化計算（處理）；求樣本相關系數(shù)矩陣R=(rij)p*p；求相關系數(shù)矩陣的特征根λi

(λ1,λ2,…,λp>0)和相應的標準正交的特征向量li；確定公共因子數(shù)；計算公共因子的共性方差hi2;對載荷矩陣進行旋轉(zhuǎn)，以求能更好地解釋公共因子；對公共因子作出專業(yè)性的解釋?，F(xiàn)在是37頁\一共有75頁\編輯于星期一9.2.3參數(shù)估計(提取因子的方法)1.主成分法（principalcomponentfactor）

每一個公共因子的載荷系數(shù)之平方和等于對應的特征根，即該公共因子的方差。估計因子荷載矩陣A=(aij)p*n和特殊方差矩陣D特殊方差矩陣D：當總體的方差矩陣未知時，用樣本協(xié)方差矩陣S代替，sii是S對角線上元現(xiàn)在是38頁\一共有75頁\編輯于星期一factor.analy1factor.analy1<-function(S,m){

p<-nrow(S);diag_S<-diag(S);sum_rank<-sum(diag_S)

rowname<-paste("X",1:p,sep="")colname<-paste("Factor",1:m,sep="")A<-matrix(0,nrow=p,ncol=m,dimnames=list(rowname,colname))eig<-eigen(S)for(iin1:m)A[,i]<-sqrt(eig$values[i])*eig$vectors[,i]h<-diag(A%*%t(A))rowname<-c("SSloadings","ProportionVar","CumulativeVar")B<-matrix(0,nrow=3,ncol=m,dimnames=list(rowname,colname))for(iin1:m){B[1,i]<-sum(A[,i]^2)B[2,i]<-B[1,i]/sum_rankB[3,i]<-sum(B[1,1:i])/sum_rank}method<-c("PrincipalComponentMethod")list(method=method,loadings=A,var=cbind(common=h,spcific=diag_S-h),B=B)}為方便，設定中間變量A=(Factor1,…,FactormX1………0………

………0………Xp)p*m現(xiàn)在是39頁\一共有75頁\編輯于星期一R實例：例9.7對55個國家和地區(qū)的男子徑賽記錄作統(tǒng)計，每位運動員記錄8項指標：100m跑(X1)、200m跑(X2)、400m跑(X3)、800m跑(X4)、1500m跑(X5)、5000m跑(X6)、10000m跑(X7)、馬拉松(X8)。8項指標的相關矩陣R如表9.4所示。取m=2，用主成分法估計因子載荷和共性方差等指標?，F(xiàn)在是40頁\一共有75頁\編輯于星期一R實現(xiàn)：x=scan('dataexample907.txt')names=c('x1','x2','x3','x4','x5','x6','x7','x8')r=matrix(0,nrow=8,ncol=8,dimnames=list(names,names))for(iin1:8){for(jin1:i){r[i,j]=x[(i-1)*i/2+j]r[j,i]=r[i,j]}}source('factor.analy1.r')fa=factor.analy1(r,m=2)Fa$method[1]"PrincipalComponentMethod"$loadingsFactor1Factor2X1-0.8171562-0.53123478X2-0.8674064-0.43231147X3-0.9151503-0.23258703X4-0.9487239-0.01184340X5-0.95935870.13153096X6-0.93764640.29276177X7-0.94395700.28715151X8-0.87992530.41074922$varcommonspcificX10.94995460.05004538X20.93928700.06071301X30.89159680.10840319X40.90021730.09978270X50.93766960.06233044X60.96489030.03510972X70.97351080.02648920X80.94298340.05701655$BFactor1Factor2SSloadings6.62258840.8775213ProportionVar0.82782360.1096902CumulativeVar0.82782360.9375137公共因子fj對X1,…,Xp的總方差貢獻現(xiàn)在是41頁\一共有75頁\編輯于星期一9.2.4因子旋轉(zhuǎn)目的：使因子負荷兩極分化，要么接近于0，要么接近于1。常用的旋轉(zhuǎn)方法：（1）方差最大正交旋轉(zhuǎn)（varimaxorthogonalrotation）基本思想：使公共因子的相對負荷（lij/hi2）的方差之和最大，且保持原公共因子的正交性和公共方差總和不變。可使每個因子上的具有最大載荷的變量數(shù)最小，因此可以簡化對因子的解釋。（2）斜交旋轉(zhuǎn)（obliquerotation）因子斜交旋轉(zhuǎn)后，各因子負荷發(fā)生了較大變化，出現(xiàn)了兩極分化。各因子間不再相互獨立，而彼此相關。各因子對各變量的貢獻的總和也發(fā)生了改變。適用于大數(shù)據(jù)集的因子分析。現(xiàn)在是42頁\一共有75頁\編輯于星期一1.理論依據(jù)設因子模型：令：(是任意m階正交矩陣).即：對任一正交矩陣Γ，Z=ΓTF也是公因子向量。相應的AΓ是公因子Z的因子載荷矩陣。為此，在因子分析的實際計算中，當求得初始因子載荷矩陣A后，反復右乘正交矩陣Γ，使得AΓ具有更明顯的意義，這種變換載荷矩陣的方法，稱為因子軸的正交旋轉(zhuǎn)?，F(xiàn)在是43頁\一共有75頁\編輯于星期一2.因子載荷方差設因子模型A=(aij)p*m為因子載荷矩陣，為變量Xi的共同度。A的每一列(因子載荷向量)數(shù)值越分散，相應的因子載荷向量的方差越大。令：第j列p個數(shù)據(jù)的方差因子載荷矩陣A的方差為：Vj越大，A的第j列值越分散，Vj趨于1或0，稱公因子Fj具有簡單化結(jié)構(gòu)。V越大越好。現(xiàn)在是44頁\一共有75頁\編輯于星期一3.方差最大的正交旋轉(zhuǎn)所謂最大方差旋轉(zhuǎn)法就是選擇正交矩陣，使得矩陣所有個列元素平方的相對方差之和達到最大。當m=2時，設已求出的因子載荷矩陣為

現(xiàn)選取正交變換矩陣Γ進行因子旋轉(zhuǎn)，Γ可以表示為:這里θ是坐標平面上因子軸按順時針方向旋轉(zhuǎn)的角度，只要求出θ，也就求出了Γ.現(xiàn)在是45頁\一共有75頁\編輯于星期一現(xiàn)在是46頁\一共有75頁\編輯于星期一m>2:當m>2時，我們可以逐次對每兩個公共因子和進行上述旋轉(zhuǎn)。對公因子Fl和Fk進行旋轉(zhuǎn)，就是對A的第l和k兩列進行正交變換，使這兩列元素平方的相對方差之和達到最大，而其余各列不變，其正交變換矩陣為現(xiàn)在是47頁\一共有75頁\編輯于星期一其中θ是因子軸Fi和Fj的旋轉(zhuǎn)角度，矩陣中其余位置上的元素全為0。m個公共因子兩兩配對旋轉(zhuǎn)共需要進行m(m-1)/2次，稱其為完成了第一次旋轉(zhuǎn)，并記第一輪旋轉(zhuǎn)后的因子載荷矩陣為A(1)。然后再重新開始，進行第二輪的Cm2次配對旋轉(zhuǎn)，新的因子載荷矩陣記為A(2)。這樣可以得到一系列的因子載荷矩陣為現(xiàn)在是48頁\一共有75頁\編輯于星期一9.2.5因子分析的計算函數(shù)Factanal：factanal(x,factors,data=NULL,covmat=NULL,n.obs=NA,subset,na.action,start=NULL,scores=c("none","regression","Bartlett"),rotation="varimax",control=NULL,...)ArgumentsX:Aformulaoranumericmatrixoranobjectthatcanbecoercedtoanumericmatrix.Factors:Thenumberoffactorstobefitted.Data:Anoptionaldataframe(orsimilar:seemodel.frame),usedonlyifxisaformula.Bydefaultthevariablesaretakenfromenvironment(formula).Covmat:Acovariancematrix,oracovariancelistasreturnedbycov.wt.Ofcourse,correlationmatricesarecovariancematrices.Scores：Typeofscorestoproduce,ifany.Thedefaultisnone,“regression”givesThompson‘sscores,“Bartlett”givenBartlett’sweightedleast-squaresscores.Partialmatchingallowsthesenamestobeabbreviated.Rotation：character."none"orthenameofafunctiontobeusedtorotatethefactors:itwillbecalledwithfirstargumenttheloadingsmatrix,andshouldreturnalistwithcomponentloadingsgivingtherotatedloadings,orjusttherotatedloadings.現(xiàn)在是49頁\一共有75頁\編輯于星期一例9.11取m=2，用factanal()函數(shù)估計例9.7因子載荷和共性方差等指標，參數(shù)選擇方差最大。x=scan('dataexample907.txt')names=c('x1','x2','x3','x4','x5','x6','x7','x8')r=matrix(0,nrow=8,ncol=8,dimnames=list(names,names))for(iin1:8){for(jin1:i){r[i,j]=x[(i-1)*i/2+j]r[j,i]=r[i,j]}}fa=factanal(factors=2,covmat=r);faUniquenesses:

x1x2x3x4x5x6x7x8

0.0810.0750.1520.1350.0820.0330.0180.087現(xiàn)在是50頁\一共有75頁\編輯于星期一Loadings:

Factor1Factor2x10.2910.914

x20.3820.882

x30.5430.744

x40.6910.622

x50.7990.529

x60.9010.393

x70.9070.399

x80.9140.278

$loadingsFactor1Factor2X1-0.8171562-0.53123478X2-0.8674064-0.43231147X3-0.9151503-0.23258703X4-0.9487239-0.01184340X5-0.95935870.13153096X6-0.93764640.29276177X7-0.94395700.28715151X8-0.87992530.41074922Factor1Factor2SSloadings4.1133.224ProportionVar0.5140.403CumulativeVar0.5140.917Thedegreesoffreedomforthemodelis13andthefitwas0.3329現(xiàn)在是51頁\一共有75頁\編輯于星期一例9.12現(xiàn)有48名應聘者應聘某公司的某職位，公司為這些應聘者的15項指標打分，某指標與得分情況見例3.17，試用因子分析的方法對15項指標作因子分析，在因子分析中選取5個因子。rt=read.table("applicant.txt")factanal(~.,factors=5,data=rt)Loadings:

Factor1Factor2Factor3Factor4Factor5FL0.1270.7220.102-0.117APP0.4510.1340.2700.2060.258

AA0.1290.686

LA0.2220.2460.827

SC0.9170.167LC0.8510.1250.279-0.420

HON0.228-0.2200.777

SMS0.8800.2660.111EXP0.7730.171DRV0.7540.3930.1990.114AMB0.9090.1870.1120.165GSP0.7830.2950.3540.148-0.181POT0.7170.3620.4460.267KJ0.4180.3990.563-0.585

SUIT0.3510.7640.148f1:外露能力f2:經(jīng)驗f3:討人喜歡f4,f5:專業(yè)能力、外貌現(xiàn)在是52頁\一共有75頁\編輯于星期一實例：表7.5是研究消費者對購買牙膏偏好的調(diào)查數(shù)據(jù)。通過市場的攔截訪問，用7級量表詢問受訪者對以下陳述的認同程度（1表示非常不同意，7表示非常同意）。

V1：購買預防蛀牙的牙膏是重要的；

V2：我喜歡使牙齒亮澤的牙膏；

V3：牙膏應當保護牙齦；

V4：我喜歡使口氣清新的牙膏；

V5：預防壞牙不是牙膏提供的一項重要利益；

V6：購買牙膏時最重要的考慮是富有魅力的牙齒?，F(xiàn)在是53頁\一共有75頁\編輯于星期一表7.5牙膏屬性評分得分表現(xiàn)在是54頁\一共有75頁\編輯于星期一將表7.5中的數(shù)據(jù)通過SPSS進行因子分析1.特征根和累計貢獻率提取兩個因子累計方差貢獻率就達到82%，第三個特征根相比下降較快，因此我們選取兩個公共因子?，F(xiàn)在是55頁\一共有75頁\編輯于星期一2.因子的含義從因子載荷陣可以看出：因子1與V1（預防蛀牙），V3（保護牙齦），V5（預防壞牙）相關性強，其中V5的載荷是負數(shù)，是由于這個陳述是反向詢問的；因子2與V2（牙齒亮澤），V4（口氣清新），V6（富有魅力）的相關系數(shù)相對較高。因此，我們命名因子1為“護牙因子”，是人們對牙齒的保健態(tài)度；因子2是“美牙因子”，說明人們“‘通過牙膏美化牙齒’影響社交活動”的重視。從這兩方面分析，對牙膏生產(chǎn)企業(yè)開發(fā)新產(chǎn)品都富有啟發(fā)意義。旋轉(zhuǎn)后因子載荷矩陣現(xiàn)在是56頁\一共有75頁\編輯于星期一R實現(xiàn)：rt=read.table("dataexample916.txt")factanal(rt,factors=2)Loadings:Factor1Factor2v10.968v20.749v30.898-0.140v40.784v5-0.887v60.830Factor1Factor2SSloadings2.5421.892ProportionVar0.4240.315CumulativeVar0.4240.739現(xiàn)在是57頁\一共有75頁\編輯于星期一典型相關性分析典型相關分析（canonicalcorrelationanalysis）是研究兩組變量之間相關關系的一種統(tǒng)計分析方法，它能夠有效地揭示兩組變量之間的相互線性依賴關系。我們知道，在一元統(tǒng)計分析中，用相關系數(shù)來衡量兩個隨機變量之間的線性相關關系；用復相關系數(shù)研究一個隨機變量和多個隨機變量的線性相關關系。然而，這些統(tǒng)計方法在研究兩組變量之間的相關關系時卻無能為力。比如要研究生理指標與訓練指標的關系，居民生活環(huán)境與健康狀況的關系，人口統(tǒng)計變量（戶主年齡、家庭年收入、戶主受教育程度）與消費變量（每年去餐館就餐的頻率、每年出外看電影的頻率）之間是否具有相關關系？閱讀能力變量（閱讀速度、閱讀才能）與數(shù)學運算能力變量（數(shù)學運算速度、數(shù)學運算才能）是否相關？這些多變量間的相關性如何分析？現(xiàn)在是58頁\一共有75頁\編輯于星期一9.3.1總體典型相關1936年霍特林（Hotelling）最早就“大學表現(xiàn)”和“入學前成績”的關系、政府政策變量與經(jīng)濟目標變量的關系等問題進行了研究，提出了典型相關分析技術。之后，Cooley和Hohnes(1971)，Tatsuoka(1971)及Mardia，Kent和Bibby(1979)等人對典型相關分析的應用進行了討論，Kshirsagar(1972)則從理論上給出了最好的分析。典型相關分析的目的是識別并量化兩組變量之間的聯(lián)系，將兩組變量相關關系的分析，轉(zhuǎn)化為一組變量的線性組合與另一組變量線性組合之間的相關關系分析。其基本思想和主成分分析非常相似。首先在每組變量中找出變量的線性組合，使得兩組的線性組合之間具有最大的相關系數(shù)。然后選取和最初挑選的這對線性組合不相關的線性組合，使其配對，并選取相關系數(shù)最大的一對，如此繼續(xù)下去，直到兩組變量之間的相關性被提取完畢為此。被選出的線性組合配對稱為典型變量，它們的相關系數(shù)稱為典型相關系數(shù)。典型相關系數(shù)度量了這兩組變量之間聯(lián)系的強度?，F(xiàn)在是59頁\一共有75頁\編輯于星期一引入：一般情況，設是兩個相互關聯(lián)的隨機向量，分別在兩組變量中選取有代表性的綜合變量U、V，使得每一個綜合變量是原變量的線性組合，即典型變量

我們希望尋找使相關系數(shù)達到最大的向量a與b，由于隨機向量乘以常數(shù)時并不改變它們的相關系數(shù)，所以，為防止結(jié)果的重復出現(xiàn)，令

現(xiàn)在是60頁\一共有75頁\編輯于星期一根據(jù)條件極值的求法引入Lagrange乘數(shù)，將問題轉(zhuǎn)化為求的極大值，其中λ，ν是Lagrange乘數(shù)。根據(jù)求極值的必要條件得

A和B具有相同的特征根，a，b則是其相應的特征向量.，最大特征根λ2對應的特征向量就是所求的典型變量的系數(shù)向量.

2.典型變量和典型相關系數(shù)的計算現(xiàn)在是61頁\一共有75頁\編輯于星期一計算過程：求A的最大特征值和相應的特征向量，令：為第1對典型相關系數(shù)，為第1對典型變量。求A的第k個最大特征值和相應的特征向量，令：為第1對典型相關系數(shù)，為第k對典型變量。現(xiàn)在是62頁\一共有75頁\編輯于星期一樣本典型相關在實際分析應用中，總體的協(xié)差陣通常是未知的，往往需要從研究的總體中隨機抽取一個樣本，根據(jù)樣本估計出總體的協(xié)差陣，并在此基礎上進行典型相關分析。設服從正態(tài)分布從該總體中抽取樣本容量為n的樣本，得到下列數(shù)據(jù)矩陣：樣本均值向量,其中，樣本協(xié)差陣現(xiàn)在是63頁\一共有75頁\編輯于星期一

計算過程：求A的第k個最大特征值和相應的特征向量，令：為第1對典型相關系數(shù)，為第k對典型變量?，F(xiàn)在是64頁\一共有75頁\編輯于星期一9.3.3典型相關分析的計算R軟件函數(shù)：cancorcancor(x,y,xcenter=TRUE,ycenter=TRUE)Argumentsxnumericmatrix(n*p1),containingthexcoordinates.ynumericmatrix(n*p2),containingtheycoordinates.xcenterlogicalornumericvectoroflengthp1,describinganycenteringtobedoneonthexvaluesbeforetheanalysis.IfTRUE(default),subtractthecolumnmeans.IfFALSE,donotadjustthecolumns.Otherwise,avectorofvaluestobesubtractedfromthecolumns.ycenteranalogoustoxcenter,butfortheyvalues.現(xiàn)在是65頁\一共有75頁\編輯于星期一表9.6生理指標和訓練指標序號X1X2X3Y1Y2Y311913650516