倔強(qiáng)的演奏者

Principal component analysis (PCA) 是一個(gè)統(tǒng)計(jì)學(xué)方法，用一組較少的不相關(guān)的變量代替大量相關(guān)變量腮鞍，同時(shí)盡可能保留初始變量的信息肮帐，這些推導(dǎo)所得的變量成為主成分甜奄。

——《R語(yǔ)言實(shí)戰(zhàn)》

介紹

主成分分析用來(lái)從多變量數(shù)據(jù)里面提取最重要的信息吆寨，一組數(shù)據(jù)的信息對(duì)應(yīng)著其總方差，所以PCA的目的就是使用一組較少不相關(guān)的變量代替大量相關(guān)變量拧廊，用principal components（下面用主成分來(lái)代指）來(lái)表示监徘，這些新變量對(duì)應(yīng)原數(shù)據(jù)的線性結(jié)合，新變量的數(shù)目少于或等于原變量數(shù)目吧碾，其中第一主成分對(duì)初始變量集的方差解釋性最大凰盔，隨后的每一個(gè)主成分都最大化它對(duì)方差的解釋程度，同時(shí)與之前所有的主成分都正交倦春。

注意：原理上理解PCA需要線性代數(shù)知識(shí)户敬，我只是學(xué)習(xí)使用R來(lái)處理數(shù)據(jù)（況且我也不懂），所以想要打破砂鍋問(wèn)到底的讀者可以查看WIKI: PCA睁本，說(shuō)的很清楚尿庐。

圖示理解

二維數(shù)據(jù)理解

image

上圖Plot A，我們可以看到這是一個(gè)簡(jiǎn)單的X-Y 坐標(biāo)系展示了一個(gè)二維數(shù)據(jù)呢堰，降維的方式在這里被展示為紅綠兩個(gè)主成分抄瑟，當(dāng)然PCA是假設(shè)對(duì)數(shù)據(jù)總方差解釋最大的成分是最重要的，也就是所謂的第一主成分枉疼。

image

圖中PC坐標(biāo)軸的PC1 就是所謂的第一主成分皮假，解釋了示例數(shù)據(jù)集最大的方差，PC2則是對(duì)方差解釋性排第二的成分骂维，同時(shí)與PC1正交惹资。這里就完成了對(duì)目標(biāo)數(shù)據(jù)的降維（等于是咱們就需要PC1就可以了，PC2對(duì)數(shù)據(jù)的解釋性不是占很大比例）航闺。

從技術(shù)上來(lái)講褪测，每個(gè)成分保留的變量比例是通過(guò)一個(gè)的叫做特征值(Eigenvalue)來(lái)計(jì)算的。因?yàn)楹芏鄷r(shí)候PCA很適合用在擁有多組相關(guān)變量的數(shù)據(jù)里面，因?yàn)橄嚓P(guān)就意味著這些變量的信息有時(shí)候是冗余的汰扭，同時(shí)PCA就可以用很少的不相關(guān)的成分來(lái)對(duì)整個(gè)數(shù)據(jù)進(jìn)行解釋?zhuān)?/p>

image

image

綜上，PCA的目的是：

• 確定數(shù)據(jù)里面的隱藏的信息
• 降維的同時(shí)去除數(shù)據(jù)噪音
• 確定相關(guān)變量

前期準(zhǔn)備

相關(guān)的R包：

• prcomp() and princomp() [built-in R stats package]
• PCA() [FactoMineR package]
• dudi.pca() [ade4 package]
• epPCA() [ExPosition package]
• principal())[psych package]

不管你選用哪些包福铅，你都需要做兩件事：

1. 數(shù)據(jù)預(yù)處理然后提取主成分
2. 可視化

本篇文章使用FactoMineR（負(fù)責(zé)分析）和factoextra（負(fù)責(zé)可視化）萝毛，首先是安裝和載入包：

> install.packages(c("FactoMineR", "factoextra"))
> library("FactoMineR")
Warning message:
程輯包‘FactoMineR’是用R版本3.5.1 來(lái)建造的 
> library("factoextra")
Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

示例數(shù)據(jù)理解

這里使用來(lái)自factoextra包的decathlon2數(shù)據(jù)集：

image

> head(decathlon2)
          X100m Long.jump Shot.put High.jump X400m X110m.hurdle Discus Pole.vault Javeline
