主成分分析(PCA)在R 及 Python中的實戰指南

大數據文摘作品，轉載要求見文末

編譯團隊|李小帥，姚佳靈

有太多不如沒有！如果一個數據集有太多變量，會怎麼樣？這裡有些可能的情況你也許會碰上——

1.你發現大部分變量是相關的。2.你失去耐心，決定在整個數據集上建模。這個模型返回很差的精度，於是你的感覺很糟糕。3.你變得優柔寡斷，不知道該做什麼。4.你開始思考一些策略方法來找出幾個重要變量。

相信我，處理這樣的情形不是像聽上去那樣難。統計技術，比如，因子分析，主成分分析有助於解決這樣的困難。在本文中，我詳細地解釋了主成分分析的概念。我一直保持說明簡要而詳實。為了操作上的理解，我也演示了在R使用這個技術並帶有解釋。

注意：

要理解本文的內容，需要有統計學的知識。

 什麼是主成分分析？

▼

簡而言之，主成分分析是一種從一個數據集的一大組可用變量中提取重要變量的方法。它從高維度數據集中提取出低維度特徵變量集合，並儘可能多地捕捉到信息。變量越少，數據可視化也變得更有意義。處理3維或者更高維度的數據集時，主成分分析方法更有效。

它總是在一個對稱相關或協方差矩陣上施行。這意味著矩陣應該是數值型的，並且有著標準化的數據。

讓我們通過一個例子來理解：

假設我們有一個300(n) X 50(p)維度的數據集。n代表著樣本集數量，p代表著預測值的數目。由於我們有個很大的p值，p = 50，因此，會有p(p-1)/2個散布圖，也就是說，有可能超過1000個散布圖需要分析變量間的關係。在這樣的數據集中做探索分析是不是一件非常繁瑣的事呀？

在這樣的情況下，選取一個捕捉到儘可能多信息的預測值子集p（p<<50）是一個明晰的辦法。接著在得到的低維度空間繪製觀察結果。

下圖顯示了利用主成分分析從高維度（三維）數據到低維度（二維）數據的轉換。請別忘了，每一個所得到的維度都是特徵p的線性組合。

圖片來源：nlpca

 什麼是主成分？

▼

主成分是數據集中的初始預測值規範化後的線性組合。在上圖中，PC1和PC2便是主成分。假設我們有一個預測值集合：X¹，X²...，Xp

主成分可以寫成：

Z¹ = Φ¹¹X¹ + Φ²¹X² + Φ³¹X³ + .... +Φp¹Xp

其中——

◇ Z¹是第一主成分

◇ Φp¹是構成第一主成分負載量（Φ¹,      Φ²…）的加載向量      。該向量被限制成模長為1。這是因為加載向量的數值巨大的模也許會導致巨大的差異。它還定義了沿著數據變化最大的主成分（Z¹）的方向。這樣一來，它使得在P維度空間中存在一條最接近n樣本集的直線。擬合的程度由歐式距離平方均值來衡量。

◇ X¹..Xp 是規範化後的預測值。規範化後的預測值的均值為0、標準差為1。

因此，

第一主成分是在數據集中捕捉最大方差的初始預測變量的線性組合。它決定了數據中最高變異性的方向。在第一主成分中，捕捉到的變異性越大，成分捕捉到的信息就越多。沒有比第一主成分有更高變異性的成分。

第一主成分形成一條最接近數據的直線，也就是說，它把數據點和該直線之間的距離平方和最小化了。

類似地，我們也能夠計算第二主成分。

第二主成分（Z²）也是捕捉到數據集中剩餘方差的線性組合，和第一主成分（Z¹）不相關。換句話說，第一主成分與第二主成分間的相關係數為0。它可以表示成：

Z² = Φ¹²X¹ + Φ²²X² + Φ³²X³ + .... + Φp2Xp

