跳转至

Variable Selection for Illness‐Death Processes Under Dual Observation Schemes

作者: Xianwei Li, Liqun Diao, Richard J. Cook
来源: Statistics in Medicine
主题: 流行病学
相关性: 4/10
机构绿灯: University of Waterloo（US News 前 50，免分进入精读）
链接: https://doi.org/10.1002/sim.70608

一、领域脉络与小综述

这个方向是什么： 这个子方向要解决的根本统计问题是：在慢性病流行病学队列中，如何对多态生存过程（特别是 illness-death 过程）进行联合建模与变量选择，当数据面临"双重观测方案"（dual observation scheme）——即疾病进展时间受区间 censoring、死亡时间受右 censoring——时，如何基于不完整数据有效估计不同转移强度（transition intensity）中的回归系数，并区分不同转移的协变量效应进行特征筛选。当前该方向在建模与算法层面已有较成熟框架，但在区间 censoring 下的 penalized likelihood 变量选择及其渐近理论（如 oracle property、效率界）仍留有口子。

发展脉络（history）： 从 intro 引用的工作可串出以下线索： - 奠基工作：Cox (1972) 建立了 multiplicative intensity 模型与部分似然框架，为连续时间多态过程的回归建模奠定基础；Aalen (1978) 提出非参数累加强度模型，给出另一条路线。 - 多态过程与 illness-death 模型：Kalbfleisch & Prentice (2002) 在其生存分析专著中系统化了多态过程的似然构造；Cook & Lawless (2018) 在 Multistate Models with Interval-Censored Data 中将 illness-death 过程与区间 censoring 结合，给出了完整的数据构造与似然理论，本文作者正是该专著的作者之一，此引用为本文设定提供了直接基石。 - 变量选择与 penalized likelihood：Tibshirani (1996) 引入 Lasso（L1 penalty）做变量选择；Fan & Li (2001) 提出 SCAD penalty 并证明其 oracle property；Zou (2006) 提出 Adaptive Lasso 以修复 Lasso 的不一致性。这些工作主要在右 censoring 生存数据或完全观测下展开。 - 多态过程的变量选择：Ambroise & McLachlan (2002) 及 Shen & Huang (2006) 等将 penalized regression 引入多态模型，但多依赖右 censoring 或当前状态（current status）数据，未处理区间 censoring 与右 censoring 并存的双重方案。 - 本文的位置：在 Cook & Lawless (2018) 的区间 censoring 多态似然框架之上，引入 Fan & Li (2001) / Zou (2006) 的分层惩罚机制，针对 illness-death 过程的三条转移路径分别设定不同 penalty，并构造 EM 箯法优化此 penalized observed data likelihood。

子线索聚类： 1. 多态过程建模与似然构造（Cook & Lawless 2018, Kalbfleisch & Prentice 2002）：聚焦于区间 censoring 下多态过程的计数过程似然与马尔可夫/半马尔可夫假设下的数据生成机制。 2. Penalized likelihood 与变量选择理论（Fan & Li 2001, Zou 2006, Tibshirani 1996）：聚焦于 penalty 函数的数学性质（凸性、无偏性、oracle property）与算法实现（LARS、Newton-Raphson、坐标下降）。 3. EM 算法在 censoring 下的应用（Dempster et al. 1977, Louis 1982, Turnbull 1976）：聚焦于不完整数据下似然优化的 EM 框架及其收敛性、观测信息矩阵提取。

这个方向在追问的核心问题： 1. 在区间 censoring 与右 censoring 并存时，如何构造 illness-death 过程的 observed data likelihood 并使其可解？ 2. 如何对不同转移路径（健康→疾病、健康→死亡、疾病→死亡）的回归系数进行分层变量选择（transition-specific selection），而非全局同一 penalty？ 3. 如何在 penalized likelihood 框架下设计算法，使得 E-step 与 M-step 可借助现有软件（如 coxph、glmnet）实现，降低实现门槛？ 4. 当前瓶颈：区间 censoring 下 penalized likelihood 的渐近理论（oracle property、收敛率、效率界）尚未严格建立；EM 算法在分层 penalty 下的收敛性质缺乏理论保证。