SEBRLE    11.04      7.58    14.83      2.07 49.81        14.69  43.75       5.02    63.19
CLAY      10.76      7.40    14.26      1.86 49.37        14.05  50.72       4.92    60.15
BERNARD   11.02      7.23    14.25      1.92 48.93        14.99  40.87       5.32    62.77
YURKOV    11.34      7.09    15.19      2.10 50.42        15.31  46.26       4.72    63.44
ZSIVOCZKY 11.13      7.30    13.48      2.01 48.62        14.17  45.67       4.42    55.37
McMULLEN  10.83      7.31    13.76      2.13 49.91        14.38  44.41       4.42    56.37
          X1500m Rank Points Competition
SEBRLE     291.7    1   8217    Decastar
CLAY       301.5    2   8122    Decastar
BERNARD    280.1    4   8067    Decastar
YURKOV     276.4    5   8036    Decastar
ZSIVOCZKY  268.0    7   8004    Decastar
McMULLEN   285.1    8   7995    Decastar

這些是用來(lái)描述運(yùn)動(dòng)員的表現(xiàn)的數(shù)據(jù)，我們可以看到有27行和13列：

• 淡藍(lán)色（Active individuals）是可用分析樣本
• 淡粉色（Active variables）是可用分析變量
• 深藍(lán)色（Supplementary individuals）是測(cè)試數(shù)據(jù)
• 最后三列則是檢驗(yàn)變量（說(shuō)白了就是測(cè)試你PCA分析做的準(zhǔn)不準(zhǔn)的）

提取分析數(shù)據(jù)集

decathlon2.active <- decathlon2[1:23, 1:10]

數(shù)據(jù)標(biāo)準(zhǔn)化

在我們接觸到的數(shù)據(jù)里面滑黔，變量擁有著不同的測(cè)量尺度的（例如單位是千米/千克/厘米等）笆包，所以我們需要對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，不然做出來(lái)的PCA分析是很不準(zhǔn)確的略荡。

最常用的就是Z-score的標(biāo)準(zhǔn)化方法了庵佣，適用于大多數(shù)數(shù)據(jù)（原理請(qǐng)自行了解），它將數(shù)據(jù)轉(zhuǎn)化為擁有方法為1汛兜，平均值為0的標(biāo)準(zhǔn)化數(shù)據(jù)巴粪。

R語(yǔ)言的基礎(chǔ)函數(shù)sacle()可以完成這件事。

而R包FactoMineR的PCA()函數(shù)是自動(dòng)會(huì)幫你的標(biāo)準(zhǔn)化的粥谬，當(dāng)然你也可以選擇自己標(biāo)準(zhǔn)化肛根。

PCA分析

res.pca <- PCA(X = decathlon2.active, scale.unit = TRUE, ncp = 5, graph = FALSE)

解釋一下PCA()函數(shù)的可選參數(shù)：

• X: a data frame. Rows are individuals and columns are numeric variables
• scale.unit: a logical value. If TRUE, the data are scaled to unit variance before the analysis. This standardization to the same scale avoids some variables to become dominant just because of their large measurement units. It makes variable comparable.
• ncp: number of dimensions kept in the final results.
• graph: a logical value. If TRUE a graph is displayed.

我們來(lái)看一下最后的分析結(jié)果：

> print(res.pca)
**Results for the Principal Component Analysis (PCA)**
The analysis was performed on 23 individuals, described by 10 variables
*The results are available in the following objects:

   name               description                          
1  "$eig"             "eigenvalues"                        
2  "$var"             "results for the variables"          
3  "$var$coord"       "coord. for the variables"           
4  "$var$cor"         "correlations variables - dimensions"
5  "$var$cos2"        "cos2 for the variables"             
6  "$var$contrib"     "contributions of the variables"     
7  "$ind"             "results for the individuals"        
8  "$ind$coord"       "coord. for the individuals"         
9  "$ind$cos2"        "cos2 for the individuals"           
10 "$ind$contrib"     "contributions of the individuals"   
11 "$call"            "summary statistics"                 
12 "$call$centre"     "mean of the variables"              
13 "$call$ecart.type" "standard error of the variables"    
14 "$call$row.w"      "weights for the individuals"        
15 "$call$col.w"      "weights for the variables"         

可視化和解讀

這里需要用到factoextra包來(lái)做這件事，無(wú)論你用的是上面哪個(gè)包做的PCA分析漏策，都可以簡(jiǎn)單快捷的對(duì)數(shù)據(jù)進(jìn)行提取并可視化派哲，參數(shù)解釋如下：

• get_eigenvalue(res.pca): Extract the eigenvalues/variances of principal components
• fviz_eig(res.pca): Visualize the eigenvalues
• get_pca_ind(res.pca), get_pca_var(res.pca): Extract the results for individuals and variables, respectively.
• fviz_pca_ind(res.pca), fviz_pca_var(res.pca): Visualize the results individuals and variables, respectively.
• fviz_pca_biplot(res.pca): Make a biplot of individuals and variables.

