R語言數據分析與挖掘(第六章):主成分分析(2)——案例講解

  • 2019 年 12 月 16 日
  • 筆記

這一講通過一個案例講解主成分分析。

在R中用於完成主成分分析的函數是princomp()，該函數有2種調用方式：

1.公式形式

基本語法為：

princomp(formula, data = NULL, subset, na.action, ...)

參數介紹：

formula：指定用於主成分的公示對象，類似於回歸分析和方差分析中的公式對象；

data：指定用於主成分分析的數據對象，一般為數據框；

subset：指定可選的向量，表示選擇的樣本子集；

na.action：一個函數，指定缺失數據的處理方法，若為NULL，則使用函數na.omit()刪除缺失數據。

2. 矩陣形式

基本語法為：

princomp(x, cor = FALSE, scores = TRUE, covmat = NULL,           subset = rep_len(TRUE, nrow(as.matrix(x))), fix_sign = TRUE, ...)

參數介紹：

X:指定用於主成分分析的數據對象，一般為數值矩陣的數據框:

Cor:邏輯值，指定用於主成分分析中採用的矩陣形式(相關矩陣或協方差陣),為TRUE表示用樣本的相關矩陣做主成分分析，為TALSE表示用樣本的協方差陣做主成分分析,默認值為FALSE;

Scores:邏輯值，指定是否計算各主成分的分量，即是否計算樣本的主成分得分，默認值TRUE;

Covmati:指定協方差陣，或者為cov.wt()提供的協方差列表，當數據不用參數x提供時，可由協方差陣提供;

Subset:指定可選向量，表示選擇的樣本子集。

函數princomp()的返回值為-個列表，包括:

sdev表示各主成分的標準差;

loadings表示載荷矩陣;

center表示各指標的樣本均值;

scale表示各指標的樣本標準差;

n.obs表示觀測樣本的個數；

scores表示主成分得分(當參數scores= TRUE時提供)。

此外，也可以利用其他函數來提取主成分分析的結果:

　　函數summary()可用於提取主成分的資訊;

　　函數loadings()可用於提取載荷矩陣;

　　函數predict()可用於計算主成分得分;

　　函數screeplot()可用於繪製主成分的碎石圖;

　　函數biplot()可用於繪製數據關於主成分的散點圖和原坐標在主成分下的方向。

綜合案例:longley數據集的變數降維及回歸