如果兩個成分是不相關的，那麼兩者應該是正交的（見下圖）。下圖是在模擬數據上用兩個預測值繪製的。需要注意的是，主成分的方向，正如預期的那樣，是正交的。這表明在這兩個主成分之間的相關係數為0。

所有接下來的主成分都跟從著相似的概念，即它們捕捉前一個主成分剩餘的變化，並與之不相關。一般而言，對於n x p維度的數據，能夠構建最小的主成分向量為（n-1, p）。

這些主成分的方向是以無監督的方式確定的，也即響應變量（Y）不是用來決定主成分方向的。因此，這是以無監督的方式。

注意：偏最小二乘法（Partial least square，簡稱PLS）是替代主成分分析的一種監督方法。偏最小二乘法分配較高的權重給與響應變量y具有強相關關係的變量，以此決定主成分。

 為什麼變量規範化是必須的？

▼

主成分是由原始預測數據規範化後提供的。這是因為原始預測數據可能具有不同的範圍尺度。例如，想像一下這麼一個數據集，在該數據集中存在很多變量的度量單位：加侖、公裡、光年等等。可以肯定的是在這些變量中的方差範圍會很大。

在沒有規範化的變量上執行主成分分析會導致帶有高方差變量近乎瘋狂的大量的負荷。反過來，這將導致一個主成分依賴於具有高方差的變量。這不是我們所希望的。

如下圖所示，主成分分析在一個數據集上執行了兩次（帶有未縮放和縮放的預測值）。該數據集有大約40個變量，正如你所見，第一主成分由變量Item_MRP所主導。同時，第二主成分由變量Item_Weight主導。這種主導普遍存在是因為變量有相關的高方差。當變量被縮放後，我們便能夠在二維空間中更好地表示變量。

 在Python & R中應用 

 主成分分析方法 

（帶有代碼註解）

▼

要選多少主成分？我可以深入研究理論，但更好是用編程實戰來回答這一問題。

作為演示示例，我將使用來自BIg Mart Prediction Challenge上的數據。

請記住，主成分分析僅能應用於數值型數據，因此，如果數據集中存在分類變量，必須將其轉換成數值型的。而且在應用這個技術前前，必須進行了基本的數據清理。讓我們快點完成原始數據的加載和清理步驟：

#目錄路徑

> path <- ".../Data/Big_Mart_Sales"

#設定工作目錄

> setwd(path)

#加載訓練和測試文件

> train <- read.csv("train_Big.csv")

> test <- read.csv("test_Big.csv")

#增加一列

> test$Item_Outlet_Sales <- 1

#合併數據集

> combi <- rbind(train, test)

#用中位值替換缺失值

> combi$Item_Weight[is.na(combi$Item_Weight)] <- median(combi$Item_Weight, na.rm = TRUE)

#用中位值替換0

> combi$Item_Visibility <- ifelse(combi$Item_Visibility == 0, median(combi$Item_Visibility),                                   combi$Item_Visibility)

#找出模式並替換

> table(combi$Outlet_Size, combi$Outlet_Type)

> levels(combi$Outlet_Size)[1] <- "Other"

至此，我們已對缺失值做了替換處理。現在剩下的都是除去了依賴性（響應）變量和其它標識符變量（如果存在的話）。正如上面所講，我們正在練習無監督學習技術，因此，響應變量必須除去。

#除去依賴性和標識符變量

> my_data <- subset(combi, select = -c(Item_Outlet_Sales, Item_Identifier,                                       Outlet_Identifier))

Let’s check the available variables ( a.k.a predictors) in the data set.

現在，檢查一下數據集中的可用變量（也即預測值）：

#檢查可用變量

> colnames(my_data)

由於主成分分析作用於數值型變量上，讓我們看看是否有不是數值型的變量。

#檢查變量的類

> str(my_data)

'data.frame': 14204 obs. of 9 variables:

$ Item_Weight : num 9.3 5.92 17.5 19.2 8.93 ...

$ Item_Fat_Content : Factor w/ 5 levels "LF","low fat",..: 3 5 3 5 3 5 5 3 5 5 ...

