R语言LCMM实战：潜在类别混合模型识别纵向数据异质性轨迹

如果你处理过纵向随访数据，大概率遇到过这个尴尬场景：把所有人的测量值汇成一条均值曲线，看上去平缓、规律、方向明确，可一旦按某个变量拆开，数据完全是另一副样子——有人快速下降，有人长期平稳，还有人缓慢爬升。更麻烦的是，基于人群平均轨迹做出的推断，可能对任何个体都不成立。这是纵向数据分析里最经典也最容易忽略的问题：群体异质性。而LCMM（潜在类别混合模型）正是针对这一问题的成熟解法。

这篇文章我打算从增长混合模型与潜在类别轨迹分析的关系切入，完整走一遍用R的lcmm包做轨迹建模的流程，包括模型结构、代码实现、类别数选择、结果解读和我实际使用中踩过的坑。适合刚接触纵向数据建模、想用轨迹分析发现潜在亚组的研究者，也适合已经在用混合模型但对LCMM细节还有疑问的同行。

1. 平均轨迹的陷阱与LCMM的定位

1.1 一条平均曲线掩盖了什么

假设你跟踪了400名受试者6年的认知得分，传统做法是拟合一个混合效应模型，得到一个总体的平均变化斜率，比如每年下降1.2分。这个数字听着很清晰，但它背后可能藏着三种完全不同的真实轨迹：一类人稳定在较高水平，一类人轻度下降，还有一类人呈现断崖式衰退。把这三组混在一起拟合，得到的“平均斜率”其实是三条截然不同轨迹的扭曲折中。

我在实际项目中见过一个典型场景：某药物临床试验的纵向数据，整体分析显示治疗组和对照组的变化趋势没有显著差异。但按LCMM把人群拆成几个潜在轨迹后，发现大约20%的人是快速进展型，治疗在这部分人群中的效果其实非常明显。这种信息如果只看均值模型，根本不会浮现。

所以轨迹分析本质上是回答一个更贴近临床的问题：这组人里面，是否存在不同变化模式的亚组？如果有，这些亚组各有什么特征，对干预或结局的响应有何不同。

1.2 GMM、LCGA、LCMM：三个术语的真实关系

很多人看到三个词就头大：增长混合模型（GMM）、潜在类别轨迹分析（LCGA或GBTM）、潜在类别混合模型（LCMM）。它们到底有什么区别？

简单说，三者的核心都是“用若干个类别的轨迹来描述异质性人群”，区别主要体现在对类别内个体差异的假设上。

• GMM（Growth Mixture Modeling）允许每个类别内部的个体有随机效应，也就是同一类别里的人并不完全同质，他们围绕类别的平均轨迹有个体波动。
• LCGA/GBTM（Latent Class Growth Analysis / Group-Based Trajectory Modeling）假设每个类别内的个体完全同质，不考虑类别内的随机效应，只估计一个平均轨迹。
• LCMM是一个更通用的框架，既可以通过设置random参数包含类别内随机效应，也可以不收随机效应从而退化为LCGA。换句话说，LCMM可以看作GMM和LCGA的统一实现框架。

在R的lcmm包里，hlme()函数配合不同的random设定，可以覆盖这几种情况。默认不指定random时，模型等价于LCGA；指定random = ~1或random = ~time后，就进入GMM的设定范围。

1.3 什么人需要学LCMM

如果你正在处理这样的问题：重复测量的纵向数据、主观上怀疑人群不是单一整体、需要识别不同的变化轨迹并进一步分析轨迹的预测因素，那LCMM就值得学。它在流行病学、临床医学、心理学、社会学、市场营销等领域都有应用，比如认知老化轨迹、抑郁症状演变、收入变化路径、顾客消费行为分型等。

不过也要提醒一点，LCMM面向的是连续型结局变量。如果你的结局是分类变量、计数变量或生存时间，lcmm包也提供了对应扩展（如Jointlcmm、multlcmm），但本文主要讨论连续结局的轨迹建模。