⚠️ 作者的 framing： - 作者将缺口 frame 为：现有变量选择方法多针对右 censoring 或单一转移路径，而慢性病队列中疾病进展受区间 censoring、死亡受右 censoring 的"双重方案"下，缺乏针对 illness-death 过程的分层 penalized likelihood 方法。这让本文的 EM + 分层 penalty 成为"显然的下一步"。 - 被淡化的竞争路线：非参数/半参数累加强度模型（Aalen 模型）下的变量选择、基于 counting process 的 penalized partial likelihood（如 Fine & Gray 2009 对竞争风险的 Lasso）、以及基于 pseudo-observation 的变量选择路线。Intro 中未引用这些直接竞争方法。 - 明显该被引却未出现的：半参数效率界理论（如 Bickel et al. 1993 Efficient and Adaptive Estimation for Semiparametric Models）、区间 censoring下的无偏估计与 influence function 工作（如 Huang 1996 对 current status 的效率界）、以及 penalized M-estimation 的渐近理论（如 Fan & Peng 2004 对高维 SCAD 的 oracle property 推广）。这些是审视本文方法渐近性质的关键参照，值得研究者去查。

张力： 未见明显对立引用。Cook & Lawless (2018) 的似然框架与 Fan & Li (2001) 的 SCAD 理论在设定上互补而非矛盾，但存在一个隐含张力：Cook & Lawless 的似然依赖马尔可夫假设（转移强度仅依赖当前状态与时间），而 Fan & Li 的 oracle property 证明依赖独立同分布样本与正则条件——区间 censoring 下的多态过程数据既非 i.i.d.（随访时间相关），又受马尔可夫约束，二者的理论前提并不直接兼容，这为后续渐近理论推导埋下难点。

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

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

  • \(k\)：个体索引，\(k = 1, \ldots, n\)，样本量。
  • \(Y_k(t)\)：个体 \(k\) 在时间 \(t\) 所处状态，取值 \(\{1, 2, 3\}\)，1=健康，2=疾病，3=死亡。
  • \(T_{12,k}\)：个体 \(k\) 从健康（1）进展到疾病（2）的随机时间（潜在变量，不可直接观测）。
  • \(T_{13,k}\)：个体 \(k\) 从健康（1）直接死亡（3）的随机时间。
  • \(T_{23,k}\)：个体 \(k\) 从疾病（2）死亡（3）的随机时间。
  • \(X_k\)：个体 \(k\) 的 \(p\) 维协变量向量（固定基线协变量）。
  • \(\beta_{12}, \beta_{13}, \beta_{23}\)：三条转移路径 \(1\to2\), \(1\to3\), \(2\to3\) 对应的 \(p\) 维回归系数向量（要估的参数 / estimand）。
  • \(\lambda_{rs}(t | X_k)\)：从状态 \(r\) 到 \(s\) 的转移强度（hazard），模型为 multiplicative intensity：\(\lambda_{rs}(t | X_k) = \lambda_{rs,0}(t) \exp(X_k^\top \beta_{rs})\)，其中 \(\lambda_{rs,0}(t)\) 是基准强度（非参数/半参数部分）。
  • \(A_k(t)\)：个体 \(k\) 在时间 \(t\) 的随访评估时间序列，\(A_k(t) = \{a_{k0}, a_{k1}, \ldots, a_{kM_k}\}\)，其中 \(a_{k0}=0\)，\(a_{kM_k} \leq C_k\)。
  • \(C_k\)：个体 \(k\) 的右 censoring 时间（死亡时间的 censoring，如失访或研究结束）。
  • \(B_k\)：个体 \(k\) 的疾病进展区间，\(B_k = (L_k, R_k]\)，其中 \(L_k\) 是最后一次确认健康的时间，\(R_k\) 是首次确认疾病的时间（若从未确认疾病，则 \(R_k = \infty\)；若首次评估即疾病，则 \(L_k = 0\)）。
  • 可观测数据：对每个个体 \(k\)，实际观测到的是 \(\{X_k, A_k, B_k, T_{D,k}, \delta_{D,k}\}\)，其中：
  • \(T_{D,k} = \min(T_{13,k}, T_{23,k}, C_k)\) 为观测到的死亡或 censoring 时间，
  • \(\delta_{D,k} = I(T_{13,k} \leq C_k \text{ 或 } T_{23,k} \leq C_k)\) 为死亡指示变量，
  • \(B_k = (L_k, R_k]\) 为疾病进展的区间 censoring 区间，
  • \(A_k\) 为评估时间序列。
  • 不可观测的潜在量：\(T_{12,k}\)（疾病进展确切时间）被区间 censoring 遮蔽，仅知其落在 \((L_k, R_k]\)；若个体在 \(L_k\) 前死亡（\(T_{D,k} < L_k\)），则 \(T_{12,k}\) 完全不可观测（既不知是否进展，也不知进展时间）。

