在多组学研究中，重复测量数据是一个宝藏。

相比于传统的单时间点静态对比，追踪同一个体在多个时间点的动态变化，能够帮我们看到更细的异质性：有的人快速上升，有的人缓慢恢复，有的人看似稳定，但其实只是在绝对水平上分层。

近期正在整理一些开源数据，发现一份代谢组学数据，它包含 5 个时间点的重复测量信息，队列是在接受新辅助化疗的晚期直肠癌患者。面对这份数据，我们产生了一个问题：

能否利用这些连续时间点，为每个重要代谢物绘制出个体的“动态响应轨迹”？进一步通过轨迹聚类，发现具有不同响应特征的亚群？

为了回答这个问题，我们进行了探索性数据分析，建立了一套工作流：

原始代谢组矩阵→ log1p 转换→ ICC + 时间效应 + 变化幅度 + 斜率异质性筛选→ 基线中心化→ latrend 单变量轨迹聚类→ BIC / APPA / OCC / Relative Entropy / 完整病例 ARI 敏感性分析

这篇文章就以这次实战为例，讲一讲多时点重复测量数据如何做轨迹聚类。

为什么不直接把所有代谢物扔进聚类？

原始数据中共有 557 个代谢物。

直觉上，我们似乎可以把每个代谢物都送进轨迹聚类模型里，然后让模型自己找模式。但这样做通常会带来三个问题：

1. 很多代谢物几乎没有明显时间变化；

2. 一些代谢物主要体现为个体间绝对水平高低，而不是随时间变化；

3. 即使有显著时间效应，也可能所有人一起升、一起降，没有足够的轨迹异质性。

所以，在聚类前，我们先做了一轮多维度过滤，筛出真正具有纵向动态潜力的候选分子。

第一步：多维度过滤，筛出真正有“动态时间变化潜力”的分子1）log1p 转换：改善右偏和异方差

代谢组丰度通常右偏明显，并且存在极端值。第一步先对原始丰度做 log 转换。

library(data.table)

# x 为某个代谢物的原始丰度向量

if (all(x >= 0, na.rm = TRUE)) {

y_log <- log1p(x)

} else {

# 如果存在负值，先整体平移再取 log

y_log <- log(x - min(x, na.rm = TRUE) + 1)

}

这里需要注意：log1p 不是“消除”右偏和异方差，而是减弱右偏、降低极端值影响、改善方差稳定性。

2）ICC 筛选：排除由个体间稳定差异主导的分子

有些代谢物看起来有差异，但差异主要来自“这个人一直高、那个人一直低”。这种情况下，轨迹聚类很容易分成“高水平组”和“低水平组”，而不是我们真正关心的“变化模式组”。

因此，我们计算了每个代谢物的组内相关系数（Intraclass Correlation Coefficient, ICC）：

ICC = between-subject variance / total variance

如果 ICC 很高，说明变异主要来自稳定的个体间差异。在这次分析中，我们将筛选阈值定义在 ICC < 0.70 ，优先保留个体内随时间变化占比更高的代谢物。

R 代码示意：

library(lme4)

fit_icc <- lmer(y_log ~ 1 + (1 | subject_id), data = df, REML = TRUE)

vc <- as.data.frame(VarCorr(fit_icc))

var_between <- vc$vcov[vc$grp == "subject_id"]

var_resid <- vc$vcov[vc$grp == "Residual"]

icc <- var_between / (var_between + var_resid)

3）LMM 时间效应：确认群体层面确实随时间变化

接下来，我们使用线性混合模型（LMM）检验每个代谢物是否存在整体时间效应。

比较两个模型：

null model: y_log ~ 1 + (1 | subject_id)

time model: y_log ~ time + (1 | subject_id)

R 代码示意：

library(lme4)

m0 <- lmer(y_log ~ 1 + (1 | subject_id), data = df, REML = FALSE)

m1 <- lmer(y_log ~ factor(timepoint) + (1 | subject_id), data = df, REML = FALSE)

lrt <- anova(m0, m1)