2. 模型结构拆解：混合在哪里，类别从哪来

2.1 从线性混合模型到“类别维度的混合”

理解LCMM最好的方式是从传统线性混合模型出发。线性混合模型通常写作：

y_it = β0 + β1 * time_it + u_i0 + ε_it

其中u_i0是个体随机截距，用于刻画个体偏离总体平均水平的程度。这个模型的所有人都共享同一组固定效应β0和β1，所以它假设人群是同质的。

LCMM则在此基础上引入一个潜变量c_i，表示个体i归属于第k个潜在类别（k = 1,…, K）。每个类别都有自己的一套固定效应参数，例如截距β0k和斜率β1k。模型变成：

y_it = β0k + β1k * time_it + u_i0 + ε_it, 如果c_i = k

这里的“混合”就体现在：总体分布不是单一的正态分布，而是K个类别的混合分布。每个个体以一定概率属于某个类别，这个概率本身也可以通过协变量建模。

2.2 三类效应的分工：fixed、mixture、random

用lcmm包拟合模型时，有三个参数最容易混淆：fixed、mixture、random。

• fixed定义的是所有类别共享的固定效应。比如y ~ time，表示时间和结局的关系是线性的，且这个关系在所有类别中都存在。
• mixture定义的是“随类别变化”的效应。比如mixture = ~time，表示不同类别的斜率可以不一样。如果某个变量在不同类别的效应不同，就应该放进mixture。
• random定义的是类别内部的随机效应。random = ~1表示每个类别的个体允许有随机截距；random = ~time表示允许个体有随机截距和随机斜率。

这三个参数的组合决定了模型的灵活度。刚入门时建议从简单开始，也就是只用fixed和mixture，随机效应先不加，等模型稳定了再考虑逐步增加复杂度。我见过不少人一开始就上全套随机效应，结果模型不收敛，反而打击信心。

2.3 类别归属概率：软分类与后验概率

LCMM不是先做聚类再拟合曲线，而是通过极大似然估计联合估计类别概率和各类别的轨迹参数。每个个体属于每一类都有一个后验概率（posterior probability），所有类别的概率之和为1。分析时通常把个体归入后验概率最大的那个类别，这就是所谓“软分类”转“硬分类”的过程。

后验概率的质量是判断模型好坏的重要依据。理想情况下，一个人应该明显属于某一类，比如概率0.92，而不是在两类之间犹豫摇摆。如果很多个体的最大后验概率只有0.5左右，说明类别界限模糊，模型可能设多了类别，或者变量选择不当。

3. R语言实操：半小时跑通第一个LCMM

3.1 为什么选lcmm包

R中可以做轨迹建模的包不少，常用的包括lcmm、traj、flexmix，以及调用Mplus或latent class相关的接口。我长期用lcmm，原因是它的函数体系完整，能从基础LCMM扩展到非线性、多元轨迹和联合模型，不需要频繁切换工具。traj包更轻量，适合快速做LCGA，但功能相对有限；flexmix是通用的有限混合模型包，功能很强但需要自己拼装纵向设计。

lcmm包由法国研究者Cécile Proust-Lima团队开发，核心函数包括hlme()（线性潜类别混合模型）、lcmm()（更通用的潜类别混合模型，支持非线性联系函数）、multlcmm()（多元轨迹）、Jointlcmm()（联合建模）。安装方式：

install.packages("lcmm")

实测在R 4.x版本下安装无问题，依赖包会自动处理。

3.2 数据准备：长格式是底线

无论什么纵向模型，数据格式都必须是长格式：每一行是一个个体在某个时间点的观测。关键变量包括个体ID、时间变量、结局变量，以及可选的协变量。用一列有代表性的模拟数据举例：

