社区首页 >专栏 >R语言ggplot2给PCA散点图结果上添加水平和垂直误差线

R语言ggplot2给PCA散点图结果上添加水平和垂直误差线

用户7010445

发布于 2021-11-08 07:01:37

1.2K01

代码可运行

文章被收录于专栏：小明的数据分析笔记本小明的数据分析笔记本

运行总次数：1

代码可运行

image.png

公众号后台有读者留言问这个图的实现办法，这个图相比于普通的PCA散点图是多了一个垂直和水平的误差线，这个如何实现之前还没有尝试过，所以查了查资料，找到了一个参考链接

https://cran.r-project.org/web/packages/SIBER/vignettes/Plot-SIA-ggplot2.html

按照这个参考链接的代码我们试试。用之前提到的小麦种子数据做示例数据

加载需要用到的包

library(ggplot2)
library(ggforce)
library(tidyverse)

读取数据

df<-read.csv("Seed_Data.csv")
df$target<-paste0("cultivar",df$target)

主成分分析

df.pca<-prcomp(df[,1:7],scale. = T)

提取主成分分析的作图数据

df.pca$x %>% as.data.frame() %>% 
  mutate(group=df$target) -> pca.result

计算PCA结果的平均值和标准差

pca.result %>% 
  group_by(group) %>% 
  summarise(pc1m=mean(PC1),
            pc1sd=sd(PC1),
            pc2m=mean(PC2),
            pc2sd=mean(PC2)) -> pca.result.a

初步的草图

ggplot()+
  geom_errorbar(data=pca.result.a,
                aes(x=pc1m,
                    ymin=pc2m-1.96*pc2sd,
                    ymax=pc2m+1.96*pc2sd,
                    color=group),
                width=0.2)+
  geom_errorbarh(data = pca.result.a,
                 aes(y=pc2m,
                     xmin=pc1m-1.96*pc1sd,
                     xmax=pc1m+1.96*pc1sd,
                     color=group))+
  stat_ellipse(data=pca.result,
               geom="polygon",
               aes(x=PC1,
                   y=PC2,
                   color=group,
                   fill=group),
               alpha=0.2) -> p
p

image.png

从图中提取作图数据

ggplot_build(p)$data[[1]] %>% 
  select(colour,x,ymin,ymax) %>% 
  pivot_longer(cols = c(ymin,ymax)) %>% 
  rename("group"="colour",
         "x" ="x",
         "y"="value") %>% 
  select(x,y,group) -> tmp.1

ggplot_build(p)$data[[2]] %>% 
  select(colour,y,xmin,xmax) %>% 
  pivot_longer(cols = c(xmin,xmax))%>% 
  rename("group"="colour",
         "y" ="y",
         "x"="value") %>% 
  select(x,y,group) -> tmp.2

tmp<-rbind(tmp.1,tmp.2)

最终作图

ggplot()+
  geom_mark_ellipse(data=tmp,
                    aes(x=x,y=y,
                        fill=group),
                    expand = unit(0,"mm"),
                    color="white")+
  geom_errorbar(data=pca.result.a,
                aes(x=pc1m,
                    ymin=pc2m-1.96*pc2sd,
                    ymax=pc2m+1.96*pc2sd),
                width=0.2)+
  geom_errorbarh(data = pca.result.a,
                 aes(y=pc2m,
                     xmin=pc1m-1.96*pc1sd,
                     xmax=pc1m+1.96*pc1sd))+
  theme_bw()+
  labs(x="PC1 (42.8%)",y="PC2(37.9%)")+
  geom_point(data=pca.result.a,
             aes(x=pc1m,y=pc2m),
             size=3)