p_time <- lrt$`Pr(>Chisq)`[2]

对所有代谢物的 p 值做 BH-FDR 校正：

screen_table[, time_FDR := p.adjust(p_time, method = "BH")]

筛选阈值：

time_FDR < 0.05

4）timepoint logFC range：不只要显著，还要有变化幅度

统计显著性不等于变化幅度有意义。

因此，我们又计算了每个代谢物在不同时间点的平均 log 丰度范围：

logFC_range = max(mean log abundance across timepoints)

- min(mean log abundance across timepoints)

R 代码示意：

fc_range <- df[, .(mean_y = mean(y_log, na.rm = TRUE)), by = timepoint][

, max(mean_y, na.rm = TRUE) - min(mean_y, na.rm = TRUE)

]

保留 logFC_range 位于全体代谢物前 25% 的分子，即：

logFC_range ≥ 75% 分位数

5）OLS 斜率异质性：寻找“每个人变化速度不同”的分子

最后，我们希望候选分子不仅有整体时间变化，还要有个体间变化趋势差异。

做法是：对每个个体单独拟合一条简单线性回归，提取 slope，然后计算 slope 的标准差。

get_subject_slope <- function(d) {

if (sum(!is.na(d$y_log)) < 2) return(NA_real_)

coef(lm(y_log ~ time_num, data = d))["time_num"]

}

slopes <- df[, .(slope = get_subject_slope(.SD)), by = subject_id]

slope_sd_OLS <- sd(slopes$slope, na.rm = TRUE)

筛选阈值：

slope_sd_OLS ≥ 75% 分位数

过滤成果：从 557 个代谢物到 63 个候选分子

在这份数据中，筛选结果如下：

原始总代谢物：557 个

ICC < 0.70 且 time FDR < 0.05：331 个

严格满足 ICC、time FDR、logFC_range ≥ q75、slope_sd_OLS ≥ q75：63 个

这 63 个候选分子 就是最终进入轨迹聚类分析的代谢物。

第二步：基线中心化，让模型关注“相对变化”

在纵向数据中，一个常见问题是：模型可能按照个体的绝对基线水平分组。

比如某个代谢物，A 组从头到尾都高，B 组从头到尾都低。图上看起来分成了两条线，但这未必说明它们有不同的动态响应模式。

为了减少这种基线水平驱动的分组，在进入 latrend 前做了基线中心化：

y_bc,it = log1p(y_it) - log1p(y_i,Baseline)

其中：

• i 表示第 i 个个体；

• t 表示时间点；

• Baseline 表示该个体自己的基线时间点。

R 代码示意：

library(data.table)

setDT(df)

setorder(df, subject_id, time_num)

# df 包含 subject_id, timepoint, time_num, raw_value

# 先 log1p

if (all(df$raw_value >= 0, na.rm = TRUE)) {

df[, y_log := log1p(raw_value)]

} else {

df[, y_log := log(raw_value - min(raw_value, na.rm = TRUE) + 1)]

}

# 提取每个个体的 baseline log 值

baseline_values <- df[timepoint == "Baseline", .(

baseline_log = y_log[1]

), by = subject_id]

# 合并回长表，并计算 baseline-centered value

df <- merge(df, baseline_values, by = "subject_id", all.x = TRUE)

df[, y_bc := y_log - baseline_log]

基线中心化后的 y_bc 可以理解为：

log-scale change from baseline ≈ 相对于个体自身基线的 log fold-change

它的优点是：

• 减少个体初始水平差异的影响；

• 聚类更关注“随时间怎么变”；

• 结果更容易解释为上升型、下降型、恢复型或波动型。

当然，它也有代价：如果一个代谢物的生物学意义主要体现在绝对水平高低，那么 基线中心化可能会削弱分类信号。

所以这一步的目的不是让所有模型都变好，而是让聚类问题从“谁高谁低”转向“谁怎么变化”。

第三步：基于 latrend 的单变量轨迹聚类

筛选出候选分子并完成基线中心化后，就可以进入轨迹聚类分析了。

