本篇目錄

主成份分析
選擇多少個主成份？
- 陡坡圖(Scree plot)-凱莎原則
- 累積解釋圖(Pareto plot)
線性組合的係數 (主成份和原變數的關係)
總結
Reference
R and packages version

主成份分析(Principal Component Analysis)

這裡拿網路上一個公開資料，2012年美國職棒MLB的資料，來進行分析，資料載點如下。

在本篇文章中，會使用到以下函式：

prcomp()：主成份分析的基本函式
plot()：繪製陡坡圖(screet plot)，選擇多少個主成份
dotchart()：繪製主成份負荷圖(PCA loadings plot)
biplot()：繪製主成份負荷圖(PCA loadings plot)

1. 主成份分析

當下載好資料後，第一步便是先讀取資料：

data <- read.csv("2012MLB.csv",  # 資料檔名 
                 header=T,          # 資料中的第一列，作為欄位名稱
                 sep=",")           # 將逗號視為分隔符號來讀取資料

head(data)

##                   Team   G   R    H H1B H2B H3B  HR RBI  BB   SO  SB   AVG
## 1        Texas Rangers 152 764 1444 941 283  32 188 738 448 1030  89 0.275
## 2   Los Angeles Angels 153 726 1424 972 252  20 180 691 423 1041 125 0.273
## 3     Colorado Rockies 152 718 1425 932 286  49 158 680 432 1129  96 0.272
## 4  St. Louis Cardinals 153 723 1447 983 277  35 152 691 503 1128  89 0.272
## 5 San Francisco Giants 153 683 1419 997 273  54  95 642 453 1032 111 0.270
## 6       Detroit Tigers 152 689 1379 922 267  37 153 663 490 1042  53 0.268
##     OBP
## 1 0.335
## 2 0.331
## 3 0.329
## 4 0.338
## 5 0.327
## 6 0.336

在這裡，選擇一壘安打、二壘安打、三壘安打、全壘打、打點、盜壘次數、四壞球，這七個變數，進行主成份分析prcomp()：

pca <- prcomp(formula = ~ H1B+H2B+H3B+HR+RBI+SB+BB,  #選擇七個變數 
              data = data,                           # 資料
              scale = TRUE)                          # 正規化資料
# 這就是我們的主成份
pca

## Standard deviations (1, .., p=7):
## [1] 1.4222856 1.3785035 1.0108522 0.9578441 0.7700729 0.7131148 0.1897347
## 
## Rotation (n x k) = (7 x 7):
##             PC1        PC2         PC3          PC4         PC5
## H1B -0.40991503 -0.4681242  0.07174689  0.056704066 -0.07882016
## H2B -0.51441491 -0.2004156  0.01669591  0.255448162 -0.46809834
## H3B  0.01853759 -0.5595940 -0.19427151 -0.004051477  0.71490431
## HR  -0.34336124  0.5417488 -0.03416307 -0.394140194  0.26396281
## RBI -0.64629912  0.1016251 -0.25396353 -0.156840299  0.20084751
## SB   0.16722000 -0.2741655 -0.52853255 -0.679207860 -0.39181790
## BB   0.05866272  0.2203673 -0.78218985  0.538744635 -0.00676713
##             PC6         PC7
## H1B -0.66701643 -0.39159028
## H2B  0.62315846 -0.14911690
## H3B  0.34921449 -0.12535585
## HR   0.10991843 -0.59189466
## RBI -0.12239863  0.65387482
## SB   0.04037564 -0.03239357
## BB  -0.12695740 -0.17252886

上面的報表是這樣解釋：

Standard deviations：特徵值開根號
Rotation：特徵向量，也就是各個主成份，所對應的線性組合(linear combination)的係數

2. 選擇多少個主成份？

當主成份算出來以後，接下來要做的是「選擇幾個主成份」！

我們可以繪製「陡坡圖Scree plot」以及「累積解釋圖Pareto plot」：

陡坡圖(Scree plot)-凱莎原則

# 使用plot()函式
plot(pca,         # 放pca
     type="line", # 用直線連結每個點
     main="Scree Plot for 2012MLB") # 主標題

# 用藍線標示出特徵值=1的地方
abline(h=1, col="blue") # Kaiser eigenvalue-greater-than-one rule

根據凱莎原則，特徵值大於1的主成份就可以選取；而且第三個以後的主成份變異趨於平緩，因此選擇前三個主成份是比較好的選擇。

累積解釋圖(Pareto plot)

這裡就比較複雜，需要進行四個步驟：

求出每個主成份的特徵值(也就是variance = std^2)

    vars <- (pca$sdev)^2  # 從pca中取出標準差(pca$sdev)後再平方，計算variance(特徵值)
    vars

## [1] 2.02289644 1.90027181 1.02182222 0.91746533 0.59301228 0.50853268
## [7] 0.03599925