下面利用第5章提到的longley數據集進行操作演練，該數據集收集了1947年至1962年16年的宏觀經濟數據，包含7個變數，分別為物價指數平減後的國民生產總值(GNP.deflator)、 國民生產總值(GNP)、(Amed.Frees)人口總量(Popultin)、年份(Year) 和就業人數(Emplyed)。失業人數(Unemployed)、 軍隊人數（Amed.Forces），人口總(Population）、年份(Year)和就業人數（Employed）。

> data(longley)  > summary(longley)    GNP.deflator         GNP          Unemployed     Armed.Forces   Min.   : 83.00   Min.   :234.3   Min.   :187.0   Min.   :145.6   1st Qu.: 94.53   1st Qu.:317.9   1st Qu.:234.8   1st Qu.:229.8   Median :100.60   Median :381.4   Median :314.4   Median :271.8   Mean   :101.68   Mean   :387.7   Mean   :319.3   Mean   :260.7   3rd Qu.:111.25   3rd Qu.:454.1   3rd Qu.:384.2   3rd Qu.:306.1   Max.   :116.90   Max.   :554.9   Max.   :480.6   Max.   :359.4     Population         Year         Employed   Min.   :107.6   Min.   :1947   Min.   :60.17   1st Qu.:111.8   1st Qu.:1951   1st Qu.:62.71   Median :116.8   Median :1954   Median :65.50   Mean   :117.4   Mean   :1954   Mean   :65.32   3rd Qu.:122.3   3rd Qu.:1958   3rd Qu.:68.29   Max.   :130.1   Max.   :1962   Max.   :70.55  > corr<-cor(longley)  > library(corrplot)  corrplot 0.84 loaded  > corrplot(corr)

該數據的相關矩陣如上圖，可以很直觀的看出，數據存在多重共線性問題，下面計算方差膨脹因子(VIF)來判斷個變數的共線性問題的嚴重性：

> library(car)  > vif(lm(Employed~.,data=longley))  GNP.deflator          GNP   Unemployed Armed.Forces   Population     135.53244   1788.51348     33.61889      3.58893    399.15102          Year     758.98060

一般認為：

當0<VIF<10，不存在多重共線性；

當10<=VIF<100,存在較強的多重共線性；

當VIF>= 100，多重共線性非常嚴重。

上訴輸出結果顯示，變數GNP.deflator 、GNP、Population和Year存在嚴重的多重共線性。 為了解決多重共線性問題，下面對數據集進行主成分分析，去掉響應變數後進行分析：

> (pr1<-princomp(longley[,-7], cor = TRUE))  Call:  princomp(x = longley[, -7], cor = TRUE)    Standard deviations:      Comp.1     Comp.2     Comp.3     Comp.4     Comp.5     Comp.6  2.14554820 1.08413122 0.45102702 0.12218125 0.05051797 0.01940897     6  variables and  16 observations.  > summary(pr1)  Importance of components:                            Comp.1    Comp.2     Comp.3      Comp.4  Standard deviation     2.1455482 1.0841312 0.45102702 0.122181253  Proportion of Variance 0.7672295 0.1958901 0.03390423 0.002488043  Cumulative Proportion  0.7672295 0.9631196 0.99702383 0.999511871                               Comp.5       Comp.6  Standard deviation     0.0505179747 1.940897e-02  Proportion of Variance 0.0004253443 6.278469e-05  Cumulative Proportion  0.9999372153 1.000000e+00

輸出結果中：Standard deviation表示主成分的標準差，即特徵值的開方；Proportion of Variance表示每個主成分的貢獻率；Cumulative Proportion表示主成分的累計貢獻率，通常情況下，參考該參數的值進行主成分的選擇，一般累計貢獻率超過85%的主成分，在本案例中，前2個主成分的累計貢獻率已經達到96%，即前2個主成分能解釋原始變數的96%的資訊，故應該選擇前2個主成分，剔除後面4個主成分，達到降維的目的。

下面利用函數loadings()輸出載荷矩陣：

> loadings(pr1)    Loadings:               Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6  GNP.deflator  0.462         0.149  0.793  0.338  0.135  GNP           0.462         0.278 -0.122 -0.150 -0.818  Unemployed    0.321 -0.596 -0.728               -0.107  Armed.Forces  0.202  0.798 -0.562  Population    0.462         0.196 -0.590  0.549  0.312  Year          0.465         0.128        -0.750  0.450                   Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6  SS loadings     1.000  1.000  1.000  1.000  1.000  1.000  Proportion Var  0.167  0.167  0.167  0.167  0.167  0.167  Cumulative Var  0.167  0.333  0.500  0.667  0.833  1.000

可以得到|成分!原的變數的線性關係:

第1主成分：

F1=0.462*CNPdnator-0.462*GNP-0.321*Unemploycd-0.202*Amed.Forces-0.462*Ppulation-0.465*Year

第2主成分

Fn=-0.596*Unemployed+0.798* Armed.Forces

下面對選樣的主成分進行命名，主成分命名的原則應注意以下幾點:

由於F1中包含不同變數的資訊，在進行命名時需費綜合考慮不同變數的資訊；需要注意線性關係中係數的符號和絕對值大小。

基於上述原則，第一主成分對應係數的符號相同，均為負號，除變最Unemployed和Amed.Forces以外，其他原始變數的係數絕對值均為0.5左有，故第1主成分反映了經濟的蕭條程度，經濟越蕭條，F1的取值越大，各原始變數的取值越小，反之，經濟越不蕭條，F1的取值越小，各原始變數的取值越大。

第二主成分要與原始變數Unemployed和Amed.Forces有關,且變Unemployed對應的係數為負，Amed.Forces對應的係數為正，係數的絕對值均大於0.5,故第二主成分反映了軍備程度，軍備程度越高，F2的取值越大，失業人數越小，且軍隊人數越多。反之，軍備程度越低，F2的取值越小，失業人數越多，軍隊人數越少。

　　儘管上述分析對兩個主成分進行了命名，但命名結果還是不能很好地解釋原版始變數的資訊，因為主成分命名本身就具有較人的難度，我們能做到的就是盡量個面地涵蓋原始變數的資訊。

下面是利用函數predict()輸出主成分的得分。

> predict(pr1)            Comp.1     Comp.2      Comp.3       Comp.4       Comp.5  1947 -3.59294203 -0.7761110  0.31804507 -0.169624216  0.009085721  1948 -3.10924187 -0.8768853  0.66329350  0.130051755  0.063565251  1949 -2.42014988 -1.5905016 -0.50961596 -0.009106066  0.005934717  1950 -2.16257276 -1.3181781 -0.11494311 -0.063265342 -0.063873551  1951 -1.48540767  1.2763227 -0.03004947  0.100660376  0.053969833  1952 -1.04339031  1.9851409 -0.16650384  0.047914580  0.038252420  1953 -0.72546382  1.9732212  0.06933769 -0.073217347  0.022425386  1954  0.03355589  0.6124972 -1.07307182 -0.067293112 -0.002331281  1955  0.10277563  0.7162362 -0.10076710 -0.104426569 -0.102049470  1956  0.46417097  0.5658113  0.30255566  0.018137849 -0.086508832  1957  0.98638447  0.4435320  0.45984009  0.123243732 -0.024470311  1958  1.87669233 -0.8914800 -0.69963411  0.193195372  0.022381351  1959  2.00361238 -0.3992485  0.27468475  0.148640632 -0.037890240  1960  2.43855615 -0.5154723  0.37766452  0.063619308 -0.016767727  1961  3.17896368 -1.0224220 -0.20859150 -0.070338720  0.058280241  1962  3.45445683 -0.1824626  0.43775561 -0.268192230  0.059996493             Comp.6  1947  0.002663801  1948  0.012373058  1949  0.005226048  1950 -0.014125087  1951 -0.044081763  1952 -0.013169060  1953  0.032810005  1954  0.017241531  1955 -0.019544921  1956  0.014604365  1957  0.028044545  1958  0.008369752  1959 -0.024302116  1960  0.004500429  1961 -0.001372125  1962 -0.009238462  > sort(predict(pr1)[,1],decreasing = T)         1962        1961        1960        1959        1958        1957   3.45445683  3.17896368  2.43855615  2.00361238  1.87669233  0.98638447         1956        1955        1954        1953        1952        1951   0.46417097  0.10277563  0.03355589 -0.72546382 -1.04339031 -1.48540767         1950        1949        1948        1947  -2.16257276 -2.42014988 -3.10924187 -3.59294203  > sort(predict(pr1)[,2],decreasing = T)        1952       1953       1951       1955       1954       1956   1.9851409  1.9732212  1.2763227  0.7162362  0.6124972  0.5658113        1957       1962       1959       1960       1947       1948   0.4435320 -0.1824626 -0.3992485 -0.5154723 -0.7761110 -0.8768853        1958       1961       1950       1949  -0.8914800 -1.0224220 -1.3181781 -1.5905016

從第一主成分看，1947、1948和1949年的得分較高，說明那幾年的經濟較為蕭條，得分較低的是1960、1961和1962年，說明那幾年的經濟形勢較好。

從第二主成分來看，1952、1953和1951年得分較高，說明那幾年的軍備情況較好，1961、1950和1949年的得分較低，說明那幾年的軍備情況較差。

下面利用函數是繪製主成分碎石圖：

screeplot(pr1,type="lines",main="")

繪圖形勢選擇折線圖，碎石圖展現的是各主成分沒有解釋道原始數據的資訊情況，折線圖越靠近x軸，說明主成分解釋原始變數的資訊越多。可以看出選擇前2個主成分是合理的。

下面利用函數biplot()繪製雙坐標圖，默認情況下，該函數繪製第一、二主成分樣本散點圖和原始坐標在第一、二主成分下的方向。

biplot(pr1)

主成分回歸 主成分回歸可以參考前面講的回歸分析看程式碼：

pre1<-predict(pr1)  longley$F1<-pre1[,1]  longley$F2<-pre1[,2]  longley$F3<-pre1[,3]  lm.pr1<-lm(Employed~F1+F2+F3,data=longley)  summary(lm.pr1)  beta<-coef(lm.pr1)  A<-loadings(pr1)[,1:3]  x.mean<-pr1$center  x.sd<-pr1$scale  coef<-A%*%beta[2:4]/x.sd  beta0<-beta[1]-x.mean%*%coef  c(beta0,coef)

VirMach 便宜 VPS

QNews

QNews