这里我们使用 R 包 latrend，并通过其统一接口调用底层的 lcmm 模型。

对每个代谢物，分别拟合：

GBTM: Group-Based Trajectory ModelGMM: Growth Mixture Model

候选类别数：

K = 2, 3, 4

时间函数：

quadratic polynomial: poly(time, 2, raw = TRUE)

核心 R 代码

下面是一段可复用的核心代码框架。

library(latrend)

library(lcmm)

library(data.table)

make_methods <- function(K) {

K <- as.integer(K)

gbtm <- lcMethodLcmmGBTM(

fixed = y_bc ~ poly(time_num, 2, raw = TRUE),

mixture = ~ poly(time_num, 2, raw = TRUE),

id = "subject_id",

time = "time_num",

nClusters = K

)

gmm <- lcMethodLcmmGMM(

fixed = y_bc ~ poly(time_num, 2, raw = TRUE),

mixture = ~ poly(time_num, 2, raw = TRUE),

random = ~ time_num,

id = "subject_id",

time = "time_num",

nClusters = K

)

methods <- list(gbtm, gmm)

names(methods) <- c(paste0("GBTM_K", K), paste0("GMM_K", K))

methods

}

methods <- unlist(lapply(2:4, make_methods), recursive = FALSE)

models <- latrendBatch(

methods,

data = df[!is.na(y_bc)],

cartesian = FALSE,

errorHandling = "pass"

)

模型质量怎么判断？不能只看轨迹图

轨迹图很直观，但不能只凭肉眼判断。

在潜类别轨迹模型里，最容易犯的错误是：看到几条线分开了，就以为“分型成立”。实际上，轨迹聚类至少要同时回答三个问题：

1. 模型整体分得清不清？

2. 每个亚群内部分类把握度够不够？

3. 换一批更严格的完整病例样本后，分型是否还能复现？

因此，我们综合使用了以下指标：

下面逐个解释这些指标。

1）group_sizes 与 min group size：先看有没有“小到不可信”的亚群

group_sizes 指的是最终模型中每个潜在轨迹类别包含的个体数。例如：

group_sizes: A:90; B:75

意思是模型把个体分成 A、B 两类，其中 A 类 90 人，B 类 75 人。

min_group_n 则是所有类别中人数最少的那一类：

min_group_n = min(number of subjects in each trajectory class)

为什么要看这个？

因为潜类别模型有时会切出一个人数极少的小类。这个小类可能看起来很“特别”，但也可能只是模型为了追求更低 BIC 而切出来的噪声亚群。

在本次工作流中，我们对 165 例全样本主模型优选 设定：

min_group_n ≥ 20

也就是说，在主分析中，每个轨迹亚群至少要有 20 个个体，才进入优先考虑范围。

但需要注意：96 例完整病例分析的定位是敏感性分析。它主要用于检验“同样的模型类型和类别数，在完整 5 时间点样本中是否能得到相似分型”，因此我们在 96 例模型中报告实际的 group_sizes、APPA、OCC 和 ARI，但没有再次强制使用 min_group_n ≥ 20 作为模型入选门槛。

R 代码示意：

assignments <- trajectoryAssignments(final_model)

class_sizes <- table(assignments)

group_sizes <- paste(names(class_sizes), as.integer(class_sizes), sep = ":", collapse = ";")

min_group_n <- min(as.integer(class_sizes))

group_sizes

min_group_n

2）APPA：每个亚群内部的“分类把握度”

APPA，全称是 Average Posterior Probability of Assignment，即平均后验分配概率。

它衡量的是：

被分到某一类的个体，平均有多大概率真的属于这一类。

举个例子，如果 A 组的 APPA = 0.92，意思是：模型把这些人分到 A 组时，平均后验概率达到 0.92，分类把握度较高。

理论上，所有类别的 APPA 都应该达到：

APPA ≥ 0.70

如果某一类 APPA 很低，说明这个亚群边界模糊，很多个体其实处在“可分可不分”的状态。