library(lcmm) set.seed(2024) n <- 400 waves <- 0:5 time <- rep(waves, n) id <- rep(seq_len(n), each = length(waves)) # 模拟两个协变量，后续用于类别归属预测 x1 <- rnorm(n) x2 <- rbinom(n, 1, 0.4) # 让类别归属依赖 x1、x2（多项式 logit 设定） eta1 <- 0.3 * x1 - 0.5 * x2 eta2 <- -0.4 * x1 + 0.6 * x2 p1 <- exp(0) / (1 + exp(eta1) + exp(eta2)) p2 <- exp(eta1) / (1 + exp(eta1) + exp(eta2)) p3 <- exp(eta2) / (1 + exp(eta1) + exp(eta2)) true_class <- apply(cbind(p1, p2, p3), 1, function(p) { sample(1:3, size = 1, prob = p) }) # 三条轨迹的真实参数 intercept <- c(48, 50, 46) slope <- c(-0.8, -2.5, 0.3) sigma <- 1.8 ri <- rnorm(n, 0, 1.2) # 类别内的个体随机截距 y <- numeric(length(id)) for (i in seq_len(n)) { idx <- id == i c <- true_class[i] y[idx] <- intercept[c] + ri[i] + slope[c] * time[idx] + rnorm(length(waves), 0, sigma) } dat <- data.frame( id = id, time = time, x1 = rep(x1, each = length(waves)), x2 = rep(x2, each = length(waves)), y = y ) head(dat)

这里模拟的三类轨迹分别是：轻度下降、快速下降、稳定维持。注意类别内加入了随机截距ri，也就是GMM设定。如果你只做LCGA，可以把这个随机截距设为0。

3.3 从单类模型出发：基准模型的意义

拟合多类别模型之前，务必先跑一个单类别模型。它有两个作用：一是作为后续多类别模型的起始值基础，二是提供一个“不划分亚组”的基准，用于比较模型改进程度。

m1 <- hlme(y ~ time, subject = "id", data = dat, ng = 1) summary(m1)

ng = 1表示只有一个类别，此时模型等价于普通的线性混合模型。summary输出的固定效应部分会给出截距和time的估计，随机效应部分给出残差方差和随机截距方差。

3.4 gridsearch与多类别模型拟合

从两类模型开始，逐步增加类别数。多类别模型的极大似然函数是非凸的，存在多个局部最优解，直接用默认起始值跑很容易陷进去。我的做法是固定使用gridsearch函数，重复搜索多个起始值。

m2 <- gridsearch( rep = 30, maxiter = 5, minfit = m1, hlme(y ~ time, subject = "id", data = dat, ng = 2, mixture = ~time) ) m3 <- gridsearch( rep = 30, maxiter = 5, minfit = m1, hlme(y ~ time, subject = "id", data = dat, ng = 3, mixture = ~time) ) m4 <- gridsearch( rep = 30, maxiter = 5, minfit = m1, hlme(y ~ time, subject = "id", data = dat, ng = 4, mixture = ~time) )

这里简单解释一下gridsearch的逻辑：rep = 30表示在基准模型估计值附近随机生成30组起始值；maxiter = 5表示每组起始值最多迭代5次，用来快速淘汰差的起始点；minfit = m1指定了搜索的基准模型。最终它会返回多次搜索中似然最优的结果。

为什么不直接多跑几次相同模型就行？因为网格搜索策略能将起始值控制在一定范围内，相比完全随机的多起点法更高效。在实际项目中，如果数据量大，rep可以降到15~20以节省时间。

4. 一个完整的轨迹发现案例：从数据生成到报告解读

4.1 用拟合指标确定类别数

模型跑完，第一步是比较BIC。BIC是潜类别模型中最受认可的模型选择指标之一，它同时惩罚了似然和参数数量，比AIC更偏好简洁模型。

summarytable(m1, m2, m3, m4, which = c("G", "npm", "AIC", "BIC", "entropy"))

模拟数据下可能得到这样一组结果（数值以你本机运行为准）：

BIC的最小值落在三类模型上，四类模型的BIC反而略有回升，说明增加第四类的获益不足以抵消参数复杂度。AIC在四类时仍然最低，但实际分析中不要单看AIC，因为AIC倾向于选更复杂的模型。以BIC为主、结合解释性来判断，是最稳妥的策略。

