在期刊《Statistics in Medicine》2024年第43卷第5748-5764页，我们关注到一篇由美国密歇根大学（University of Michigan）的研究团队发表的重要原创性研究。该研究的第一作者兼通讯作者是生物统计学系的Mingyan Yu，共同作者包括同系的Zhenke Wu和Michael R. Elliott教授，以及来自社会研究所（Institute for Social Research）的Margaret Hicken教授。这项工作得到了美国国立卫生研究院（NIH）两项基金（R56AG-066693和R01MD-016244）的支持，旨在解决一个日益重要的方法论挑战：如何利用高频率、高强度的纵向生物标志物数据中的“变异性”信息来预测健康结局。

研究的学术背景与目标

在过去的几十年里，利用纵向测量数据（如生物标志物、调查问卷评分）来预测健康结局的研究已形成大量文献。传统方法（如线性混合效应模型与生存模型的联合建模）主要关注生物标志物轨迹的“均值”，并将其与结局相关联，而将“方差”（variability）视为需要被控制的干扰参数。然而，越来越多的证据表明，个体内部生物标志物的变异性本身可能蕴含着与健康结局密切相关的宝贵信息。例如，情绪变异性与自杀意念、睡眠时间变异性与抑郁症状、认知功能变异性与痴呆症发病风险之间均存在显著关联。这一认识催生了一个虽小但快速增长的、专注于探索变异性预测价值的研究领域。

随着可穿戴设备等技术的发展，生态瞬时评估（Ecological Momentary Assessment, EMA）数据或高强度纵向生物标志物数据变得越来越普遍。这类数据采集频率极高（如每秒一次的心率），导致每个个体在数小时或数天内即可产生成百上千个观测值。这为研究短期效应和更精准地预测长期结局提供了可能，但也对统计分析方法提出了特殊要求。现有研究在处理EMA数据与结局的关系时，大多仍聚焦于均值结构，或采用两阶段法（先估计变异性，再将其作为预测因子纳入结局模型），后者会因未能正确传导第一阶段估计的不确定性而导致有偏的推断和不可控的方差。

具体到心率数据，传统的心率变异性（Heart Rate Variability, HRV）分析通常基于心搏间期计算时间域或频域指标，并在特定时间段内（如5分钟）进行平均，得到一个单一的HRV汇总值。本研究提出了一种不同的、据信是新颖的方法，旨在从标准的HRV分析中挖掘出更多信息。

因此，本研究的核心目标在于提出并验证一个新颖的贝叶斯分层联合模型，用于将高强度纵向生物标志物的变异性作为预测因子，与横截面健康结局进行联合建模。研究的创新之处在于：1）同时建模生物标志物的残差变异性（“秒到秒”短期波动）和随机效应变异性（“分钟到分钟”长期趋势波动）；2）通过主体层面的三次B样条和跨个体的信息共享来处理数据的“高强度”特性；3）采用完全贝叶斯联合建模框架，而非两阶段方法，以确保推断的无偏性和有效性。研究者通过一个涉及社交压力研究的、具有赫兹级分辨率心率数据的生物监测应用，展示了所提模型的实用性。

详细的工作流程

本研究的工作流程严谨而系统，主要包括方法学模型构建、模拟研究验证和实际数据应用三个核心部分。

第一部分：统计模型构建 研究者构建了一个由两个子模型组成的联合模型框架。

1. 预测因子子模型（纵向生物标志物模型）：用于建模高强度纵向生物标志物数据 (X{ij})。核心是使用三次B样条基函数来拟合随时间变化的非线性的均值轨迹。模型形式为 (X{ij} = f(t_{ij}; \beta_0, \mathbf{b}i) + \epsilon{ij})。其中，(f(\cdot)) 是由总体固定截距 (\beta_0) 和个体特定随机效应 (\mathbf{b}_i) 定义的样条函数。(\mathbf{b}i) 被假设服从均值为零、方差为 (\sigma^2{bi}) 的正态分布，此方差被定义为个体的“长期变异性”，反映了心率分钟级别的波动趋势（即样条曲线的“弯曲度”）。残差项 (\epsilon{ij}) 则捕捉“短期变异性”，反映了秒级别的瞬时波动。为了处理高频率数据中必然存在的自相关，研究者假设残差具有一阶自回归结构，即 (\epsilon_{ij} = \rhoi \epsilon{i(j-1)} + w{ij})，其中 (w{ij}) 为独立同分布的白噪声，方差为 (\sigma^2_{wi})。因此，(\sigma^2{wi}) 被定义为个体的“短期变异性”。模型通过设置分层先验（如 (\log(\sigma^2{b_i}) \sim N(v_b, \psi^2_b))），实现了不同个体间在长期和短期变异性上的信息共享。