在文章结果中，我们展示的是：

min APPA

也就是所有类别中最差的那个 APPA：

min APPA = min(APPA across trajectory classes)

为什么展示最小值？因为平均值可能掩盖问题。一个模型可能大组 APPA 很高，小组 APPA 很低，平均看起来还不错，但小组其实并不可靠。

3）OCC：分类几率相对于随机分配提升了多少

OCC，全称是 Odds of Correct Classification，即正确分类几率比。

相比 APPA，OCC 的优势是：它考虑了类别样本量不均衡的影响。

直观地说，OCC 衡量的是：

模型把个体分到某一类的正确分类几率，相比随机分配提高了多少倍。

常见经验标准是：

OCC > 2.0：通常可认为分类优于随机分配OCC > 5.0：更严格，说明错分风险更低

本次分析中，我们在模型优选时采用了更严格的标准：

所有类别 OCC > 5

和 APPA 一样，文章中展示的是：

min OCC

也就是所有类别中最差的那个 OCC。这样可以避免某个边缘亚群分类质量很差，却被整体结果掩盖。

4）APPA / OCC / group_sizes 的 R 代码

APPA 和 OCC 可以从 lcmm 模型的后验概率表中手动计算。

library(data.table)

calc_quality <- function(model, model_name) {

# latrend 模型底层存放的是 lcmm 对象

lc <- slot(model, "model")

# pprob 包含：个体 ID、最终分类 class、以及 prob1/prob2/... 后验概率

pp <- as.data.table(lc$pprob)

n_total <- nrow(pp)

class_quality <- rbindlist(lapply(sort(unique(pp$class)), function(k) {

prob_col <- paste0("prob", k)

# 该类别人数

n_k <- sum(pp$class == k)

# 该类别占比

pi_k <- n_k / n_total

# APPA：被分到第 k 类个体，其第 k 类后验概率的平均值

appa <- mean(pp[[prob_col]][pp$class == k], na.rm = TRUE)

# OCC：正确分类 odds 相对于随机分配 odds 的提升倍数

occ <- (appa / (1 - appa)) / (pi_k / (1 - pi_k))

data.table(

model = model_name,

class = k,

n = n_k,

proportion = pi_k,

APPA = appa,

OCC = occ

)

}))

# 文章中展示的是最保守的 min APPA / min OCC

summary_quality <- class_quality[, .(

group_sizes = paste0("class", class, ":", n, collapse = ";"),

min_group_n = min(n),

min_APPA = min(APPA, na.rm = TRUE),

min_OCC = min(OCC, na.rm = TRUE)

)]

list(

by_class = class_quality,

summary = summary_quality

)

}

5）Relative Entropy：模型整体边界清不清

相对熵（Relative Entropy，也常写作 Entropy）是一个全局分类不确定性指标。

它的取值范围通常在：

0 到 1

解释上可以这样理解：

Relative Entropy 越接近 1：类别之间边界越清晰，后验概率越集中Relative Entropy 越接近 0：类别之间重叠越严重，分类越不确定

在轨迹模型评估中，常见经验判断是：

Relative Entropy ≥ 0.80：分类边界较清晰，是比较理想的结果0.60–0.80：中等清晰，需要结合 APPA/OCC/ARI 判断< 0.60：整体分类边界模糊，需要谨慎解释

在 latrend 中，可以直接使用 metric() 提取相对熵：

library(latrend)

relative_entropy <- metric(final_model, "relativeEntropy")

relative_entropy

如果希望自己从后验概率表中复核，也可以根据 posterior probability 手动计算。下面给出一个常见计算方式：

calc_relative_entropy <- function(model) {

lc <- slot(model, "model")

pp <- as.data.table(lc$pprob)

# 提取 prob1, prob2, ..., probK

prob_cols <- grep("^prob[0-9]+$", names(pp), value = TRUE)

P <- as.matrix(pp[, ..prob_cols])

n <- nrow(P)

K <- ncol(P)

# 避免 log(0)

P[P <= 0] <- .Machine$double.eps

# classification entropy

entropy_raw <- -sum(P * log(P))

# relative entropy, normalized to 0-1

relative_entropy <- 1 - entropy_raw / (n * log(K))

relative_entropy

}

