多元統(tǒng)計分析-基于R（第3版） 課件 費宇 第9、10章-主成分分析、因子分析

第9章主成分分析2024/9/52

主成分分析（principalcomponentanalysis，PCA）是將具有相關關系的多個變量有效地轉化為少數(shù)幾個綜合變量來處理，從而簡化相關統(tǒng)計分析的一種多元統(tǒng)計方法．主成分分析也稱主分量分析，是由Pearson于1901年首先提出，到1933年由Hotelling加以完善后發(fā)展起來的．目前，在涉及到高維數(shù)據(jù)分析處理的諸多領域主成分分析都有廣泛的應用．本章介紹主成分的基本思想、總體主成分和樣本主成分的定義、性質、主成分的計算和解釋、主成分回歸等．2024/9/53

基本思想：用個數(shù)較少，但是保留了原始變量的大部分信息幾個不相關的綜合變量(即主成分)來代替原來較多的變量，從而可以簡化數(shù)據(jù)，對原來復雜的數(shù)據(jù)關系進行簡明有效的統(tǒng)計分析．

主成分分析的本質就是“有效降維”，既要減少變量個數(shù)，又不能損失太多信息．換句話說，就是“降噪”、“擠水分”或者說“冗余消除”，將高維數(shù)據(jù)有效地轉化為低維數(shù)據(jù)來處理．揭示變量之間的內在聯(lián)系，進而分析解決實際問題.

第9章

主要內容9.1

主成分分析的基本思想9.2

總體主成分

主成分的含義、計算、主要性質、主成分個數(shù)的確定、變量的標準化．9.3

樣本主成分

樣本主成分的性質和計算；主成分分析的步驟和相關R函數(shù)．9.4案例：主成分回歸分析2024/9/54

9.1

主成分分析的基本思想2024/9/55

實際統(tǒng)計分析中，經(jīng)常處理多變量、大維數(shù)數(shù)據(jù)分析問題，分析過程較復雜，難度較大.全部變量中可能存在信息的重疊．為去除這些信息重疊，希望用個數(shù)較少，但是保留了原始變量的大部分信息幾個不相關的綜合變量(即主成分)來代替原來較多的變量，達到降維的目的，從而進行簡明有效的統(tǒng)計分析.主成分分析中的信息，就是指變量的變異性，常用標準差或方差來表示它.9.2總體主成分

9.2.1

主成分的含義

以二維正態(tài)分布樣本點來直觀說明：如圖，設有n個樣品點大致分布在平面上一個橢圓內：2024/9/56

圖9－1二維情形主成分示意圖二維情形主成分的解釋

樣本點之間的差異是由

x1和

x2的變化引起的，兩者變動的相差不大，但如果用新坐標

y1和

y2來代替，易見，這些樣本點的差異主要體現(xiàn)在

y1軸上，n個點在

y1軸方向上的方差達到最大，即在此方向上包含了有關n個樣品的最多的信息．

將這些點投影到

y1軸方向能使信息的損失最小，如果

y1軸方向的差異占了全部樣本點差異的絕大部分，那么將

y2忽略是合理的，這樣就把兩個變量簡化為一個，顯然這里的

y1軸代表了數(shù)據(jù)變化最大的方向，稱之為第一主成分．y2稱為第二主成分，并要求已經(jīng)包含在

y1中的信息不出現(xiàn)在

y2中,即2024/9/57

二維情形主成分的解釋(續(xù))

注意兩個主成分

y1和

y2都是

x1和

x2的線性組合：其中P為旋轉變換矩陣，它是正交矩陣．

極端情形1：橢圓變成圓，第一主成分

y1只體現(xiàn)了原始二維樣品點約一半的信息，若此時將

y2忽略，則將損失約50％的信息，這顯然是不可取的．

極端情形2：橢圓扁平到了極限，變成

y1軸上一條線段，第一主成分

y1幾乎包含有二維樣品點的全部信息，僅用

y1代替原始數(shù)據(jù)幾乎不會有任何的信息損失，此時降維效果是非常理想的.2024/9/58

利用R程序來模擬這一過程（先下載安裝mvtnorm）2024/9/59

>library(mvtnorm)

#先加載多元正態(tài)及t分布程序包mvtnorm>set.seed(8)

#設置隨機數(shù)種子>sigma<-matrix(c(1,0.9,0.9,1),ncol=2)#設定協(xié)方差矩陣,相關系數(shù)為0.9>mnorm<-rmvnorm(n=200,mean=c(0,0),sigma=sigma)>plot(mnorm)#產(chǎn)生200個二維正態(tài)分布隨機數(shù)并畫散點圖（見圖9-2）>abline(a=0,b=1);abline(a=0,b=-1)#畫坐標軸旋轉45度后的二條直線>eig<-eigen(cor(mnorm));eig#求特征值和特征向量>vector1<-eig$vectors[,1];vector2<-eig$vectors[,2]>y1<-scale(mnorm)%*%vector1;y2<-scale(mnorm)%*%vector2#將數(shù)據(jù)標準化>plot(y1,y2,ylim=c(-2,2));abline(h=0,v=0)

#見圖9-3>cbind(var(y1),var(y2),cor(y1,y2))

[,1][,2][,3][1,]1.885428