計算每個主成分的解釋比例 = 各個主成份的特徵值/總特徵值

    # 計算每個主成分的解釋比例 = 各個主成分的特徵值/總特徵值
    props <- vars / sum(vars)    
    props

## [1] 0.288985205 0.271467401 0.145974603 0.131066475 0.084716040 0.072647526
## [7] 0.005142749

累加每個主成份的解釋比例(aggregated effects)

    cumulative.props <- cumsum(props)  # 累加前n個元素的值
    cumulative.props

## [1] 0.2889852 0.5604526 0.7064272 0.8374937 0.9222097 0.9948573 1.0000000

把累積解釋比例畫成圖：

    #當我們取前三個主成份，可以解釋 70.64% 的變異
    cumulative.props[3]

## [1] 0.7064272

    # 累積解釋比例圖
    plot(cumulative.props)

所以原本的資料集，經過主成份分析後，會轉換成新的以主成份代替的資料集(pca$x)。
以下步驟是取前三個主成份，作為新的資料集：

# pca$rotation 
top3_pca.data <- pca$x[, 1:3]
top3_pca.data

##            PC1         PC2         PC3
## 1  -2.65536140  0.04641055  0.05124254
## 2  -1.37712847  0.01360254  0.24495540
## 3  -1.57754875 -1.72554295  0.14807629
## 4  -1.76751032 -0.84074064 -0.75258974
## 5  -0.44097214 -3.66454431 -0.36109275
## 6  -1.08166263 -0.10861369  0.27429500
## 7  -0.45474791 -2.65730709  1.13595923
## 8  -2.46449798  0.03093137  1.32250201
## 9  -0.11216989 -1.29948488 -0.83784413
## 10 -1.01441489  0.55336109  0.51565842
## 11 -1.53519975  3.14421085 -1.02043498
## 12 -1.20736887 -0.28285945 -1.24157398
## 13 -1.06901074  0.22513361 -0.72630319
## 14  0.01063384 -0.25128338  0.55678238
## 15 -0.25896225  1.04428266  0.46364492
## 16 -0.32397454  0.64019528  0.36177440
## 17  0.46087245  0.61660868  0.67675736
## 18  0.41732553  0.71111234 -1.21591003
## 19  1.04120451  0.11986429 -0.63306531
## 20  1.39709979 -0.53060810  0.67003668
## 21  1.64046682 -1.54053971 -1.66879107
## 22 -0.51860317  2.45956726  1.36247833
## 23  1.10286061  0.60454851  1.70393647
## 24  1.91924719 -1.14973539 -0.64690007
## 25  0.84451316  1.55315562  0.06346943
## 26  1.74614534 -0.60978067  1.36484924
## 27  1.23431051  0.91240050 -2.16643033
## 28  2.60386762  0.08178357  0.87401532
## 29  1.00595356  1.50672401 -1.44179251
## 30  2.43463278  0.39714752  0.92229467

3. 線性組合的係數 (主成份和原變數的關係)

每一個主成份，都是原變數經過線性組合後產生的值。

若有 n 個原變數(columns)，則可以投影出 n 個主成份，其線性組合的格式如下：

其中，觀察原變數在線性組合中的係數(Φ；特徵向量)，可以知道主成份和原變數之間的關係(例如：正影響/負影響；影響程度多大)。

從計算出的 pca 中可以取出 rotation 資訊，其對應的意義為 eigen vector (線性組合中的係數)

【註】在 R 的官方文件中，描述此資訊為「變數負荷」(variable loadings)，容易跟「主成份負荷」(loadings of PC; correlation coefficient between PCs and original varialbes)混淆，此兩者在學術上的意義不同，特此提醒。

# 特徵向量(原變數的線性組合)
pca$rotation

##             PC1        PC2         PC3          PC4         PC5
## H1B -0.40991503 -0.4681242  0.07174689  0.056704066 -0.07882016
## H2B -0.51441491 -0.2004156  0.01669591  0.255448162 -0.46809834
## H3B  0.01853759 -0.5595940 -0.19427151 -0.004051477  0.71490431
## HR  -0.34336124  0.5417488 -0.03416307 -0.394140194  0.26396281
## RBI -0.64629912  0.1016251 -0.25396353 -0.156840299  0.20084751
## SB   0.16722000 -0.2741655 -0.52853255 -0.679207860 -0.39181790
## BB   0.05866272  0.2203673 -0.78218985  0.538744635 -0.00676713
##             PC6         PC7
## H1B -0.66701643 -0.39159028
## H2B  0.62315846 -0.14911690
## H3B  0.34921449 -0.12535585
## HR   0.10991843 -0.59189466
## RBI -0.12239863  0.65387482
## SB   0.04037564 -0.03239357
## BB  -0.12695740 -0.17252886

取前三個主成份的特徵向量：

top3.pca.eigenvector <- pca$rotation[, 1:3]
top3.pca.eigenvector

