其他
R 语言笔记：绘制 3D 的主成分分析图（PCA），常用的几种方法
有时候，3D 的 PCA 图是很有用的，所以 PCA 显示 3 个轴在论文中（包括不少好论文）也偶尔可以看到。
例如：Plant chemical traits define functional and phylogenetic axes of plant biodiversity
https://onlinelibrary.wiley.com/doi/full/10.1111/ele.14262
怎么制作一个 3D 的 PCA 图？
第一步，使用 vegan 包运行一个 PCA 分析，
#################################################################
## PCA 分析 ##
#################################################################
# 包
library(tidyverse)
library(vegan)
library(lattice)
library(rayrender)
# 执行 PCA 分析
pca <- iris[, 1:4] |>
scale() |>
rda()
第二步，提取 PCA 得分数据（即获取 PC1、PC2、PC3 的坐标）
lab.pca1 <- sprintf("Axis 1 (%.2f%%)", pca_smr$cont$importance[2, 1] * 100)
lab.pca2 <- sprintf("Axis 2 (%.2f%%)", pca_smr$cont$importance[2, 2] * 100)
lab.pca3 <- sprintf("Axis 3 (%.2f%%)", pca_smr$cont$importance[2, 3] * 100)
# 坐标提取：样点（花朵种类）
st_coords <- summary(pca, scaling = 2)$sites |>
as.data.frame() |>
mutate(type = iris$Species)
# 坐标提取：变量数据（花朵变量）
sp_coords <- (summary(pca, scaling = 2)$species / 2) |>
as.data.frame() |>
rownames_to_column(var = "var")
# 绘图
ggplot(sp_coords, aes(PC1, PC2, color = var)) +
geom_point(data = st_coords, aes(PC1, PC2), inherit.aes = FALSE) +
geom_segment(aes(x = 0, y = 0, xend = PC1, yend = PC2), arrow = grid::arrow(45, unit(10, "pt"))) +
scale_color_viridis_d() +
labs(title = "Correlation plot (scaling = 2)") +
coord_fixed() +
theme_minimal(base_size = 16) +
theme(panel.grid.minor = element_blank())
ggsave("ggplot2-PCA-2d.png", width = 5, height = 4, dpi = 600, bg = "white")
第三步，第一种方法，使用 scatterplot3d 包绘制一个 3D PCA 图（显示 3 个轴），
#################################################################
## 常规的 3D 绘图 ##
#################################################################
# 绘 3D 图（使用 scatterplot3d 或 lattice 包）
st.colors <- setNames(RColorBrewer::brewer.pal(nlevels(iris$Species), "Dark2"), levels(iris$Species))
png("scatterplot3d-PCA-3d.png", width = 5, height = 5, units = "in", res = 600)
scatterplot3d::scatterplot3d(
st_coords$PC1,
st_coords$PC2,
st_coords$PC3,
pch = 19,
color = st.colors[st_coords$type],
xlab = lab.pca1,
ylab = lab.pca2,
zlab = lab.pca3
)
legend("top", legend = levels(iris$Species), pch = 19, col = st.colors)
dev.off()
png("lattice-PCA-3d.png", width = 3, height = 3, units = "in", res = 600)
lattice::cloud(PC3 ~ PC1 * PC2, group = type, pch = 19, data = st_coords)
dev.off()
#################################################################
## 立体感 3D 绘图 ##
#################################################################
# 物种颜色
st.colors <- setNames(RColorBrewer::brewer.pal(nlevels(iris$Species), "BrBG"), levels(iris$Species))
spheres <- pmap_dfr(
list(st_coords$PC1, st_coords$PC2, st_coords$PC3, st_coords$type),
function(PC1, PC2, PC3, type) {
sphere(
PC1,
PC2,
PC3,
radius = 0.1 / 1.5,
material = diffuse(color = st.colors[type])
)
}
)
# 准备线条
l.axis <- 10
l.width <- 0.1 / 10
path.1 <- list(c(0, 0, 0), c(l.axis, 0, 0), c(-l.axis, 0, 0))
path.2 <- list(c(0, 0, 0), c(0, l.axis, 0), c(0, -l.axis, 0))
path.3 <- list(c(0, 0, 0), c(0, 0, l.axis), c(0, 0, -l.axis))
# 变量（及其标签）颜色
sp.colors <- setNames(RColorBrewer::brewer.pal(nlevels(factor(sp_coords$var)), "Set1"), levels(factor(sp_coords$var)))
paths <- pmap_dfr(
list(sp_coords$PC1, sp_coords$PC2, sp_coords$PC3, sp_coords$var),
function(PC1, PC2, PC3, var) {
path(
list(c(PC1, PC2, PC3), c(0, 0, 0)),
width = l.width * 1.5,
material = light(color = sp.colors[var], intensity = 2)
)
}
)
# 准备文字
h.text <- 0.1
lab.pca1 <- sprintf("Axis 1 (%.2f%%)", pca_smr$cont$importance[2, 1] * 100)
lab.pca2 <- sprintf("Axis 2 (%.2f%%)", pca_smr$cont$importance[2, 2] * 100)
lab.pca3 <- sprintf("Axis 3 (%.2f%%)", pca_smr$cont$importance[2, 3] * 100)
texts <- pmap_dfr(
list(sp_coords$var, sp_coords$PC1, sp_coords$PC2, sp_coords$PC3),
function(lab, PC1, PC2, PC3) {
text3d(
label = str_to_title(str_replace(lab, "[.]", " ")),
PC1 * 1.1,
PC2 * 1.1,
PC3 * 1.1,
text_height = h.text * 1.3,
material = light(color = sp.colors[lab], intensity = 2)
)
}
)
# 构造 3D 物体
f.scene <- function(angle1 = 0, angle2 = 0, angle3 = 0) {
generate_ground(depth = -2, material = diffuse(checkercolor = "grey20")) |>
add_object(
group_objects(
rbind(spheres, paths, texts) |>
# add_object(path(path.1, width = l.width / 3, type = "flat", material = light(intensity = 2))) |>
# add_object(path(path.2, width = l.width / 3, type = "flat", material = light(intensity = 2))) |>
# add_object(path(path.3, width = l.width / 3, type = "flat", material = light(intensity = 2))) |>
add_object(text3d(-1, 0, 0, text_height = h.text, labe = lab.pca1, material = light(intensity = 2))) |>
add_object(text3d(0, 1, 0, text_height = h.text, labe = lab.pca2, material = light(intensity = 2))) |>
add_object(text3d(0, 0, 1, text_height = h.text, labe = lab.pca3, material = light(intensity = 2))),
angle = c(angle1, angle2, angle3)
)
) |>
add_object(sphere(x = -5, y = 7, z = 8, radius = 2, material = light(intensity = 20))) |>
add_object(sphere(x = 5, y = 7, z = 8, radius = 2, material = light(intensity = 10)))
}
render_preview(f.scene(angle1 = 0, angle2 = 0, angle3 = 0), fov = 20, lookfrom = c(1, 1, 6), exponent = Inf)
# 渲染 3D 散点图
f.width <- f.height <- 1000 / 1; f.width
render_scene(
f.scene(angle1 = 0, angle2 = 30, angle3 = 0),
filename = "pca-3D.png",
fov = 17,
lookfrom = c(1, 4, 10),
aperture = 0,
width = f.width,
height = f.height,
samples = 300,
# clamp_value = 10, # 此处设置，增加渲染画质，耗时更长
# min_variance = 5e-6, # 此处设置，极大增加渲染画质，耗时更长（可能长达 30 分钟）
parallel = TRUE
)
效果如下，
记录之。