第二步：最小内核——最简特例：马尔可夫 illness-death 过程 + 区间 censoring + 分层 L1 penalty

剥掉一般性设定（多条转移、不同 penalty 函数、非参数基准强度），取最简特例： - 特例设定：\(p=1\)（单协变量），三条转移各有系数 \(\beta_{12}, \beta_{13}, \beta_{23}\)（各为标量），基准强度 \(\lambda_{rs,0}(t)\) 取常数（参数模型），penalty 取 L1（\(\|\beta_{rs}\|_1\)），马尔可夫假设成立（转移强度仅依赖当前状态）。 - 要估的参数：\(\theta = (\lambda_{12,0}, \lambda_{13,0}, \lambda_{23,0}, \beta_{12}, \beta_{13}, \beta_{23})\)。 - Observed data likelihood（个体 \(k\)）： 若 \(\delta_{D,k}=1\)（观测到死亡）且死亡前确诊疾病（\(R_k < T_{D,k}\)），则似然因子为：

\[L_k(\theta) = P(Y_k(L_k)=1) \cdot P(Y_k(R_k)=2 | Y_k(L_k)=1) \cdot P(T_{D,k} \in [t, t+dt], Y_k(t)=2 | Y_k(R_k)=2)\]
在马尔可夫假设下，这些概率可由转移强度 \(\lambda_{rs}\) 的指数积分算出（如 \(P(Y_k(R_k)=2 | Y_k(L_k)=1) = \int_{L_k}^{R_k} \lambda_{12}(s|X_k) \exp\{-\int_{L_k}^s [\lambda_{12}(u|X_k)+\lambda_{13}(u|X_k)]du\} ds\)）。 若死亡前未确诊疾病（\(T_{D,k} < L_k\) 或 \(L_k < T_{D,k} < R_k\)），则似然需积分掉 \(T_{12,k}\) 的不确定性。 - Penalized likelihood：
\[Q(\theta) = \log L(\theta) - n(\gamma_{12}|\beta_{12}| + \gamma_{13}|\beta_{13}| + \gamma_{23}|\beta_{23}|)\]
其中 \(\gamma_{rs}\) 是分层 penalty 系数（不同转移用不同权重）。 - 核心数学困难：\(T_{12,k}\) 被区间 censoring 遮蔽，observed data likelihood 中包含对 \(T_{12,k}\) 的积分（或条件概率的指数积分），导致 \(Q(\theta)\) 无法直接对 \(\beta_{rs}\) 求导并套用坐标下降/LARS——因为似然不是 Cox 部分似然那样的线性累加形式，而是状态转移概率的复杂函数。 - 本文的关键想法：用 EM 算法绕过积分困难——E-step 将 \(T_{12,k}\) 的条件分布（给定观测数据与当前参数估计）离散化到评估时间点上，计算"伪个体"的期望计数过程数据（期望转移次数、期望风险集大小）；M-step 将这些期望数据代入 Cox 部分似然形式，此时 penalized M-step 退化为分层 penalized Cox 回归（不同转移路径的系数用不同 penalty），可直接用现有软件（如 coxph + penalty 或 glmnet 的分层 Lasso）求解。 - 为什么成立：马尔可夫假设下，给定 \(T_{12,k}\) 的确切时间，多态过程的似然退化为三条独立 Cox 部分似然的乘积（每条转移一个）；EM 的 E-step 补全了 \(T_{12,k}\) 的期望，M-step 就回到了标准 penalized Cox 框架。分层 penalty 的引入不影响 EM 的收敛性框架（因为 M-step 仍是一个凸优化），但不同 penalty 权重 \(\gamma_{rs}\) 的选择需外层调参（如 BIC/CV）。