relative_entropy_manual <- calc_relative_entropy(final_model)

relative_entropy_manual

这也是为什么在实战结果中，除了 APPA/OCC，我们还展示了每个模型的 relative entropy。

6）ARI：全样本与完整病例分型是否可重复

完整病例敏感性分析比较的是：

全样本模型：165 个个体，允许部分时间点缺失完整病例模型：96 个个体，均具有 5 个时间点

核心问题是：全样本模型和完整病例模型是否给出了相似的个体分型？

这里使用 Adjusted Rand Index（ARI） 比较两套聚类标签。

ARI 的取值通常可以这样理解：

ARI ≥ 0.90：高度稳健与可重复0.40–0.90：存在一定一致性，但需要结合图形和其他指标判断ARI < 0.40：分类存在高度不稳定性或随机性，需要警惕ARI 接近 0：一致性接近随机水平ARI < 0：比随机一致性还差

R 代码示意：

library(mclust)

library(data.table)

# all_assign: 165 例全样本模型的分类结果

# complete_assign: 96 例完整病例模型的分类结果

# 两个表都至少包含 id 和 traj_group 两列

all_assign <- fread("final_assignments_all_165_baseline_centered.csv")

complete_assign <- fread("final_assignments_complete_96_baseline_centered.csv")

# 只比较两套结果共有的 96 个完整病例个体

ari_dt <- merge(

all_assign,

complete_assign,

by = "id",

suffixes = c("_all165", "_complete96")

)

ari <- adjustedRandIndex(

ari_dt$traj_group_all165,

ari_dt$traj_group_complete96

)

ari

这里还有一个技术细节：潜类别模型可能发生 label switching，也就是 A/B/C 标签名称在不同模型里被置换。

ARI 本身不受标签名称置换影响，因此适合用于比较两套聚类结果的一致性；但如果要并排画图展示，就需要根据最大重叠关系，把完整病例模型的 A/B/C 标签重标到全样本模型的语义上。

7）这些指标在工作流中如何协同？

我们的实际判断顺序是这样的。

第一步：先看全局 Relative Entropy。

Relative Entropy 相当于先给模型做一个“全身检查”：

Entropy ≥ 0.80：整体分类边界清晰，可以进入下一步局部检查Entropy < 0.60：整体边界模糊，即使图形好看，也要谨慎解释0.60–0.80：中间地带，需要结合 APPA、OCC、ARI 和轨迹形态综合判断

第二步：再看分亚群的 APPA 和 OCC。

在整体边界还可以的前提下，进一步排查每个亚群：

每个亚群 APPA 是否 ≥ 0.70？每个亚群 OCC 是否至少 > 2.0？如果采用更严格标准，是否能达到 OCC > 5.0？

这一步是为了防止某些“人数很少、边界很模糊”的亚群被整体指标掩盖。

第三步：看 group_sizes。

即使 APPA/OCC 达标，如果某个类别只有几个人，也不适合作为稳定亚型解释。本次设置了较保守的：

min_group_n ≥ 20

第四步：最后看 complete-case ARI。

如果全样本模型看起来很好，但 96 例完整病例模型给出完全不同的划分，那说明结果可能受到缺失时间点或样本组成影响。

因此，ARI 是最后的复现性检查：

ARI ≥ 0.90：非常稳健ARI < 0.40：高度不稳定，需要作为警示案例

一句话总结：

Relative Entropy 看整体边界，APPA/OCC 看每个亚群的分类把握度，group_sizes 看亚群是否太小，ARI 看全样本和完整病例之间是否可重复。四者结合，才比单纯看轨迹图可靠得多。

实战结果展示

从下面的实战结果开始，每个代谢物都并排展示两张图：左侧为全样本 165 例模型，右侧为完整 5 时间点的 96 例敏感性分析模型。这样可以同时观察主分析轨迹形态和完整病例中的稳定性。

