一、主成分分析是利用降維的方法,在損失很少信息量很少的前提下,把多個指標轉換為幾個綜合指標的多元統計方法。通常把轉化的綜合指標稱為主成分。
二、基本原理
在對某一事物進行研究時,為了更全面、准確地反應事物的特征及其發展規律人們通常考慮一起有關系的多個指標,也叫變量。
三、主成分分析步驟
1、根據問題選取初始變量
2、根據初始變量特性判斷由協方差矩陣求主成分還是由相關陣求主成分
3、求協方差矩陣或相關矩陣的特征根和相應的標准特征向量
4、判斷是否存在明顯的多重共線性,若存在則回到第一步。
5、得到主成分的表達式並確定主成分個數,選取主成分
6、結合主成分對研究問題進行深入的分析
四、主成分上機實現
1、princomp函數為主成分分析的主要函數
其用法為princomp(formula,data=NULL,subset,na.action.....)
其中formula是沒有響應變量的公式,data是數據框,或者
princomp(x,cor=FALSE,scores=TURE,covmat=NULL,subset=rep(TURE,nrow(as.matrix(x))),....)
其中x是主要分析的數據,以數值矩陣或數據框的形式給出,cor為邏輯變量,當cor為TURE表示用樣本的協方差陣做主成分分析,反之用協方差陣,covmat是協方差陣,如果數據不由x提供就是由其提供。
2、summary函數提取出相應的信息
summary(object,loadings=FALSE,cutoff=0,1...)
objict表示princom()得到的對象,loadings為TRUE時,表示顯示loadings的內容,反之顯示。
3、loadings()顯示主成分分析或因子分析中的內容
4、predict()
predict(objict)用來預測主成分的值。
5、screeplot()函數用於畫碎石圖,使用格式為screeplot(x,npcs=min(10,length(x$sdev)),type=c("barplot","lines"))
npcs是主成分的個數,type描述其類型,分別為直方圖和直線圖的類型。
6、biplot()函數主要畫出主成分的散點圖和原坐標在主成分下的方向,使用格式為
biplot(x,choice=1:2,scale=1,pc.boplot=FALSE....)
choice是選擇主成分,默認為1,2主成分,bc.plot為邏輯變量默認為FALSE。
7、實例
下表是某工業部門13個行業,分別是冶金(1)、電力(2)、煤炭、化學、機械、建材、森工、食品、紡織、縫紉、皮革、造紙、和文化教育藝術品,8個指標分別為
年末固定資產凈產值(x1,萬元)、職工人數、工業總產值、全員勞動生產率、百元固定原值實現產值、資金利潤率(x6)、標准燃料消費量和能源利用效果。
(1)試用主成分分析方法確定8個指標的幾個主成分,並對其進行解釋
(2)利用主成分得分對13個行業進行排序和分類。
X1 X2 X3 X4 X5 X6 X7 X8 1 90342 52455 101091 19272 82.0 16.1 197435 0.172 2 4903 1973 2035 10313 34.2 7.1 592077 0.003 3 6735 21139 3767 1780 36.1 8.2 726396 0.003 4 49454 36241 81557 22504 98.1 25.9 348226 0.985 5 139190 203505 215898 10609 93.2 12.6 139572 0.628 6 12215 16219 10351 6382 62.5 8.7 145818 0.066 7 2372 6572 8103 12329 184.4 22.2 20921 0.152 8 11062 23078 54935 23804 370.4 41.0 65486 0.263 9 17111 23907 52108 21796 221.5 21.5 63806 0.276 10 1206 3930 6126 15586 330.4 29.5 1840 0.437 11 2150 5704 6200 10870 184.2 12.0 8913 0.274 12 5251 6155 10383 16875 146.4 27.5 78796 0.151 13 14341 13203 19396 14691 94.6 17.8 6354 1.574
1 1 w=read.csv("afh.csv") 2 2 w 3 3 attach(w) 4 4 X1 5 5 #### 作主成分分析 6 6 w.pr<-princomp(w, cor=TRUE)#利用協方差作主成分分析 7 7 w.pr 8 8 #### 並顯示分析結果
Standard devi....其中可以看出其標准差,即主成分的方差開方,也就是相應特征值lambda1、lambda2、lambda3到lambda8的特征值的開方
1 summary(w.pr, loadings=TRUE) 2 3 #### 作預測 4 predict(w.pr) 5 6 #### 畫碎石圖 7 screeplot(w.pr) 8 screeplot(w.pr,type="lines",pch=25,col=4) 9 10 biplot(w.pr) 11 12 princomp(~X1+X2+X3+X4+X5+X6+X7+X8,data=w, cor=T)
在這里我們也同樣可以看出一些關鍵的信息。
其中Standard devi....其中可以看出其標准差,即主成分的方差開方,也就是相應特征值lambda1、lambda2、lambda3到lambda8的特征值的開方。
proportion of Var....表示方差的貢獻率。
Cumulative proportion表示方差的累積貢獻率
Loadings=TURE表示列出載荷的內容,即有8個主成分,因此可以根據方差累積貢獻率大於85%可以選取主成分,以達到降維的目的。在這可以選取comp.1、comp.2、comp.2
第一主成分和第二主成分的散點圖
預測值:
(2)綜合得分:以個主成分的貢獻率為權,將其線性組合得到綜合評價函數:
C=(lambda1*C1+lambda2*C2+….+lambdaiCi)/(lambda1+lambda2+….+lambdai)
在這我們選取的是第一第二主成分,這其主成分得分為:w8
A<-c(1.7620762 ,1.7021873,0.96447683 )####默認為列
B<-sqrt(A)
B[1]
B[2]
B[3]
X1<-w.pr$scores[,1]
X1
X2<-w.pr$scores[,2]
X3<-w.pr$scores[,3]
X3
w2<-data.frame(X1,X2,X3)
w4<-matrix(X1,nrow=13)
w4
w41<-matrix(X2,nrow=13)
w41
w42<-matrix(X3,nrow=13)
w42
w5<-w4%*%B[1]/(B[1]+B[2]+B[3])
w5
w6<-w41%*%B[2]/(B[1]+B[2]+B[3])
w6
w7<-w42%*%B[3]/(B[1]+B[2]+B[3])
w7
w8=w5+w6+w7
w8
rank(w8)
或者這樣也可以調用mvstats包運用princomp.rank(m=3,w.pr)函數就可以
princomp.rank(m=3,w.pr)
(2)分類:
w=scale(w)####數據標准化
d=dist(w)####聚類距離
d
plot(h1,hang=-1,main=c("聚類分析"),sub=c("最短距離或者最近鄰法"))###使用labels=iris$Species[iris]加下標
hh1=rect.hclust(h1,k=3,border=c(2:3))###將h1分類加顏色
######表示最短距離或者最近鄰法
二、主成分回歸
在回歸分析中,曾經講過,當自變量出現多重共線性時,經典回歸方法作回歸系數的最小二乘法估計一般效果會比較差,而采取主成分回歸能夠客服直接回歸的不足,下面用例子進行說明。
對某地區某類消費品的銷售量y進行調查,它與下面四個變量有關,居民可支配收入(x1)、該類消費品平均價格指數(x2)、社會該消費品保有量(x3)、其他消費品平均價格指數(x4)。
用主成分回歸方法建立營銷量y與其他變量的關系。
數據如下:
x1=c(82.9,88.0,99.9,105.3,117.7,131.0,148.2,161.8,174.2,184.7)
x2=c(92.0,93.0,96.0,94.0,100,101,105,112,112,112)
x3=c(17.1,21.3,25.1,29.0,34,40,44,49,51,53)
x4=c(94,96,97,99,100,101,104,109,111,111)
y1=c(8.4,9.6,10.4,11.4,12.2,14.2,15.8,17.9,29.6,20.8)
用R編寫程序可得:
1 x1=c(82.9,88.0,99.9,105.3,117.7,131.0,148.2,161.8,174.2,184.7) 2 x2=c(92.0,93.0,96.0,94.0,100,101,105,112,112,112) 3 x3=c(17.1,21.3,25.1,29.0,34,40,44,49,51,53) 4 x4=c(94,96,97,99,100,101,104,109,111,111) 5 y1=c(8.4,9.6,10.4,11.4,12.2,14.2,15.8,17.9,29.6,20.8) 6 length(x3) 7 w=data.frame(x1,x2,x3,x4,y1) 8 lm.sol<-lm(y1~x1+x2+x3+x4) 9 summary(lm.sol)
運行結果如下:
在這可以看出,按4個變量得到的回歸方程為
y=-114.8655+0.02438x1-0.215x2-0.233x3+1.53x4
cor(w)
可以得到:相關系數矩陣
仔細分析得到的結果,可以看出計算出來的回歸系數並不顯著,可能存在多重共線性。
因此,用主成分回歸進行分析,其代碼如下:
1 conomy<-princomp(~x1+x2+x3+x4,cor=T)####主成分分析 2 summary(conomy,loadings=T)
得到結果如下:
首先我們可以看出前兩個主成分的貢獻率以達到99%,在這我們可以看到第一主成分是關於居民可支配收入(x1)、該類消費品平均價格指數(x2)、社會該消費品保有量(x3)、其他消費品平均價格指數(x4)。在這我們可以看到第二主成分也是如此。
但是我們將兩個主成分帶入原公式,這是不符合經濟意義的,所以我們的多重共線性並未消。
下面作主成分回歸
首先計算主成分的預測值,並且將第一主成分和第二主成分的預測值放入數據框中,然后做主成分回歸。
pre<-predict(conomy)
pre
w$z1<-pre[,1]
w$z2<-pre[,2]
lm.sol2<-lm(y1~z1+z2,data=w)
summary(lm.sol2)
可以得到回歸方程為:y1=15.03-2.738z1+4.4458z2
這樣看不出與變量之間的關系於是我們變換一下表達式
ba<-coef(lm.sol2)###提取回歸系數
ba[2]
A<-loadings(conomy)
x.bar<-conomy$center####s數據中心,均值
x.sd<-conomy$scale###標准差
coef<-(ba[2]*A[,1]+ba[3]*A[,2])/x.sd###x1x2x3x4的系數
coef
ba0<-ba[1]-sum(x.bar*coef)###常數
ba0
即回歸方程為:y1= -641.2539+0.01315832x1+0.49188577x2 -0.14875597 x3+ 0.51413117x4
可見這樣較為合理。