$ Item_Visibility : num 0.016 0.0193 0.0168 0.054 0.054 ...

$ Item_Type : Factor w/ 16 levels "Baking Goods",..: 5 15 11 7 10 1 14 14 6 6 ...

$ Item_MRP : num 249.8 48.3 141.6 182.1 53.9 ...

$ Outlet_Establishment_Year: int 1999 2009 1999 1998 1987 2009 1987 1985 2002 2007 ...

$ Outlet_Size : Factor w/ 4 levels "Other","High",..: 3 3 3 1 2 3 2 3 1 1 ...

$ Outlet_Location_Type : Factor w/ 3 levels "Tier 1","Tier 2",..: 1 3 1 3 3 3 3 3 2 2 ...

$ Outlet_Type : Factor w/ 4 levels "Grocery Store",..: 2 3 2 1 2 3 2 4 2 2 ... -

可惜，9個變量中有6個本質上是分類變量。現在，我們得做些額外工作。我們將使用一位有效編碼將分類變量轉換成數值型。

#載入庫

> library(dummies)

#創建一個虛擬的數據幀

> new_my_data <- dummy.data.frame(my_data, names = c("Item_Fat_Content","Item_Type",

                                "Outlet_Establishment_Year","Outlet_Size",

                                "Outlet_Location_Type","Outlet_Type"))

檢查一下，現在我們是否有了一個整數值數據集，只要簡單地寫：

 #檢查數據集

> str(new_my_data)

是的，現在數據類型全部為數值型的。現在我們能夠繼續工作，應用主成分分析了。

基本R函數prcomp()用來實施主成分分析。默認情況下，它讓變量集中擁有等於0的均值。用上參數scale. = T，我們規範化變量使得標準偏差為1。

#主成分分析

> prin_comp <- prcomp(new_my_data, scale. = T)

> names(prin_comp)

[1] "sdev"     "rotation" "center"   "scale"    "x"

prcomp()函數形成5種實用操作：

1. 中心和規模是指在實施主成分分析前用於標準化變量的各均值和標準偏差

#輸出變量的均值

prin_comp$center

#輸出變量的標準偏差

prin_comp$scale

2.旋轉措施提供主成分的負載。旋轉矩陣的每一列包含主成分負載向量。這是我們應該感興趣的最重要措施。

它返回44個主成分負載。正確嗎？當然。在一個數據集中，主成分負載的最大值至少為(n-1, p)。下面我們看看前4個主成分和前5行。

> prin_comp$rotation[1:5,1:4]

                                PC1            PC2            PC3             PC4

Item_Weight                0.0054429225   -0.001285666   0.011246194   0.011887106

Item_Fat_ContentLF        -0.0021983314    0.003768557  -0.009790094  -0.016789483

Item_Fat_Contentlow fat   -0.0019042710    0.001866905  -0.003066415  -0.018396143

Item_Fat_ContentLow Fat    0.0027936467   -0.002234328   0.028309811   0.056822747

Item_Fat_Contentreg        0.0002936319    0.001120931   0.009033254  -0.001026615

3.為了計算主成分評價向量，我們不必將數據和負載相乘。相反，矩陣X具有14204 x 44 維度的主成分評價向量。

> dim(prin_comp$x)

[1] 14204    44

讓我們來繪製產生的主成分——

> biplot(prin_comp, scale = 0)

參數scale = 0確保上圖中箭頭的縮放代表負載。為了從上圖中作出推斷，關注圖中的最末端（上、下、左、右）。

我們推斷第一主成分與Outlet_TypeSupermarket、Outlet_Establishment_Year2007的量度對應。類似的，第二成分和Outlet_Location_TypeTier1、Outlet_Sizeother的量度對應。對應一個成分裡的某個變量的精確量度，應該再看看旋轉後的矩陣（如上）。

4.prcomp函數也提供了計算每一個主成分標準偏差的便利。sdev是指主成分標準偏差。