第一个展示的是 Phenylalanine

全样本模型，165 例：

final_model: GMM_K2

group_sizes: A:90; B:75

min APPA: 0.850

min OCC: 5.199

relative entropy: 0.550

完整病例模型，96 例，均有 5 个时间点：

complete_model: GMM_K2

group_sizes: A:42; B:54

min APPA: 0.866

min OCC: 5.026

relative entropy: 0.584

全样本 vs 完整病例标签一致性：

complete-case ARI: 0.558

-敏感性校验

全样本模型与完整病例模型的 ARI 为 0.558，提示两套分型之间存在中等一致性，但未达到高度稳健水平。

因此，Phenylalanine 的两类轨迹结构具有一定可重复性，但仍可能受到样本完整性或缺失时间点的影响。

-结果说明

在基线中心化后，Phenylalanine 的全样本轨迹可见两类相对变化模式。两组在时间轴上的均值曲线存在差异，提示该代谢物的动态变化并非完全由个体基线水平驱动。

APPA和OCC均达到预设标准，但Relative entropy < 0.60，提示全局分类边界并不十分清晰。该结果更适合作为探索性发现，而非直接作为稳定亚型结论。

第二个例子是 12alpha-Hydroxy-3-oxochola-4,6-dienoate

全样本模型，165 例：

final_model: GBTM_K2

group_sizes: A:42; B:123

min APPA: 0.951

min OCC: 8.485

relative entropy: 0.842

完整病例模型，96 例，均有 5 个时间点：

complete_model: GBTM_K2

group_sizes: A:20; B:76

min APPA: 0.915

min OCC: 10.351

relative entropy: 0.856

全样本 vs 完整病例标签一致性：

complete-case ARI: 1.000

-敏感性校验

全样本模型与完整病例模型的 ARI 为 1.000，说明两套模型在共有个体上的分型完全一致。该结果支持该代谢物轨迹分型具有较高稳健性，受缺失时间点影响较小。

-结果解释

12alpha-Hydroxy-3-oxochola-4,6-dienoate 的全样本轨迹显示出两类较明确的相对变化模式。完整病例图中，两类轨迹的主要形态与全样本结果保持一致，说明该分型结构在样本限制后仍然可见。

APPA和OCC均达到预设标准，Relative entropy > 0.80，提示类别边界清晰。该代谢物是较适合后续分析的候选对象。

其分型同时具备清晰的视觉分离、较高的分类纯度、较好的全局 entropy 以及完整病例复现性，可优先用于后续与疗效、病理反应或代谢通路变化的关联分析。

第三个例子是 N-Acetylisoniazid

全样本模型，165 例：

final_model: GMM_K3

group_sizes: A:122; B:21; C:22

min APPA: 0.871

min OCC: 11.440

relative entropy: 0.888

完整病例模型，96 例，均有 5 个时间点：

complete_model: GMM_K3

原始标签 group_sizes: A:17; B:64; C:15

按全样本语义重标后 group_sizes: A:64; B:15; C:17

min APPA: 0.891

min OCC: 27.438

relative entropy: 0.906

全样本 vs 完整病例标签一致性：

complete-case ARI: 0.966

-敏感性校验

全样本模型与完整病例模型的 ARI 为 0.966，说明在共有个体中，三类分型具有高度一致性。

需要强调的是，ARI 支持分型结构可复现，但不能抵消完整病例模型中小类样本量不足的问题。

-结果解释

N-Acetylisoniazid 的全样本轨迹呈现三类相对变化模式，而非简单的二分类结构。三类轨迹形态在完整病例中仍可观察到。

APPA和OCC均达到预设标准，Relative entropy > 0.80，提示三类结构具有较高分类确定性和较清晰的全局边界。在后续临床关联分析中，应特别注意小亚群样本量对效应估计和置信区间稳定性的影响。

第四个例子是 5-Hydroxyectoine

全样本模型，165 例：

final_model: GBTM_K2