##             PC1        PC2         PC3
## H1B -0.40991503 -0.4681242  0.07174689
## H2B -0.51441491 -0.2004156  0.01669591
## H3B  0.01853759 -0.5595940 -0.19427151
## HR  -0.34336124  0.5417488 -0.03416307
## RBI -0.64629912  0.1016251 -0.25396353
## SB   0.16722000 -0.2741655 -0.52853255
## BB   0.05866272  0.2203673 -0.78218985

我們可以繪製主成份負荷圖，觀察原變數和主成份之間的關係：

first.pca <- top3.pca.eigenvector[, 1]   #  第一主成份
second.pca <- top3.pca.eigenvector[, 2]  #  第二主成份
third.pca <- top3.pca.eigenvector[, 3]   #  第三主成份

(以下有用到排序的技巧，可以參考前篇筆記內的order(), sort())

第一主成份

SB(盜壘)、BB(四壞球)與PC-1呈現正相關，看起來和「上壘」有關。

# 第一主成份：由小到大排序原變數的係數
first.pca[order(first.pca, decreasing=FALSE)]

##         RBI         H2B         H1B          HR         H3B          BB 
## -0.64629912 -0.51441491 -0.40991503 -0.34336124  0.01853759  0.05866272 
##          SB 
##  0.16722000

# 使用dotchart，繪製主成份負荷圖
dotchart(first.pca[order(first.pca, decreasing=FALSE)] ,   # 排序後的係數
         main="Loading Plot for PC1",                      # 主標題
         xlab="Variable Loadings",                         # x軸的標題
         col="red")                                        # 顏色

第二主成份

HR(全壘打)、BB(四壞球)、RBI(打點)與PC-2呈現正相關，看起來和「打擊者」有關。

# 第二主成份：由小到大排序原變數的係數
second.pca[order(second.pca, decreasing=FALSE)]

##        H3B        H1B         SB        H2B        RBI         BB 
## -0.5595940 -0.4681242 -0.2741655 -0.2004156  0.1016251  0.2203673 
##         HR 
##  0.5417488

# 使用dotchart，繪製主成份負荷圖
dotchart(second.pca[order(second.pca, decreasing=FALSE)] ,  # 排序後的係數
         main="Loading Plot for PC2",                       # 主標題
         xlab="Variable Loadings",                          # x軸的標題
         col="blue")                                        # 顏色

第三主成份

H1B(一壘安打)、H2B(二壘安打)與PC-3呈現正相關，看起來和「安打」有關。

# 第三主成份：由小到大排序原變數的係數
third.pca[order(third.pca, decreasing=FALSE)]

##          BB          SB         RBI         H3B          HR         H2B 
## -0.78218985 -0.52853255 -0.25396353 -0.19427151 -0.03416307  0.01669591 
##         H1B 
##  0.07174689

# 使用dotchart，繪製主成份負荷圖
dotchart(third.pca[order(third.pca, decreasing=FALSE)] ,   # 排序後的係數
         main="Loading Plot for PC3",                      # 主標題
         xlab="Variable Loadings",                         # x軸的標題
         col="purple")                                     # 顏色

我們也可以繪製另一種主成份負荷圖，觀察每個球隊擅長的特性是什麼：

右邊的球隊適合上壘，多以盜壘(SB)和四壞球(BB)保送見長，全壘打表現中等(e.g.編號19)
左上方的球隊以力量取勝，在全壘打(HR)和打點(RBI)上有顯著的優勢(e.g.編號11)
下方的球隊不擅長全壘打，但在安打上的表現遠勝於其他球隊，盜壘也有一定水準(e.g.編號5)

# 選取 PC1 和 PC2 繪製主成份負荷圖
biplot(pca, choices=1:2)

(註)

根據2012年美國職棒的統計：

編號19的洛杉磯道奇隊，累積97次盜壘，只有101支全壘打。
編號11的紐約洋基隊，全壘打產量累積224支，高居全聯盟之冠。
編號5的舊金山巨人隊，全壘打累積95支，是所有球隊唯一未破百的隊伍；但一壘安打卻將近1000支，盜壘109次。

總結

本篇筆記以2012美國的職棒資料進行主成份分析，練習如何挑選幾個主成份，並且根據主份負荷進行解釋上的探討。

主成份分析，因為會對原始資料進行轉軸，因此有時候會比較難解釋；不過換個角度思考，只要能從主成份中找到有趣的故事，那麼這一次分析往往就能發現有價值的策略或insight。

It’s still a long way to go~

Reference

本篇筆記參考R的世界-主成分分析製作而成。

6. R and packages version

這是本篇筆記使用R跟套件版本：

pkgs = c("base")
r_version = paste("R", getRversion())
pkg_version = c()
for (package_name in pkgs) {
    pkg_version = c(pkg_version, 
                    paste(package_name,packageVersion(package_name)))
}
c(r_version, pkg_version)

## [1] "R 3.5.3"    "base 3.5.3"

R筆記–(7)主成份分析(2012美國職棒MLB)

skydome20

2019/03/16