#計算每個主成分的標準偏差

> std_dev <- prin_comp$sdev

#計算方差

> pr_var <- std_dev^2

#查看前10個成分的方差

> pr_var[1:10]

[1] 4.563615 3.217702 2.744726 2.541091 2.198152 2.015320 1.932076 1.256831

[9] 1.203791 1.168101

我們的目標是尋找能夠說明最大方差的成分，這是因為，我們想在使用這些成分時儘可能多地保留信息。因此，如果用來說明的方差越大，那麼這些成分包含的信息也就越多。

為計算被每個主成分解釋的方差的佔比，我們簡單地將該方差除以方差總和。結果如下：

#被解釋的方差的佔比

> prop_varex <- pr_var/sum(pr_var)

> prop_varex[1:20]

[1] 0.10371853 0.07312958 0.06238014 0.05775207 0.04995800 0.04580274

[7] 0.04391081 0.02856433 0.02735888 0.02654774 0.02559876 0.02556797

[13] 0.02549516 0.02508831 0.02493932 0.02490938 0.02468313 0.02446016

[19] 0.02390367 0.02371118

如上結果顯示，第一主成分能夠說明10.3%的方差。第二成分能夠說明7.3%的方差，第三成分說明了6.2%的方差等等。那麼，對於建模階段，我們究竟需要選用多少成分呢？

用碎石圖可以解決上面的問題。碎石圖用來訪問成分或說明數據中最可變性的因素。它代表降序排列值。

#碎石圖

> plot(prop_varex, xlab = "Principal Component",

             ylab = "Proportion of Variance Explained",

             type = "b")

上圖顯示：約30個成分說明了數據集中98.4%的方差。換句話說，利用主成分分析算法，我們將預測值從44個降到30個，而不影響說明的方差。這就是主成分分析算法的強大之處。讓我們通過繪製一個累計方差圖做確認核查。它將向我們展示成分數量的清晰畫面。

#cumulative scree plot

> plot(cumsum(prop_varex), xlab = "Principal Component",

              ylab = "Cumulative Proportion of Variance Explained",

              type = "b")

上圖顯示，30個成分對方差的影響接近98%。因此，在這個案例中，我們選擇30種成分（PC1到PC30），並且用在建模階段。這個使得在訓練集上實施主成分分析的步驟變得完整了。對於建模，我們將使用30個成分作為預測變量並按照正常的過程進行。

 用主成分分析成分預測建模 

▼

我們在訓練集上完成主成分計算之後，現在讓我們理解利用這些成分在測試數據上做預測的過程。這個過程是簡單的。就像我們已經在訓練集上獲得主成分分析成分那樣，我們將在測試集上取另外一組成分。最後，我們訓練模型。

但是，要理解幾個要點：

我們不應該把訓練集和測試集合在一起來一次性地獲得整個數據的主成分分析成分。因為，由於測試數據會「洩露」到訓練集中，這會違背整個概括假設。換句話說，測試數據集不再保持「沒看見」的狀態。最終，這會打擊模型的泛化能力。

我們不應該在測試和訓練數據集上分開進行主成分分析。因為，來自訓練和測試的主成分的組合向量將有不同的方向（方差不同的緣故）。由於這個原因，我們最終會比較在落在不同軸上的數據。這樣，來自訓練和測試數據的結果向量應該有相同的軸。

 那麼，我們應該做什麼？

▼

我們應該像我們在訓練集上所做的一樣，在測試集上做相同的轉換，包括集中和度量特徵。讓我們在R中做一下：

#加上帶主成分的訓練集

> train.data <- data.frame(Item_Outlet_Sales = train$Item_Outlet_Sales, prin_comp$x)

#we are interested in first 30 PCAs

#我們對前30個主成分感興趣

> train.data <- train.data[,1:31]

#運行決策樹

> install.packages("rpart")

> library(rpart)

> rpart.model <- rpart(Item_Outlet_Sales ~ .,data = train.data, method = "anova")

