如何在含有序/混合变量的PCoA中实现特征向量可视化？

问题 #

我有一个包含多个物种有序变量的数据集，希望使用主坐标分析（PCoA）进行可视化。当将数据视为连续型（数值型）时，可直接使用vegan::vegdist创建Bray相异度指数，ape::pcoa计算主坐标分解，再通过biplot可视化变量：

library(ape)
library(vegan)
library(FD)

df <- data.frame(a = sample.int(4, 20, replace=TRUE),
b = sample.int(4, 20, replace=TRUE),
c = sample.int(4, 20, replace=TRUE),
d = sample.int(4, 20, replace=TRUE),
e = sample.int(4, 20, replace=TRUE))
rownames(df) <- paste0("species_", letters[1:20])

df.distance <- vegdist(df, "bray")
res <- pcoa(df.distance)
#biplot(res)
biplot(res, df)

但由于变量是有序型，vegdist无法处理，因此我改用FD::gowdis计算适用于混合变量的Gower相异度：

df.ordinal <- df
df.ordinal$a <- factor(df.ordinal$a,levels=1:4,labels = c("low","medium","high","veryhigh"),ordered=T)
df.ordinal$b <- factor(df.ordinal$b,levels=1:4,labels = c("low","medium","high","veryhigh"),ordered=T)
df.ordinal$c <- factor(df.ordinal$c,levels=1:4,labels = c("low","medium","high","veryhigh"),ordered=T)
df.ordinal$d <- factor(df.ordinal$d,levels=1:4,labels = c("low","medium","high","veryhigh"),ordered=T)
df.ordinal$e <- factor(df.ordinal$e,levels=1:4,labels = c("low","medium","high","veryhigh"),ordered=T)

df.distance.gower <- gowdis(df.ordinal, ord="podani")
res <- pcoa(df.distance.gower)
biplot(res)

考虑到有序数据后排序结果有所不同，但我无法将叠加变量可视化为特征向量，执行biplot(res.ordinal, df.ordinal)时出现错误：

> biplot(res.ordinal, df.ordinal)
Error in cov(Y, points.stand) : 
  is.numeric(x) || is.logical(x) is not TRUE

推测这是因为变量现在是有序型数据，而非vegdist示例中的连续型。请问是否存在可在混合数据集的PCoA中可视化特征向量/载荷的方法，或是存在理论上无法实现的原因？

解决方案 #

核心原因 #

ape::biplot.pcoa函数要求传入的变量矩阵必须是数值型，因为它内部会计算变量与PCoA轴之间的协方差，而有序因子无法直接参与协方差计算，这就是报错的直接原因。

可行实现方法 #

1. 将有序因子转为数值型适配biplot

利用有序因子的内在顺序信息，将其转换为整数编码的数值矩阵，再传入biplot函数：

# 提取有序因子的整数编码，转为数值矩阵
df.num <- sapply(df.ordinal, function(x) as.integer(x))
rownames(df.num) <- rownames(df.ordinal)

# 用之前计算好的Gower距离PCoA结果绘制biplot
res <- pcoa(df.distance.gower)
biplot(res, df.num)

这种方法保留了有序变量的顺序关系，同时满足biplot的数值输入要求。

2. 使用vegan包的capscale替代方案

vegan::capscale支持基于Gower距离的排序分析，能自动处理有序变量，直接生成带变量载荷的可视化图：

library(vegan)
# 基于Gower距离执行CAP分析（与PCoA逻辑高度关联）
cap_res <- capscale(df.distance.gower ~ ., data = df.ordinal)
# 绘制包含样本点和变量箭头的排序图
plot(cap_res, display = c("sites", "species"))

3. 手动计算载荷并用ggplot2定制可视化

如果需要更灵活的图控，可以手动计算变量与PCoA轴的相关性，再用ggplot2绘制：

library(ggplot2)
# 提取PCoA前两轴的样本得分
pcoa_scores <- as.data.frame(res$vectors[,1:2])
colnames(pcoa_scores) <- c("PC1", "PC2")

# 将有序变量转为数值，计算与PCoA轴的相关性作为载荷
df.num <- sapply(df.ordinal, as.integer)
loadings <- cor(df.num, pcoa_scores)
loadings <- as.data.frame(loadings)
loadings$var <- rownames(loadings)

# 绘制散点图+变量箭头
ggplot(pcoa_scores, aes(x=PC1, y=PC2)) +
  geom_point(size=2) +
  geom_segment(data=loadings, aes(x=0, xend=PC1*2, y=0, yend=PC2*2), 
               arrow=arrow(length=unit(0.2,"cm")), color="#E63946") +
  geom_text(data=loadings, aes(x=PC1*2.2, y=PC2*2.2, label=var), color="#E63946", size=4) +
  theme_bw()

注：代码中的系数2和2.2用于调整箭头长度和标签位置，可根据图的范围自行修改。

理论说明 #

混合数据PCoA的变量载荷计算，本质是量化变量与排序轴的关联程度。只要能将有序/分类变量转换为保留其内在结构的数值（比如有序因子的整数编码），就可以计算这种关联，不存在理论上的不可行性，核心是要适配工具函数的数据类型要求。

内容的提问来源于stack exchange，提问作者Thomas Moore