0.11457184.418324e-16#可見y1方差很大為1.885，y2方差相對很小為0.115，且y1和y2不相關

#可以認為原來二維數(shù)據(jù)的變化都體現(xiàn)在y1一個維度上了2024/9/510

圖9－2二維正態(tài)分布模擬數(shù)據(jù)的主成分示意圖2024/9/511

圖9－3坐標軸旋轉以后的散點圖一般,設總體的p個主成分為：

第i

個主成分yi的方差為y1

是

X的一切線性組合中方差最大者；y2

是與y1

不相關是

X的一切線性組合中方差最大者；······2024/9/512

9.2.2主成分的計算2024/9/513

9.2.3主成分的主要性質2024/9/514

設總體X的p個主成分所成向量為：性質1

：性質2：性質3：載荷與載荷矩陣：2024/9/515

9.2.4主成分個數(shù)的確定2024/9/516

第i

個主成分yi的方差貢獻率為：通常取使得累積貢獻率滿足的最小的k為主成分個數(shù)．有的文獻取累積貢獻率首次超過85%的k．

9.2.5變量的標準化及意義從總體協(xié)方差矩陣Σ出發(fā)做主成分分析傾向于反映方差大的變量的信息，會出現(xiàn)“大數(shù)吃小數(shù)”的現(xiàn)象．為了均等地對待每一個原始變量，常常將各原始變量作標準化處理：標準化后的總體

的協(xié)方差矩陣恰好是原總體X的相關系數(shù)矩陣

ρ．綜上討論，既可從Σ出發(fā)，也可以從ρ出發(fā)做主成分分析，考慮到現(xiàn)實經(jīng)濟意義，后者用得更多.2024/9/517

9.3樣本主成分實際問題中Σ和ρ往往是未知的,需要用樣本的協(xié)方差矩陣

S和樣本的相關系數(shù)矩陣

R來估計：2024/9/518

9.3.1樣本主成分的性質和計算設S的p個特征值為，對應的單位正交特征向量為，則樣本的第

i個主成分為性質1

：性質2

：性質3

：性質4：2024/9/519

9.3.2主成分分析的步驟

實際問題中更常用的是從樣本相關系數(shù)矩陣R

出發(fā)進行主成分分析,方法是用

R

替換

S,其余操作不變,其步驟可歸納為：將原始樣本標準化求樣本的相關系數(shù)矩陣

R

求R的特征值以及對應的單位正交特征向量按主成分累積貢獻率超過80%確定主成分的個數(shù)k,并寫出主成分表達式為(5)對分析結果做統(tǒng)計意義和實際意義兩方面的解釋.2024/9/520

主成分分析特別說明：2024/9/521

9.3.2

(續(xù))主成分相關的R函數(shù)1.

princomp函數(shù)(這是主成分分析最常用的函數(shù))princomp(x,cor=FALSE,score=TRUE,…)2.summary函數(shù)(提取主成分的信息)summary(object,loadings=FALSE,…)3.loadings函數(shù)(顯示主成分的載荷陣)loadings(object)4.

predict函數(shù)(預測主成分的值)predict(object,newdata,…)5.

