詳解R語言中的PCA分析與視覺化

2021-03-10 16:02:24

1. 常用術語

（1）標準化（Scale）

如果不對資料進行scale處理，本身數值大的基因對主成分的貢獻會大。如果關注的是變數的相對大小對樣品分類的貢獻，則應SCALE，以防數值高的變數匯入的大方差引入的偏見。但是定標(scale)可能會有一些負面效果，因為定標後變數之間的權重就是變得相同。如果我們的變數中有噪音的話，我們就在無形中把噪音和資訊的權重變得相同，但PCA本身無法區分訊號和噪音。在這樣的情形下，我們就不必做定標。

（2）特徵值 (eigen value)

特徵值與特徵向量均為矩陣分解的結果。特徵值表示標量部分，一般為某個主成分的方差，其相對比例可理解為方差解釋度或貢獻度；特徵值從第一主成分會逐漸減小。

（3）特徵向量(eigen vector)

特徵向量為對應主成分的線性轉換向量（線性迴歸係數），特徵向量與原始矩陣的矩陣積為主成分得分。特徵向量是單位向量，其平方和為1。特徵向量主要起轉換作用，其數值不能說明什麼問題，解釋力更強的是載荷loadings，但很多R輸出中經常混用，engen vector與loadings。

（4）載荷(loading)

因子載荷矩陣並不是主成分的特徵向量，即不是主成分的係數。主成分系數的求法：各自因子載荷向量除以各自因子特徵值的算數平方根。

特徵向量是單位向量，特徵向量乘以特徵值的平方根構造了載荷loading。
列上看，不同變數對某一PC的loadings的平方和等於其徵值，因此每個變數的loadings值可表徵其對PC的貢獻。
行上看，同一變數對不同PCs的loadings行平方和為1，表徵不同PCs對某一變數方差的解釋度。

（5）得分(score)

指主成分得分，矩陣與特徵向量的積。·

2. PCA分析過程

2.0 手動計算

#特徵分解
dat_eigen<-scale(iris[,-5],scale=T)%>%cor()%>%eigen()
#特徵值提取
dat_eigen$values 
#特徵向量提取
dat_eigen$vectors
#求loadings=eigen vector*sqrt(eigen value)，與princomp不同
#主成分載荷表示各個主成分與原始變數的相關係數。
sweep(dat_eigen$vectors,2,sqrt(dat_eigen$values),"*")
#將中心化的變數矩陣得到每個觀測值的得分
scale(iris[,-5],scale=T)%*%dat_eigen$vectors%>%head()

2.1 prcomp函數

prcomp函數使用較為簡單，但是不同於常規的求取特徵值和特徵向量的方法，prcomp函數是對變數矩陣(相關矩陣)採用SVD方法計算其奇異值（原理上是特徵值的平方根），函數幫助中描述為函數結果中的sdev。
prcomp函數輸入引數為變數矩陣（x）,中心化（center,預設為true），標準化（scale，預設為false，建議改為true），主成份個數（rank）。
prcomp函數輸出有sdev(各主成份的奇異值)，rotation（特徵向量，迴歸係數），x（score得分矩陣）。

iris.pca<-prcomp(iris[,-5],scale=T,rank=4,retx=T) #相關矩陣分解
	#retx表四返回score，scale表示要標準化
	summary(iris.pca) #方差解釋度
	iris.pca$sdev #特徵值的開方
	iris.pca$rotation #特徵向量，迴歸係數
	iris.pca$x #樣本得分score

2.2 princomp函數

princomp以計算相關矩陣或者協方差矩陣的特徵值為主要手段。
princomp函數輸出有主成份的sd,loading,score,center,scale.
prcomp函數使用較為簡單，但是不同於常規的求取特徵值和特徵向量的方法，prcomp函數是對變數矩陣(相關矩陣)採用SVD方法計算其奇異值（原理上是特徵值的平方根），函數幫助中描述為函數結果中的sdev。
prcomp函數輸入引數為變數矩陣（x）,中心化（center,預設為true），標準化（scale，預設為false，建議改為true），主成份個數（rank）。
prcomp函數輸出有sdev(各主成份的奇異值及其方差累積)，rotation（載荷矩陣），x（得分矩陣），center(變數的均值),scale（變數的標準偏差）

