微信公眾號:生信小知識
關注可了解更多的教程及生信知識。問題或建議,請公眾號留言;
前言統計學背景知識協方差相關係數函數總結實例講解1.載入原始數據2.作主成分分析3.結果解讀4.畫主成分的碎石圖並預測5.PCA結果繪製後記
前言PCA分析大家肯定經常看到,但是你真的懂PCA分析的結果嗎?
圖我也會看,我只是不是很清楚PCA背後輸出結果的解讀而已。正好看到一篇不錯的博客,就把主要的知識點記錄下 。
reference:
統計學背景知識協方差可以通俗的理解為:兩個變量在變化過程中是同方向變化?還是反方向變化?同向或反向程度如何?
從數值來看,協方差的數值越大,兩個變量同向程度也就越大。反之亦然。
從公式出發來理解一下:
公式簡單翻譯一下是:如果有X,Y兩個變量,每個時刻的「X值與其均值之差」乘以「Y值與其均值之差」得到一個乘積,再對這每時刻的乘積求和並求出均值(其實是求「期望」,但就不引申太多新概念了,簡單認為就是求均值了)。
具體例子可以去知乎詳細查看:
https://www.zhihu.com/question/20852004
相關係數對於相關係數,我們從它的公式入手。一般情況下,相關係數的公式為:
翻譯一下:就是用X、Y的協方差除以X的標準差和Y的標準差。
所以,相關係數也可以看成協方差:一種剔除了兩個變量量綱影響、標準化後的特殊協方差。
既然是一種特殊的協方差,那它:
1、也可以反映兩個變量變化時是同向還是反向,如果同向變化就為正,反向變化就為負。
2、由於它是標準化後的協方差,因此更重要的特性來了:它消除了兩個變量變化幅度的影響,而只是單純反應兩個變量每單位變化時的相似程度。
具體例子可以去知乎詳細查看:
https://www.zhihu.com/question/20852004
函數總結注意:這裡的輸入數據,rownames是樣本名,colnames是樣本的特徵。(與正常數據正好相反,需要用t()來轉置數據)
princomp()主成分分析 可以從相關陣或者從協方差陣做主成分分析
fviz_pca_ind對princomp()結果進行展示
summary()提取主成分信息
loadings()顯示主成分分析或因子分析中載荷的內容
predict()預測主成分的值
screeplot()畫出主成分的碎石圖
biplot()畫出數據關於主成分的散點圖和原坐標在主成分下的方向
實例講解現有30名中學生身高、體重、胸圍、坐高數據,對身體的四項指標數據做主成分分析。
1.載入原始數據# 清空環境
rm(list = ls())
options(stringsAsFactors = F)
# 準備數據
if (T) {
test<-data.frame(
X1=c(148, 139, 160, 149, 159, 142, 153, 150, 151, 139,
140, 161, 158, 140, 137, 152, 149, 145, 160, 156,
151, 147, 157, 147, 157, 151, 144, 141, 139, 148),
X2=c(41, 34, 49, 36, 45, 31, 43, 43, 42, 31,
29, 47, 49, 33, 31, 35, 47, 35, 47, 44,
42, 38, 39, 30, 48, 36, 36, 30, 32, 38),
X3=c(72, 71, 77, 67, 80, 66, 76, 77, 77, 68,
64, 78, 78, 67, 66, 73, 82, 70, 74, 78,
73, 73, 68, 65, 80, 74, 68, 67, 68, 70),
X4=c(78, 76, 86, 79, 86, 76, 83, 79, 80, 74,
74, 84, 83, 77, 73, 79, 79, 77, 87, 85,
82, 78, 80, 75, 88, 80, 76, 76, 73, 78)
)
rownames(test) <- paste0("student_",1:30)
colnames(test) <- c("height","weight","chest","sit-h")
}
# PCA分析
if (T) {
test.pr<-princomp(test,cor=TRUE)
summary(test.pr,loadings=TRUE)
}
# Importance of components:
# Comp.1 Comp.2 Comp.3 Comp.4
# Standard deviation 1.8817805 0.55980636 0.28179594 0.25711844
# Proportion of Variance 0.8852745 0.07834579 0.01985224 0.01652747
# Cumulative Proportion 0.8852745 0.96362029 0.98347253 1.00000000
#
# Loadings:
# Comp.1 Comp.2 Comp.3 Comp.4
# height 0.497 0.543 0.450 0.506
# weight 0.515 -0.210 0.462 -0.691
# chest 0.481 -0.725 -0.175 0.461
# sit-h 0.507 0.368 -0.744 -0.232
結果解讀:
Standard deviation 標準差 其平方為方差=特徵值
Proportion of Variance 方差貢獻率
Cumulative Proportion 方差累計貢獻率
由結果顯示:前兩個主成分的累計貢獻率已經達到96.36%,可以捨去另外兩個主成分,達到降維的目的。
因此可以得到函數表達式:
注意要點:
cor是邏輯變量,當cor=TRUE表示用樣本的相關矩陣R做主成分分析,當cor=FALSE表示用樣本的協方差陣S做主成分分析
loading是邏輯變量,當loading=TRUE時表示顯示loading 的內容,loadings的輸出結果為載荷是主成分對應於原始變量的係數,即Q矩陣
3.結果解讀這裡我們可以看一看得到的test.pr變量的結構:
sdev是標準偏差
center是每列計算是減去的均值
scores即降維之後的結果
我們可以利用函數來驗證下scores的結果到底是什麼意思:
library(factoextra)
# PCA結果圖
fviz_pca_ind(test.pr)
# 手動畫散點圖
ggplot(as.data.frame(test.pr$scores),aes(Comp.1,Comp.2)) + geom_point()
PCA結果圖:
手動畫散點圖:
可以看到,這兩者的結果圖是一樣的!
4.畫主成分的碎石圖並預測screeplot(test.pr,type="lines")
主要用到的函數是fviz_pca_ind,這個函數來自factoextraR包,所以需要先安裝&加載才可使用,下面記錄下關於這個函數最常用的幾個選項:
Usage
fviz_pca_ind(X, axes = c(1, 2), geom = c("point", "text"),
geom.ind = geom, repel = FALSE, habillage = "none", palette = NULL,
addEllipses = FALSE, col.ind = "black", fill.ind = "white",
col.ind.sup = "blue", alpha.ind = 1, select.ind = list(name = NULL, cos2
= NULL, contrib = NULL), ...)
Arguments
# geom——指定圖形上是只顯示點,還是同時也顯示標籤。默認同時顯示。
# palette——自行指定顏色
# addEllipses——加95%置信橢圓
# col.ind——每個點的顏色
# legend.title——指定legend的名字
下面看實例:
fviz_pca_ind(test.pr,
geom.ind = "point",
col.ind = as.character(c(rep("Normal",15),rep("Tumor",15))),
palette = c("red", "black"),
addEllipses = T,
legend.title = "Groups")
是有點醜了,不過也是為了方便理解這個函數每個參數的意義。
後記稍微整理了下,感覺對PCA怎麼畫有了更多了解,雖然之前畫過,但是都是跑流程,從沒有關注具體結果,所以,看似簡單,但是卻不熟悉。