特征值/方差解釋性 Eigenvalues / Variances

就像上面說(shuō)的那樣，特征值可以衡量各個(gè)成分保留的數(shù)據(jù)方差的多少掺喻，從第一主成分開(kāi)始依次降低芭届，我們通過(guò)觀察特征值來(lái)決定我們需要考慮保留幾個(gè)成分。

> library("factoextra")
> eig.val <- get_eigenvalue(res.pca)
> eig.val
       eigenvalue variance.percent cumulative.variance.percent
Dim.1   4.1242133        41.242133                    41.24213
Dim.2   1.8385309        18.385309                    59.62744
Dim.3   1.2391403        12.391403                    72.01885
Dim.4   0.8194402         8.194402                    80.21325
Dim.5   0.7015528         7.015528                    87.22878
Dim.6   0.4228828         4.228828                    91.45760
Dim.7   0.3025817         3.025817                    94.48342
Dim.8   0.2744700         2.744700                    97.22812
Dim.9   0.1552169         1.552169                    98.78029
Dim.10  0.1219710         1.219710                   100.00000

eigenvalue的值cumulative之后是10感耙，對(duì)反差的解釋程度在第二列褂乍，例如此處的第一主成分解釋了數(shù)據(jù)41%的方差，最后一列表明了前幾個(gè)主成分解釋了多少的方差即硼，那我們可以看到一二主成分就解釋了大約60%的差異树叽，而在eigenvalue的選擇上：

• >1的eigenvalue說(shuō)明了這個(gè)主成分不止解釋了一個(gè)原始變量。
• 同時(shí)你還可以限制解釋比例來(lái)獲取主成分谦絮，例如cumulative以后超過(guò)70%的题诵。

目前還沒(méi)有一個(gè)統(tǒng)一的標(biāo)準(zhǔn)來(lái)定義如何去選去主成分?jǐn)?shù)目，所以這很大程度上依賴(lài)于你對(duì)數(shù)據(jù)分析的目的和你對(duì)數(shù)據(jù)的認(rèn)知程度层皱。

在這個(gè)示例數(shù)據(jù)里面性锭，前三個(gè)主成分解釋了72%的方差差異性，這就是一個(gè)可以接受的比例了叫胖。

在R語(yǔ)言實(shí)戰(zhàn)里面草冈，psych包可以畫(huà)一個(gè)碎石圖，根據(jù)碎石檢驗(yàn)?zāi)玫轿覀冃枰Ａ舻財(cái)?shù)據(jù)，也不失為以一個(gè)可取的方法怎棱。

當(dāng)然factoextra包也是可以做這個(gè)事情的fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50))：

image

對(duì)變量作圖

factoextra包自帶了提取變量的分析結(jié)果→get_pca_var

> var <- get_pca_var(res.pca)
> var
Principal Component Analysis Results for variables
 ===================================================
  Name       Description                                    
1 "$coord"   "Coordinates for the variables"                
2 "$cor"     "Correlations between variables and dimensions"
3 "$cos2"    "Cos2 for the variables"                       
4 "$contrib" "contributions of the variables"  

可以看到相關(guān)的描述哩俭，其中var.cos2 = var.coord * var.coord，而根據(jù)這些數(shù)據(jù)我們就可以進(jìn)行作圖：

> # Coordinates
> head(var$coord)
                  Dim.1       Dim.2      Dim.3       Dim.4      Dim.5
X100m        -0.8506257 -0.17939806  0.3015564  0.03357320 -0.1944440
Long.jump     0.7941806  0.28085695 -0.1905465 -0.11538956  0.2331567
Shot.put      0.7339127  0.08540412  0.5175978  0.12846837 -0.2488129
High.jump     0.6100840 -0.46521415  0.3300852  0.14455012  0.4027002
X400m        -0.7016034  0.29017826  0.2835329  0.43082552  0.1039085
X110m.hurdle -0.7641252 -0.02474081  0.4488873 -0.01689589  0.2242200
> # Cos2: quality on the factore map
> head(var$cos2)
                 Dim.1        Dim.2      Dim.3        Dim.4      Dim.5
X100m        0.7235641 0.0321836641 0.09093628 0.0011271597 0.03780845
Long.jump    0.6307229 0.0788806285 0.03630798 0.0133147506 0.05436203
Shot.put     0.5386279 0.0072938636 0.26790749 0.0165041211 0.06190783
High.jump    0.3722025 0.2164242070 0.10895622 0.0208947375 0.16216747
X400m        0.4922473 0.0842034209 0.08039091 0.1856106269 0.01079698
X110m.hurdle 0.5838873 0.0006121077 0.20149984 0.0002854712 0.05027463
> # Contributions to the principal components
> head(var$contrib)
                 Dim.1      Dim.2     Dim.3       Dim.4     Dim.5