data(wine) #三種葡萄釀造的紅酒品質分析資料集
wine.pca<-princomp(wine,cor=T,scores=T) 
#預設方差矩陣(cor=F),改為cor=T則結果與prcomp相同
summary(wine.pca) #各主成份的SVD值以及相對方差
wine.pca$loading #特徵向量，迴歸係數
wine.pca$score
screenplot(wine.pca) #方差分佈圖
biplot(wine.pca,scale=F) #碎石圖,直接把x與rotation繪圖，而不標準化

2.3 psych::principal

實際上該principal主要用於因子分析。

model_pca<-psych::principal(iris[,-5],nfactors=4,rotate="none")
model_pca$values # 特徵值=sdev^2
# 此處loadings與手動計算相同=prcomp的rotation*sdev 
model_pca%>%.$loadings #載荷，不是特徵向量
#此處score=prcomp的score/sdev
model_pca$scores[1:5,] #此處為因子得分，不是主成分得分
model_pca$weights #此處為上面loadings/特徵值，也稱成份得分系數或者因子係數

3. PCA結果解釋

下文參照chentong的內容

prcomp函數會返回主成分的標準差、特徵向量和主成分構成的新矩陣。
不同主成分對資料差異的貢獻和主成分與原始變數的關係。
1. 主成分的平方為為特徵值,其含義為每個主成分可以解釋的資料差異，計算方式為 eigenvalues = (pca$sdev)^2
2. 每個主成分可以解釋的資料差異的比例為 percent_var = eigenvalues*100/sum(eigenvalues)
3. 可以使用summary(pca)獲取以上兩條資訊。

這兩個資訊可以判斷主成分分析的質量：
 成功的降維需要保證在前幾個為數不多的主成分對資料差異的解釋可以達到80-90%。

指導選擇主成分的數目：
 1. 選擇的主成分足以解釋的總方差大於80% (方差比例碎石圖)
 2. 從前面的協方差矩陣可以看到，自動定標(scale)的變數的方差為1 (協方差矩陣對角線的值)。待選擇的主成分應該是那些方差大於1的主成分，即其解釋的方差大於原始變數（特徵值碎石圖，方差大於1，特徵值也會大於1，反之亦然）。

鑑定核心變數和變數間的隱性關係:
 原始變數與主成分的相關性 Variable correlation with PCs (var.cor) = loadings * sdev
 原始資料對主成分的貢獻度 var.cor^2 / (total var.cor^2)

4. PCA結果視覺化

4.1 ggbiplot包

devtools::install_github("vqv/ggbiplot")
library(ggbiplot)
ggscreeplot(wine.pca) #碎石圖

碎石圖

biplot

ggbiplot(wine.pca, obs.scale = 1, var.scale = 1,
   groups = wine.class, ellipse = TRUE, circle = TRUE) +
 scale_color_discrete(name = '') +
 theme(legend.direction = 'horizontal', legend.position = 'top')

4.2 ggfortify包

ggfortify包中autoplot函數可自動繪製。

library(ggfortify)
pca1<-iris%>%select(-5)%>%prcomp()
autoplot(pca1,data = iris,col= 'Species',size=2,
   loadings =T,loadings.label = TRUE,
   frame = TRUE,frame.type='norm',
   label = TRUE, label.size = 3
)+ theme_classic()

4.3 factoextra包視覺化

FactoMineR與factoextra分別進行PCA分析與視覺化，當然factoextra包中函數也可對prcomp、princomp函數結果進行視覺化。

library(factoextra)
library(FactoMineR)
# 利用FactoMineR包中PCA函數進行PCA分析
> wine.pca2 <- PCA(wine,scale.unit = T,ncp=5,graph = T) #

wine.pca2中內容

> print(wine.pca2)
**Results for the Principal Component Analysis (PCA)**
The analysis was performed on 178 individuals, described by 13 variables
*The results are available in the following objects:

 name    description       
1 "$eig"    "eigenvalues"      
2 "$var"    "results for the variables"   
3 "$var$coord"  "coord. for the variables"   
4 "$var$cor"   "correlations variables - dimensions"
5 "$var$cos2"  "cos2 for the variables"    
6 "$var$contrib"  "contributions of the variables"  
7 "$ind"    "results for the individuals"  
8 "$ind$coord"  "coord. for the individuals"   
9 "$ind$cos2"  "cos2 for the individuals"   
10 "$ind$contrib"  "contributions of the individuals" 
11 "$call"   "summary statistics"     
12 "$call$centre"  "mean of the variables"    
13 "$call$ecart.type" "standard error of the variables" 
14 "$call$row.w"  "weights for the individuals"  
15 "$call$col.w"  "weights for the variables"