group_sizes: A:82; B:83

min APPA: 0.880

min OCC: 7.275

relative entropy: 0.651

完整病例模型，96 例，均有 5 个时间点：

complete_model: GBTM_K2

group_sizes: A:50; B:46

min APPA: 0.872

min OCC: 7.396

relative entropy: 0.669

全样本 vs 完整病例标签一致性：

complete-case ARI: 0.958

-敏感性校验

全样本模型与完整病例模型的 ARI 为 0.958，提示两套分型在共有个体中具有高度一致性。结合均衡 group sizes，该代谢物的二分类结果具有较好的稳健性。

-结果解释

5-Hydroxyectoine 的基线中心化轨迹可分为两类，且两类在全样本和完整病例图中均呈现相似的分离趋势。相较于极端不均衡分组，该代谢物的两类轨迹在人群规模上更接近。

APPA和OCC均达到预设标准，Relative entropy在0.60–0.80之间，说明各类别后验分配质量达标，但全局边界清晰度处于中等水平。算是适合后续比较分析的候选分子。

其优势在于两类样本量均衡、分类质量达标、完整病例复现性较好；但 Relative entropy 未超过 0.80，因此仍应避免将其解读为边界极其清晰的分型。

最后一个例子是 Flavone

全样本模型，165 例：

final_model: GMM_K2

group_sizes: A:39; B:126

min APPA: 0.794

min OCC: 2.041

relative entropy: 0.519

完整病例模型，96 例，均有 5 个时间点：

complete_model: GMM_K2

group_sizes: A:43; B:53

min APPA: 0.846

min OCC: 5.493

relative entropy: 0.527

全样本 vs 完整病例标签一致性：

complete-case ARI: 0.309

-敏感性校验

全样本模型与完整病例模型的 ARI 为 0.309，低于 0.40，提示两套分型一致性较差，存在明显不稳定性。

-结果解释

Flavone 的轨迹图在全样本和完整病例中均可见两类曲线，单纯从图形上看，该代谢物似乎可以被划分为两类，但视觉分离并不能直接证明分型可靠。

全样本 min APPA 为 0.794，达到 APPA ≥0.70 的最低要求；但 min OCC 为 2.041，仅接近宽松阈值，未达到本次主分析优先采用的 OCC >5 标准，Relative entropy < 0.60，提示整体分类边界不清晰。

96 例完整病例模型的 min APPA 为 0.846，min OCC 为 5.493，Relative entropy < 0.60，说明完整病例模型内部虽然局部指标改善，全局分类边界仍不清晰。

这套流程还能怎么扩展？

本次分析使用的是单变量轨迹聚类。也就是说，每次只针对一个代谢物建模。

在真实的前瞻性临床研究中，还可以继续往下做：

1）把轨迹亚群作为暴露因素

例如，将某个代谢物的轨迹亚型作为核心暴露变量，进一步分析它与疗效、预后或不良事件之间的关系。

可选模型包括：

Cox proportional hazards modellogistic regressionlinear mixed modelGEE

2）引入临床协变量

如果有年龄、性别、BMI、疾病分期、用药、生活方式等元数据，可以在后续关联模型中进行校正。

3）从单分子走向模块层面

单个代谢物的轨迹有时噪声较大。可以先通过 WGCNA、通路富集或 PCA 得到代谢模块得分，再对模块得分做轨迹聚类。

这种思路可能更接近系统生物学层面的动态异质性。

4）尝试多变量潜在类别模型

如果目标是联合多个代谢物的动态变化，可以尝试 Joint / Multivariate LCMM 或其他多变量纵向聚类方法。

本文转自：谷禾健康

转载本文请联系原作者获取授权，同时请注明本文来自牛耀芳科学网博客。
链接地址：https://blog.sciencenet.cn/blog-2040048-1543431.html

上一篇：[转载]儿童发育迟缓，不只是“吃得少”：营养、肠道菌群与生长的复杂关系


欢迎参加科学网十佳博文评选活动！

主办单位：    支持单位：