2. 结局子模型（二分类结局模型）：用于建模横截面二分类结局 (Y_i)。研究者采用了一个概率单位回归模型：(P(Y_i=1) = \Phi(\eta_i))，其中 (\Phi) 为标准正态累积分布函数，线性预测项 (\eta_i) 的关键部分由个体特异性变异性的对数构成，即 (\eta_i = \alpha_1 + \alpha2 \times \log(\sigma{w_i}) + \alpha3 \times \log(\sigma{b_i}) + \boldsymbol{\gamma} \mathbf{Z}_i)。这样，从纵向子模型提取出的长期和短期变异性被直接作为预测因子纳入了结局模型。

3. 贝叶斯推断与实现：由于模型包含了多个非共轭先验，后验分布没有闭合解。因此，研究采用马尔可夫链蒙特卡洛方法进行后验推断，具体使用JAGS软件并通过R语言中的r2jags包进行调用。JAGS能够识别纵向子模型中的广义线性混合模型结构，从而对固定效应和随机效应进行块更新，显著提高了抽样效率。

第二部分：模拟研究 在将模型应用于实际数据之前，研究者进行了深入的模拟研究，有两个主要目的：一是检验所提模型恢复真实参数值的能力；二是与常见的两阶段方法进行比较。

1. 模拟设置：模拟了N=150名个体，每人有(ni)=600个观测值，进行了R=200次重复。模拟参数基于从实际应用数据集（密歇根工作-生活研究）中随机选取的40名参与者的心率数据拟合得到，以确保模拟数据具有真实心率数据的特征。研究者设置了两种场景来考察测量误差大小的影响：场景1（低测量误差） 和 场景2（高测量误差），主要通过调整白噪声方差(\sigma^2{wi})的先验分布均值来实现。在每种场景下，首先根据纵向子模型参数生成模拟心率数据，然后基于生成的真实(\sigma{bi})和(\sigma{w_i})，按照设定的回归系数(\boldsymbol{\alpha})生成二分类结局。

2. 对比方法：作为对比的两阶段方法是：第一阶段，使用R语言nlme包中的lme()函数，对所有人的模拟心率数据拟合一个具有AR(1)残差结构的线性混合效应模型，然后通过计算随机效应的经验方差和残差白噪声的经验方差，分别得到每个个体的长期和短期变异性估计值；第二阶段，将这些估计值作为已知预测因子，使用glm()函数拟合一个概率单位回归模型来预测结局。

3. 评价指标：通过覆盖率、偏差、平均区间长度和均方根误差四个指标来评估和比较两种方法的性能。

第三部分：实际应用（密歇根工作-生活研究） 研究者将提出的联合模型应用于“密歇根工作-生活研究”的实证数据中，以探讨赫兹级心率变异性与社会压力（通过实验分组反映）之间的关联。

1. 研究描述：该研究是一项随机对照实验，旨在调查与潜在歧视和偏见情境相关的预期性和持续性压力。232名健康的年轻黑人女性被随机分配到四个实验组之一，接受为期2小时的实验室社交压力测试（基于特里尔社会压力测试，TSST）。四个组分别是：组1（WH：白人敌对观众）、组2（WF：白人友好观众）、组3（BF：黑人友好观众）和组4（C：对照组，无观众）。整个实验过程连续采集心率数据。在排除心率数据大量缺失的参与者后，最终分析样本为165人。

2. 数据分析框架调整：

   • 纵向子模型扩展：为了更精细地刻画压力反应过程，研究者将实验时间段（从标记1到标记6）划分为五个时期：预期期、第一次测量间歇期、压力期、第二次测量间歇期、恢复期。并假设长期变异性(\sigma^2_{bi})和短期变异性(\sigma^2{w_i})在这五个时期是不同的，即为每个时期估计一组变异参数。
   • 结局子模型扩展：研究关注的是组别分配这个二分类结局，但进行了多组两两比较（如WH vs WF， WF vs BF， WH vs C等），以对应不同的科学假设。每次比较拟合一个独立的联合模型。线性预测项(\eta_i)扩展为包含五个时期的长期和短期变异性对数，共11个回归系数。

3. 模型拟合与诊断：对每个配对比较模型，运行3条链，每条链18000次迭代（前9000次作为退火）。使用Gelman-Rubin统计量评估收敛性，所有参数R值均小于1.1，表明收敛良好。此外，研究者进行了后验预测检查，通过计算卡方差异统计量的后验预测p值来评估模型对纵向心率数据和二分类结局的拟合优度，结果均显示模型拟合良好。

主要结果

