【R画图学习8.1】PCA图

主成分分析，即Principle Component Analysis (PCA)，是一种传统的统计学方法，被机器学习领域引入后，通常被认为是一种特殊的非监督学习算法，其可以对复杂或多变量的数据做预处理，以减少次要变量，便于进一步使用精简后的主要变量进行数学建模和统计学模型的训练，所以PCA又被称为主变量分析。

R语言中，实现这个方法的函数主要有 princomp ， prcomp 等。查资料后，这两种方法的不同主要体现如下：

具体细节大家可以自行查资料学习。

# 使用R自带iris数据集（150*5，行为样本，列为特征） 可以自己对比一下2个函数结果的异同

############ prcomp函数 ################

# retx表示返回score，scale表示要标准化

iris.pca <- prcomp(iris[,-5],scale=T,rank=4,retx=T) #相关矩阵分解

summary(iris.pca) #方差解释度

iris.pca$sdev #特征值的开方

iris.pca$rotation #特征向量，回归系数

iris.pca$x #样本得分score

############ princomp函数 ################

library(stats)

# 默认方差矩阵(cor=F),改为cor=T则结果与prcomp相同

iris.pca<-princomp(iris[,-5],cor=T,scores=T)

# 各主成份的SVD值以及相对方差

summary(iris.pca)

# 特征向量，回归系数

iris.pca$loading

iris.pca$score

下面我们就用R这个自带的iris数据来画一下PCA图。

从这个数据来看，包含150行，5列，最后一列是按species分成了三类。

我们用prcomp函数来计算各个主成分，但是因为最后一列是Species，所以计算的时候需要先排除。

pca1 <- prcomp(iris[,-ncol(iris)],center = TRUE,scale= TRUE)

可以看出PC的信息存在x中。

df1 <- pca1$x # 提取PC score

df1 <- as.data.frame(df1) # 注意：如果不转成数据框形式后续绘图时会报错

可以通过summary函数查看每个主成分的贡献度。

summ1 <- summary(pca1)

summ1

所以我们可以summary里面的信息把PC1 PC2提取出来作为X和Y轴的lab。summ1$importance[2,1]和summ1$importance[2,2]是PC1和PC2。round函数是取小数点的位数。

xlab1 <- paste0("PC1(",round(summ1$importance[2,1]*100,2),"%)")

ylab1 <- paste0("PC2(",round(summ1$importance[2,2]*100,2),"%)")

df1$Species <- iris_input$Species

把原先的Specis类别信息加进来。

首先还是用ggplot2的来画，先画个基本版本的二维PCA图。

ggplot(data = df1,aes(x = PC1,y = PC2,color = Species))+

geom_point(size = 2.5)

下面我们添加上椭圆形的置信区间。

ggplot(data = df1,aes(x = PC1,y = PC2,color = Species))+

  geom_point(size = 2.5)+

  # 绘制置信椭圆：

  stat_ellipse(aes(fill = Species),type = "norm",geom = "polygon",alpha = 0.25,color = NA)

添加上x，y的lab和title。

ggplot(data = df1,aes(x = PC1,y = PC2,color = Species))+

  geom_point(size = 2.5)+

  # 绘制置信椭圆：

  stat_ellipse(aes(fill = Species),type = "norm",geom = "polygon",alpha = 0.25,color = NA)+

  labs(x = xlab1,y = ylab1,title = "PCA Scores Plot")

换个背景模式。

ggplot(data = df1,aes(x = PC1,y = PC2,color = Species))+

  geom_point(size = 2.5)+

  # 绘制置信椭圆：

  stat_ellipse(aes(fill = Species),type = "norm",geom = "polygon",alpha = 0.25,color = NA)+

  labs(x = xlab1,y = ylab1,title = "PCA Scores Plot")+

  theme_bw()

下面我们尝试修改色系。

ggplot(data = df1,aes(x = PC1,y = PC2,color = Species))+

  geom_point(size = 2.5)+

  # 绘制置信椭圆：

  stat_ellipse(aes(fill = Species),type = "norm",geom = "polygon",alpha = 0.25,color = NA)+

  labs(x = xlab1,y = ylab1,title = "PCA Scores Plot")+

  theme_bw()+

  scale_colour_manual(values = c("purple","orange","pink"))   #改变点的颜色

ggplot(data = df1,aes(x = PC1,y = PC2,color = Species))+

  geom_point(size = 2.5)+

  # 绘制置信椭圆：

  stat_ellipse(aes(fill = Species),type = "norm",geom = "polygon",alpha = 0.25,color = NA)+

  labs(x = xlab1,y = ylab1,title = "PCA Scores Plot")+

  theme_bw()+

  scale_colour_manual(values = c("purple","orange","pink"))+

  scale_fill_manual(values = c("purple","orange","pink"))   #可以改变填充色

下面我们试着调整文字等外观。

ggplot(data = df1,aes(x = PC1,y = PC2,color = Species))+

  geom_point(size = 2.5)+

  # 绘制置信椭圆：

  stat_ellipse(aes(fill = Species),type = "norm",geom = "polygon",alpha = 0.25,color = NA)+

  labs(x = xlab1,y = ylab1,title = "PCA Scores Plot")+

  theme_bw()+

  scale_colour_manual(values = c("purple","orange","pink"))+

  scale_fill_manual(values = c("purple","orange","pink"))+

  theme(plot.title = element_text(hjust = 0.5,size = 15), #修改title的位置

        axis.text = element_text(size = 11),axis.title = element_text(size = 13), #分别修改text和title文字的大小

        legend.text = element_text(size = 11),legend.title = element_text(size = 13),#分别修改legend文字的大小

        plot.margin = unit(c(0.4,0.4,0.4,0.4),'cm'))

下面，我们再测试另外一款经常用来画PCA的包ggbiplot，这个我经常在宏基因组领域看见，其实也是基于ggplot之上开发的包，只不过替用户省了很多事。

library(devtools)

install_github("vqv/ggbiplot")

library(ggbiplot)

另外还可以把ggplot上面的有些属性望上面叠加的。

ggbiplot(pca1, obs.scale = 1, var.scale = 1,

        #是否显示与变量的箭头

        var.axes=1,

        #改变箭头对应的变量的字体大小

        varname.size=4,

        # 分组信息：

        groups = iris$Species,

        # 是否显示置信椭圆：

        ellipse = T,

        # 是否显示中心的圆：

        circle = F,

        labels.size = 5) +

  theme_bw()+

  theme(legend.direction = 'horizontal', legend.position = 'top')+

  scale_color_discrete(name = 'Species') +

  theme(axis.text = element_text(size = 15),axis.title = element_text(size = 15),

        legend.text = element_text(size = 15),legend.title = element_text(size = 15))

其中箭头代表原始变量，方向代表原始变量和主成分的相关性，长度代表原始变量对于主成分的贡献度。

当然，还有很多其它的包也可以实现简便的PCA作图。

library(factoextra)

library(FactoMineR)

iris.pca <- PCA(iris[,-5], graph = T)

PC1,2...柱状图

fviz_screeplot(iris.pca, addlabels = TRUE, ylim = c(0, 80))

fviz_pca_biplot(iris.pca, label = "var",

                habillage=iris$Species,

                addEllipses=TRUE,

                ellipse.level=0.95,

                ggtheme = theme_minimal())