1 背景
2 拆解主成分分析步驟
2.1 測試數據
2.2 為什么要做主成分分析
2.3 step1:數據標準化(中心化)
2.4 step2:求相關系數矩陣
2.5 step3:計算特征值和特征向量
2.6 step4:崖低碎石圖和累積貢獻圖
2.7 step5:主成分載荷
2.8 step6:主成分得分計算和圖示
3 實戰一
4 實戰二
5 進階的主成分分析-psych包
6 推薦一個R包factoextra
7.主成分分析的生物信息學應用
8. 主成分分析的其它可視化方法
9.其它學習資料
主成分分析法是數據挖掘中常用的一種降維算法,是Pearson在1901年提出的,再后來由hotelling在1933年加以發展提出的一種多變量的統計方法，其最主要的用途在于“降維”，通過析取主成分顯出的最大的個別差異,也可以用來削減回歸分析和聚類分析中變量的數目，與因子分析類似。
所謂降維，就是把具有相關性的變量數目減少，用較少的變量來取代原先變量。如果原始變量互相正交，即沒有相關性，則主成分分析沒有效果。
在生物信息學的實際應用情況中，通常是得到了成百上千個基因的信息，這些基因相互之間會有影響，通過主成分分析后，得到有限的幾個主成分就可以代表它們的基因了。也就是所謂的降維。
R語言有非常多的途徑做主成分分析，比如自帶的princomp()和psych包的principal()函數，還有gmodels包的fast.prcomp函數。
實際應用時我們通常會選擇主成分分析函數，直接把input數據一步分析到位，只需要看懂輸出結果即可。但是為了加深理解，這里一步步拆解主成分分析步驟，講解原理。
數據集USJudgeRatings包含了律師對美國高等法院法官的評分。數據框包含43個樣本，12個變量!
下面簡單看一看這12個變量是什么，以及它們的相關性。
library(knitr)
kable(head(USJudgeRatings))
這12個變量的介紹如下：
[,1] CONT Number of contacts of lawyer with judge.
[,2] INTG Judicial integrity.司法誠實性
[,3] DMNR Demeanor.風度；舉止；行為
[,4] DILG Diligence.勤奮，勤勉；注意的程度
[,5] CFMG Case flow managing.
[,6] DECI Prompt decisions.
[,7] PREP Preparation for trial.
[,8] FAMI Familiarity with law.
[,9] ORAL Sound oral rulings.
[,10] WRIT Sound written rulings.
[,11] PHYS Physical ability.
[,12] RTEN Worthy of retention.
這些是專業領域的用詞，大家可以先不用糾結具體細節。
不管三七二十一就直接套用統計方法都是耍流氓，做主成分分析并不是拍腦袋決定的。 在這個例子里面，我們拿到了這43個法官的12個信息，就可以通過這12個指標來對法官進行分類，但也許實際情況下收集其他法官的12個信息比較麻煩，所以我們希望只收集三五個信息即可，然后也想達到比較好的分類效果。或者至少也得剔除幾個指標吧，這個時候主成分分析就能派上用場啦。這12個變量能得到12個主成分，如果前兩個主成分可以揭示85%以上的變異度，也就是說我們可以用兩個主成分來代替那12個指標。
在生物信息學領域，比如我們測了1000個病人的2萬個基因的表達矩陣，同時也有他們的健康狀態信息。那么我們想仔細研究這些數據，想得到基因表達與健康狀態的某種關系。這樣我就可以對其余幾十億的人檢測基因表達來預測其健康狀態。如果我們進行了主成分分析，就可以選擇解釋度比較高的主成分對應的基因，可能就幾十上百個而已，大幅度的降低廣泛的基因檢測成本。
dat_scale=scale(USJudgeRatings,scale=F)
options(digits=4, scipen=4)
kable(head(dat_scale))
scale()是對數據中心化的函數，當參數scale=F時，表示將數據按列減去平均值，scale=T表示按列進行標準化，公式為（x-mean(x)）/sd(x)，本例采用中心化。
dat_cor=cor(dat_scale)
options(digits=4, scipen=4)
kable(dat_cor)
從相關系數看，只有 CONT 變量跟其它變量沒有相關性。
當然，這樣的矩陣不易看清楚規律，很明顯，畫個熱圖就一眼看出來了。
利用eigen函數計算相關系數矩陣的特征值和特征向量。這個是主成分分析法的精髓。
dat_eigen=eigen(dat_cor)
dat_var=dat_eigen$values ## 相關系數矩陣的特征值
options(digits=4, scipen=4)
dat_var
## [1] 10.133504 1.104147 0.332902 0.253847 0.084453 0.037286 0.019683
## [8] 0.015415 0.007833 0.005612 0.003258 0.002060
pca_var=dat_var/sum(dat_var)
pca_var
## [1] 0.8444586 0.0920122 0.0277418 0.0211539 0.0070377 0.0031072 0.0016402
## [8] 0.0012846 0.0006528 0.0004676 0.0002715 0.0001717
pca_cvar=cumsum(dat_var)/sum(dat_var)
pca_cvar
## [1] 0.8445 0.9365 0.9642 0.9854 0.9924 0.9955 0.9972 0.9984 0.9991 0.9996
## [11] 0.9998 1.0000
可以看出，PC1(84.4%)和PC2(9.2%)共可以解釋這12個變量的93.6的程度，除了CONT外的其他的11個變量與PC1都有較好的相關性，所以PC1與這11個變量基本斜交，而CONT不能被PC1表示，所以基本與PC1正交垂直，而PC2與CONT基本平行，表示其基本可以表示CONT。
library(ggplot2)
p=ggplot(,aes(x=1:12,y=pca_var))
p1=ggplot(,aes(x=1:12,y=pca_cvar))
p+geom_point(pch=2,lwd=3,col=2)+geom_line(col=2,lwd=1.2)
p1+geom_point(pch=2,lwd=3,col=2)+geom_line(col=2,lwd=1.2)
崖低碎石圖（scree plot）即貢獻率圖，是希望圖形一開始很陡峭，如懸崖一般，而剩下的數值都很小，如崖底的碎石一樣。
崖低碎石圖和累積貢獻率圖是對主成分貢獻率和累積貢獻率的一種直觀表示，用以作為選擇主成分個數的參考。本例中第一個主成分解釋總變異的84.4%,可以只選擇第一個主成分，但第二主成分也不小，因此選擇前兩個主成分。
主成分的個數選擇沒有一定之規，需按實際情況具體分析，一般要求累積貢獻率大于85%或特征值大于1.
但是在實際的生物信息學應用中，通常達不到這個要求。
主成分載荷表示各個主成分與原始變量的相關系數。
pca_vect= dat_eigen$vector ## 相關系數矩陣的特征向量
loadings=sweep(pca_vect,2,sqrt(pca_var),"*")
rownames(loadings)=colnames(USJudgeRatings)
options(digits=4, scipen=4)
kable(loadings[,1:2])
| CONT | 0.0028 | 0.2830 |
| INTG | -0.2652 | -0.0552 |
| DMNR | -0.2636 | -0.0599 |
| DILG | -0.2797 | 0.0110 |
| CFMG | -0.2780 | 0.0511 |
| DECI | -0.2774 | 0.0388 |
| PREP | -0.2843 | 0.0098 |
| FAMI | -0.2818 | -0.0004 |
| ORAL | -0.2874 | -0.0011 |
| WRIT | -0.2858 | -0.0095 |
| PHYS | -0.2580 | 0.0270 |
| RTEN | -0.2847 | -0.0119 |
結果表明，CONT指標跟其它指標表現完全不一樣，第一個主成分很明顯跟除了CONT之外的所有其它指標負相關，而第二個主成分則主要取決于CONT指標。
將中心化的變量dat_scale乘以特征向量矩陣即得到每個觀測值的得分。
pca_score=as.matrix(dat_scale)%*%pca_vect
head(pca_score[,1:2])
## [,1] [,2]
## AARONSON,L.H. 0.5098 -1.7080
## ALEXANDER,J.M. -2.2676 -0.8508
## ARMENTANO,A.J. -0.2267 -0.2903
## BERDON,R.I. -3.4058 -0.5997
## BRACKEN,J.J. 6.5937 0.2478
## BURNS,E.B. -2.3336 -1.3563
將兩個主成分以散點圖形式展示,看看這些樣本被這兩個主成分是如何分開的
p=ggplot(,aes(x=pca_score[,1],y=pca_score[,2]))+geom_point(color=USJudgeRatings[,1],pch=USJudgeRatings[,2])
p
對于主成分分析，不同數據會有不同的分析方法，應具體情況具體分析。
總結一下PCA的算法步驟：
設有m條n維數據。
1）將原始數據按列組成n行m列矩陣X
2）將X的每一行（代表一個屬性字段）進行零均值化，即減去這一行的均值
3）求出協方差矩陣
4）求出協方差矩陣的特征值及對應的特征向量
5）將特征向量按對應特征值大小從上到下按行排列成矩陣，取前k行組成矩陣P
6）Y=PX即為降維到k維后的數據
PCA本質上是將方差最大的方向作為主要特征，并且在各個正交方向上將數據“離相關”，也就是讓它們在不同正交方向上沒有相關性。
PCA也存在一些限制，例如它可以很好的解除線性相關，但是對于高階相關性就沒有辦法了，對于存在高階相關性的數據，可以考慮Kernel PCA，通過Kernel函數將非線性相關轉為線性相關，關于這點就不展開討論了。另外，PCA假設數據各主特征是分布在正交方向上，如果在非正交方向上存在幾個方差較大的方向，PCA的效果就大打折扣了。
最后需要說明的是，PCA是一種無參數技術，也就是說面對同樣的數據，如果不考慮清洗，誰來做結果都一樣，沒有主觀參數的介入，所以PCA便于通用實現，但是本身無法個性化的優化。
比如你要做一項分析人的糖尿病的因素有哪些，這時你設計了10個你覺得都很重要的指標，然而這10個指標對于你的分析確實太過繁雜，這時你就可以采用主成分分析的方法進行降維。10個指標之間會有這樣那樣的聯系，相互之間會有影響，通過主成分分析后，得到三五個主成分指標。此時,這幾個主成分指標既涵蓋了你10個指標中的絕大部分信息，這讓你的分析得到了簡化（從10維降到3、5維）。
下面是442個糖尿病人相關的數據集，具體如下：
x a matrix with 10 columns (自變量)
y a numeric vector (因變量)
library(lars)
library(glmnet)
data(diabetes)
attach(diabetes)
summary(x)
## age sex bmi
## Min. :-0.10723 Min. :-0.0446 Min. :-0.09028
## 1st Qu.:-0.03730 1st Qu.:-0.0446 1st Qu.:-0.03423
## Median : 0.00538 Median :-0.0446 Median :-0.00728
## Mean : 0.00000 Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.03808 3rd Qu.: 0.0507 3rd Qu.: 0.03125
## Max. : 0.11073 Max. : 0.0507 Max. : 0.17056
## map tc ldl
## Min. :-0.11240 Min. :-0.12678 Min. :-0.11561
## 1st Qu.:-0.03666 1st Qu.:-0.03425 1st Qu.:-0.03036
## Median :-0.00567 Median :-0.00432 Median :-0.00382
## Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.03564 3rd Qu.: 0.02836 3rd Qu.: 0.02984
## Max. : 0.13204 Max. : 0.15391 Max. : 0.19879
## hdl tch ltg
## Min. :-0.10231 Min. :-0.07639 Min. :-0.12610
## 1st Qu.:-0.03512 1st Qu.:-0.03949 1st Qu.:-0.03325
## Median :-0.00658 Median :-0.00259 Median :-0.00195
## Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.02931 3rd Qu.: 0.03431 3rd Qu.: 0.03243
## Max. : 0.18118 Max. : 0.18523 Max. : 0.13360
## glu
## Min. :-0.13777
## 1st Qu.:-0.03318
## Median :-0.00108
## Mean : 0.00000
## 3rd Qu.: 0.02792
## Max. : 0.13561
dim(x)
## [1] 442 10
length(y)
## [1] 442
summary(y)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 25 87 140 152 212 346
boxplot(y)
其中x矩陣含有10個變量，分別是：“age” “sex” “bmi” “map” “tc” “ldl” “hdl” “tch” “ltg” “glu” 它們都在一定程度上或多或少的會影響個體糖尿病狀態。
數據的詳細介紹見 Efron, Hastie, Johnstone and Tibshirani (2003) “Least Angle Regression” (with discussion) Annals of Statistics;
一步法主成分分析
上面我們把整個主成分分析步驟拆解開來講解具體原理，但是實際數據處理過程中我們通常是用現成的函數一步法完成主成分分析，而且這個是非常高頻的分析，所以R里面自帶了一個函數princomp()來完成主成分分析，如下：
data=x ## 這里的x是上面的 diabetes 數據集里面的 442個病人的10個生理指標
pca<-princomp(data, cor=FALSE)
cor是邏輯變量,當cor=TRUE表示用樣本的相關矩陣R做主成分分析,當cor=FALSE表示用樣本的協方差陣S做主成分分。
summary(pca)
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
## Standard deviation 0.09542 0.05811 0.05223 0.04649 0.03871 0.03693
## Proportion of Variance 0.40242 0.14923 0.12060 0.09555 0.06622 0.06027
## Cumulative Proportion 0.40242 0.55165 0.67225 0.76780 0.83402 0.89429
## Comp.7 Comp.8 Comp.9 Comp.10
## Standard deviation 0.03484 0.03132 0.013311 0.0044009
## Proportion of Variance 0.05366 0.04337 0.007832 0.0008561
## Cumulative Proportion 0.94794 0.99131 0.999144 1.0000000
可以看到前三個主成份的信息量也只有67.2%，達不到我們前面說到85%，所以很難說可以用這3個主成分去代替這10個生理指標來量化病人的狀態。
值得一提的是，如果你看懂了前面的主成分分析的拆解步驟，就應該明白有多少個變量就有多少個主成分，只是并不是所有的主成分都有意義，理想狀態下我們希望有限的幾個主成分就可以代替數量繁多的變量，尤其是生物信息學里面的基因表達矩陣，兩三萬個基因的表達情況我們希望幾十個基因就可以替代它們，因為那些基因之間是相互關聯的。
碎石圖
也可以畫出主成分的碎石圖，來決定選幾個主成分。
screeplot(pca, type='lines')
由碎石圖可以看出，第5個主成分之后，圖線變化趨于平穩，因此可以選擇前5個主成分做分析。
樣本分布的散點圖
根據前兩個主成分畫出樣本分布的散點圖。
biplot(pca)
上面這個圖似乎意義不大，因為大部分情況下都是需要結合樣本的分組信息來看看這些主成分是否可以把樣本比較不錯的分開。
** 獲得訓練數據前4個主成份的值 **
kable(predict(pca, data)[1:4,])
kable(data[1:4,])
預測主成份的值，這里用的就是訓練數據，所以得出訓練數據主成分的值。
R中自帶數據集data(Harman23.cor)數據集中包含305名受試者的8個身體測量指標
data(Harman23.cor)
kable(Harman23.cor[1:5])
## Warning in kable_markdown(x = structure(c("0", "0", "0", "0", "0", "0", :
## The table should have a header (column names)
## Warning in kable_markdown(x = structure("305", .Dim = c(1L, 1L), .Dimnames
## = list(: The table should have a header (column names)
正文中的princomp()函數為基礎包中進行主成分分析的函數。 另外，R中psych包中提供了一些更加豐富有用的函數，這里列出幾個相關度較高的函數，以供讀者了解。
還有很多主成分分析結果可視化包，在直播我的基因組里面都提到過。
factoextra是一個R包，易于提取和可視化探索性多變量數據分析的輸出，包括：
主成分分析（PCA），用于通過在不丟失重要信息的情況下降低數據的維度來總結連續（即定量）多變量數據中包含的信息。
對應分析（CA）是適用于分析由兩個定性變量（或分類數據）形成的大型應變表的主成分分析的擴展。
多重對應分析（MCA），它是CA對包含兩個以上分類變量的數據表的適應。
多因素分析（MFA）專門用于數據集，其中變量被組織成組（定性和/或定量變量）。
分層多因素分析（HMFA）：在將數據組織成層次結構的情況下，MFA的擴展。
混合數據因子分析（FAMD）是MFA的一個特例，專門用于分析包含定量和定性變量的數據集。
有許多R包實現主要組件方法。這些包包括：FactoMineR，ade4，stats，ca，MASS和ExPosition。然而，根據使用的包，結果呈現不同。為了幫助解釋和多變量分析的可視化（如聚類分析和維數降低分析），所以作者開發了一個名為factoextra的易于使用的R包。
比如對千人基因組計劃的對VCF突變數據進行主成分分析來區分人種：https://www.biostars.org/p/185869/
再比如每次做表達矩陣都需要根據樣本分組信息PCA分析及可視化看看分組是否可靠。
詳細例子見：http://github.com/jmzeng1314/GEO/
然后就是單細胞轉錄組數據也經常會PCA看看分群，或者PCA來去除前幾個主成分因素來抹掉某些影響等等。
看到一個包ropls可視化做的不錯，本來以為ropls肯定是一個正常的r包，應該是在cran里面，結果
> install.packages('ropls')
Warning in install.packages :
package 'ropls' is not available (for R version 3.4.3)
Warning in install.packages :
cannot open URL 'https://cran.rstudio.com/bin/macosx/el-capitan/contrib/3.4/PACKAGES.rds': HTTP status was '404 Not Found'
> BiocInstaller::biocLite('ropls')
BioC_mirror: https://bioconductor.org
Using Bioconductor 3.6 (BiocInstaller 1.28.0), R 3.4.3 (2017-11-30).
Installing package(s) 'ropls'
trying URL 'https://bioconductor.org/packages/3.6/bioc/bin/macosx/el-capitan/contrib/3.4/ropls_1.10.0.tgz'
Content type 'application/x-gzip' length 1223650 bytes (1.2 MB)
==================================================
downloaded 1.2 MB
現在什么包都往bioconductor里面丟真奇怪。但是作圖顏值還可以，大家可以看看：
后來仔細看標題，終于明白了。
ropls: PCA, PLS(-DA) and OPLS(-DA) for multivariate analysis and feature selection of omics data
構建的就是組學數據，所以放在bioconductor也是正常。
http://www.cs.otago.ac.nz/cosc453/student_tutorials/principal_components.pdf
https://www.cs.princeton.edu/picasso/mats/PCA-Tutorial-Intuition_jp.pdf
http://www.cs.umd.edu/~samir/498/PCA.pdf
http://www.yale.edu/ceo/Documentation/PCA_Outline.pdf
http://people.tamu.edu/~alawing/materials/ESSM689/pca.pdf （R相關）
http://www2.dc.ufscar.br/~cesar.souza/publications/pca-tutorial.pdf （2012）