摘要資訊

> summary(wine.pca2)

Call:
PCA(X = wine, scale.unit = T, ncp = 5, graph = T) 


Eigenvalues
      Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6 Dim.7 Dim.8 Dim.9 Dim.10
Variance    4.706 2.497 1.446 0.919 0.853 0.642 0.551 0.348 0.289 0.251
% of var.    36.199 19.207 11.124 7.069 6.563 4.936 4.239 2.681 2.222 1.930
Cumulative % of var. 36.199 55.406 66.530 73.599 80.162 85.098 89.337 92.018 94.240 96.170
      Dim.11 Dim.12 Dim.13
Variance    0.226 0.169 0.103
% of var.    1.737 1.298 0.795
Cumulative % of var. 97.907 99.205 100.000

...省略...

還輸出了簡易的圖

4.3.1 特徵值視覺化提取特徵值

> get_eigenvalue(wine.pca2) #標準化資料中特徵值>1的變數解釋能力較強
  eigenvalue variance.percent cumulative.variance.percent
Dim.1 4.7058503  36.1988481     36.19885
Dim.2 2.4969737  19.2074903     55.40634
Dim.3 1.4460720  11.1236305     66.52997
...省略部分

碎石圖

fviz_eig(wine.pca2,addlabels = TRUE) #碎石圖,展示方差解釋度

4.3.2 變數資訊視覺化

變數提取主要有get_pca_var()函數,輸出變數在主成分投影上的座標，變數與主成分PC的相關係數，相關係數的平方，變數對某一PC的相關貢獻

#get_pca_var()等同於get_pca(element="var")
> get_pca_var(wine.pca2)#coord cor cos2 contribution
Principal Component Analysis Results for variables
 ===================================================
 Name  Description         
1 "$coord" "Coordinates for the variables"    
2 "$cor"  "Correlations between variables and dimensions"
3 "$cos2" "Cos2 for the variables"      
4 "$contrib" "contributions of the variables" 
> wine.pca2$var #輸出上述資料
> get_pca_var(wine.pca2)$coord
> get_pca_var(wine.pca2)$cos2

變數座標(coord)與相關性(cor)視覺化

coord是座標(實際的loading)，與cor數值相同
coord=eigen vector * stdev
相關圖中，靠近的變數表示正相關；對向的是負相關。
箭頭越遠離遠原點、越靠經圓圈表明PC對其的代表性高（相關性強）

fviz_pca_var(wine.pca2) #變數相關性視覺化圖

cos2視覺化

cos2代表不同主成分對變數的代表性強弱，對特定變數，所有組成份上的cos2之和為1，因為cos2為cor的平方，所以也認為是相關性越強，其結果與cor類似。

#cos2是coord的平方，表徵特定變數在所有PC上的代表性，某個變數的所有cos2總和為1
library("corrplot")
corrplot(get_pca_var(wine.pca2)$cos2, is.corr=FALSE)

下圖中PC1對Phenols變數的代表性最好

fviz_cos2(wine.pca2, choice = "var", axes = 1:2) 
# cos2在主成分上的加和，並排序

#不同按照cos2大小設定顏色梯度，也可以設定alpha梯度
fviz_pca_var(wine.pca2,axes=c(1,2),
 col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
    repel = TRUE) # Avoid text overlapping

contrib視覺化

contrib是每個變數對某一特定PC的貢獻
contrib=(var.cos2 * 100) / (total cos2 of the PC component)
多個變數的contrib = [(Contrb1 * Eig1) + (Contrib2 * Eig2)]/(Eig1 + Eig2)

> get_pca_var(wine.pca2)$contrib
      Dim.1  Dim.2  Dim.3  Dim.4  Dim.5
