1. 项目概述与核心思路做纵向数据分析的人十有八九都会遇到这样一个问题个体随着时间的变化轨迹到底能不能分类传统的混合效应模型会给出一个“平均轨迹加上随机偏移”但如果你关心的是“人群是否存在不同发展模式”普通线性混合模型就不太够用了。这时就需要用到组轨迹模型也就是标题里的GBTMGroup-Based Trajectory Models。这个模型最早是Nagin提出的在犯罪学、心理学、医学等领域用得非常多比如青少年攻击行为的发展轨迹、慢性病患者的用药依从性变化、认知衰退的异质性路径等等。trajeR包就是R语言里专做GBTM分析的工具。第一次接触这个包的人多半是被“trajeR”这个名字吸引进来的——它对标的是SAS里的PROC TRAJ如果你以前用过SAS那套流程转过来会觉得很亲切。但和SAS的闭源商业环境不同trajeR直接在R生态里跑能和ggplot2、dplyr这些常用包无缝配合出图、清洗数据、做后续统计分析都在同一个工作流里完成。这篇博文适合谁正在做纵向数据分析、想识别异质性发展轨迹的研究者尤其是手里已经有面板数据或者随访数据但不知道从哪下手的人。我会从模型原理讲到R语言实操再到结果解读和常见坑尽量把路径讲完整。trajeR包的优势在于它实现了Nagin提出的完整GBTM框架包括Censored Normal模型、Zero-Inflated Poisson模型和Logit模型三种分布假设覆盖了连续型、计数型和二元型三种数据类型。它不仅能拟合模型还提供了模型比较指标BIC、后验概率计算、分类诊断等功能等于把整个分析闭环都做完了。1.1 核心需求解析要理解GBTM首先得跳出“回归一个总体平均轨迹”的思维。GBTM的核心假设是总体由若干潜在子群体组成每个子群体有一条属于自己的平均轨迹个体在子群体内部围绕该平均轨迹随机波动。它不要求预先知道个体属于哪个组而是通过模型估计出每个个体属于各组的后验概率再把个体分配到概率最高的那一组。这和我们熟悉的K-means聚类有本质区别。K-means是先计算个体间的距离再聚类而GBTM是基于每个个体完整的纵向轨迹来建模考虑了时间结构。所以它更适合处理“同一个体多次测量”的数据而不是横截面数据。trajeR包的基本流程可以概括为准备数据、指定模型类型、设定组数和多项式阶数、拟合、比较BIC、检查诊断指标、输出结果。看起来不复杂但每一步都有很多细节坑比如数据格式不对、模型不收敛、组数选多了导致某组样本量太少等等。我在后面的章节会把这些坑一个个讲清楚。2. 环境搭建与数据准备2.1 R环境与trajeR包安装在跑trajeR之前先把R环境准备好。这篇内容假设你已经安装了R语言和RStudio如果还没装先去官网下载R和RStudio版本建议R 4.0以上。trajeR包对R版本有一定要求我自己在R 4.2和R 4.3上都跑过都没问题。trajeR包目前没有发布到CRAN官方仓库需要从GitHub安装这一步是新手最容易卡住的地方# 先安装devtools包 install.packages(devtools) # 从GitHub安装trajeR devtools::install_github(hadjiakbarh/trajeR)安装过程中会自动下载依赖包包括Rcpp、MASS、ggplot2这些。如果在Windows上编译报错多半是缺少Rtools去Rtools官网下载对应版本装上就行。注意Rtools版本必须和你的R版本匹配这一点非常容易踩坑。2.2 数据格式与结构要求trajeR对数据格式有比较严格的要求这一节说清楚能帮你省半天的时间。GBTM分析需要纵向数据每个个体在多个时间点有观测值。trajeR要求的输入格式是宽格式也就是说每一行代表一个个体每一列代表一个时间点的测量值。举例来说如果你有100个个体每个个体在5个时间点有测量那么数据框应该有100行和至少5列测量值。# 模拟一个适合trajeR的宽格式数据集 set.seed(123) n - 300 # 300个个体 time_points - 5 # 每个个体5次测量 # 生成时间变量 time - seq(1, 5, length.out time_points) # 真实分组Group 1: 低水平平稳, Group 2: 中水平上升, Group 3: 高水平下降 group - sample(1:3, n, replace TRUE, prob c(0.4, 0.35, 0.25)) data - data.frame(id 1:n) # 为每个组生成不同的轨迹 for (t in 1:time_points) { col_name - paste0(Y, t) data[[col_name]] - NA for (i in 1:n) { if (group[i] 1) { data[[col_name]][i] - 2 rnorm(1, sd 0.5) } else if (group[i] 2) { data[[col_name]][i] - 2 1.2 * time[t] rnorm(1, sd 0.5) } else { data[[col_name]][i] - 8 - 1.5 * time[t] rnorm(1, sd 0.5) } } } head(data)数据里有时间点但时间变量不需要放在测量值列里。trajeR通过一个独立的参数来指定时间变量向量后面实操部分会讲到。还有一点trajeR不支持时间点不规则的个体也就是说每个个体必须都有完整的时间点测量如果你只有部分时间点的数据需要在拟合之前在数据框里把缺失的时间点设置为NA。2.3 缺失值处理建议GBTM对缺失值的处理并不复杂但有一个原则如果某个个体在某一个时间点缺失它在那一列依然是NA但在拟合时trajeR会删除该个体在那一列的记录。如果缺失比例过高我建议在拟合前就直接删除这些个体否则模型会基于很少的记录来估计那条轨迹很不稳定。更稳妥的做法是在建模前先做一次缺失值数据结构分析# 计算每个个体的缺失情况 data$miss_count - rowSums(is.na(data[, paste0(Y, 1:5)])) table(data$miss_count) # 如果某个体缺失超过2个时间点建议删除 data_clean - data[data$miss_count 2, ]这样处理之后剩下的数据再去跑trajeR模型稳定性和结果可解释性都会好很多。3. 实操核心——拟合组轨迹模型3.1 trajeR函数基本结构trajeR的核心函数就叫trajeR语法格式如下trajeR(Y, A, deg, Model, ng, ...)其中Y宽格式的测量值数据框不含时间变量A时间变量向量deg每组轨迹的多项式阶数向量Model模型类型可选CNORM正态截断、ZIP零膨胀泊松、LOGIT二分类ng组的数量初次接触时最容易犯的错误是deg参数。比如说你打算拟合3组每组用二次多项式抛物线那么deg应该写成c(2,2,2)而不是一个单独的2。trajeR会把这个参数和ng一一对应起来如果长度不匹配直接报错。3.2 选择模型类型选择Model参数是整个分析中最关键的决策之一它直接决定了你对数据的假设。如果你的测量指标是连续性分数比如量表得分、体重、血压用CNORMCensored Normal模型。这个模型假设数据服从正态分布但考虑到可能存在天花板或地板效应允许设置阈值。trajeR的CNORM模型可以设置一个最小值和最大值避免预测值超出合理范围。如果你的测量指标是计数数据比如不良事件次数、犯罪次数、住院次数用ZIP模型。ZIP模型的特别之处在于它把“零”作为一个特殊群体来处理因为计数数据里常出现大量零值——比如多数人没有犯罪记录只有少数人有这时候普通的泊松模型会低估零的比例。如果你关心的是某个二分类结果比如是否患病、是否就业、是否抑郁用Logit模型。它拟合的是事件发生的概率每个组的轨迹表示该组随时间变化的发生概率。从我的使用经验来看CNORM是大多数人的选择因为很多人用的是连续性结局指标。我在这篇博文里以CNORM模型为主来演示但ZIP和Logit的调用方式几乎完全相同只是结果解释时注意概率和期望值的差异。3.3 多项式阶数选择deg参数决定了每条轨迹的形态。通常的做法是先试探性拟合一个简单的模型比如所有组都用一阶线性然后逐步提高阶数用BIC来评估。# 加载trajeR包 library(trajeR) # 准备数据 Y - data_clean[, paste0(Y, 1:5)] # 测量值 A - c(1, 2, 3, 4, 5) # 时间点 # 拟合一个简单的3组模型所有组线性 model_3g_linear - trajeR(Y Y, A A, deg c(1,1,1), Model CNORM, ng 3) summary(model_3g_linear)这里的summary输出包括BIC值、每组样本量占比、参数估计等。BIC越低模型越好但要注意BIC不能无限制地靠增加组数来优化否则会过拟合。3.4 组数选择策略组数ng的选择没有绝对正确的答案需要结合统计指标和专业理论来综合判断。我常用的思路是先设定一个比较宽的搜索范围比如从2组到5组每组用一阶或二阶多项式拟合然后记录每个模型的BIC。results - data.frame(ng 2:5, BIC NA) for (g in 2:5) { m - trajeR(Y Y, A A, deg rep(2, g), Model CNORM, ng g) results[results$ng g, BIC] - m$BIC } print(results)然后比较BIC的变化趋势。BIC通常会随着组数增加而下降但下降幅度会越来越小。一般的经验法则是选“BIC下降幅度明显变缓”的那个点也就是所谓的“拐点”。还有一个辅助指标是平均后验概率每个组的平均后验概率最好都在0.7以上至少不低于0.6。我在实际分析中发现如果某组的人数占比低于5%这组可能没有实际意义哪怕BIC显示它应该存在。这种情况在医学研究中很常见——一个极小的亚组在统计上显著但临床意义存疑。3.5 实例演示3组CNORM模型拟合回到我们的模拟数据假设我们决定拟合3组、每组二阶多项式的模型# 拟合3组模型每组二次多项式 model_3g - trajeR(Y Y, A A, deg c(2,2,2), Model CNORM, ng 3) summary(model_3g)summary输出会包含模型的对数似然值、AIC、BIC、每组的参数估计Intercept、Linear、Quadratic、每组的样本量等。trajeR还提供了confint方法可以查看参数估计的置信区间confint(model_3g)这里有一个值得关注的细节即使你设置了每组二次多项式有些组的二次项系数也可能不显著。这意味着该组的真实轨迹可能是线性的过拟合会引入不必要的参数增加模型复杂度。所以一个常用的策略是先拟合一个全部为高阶的模型然后根据显著性检验把不显著的阶数降下来再重新拟合比较BIC。这种做法在trajeR里可以通过手动调整deg向量来完成是一个很灵活的模型优化方式。3.6 后验概率与分类拟合完成后trajeR会自动计算每个个体属于每个组的后验概率。你可以用postprob方法提取后验概率矩阵posterior - postprob(model_3g) head(posterior)每个个体会被分配到后验概率最高的那一组。你可以把这个分组信息添加到原始数据里方便后续分析data_clean$assigned_group - posterior$group然后你可以用这个分组变量做后续的交叉分析比如各组在基线特征上的差异比较或者和外部结局变量的关联分析。这就是GBTM的灵活之处——它为后续分析提供了一个“潜在分类变量”这是传统方法做不到的。4. 结果可视化与输出4.1 绘制轨迹图模型拟合完最重要的一步就是把轨迹画出来。一张好的轨迹图能直观传达“这几组在怎么发展变化”的信息比任何表格都有说服力。trajeR提供了一个绘图函数trajPlottrajPlot(model_3g, Y Y, A A, Model CNORM, main 3-Group Trajectory Model, xlab Time, ylab Y Value)这个函数会画出三条曲线分别是每组估计的平均轨迹并附带95%置信区间。如果你想把观测数据点和轨迹叠放在一起可以在调用时加上ciTRUE参数或者用ggplot2自己画。4.2 自己用ggplot2画图trajeR自带的绘图功能够用但我更推荐自己用ggplot2画这样图表样式可以完全自定义配色、字体、尺寸都能和论文或报告的要求匹配。画图前需要先提取每组的预测轨迹library(ggplot2) # 提取参数 params - model_3g$params # params是一个矩阵每行对应一组列包括Intercept, Linear, Quadratic等 # 生成预测值 time_seq - seq(1, 5, length.out 100) predictions - data.frame() for (g in 1:3) { pred - params[g, Intercept] params[g, Linear] * time_seq params[g, Quadratic] * time_seq^2 temp - data.frame(Time time_seq, Y pred, Group as.factor(g)) predictions - rbind(predictions, temp) } # 绘制轨迹 ggplot(predictions, aes(x Time, y Y, color Group)) geom_line(size 1.2) labs(title 3-Group Trajectory Model, x Time, y Predicted Y) theme_minimal()这里的params提取方式是一个通用思路实际使用时你需要查看model_3g$params的具体结构来确定列名称。如果结果显示某些组的“Quadratic”为0那说明该组在优化中被降阶了如果你在trajeR调用中设置了不同的deg画图时也会只画线性趋势。4.3 平均后验概率和组间差异除了预测轨迹图还有一个重要的诊断指标叫平均后验概率Average Posterior Probability简称AvePP。每组应该计算该组所有个体的后验概率均值理想情况下应该在0.7以上。如果某组的AvePP很低说明这组区分度不够组间存在大量模糊分类的个体。# 计算每组的平均后验概率 assignments - posterior$posterior # 后验概率矩阵 groups - posterior$group # 分组结果 avepp - sapply(1:3, function(g) mean(assignments[groups g, g])) names(avepp) - paste0(Group , 1:3) print(avepp)当发现某组AvePP低时处理方法有几个方向一是减少组数让模糊的个体被合并到其他组二是增加多项式阶数看是否能通过更灵活的轨迹形状提高区分度三是检查是否有个体存在极端观测值考虑是否删除后重新拟合。4.4 组间分布与后续分析最后你可能会想比较各个潜在组在基线特征上的差异或者评估这些组对远期结局的预测能力。这些分析在获得分组变量后直接进行即可。我这里举一个简单的例子比较各组的某个基线变量# 添加分组标签 data_clean$group_label - factor(data_clean$assigned_group, labels c(Low Stable, Moderate Increase, High Decrease)) # 比较基线变量假设有一个基线变量baseline # data_clean$baseline - rnorm(nrow(data_clean), mean 5, sd 2) # 例如用线性模型看组间差异 lm_model - lm(baseline ~ group_label, data data_clean) summary(lm_model)这样整个GBTM分析流程就闭环了从数据准备到模型拟合再到可视化与后续分析trajeR都提供了配套工具。5. 常见问题与排查技巧实录5.1 安装报错Rtools问题这是最常遇到的问题。Windows用户在执行install_github时如果出现“compilation failed”或者“SystemRequirements: C11”之类的报错几乎都是Rtools没装或者不匹配。解决办法是查看你的R版本到Rtools官网下载对应版本在安装时勾选“Add Rtools to PATH”这个选项然后重启RStudio再装一次。我这里强调一下R 4.2对应Rtools42R 4.3对应Rtools43版本不匹配确实会白折腾一场。还有一个小技巧如果GitHub下载实在太慢或者连接不稳定可以先把仓库打包下载到本地再用install.packages(trajeR_0.1.0.tar.gz, repos NULL, type source)本地安装。不过需要注意本地安装依然需要Rtools编译源码。5.2 模型不收敛或参数异常使用trajeR时偶尔会遇到模型迭代不收敛或者收敛后参数估计出现极端值。这种情况我遇到过几次原因多半是初始值不当或者组数过多导致参数空间太大。trajeR允许你手动设定迭代的最大次数默认值是1000次。如果模型在1000次迭代后还没收敛可以尝试把这个值调大。可以通过options参数来设置model_3g - trajeR(Y Y, A A, deg c(2,2,2), Model CNORM, ng 3, itermax 2000)还有一个办法是先拟合一个较简单的模型比如所有组只取线性deg c(1,1,1)得到参数估计后作为初始值再拟合更复杂的模型。trajeR的文档中没有直接提供设置初始值的接口但通过逐步增加复杂度可以缓解这类问题。还有一种参数异常的情况某组的样本量占比特别小比如只占总体的1%~2%。这时候参数估计往往很不稳定置信区间极宽。我的建议是如果某组占比过低减少组数重新拟合或者考虑把它和其他相似的组合并。5.3 数据格式错误trajeR要求Y参数是宽格式数据框A是时间变量的数值向量。一个常见错误是数据里包含了个体ID列但没有从Y中剔除。比如# 错误示例Y里包含了ID列 trajeR(Y data_clean[, c(id, Y1, Y2, Y3, Y4, Y5)], A c(1, 2, 3, 4, 5), deg c(2,2,2), Model CNORM, ng 3)这样会导致trajeR把ID列也当作一个测量时间点来建模结果会非常离谱。务必确保Y只包含各时间点的测量值。还有一类问题出现在A的长度和Y的列数不一致时。如果你有5个时间点A就必须是长度为5的数值向量。如果时间点之间间隔不均匀比如第1、第3和第10个月有测量A应该写成c(1,3,10)trajeR会在多项式拟合时自动处理这个不等距设计。5.4 组数选择的纠结这是GBTM分析中最让人头疼的问题。有时4组模型BIC确实比3组低但4组中的某一组在专业上很难解释。我的态度是统计指标是参考领域理论是根本依据。一个在专业上难以解释的组即使统计指标更优也不应该强行保留。我有一次在分析慢性病患者用药依从性时5组模型的BIC明显优于4组但第5组只包含了6个人总样本量800多而且这6个人的特征非常离散找不到共同点。最终我选择了4组模型因为在临床场景里这4组的分类更有指导意义。后来在审稿时审稿人也比较认可这种做法认为“组的选择应当具有实质性的可解释性”。5.5 ZIP模型的零膨胀问题如果你用的是ZIP模型要注意trajeR的输出参数里包含零膨胀概率的参数。ZIP模型假设一部分个体的计数是结构性的零值另一部分个体在一个潜在的泊松过程中随机变化。这个零膨胀概率会随时间变化。对于ZIP模型输出的参数有每组轨迹的参数通常是对数尺度的还有一个零膨胀概率的参数。解读时需要注意轨迹曲线画出来的可能是期望计数而不是某个特定的发生率。如果你的数据中零的比例过高比如超过50%ZIP模型可能比重力模型更好但前提是你对这种“结构性零”有理论上的预期。6. 实操心得体会在多个项目里用过trajeR之后我总结了几条亲测有效的经验。第一建模前花时间在数据探索上绝对值得。先看看数据里有多少人完成了全部时间点的测量、每个时间点的分布形态、是否存在明显的天花板或地板效应这些都会直接影响Model和deg的初始设定。第二不要在BIC上做“死磕”。很多刚开始用GBTM的人会陷入一个误区不断加组、加阶数直到BIC降到最低。最终得到的是一个过度拟合、无法解释的结果。模型只是工具文章的核心是领域故事。第三如果条件允许结合多重插补法处理缺失值后再跑trajeR要比直接listwise deletion稳健得多。我通常用mice包做多重插补然后在每个插补数据集上分别跑GBTM再汇总结果。这样能更充分地利用数据信息得到的组分配也更稳定。第四trajeR的输出对象结构不算特别直观建议在拟合完模型后立刻用str()查看一下对象内部结构了解哪些元素对应参数估计、哪些对应BIC、哪些对应后验概率。在自动化分析时这种对对象结构的熟悉会让你省很多时间。最后再分享一个小技巧如果论文或报告里需要展示分组后的个体归属概率不妨画一个“全样本后验概率热图”把每个个体在每组的后验概率用颜色深浅表示出来。这张图能让读者直观看到分类的清晰度。实现方式很简单就是用ggplot2的geom_tile绘制这个矩阵。我做了几次之后发现它比任何文字描述都更有说服力。 R语言trajeR包实现GBTM组轨迹模型：从原理到实战
1. 项目概述与核心思路做纵向数据分析的人十有八九都会遇到这样一个问题个体随着时间的变化轨迹到底能不能分类传统的混合效应模型会给出一个“平均轨迹加上随机偏移”但如果你关心的是“人群是否存在不同发展模式”普通线性混合模型就不太够用了。这时就需要用到组轨迹模型也就是标题里的GBTMGroup-Based Trajectory Models。这个模型最早是Nagin提出的在犯罪学、心理学、医学等领域用得非常多比如青少年攻击行为的发展轨迹、慢性病患者的用药依从性变化、认知衰退的异质性路径等等。trajeR包就是R语言里专做GBTM分析的工具。第一次接触这个包的人多半是被“trajeR”这个名字吸引进来的——它对标的是SAS里的PROC TRAJ如果你以前用过SAS那套流程转过来会觉得很亲切。但和SAS的闭源商业环境不同trajeR直接在R生态里跑能和ggplot2、dplyr这些常用包无缝配合出图、清洗数据、做后续统计分析都在同一个工作流里完成。这篇博文适合谁正在做纵向数据分析、想识别异质性发展轨迹的研究者尤其是手里已经有面板数据或者随访数据但不知道从哪下手的人。我会从模型原理讲到R语言实操再到结果解读和常见坑尽量把路径讲完整。trajeR包的优势在于它实现了Nagin提出的完整GBTM框架包括Censored Normal模型、Zero-Inflated Poisson模型和Logit模型三种分布假设覆盖了连续型、计数型和二元型三种数据类型。它不仅能拟合模型还提供了模型比较指标BIC、后验概率计算、分类诊断等功能等于把整个分析闭环都做完了。1.1 核心需求解析要理解GBTM首先得跳出“回归一个总体平均轨迹”的思维。GBTM的核心假设是总体由若干潜在子群体组成每个子群体有一条属于自己的平均轨迹个体在子群体内部围绕该平均轨迹随机波动。它不要求预先知道个体属于哪个组而是通过模型估计出每个个体属于各组的后验概率再把个体分配到概率最高的那一组。这和我们熟悉的K-means聚类有本质区别。K-means是先计算个体间的距离再聚类而GBTM是基于每个个体完整的纵向轨迹来建模考虑了时间结构。所以它更适合处理“同一个体多次测量”的数据而不是横截面数据。trajeR包的基本流程可以概括为准备数据、指定模型类型、设定组数和多项式阶数、拟合、比较BIC、检查诊断指标、输出结果。看起来不复杂但每一步都有很多细节坑比如数据格式不对、模型不收敛、组数选多了导致某组样本量太少等等。我在后面的章节会把这些坑一个个讲清楚。2. 环境搭建与数据准备2.1 R环境与trajeR包安装在跑trajeR之前先把R环境准备好。这篇内容假设你已经安装了R语言和RStudio如果还没装先去官网下载R和RStudio版本建议R 4.0以上。trajeR包对R版本有一定要求我自己在R 4.2和R 4.3上都跑过都没问题。trajeR包目前没有发布到CRAN官方仓库需要从GitHub安装这一步是新手最容易卡住的地方# 先安装devtools包 install.packages(devtools) # 从GitHub安装trajeR devtools::install_github(hadjiakbarh/trajeR)安装过程中会自动下载依赖包包括Rcpp、MASS、ggplot2这些。如果在Windows上编译报错多半是缺少Rtools去Rtools官网下载对应版本装上就行。注意Rtools版本必须和你的R版本匹配这一点非常容易踩坑。2.2 数据格式与结构要求trajeR对数据格式有比较严格的要求这一节说清楚能帮你省半天的时间。GBTM分析需要纵向数据每个个体在多个时间点有观测值。trajeR要求的输入格式是宽格式也就是说每一行代表一个个体每一列代表一个时间点的测量值。举例来说如果你有100个个体每个个体在5个时间点有测量那么数据框应该有100行和至少5列测量值。# 模拟一个适合trajeR的宽格式数据集 set.seed(123) n - 300 # 300个个体 time_points - 5 # 每个个体5次测量 # 生成时间变量 time - seq(1, 5, length.out time_points) # 真实分组Group 1: 低水平平稳, Group 2: 中水平上升, Group 3: 高水平下降 group - sample(1:3, n, replace TRUE, prob c(0.4, 0.35, 0.25)) data - data.frame(id 1:n) # 为每个组生成不同的轨迹 for (t in 1:time_points) { col_name - paste0(Y, t) data[[col_name]] - NA for (i in 1:n) { if (group[i] 1) { data[[col_name]][i] - 2 rnorm(1, sd 0.5) } else if (group[i] 2) { data[[col_name]][i] - 2 1.2 * time[t] rnorm(1, sd 0.5) } else { data[[col_name]][i] - 8 - 1.5 * time[t] rnorm(1, sd 0.5) } } } head(data)数据里有时间点但时间变量不需要放在测量值列里。trajeR通过一个独立的参数来指定时间变量向量后面实操部分会讲到。还有一点trajeR不支持时间点不规则的个体也就是说每个个体必须都有完整的时间点测量如果你只有部分时间点的数据需要在拟合之前在数据框里把缺失的时间点设置为NA。2.3 缺失值处理建议GBTM对缺失值的处理并不复杂但有一个原则如果某个个体在某一个时间点缺失它在那一列依然是NA但在拟合时trajeR会删除该个体在那一列的记录。如果缺失比例过高我建议在拟合前就直接删除这些个体否则模型会基于很少的记录来估计那条轨迹很不稳定。更稳妥的做法是在建模前先做一次缺失值数据结构分析# 计算每个个体的缺失情况 data$miss_count - rowSums(is.na(data[, paste0(Y, 1:5)])) table(data$miss_count) # 如果某个体缺失超过2个时间点建议删除 data_clean - data[data$miss_count 2, ]这样处理之后剩下的数据再去跑trajeR模型稳定性和结果可解释性都会好很多。3. 实操核心——拟合组轨迹模型3.1 trajeR函数基本结构trajeR的核心函数就叫trajeR语法格式如下trajeR(Y, A, deg, Model, ng, ...)其中Y宽格式的测量值数据框不含时间变量A时间变量向量deg每组轨迹的多项式阶数向量Model模型类型可选CNORM正态截断、ZIP零膨胀泊松、LOGIT二分类ng组的数量初次接触时最容易犯的错误是deg参数。比如说你打算拟合3组每组用二次多项式抛物线那么deg应该写成c(2,2,2)而不是一个单独的2。trajeR会把这个参数和ng一一对应起来如果长度不匹配直接报错。3.2 选择模型类型选择Model参数是整个分析中最关键的决策之一它直接决定了你对数据的假设。如果你的测量指标是连续性分数比如量表得分、体重、血压用CNORMCensored Normal模型。这个模型假设数据服从正态分布但考虑到可能存在天花板或地板效应允许设置阈值。trajeR的CNORM模型可以设置一个最小值和最大值避免预测值超出合理范围。如果你的测量指标是计数数据比如不良事件次数、犯罪次数、住院次数用ZIP模型。ZIP模型的特别之处在于它把“零”作为一个特殊群体来处理因为计数数据里常出现大量零值——比如多数人没有犯罪记录只有少数人有这时候普通的泊松模型会低估零的比例。如果你关心的是某个二分类结果比如是否患病、是否就业、是否抑郁用Logit模型。它拟合的是事件发生的概率每个组的轨迹表示该组随时间变化的发生概率。从我的使用经验来看CNORM是大多数人的选择因为很多人用的是连续性结局指标。我在这篇博文里以CNORM模型为主来演示但ZIP和Logit的调用方式几乎完全相同只是结果解释时注意概率和期望值的差异。3.3 多项式阶数选择deg参数决定了每条轨迹的形态。通常的做法是先试探性拟合一个简单的模型比如所有组都用一阶线性然后逐步提高阶数用BIC来评估。# 加载trajeR包 library(trajeR) # 准备数据 Y - data_clean[, paste0(Y, 1:5)] # 测量值 A - c(1, 2, 3, 4, 5) # 时间点 # 拟合一个简单的3组模型所有组线性 model_3g_linear - trajeR(Y Y, A A, deg c(1,1,1), Model CNORM, ng 3) summary(model_3g_linear)这里的summary输出包括BIC值、每组样本量占比、参数估计等。BIC越低模型越好但要注意BIC不能无限制地靠增加组数来优化否则会过拟合。3.4 组数选择策略组数ng的选择没有绝对正确的答案需要结合统计指标和专业理论来综合判断。我常用的思路是先设定一个比较宽的搜索范围比如从2组到5组每组用一阶或二阶多项式拟合然后记录每个模型的BIC。results - data.frame(ng 2:5, BIC NA) for (g in 2:5) { m - trajeR(Y Y, A A, deg rep(2, g), Model CNORM, ng g) results[results$ng g, BIC] - m$BIC } print(results)然后比较BIC的变化趋势。BIC通常会随着组数增加而下降但下降幅度会越来越小。一般的经验法则是选“BIC下降幅度明显变缓”的那个点也就是所谓的“拐点”。还有一个辅助指标是平均后验概率每个组的平均后验概率最好都在0.7以上至少不低于0.6。我在实际分析中发现如果某组的人数占比低于5%这组可能没有实际意义哪怕BIC显示它应该存在。这种情况在医学研究中很常见——一个极小的亚组在统计上显著但临床意义存疑。3.5 实例演示3组CNORM模型拟合回到我们的模拟数据假设我们决定拟合3组、每组二阶多项式的模型# 拟合3组模型每组二次多项式 model_3g - trajeR(Y Y, A A, deg c(2,2,2), Model CNORM, ng 3) summary(model_3g)summary输出会包含模型的对数似然值、AIC、BIC、每组的参数估计Intercept、Linear、Quadratic、每组的样本量等。trajeR还提供了confint方法可以查看参数估计的置信区间confint(model_3g)这里有一个值得关注的细节即使你设置了每组二次多项式有些组的二次项系数也可能不显著。这意味着该组的真实轨迹可能是线性的过拟合会引入不必要的参数增加模型复杂度。所以一个常用的策略是先拟合一个全部为高阶的模型然后根据显著性检验把不显著的阶数降下来再重新拟合比较BIC。这种做法在trajeR里可以通过手动调整deg向量来完成是一个很灵活的模型优化方式。3.6 后验概率与分类拟合完成后trajeR会自动计算每个个体属于每个组的后验概率。你可以用postprob方法提取后验概率矩阵posterior - postprob(model_3g) head(posterior)每个个体会被分配到后验概率最高的那一组。你可以把这个分组信息添加到原始数据里方便后续分析data_clean$assigned_group - posterior$group然后你可以用这个分组变量做后续的交叉分析比如各组在基线特征上的差异比较或者和外部结局变量的关联分析。这就是GBTM的灵活之处——它为后续分析提供了一个“潜在分类变量”这是传统方法做不到的。4. 结果可视化与输出4.1 绘制轨迹图模型拟合完最重要的一步就是把轨迹画出来。一张好的轨迹图能直观传达“这几组在怎么发展变化”的信息比任何表格都有说服力。trajeR提供了一个绘图函数trajPlottrajPlot(model_3g, Y Y, A A, Model CNORM, main 3-Group Trajectory Model, xlab Time, ylab Y Value)这个函数会画出三条曲线分别是每组估计的平均轨迹并附带95%置信区间。如果你想把观测数据点和轨迹叠放在一起可以在调用时加上ciTRUE参数或者用ggplot2自己画。4.2 自己用ggplot2画图trajeR自带的绘图功能够用但我更推荐自己用ggplot2画这样图表样式可以完全自定义配色、字体、尺寸都能和论文或报告的要求匹配。画图前需要先提取每组的预测轨迹library(ggplot2) # 提取参数 params - model_3g$params # params是一个矩阵每行对应一组列包括Intercept, Linear, Quadratic等 # 生成预测值 time_seq - seq(1, 5, length.out 100) predictions - data.frame() for (g in 1:3) { pred - params[g, Intercept] params[g, Linear] * time_seq params[g, Quadratic] * time_seq^2 temp - data.frame(Time time_seq, Y pred, Group as.factor(g)) predictions - rbind(predictions, temp) } # 绘制轨迹 ggplot(predictions, aes(x Time, y Y, color Group)) geom_line(size 1.2) labs(title 3-Group Trajectory Model, x Time, y Predicted Y) theme_minimal()这里的params提取方式是一个通用思路实际使用时你需要查看model_3g$params的具体结构来确定列名称。如果结果显示某些组的“Quadratic”为0那说明该组在优化中被降阶了如果你在trajeR调用中设置了不同的deg画图时也会只画线性趋势。4.3 平均后验概率和组间差异除了预测轨迹图还有一个重要的诊断指标叫平均后验概率Average Posterior Probability简称AvePP。每组应该计算该组所有个体的后验概率均值理想情况下应该在0.7以上。如果某组的AvePP很低说明这组区分度不够组间存在大量模糊分类的个体。# 计算每组的平均后验概率 assignments - posterior$posterior # 后验概率矩阵 groups - posterior$group # 分组结果 avepp - sapply(1:3, function(g) mean(assignments[groups g, g])) names(avepp) - paste0(Group , 1:3) print(avepp)当发现某组AvePP低时处理方法有几个方向一是减少组数让模糊的个体被合并到其他组二是增加多项式阶数看是否能通过更灵活的轨迹形状提高区分度三是检查是否有个体存在极端观测值考虑是否删除后重新拟合。4.4 组间分布与后续分析最后你可能会想比较各个潜在组在基线特征上的差异或者评估这些组对远期结局的预测能力。这些分析在获得分组变量后直接进行即可。我这里举一个简单的例子比较各组的某个基线变量# 添加分组标签 data_clean$group_label - factor(data_clean$assigned_group, labels c(Low Stable, Moderate Increase, High Decrease)) # 比较基线变量假设有一个基线变量baseline # data_clean$baseline - rnorm(nrow(data_clean), mean 5, sd 2) # 例如用线性模型看组间差异 lm_model - lm(baseline ~ group_label, data data_clean) summary(lm_model)这样整个GBTM分析流程就闭环了从数据准备到模型拟合再到可视化与后续分析trajeR都提供了配套工具。5. 常见问题与排查技巧实录5.1 安装报错Rtools问题这是最常遇到的问题。Windows用户在执行install_github时如果出现“compilation failed”或者“SystemRequirements: C11”之类的报错几乎都是Rtools没装或者不匹配。解决办法是查看你的R版本到Rtools官网下载对应版本在安装时勾选“Add Rtools to PATH”这个选项然后重启RStudio再装一次。我这里强调一下R 4.2对应Rtools42R 4.3对应Rtools43版本不匹配确实会白折腾一场。还有一个小技巧如果GitHub下载实在太慢或者连接不稳定可以先把仓库打包下载到本地再用install.packages(trajeR_0.1.0.tar.gz, repos NULL, type source)本地安装。不过需要注意本地安装依然需要Rtools编译源码。5.2 模型不收敛或参数异常使用trajeR时偶尔会遇到模型迭代不收敛或者收敛后参数估计出现极端值。这种情况我遇到过几次原因多半是初始值不当或者组数过多导致参数空间太大。trajeR允许你手动设定迭代的最大次数默认值是1000次。如果模型在1000次迭代后还没收敛可以尝试把这个值调大。可以通过options参数来设置model_3g - trajeR(Y Y, A A, deg c(2,2,2), Model CNORM, ng 3, itermax 2000)还有一个办法是先拟合一个较简单的模型比如所有组只取线性deg c(1,1,1)得到参数估计后作为初始值再拟合更复杂的模型。trajeR的文档中没有直接提供设置初始值的接口但通过逐步增加复杂度可以缓解这类问题。还有一种参数异常的情况某组的样本量占比特别小比如只占总体的1%~2%。这时候参数估计往往很不稳定置信区间极宽。我的建议是如果某组占比过低减少组数重新拟合或者考虑把它和其他相似的组合并。5.3 数据格式错误trajeR要求Y参数是宽格式数据框A是时间变量的数值向量。一个常见错误是数据里包含了个体ID列但没有从Y中剔除。比如# 错误示例Y里包含了ID列 trajeR(Y data_clean[, c(id, Y1, Y2, Y3, Y4, Y5)], A c(1, 2, 3, 4, 5), deg c(2,2,2), Model CNORM, ng 3)这样会导致trajeR把ID列也当作一个测量时间点来建模结果会非常离谱。务必确保Y只包含各时间点的测量值。还有一类问题出现在A的长度和Y的列数不一致时。如果你有5个时间点A就必须是长度为5的数值向量。如果时间点之间间隔不均匀比如第1、第3和第10个月有测量A应该写成c(1,3,10)trajeR会在多项式拟合时自动处理这个不等距设计。5.4 组数选择的纠结这是GBTM分析中最让人头疼的问题。有时4组模型BIC确实比3组低但4组中的某一组在专业上很难解释。我的态度是统计指标是参考领域理论是根本依据。一个在专业上难以解释的组即使统计指标更优也不应该强行保留。我有一次在分析慢性病患者用药依从性时5组模型的BIC明显优于4组但第5组只包含了6个人总样本量800多而且这6个人的特征非常离散找不到共同点。最终我选择了4组模型因为在临床场景里这4组的分类更有指导意义。后来在审稿时审稿人也比较认可这种做法认为“组的选择应当具有实质性的可解释性”。5.5 ZIP模型的零膨胀问题如果你用的是ZIP模型要注意trajeR的输出参数里包含零膨胀概率的参数。ZIP模型假设一部分个体的计数是结构性的零值另一部分个体在一个潜在的泊松过程中随机变化。这个零膨胀概率会随时间变化。对于ZIP模型输出的参数有每组轨迹的参数通常是对数尺度的还有一个零膨胀概率的参数。解读时需要注意轨迹曲线画出来的可能是期望计数而不是某个特定的发生率。如果你的数据中零的比例过高比如超过50%ZIP模型可能比重力模型更好但前提是你对这种“结构性零”有理论上的预期。6. 实操心得体会在多个项目里用过trajeR之后我总结了几条亲测有效的经验。第一建模前花时间在数据探索上绝对值得。先看看数据里有多少人完成了全部时间点的测量、每个时间点的分布形态、是否存在明显的天花板或地板效应这些都会直接影响Model和deg的初始设定。第二不要在BIC上做“死磕”。很多刚开始用GBTM的人会陷入一个误区不断加组、加阶数直到BIC降到最低。最终得到的是一个过度拟合、无法解释的结果。模型只是工具文章的核心是领域故事。第三如果条件允许结合多重插补法处理缺失值后再跑trajeR要比直接listwise deletion稳健得多。我通常用mice包做多重插补然后在每个插补数据集上分别跑GBTM再汇总结果。这样能更充分地利用数据信息得到的组分配也更稳定。第四trajeR的输出对象结构不算特别直观建议在拟合完模型后立刻用str()查看一下对象内部结构了解哪些元素对应参数估计、哪些对应BIC、哪些对应后验概率。在自动化分析时这种对对象结构的熟悉会让你省很多时间。最后再分享一个小技巧如果论文或报告里需要展示分组后的个体归属概率不妨画一个“全样本后验概率热图”把每个个体在每组的后验概率用颜色深浅表示出来。这张图能让读者直观看到分类的清晰度。实现方式很简单就是用ggplot2的geom_tile绘制这个矩阵。我做了几次之后发现它比任何文字描述都更有说服力。