> rpart.model

#把測試轉換成主成分分析

> test.data <- predict(prin_comp, newdata = pca.test)

> test.data <- as.data.frame(test.data)

#選擇前30個成分

> test.data <- test.data[,1:30]

#在測試數據上做出預測

> rpart.prediction <- predict(rpart.model, test.data)

#最後查看你的分數排行榜，這只是為了好玩

> sample <- read.csv("SampleSubmission_TmnO39y.csv")

> final.sub <- data.frame(Item_Identifier = sample$Item_Identifier, Outlet_Identifier = sample$Outlet_Identifier, Item_Outlet_Sales = rpart.prediction)

> write.csv(final.sub, "pca.csv",row.names = F)

這是提取主成分之後完整的建模過程。我保證你在上傳解決方案後不會對你的分數排行榜感到高興。試試用下隨機森林。

對於Python用戶：為了在Python中運行主成分分析，只需從sklearn庫導入主成分分析。和上文提到的對R用戶的解釋是一樣的。當然，用Python的結果是用R後派生出來的。Python中所用的數據集是清洗後的版本，缺失值已經被補上，分類變量被轉換成數值型。建模過程保持不變，和上面對R用戶所說的一樣。

import numpy as np

from sklearn.decomposition import PCA

import pandas as pd

import matplotlib.pyplot as plt

from sklearn.preprocessing import scale

%matplotlib inline

#加載數據集

data = pd.read_csv('Big_Mart_PCA.csv')

#轉換成數字型數組

X=data.values

#Scaling the values

X = scale(X)

pca = PCA(n_components=44)

pca.fit(X)

#每一個PC說明的方差數量

var= pca.explained_variance_ratio_

#Cumulative Variance explains

var1=np.cumsum(np.round(pca.explained_variance_ratio_, decimals=4)*100)

print var1

[  10.37   17.68   23.92   29.7    34.7    39.28   43.67   46.53   49.27

51.92   54.48   57.04   59.59   62.1    64.59   67.08   69.55   72.

74.39   76.76   79.1    81.44   83.77   86.06   88.33   90.59   92.7

94.76   96.78   98.44  100.01  100.01  100.01  100.01  100.01  100.01

100.01  100.01  100.01  100.01  100.01  100.01  100.01  100.01]

plt.plot(var1)

#看看我採用了30個變量的輸出圖

pca = PCA(n_components=30)

pca.fit(X)

X1=pca.fit_transform(X)

print X1

要點回顧——

◇主成分分析被用來克服數據集中的冗餘。

◇這些特徵具有低維的性質。

◇這些特徵（也即成分）是原始預測變量規範化線性組合形成的結果。

◇這些成分旨在用高可釋方差抓取儘可能多的信息。

◇第一成分有最高的方差，接下來是第二和第三成分，以此類推。

◇這些成分必須是不相關的，請參考前文所述。

◇規範化數據在預測值用不同單位測量時變得極其重要。

◇主成分分析在3維及以上維度的數據集中最有成效。因為，維度越高，就越難從最終的數據云做出解釋。

◇主成分分析應用於數值型變量的數據集上。

◇主成分分析是一個工具，有助於生成高維度數據的更好的可視化。

來源：http://www.analyticsvidhya.com/blog/2016/03/practical-guide-principal-component-analysis-python/

如需轉載，請在開篇顯著位置註明作者和出處（轉自：大數據文摘 |bigdatadigest），並在文章結尾放置大數據文摘醒目二維碼。無原創標識文章請按照轉載要求編輯，可直接轉載，轉載後請將轉載連結發送給我們；有原創標識文章，請發送【文章名稱-待授權公眾號名稱及ID】給我們申請白名單授權。未經許可的轉載以及改編者，我們將依法追究其法律責任。聯繫郵箱：[e-mail]。
回復「志願者」了解如何加入我們


Kaggle | 使用Python和R繪製數據地圖的十七個經典案例（附資源）