详解R语言中的PCA分析与可视化

2023-07-04 11:12:03 227

1.常用术语

（1）标准化（Scale）

如果不对数据进行scale处理，本身数值大的基因对主成分的贡献会大。如果关注的是变量的相对大小对样品分类的贡献，则应SCALE，以防数值高的变量导入的大方差引入的偏见。但是定标(scale)可能会有一些负面效果，因为定标后变量之间的权重就是变得相同。如果我们的变量中有噪音的话，我们就在无形中把噪音和信息的权重变得相同，但PCA本身无法区分信号和噪音。在这样的情形下，我们就不必做定标。

（2）特征值(eigenvalue)

特征值与特征向量均为矩阵分解的结果。特征值表示标量部分，一般为某个主成分的方差，其相对比例可理解为方差解释度或贡献度；特征值从第一主成分会逐渐减小。

（3）特征向量(eigenvector)

特征向量为对应主成分的线性转换向量（线性回归系数），特征向量与原始矩阵的矩阵积为主成分得分。特征向量是单位向量，其平方和为1。特征向量主要起转换作用，其数值不能说明什么问题，解释力更强的是载荷loadings，但很多R输出中经常混用，engenvector与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=eigenvector*sqrt(eigenvalue)，与princomp不同
#主成分载荷表示各个主成分与原始变量的相关系数。
sweep(dat_eigen$vectors,2,sqrt(dat_eigen$values),"*")
#将中心化的变量矩阵得到每个观测值的得分
scale(iris[,-5],scale=T)%*%dat_eigen$vectors%>%head()

2.1prcomp函数

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.2princomp函数

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.3psych::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，反之亦然）。

鉴定核心变量和变量间的隐性关系:
原始变量与主成分的相关性VariablecorrelationwithPCs(var.cor)=loadings*sdev
原始数据对主成分的贡献度var.cor^2/(totalvar.cor^2)

4.PCA结果可视化

4.1ggbiplot包

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.2ggfortify包

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.3factoextra包可视化

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)
**ResultsforthePrincipalComponentAnalysis(PCA)**
Theanalysiswasperformedon178individuals,describedby13variables
*Theresultsareavailableinthefollowingobjects:

namedescription
1"$eig""eigenvalues"
2"$var""resultsforthevariables"
3"$var$coord""coord.forthevariables"
4"$var$cor""correlationsvariables-dimensions"
5"$var$cos2""cos2forthevariables"
6"$var$contrib""contributionsofthevariables"
7"$ind""resultsfortheindividuals"
8"$ind$coord""coord.fortheindividuals"
9"$ind$cos2""cos2fortheindividuals"
10"$ind$contrib""contributionsoftheindividuals"
11"$call""summarystatistics"
12"$call$centre""meanofthevariables"
13"$call$ecart.type""standarderrorofthevariables"
14"$call$row.w""weightsfortheindividuals"
15"$call$col.w""weightsforthevariables"

摘要信息

>summary(wine.pca2)

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


Eigenvalues
Dim.1Dim.2Dim.3Dim.4Dim.5Dim.6Dim.7Dim.8Dim.9Dim.10
Variance4.7062.4971.4460.9190.8530.6420.5510.3480.2890.251
%ofvar.36.19919.20711.1247.0696.5634.9364.2392.6812.2221.930
Cumulative%ofvar.36.19955.40666.53073.59980.16285.09889.33792.01894.24096.170
Dim.11Dim.12Dim.13
Variance0.2260.1690.103
%ofvar.1.7371.2980.795
Cumulative%ofvar.97.90799.205100.000

...省略...

还输出了简易的图

4.3.1特征值可视化提取特征值

>get_eigenvalue(wine.pca2)#标准化数据中特征值>1的变量解释能力较强
eigenvaluevariance.percentcumulative.variance.percent
Dim.14.705850336.198848136.19885
Dim.22.496973719.207490355.40634
Dim.31.446072011.123630566.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)#coordcorcos2contribution
PrincipalComponentAnalysisResultsforvariables
===================================================
NameDescription
1"$coord""Coordinatesforthevariables"
2"$cor""Correlationsbetweenvariablesanddimensions"
3"$cos2""Cos2forthevariables"
4"$contrib""contributionsofthevariables"
>wine.pca2$var#输出上述数据
>get_pca_var(wine.pca2)$coord
>get_pca_var(wine.pca2)$cos2

变量坐标(coord)与相关性(cor)可视化

coord是坐标(实际的loading)，与cor数值相同
coord=eigenvector*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)#Avoidtextoverlapping

contrib可视化

contrib是每个变量对某一特定PC的贡献
contrib=(var.cos2*100)/(totalcos2ofthePCcomponent)
多个变量的contrib=[(Contrb1*Eig1)+(Contrib2*Eig2)]/(Eig1+Eig2)

>get_pca_var(wine.pca2)$contrib
Dim.1Dim.2Dim.3Dim.4Dim.5
Alcohol2.083097e+0023.3918819714.30075530.031884757.0577176
MalicAcid6.011695e+005.0593925350.792329428.825117660.1240000
Ash4.206853e-049.98994952039.21563744.587117142.0456286
AlcAsh5.727426e+000.01121587437.46423550.370386830.4369599
Mg2.016174e+008.9780535901.709737612.3760833852.8599543
Phenols1.557572e+010.4230138102.13682893.923107042.2295987
Flav1.788734e+010.0011288342.27050352.319370451.1886633
NonFlavPhenols8.912201e+000.0828258942.90253114.1331304725.0703477
Proa9.823804e+000.1544625372.233659115.924611641.8730611
Color7.852920e-0128.0895412411.88529960.434619550.5842581
Hue8.803953e+007.7972267840.726277618.298837873.0142002
OD1.415019e+012.7058997462.75575233.390044791.0233546
Proline8.222684e+0013.3154076651.60645285.385688422.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)#coordcorcos2contribution
get_pca(element="ind)

get_pca_ind(wine.pca2)#coordcorcos2contribution
wine.pca2$ind#coordcos2contribdist
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,
#oneof"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.4biplot

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分析与可视化内容请搜索毛票票以前的文章或继续浏览下面的相关文章希望大家以后多多支持毛票票！

声明：本文内容来源于网络，版权归原作者所有，内容由互联网用户自发贡献自行上传，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任。如果您发现有涉嫌版权的内容，欢迎发送邮件至：czq8825#qq.com（发邮件时，请将#更换为@）进行举报，并提供相关证据，一经查实，本站将立刻删除涉嫌侵权内容。