三、这篇论文做了什么

三句话： ①研究了 illness-death 过程在双重观测方案（疾病进展区间 censoring + 死亡右 censoring）下的变量选择问题； ②核心工具是分层 penalized observed data likelihood + EM 箯法（E-step 离散化条件分布补全期望数据，M-step 分层 penalized Cox 回归）； ③主要结论是该方法可借助现有软件实现，模拟显示有限样本下变量选择一致性较好，NACC 痴呆数据应用揭示了不同转移路径的协变量差异。

关键设定与假设： - 马尔可夫假设：转移强度 \(\lambda_{rs}(t|X_k)\) 仅依赖当前状态 \(r\) 与时间 \(t\)，不依赖历史路径（如 \(T_{12,k}\) 的具体值不影响 \(\lambda_{23}\)）。这是似然可分解为 Cox 形式的基础，也是 EM 补全数据后 M-step 可用 Cox 软件的必要条件。相比 Cook & Lawless (2018) 的半马尔可夫扩展，本文强化了马尔可夫假设以简化算法。 - Multiplicative intensity 模型：\(\lambda_{rs}(t|X_k) = \lambda_{rs,0}(t) \exp(X_k^\top \beta_{rs})\)，每条转移有独立的基准强度 \(\lambda_{rs,0}(t)\) 与系数 \(\beta_{rs}\)。这是 Cox 模型的直接推广，相比 Aalen 累加模型限制了效应的乘法结构。 - 独立 censoring 假设：评估时间序列 \(A_k\) 与右 censoring \(C_k\) 的生成机制独立于多态过程 \(Y_k(t)\) 的转移时间（给定协变量）。这是 observed data likelihood 可忽略评估机制的条件。 - 分层 penalty 设定：不同转移路径的系数 \(\beta_{12}, \beta_{13}, \beta_{23}\) 使用不同 penalty 函数（如 \(\beta_{12}\) 用 SCAD，\(\beta_{13}\) 用 Adaptive Lasso，\(\beta_{23}\) 用 L1）或不同权重 \(\gamma_{rs}\)。这允许不同转移有不同的稀疏度与变量集，是本文与全局同一 penalty 的关键区别。