X100m        17.544293  1.7505098  7.338659  0.13755240  5.389252
Long.jump    15.293168  4.2904162  2.930094  1.62485936  7.748815
Shot.put     13.060137  0.3967224 21.620432  2.01407269  8.824401
High.jump     9.024811 11.7715838  8.792888  2.54987951 23.115504
X400m        11.935544  4.5799296  6.487636 22.65090599  1.539012
X110m.hurdle 14.157544  0.0332933 16.261261  0.03483735  7.166193

相關(guān)曲線作圖 Correlation circle

一個(gè)變量和一個(gè)主成分之間的關(guān)系的代表著在PC坐標(biāo)系里面該變量的坐標(biāo)拳恋，對(duì)變量作圖我們可以用fviz_pca_var函數(shù)：

fviz_pca_var(res.pca, col.var = "black")

image

這張圖也可以稱(chēng)為變量相關(guān)圖凡资，它展示了變量組內(nèi)包括和主成分之間的關(guān)系，正相關(guān)的變量是彼此靠近的谬运，負(fù)相關(guān)的變量師南轅北轍的隙赁，而從中心點(diǎn)到變量的長(zhǎng)度則代表著變量在這個(gè)維度所占的比例（也可以理解為質(zhì)量，quality）梆暖。

代表質(zhì)量作圖 Quality of representation

變量在PCA結(jié)果里面的質(zhì)量（quality）稱(chēng)為cos2 (square cosine, squared coordinates) 伞访，可以用corrplot包進(jìn)行可視化：

library("corrplot")
corrplot(var$cos2, is.corr=FALSE)

image

或者使用factoextra里面的fviz_cos()函數(shù)也可以搞定，只不過(guò)沒(méi)有上圖好看而已：

# Total cos2 of variables on Dim.1 and Dim.2
fviz_cos2(res.pca, choice = "var", axes = 1:2)

image

需要注意的是：

• 一個(gè)較高的cos2值代表著這個(gè)變量對(duì)該主成分有較大的貢獻(xiàn)值轰驳，在這個(gè)示例數(shù)據(jù)里面就表現(xiàn)在相關(guān)曲線作圖里面的約靠近圓的邊緣厚掷。
• 一個(gè)較低的cos2值代表著這個(gè)變量并沒(méi)有很好的被主成分所代表，在相關(guān)曲線作圖里面就約靠近圓心级解。
• cos2值就是為了衡量一個(gè)變量的有用程度蝗肪。
• 越高就代表著這個(gè)變量在主成分分析里面約重要。
• 第一張圖每一行加在一起等于1蠕趁。

由于該包的作圖是基于ggplot2的薛闪，所以在作圖的時(shí)候也可以選擇的color = 'cos2'來(lái)進(jìn)行上色：

# Color by cos2 values: quality on the factor map
fviz_pca_var(res.pca, col.var = "cos2",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
             repel = TRUE # Avoid text overlapping
             )

image

當(dāng)然，我們也可以通過(guò)改變變量的透明度來(lái)說(shuō)明其重要性：

# Change the transparency by cos2 values
fviz_pca_var(res.pca, alpha.var = "cos2")

image

對(duì)主成分的貢獻(xiàn)作圖 Contributions of variables to PCs

每個(gè)變量對(duì)主成分的貢獻(xiàn)比例是不一樣的俺陋，而對(duì)第一二主成分貢獻(xiàn)最大的變量是最重要的豁延，而對(duì)主成分貢獻(xiàn)很低的變量甚至可以移出數(shù)據(jù)的簡(jiǎn)化分析，可以通過(guò)類(lèi)似質(zhì)量作圖的模式來(lái)進(jìn)行可視化：

library("corrplot")
corrplot(var$contrib, is.corr=FALSE) 

image

factorextra包里面的fviz_contrib函數(shù)也可以用于可視化：

# Contributions of variables to PC1
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10)
# Contributions of variables to PC2
fviz_contrib(res.pca, choice = "var", axes = 2, top = 10)

image

image

其中axes指示坐標(biāo)腊状，top指示畫(huà)前多少個(gè)變量诱咏，當(dāng)然也可以畫(huà)在一張圖上：

fviz_contrib(res.pca, choice = "var", axes = 1:2, top = 10)

image

紅色虛線代表著平均貢獻(xiàn)，高于平均值的可以被認(rèn)為算是重要變量缴挖。

更新截至日期袋狞，2018-10-27.
待續(xù)。