screeplot函數(shù)(畫出主成分的碎石圖)screeplot(object,type=c(“barplot”,“l(fā)ines”,…)2024/9/522

2024/9/523表9-1給出了52名學生的數(shù)學

(x1)、物理

(x2)、化學

(x3)、語文

(x4)、歷史

(x5)和英語

(x6)成績，對其進行主成分分析．例9.1學生六門課成績數(shù)據(jù)的主成分分析>setwd("C:/data")#設定工作路徑>d9.1<-read.csv(“exam9.1.csv”,header=T)#讀入數(shù)據(jù)>R=round(cor(d9.1),3);R#樣本相關系數(shù)陣保留三位小數(shù)

x1x2x3x4x5x6x11.0000.6470.696-0.561-0.456-0.439x2

0.6471.0000.573-0.503-0.351-0.458x30.6960.5731.000

-0.380-0.274-0.244x4-0.561-0.503-0.3801.0000.8130.835x5-0.456-0.351-0.2740.8131.0000.819x6-0.439-0.458-0.244

0.8350.8191.000解

先讀取數(shù)據(jù)，求樣本相關系數(shù)矩陣，R程序如下：表6-017個地質勘探點樣品的標準化數(shù)據(jù)2024/9/524

學號x1x2x3x4x5x6165617284817927777766470553676349656757478847562716456671675265576831007941675078694975163558678453586656……………………………………4599100995363604678685275746647729073768079486964606874804952626510096100507072567482745172747588918652687470878783表

9－152名學生六門課程成績數(shù)據(jù)2024/9/525

易見，文科三門課程語文(x4)、歷史(x5)和英語(x6)相關性較強；理科三門課程數(shù)學(x1)、物理(x2)和化學(x3)相關性也較強．可以進一步作主成分分析，求樣本相關矩陣的特征值和主成分載荷．

由下面的R程序運行結果可知主成分的標準差，即相關系數(shù)矩陣的六個特征值開方各為：同時前兩個主成分的累積貢獻率為0.618+0.210=0.829，已經(jīng)超過80%，所以取兩個主成分就可以了．2024/9/526>PCA9.1=princomp(d9.1,cor=T)

#用樣本相關系數(shù)陣做主成分分析>PCA9.1Call:princomp(x=d9.1,cor=T)Standarddeviations:Comp.1Comp.2Comp.3Comp.4Comp.5Comp.61.926

1.1240.6640.5200.4120.3836variablesand52observations.>summary(PCA9.1,loadings=T)

#列出主成分分析結果Importanceofcomponents:

Comp.1Comp.2Comp.3Comp.4Comp.5Comp.6Standarddeviation1.926

1.124

0.6640.5200.4120.383ProportionofVariance

0.6180.2100.0730.0450.0280.024CumulativeProportion

0.6180.8290.9020.9470.9761.000Loadings:

Comp.1Comp.2Comp.3Comp.4Comp.5Comp.6x1-0.412-0.3760.2160.788-0.145x2-0.381-0.357-0.806

-0.1180.212

0.141x3-0.332-0.563

0.467-0.588x40.461-0.279

0.599-0.590x50.421-0.415-0.250-0.738

-0.205x60.430-0.4070.1460.1340.2220.749第一主成分和第二主成分分別為：參見教材分析，可將將它們分別理解為“課程差異主成分”

和“課程均衡主成分”.2024/9/5

27>round(predict(PCA9.1),3)

#作預測,計算主成分得分并解釋>screeplot(PCA9.1,type=“l(fā)ines”)

#畫線型碎石圖(見圖9-4)>load=loadings(PCA9.1)

#提取主成分載荷矩陣為load>plot(load[,1:2],xlim=c(-0.6,0.6),ylim=c(-0.6,0.6))#作散點圖>rnames=c(“數(shù)學”,“物理”,“化學”,“語文”,“歷史”,“英語”)#命名>text(load[,1],load[,2],labels=rnames,adj=c(-0.3,1.5))#用中文為散點圖標注>abline(h=0,v=0,lty=3)

#用虛線劃分四個象限(見圖9-5)2024/9/5

28圖9－452名學生六門課程成績的主成分線型碎石圖2024/9/5

29圖9－5前兩個主成分的載荷散點圖2024/9/5

30圖9－652名學生成績數(shù)據(jù)的雙坐標散點圖>biplot(PCA9.1,scale=0.5)

#繪制52個樣本點關于前兩個主成分的散點圖2024/9/531

由于第一主成分是文理課程差異因子，理科課程在第一主成分上的載荷絕對值大且取負值，文科課程在第一主成分上的載荷絕對值大且取正值，因此圖中Comp.1軸方向靠左的樣本點，如6，7和45號樣本點，對應理科成績好、文科成績差的學生；相對的Comp.1軸方向靠右的樣本點，如30和49號樣本點，對應文科成績好、理科成績差的學生．又第二主成分表示課程均衡因子，在圖中Comp.2軸方向靠下的樣本點，如26，33號樣本點，對應各科成績都較好學生，相對的Comp.2軸方向靠上的樣本點，如3，5和8號樣本點，對應各科成績都較差的學生，而居中的樣本點，如42，24和39號樣本點，對應各科成績都屬于中等且差異不大的學生．這樣就可以對52名學生按對應樣本點所在的位置進行大致分類．2024/9/532

在某沉積盆地一坳陷區(qū)的17個取樣點經(jīng)勘探測定了六個地質變量：x1為有機碳（%）；x2為生油層埋深（米）；x3油層孔隙度（%）；x4為儲層厚度（米）；x5為地下水含碘量（p.p.m）；x6為地下水礦化度(克/升)，見表9-2．要求根據(jù)這些數(shù)據(jù)進行主成分分析.例9.2石油勘探樣品數(shù)據(jù)分析(數(shù)據(jù)exam9.2)表6-017個地質勘探點樣品的標準化數(shù)據(jù)2024/9/533

點號x1x2x3x4x5x61-0.9142-0.7119-0.9293-0.4385-0.57100.73612-0.3095-0.5206-1.3309-0.2764-0.57100.57143-1.0654-0.71190.2756-0.7626-1.09570.90074-1.3073-0.95111.25740.3718-1.09571.394650.1743-0.47270.3203-0.9895-0.0463-0.25186-0.8235-0.59230.40951.3441-0.83330.406870.90002.1583-0.1260-0.85981.7901-1.89838-0.0071-0.3532-1.4201-1.0219-0.0463-0.581191.20231.6799-0.7508-0.60052.3148-1.2397100.1743-0.3532-0.97391.3441-0.04630.2421112.26061.44070.72192.64050.7407-1.075012-1.4282-0.95110.0079-0.7950-1.09571.065313-0.3397-0.52062.1499-0.1144-0.57100.4068140.7790-0.23361.19700.69590.21610.9104150.41620.72321.0789-0.30880.47840.745716-0.6118-0.71190.36490.0477-0.57101.5593170.90001.08200.1418-0.27641.0031-0.5811表

9－217個地質勘探點樣品的標準化數(shù)據(jù)2024/9/534>setwd("C:/data")

#設定工作路徑>d9.2<-read.csv("exam9.2.csv",header=T)#將exam9.2數(shù)據(jù)讀入到d9.2中>R=round(cor(d9.2),3);R

#求樣本相關系數(shù)矩陣

x1x2x3x4

x5x6x11.0000.8400.003

0.3470.839-0.747x20.8401.000-0.051

0.077

0.939-0.839x30.003-0.0511.0000.259-0.1640.285x40.3470.0770.259

1.000-0.0370.022x50.8390.939-0.164-0.0371.000-0.827x6-0.747-0.8390.2850.022-0.8271.000

易見，x2與x5相關性最強，其絕對值在0.9~0.95，x1與x2，x1與x5，x2與x6，x5與x6的相關性較強，其絕對值在0.8~0.9，說明六個變量之間確實存在較強的相關關系，應當進行“降維”處理，可以作主成分分析．2024/9/535>options(digits=3)

#設置小數(shù)點位數(shù)為3>PCA9.2=princomp(d9.2,cor=T,scores=T);PCA9.2#作主成分分析Call:princomp(x=d9.2,cor=T,scores=T)Standarddeviations:Comp.1Comp.2Comp.3Comp.4Comp.5Comp.6

1.8851.1700.8600.4300.3400.197>summary(PCA9.2,loadings=T)

#列出主成分分析結果Importanceofcomponents:

Comp.1Comp.2Comp.3Comp.4Comp.5Comp.6Standarddeviation

1.8851.1700.8600.43010.33990.19653ProportionofVariance0.5920.2280.1230.03080.01930.00644CumulativeProportion0.592

0.8200.9430.97430.99361.00000Loadings:

Comp.1Comp.2Comp.3Comp.4Comp.5Comp.6x10.4850.2390.2910.7350.274x20.510-0.166-0.587

0.600x30.646-0.728-0.181x40.7020.640-0.254-0.153x50.509-0.1540.409-0.187-0.713x6-0.484

0.1590.837-0.1180.155前兩個主成分的累積貢獻率為0.592+0.228=0.82，已經(jīng)超過80%，所以只需取兩個主成分．第一主成分和第二主成分各為（為簡明起見，樣本主成分表達式中的所有“*”省略，以下同）：2024/9/536

四個變量x1

(有機碳)，x2(生油層埋深)，x5(地下水含碘量)和x6(地下水礦化度)在主成分z1上載荷較大，故第一主成分z1可解釋為“生油條件”主成分；第二主成分z2與x3(油層孔隙度)和x4(儲層厚度)這兩個變量關系特別密切，可解釋為“儲油條件”主成分．這樣的分析結果與石油地質理論是相符合的．2024/9/537>screeplot(PCA6.1,type=“l(fā)ines”)#畫碎石圖，用直方圖類型（見圖9-7）圖9-717個石油地質勘測點樣本數(shù)據(jù)的主成分碎石圖2024/9/538用主成分載荷矩陣前兩列數(shù)據(jù)作主成分載荷散點圖(見圖9-8)，R程序如下：．load=loadings(PCA9.2)#提取主成分載荷矩陣plot(load[,1:2],xlim=c(-0.5,1),ylim=c(-0.2,0.8))

#作散點圖rnames=c(“x1有機碳”,“x2生油層埋深”,“x3油層孔隙度”,“x4儲層厚度","x5地下水含碘量","x6地下水礦化度")#見圖9－8text(load[,1],load[,2],labels=rnames,cex=0.8,adj=c(-0.1,0.6))

#用中文為散點標號abline(h=0,v=0,lty=3)

#用虛線劃分象限六個變量在主成分z1和z2坐標面上的載荷散點圖表明了兩個主成分z1和z2具有明顯的“生油”和“儲油”傾向特征．2024/9/539

圖9-8兩個主成分的載荷散點圖2024/9/540>A=round(PCA9.2$scores,3)

#計算主成分得分，取3位小數(shù)>B=round(apply(A[,1:2],1,crossprod),2)

#按行加總前2個主成分上的載荷平方>cbind(A,"綜合得分“=B,"排名“=rank(B))

#按列合并主成分得分、綜合得分和排名

Comp.1Comp.2Comp.3Comp.4Comp.5Comp.6綜合得分

排名

[1,]

-1.333

-1.1500.6350.136-0.150

-0.015

3.10

9

[2,]

-0.798

-1.1751.0730.252

0.1530.261

2.024

[3,]

-1.901

-0.540

-0.429

-0.187-0.047

0.287

3.9112

[4,]

-2.421

0.956

-0.424

-0.124-0.402-0.125

6.7814

[5,]

-0.037

-0.640

-0.669

-0.357

0.749

-0.144

0.412

[6,]

-1.2091.0340.836

-0.626-0.472-0.190

2.536

[7,]3.560

-1.001

-0.862

-0.504-0.495

0.093

13.6816

[8,]0.265

-1.930

0.606

-0.359

0.514

0.0033.8011

[9,]3.474

-1.110

-0.2160.441-0.255-0.37513.30

15[10,]0.0420.2641.8140.143-0.052-0.2290.071[11,]3.0332.6361.183

-0.323

0.189

0.30416.1517[12,]

-2.269

-0.814

-0.245

-0.119-0.2010.0955.8113[13,]

-1.0471.250

-1.458

-0.573

0.258-0.144

2.667[14,]

-0.0401.500

-0.2940.678

0.464-0.181

2.255[15,]0.4510.589

-1.0980.725-0.161

0.258

0.553[16,]

-1.6920.285

-0.0300.668-0.095

0.000

2.948[17,]1.922

-0.152

-0.4200.129

0.003

0.102

3.7210

11號樣本點綜合排名最高，為17分；7號和9號排名次之，分別為16分和15分；之后樣本點得分排名從高到低依次為4、12、3、8和17號．

利用函數(shù)biplot來繪制它們在z1和z2構成的坐標面

z1Oz2上的散點圖，并且加入六個變量在同一坐標面

z1Oz2上的載荷散點圖，得到所謂的“雙坐標”散點圖（見圖9-9）．

借助該圖可以對17個勘測樣本點進行大致分類：11號樣本點獨居右上，它在

“生油”主成分z1和“儲油”主成分z2上得分均高，應該首先重點關注．7、9號樣本點相鄰且最靠右，且在z1得分很高，可合為一類，次重點考慮；此外，在z1和z2上至少有一個得分較高的3、4、8、12和17號樣本點也應該重點考察．這與上面的綜合得分和排名一致．2024/9/541

2024/9/542

圖9－917個石油地質勘測點樣本數(shù)據(jù)的雙坐標散點圖>biplot(PCA9.2,scale=0.5)#繪制17個樣本點和6個變量對z1和z2的散點圖9.4案例:主成分回歸分析

案例9.1(數(shù)據(jù)文件為case9.1)表9-3給出了2019年全國31個地區(qū)相關數(shù)據(jù)．它們分別為：貨運量x1(萬噸)，貨物周轉量x2(億噸公里)，GDPx3(億元)，人均GDPx4(元)，城鎮(zhèn)居民人均可支配收入y(元)．根據(jù)這些數(shù)據(jù)做線性回歸分析和主成分回歸分析，并比較它們的異同.2024/9/543

2024/9/544

表9-32019年全國部分地區(qū)貨運量、GDP、人均可支配收入等數(shù)據(jù)城市x1x2x3x4y北京22808108935371.316422073848.5天津50093266214104.39037146118.9河北2424451356335104.54634835737.7山西192192546617026.74572433262.4內蒙古182702458717212.56785240782.5遼寧178253892124909.55719139777.2吉林43193180311726.84347532299.2黑龍江50475161513612.73618330944.6………………………………貴州83402123516769.34643334404.2云南122727155223223.84794436237.7西藏40251541697.84890237410.0陜西154749348225793.26664936098.2甘肅6361024968718.33299532323.4青海149453982966.04898133830.3寧夏425116513748.55421734328.5新疆844231948135977解(1)先做線性回歸分析，R程序及結果如下：>setwd("C:/data")

#設定工作路徑>c9.1<-read.csv("case9.1.csv",header=T)

#將數(shù)據(jù)讀入到c9.1中>options(digits=3)

#取三位有效數(shù)字>lmc9.1<-lm(y~1+x1+x2+x3+x4,data=c9.1)>summary(lmc9.1)從輸出結果（見教材）可以看出，回歸方程是非常顯著的，R2為0.923，模型擬合效果很好，但x1、x2和x3的回歸系數(shù)沒有通過顯著性檢驗(在0.05的顯著性水平下)．回歸方程為：然后作逐步回歸，R程序及結果如下：>summary(step(lmc9.1))

回歸方程和回歸系數(shù)均顯著，R2為0.921，逐步回歸方程為：2024/9/545

(2)再作主成分回歸分析,先求樣本相關系數(shù)陣：>R=round(cor(c9.1[,2:6]),3);R

#求樣本相關系數(shù)矩陣發(fā)現(xiàn)x4與

y高度相關，

x1

,

x2

,

x3相關性較強，可用主成分降維>c9.1pr<-princomp(~x1+x2+x3+x4,data=c9.1,cor=T)#使用公式法>summary(c9.1pr,loadings=T)

前兩個主成分累積貢獻率已達88%，故選擇前兩個主成分>pre<-predict(c9.1pr)

#計算主成分得分>c9.1$z1<-pre[,1];c9.1$z2<-pre[,2]>lmpr<-lm(y~z1+z2,data=c9.1)

#做y關于主成分z1和z2的回歸>summary(lmpr)輸出結果顯示：y關于兩個主成分z1和z2的回歸方程和三個回歸系數(shù)均是非常顯著的，R2為0.887，主成分回歸方程為：2024/9/546

(3)還利用主成分與原來自變量間的關系將主成分還原為原來的自變量，參見參考文獻[4]．R程序及結果如下：>beta<-coef(lmpr);A<-loadings(c9.1pr)[,1:2]>x.bar<-c9.1pr$center;x.sd<-c9.1pr$scale>coef<-A%*%beta[2:3]/x.sd>beta0<-beta[1]-x.bar%*%coef>c(beta0,coef)[1]2.36e+04-2.86e-023.98e-014.64e-022.52e-01

由輸出結果知主成分z1和z2還原為原始變量后的回歸方程為：可將它和最初得到的回歸方程進行比較．前者是從主成分回歸方程（方程和回歸系數(shù)均顯著）變形而來，更合理，預測效果也更好．2024/9/547

2024/9/548主編：費宇，魯筠中國人民大學出版社，2024年9/5/2024主編：費宇50第10章因子分析因子分析（factoranalysis）最早起源于KarlPearson和ChalesSpearman等人關于智力的定義和測量工作，因子分析的基本目的是，只要可能，就用少數(shù)幾個潛在的不能觀察的隨機變量（稱為因子）去描述許多個隨機變量之間的協(xié)方差關系．從這點上看，因子分析與主成分分析有相似之處，但因子分析中的因子是不可觀察的，也不必是相互正交的變量．因子分析可以視為主成分分析的一種推廣，它的基本思想是：根據(jù)相關性大小把變量分組，使得組內的變量相關性較高，但不同組的變量相關性較低，則每組變量可以代表一個基本結構，稱為因子，它反映已經(jīng)觀測到的相關性．因子分析可以用來研究變量之間的相關關系，稱為R型因子分析；也可以用來研究樣品之間的相關關系，稱為Q型因子分析．二者雖然形式上有所不同，但數(shù)學處理上是一樣的，所以本章只介紹R型因子分析．9/5/2024主編：費宇51第10章因子分析10.1正交因子模型10.2因子模型的估計10.3因子正交旋轉10.4因子得分10.5因子分析小結10.6案例分析9/5/2024主編：費宇5210.1正交因子模型1.模型定義設p維隨機向量的期望為,方差-協(xié)方差矩陣為Σ,

假定X線性地依賴于少數(shù)幾個不可觀測的隨機變量f1,…,fm(m<p)和p個附加的方差源

ε1,…,εp,

一般稱f1,…,fm為公因子，稱ε1,…,εp為特殊因子，或誤差.

1.模型定義9/5/2024主編：費宇53那么，因子模型為1.模型定義引入矩陣符號，記

那么因子模型(10.1)可以寫為9/5/2024主編：費宇541.

模型定義其中aij稱為第i個變量在第j個因子上的載荷,矩陣A稱為載荷矩陣.我們假定

如果模型(10.2)滿足假定(10.3)，則稱該模型為正交因子模型，如果F的各個分量相關，即

不是單位陣，則相應的模型稱為斜交因子模型，本書只討論正交因子模型.9/5/2024主編：費宇551.

模型定義從正交因子模型容易求得X的協(xié)方差9/5/2024主編：費宇561.模型定義同樣，容易求得由(10.4)可得該式說明xi的方差由兩部分構成：m個公因子和一個特殊因子，其中表示第j個公因子對xi的方差貢獻，而Фi是第i個特殊因子對xi的方差貢獻，稱之為特殊度.記

hi2=ai12+ai22+…+aim2，它表示m個公因子對變量xi的方差貢獻總和，稱之為第i個共同度，它是載荷矩陣A的第i行元素平方和.9/5/2024主編：費宇571.模型定義由(10.5)可得上式說明aij表示變量xi與公因子fj的協(xié)方差.另一方面，我們也可以考慮某個公因子fj對各個變量x1,…,xp的影響，采用來度量這個影響的大小，bj2是載荷矩陣A第j列元素的平方和，稱之為公因子fj對各p個變量的方差貢獻，bj2越大，表示fj對各p個變量的影響越大，它可以作為公因子fj重要性的一個度量.9/5/2024主編：費宇581.模型定義需要指出的是，當

m>1時，因子模型是不唯一的,設T為m×m正交矩陣，即TTT=TTT=I，模型(10.2)可改寫為

式中，A*=AT,F*=TTF9/5/2024主編：費宇591.模型定義注意到即F*也滿足(10.3)，顯然因子F與F*有相同的統(tǒng)計性質，但相應的載荷矩陣A與A*是不相同的，但它們產(chǎn)生相同的方差-協(xié)方差矩陣Σ

，即9/5/2024主編：費宇601.模型定義一方面，因為F*=TTF，即F*是由F經(jīng)正交變換得到,而A*=AT，即A*=(a*ij)是由A=(aij)經(jīng)正交變換得到,另一方面，由(10.11)易知，變量xi的共同度為即正交變換不改變公因子的共同度.9/5/2024主編：費宇6110.2

因子模型的估計建立因子模型首先要估計載荷矩陣及特殊方差，常用的估計方法有主成分法、主因子法和極大似然法.1.

主成份法設Σ的特征值為λ1,

λ2,…,λp(λ1≥λ2≥…≥λp≥0)，e1,e2,…,ep為對應的標準正交化特征向量，那么Σ可以寫為9/5/2024主編：費宇621.

主成份法這個分解是公因子個數(shù)為p，特殊因子方差為0的因子模型的方差-協(xié)方差矩陣結構形式，即雖然上式給出的Σ因子分析表達式是精確的，但實際應用中沒有價值，因為因子分析的目的是要尋找少數(shù)m(m<p)個公因子解釋原來p個變量的協(xié)方差結構，所以，采用主成分分析的思想，如果Σ的最后p-m個特征值很小，在(10.13)中略去λm+1em+1eTm+1+…+λpepeTp對Σ的貢獻，9/5/2024主編：費宇631.

主成份法于是得這里假定了(10.2)中的特殊因子是可以在Σ的分解中忽略的，如果特殊因子不能忽略，那么它們的方差可以取Σ-AAT的對角元，9/5/2024主編：費宇641.

主成份法此時有其中9/5/2024主編：費宇651.

主成份法實際應用中Σ是未知的，通常用它的估計，即樣本協(xié)方差矩陣S來代替，考慮到變量的量綱差別，往往需要將數(shù)據(jù)標準化，這樣求得的樣本方差-協(xié)方差矩陣就是原來數(shù)據(jù)的相關系數(shù)矩陣R，所以可以從R出發(fā)來估計因子載荷矩陣和特殊因子的方差.9/5/2024主編：費宇661.

主成份法設R的特征值為 ,

為對應的標準正交化特征向量，設

m<p，則由R出發(fā)因子模型的載荷矩陣的估計為特殊因子的方差фi的估計為9/5/2024主編：費宇671.

主成份法這時，共同度hi2的估計為變量xi與公因子fj協(xié)方差的估計為，公因子fj對各個變量的貢獻bj2的估計為9/5/2024主編：費宇681.

主成份法那么，如何確定公因子數(shù)目m呢？可以仿照主成分分析的思想，比如尋找m使得來確定公因子數(shù)m.9/5/2024主編：費宇692.

主因子法假定原始向量X的各分量已作了標準化變換.如果其滿足正交因子模型，則有9/5/2024主編：費宇70則稱為X的約相關矩陣.其中，

R為X的相關矩陣.令2.

主因子法9/5/2024主編：費宇71

中的對角線元素是

，而不是1，非對角線元素和R中是完全一樣的，并且是一個非負定矩陣.設是特殊方差的一個合適的初始估計，則約相關矩陣可估計為：2.

主因子法9/5/2024主編：費宇72其中，，是的初始估計.又設的前個特征值依次為，相應的正交單位特征向量為，則A的主因子解為：2.

主因子法9/5/2024主編：費宇73由此我們可以重新估計特殊方差，的最終估計為：如果我們希望求得擬合程度更高的解，則可以采用迭代的方法，即利用式(10.26)中的再作為特殊方差的初始估計，重復上述步驟，直至解穩(wěn)定為止．2.

主因子法9/5/2024主編：費宇74特殊(或共性)方差的常用初始估計方法有：(1)取

，其中是的第個對角線元素，此時共性方差的估計為

，它是

和其他個變量間樣本復相關系數(shù)的平方，該初始估計方法最為常用．(2)取

，此時.(3)取

，此時

，得到的是一個主成分解.3.

極大似然法9/5/2024主編：費宇75設公共因子，特殊因子，且相互獨立，則必然有原始向量．由樣本計算得到的似然函數(shù)是和的函數(shù)．由于，故似然函數(shù)可更清楚地表示為．記的極大似然估計為，即有可以證明，而和滿足方程組：3.

極大似然法9/5/2024主編：費宇76式中，由于A的解是不唯一的，故為了得到唯一解，可附加計算上方便的唯一性條件：是對角矩陣.3.

極大似然法9/5/2024主編：費宇77方程組（10.28）的和一般可用迭代方法解得.對極大似然解，當因子數(shù)增加時，原來因子的估計載荷及對x的貢獻將發(fā)生變化，這與主成分解及主因子解不同．10.3

因子正交旋轉在第10.1節(jié)我們已經(jīng)看到，滿足方差結構Σ

=AAT+Ф的因子模型并不惟一，模型的公因子與載荷矩陣不惟一.如果F是模型的公因子，A是相應的載荷矩陣，而T是m×m正交矩陣，則F*=TTF也是公因子，相應的載荷矩陣為A*=AT，A*也滿足Σ=A*A*T+Ф這說明，公因子和因子載荷矩陣作正交變換后，并不改變共同度，我們稱因子載荷的正交變換和伴隨的因子正交變換為因子正交旋轉.9/5/2024主編：費宇7810.3

因子正交旋轉設

是用某種方法（比如主成分法）得到的因子載荷矩陣的估計，T為

m×m正交陣,則是旋轉載荷矩陣.問題是：為什么要進行因子旋轉？其目的是什么？9/5/2024主編：費宇7910.3

因子正交旋轉如果初始載荷不易解釋時，就需要對載荷作旋轉，以便得到一個更簡單的結構.最理想的情況是這樣的載荷結構，每個變量僅在一個因子上有較大的載荷，而在其余因子上的載荷比較小，至多是中等大小，這樣公因子fi的具體含義可由載荷較大的變量根據(jù)具體問題加以解釋.如何進行因子旋轉尋找一個簡單結構的載荷矩陣，這里不作詳細介紹.9/5/2024主編：費宇8010.4

因子得分在因子分析中，雖然我們關心模型中載荷矩陣的估計和對公因子的解釋，但對于公因子的估計，即因子得分，有時也是需要的.但是因子得分的計算并不同于通常意義下的參數(shù)估計，而是對不可觀測的因子fj取值的估計，下面介紹用加權最小二乘法估計因子得分.9/5/2024主編：費宇811.

加權最小二乘法給定因子模型X=μ+AF+ε,假定均值向量μ,載荷矩陣A和特殊方差陣Ф已知，把特殊因子ε看作誤差，因為Var(εi)=фi(i=1,2,…,p)未必相等,所以我們用加權最小二乘法估計公因子F.首先將因子模型

(10.2)改寫為9/5/2024主編：費宇821.

加權最小二乘法兩邊左乘Ф-1/2得記X*=Ф-1/2(X-μ),A*=Ф-1/2A,ε*=Ф-1/2ε,則上式可以寫成注意到E(ε*)=Ф-1/2E(ε)=0,

cov(ε*)=E(ε*ε*T)=Ф-1/2E(εεT)Ф-1/2=I9/5/2024主編：費宇831.

加權最小二乘法所以(10.32)是經(jīng)典的回歸模型，由最小二乘法知F的估計為

實際中,A,

Ф和μ都是未知的,通常用它們的某種估計來代替,比如我們采用正交旋轉后的載荷矩陣A的估計,和樣本均值 ,分別代替A,Ф和μ9/5/2024主編：費宇841.

加權最小二乘法于是可得對應于xj的因子得分9/5/2024主編：費宇852.

回歸法在正交因子模型中，假設服從(m+p)元正態(tài)分布，用回歸預測方法可將估計為：9/5/2024主編：費宇862.

回歸法在實際應用中，可用，和分別代替上式中的，和來得到因子得分．樣品的因子得分9/5/2024主編：費宇873.

綜合因子得分9/5/2024主編：費宇88個因子任意若干個取相反符號，特別是全部取相反符號仍然滿足因子分析模型，所以仍然可以作為因子。3.

綜合因子得分以各因子的方差貢獻率為權重，由各因子的線性組合得到綜合評價指標函數(shù)：式中，9/5/2024主編：費宇893.

綜合因子得分9/5/2024主編：費宇90那么這樣的因子得分函數(shù)將會有種不同的組合。所以這樣的因子得分實際上是不好解釋的，此外，使用不同的因子旋轉會得到不同的因子，從而綜合評價函數(shù)也就不同，哪一個才是對的呢？還有，綜合起來表示的是什么因子呢？所以，因子綜合得分是沒有合理的解釋的。例10.1數(shù)據(jù)文件為eg9.1前面第9章例9.1表9-1給出了52名學生的數(shù)學（x1）、物理（x2）、化學（x3）、語文（x4）、歷史（x5）和英語（x6）成績，試進行學生成績的因子分析.解:采用R軟件對樣本數(shù)據(jù)進行因子分析，首先計算樣本數(shù)據(jù)的相關系數(shù)矩陣，觀察各變量之間的相關性.

R程序及結果如下：9/5/2024主編：費宇91例10.1數(shù)據(jù)文件為eg9.1#假設已經(jīng)讀取了52名學生成績數(shù)據(jù)cor(X)#計算樣本數(shù)據(jù)的相關系數(shù)矩陣x1x2x3x4x5x6x11.000.650.70-0.56-0.46-0.44x20.651.000.57-0.50-0.35-0.46x30.700.571.00-0.38-0.27-0.24x4-0.56-0.50-0.381.000.810.83x5-0.46-0.35-0.270.811.000.82x6-0.44-0.46-0.240.830.821.009/5/2024主編：費宇92例10.1數(shù)據(jù)文件為eg9.1從樣本數(shù)據(jù)各變量的相關系數(shù)上可以看出,x4、x5和x6之間存在較強的相關性.為了消除各變量之間的相關性，下面分別采用R軟件中基于極大似然法的因子分析函數(shù)factanal()和基于主成分法的因子分析函數(shù)factpc()對數(shù)據(jù)進行因子分析提取因子.

R程序及結果如下：9/5/2024主編：費宇93例10.1數(shù)據(jù)文件為eg9.1#極大似然法做因子分析factanal(X,factors=2,rotation="none")Call:factanal(x=X,factors=2,rotation="none")Uniquenesses:x1x2x3x4x5x60.230.460.330.150.210.15Loadings:Factor1Factor2x1-0.680.56x2-0.600.43x3-0.490.66x40.920.10x50.860.24x60.880.279/5/2024主編：費宇94例10.1數(shù)據(jù)文件為eg9.1

Factor1Factor2SSloadings3.401.07ProportionVar0.570.18CumulativeVar0.570.74

Testofthehypothesisthat2factorsaresufficient.Thechisquarestatisticis3.6on4degreesoffreedom.Thep-valueis0.46#主成分法做因子分析library(mvstats)#加載mvstats包fac=factpc(X,2)fac9/5/2024主編：費宇95例10.1數(shù)據(jù)文件為eg9.1$VarsVarsVars.PropVars.CumFactor13.7100.618361.83Factor21.2620.210482.87$loadingsFactor1Factor2X1-0.79370.4224x2-0.73420.4008x3-0.63970.6322x40.88830.3129x50.81010.4661x60.82850.45679/5/2024主編：費宇96例10.1

數(shù)據(jù)文件為eg9.1從上述極大似然法和主成分法得出的因子分析結果上可以看出，極大似然法前兩個因子累計貢獻率只有74%，而主成分法累計貢獻率達到了82.87%，說明主成分法效果比極大似然分析法效果好，其原因在于，極大似然法做因子分析要求數(shù)據(jù)樣本要服從多元正態(tài)分布，但在實際中大多數(shù)數(shù)據(jù)都很難滿足多元正態(tài)要求。接下來為了更好地解釋因子的含義，我們基于主成分法采用方差最大化作因子正交旋轉。R程序及結果如下：9/5/2024主編：費宇97例10.1

數(shù)據(jù)文件為eg9.1fac1=factpc(X,2,rotation="varimax")#用主成分法采用方差最大化作因子正交旋轉Fac