还有一个需要注意的点：entropy（熵）用来衡量分类的清晰度，取值在0到1之间，越接近1越好。上面表格里三类的entropy达到0.85，是三个模型中最高的，说明三类划分的边界最干净。

4.2 用后验概率评估分类质量

仅仅BIC最低还不够，还需要检查每个类别的平均后验概率和样本量占比。

postprob(m3)

输出中有一张表，展示每个类别中个体的平均后验概率。经验法则是各类别的平均后验概率最好在0.7以上，如果低于0.7，说明这个类别很可能不稳定，需要通过减少类别数或调整模型结构来改善。模拟数据的输出一般会比较理想，比如三类平均后验概率分别接近0.90、0.92、0.88。

同时要检查每个类别的样本量占比。类别占比没有绝对标准，但一个类别如果占比小于5%，它的轨迹估计会非常不稳定，后续关联分析也很难有统计效能。此时要考虑减少类别数，或者看是否该类别由极端值驱动。

4.3 轨迹图：预测均值如何解读

确定三类模型后，画轨迹图是必不可少的环节。

plot(m3, which = "fit", legend.loc = "topright", lwd = 2)

这张图会展示三条估计的平均轨迹。如果你是用模拟数据，可以对照真实类别计算个体轨迹的均值，你会发现模型恢复得相当好。

读图时有个容易忽视的细节：图上的轨迹是模型预测的均值轨迹，不是每个类别内的个体观测均值。当类别内随机效应较大时，个体观测点会明显散落在预测轨迹两侧。汇报结果时，应该同时报告预测轨迹的参数估计和类别内个体变异的大小，而不是只丢一张图。

4.4 把类别当作新变量：后续关联分析怎么做

轨迹模型跑完后，最自然的后续问题就是：哪些因素决定了这个人属于快速下降组而不是稳定组？这时可以使用classmb参数，在建模的同时直接估计类别归属的预测因素。

m3_cov <- gridsearch( rep = 30, maxiter = 5, minfit = m3, hlme(y ~ time, subject = "id", data = dat, ng = 3, mixture = ~time, classmb = ~ x1 + x2) ) summary(m3_cov)

classmb指定的是一个多项式logit模型，类似多分类逻辑回归，用于预测个体归属每个类别的概率。默认以第一个类别为参照，输出的系数表示相对于参照类别，该变量增加一个单位时，个体更可能归属哪个类别。

也可以先保存每个个体的最大后验概率类别，再用传统的多分类逻辑回归做分析。两种思路结果接近，但classmb一次性联合估计，统计效率更高，标准误也更可信。

5. 拟合LCMM时最常见的五个“翻车场景”

5.1 收敛失败与“怪异的类别”

模型不收敛，或者收敛后某个类别的轨迹异常（比如斜率是其他类别的十倍），在LCMM里非常常见。我遇到过最典型的情况是：指定mixture = ~time后，某个类别在时间上完全交错，画出来像是噪声。

处理顺序应该是：先降低模型复杂度，比如把random从~time简化到~1甚至去掉；再检查时间变量是否标准化。随访时间跨度很大时（比如0到20年），斜率数量级差异大会让优化过程很吃力，建议把时间缩放到0到1或均值中心化。最后再考虑增加gridsearch的重复次数。

5.2 局部最优解：跑两次结果不一样

我早期用lcmm时犯过一个错：同样的数据、同样的类别数，换了一组起始值，BIC差了几十，轨迹形态完全变了。这就是局部最优解的典型表现。函数默认使用的起始值并不是全局最优的保证。尤其类别数超过3个以后，似然面的复杂程度会急剧上升。

解决方案就是坚持用gridsearch，并且在重要分析中至少跑两次不同随机种子的搜索，确认最优结果一致。如果两次结果不一样，以BIC最小的为准，同时报告这一不稳定性。审稿人如果问起来，这反而是分析严谨的体现。