主要结果： 1. Penalized observed data likelihood 的构造（Section 2-3）： - 在马尔可夫假设与独立 censoring 下，个体 \(k\) 的 observed data likelihood 可写为状态转移概率 \(P(Y_k(a_{kj})=s | Y_k(a_{kj-1})=r, X_k)\) 的乘积（在评估时间点序列上），这些概率由转移强度的指数积分算出。 - Penalized 目标函数：\(Q(\theta) = \log L(\theta) - n \sum_{r<s} p_{\gamma_{rs}}(\beta_{rs})\)，其中 \(p_{\gamma_{rs}}\) 是 penalty 函数（SCAD/Adaptive Lasso/L1），\(\gamma_{rs}\) 是转移特定的 tuning 参数。 - 统计含义：分层 penalty 允许健康→疾病的协变量（如认知测试分数）与健康→死亡的协变量（如心血管风险因子）有不同的稀疏结构，避免全局 penalty 强制同一稀疏度。

  1. EM 箯法的构造与实现（Section 4，核心算法结果）：
  2. E-step：给定当前参数估计 \(\theta^{(m)}\)，计算 \(T_{12,k}\) 在每个评估区间 \((a_{kj-1}, a_{kj}]\) 内的条件概率 \(w_{kj}^{(m)} = P(T_{12,k} \in (a_{kj-1}, a_{kj}] | \text{observed data}, \theta^{(m)})\)，以及个体在死亡时所处状态的条件概率。这些权重通过马尔可夫转移概率的递推计算（前向-后向算法类似物）得到。
  3. M-step：将 E-step 的权重 \(w_{kj}^{(m)}\) 视为"伪个体"的期望转移次数，构造分层 penalized Cox 部分似然：
    \[Q^{(m)}(\theta) = \sum_{r 其中 \(w_{rs,kj}^{(m)}\) 是期望计数，\(R_{rs,kj}\) 是期望风险集。此目标函数可按转移路径分层优化，每层用现有 penalized Cox 软件（如 R 的 coxph + penalty 参数）求解。

  4. 收敛性：EM 框架保证 \(Q(\theta)\) 单调上升，但分层 penalty 的非凸性（SCAD）可能使 M-step 有多个局部极值；作者未严格证明全局收敛，仅依赖模拟观察收敛。

  5. Tuning 参数选择（Section 5）：

  6. 用修正 BIC（考虑区间 censoring 的有效样本量）选择 \(\gamma_{rs}\)：\(\text{BIC} = -2 \log L(\hat{\theta}) + \log(n_{\text{eff}}) \cdot df\)，其中 \(n_{\text{eff}}\) 是期望事件数，\(df\) 是非零系数数。
  7. 模拟显示 BIC 在区间 censoring 下比 AIC 更稳定，但有效样本量的定义依赖 E-step 的期望计数，理论性质未严格建立。

证明路线与技术技巧： - 整体路线： 1. 构造 observed data likelihood（基于马尔可夫转移概率的指数积分）→ 2. 引入分层 penalty 构造 penalized 目标函数 → 3. 识别直接优化的困难（似然中的积分不可解析求导）→ 4. 构造 EM 算法：E-step 计算条件概率权重（离散化到评估时间点），M-step 转化为分层 penalized Cox 回归 → 5. 用修正 BIC 选 tuning 参数 → 6. 模拟验证有限样本表现。 - 关键跳跃点： - E-step 的条件概率计算：\(P(T_{12,k} \in (a_{kj-1}, a_{kj}] | \text{observed data}, \theta^{(m)})\) 的递推计算是核心难点——观测数据包含区间 censoring（\(B_k\)）与右 censoring（\(T_{D,k}, \delta_{D,k}\)），条件概率需在马尔可夫假设下用转移概率的乘积积分与贝叶斯规则递推。作者用前向算法（类似隐马尔可夫模型的概率计算）逐评估点递推，避免了直接积分。 - M-step 的分层 penalized Cox 构造：将期望计数 \(w_{rs,kj}^{(m)}\) 代入 Cox 部分似然形式，使得 M-step 可用现有软件求解——这是算法可行性的关键跳跃，因为直接优化 penalized observed data likelihood 需自定义数值方法（如 Newton-Raphson + 数值积分），而 EM + 分层 Cox 回归可复用 coxph 的成熟实现。 - 技术技巧点名： - EM 箯法（Dempster et al. 1977）：用于处理 \(T_{12,k}\) 的缺失数据，将 observed data likelihood 优化转化为完整数据似然的期望优化。 - 前向-后向算法（类似 HMM 的概率递推）：用于 E-step 中计算 \(T_{12,k}\) 在各评估区间的条件概率，避免直接积分。 - 分层 penalized Cox 回归：M-step 中不同转移路径的系数用不同 penalty 函数/权重优化，复用 coxph 的 penalty 参数或 glmnet 的分层坐标下降。 - 修正 BIC：用期望事件数 \(n_{\text{eff}}\) 替代样本量 \(n\)，适应区间 censoring 下的有效信息量。

真实例子与应用： - 数据：National Alzheimer's Coordinating Center (NACC) 的统一数据集（UDS），包含 2005-2020 年的纵向评估数据，约 30,000+ 个体，评估间隔约 1 年，记录认知状态（正常/MCI/痴呆）与死亡。 - 如何用上去：将认知状态简化为健康（正常/MCI）与疾病（痴呆），构建 illness-death 过程（健康→痴呆、健康→死亡、痴呆→死亡），协变量包括年龄、性别、教育、APOE4 基因、认知测试分数等。疾病进展时间受区间 censoring（每年评估一次），死亡时间受右 censoring（失访或研究结束）。 - 得到什么结果： - 变量选择结果显示：健康→痴呆的系数选中了 APOE4、认知测试分数（逻辑记忆、动物命名）、年龄；健康→死亡的系数选中了年龄、性别、心血管风险因子；痴呆→死亡的系数选中了年龄、基线认知分数。不同转移的选中变量集不同，验证了分层 penalty 的必要性。 - 估计的系数方向与流行病学文献一致（APOE4 增加痴呆风险、年龄增加死亡风险）。 - 这个例子想说明什么：展示分层 penalty 可揭示不同转移路径的协变量差异（全局 penalty 会强制同一稀疏度，可能漏选或误选），以及 EM 箯法在真实区间 censoring 数据上的可行性。

🔎 结论是否比证明窄： - Oracle property 的缺失：文中 claim 分层 penalty 可实现变量选择一致性（模拟显示），但未在定理中严格证明 SCAD/Adaptive Lasso 在区间 censoring 多态过程下的 oracle property——这是条件 X（马尔可夫 + 独立 censoring + penalty 函数性质）下应可证明但被泛泛 claim 的结论，扎根在 Section 6 的模拟讨论与 Section 7 的应用分析中。 - EM 算法的收敛速率：文中 claim EM 收敛（单调上升），但未给出收敛速率（如线性收敛的条件）或有限步后的误差界——这是 Dempster et al. (1977) 框架下的已知困难，本文回避了。 - 修正 BIC 的理论保证：BIC 的有效样本量 \(n_{\text{eff}}\) 定义依赖 E-step 的期望计数，其模型选择一致性未严格证明，仅在模拟中验证。

四、开放问题（点到为止，扎根具体语句）

  1. Oracle property 的严格证明：在马尔可夫 illness-death 过程 + 区间 censoring + 分层 SCAD/Adaptive Lasso penalty 下，证明 penalized estimator 的 oracle property（变量选择一致性 + 估计的渐近正态性）。扎根在 Section 6 末尾的讨论："Simulation studies demonstrate the finite sample performance of the method"——模拟验证了但理论未建立，需引用 Fan & Li (2001) 的 oracle property 证明路线并适配区间 censoring 似然的非标准结构。

  2. 半马尔可夫扩展下的 EM 算法：若放宽马尔可夫假设（\(\lambda_{23}\) 依赖 \(T_{12,k}\) 的具体值），E-step 的条件概率计算需联合 \(T_{12,k}\) 与 \(T_{23,k}\)，M-step 的 Cox 部分似然分解不再成立。扎根在 Section 2 的假设陈述："We assume a Markov illness-death process"——这是为算法可行性做的简化，放宽后需重新设计 M-step。

  3. Penalized likelihood 的渐近效率界：在区间 censoring + 右 censoring 下，未 penalized 的 observed data estimator 的半参数效率界是什么？penalized estimator 是否达到此界（在非零系数子模型上）？扎根在 intro 对 Fan & Li (2001) 的引用——他们证明了 SCAD 的 oracle property 包含渐近正态性，但未讨论效率界；本文亦未触及，需引用 Bickel et al. (1993) 的半参数效率理论。

  4. 高维设定下的变量选择：本文模拟的 \(p\) 较小（\(p \leq 10\)），当 \(p \gg n\) 时，分层 penalized Cox 的 M-step 是否仍可行（coxph + penalty 在高维下不稳定），是否需换用 glmnet 的坐标下降？扎根在 Section 5 的 tuning 参数选择讨论——仅考虑 \(p\) 较小的 BIC，未讨论高维下的 tuning 策略（如 CV-1se 或 stability selection）。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论