虽然PCA和RDA分析及绘图都写过教程，但是对于结果的解释都没有写的很详细，刚好最近有人询问怎样使用FactoMineR factoextra包进行PCA分析。所以使用R统计绘图-环境因子相关性热图中的不同土壤环境因子数据进行PCA绘图和结果解读推文。
一、 PCA分析过程
1.1 数据准备
# 1.1.1 设置工作路径
setwd("D:\\EnvStat\\PCA\\3DPCA")
getwd()#获取工作路径
# 1.1.2 调用R包
library(factoextra)
library(scales)
library(ggsci)
library(pca3d)
library(rgl)
# 1.1.3 导入数据
data = read.csv("env.csv",
header = TRUE,row.names=1,
stringsAsFactors = FALSE,
check.names = FALSE)
##设置check.names=FALSE，读入数据是不会更改行名和列名中的特殊符号(-、空格等符号)
head(data)
dim(data) # 1-3列为分组信息，4-14列为分析变量信息。
图1|环境因子及分组信息表，env.csv。行为样品名称，列为环境因子名称和分组信息，共有11个环境变量，3个分组信息。
1.2 主成分分析(PCA)
使用R中能进行PCA分析的R函数(vegan::rda(),prcomp()或princomp)就行，pca3D需要prcomp()输出结果，所以这里使用prcomp()进行PCA分析。
## 1.2.1 使用prcomp()进行PCA分析
#library(vegan)
#pca=rda(data[4:14],scale=TRUE)
#pca = prcomp(data[4:14], scale.= TRUE)
pca <- princomp(data[4:14],cor=TRUE,scores=TRUE) # cor = TRUE使用相关性矩阵进行相关性分析。
## 1.2.2 提取样本PC1-3坐标数据
ind = data.frame(get_pca_ind(pca)$coord[,1:3])
colnames(ind)= c("PC1","PC2","PC3")
ind = cbind(ind,data[2:3])
ind
####设置因子水平---基于tillage和depth分组信息绘图
ind$tillage = factor(ind$tillage,levels = c("CK","NT","PT","WL"))
ind$depth = factor(ind$depth,levels = c("0-10 cm","10-20 cm","20-30 cm"))
## 1.2.3 提取变量PC1-PC3数据
var = data.frame(get_pca_var(pca)$coord[,1:3])
colnames(var)= c("PC1","PC2","PC3")
var注：这里只展示绘图过程，PCA分析的详细分析及结果解释流程都写在R统计绘图-PCA分析绘图及结果解读(误差线，多边形，双Y轴图、球形检验、KMO和变量筛选等)中。
图2|样本坐标数据及分组信息表,ind。
图3|变量坐标数据表,var。
二、绘制3D PCA图
3D PCA样本图其实就是3D散点图，如果要绘制3D PCA双序图，可以以箭头等其他形式展示变量。绘图原理都是差不多的，本质上都需要提供坐标信息，向图中添加各种需要的元素并设置图中的各种参数美化图形。可以绘制3D图形的R包挺多的，我觉得rgl是最强大的，其他很多可以绘制3D图的函数，很多都是调用的rgl包。另外还有一个pca3d专门用于绘制3D PCA图。这里主要介绍使用rgl或pca3d绘制3D PCA图以及输出pdf格式图。
2.1 pca3d绘制3D PCA图
样本直接用散点显示(样本名称太长，重叠情况严重，所以不添加样本标签)，使用分组信息分别着色。变量使用红色箭头显示。
## 2.1.1 准备绘图颜色
###数据分析时，最好统一绘图结果的色系，所以一般同一篇文章的数据分析结果，我都先根据需要的颜色数量，选择ggsci中的期刊颜色。之后都统一使用这一颜色集。这里选择d3颜色集
cols = pal_d3("category10")(10)
show_col(cols) # 查看颜色的序号、十六进制名称和颜色。可以自己手动选择自己想要的颜色。
### 如果有想要使用的RGB颜色，也可以使用rgb()函数转为十六进制颜色值
### col2rgb()可以将颜色名称或十六进制颜色转换为RGB值。
### grDevices::convertColor() # 可以更灵活的在各个颜色空间进行转换。
### 例：RGB(80,196,143);RGB(54.133,254);SD:RGB(245,97,111)
#col2rgb('blue', alpha=T) # 颜色名称转为RGB值。
#col2tgb("hexadecimal",alpha = TRUE) # 16进制转为RGB值。
#cols = c(rgb(red=80, green=196, blue=143, alpha=255, max=255),
#rgb(red=54, green=133, blue=254, alpha=255, max=255),
#rgb(red=245, green=97, blue=111, alpha=255, max=255))
## 2.1.2 pca3d绘制3D PCA图
#browseVignettes("pca3d")# 查看帮助手册
#??pca3d
### 默认参数绘图
pca3d(pca) # 图5
### 添加分组信息和自定义形状和颜色
listShapes() # 可以查看pca2d()和pca3d可以使用的形状
pca3d(pca, # 图6
components = 1:3, # 选择PC1-PC3用于绘图
group=ind[,4],# 选择的分组信息,不能使用ind[4]会报“Error in matchShapes(shape) : Incorrect shapes: NA”错误
palette=cols,# 设置色板
shape ="s", # 需要同时设置形状，不然默认根据分组信息分配不同的形状(tetrahedrons, cubes, spheres和octahedrons)。
radius = 2 # 类似于cex参数，设置形状大小。
) ## 不使用分组信息，向参数palette和shape提供一个与样本数等长的颜色和性状矢量也可指定颜色和形状。
### 添加分组椭圆
pca3d(pca, # 图7
components = 1:3,
group=ind[,4],
palette=cols,
shape ="s",
show.plane = TRUE, # 在y=0处添加水平平面
show.shadows = TRUE, # 绘制棒棒图，点到y=0平面的垂直线和阴影。
show.ellipses = TRUE, # 添加分组椭圆
ellipse.ci = 0.95,# 椭圆集合区间水平，默认0.95。
show.axes = TRUE,
show.scale =TRUE, # 显示轴刻度
show.centroids=TRUE, # 显示聚类质心
show.group.labels=TRUE,# 显示分组标签，可以自定义字符，显示在质心处。
axes.color= "black",#设置轴线颜色
radius = 2,
bg = "white",# 设置图背景颜色
axe.titles = c("PC1","PC2","PC3") # 可用于更改轴标题。
) # 图6
### 绘制3D PCA双序图---样本以形状显示
pca3d(pca, # 图8
components = 1:3,
group = ind[,4],
biplot = TRUE, # 绘制双序图 ，变量以红色箭头显示。
biplot.vars = 5, # 需要展示的变量数量，默认是5。
palette = cols, # 设置色板
shape ="s",
new= TRUE, # 在新的窗口展示图，不覆盖之前的绘图窗口。
#fancy = TRUE,# 相当于设置同时展示样本标签，质心、样本阴影和分组标签
bg= "white",
show.axe.titles=TRUE,
legend="right",
show.axes = TRUE,
show.scale = TRUE,
axes.color= "black")# 设置图例位置，
### 绘制3D PCA双序图---样本以样本名展示
pca3d(pca, # 图9
components = 1:3,
biplot=TRUE,
biplot.vars=6,
group=ind[,4],
palette=cols[1:4], # 会按照分组因子水平按序分配颜色。
show.ellipses = TRUE,
show.shapes = FALSE,
show.labels=TRUE, # 显示样本标签
labels.col = c(rep(cols[1],9),
rep(cols[2],9),
rep(cols[3],9),
rep(cols[4],9)), # 设置样本标签颜色
#shape ="s",
bg= "white",
legend="right",
show.axes = TRUE,
axes.color= "black",
show.scale = FALSE,
new= TRUE
)
###绘图后可用鼠标将图调整到需要的角度，然后保存该角度到本地
snapshotPCA3d(file="pca3d.png")
图4|ggsci中d3的category10颜色集。
图5|pca3d默认参数绘图结果。用鼠标可以调整图形展示角度。
图6|添加分组信息和自定义形状和颜色。
图7|设置置信椭圆、质心、轴刻度标签等信息。
图8|3D PCA双序图。变量以红色箭头显示。
图9|3D PCA双序图。样本以字符展示。样本太多，就不适合用字符标签显示，可用于探索数据。
2.2 rgl包绘制3D PCA图
rgl是一个很强大的绘制3D图的R包。许多绘制3D图的函数都需要调用rgl包。比如car::scatter3d();compositions::plot3D等函数。使用rgl可以绘制多种类型的3D图形(hist3D(),points3d(),lines3d(),segments3d(),triangles3d()和quads3d()等),可以更方便的添加绘图元素和调整图形参数。使用说明：http://127.0.0.1:16323/library/rgl/doc/rgl.html
## 2.2.1 rgl包绘制3D PCA双序图---单个分组因子
###browseVignettes("rgl") # 查看帮助信息
## 设置颜色
cols = pal_d3("category10")(10)
col = cols[c(7,1,5,2)] # 选择四个颜色，为每个分组指定特定颜色。
get_colors <- function(groups, cols = palette()){
groups <- as.factor(groups)
ngrps <- length(levels(groups))
if(ngrps > length(cols))
cols <- rep(cols, ngrps)
color <- cols[as.numeric(groups)]
names(color) <- as.vector(groups)
return(color)
}
labels.col = get_colors(ind[,4],col)
labels.col
### 或者
ind[,4] = factor(ind[,4],levels = c("CK","NT","PT","WL")) # 设置分组的因子水平
labels.col = col[as.numeric(ind[,4])]
names(labels.col) <- as.vector(ind[,4])
labels.col
### 或者，三种设置颜色方法，效果是一致的
labels.col = c(rep(cols[7],9),
rep(cols[2],9),
rep(cols[1],9),
rep(cols[5],9)) # 根据tillage分组为样本赋予颜色，这个要根据样本的排布顺序设置颜色。
open3d() # 打开一个绘图窗口；图10.
par3d(family="serif",
cex=1.2,
font=1) # 设置图形参数，应用于后面的所有绘图对象,需要与绘图参数一起运行。
### 绘制样本散点3D图
plot3d(ind[,1:3],
type="p",# "p"表示点，"s"表示球，"l"表示线,"h"表示线段
col=labels.col,
lwd = 1,
size=12,
ticktype = "detailed",
#theta = 20,
#phi = 10, # theta和phi是极坐标
#fov = 60, # 以度为单位的视场角
#zoom = 1 # zoom表示缩放系数,这些都可以使用绘图后，使用鼠标调整。
)
## 3D图在屏幕的视图也可以用view3d()设置
#view3d(theta = 20,
#phi = 10, # theta和phi是极坐标
#fov = 90, # 以度为单位的视场角
#zoom = 1 )
### 默认图形以立方体形式展示，不设置网格线
axes3d() # 只显示x y和z 3个轴
grid3d(c("x","y","z")) #x y和z 3个平面添加网格线
### 使用arrow3d()添加变量箭头
####需要准备箭头的起始点(0,0,0)数据表
origin = data.frame(PC1= rep(0,11),
PC2= rep(0,11),
PC3= rep(0,11),
row.names = rownames(var))
#### 绘制变量箭头
####11个环境因子数量有点多，这里只在图中展示5个环境因子可以使用vegan包先筛选环境因子。
####为了图美观，环境因子坐标都*3.
for (i in 1:nrow(var[1:5,])) {
arrow3d(p0 = origin[i,],
p1 = var[i,]*3,
type = "lines", #箭头类型可选 "extrusion"，"lines","rotation","flat"
width = 3,# 箭头末端倒钩的宽度
thickness = 0.618 * width, # 箭头倒钩的厚度
theta = 0.45, # 箭头末端倒钩的角度
s = 0.2,# 箭头末端倒钩的长度
n = 2, # 箭头末端倒钩的个数，默认为3
col = "red",
plot = TRUE) # plot = FALSE则返回绘图参数，但不绘图
}
#### 添加变量标签
text3d(var[1:5,]*3,
texts=rownames(var[1:5,]),
#adj = c(1,1,1), # 调整字符标签相对于坐标的位置.(0,0)将标签放在左下角落，(0.5,0.5) #中心，(1,1)右上
offset = 0.5, # 与pos联用,设置标签与坐标点的距离。
pos = 3, # 设置标签位置。0=坐标所在位置, 1=below,2= left,3=above,4=right,5=in front of,6=behind
col="black",
#family = "serif",
cex=1.5)
#### 添加图例
legend3d("right",
levels(ind$tillage),
pch=16, # 设置形状
col=cols[1:4] # 设置颜色
)
#### 从open3d()开始运行完命令后，调整为自己想要的角度和缩放程度后，保存pdf格式图片。图例没有输出
rgl.postscript("rgl.pca1.pdf", fmt = "pdf", drawText = T)
close3d() # 关闭3D图的屏幕显示
图10|rgl包绘制3D PCA双序图---单个分组因子，rgl.pca1.pdf。
上面只使用颜色显示了tillage的分组信息，还要一个分组信息depth可以使用形状显示。下面介绍使用pch3d()根据tillage设置样本点颜色，根据depth设置样本点形状。
## 2.2.2 rgl包绘制3D PCA双序图---双分组因子
### 样本根据tillage和depth分组信息设置不同的颜色和形状
#### 设置形状：?points，可以查看可以使用的平面点形状(原点，三角等)。
#### 立体图形有：cube3d();tetrahedron3d();octahedron3d();icosahedron3d();dodecahedron3d();cuboctahedron3d()。
## 设置颜色：颜色上面已经设置好，就不再重复了。
open3d() # 图11
pch3d(ind[,1:3],
pch = rep(c(15,16,17),# ?points,可以查看可使用的shape
time = 12),# 根据depth设置不同的形状
col=labels.col, # 根据tillage设置形状颜色
cex = 0.5, # 形状相对于图形的大小
radius = 0.25 # 坐标中绘制图形大小
)
bg3d("lightgrey") # 设置图形背景颜色为灰色
axes3d(edges = "bbox",
box = TRUE, # box = FALSE，则只显示立方体三面的框线。
col = "black",# 设置轴为黑色
lwd = 2,# 设置轴线宽度
family = "serif",
font = 1,
marklen = 0.3,
marklen.rel = FALSE) # 设置刻度线宽度,marklen.rel = TREU,则刻度线宽度 = 1/marklen * axis length。
aspect3d(1, 1, 1) # 设置x,y和z 轴的长度比率为1:1:1，图形近似一个正方体。
title3d(xlab = "PC1",
ylab = "PC2",
zlab = "PC3",
#main = "3D PCA biplot", # 设置图标题
col = "black")
for (i in 1:nrow(var[1:5,])) {
arrow3d(p0 = origin[i,],
p1 = var[i,]*3,
type = "lines",
width = 3,
thickness = 0.618 * width,
theta = 0.45,
s = 0.2,
n = 2,
lwd = 2,
col = "red",
plot = TRUE
)
}
text3d(var[1:5,]*3,
texts=rownames(var[1:5,]),
offset = 0.5,
pos = 3,
col="black",
family = "serif",
cex=1.2)
## 调整好3D图形展示角度之后，输出pdf图。设置的灰色图形背景没有输出到图中。
rgl.postscript("rgl.pca2.pdf", fmt = "pdf", drawText = T)
图11|rgl包绘制3D PCA双序图---双分组因子，rgl.pca2.pdf。样本点用形状和颜色显示不同分组信息。设置刻度线长度，轴标签等图形参数。
## 2.2.3 rgl包绘制3D PCA双序图---设置xyz坐标轴在(0,0,0)相交。
### rgl实现类似pca3d()的绘图结果
open3d() # 图12
#par3d() # 3d绘图参数设置
#layout3d(matrix(1:4, 2,2)) # 3D图排版设置
###heights表示图布局中行的高度，widths表示图布局中列的宽度，matrix()设置多张图依次在屏幕上的显示位置。
#mfrow3d(2,2,byrow = TRUE,sharedMouse = TRUE) # 设置绘图参数，画布被按行分为4份。
par3d(windowRect = 50 + c( 0, 0, 640, 640)) # 设置绘图窗口大小
pch3d(ind[,1:3],
pch = rep(c(15,16,17),
time = 12),
col=labels.col,
cex = 0.5,
radius = 0.25
)
### 设置xyz轴线在(0,0,0)相交，轴线颜色设置为灰色
axis3d("x",pos = c(NA,0,0),labels = TRUE,color = "lightgrey") # pos决定了轴线在另两个轴的位置，不与level同时使用。
axis3d("y",pos = c(0,NA,0),labels = TRUE,color = "lightgrey")
axis3d("z",pos = c(0,0,NA),labels = TRUE,color = "lightgrey")
### 添加轴标签
#### 与轴线相同方向的平行线，刻度线从line = 0到line=1，刻度标签放在line=2。level表示与轴线正交的线，
title3d(xlab = "PC1",
ylab = "PC2",
zlab = "PC3",
col = "black")
### 绘制变量箭头与名称
for (i in 1:nrow(var[1:5,])) {
arrow3d(p0 = origin[i,],
p1 = var[i,]*3,
type = "lines",
width = 3,
thickness = 0.618 * width,
theta = 0.45,
s = 0.2,
n = 2,
lwd = 2,
col = "black", # 变量箭头设置为黑色
plot = TRUE
)
}
text3d(var[1:5,]*3,
texts=rownames(var[1:5,]),
offset = 0.5,
pos = 3,
col="black",
family = "serif",
cex=1.2)
## 调整好3D图形展示角度之后，输出pdf图。
#rgl.postscript("rgl.pca3.pdf", fmt = "pdf", drawText = T)
### 添加框面
bbox3d(draw_front = FALSE, # 只绘制立方体的三个边界
marklen = 0.3,
marklen.rel = FALSE,
color = "#333377", # 设置边界框颜色
emission="#333377", #emission, specular, shininess: properties for lighting calculation
specular="#3333FF", shininess=5, alpha=0.8
) # bbox3d()与aexs3d()功能有重叠，先运行的会被后运行的结果覆盖。
#highlevel(integer()) # 以rglwidget显示，运行命令前后，显示的图没什么差别。
## 调整好3D图形展示角度之后，输出pdf图。设置的狂免颜色亮度没有输出到图中。
rgl.postscript("rgl.pca3.pdf", fmt = "pdf",drawText = TRUE)
图12|rgl模仿pca3d()绘图结果，rgl.pca3.pdf。
## 2.2.4 rgl包绘制3D PCA双序图---根据分组信息添加集合椭圆/球
names(col) = c("CK","NT","PT","WL") # 因子水平与样本排列顺序不一致，在设置颜色和形状时要特别注意。样本点的颜色与集合椭圆/球的颜色保持一致。
open3d() #图13
pch3d(ind[,1:3],
pch = rep(c(15,16,17),# ?points,可以查看可使用的shape
time = 12),# 根据depth设置不同的形状
col=labels.col, # 根据tillage设置形状颜色
cex = 0.5, # 形状相对于图形的大小
radius = 0.25 # 坐标中绘制图形大小
)
axes3d(edges = "bbox",labels = TRUE,
box = FALSE, # box = TRUE，则只显示立方体六面的框线。
col = "black",# 设置轴为黑色
lwd = 2,# 设置轴线宽度
family = "serif",
font = 1,
marklen = 0.3,
marklen.rel = FALSE) # 设置刻度线宽度,marklen.rel = TREU,则刻度线宽度 = 1/marklen * axis length。
aspect3d(1, 1, 1)
### 添加轴标签
#### 与轴线相同方向的平行线，刻度线从line = 0到line=1，刻度标签放在line=2。level表示与轴线正交的线，
title3d(xlab = "PC1",
ylab = "PC2",
zlab = "PC3",
col = "black")
### 绘制变量箭头与名称
for (i in 1:nrow(var[1:5,])) {
arrow3d(p0 = origin[i,],
p1 = var[i,]*3,
type = "lines",
width = 3,
thickness = 0.618 * width,
theta = 0.45,
s = 0.2,
n = 2,
lwd = 2,
col = "black", # 变量箭头设置为黑色
plot = TRUE
)
}
text3d(var[1:5,]*3,
texts=rownames(var[1:5,]),
offset = 0.5,
pos = 3,
col="black",
family = "serif",
cex=1.2)
### 根据tillage绘制分组集合区域
for (i in 1:length(levels(ind[,4]))){
group <- levels(ind[,4])[i]
xyz <- ind[which(ind[,4] == group),]
ellips <- ellipse3d(cov(cbind(xyz[,1:3])),
centre=c(mean(xyz[,1]), xyz[,2], mean(xyz[,3])), level = 0.95)
#shade3d(ellips, col = col[i], alpha = 0.2, lit = FALSE) # 阴影椭圆形式展示集合区间
wire3d(ellips, col = col[i], lit = FALSE) # 网格椭球形式
# show group labels
texts3d(apply(xyz[1:3], 2, mean),
text = group, # 在椭圆中心添加分组标签
col= col[i], cex = 2)
}
## 调整好3D图形展示角度之后，输出pdf图。设置的框面颜色亮度没有输出到图中。
rgl.postscript("rgl.pca4.pdf", fmt = "pdf",drawText = TRUE)
## 添加单个集合区间
### 获得一个mesh3d对象，用于后续绘图。
ellips <- ellipse3d(cov(ind[,1:3]), # 相关性或协方差矩阵
centre=c(mean(ind[,1]), mean(ind[,2]), mean(ind[,3])),# 椭圆中心
level = 0.95) # level设置集合区间水平，用于控制椭球的大小。
### 绘制单个集合区间，阴影或网格形式
shade3d(ellips, col = "#D95F02", alpha = 0.1, lit = FALSE) #阴影平面
wire3d(ellips, col = "#D95F02", lit = FALSE) # 网格椭球
#### shade3d()与wire3d()同时运行，会形成网格椭圆集合。
图13|添加集合区间3D PCA图，rgl.pca4.pdf。
数据表和代码可以QQ交流群文件夹中下载，或微信公众号后台发送“PCA3D”获取。
参考资料
http://www.sthda.com/english/wiki/a-complete-guide-to-3d-visualization-device-system-in-r-r-software-and-data-visualization#d-scatter-plothttps://dmurdoch.github.io/rgl/dev/reference/material.htmlhttps://dmurdoch.github.io/rgl/dev/reference/polygon3d.htmlhttps://dmurdoch.github.io/rgl/dev/reference/rglwidget.html
推荐阅读
R绘图-RDA排序分析
R统计绘图-VPA(方差分解分析)
R统计绘图-RDA分析、Mantel检验及绘图
R统计-PCA/PCoA/db-RDA/NMDS/CA/CCA/DCA等排序分析教程R统计绘图-PCA分析绘图及结果解读(误差线，多边形，双Y轴图、球形检验、KMO和变量筛选等)
R统计-微生物群落结构差异分析及结果解读
R统计绘图-PCA分析及绘制双坐标轴双序图
R统计绘图-分子生态相关性网络分析
R中进行单因素方差分析并绘图R统计-多变量单因素参数、非参数检验及多重比较
R绘图-相关性分析及绘图
R绘图-相关性系数图
R统计绘图-环境因子相关性热图
R统计绘图-corrplot绘制热图及颜色、字体等细节修改
R统计绘图-corrplot热图绘制细节调整2(更改变量可视化顺序、非相关性热图绘制、添加矩形框等)
R数据可视化之美-节点链接图R统计绘图-rgbif包下载GBIF数据及绘制分布图
R统计-正态性分布检验[Translation]
R统计-数据正态分布转换[Translation]
R统计-方差齐性检验[Translation]
R统计-Mauchly球形检验[Translation]R统计绘图-单、双、三因素重复测量方差分析[Translation]
R统计绘图-混合方差分析[Translation]
R统计绘图-协方差分析[Translation]
R统计绘图-One-Way MANOVA
