在R中处理KML3D结果：获取簇标签、绘图及最优簇数选择

1. 获取模型生成的簇标签 #

首先必须保存kml3d()的运行结果，后续操作都基于这个结果对象：

# 重新运行并保存KML3d分析结果
kml_results <- kml3d(trajectory_knee, nbClusters = 2:5, nbRedrawing = 4, toPlot = "both")

针对你选定的簇数（比如3簇），提取每个观测的簇标签：

# 提取3簇模型的簇标签，将数字替换为你最终选择的最优簇数
cluster_labels <- partition(kml_results, nbClusters = 3)
# 将簇标签合并到原始数据集中
data$cluster <- cluster_labels

如果需要查看所有簇数对应的标签，可使用循环批量提取：

# 提取2-5簇所有模型的标签
all_clusters <- lapply(2:5, function(n) partition(kml_results, nbClusters = n))
names(all_clusters) <- paste0("cluster_", 2:5)

2. 绘制簇的轨迹图 #

使用KML包自带的plot()方法，可直接绘制不同簇数下的轨迹，支持同时展示两个指标或单独展示：

# 同时展示Oxford和Eq5d的簇轨迹图（3簇为例）
plot(kml_results, nbClusters = 3, type = "both")

# 单独绘制Oxford评分的簇轨迹
plot(kml_results, nbClusters = 3, varName = "Oxford score")

# 单独绘制Eq5d评分的簇轨迹
plot(kml_results, nbClusters = 3, varName = "Eq5d score")

若需要更精细的可视化，可提取簇中心数据后用ggplot2自定义绘制：

library(ggplot2)
# 提取3簇模型的簇中心数据
cluster_centers <- getCenters(kml_results, nbClusters = 3)
# 整理为长格式用于ggplot
centers_long <- reshape2::melt(cluster_centers, id.vars = "time", 
                               variable.name = "indicator", value.name = "score")
# 绘制自定义轨迹图
ggplot(centers_long, aes(x = time, y = score, color = factor(cluster), group = cluster)) +
  geom_line(size = 1.2) +
  facet_wrap(~indicator) +
  labs(x = "时间点", y = "评分", color = "簇") +
  theme_bw()

3. 结合BIC指标和轨迹图选择最优簇数 #

提取并可视化BIC值 #

KML3d模型的BIC值可直接从结果对象中提取，BIC通常越大越好，同时关注拐点（当簇数增加后BIC提升变缓时，即为候选最优簇数）：

# 提取2-5簇模型的BIC值
bic_values <- sapply(kml_results@models, function(model) model@bic)
# 整理为数据框查看
bic_df <- data.frame(簇数 = 2:5, BIC值 = bic_values)
print(bic_df)

# 绘制BIC随簇数变化的折线图
plot(bic_df$簇数, bic_df$BIC值, type = "b", pch = 16, 
     xlab = "簇数", ylab = "BIC值", main = "BIC与簇数的关系")

选择最优簇数的判断逻辑 #

1. 优先看BIC值：选择BIC最大或提升幅度明显变缓的簇数；
2. 结合Calinski-Harabasz图：该指标越大说明簇内差异越小、簇间差异越大，优先选指标较高的簇数；
3. 验证轨迹图：确保簇的轨迹具有临床/统计上的解释性，避免选择轨迹区分度极低的簇数。

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