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函數。
2 拆解主成分分析步驟實際應用時我們通常會選擇主成分分析函數，直接把input數據一步分析到位，只需要看懂輸出結果即可。但是為了加深理解，這裡一步步拆解主成分分析步驟，講解原理。
2.1 測試數據數據集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.
這些是專業領域的用詞，大家可以先不用糾結具體細節。
2.2 為什麼要做主成分分析不管三七二十一就直接套用統計方法都是耍流氓，做主成分分析並不是拍腦袋決定的。 在這個例子裡面，我們拿到了這43個法官的12個信息，就可以通過這12個指標來對法官進行分類，但也許實際情況下收集其他法官的12個信息比較麻煩，所以我們希望只收集三五個信息即可，然後也想達到比較好的分類效果。或者至少也得剔除幾個指標吧，這個時候主成分分析就能派上用場啦。這12個變量能得到12個主成分，如果前兩個主成分可以揭示85%以上的變異度，也就是說我們可以用兩個主成分來代替那12個指標。
在生物信息學領域，比如我們測了1000個病人的2萬個基因的表達矩陣，同時也有他們的健康狀態信息。那麼我們想仔細研究這些數據，想得到基因表達與健康狀態的某種關係。這樣我就可以對其餘幾十億的人檢測基因表達來預測其健康狀態。如果我們進行了主成分分析，就可以選擇解釋度比較高的主成分對應的基因，可能就幾十上百個而已，大幅度的降低廣泛的基因檢測成本。
2.3 step1:數據標準化(中心化)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)，本例採用中心化。
2.4 step2:求相關係數矩陣dat_cor=cor(dat_scale)
options(digits=4, scipen=4)
kable(dat_cor)
從相關係數看，只有 CONT 變量跟其它變量沒有相關性。
當然，這樣的矩陣不易看清楚規律，很明顯，畫個熱圖就一眼看出來了。
2.5 step3:計算特徵值和特徵向量利用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。
2.6 step4:崖低碎石圖和累積貢獻圖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.
但是在實際的生物信息學應用中，通常達不到這個要求。
2.7 step5:主成分載荷主成分載荷表示各個主成分與原始變量的相關係數。
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指標跟其它指標表現完全不一樣，第一個主成分很明顯跟除了CONT之外的所有其它指標負相關，而第二個主成分則主要取決於CONT指標。
2.8 step6:主成分得分計算和圖示將中心化的變量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維數據。
PCA本質上是將方差最大的方向作為主要特徵，並且在各個正交方向上將數據「離相關」，也就是讓它們在不同正交方向上沒有相關性。
PCA也存在一些限制，例如它可以很好的解除線性相關，但是對於高階相關性就沒有辦法了，對於存在高階相關性的數據，可以考慮Kernel PCA，通過Kernel函數將非線性相關轉為線性相關，關於這點就不展開討論了。另外，PCA假設數據各主特徵是分布在正交方向上，如果在非正交方向上存在幾個方差較大的方向，PCA的效果就大打折扣了。
最後需要說明的是，PCA是一種無參數技術，也就是說面對同樣的數據，如果不考慮清洗，誰來做結果都一樣，沒有主觀參數的介入，所以PCA便於通用實現，但是本身無法個性化的優化。
3 實戰一比如你要做一項分析人的糖尿病的因素有哪些，這時你設計了10個你覺得都很重要的指標，然而這10個指標對於你的分析確實太過繁雜，這時你就可以採用主成分分析的方法進行降維。10個指標之間會有這樣那樣的聯繫，相互之間會有影響，通過主成分分析後，得到三五個主成分指標。此時,這幾個主成分指標既涵蓋了你10個指標中的絕大部分信息，這讓你的分析得到了簡化（從10維降到3、5維）。
下面是442個糖尿病人相關的數據集，具體如下：
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)
#
length(y)
#
summary(y)
#
#
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
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,])
預測主成份的值，這裡用的就是訓練數據，所以得出訓練數據主成分的值。
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)
#
#
正文中的princomp()函數為基礎包中進行主成分分析的函數。 另外，R中psych包中提供了一些更加豐富有用的函數，這裡列出幾個相關度較高的函數，以供讀者了解。
還有很多主成分分析結果可視化包，在直播我的基因組裡面都提到過。
6 推薦一個R包factoextrafactoextra是一個R包，易於提取和可視化探索性多變量數據分析的輸出，包括：
主成分分析（PCA），用於通過在不丟失重要信息的情況下降低數據的維度來總結連續（即定量）多變量數據中包含的信息。
對應分析（CA）是適用於分析由兩個定性變量（或分類數據）形成的大型應變表的主成分分析的擴展。
多重對應分析（MCA），它是CA對包含兩個以上分類變量的數據表的適應。
多因素分析（MFA）專門用於數據集，其中變量被組織成組（定性和/或定量變量）。
分層多因素分析（HMFA）：在將數據組織成層次結構的情況下，MFA的擴展。
混合數據因子分析（FAMD）是MFA的一個特例，專門用於分析包含定量和定性變量的數據集。
有許多R包實現主要組件方法。這些包包括：FactoMineR，ade4，stats，ca，MASS和ExPosition。然而，根據使用的包，結果呈現不同。為了幫助解釋和多變量分析的可視化（如聚類分析和維數降低分析），所以作者開發了一個名為factoextra的易於使用的R包。
7.主成分分析的生物信息學應用比如對千人基因組計劃的對VCF突變數據進行主成分分析來區分人種：https://www.biostars.org/p/185869/
再比如每次做表達矩陣都需要根據樣本分組信息PCA分析及可視化看看分組是否可靠。
詳細例子見：http://github.com/jmzeng1314/GEO/
然後就是單細胞轉錄組數據也經常會PCA看看分群，或者PCA來去除前幾個主成分因素來抹掉某些影響等等。
8. 主成分分析的其它可視化方法看到一個包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也是正常。
9.參考資料：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）
生信基礎知識大全系列：生信基礎知識100講
(包括統計學，測序原理，生物學，計算機基礎，繪圖及可視化，醫學)
生信小技巧系列B站高清視頻：
https://www.bilibili.com/video/av25131640/
生信小技巧第09課-生信基礎資料大全 這些資料長期有效。
史上最強的生信自學環境準備課來啦！！ 7次改版，11節課程，14K的講稿，30個夜晚打磨，100頁PPT的課程。
還有GEO課程：請你師妹喝奶茶，不如一起看GEO 感興趣的可以推薦給自己的師弟師妹哦。
當然了，閱讀原文有驚喜！！！