模拟研究结果清晰地证明了联合建模方法的优越性。 在场景1（低测量误差） 下，联合模型与两阶段方法的表现相近，联合模型的参数覆盖率在94%-96.5%之间，偏差很小。两阶段方法虽覆盖率略低（92%-94%），但偏差也较小。 然而，在场景2（高测量误差） 下，联合模型的优势变得非常明显。联合模型依然保持了优异的性能，覆盖率在92.5%-95.5%之间，偏差可控。相反，两阶段方法彻底失效：对于参数(\alpha_1)、(\alpha_2)和(\alpha_3)，其覆盖率分别暴跌至11.5%、24%和61%，并且出现了严重的偏差（例如(\alpha_2)的偏差达1.01）。这充分验证了理论预期：当纵向标记的测量误差较大时，两阶段方法因无法将第一阶段估计的不确定性传递到第二阶段，会导致有偏且不可靠的推断。而联合模型通过同时估计所有参数，恰当地考虑了这种不确定性，因此无论测量误差水平如何，都能提供有效推断。

实际应用结果揭示了心率变异性与社会压力之间复杂而有趣的关系模式。 研究者根据TSST和压力生物学文献，对不同实验组在各个时期的压力反应差异提出了具体假设。应用模型的结果表明，短期（秒到秒）变异性和长期（分钟到分钟）变异性可能捕捉了压力反应过程的不同方面。 * 长期变异性的模式与传统HRV研究的发现更一致。通常，更高的平均HRV（对应于更低的心率和更低的急性压力）与更低的压力相关。在本研究中，对照组（理论上压力暴露最低）在多个时期表现出显著更高的长期变异性。例如，与BF组相比，C组在预期期和压力期的长期变异性更高；与WF组相比，C组在两次测量间歇期的长期变异性更高；与WH组相比，C组在第一次测量间歇期的长期变异性更高。这支持了“观众存在本身即构成压力源”的假设。然而，这种模式并非绝对。例如，在压力期，WF组的长期变异性反而高于C组，这可能意味着与一位不苟言笑的黑人工作人员互动所产生的压力，与面对一个不苟言笑的白人友好观众团体进行演讲的压力水平相似。 * 短期变异性则可能呈现出与压力体验的正相关关系。例如，在压力期，暴露于白人敌对观众（WH）的参与者，其短期变异性有高于暴露于白人友好观众（WF）参与者的趋势（系数为负但95%置信区间包含0）。更重要的是，在那些观众组与对照组存在差异的情况下，进行演讲的观众组通常表现出更高的短期变异性。例如，在预期期和第二次测量间歇期，WF组和BF组的短期变异性显著高于C组。这提示短期变异性可能是急性压力反应的一个敏感指标。

后验预测检查的结果为模型的有效性提供了支持。对于所有五个配对比较模型，纵向心率数据的后验预测p值中位数均在0.36左右，二分类结局的后验预测p值在0.34-0.46之间，表明模型生成的复制数据与观测数据吻合良好，模型设定是合理的。

结论与价值

本研究的主要结论是，成功开发并验证了一个贝叶斯分层联合模型，能够有效地将高强度纵向生物标志物（如心率）中不同时间尺度的变异性提取出来，并将其作为预测因子与健康结局进行联合建模。该方法在方法论上优于传统的两阶段方法，特别是在数据存在较大测量误差时，能保证推断的无偏性和有效性。

其科学价值体现在两方面：方法论上，它将联合建模技术扩展到了高强度纵向数据场景，并创新性地同时利用残差和随机效应两种变异性作为预测因子，为分析密集采样数据提供了新的强大工具。应用上，通过对密歇根工作-生活研究的分析，揭示了心率变异性的不同成分（秒级与分钟级）可能分别对应压力反应的不同生理心理过程，深化了我们对社交压力下自主神经系统反应的理解，为利用可穿戴设备数据进行精准的心理生理状态评估提供了新思路。

研究亮点

1. 方法创新性：提出了首个针对高强度纵向生物标志物数据的、同时包含残差和随机效应变异性作为预测因子的贝叶斯联合模型。通过B样条和分层先验巧妙处理了高维度、高自相关的数据挑战。
2. 理论验证扎实：通过精心设计的模拟研究，有力证明了联合建模相对于两阶段方法在处理测量误差方面的绝对优势，为方法的选择提供了实证依据。
3. 应用发现新颖：在实际应用中，区分并解释了心率中“短期”与“长期”变异性与社交压力关联的不同模式，这一发现超越了传统HRV平均值的分析，提供了更精细的洞察。
4. 解决实际问题：直接针对可穿戴设备产生的海量高频数据进行分析，具有很高的现实应用价值，为行为医学、健康心理学、流行病学等领域利用此类数据预测健康结局树立了方法学典范。

其他有价值内容

研究者在讨论部分也坦诚指出了本工作的若干局限性，如样本量相对较小（尤其在对照组）、对自相关结构做了较强的AR(1)假设、模型计算耗时较长等。同时，他们也提出了有前景的未来扩展方向，例如考虑生物标志物的其他分布形式、在模型中纳入更多生物信号（如皮电活动）并处理其相关性、以及将个体特异性变异参数进一步建模为时间的连续函数以更动态地捕捉变化等。这些思考为后续研究指明了道路。