Alcohol  2.083097e+00 23.391881971 4.3007553 0.03188475 7.0577176
MalicAcid  6.011695e+00 5.059392535 0.7923294 28.82511766 0.1240000
Ash   4.206853e-04 9.989949520 39.2156374 4.58711714 2.0456286
AlcAsh   5.727426e+00 0.011215874 37.4642355 0.37038683 0.4369599
Mg    2.016174e+00 8.978053590 1.7097376 12.37608338 52.8599543
Phenols  1.557572e+01 0.423013810 2.1368289 3.92310704 2.2295987
Flav   1.788734e+01 0.001128834 2.2705035 2.31937045 1.1886633
NonFlavPhenols 8.912201e+00 0.082825894 2.9025311 4.13313047 25.0703477
Proa   9.823804e+00 0.154462537 2.2336591 15.92461164 1.8730611
Color   7.852920e-01 28.089541241 1.8852996 0.43461955 0.5842581
Hue   8.803953e+00 7.797226784 0.7262776 18.29883787 3.0142002
OD    1.415019e+01 2.705899746 2.7557523 3.39004479 1.0233546
Proline  8.222684e+00 13.315407665 1.6064528 5.38568842 2.4922558
fviz_contrib(wine.pca2, choice = "var", axes = 1:2)

corrplot(get_pca_var(wine.pca2)$contrib, is.corr=FALSE)

fviz_contrib(wine.pca2, choice = "var", axes = 1:2)

根據contribution將變數顏色分類

fviz_pca_var(wine.pca2,col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

變數分組

#人為分組
bb<-as.factor(c(rep(c("soil","micro","plant"),4),"soil"))
names(bb)<-row.names(wine.pca2$var$contrib)
fviz_pca_var(wine.pca2, col.var = bb, palette = c("#0073C2FF", "#EFC000FF", "#868686FF"),
    legend.title = "Cluster")

4.3.3 樣本視覺化scores 樣本座標視覺化

get_pca_ind(wine.pca2) #coord cor cos2 contribution
get_pca(element="ind)

get_pca_ind(wine.pca2) #coord cor cos2 contribution
wine.pca2$ind #coord cos2 contrib dist
fviz_pca_ind(wine.pca2)#score 視覺化coord

fviz_pca_ind(wine.pca2, geom=c("point","text"),
   addEllipses = T,
   pointshape=21,col.ind="black",pointsize="cos2",
   fill.ind = wine.class,palette = "npg",
   #col.ind="cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
   #col.ind="contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
   # label=wine.class,
   repel = TRUE)

fviz_pca_ind(wine.pca2, axes = c(1, 2),label="none", #habillage只能是分類變數
    addEllipses = TRUE, ellipse.type="norm",ellipse.level=0.9,
    # one of "confidence","t","norm","euclid","convex"
    habillage = wine.class,palette = "jco",
    mean.point=F
    )

樣本的cos2與contrib圖

fviz_cos2(wine.pca2, choice = "ind",axes=1:2)

fviz_contrib(wine.pca2, choice = "ind", axes = 1:2)

4.3.4 biplot

biplot不需要關注具體數值，只需要關注方向與位置
樣本在變數同側是具有高數值，反之則值低

fviz_pca_biplot(wine.pca2, axes = c(1,2),repel = F,
    addEllipses = T,ellipse.alpha=0.15,
    geom=c("point"),geom.var=c("arrow","text"),
    arrowsize=1.5,labelsize=5, #arrow與text大小
    pointshape=21,pointsize=5,fill.ind = wine.class,col.ind="black", #point
    col.var=factor(c(rep(c("soil","plant"),6),"plant"))
    
)%>%ggpar(xlab="PC1",ylab="PC2",title="PCA-Biplot",
   font.x=14,font.y=14,font.tickslab = 14,ggtheme=theme_bw(),ylim=c(-4.5,4.5),
   legend.title = list(color="Variable",fill="Class"),font.legend = 12,
   )+
  ggpubr::fill_palette("jco")+ggpubr::color_palette("npg")+
 theme(axis.ticks.length= unit(-0.25, 'cm'), #設定y軸的刻度長度為負數，即向內
  axis.text.y.left = element_text(margin = unit(c(0.5, 0.5, 0.5, 0.05), 'cm')),
  axis.text.x.bottom = element_text(margin = unit(c(0.5, 0.5, 0.05, 0.5), 'cm'))
  )

到此這篇關於詳解R語言中的PCA分析與視覺化的文章就介紹到這了,更多相關R語言PCA分析與視覺化內容請搜尋it145.com以前的文章或繼續瀏覽下面的相關文章希望大家以後多多支援it145.com！