5.3 类别太小或太大：是该合并还是该拆

模拟数据里三类占比大约40%、35%、25%，分布均匀。但真实数据里经常出现占3%的极小类别，而且这个类别的轨迹往往非常极端——比如全部快速下降且终点很低。这时不要急着高兴。

小类别可能代表真实但罕见的亚组，也可能只是少数极端值“拉”出来的伪类别。常规建议是：先看该类别个体的原始数据，确认不是数据录入错误；然后把类别数减一，比较BIC的变化；如果BIC依然支持小类别模型，且该类别符合领域理论预期，才考虑保留。这个判断过程务必写进报告，因为读者也会关心。

5.4 数据稀疏时怎么办

随访次数太少会直接影响轨迹估计。每类至少要有足够的重复测量次数，我一般建议每个体至少3次随访，类别数不超过4个；如果只有2次随访，基本不要考虑随机斜率，模型极容易退化。缺失数据比例高时，lcmm基于全信息极大似然可以处理缺失结局，前提是缺失机制满足随机缺失假设。如果缺失与轨迹本身有关，模型估计就会有偏。

5.5 常见解读误区：类别不是“真实实体”

这一点必须反复提醒：LCMM估计出的类别，是对异质性的一种近似描述，不等同于生物学或社会学意义上真实存在的亚型。两个类别之间可能只是量的差异（比如程度不同），而非质的差异（机制不同）。在论文里谨慎的表述是“识别出三个潜在轨迹组”而不是“发现三种疾病亚型”。很多研究因为措辞不当，让读者对类别产生了过度实体化的理解，这是潜在类别分析被批评最多的点。

6. 从LCMM继续延伸：复杂轨迹建模的可选路线

6.1 非线性轨迹与更灵活的模型

线性假设在很多实际场景下并不成立，比如认知功能往往在某个时间点后加速下降。这时你可以直接加入多项式项，或者用lcmm()函数配合样条。lcmm()是hlme()的泛化版本，支持更灵活的联系函数和非线性轨迹设定：

m_spline <- gridsearch( rep = 30, maxiter = 5, minfit = m1, lcmm(y ~ time, subject = "id", data = dat, ng = 3, mixture = ~time, link = "3-equi-splines") )

使用样条会增加参数数量，对数据量有一定要求。如果样本量不足，先试二次项通常更实际。

6.2 多元轨迹与联合模型

当你同时关注多个纵向指标时，比如认知得分和抑郁得分，可以分别建模，再分析各类别的交叉关系，也可以直接用multlcmm()做多元轨迹联合建模。它允许不同指标拥有不同的类别结构，分析它们之间的关联。

更进一步的玩法是Jointlcmm()，把纵向轨迹和生存事件联合建模。一个典型问题是：认知快速下降轨迹是否与痴呆发病风险相关。联合模型通过共享潜在类别结构，同时估计轨迹和风险，避免了两步分析带来的偏差。

6.3 一个重要的提醒：类别数不是越多越好

最后一件事，也是我每次和同行交流都会强调的：确定类别数，最终依赖的是统计指标、领域解释性、外部验证三者共同判断，而不是机械地选BIC最小的那个。有些场景下，3类模型和4类模型的BIC只差很小的数值，但4类中有一个类别的轨迹形态缺乏实际意义，这时我宁愿选3类。统计模型是为研究问题服务的，不是让研究问题反过来迁就模型的输出。

在我自己的项目里，这套流程已经跑过不少轮。最开始接触LCMM时，我也试图在一个模型里塞进尽可能多的随机效应和协变量，结果换来一堆收敛警告。后来养成了从简单到复杂、从单类到多类的习惯，反而每一步的检查点都清晰了。如果你正准备用LCMM分析自己的纵向数据，我建议你也按这个顺序走：先明确研究问题是否需要识别亚组，再用一两个协变量跑通全流程，最后才逐步增加模型复杂性。轨迹分析真正的价值，不在于把人群拆成几组，而在于拆出的这几组能否带来对问题更深入的理解。