纵向人群研究常对同一对象进行重复测量(如血压、血糖、认知评分及功能状态等),时间跨度可达数十年。不同临床问题需采用相应的纵向数据分析方法。例如,联合模型通过共享随机效应整合纵向过程与生存数据,适用于探索标志物长期变化对结局的影响机制
[1];而超级界标模型兼具灵活性与计算效率,更适合多时间点动态风险分层,尤其在处理多个纵向指标或竞争风险数据时优势显著
[2]。
当研究问题聚焦于人群纵向变化的总体模式时,线性混合效应模型(linear mixed model,LMM)足以刻画总体平均轨迹及个体围绕该轨迹的随机波动。然而,若人群存在异质性,例如:部分个体具有不同的起始水平、变化方向或变化速度,单一平均轨迹难以直接揭示这些潜在模式。群组轨迹模型(group-based trajectory model,GBTM)正是解决这类问题的方法学工具,它以个体完整纵向观测序列为分析单位,通过有限亚组的特异性轨迹近似描述总体中的纵向异质性
[3],可将复杂的个体轨迹转化为可描述、可比较的群体模式,为识别早期影响因素及预测后续结局提供分析依据。该方法源于犯罪学与心理学,现已广泛应用于临床医学和公共卫生领域
[4-6]。
经典GBTM假定总体由若干不可直接观测的轨迹组构成,同一轨迹组内个体共享同组特异性平均轨迹
[3,7]。本文旨在比较GBTM与相关纵向轨迹模型的核心差异,阐述经典GBTM的基本原理、模型选择与报告规范,并结合实例演示其标准分析流程。
基本思想与方法学框架
GBTM属于有限混合模型(finite mixture model)在纵向资料中的一种专门应用
[3]。经典GBTM假定研究人群并非同质,总体由G个不可直接观测的轨迹组构成。每一亚组内的个体共享同一条平均发展轨迹,个体间的残余差异由组内随机误差刻画。设个体
i在时间
接受第
次测量,结局为
,全部纵向观测的完整向量为
;潜在轨迹组数为G,π
g为随机抽取一名个体属于第g组的概率,且各组概率之和为1;
为给定第
组下该个体观测序列的条件概率。个体数据的边际分布可写成各组条件分布的加权和,简洁表示为:
在经典GBTM中,假定给定组别后各次测量条件独立,等于各时点观测概率的乘积,故 ,其中由结局分布决定如正态分布、Bernoulli分布或Poisson分布等。
GBTM不是先按观察曲线聚类再拟合,而是在同一个概率模型中同时估计组比例、组特异性轨迹参数和结局分布参数。经典GBTM组内通常不设置个体随机截距或随机斜率,模型主要用“离散的潜在轨迹组”来描述个体间异质性
[8];同组成员共享一条均值轨迹,组内个体围绕组平均轨迹的偏离由结局分布和误差项吸收。组内同质和条件独立假设使模型结构简明,曲线易于解释。当真实资料存在较强的组内连续差异或个体内相关结构时,GBTM可能用额外类别去吸收本应由随机效应解释的变异,类别数量或形状可能受到影响
[9-10]。
个体的完整纵向观测由所在组特异性轨迹共同解释,每个轨迹组都有自己的平均变化函数,第g组的期望轨迹可写为:
其中,p为该组轨迹的多项式阶数,各组可以不同。
组特异性均值轨迹常用时间的低阶多项式表示
[11]。因而,时间尺度是GBTM建模的第一构件选择。当以实际年龄为时间轴,轨迹描述指标随年龄增长的变化;当以距诊断、距手术或距治疗开始的时间为轴,轨迹描述共同临床起点后的变化。为便于数值计算和参数解释,可对时间变量进行中心化和缩放。例如,以实际年龄为时间轴时,将年龄转换为(age-65)/10,即以65岁为时间原点、10岁为一个时间单位。该变换可减小年龄及其高次项的数值量级,使截距对应65岁时的模型水平。需要注意,以实际年龄拼接不同入组年龄对象的数据时,曲线同时利用个体内随访变化和不同年龄入组对象之间的差异;若存在出生队列效应、时期效应或选择性存活,整条年龄曲线不能简单理解为同一批人从低龄到高龄的自然史
[12]。
另一个GBTM模型构件即结局分布的选择。结局分布应依据数据特征选择,并检验模型拟合。连续、二分类和计数结局可分别采用正态、Bernoulli和Poisson分布,后二者常分别使用logit和对数连接;比例结局可考虑Beta分布,零值较多时可考虑零膨胀模型。具有明确上下界且边界值较多的量表可采用删失正态分布(censored normal,CNORM)
[3,13],如取值为0~30分的简易精神状态检查量表(mini-mental state examination,MMSE)。CNORM将观测值视为潜在连续变量在量表边界处删失后的结果,此处删失不同于生存分析中的失访或右删失。其潜变量均值不等于观测量表均值,报告期望轨迹时须转换至原始量表尺度。
总体似然为各个体边际似然的乘积。模型采用最大似然估计,可通过期望最大化(expectation-maximization,EM)等算法求解。根据估计参数计算每名个体属于各组的后验概率,再按最大后验概率分配轨迹组。
与其他纵向轨迹模型的比较
不同纵向轨迹模型的核心差异(
表1),在于对研究对象之间的异质性如何刻画:异质性被看作总体平均轨迹周围的连续个体差异,还是若干潜在类别之间的离散差异,抑或二者并存。线性混合效应模型用连续随机效应刻画个体围绕同一总体均值的偏离,不产生离散亚组
[14]。GBTM与不设组内随机效应的潜在类别增长分析(latent class growth analysis,LCGA)高度相近
[8]。增长混合模型(growth mixture model,GMM)和潜在类别混合模型(latent class mixed model,LCMM)则在类别间差异的基础上进一步保留随机效应
[9],因而能够同时描述离散与连续异质性,但参数更多、可识别性和收敛更具挑战。因此,模型选择应基于研究问题、样本量、测量次数、可识别性和计算稳定性
[15]。
模型构建、选择与报告
规范的GBTM分析应从明确的研究问题开始,而不是由软件自动寻找类别。研究方案首先应明确结局指标、时间原点、时间尺度和纳入标准;其次检查每名对象的测量次数、观测时间分布、结局上下界、极端值和缺失;最后依据结局类型选择概率分布,并设定合理的时间函数。多项式阶数通常从一次或二次开始,高阶项虽然可能改善样本内拟合,却容易在随访两端产生缺乏数据支持的弯曲。
候选模型比较通常涵盖组数、时间函数阶数、结局分布和方差约束等维度
[16-17]。赤池信息量准则(Akaike information criterion,AIC)和贝叶斯信息准则(Bayesian information criterion,BIC)用于权衡模型拟合与复杂度,比较须基于相同样本、结局定义和似然尺度。通常AIC、BIC越小越优;部分GBTM软件采用“对数似然减去惩罚项”的等价形式,方向相反,须先核实输出定义。
分类质量常用平均后验概率(average posterior probability,AvePP)、正确分类优势比(odds of correct classification,OCC)、模型估计组比例和最大后验分类比例报告。AvePP≥0.70、OCC≥5和组比例≥5%可作为经验参照
[11,18],模型选择还需结合类别规模、稳定性和可解释性。
有限混合似然可能存在多个局部极值,应采用不同初值重复拟合,并在标签对齐后比较似然值、参数及分类结果
[19]。除软件收敛标志外,必要时还应核查观测信息矩阵的正定性及标准误、置信区间的稳定性;若结果随初值改变,应报告其敏感性。
建模工具
R语言中,trajeR包适用于经典GBTM,支持CNORM、二分类logit、Poisson、零膨胀Poisson和Beta等分布,可计算个体后验概率及常用分类质量指标
[3];lcmm包则基于潜在类别混合效应框架,支持含随机效应或复杂链接函数的扩展模型
[19]。结果报告中须注明R与软件包版本、结局分布、时间变换、模型评价指标与各组样本量。
应用实例
数据来源与研究问题
本文使用lcmm包所附的PAQUID(Personnes Agées QUID)数据作为分析实例,使用trajeR包进行分析。PAQUID是一项基于法国西南部老年人群的前瞻性队列研究,旨在观察认知和功能老化及痴呆的发生
[20]。lcmm包所附的paquid数据为该队列的一个示例子样本,包含500名研究对象、2 250条随访记录,最长观察时间约20年,涵盖MMSE及其他认知与功能测量指标,并提供性别、教育和入组年龄等基线信息
[19]。
本例描述MMSE随年龄的轨迹,探索性别、教育和入组年龄与轨迹组归属的关系,以及轨迹组与随访期痴呆诊断的关联。示例子样本的组比例和关联估计不能直接外推至PAQUID总体或一般老年人群。
以下展示应用GBTM的标准四步流程:数据整理与模型设定→候选模型构建与选择→轨迹组确定与分类质量评价→相关因素探索性分析。
数据整理与模型设定
将MMSE在0~30分范围内且同时具有有效年龄记录的测量定义为有效观察,并纳入至少完成3次有效MMSE测量的个体。该纳入标准旨在提高本教学实例中轨迹形态估计和后验分类的稳定性,并非GBTM的一般性硬性要求。最终纳入342名研究对象,共1 977次有效观察,年龄范围为66.3~100.1岁;每人测量3~9次,中位数为5次。入组年龄为(73.3±5.9)岁,其中男性131人(38.3%),小学毕业252人(73.7%)。
以实际年龄作为时间尺度,并构建中心化变量age65=(age-65)/10,以65岁为时间原点、10岁为一个时间单位,以改善数值稳定性。CNORM模型的截距对应65岁时潜在正态变量的均值;观测尺度上的期望MMSE需结合量表上下界转换得到。将每名个体的MMSE和age65分别整理为结局矩阵Y和时间矩阵A,间断缺失位置保留为空值。MMSE具有0分和30分2个明确边界,故采用CNORM模型,并在工作模型中令各轨迹组共享残差标准差。
图1显示个体MMSE随年龄变化具有明显异质性;总体LOESS曲线呈非线性下降,但单一平滑曲线不能揭示潜在的不同变化模式。
候选模型构建与选择
组数初筛时分别拟合1~6组模型(
表2),各轨迹均采用二次多项式。在分类指标比较中,6组模型BIC最低,但6组模型的最小组比例仅3.5%,最低OCC为4.56,提示存在小类别和较弱的分类清晰度。3组模型较简约,各项分类指标均达到常用经验参考值。4组模型的最小组比例为9.1%,最低AvePP为0.789,最低OCC为4.95,略低于经验参考值5,但偏离幅度很小,且各组规模相对均衡。
为评估局部极值和初值敏感性,以原始初值及5个轻度扰动初值重复拟合4组模型,6次拟合的对数似然、最小模型估计组概率和共享残差标准差基本一致。进一步扩大初值范围后,4~6组模型均得到对数似然更高的其他解,但其最小组比例分别为3.8%、1.8%和1.5%,形成较小类别,提示该解可能反映对少数个体的进一步拆分,并对初值较为敏感。
综合模型拟合、分类质量、类别规模、数值稳定性及解释需要,本文未将单一指标作为机械判定标准,最终选择4组作为示范工作模型,并开展轨迹描述及后续探索性关联分析。这并不等于4组模型在统计意义上具有唯一正确性,而是针对本实例的综合权衡。
所选模型的轨迹及分类质量
按85岁时CNORM观测尺度上的模型期望MMSE由高至低重新排列4个轨迹组,并依据相对水平和变化速度命名。
图2显示模型拟合轨迹;90岁后的虚线表示数据支持较稀疏的区间,超过各组最大观察年龄的部分属于外推。
按最大后验概率分类,相对保持型207人(60.5%),70岁和90岁期望MMSE分别为28.43和26.65分;加速下降型62人(18.1%),约80岁后下降加快,90岁为19.38分;较低水平缓慢下降型42人(12.3%),70岁和90岁分别为24.32和22.04分;早发快速下降型31人(9.1%),70岁为27.13分,85岁和90岁分别降至16.55和6.58分。组名仅描述样本内轨迹,不对应临床诊断或病理亚型。
各组AvePP为0.789~0.947,OCC为4.95~144.26,其中相对保持型OCC接近经验参考值5。模型估计比例与后验分类比例并不完全相同,例如相对保持型分别为54.0%和60.5%,反映部分个体在相邻轨迹间仍存在分类不确定性。观察均值与拟合曲线在主要数据支持范围内总体一致,但高龄区间应谨慎解释:早发快速下降型90岁及以上仅有11次观察,最大观察年龄为94.3岁;因此曲线末端趋近0分主要反映模型外推和量表下界作用,不能解释为所有个体在相应年龄均达到0分。
轨迹组归属相关因素的探索性分析
在选定4组工作模型后,按最大后验概率分类描述各轨迹组的基线特征。各组入组年龄、性别和教育构成及未经调整的总体比较见
表3。
随后,以相对保持型为结局参照,采用多项Logistic回归分析性别、教育水平和入组年龄与轨迹组归属的关系。入组年龄按每增加5岁编码,模型同时纳入3个基线变量,结果见
表4。各预测因素的总体
P值采用似然比检验。
未经调整的组间比较显示,4组入组年龄、性别和教育构成不同(
表3)。多项Logistic回归中,仅教育水平的总体检验有统计学意义(
P<0.001);小学毕业者归入3个下降型而非相对保持型的优势比较低(
表4)。
上述关联不能作因果解释。PAQUID数据中的CEP(Certificat d’études primaires,法国小学毕业文凭)仅区分是否小学毕业,不能反映具体受教育年限或教育质量;教育水平也与MMSE表现相关
[20]。此外,最大后验分类后的回归未校正类别分类误差,结果仅用于方法演示和探索性解释。
轨迹组与随访期诊断为痴呆的关联探索
以入组年龄至诊断为痴呆或末次随访年龄的差值作为随访时间,将诊断为痴呆作为事件。
图3以Kaplan-Meier曲线估计描述各轨迹组的痴呆累积发生率,4组的痴呆累积发生率呈明显分离,组间差异具有统计学意义(log-rank
P<0.001)。其中,早发快速下降型上升最快,加速下降型和较低水平缓慢下降型次之,相对保持型最低,提示不同MMSE变化轨迹与随访期诊断为痴呆存在显著差异。
调整入组年龄、性别和教育水平后,与相对保持型相比,3个下降型均与更高的痴呆发生风险相关(
P均<0.001,
表5)。入组年龄的关联有统计学意义,性别和教育水平的关联未达统计学显著性。尤其需要注意,114例事件中有112例在诊断为痴呆当日或之后仍有MMSE观察参与轨迹估计,存在时间信息重叠和反向因果风险。因此,本分析仅用于演示探索性描述关联,不能用于评价轨迹对未来痴呆的预测效能,也不能作因果解释。
结语 医学纵向研究关注群体特征随时间的演变。总体平均轨迹可概括群体变化,但难以直接揭示不同亚组的变化模式。GBTM以有限条可解释的代表性轨迹描述人群异质性,为人群分层、识别高危亚组及探索相关因素提供直观框架。
实际应用中,研究者应根据研究问题确定时间尺度、结局分布和合理的时间函数,系统比较候选模型,并通过不同初值重复拟合,检查数值稳定性。模型选择需综合拟合优度、分类质量、类别规模和可解释性,同时报告模型估计组比例与后验分类比例,避免仅凭单一指标确定组数。
GBTM适用于多种结局分布,但其潜在类别是统计近似,不应直接视为天然存在的生物亚型。类别数量和轨迹形状受样本构成、多项式阶数、残差结构及初值影响。未来可结合生存模型、联合潜在类别模型及多轨迹模型,研究纵向变化与事件结局及多个相关指标的关系,并在适当模型框架下处理信息性脱落;也可采用样条或贝叶斯方法拓展非线性建模。进一步应用应关注轨迹的可重复性,并在后续关联分析中考虑分类误差。
国家重点研发计划(2023YFC2411403)