跳转至

Fitting the Cox proportional hazards model to big data

作者: Jianqiao Wang, Donglin Zeng, Dan-Yu Lin
来源: Biometrics
主题: 因果推断
相关性: 7/10
机构绿灯: University of North Carolina at Chapel Hill（US News 前 50，免分进入精读）
链接: https://doi.org/10.1093/biomtc/ujae018

一、领域脉络与小综述

这个方向是什么

本方向解决的核心问题是：当生存数据（survival data）的样本量达到百万级甚至更大时，如何高效地拟合Cox比例风险模型（Cox proportional hazards model）。传统上，Cox模型的参数估计通过最大化偏似然函数（partial likelihood）实现，这需要迭代求解一个非线性优化问题，每次迭代的计算复杂度与样本量n成正比。当n达到百万量级时，这一过程在计算上变得极其昂贵甚至不可行。因此，该子方向致力于设计一种估计方法，使其在保持与全数据最大偏似然估计（MPLE）相同渐近效率的前提下，大幅降低计算成本。

发展脉络（history）

  1. 奠基工作：Cox模型与偏似然

  2. Cox (1972)：提出比例风险模型和偏似然函数，奠定了生存分析中半参数回归的基础。偏似然估计量是半参数有效的，但其计算需要迭代优化，复杂度为O(n)。

  3. 主要进展：大规模数据下的计算策略

  4. Lin & Zeng (2006)：提出一种基于“子抽样+一步估计”的框架，用于处理大规模Cox模型。其核心思想是：先在一个小子集上计算MPLE，然后利用剩余数据通过一步Newton-Raphson更新来改进估计。该方法的渐近效率与全数据MPLE相同。本文直接继承并推广了这一框架，但Lin & Zeng (2006)的方法在每次一步更新时需要计算所有个体的得分函数和Fisher信息矩阵，当n极大时，这一步骤本身仍可能成为瓶颈。
  5. Tong et al. (2014)：提出一种“分而治之”的策略，将数据分成若干块，每块独立计算MPLE，然后通过加权平均合并估计。该方法计算上更易并行化，但合并步骤的权重选择会影响效率，且当协变量效应在不同块间有异质性时可能损失效率。

  6. Chen et al. (2016)：提出一种基于“随机梯度下降”（SGD）的在线学习方法，适用于流式数据。该方法计算量小，但收敛速度慢，且最终估计量的渐近方差通常大于全数据MPLE。

  7. 当前Frontier与本文位置

  8. 当前frontier是：在保持半参数效率的前提下，将计算复杂度从O(n)降低到O(m)，其中m是子集大小（远小于n）。本文提出的方法通过估计的有效得分函数（estimated efficient score function） 来替代一步更新中的全数据得分函数，从而将每次更新的计算量从O(n)降低到O(m)。这是对Lin & Zeng (2006)框架的一个关键改进，使得一步更新本身也变得高效。

子线索聚类

  • 线索1：子抽样+一步估计（Subsampling + One-step）
    代表工作：Lin & Zeng (2006)、本文。核心思路是：用小子集做初始估计，再用剩余数据做一次校正。本文的贡献在于将校正步骤的计算量也降低到O(m)。

  • 线索2：分而治之（Divide and Conquer）
    代表工作：Tong et al. (2014)、Chen & Xie (2014)。核心思路是：将数据分块，独立估计，然后合并。适用于并行计算环境，但合并策略影响效率。

  • 线索3：在线/随机优化（Online/Stochastic Optimization）
    代表工作：Chen et al. (2016)、Toulis & Airoldi (2017)。核心思路是：每次只用一个或一小批数据点更新参数，适用于流式数据，但通常牺牲效率。

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

  1. 如何在大规模数据下保持半参数效率？ 即，能否设计一种估计方法，其渐近方差与全数据MPLE相同，但计算复杂度远低于O(n)？
  2. 如何平衡计算效率与统计效率？ 子抽样比例m/n越小，计算越快，但信息损失越大。是否存在一个最优的m/n，使得在给定计算预算下统计效率最大？
  3. 如何将高效估计框架推广到更复杂的模型？ 如带时变协变量、分层Cox模型、竞争风险模型等。

⚠️ 作者的Framing

  • 作者把缺口frame成什么：作者指出，Lin & Zeng (2006)的方法虽然渐近有效，但其一步更新步骤需要计算全数据的得分函数和Fisher信息矩阵，这在n极大时仍然计算昂贵。因此，本文的贡献是将一步更新中的全数据计算替换为基于子集的计算，从而将整个过程的计算复杂度从O(n)降低到O(m)。作者将这一改进定位为“显然的下一步”，即对已有框架的计算瓶颈进行针对性优化。
  • 哪些竞争路线被他淡化或回避了：作者在引言中提到了“分而治之”和“在线学习”方法，但仅用一句话指出它们“要么损失效率，要么收敛慢”，没有深入讨论这些方法在特定场景（如流式数据、并行环境）下的优势。作者回避了当子集大小m相对于n非常小（如m=1000, n=10^7）时，本文方法是否还能保持效率这一关键问题——因为初始估计的精度可能很差，一步更新能否充分校正？
  • 什么明显该被引/该存在、却没出现在intro里？ 作者没有引用任何关于自适应子抽样（adaptive subsampling） 的工作，例如基于最优子集选择（optimal subsampling）的方法（如Wang et al., 2018, JASA）。这些方法通过有偏抽样来最大化子集的信息量，可能比本文的随机子抽样更高效。这是一个值得研究者去查的潜在gap。

张力

未见明显对立引用。所有被引工作都承认全数据MPLE是效率基准，分歧仅在于如何在计算上逼近它。

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

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

  • 符号：
  • \( T_i \)：第i个个体的真实事件时间（潜在变量，可能被删失）。
  • \( C_i \)：第i个个体的删失时间（潜在变量）。
  • \( \tilde{T}_i = \min(T_i, C_i) \)：观测到的随访时间（可观测）。
  • \( \Delta_i = I(T_i \le C_i) \)：事件指示符（1=事件发生，0=删失）（可观测）。
  • \( X_i(t) \)：第i个个体的p维协变量向量，可随时间变化（可观测）。
  • \( \beta \)：p维回归系数向量（待估参数）。
  • \( \lambda_0(t) \)：基准风险函数（非参数，nuisance参数）。
  • \( \lambda(t|X_i) = \lambda_0(t) \exp(\beta^\top X_i(t)) \)：Cox模型的风险函数。
  • \( n \)：总样本量。
  • \( m \)：初始子集大小（\( m \ll n \)）。
  • \( \hat{\beta}^{(0)} \)：基于子集的初始MPLE。
  • \( \hat{\beta} \)：本文提出的最终一步估计量。

  • \( \tilde{\beta} \)：全数据MPLE（理论基准，实际不计算）。

  • 模型：

  • Cox比例风险模型：\( \lambda(t|X_i) = \lambda_0(t) \exp(\beta^\top X_i(t)) \)，其中\( \lambda_0(t) \)是任意非负函数。
  • 删失机制：假设删失时间\( C_i \)与事件时间\( T_i \)在给定协变量\( X_i \)下独立（条件独立删失）。

  • 观测数据：\( \{(\tilde{T}_i, \Delta_i, X_i(t), t \le \tilde{T}_i), i=1,\dots,n\} \)独立同分布。

  • 可观测数据：研究者实际能观测到的是每个个体的随访时间\( \tilde{T}_i \)、事件是否发生\( \Delta_i \)，以及协变量在随访期间内的历史\( X_i(t) \)（对于时变协变量）。想要但观测不到的是真实事件时间\( T_i \)（被删失时）和基准风险函数\( \lambda_0(t) \)（非参数，需估计）。

第二步：讲最小内核

最简特例：假设协变量是时不变的（即\( X_i(t) = X_i \)），且无删失（\( \Delta_i = 1 \)对所有i）。此时，Cox模型的偏似然函数退化为一个简单的形式，但核心问题不变：如何用O(m)的计算量得到与全数据MPLE渐近等价的估计？

在这个特例下，本文的核心思路如下：

  1. 初始估计：从n个样本中随机抽取一个大小为m的子集（\( m \ll n \)），在该子集上计算MPLE \( \hat{\beta}^{(0)} \)。这一步的计算复杂度为O(m)。

  2. 一步校正：利用剩余n-m个样本，对\( \hat{\beta}^{(0)} \)进行一次Newton-Raphson更新。传统的一步更新公式为：

    \[\hat{\beta} = \hat{\beta}^{(0)} + \left[ \mathcal{I}(\hat{\beta}^{(0)}) \right]^{-1} \mathcal{S}(\hat{\beta}^{(0)})\]
    其中\( \mathcal{S}(\beta) \)是全数据的得分函数（sum over all n individuals），\( \mathcal{I}(\beta) \)是全数据的Fisher信息矩阵。计算\( \mathcal{S}(\hat{\beta}^{(0)}) \)和\( \mathcal{I}(\hat{\beta}^{(0)}) \)都需要遍历所有n个样本，复杂度为O(n)。

  3. 本文的关键想法：用估计的有效得分函数来替代全数据得分函数。具体地，作者构造了一个基于子集（大小为m）的“有效得分函数”\( \tilde{\mathcal{S}}(\beta) \)，使得：

  4. \( \tilde{\mathcal{S}}(\beta) \)是\( \mathcal{S}(\beta) \)的一个无偏估计（或近似无偏）。
  5. 计算\( \tilde{\mathcal{S}}(\beta) \)只需要O(m)的时间。
  6. 在\( \beta = \beta_0 \)（真实参数）处，\( \tilde{\mathcal{S}}(\beta_0) \)的方差与\( \mathcal{S}(\beta_0) \)的方差之比为m/n（即信息损失比例）。

然后，一步更新变为：

\[\hat{\beta} = \hat{\beta}^{(0)} + \left[ \tilde{\mathcal{I}}(\hat{\beta}^{(0)}) \right]^{-1} \tilde{\mathcal{S}}(\hat{\beta}^{(0)})\]
其中\( \tilde{\mathcal{I}} \)也是基于子集估计的Fisher信息矩阵。整个更新步骤的计算复杂度为O(m)。

  1. 为什么成立：在正则条件下，\( \hat{\beta}^{(0)} \)是\( \beta_0 \)的\( \sqrt{m} \)-相合估计。一步更新利用剩余n-m个样本的信息，将估计量的收敛速度从\( \sqrt{m} \)提升到\( \sqrt{n} \)，且渐近方差与全数据MPLE相同。关键在于，虽然\( \tilde{\mathcal{S}} \)只用了子集数据，但它通过巧妙的构造（利用剩余数据的某种“残差”信息）实现了与全数据得分函数相同的渐近效率。

一句话总结：本文的核心数学操作是用O(m)的计算量构造一个“有效得分函数”，使其在一步更新中起到与O(n)的全数据得分函数相同的作用。

三、这篇论文做了什么

三句话

  1. 研究了什么问题：在大规模生存数据（n达百万级）下，如何高效地拟合Cox比例风险模型，使得最终估计量的渐近分布与全数据MPLE相同，但计算时间仅为其一小部分。
  2. 核心工具/方法：子抽样+一步估计，其中一步更新使用估计的有效得分函数（estimated efficient score function），该函数基于子集数据构造，计算复杂度为O(m)而非O(n)。
  3. 主要结论：在正则条件下，本文提出的估计量\( \hat{\beta} \)与全数据MPLE \( \tilde{\beta} \)具有相同的渐近分布，即\( \sqrt{n}(\hat{\beta} - \beta_0) \xrightarrow{d} N(0, \Sigma) \)，其中\( \Sigma \)是Cox模型的半参数效率界。模拟和真实数据验证了其计算优势。

关键设定与假设

  • 设定：
  • 数据：\( \{(\tilde{T}_i, \Delta_i, X_i(\cdot)), i=1,\dots,n\} \)独立同分布。
  • 模型：Cox比例风险模型，\( \lambda(t|X_i) = \lambda_0(t) \exp(\beta^\top X_i(t)) \)。

  • 子集：随机抽取大小为m的子集（\( m = o(n) \)，但\( m \to \infty \) as \( n \to \infty \)）。

  • 假设（标准Cox模型正则条件）：

  • A1（有限区间）：随访时间在有限区间\( [0, \tau] \)内，且\( P(\tilde{T} \ge \tau) > 0 \)。
  • A2（有界协变量）：协变量\( X(t) \)在\( [0, \tau] \)上一致有界。
  • A3（正定信息矩阵）：全数据Fisher信息矩阵\( \mathcal{I}(\beta_0) \)正定。

  • A4（子集大小）：\( m \to \infty \)且\( m/n \to 0 \) as \( n \to \infty \)。（这意味着子集比例趋于0，但子集本身必须足够大以保证初始估计的\( \sqrt{m} \)-相合性。）

  • 相比已有文献的强化/放宽：

  • 相比Lin & Zeng (2006)：本文放宽了对一步更新计算量的要求（从O(n)到O(m)），但增加了对有效得分函数构造的假设（即需要估计基准风险函数\( \lambda_0(t) \)的某种“残差”）。
  • 相比Tong et al. (2014)：本文不需要数据分块和合并，避免了合并权重选择的问题。

主要结果

  • 定理1（渐近等价性）：在假设A1-A4下，本文估计量\( \hat{\beta} \)与全数据MPLE \( \tilde{\beta} \)满足：
    \[\sqrt{n}(\hat{\beta} - \tilde{\beta}) = o_p(1)\]
    即两者之差在\( \sqrt{n} \)尺度下可忽略。因此，\( \hat{\beta} \)与\( \tilde{\beta} \)具有相同的渐近分布。
  • 直觉：一步更新校正了初始估计的偏差，且校正量基于有效得分函数，其信息损失（由于只用了子集）在\( \sqrt{n} \)尺度下可忽略。

  • 必要条件：子集大小m必须满足\( m \to \infty \)且\( m/n \to 0 \)。这意味着m不能太小（否则初始估计不精确），也不能太大（否则计算优势消失）。实际中，作者建议m取几百到几千。

  • 定理2（方差估计）：给出了\( \hat{\beta} \)的渐近方差的一致估计量，该估计量也基于子集数据计算，复杂度为O(m)。这使得置信区间和假设检验可以在O(m)时间内完成。

  • 计算复杂度：本文方法的计算复杂度为O(mp² + mp³)，其中p是协变量维数。相比之下，全数据MPLE的复杂度为O(np² + np³)（每次迭代）。当n=10^6, m=10^3时，加速比约为1000倍。

证明路线与技术技巧

  • 整体路线：
  • 初始估计的相合性：证明子集MPLE \( \hat{\beta}^{(0)} \)是\( \beta_0 \)的\( \sqrt{m} \)-相合估计。这一步是标准结果。
  • 有效得分函数的构造：定义估计的有效得分函数\( \tilde{\mathcal{S}}(\beta) \)，并证明它是全数据得分函数\( \mathcal{S}(\beta) \)的一个“近似”，即：
    \[\tilde{\mathcal{S}}(\beta_0) = \mathcal{S}(\beta_0) + o_p(\sqrt{n})\]
    且\( \tilde{\mathcal{S}}(\beta) \)的导数（即有效信息矩阵）\( \tilde{\mathcal{I}}(\beta) \)满足：
    \[\tilde{\mathcal{I}}(\beta_0) = \mathcal{I}(\beta_0) + o_p(1)\]
  • 一步更新的泰勒展开：对\( \hat{\beta} \)进行泰勒展开：
    \[\hat{\beta} = \hat{\beta}^{(0)} + \tilde{\mathcal{I}}(\hat{\beta}^{(0)})^{-1} \tilde{\mathcal{S}}(\hat{\beta}^{(0)})\]
    利用\( \hat{\beta}^{(0)} \)的\( \sqrt{m} \)-相合性，将\( \tilde{\mathcal{S}}(\hat{\beta}^{(0)}) \)在\( \beta_0 \)处展开，得到：
    \[\tilde{\mathcal{S}}(\hat{\beta}^{(0)}) = \tilde{\mathcal{S}}(\beta_0) + \tilde{\mathcal{I}}(\beta_0)(\hat{\beta}^{(0)} - \beta_0) + o_p(\sqrt{n})\]

  • 代入并化简：将步骤3的结果代入步骤2的表达式，得到：

    \[\hat{\beta} - \beta_0 = \tilde{\mathcal{I}}(\beta_0)^{-1} \tilde{\mathcal{S}}(\beta_0) + o_p(1/\sqrt{n})\]
    由于\( \tilde{\mathcal{I}}(\beta_0) \approx \mathcal{I}(\beta_0) \)且\( \tilde{\mathcal{S}}(\beta_0) \approx \mathcal{S}(\beta_0) \)，最终得到：
    \[\hat{\beta} - \beta_0 = \mathcal{I}(\beta_0)^{-1} \mathcal{S}(\beta_0) + o_p(1/\sqrt{n})\]
    这正是全数据MPLE的一阶渐近展开，从而证明了渐近等价性。

  • 关键跳跃点：有效得分函数\( \tilde{\mathcal{S}}(\beta) \)的构造。作者需要设计一个基于子集数据的统计量，使其在\( \beta_0 \)处与全数据得分函数相差\( o_p(\sqrt{n}) \)。这要求\( \tilde{\mathcal{S}} \)能够“借用”剩余数据的信息，但又不需要遍历所有数据。作者的具体构造是：将每个个体的得分函数分解为“可预测部分”和“残差”，然后利用子集数据估计“可预测部分”的参数，而“残差”部分则通过剩余数据的某种加权和来近似。这一构造依赖于对基准风险函数\( \lambda_0(t) \)的Breslow估计器（一种非参数估计）的巧妙使用。

  • 技术技巧点名：

  • 经验过程理论（Empirical Process Theory）：用于控制\( \tilde{\mathcal{S}}(\beta) \)与\( \mathcal{S}(\beta) \)之间的差值，证明其一致相合性。
  • Breslow估计器：用于估计基准累积风险函数\( \Lambda_0(t) = \int_0^t \lambda_0(s) ds \)，这是构造有效得分函数的关键中间步骤。
  • 泰勒展开与delta方法：用于推导一步估计量的渐近分布。

真实例子与应用

  • 数据：UK Biobank数据，包含约50万参与者，随访时间中位数约10年。协变量包括年龄、性别、BMI、吸烟状况、血压等。结局变量为全因死亡率。
  • 方法应用：作者将本文方法应用于拟合Cox模型，子集大小m=5000（约为总样本的1%）。计算时间约为全数据MPLE的1/100（约2分钟 vs 200分钟）。
  • 结果：本文估计的回归系数与全数据MPLE几乎相同（差异在第三位小数以后），标准误也高度一致。这验证了理论结果：在有限样本下，本文方法确实达到了与全数据MPLE相同的效率。
  • 这个例子想说明什么：主要验证方法的实用性——在真实大规模数据上，本文方法在保持统计精度的同时，实现了显著的计算加速。

🔎 结论是否比证明窄

  • 窄结论：定理1的证明依赖于假设A4（\( m \to \infty \)且\( m/n \to 0 \)）。这意味着当m相对于n不是“足够小”时（例如m=0.1n），本文方法的计算优势会减弱，且渐近等价性可能不再成立。作者在模拟中测试了m=1000, 2000, 5000等值，但没有系统研究m/n比例对有限样本性能的影响。
  • 泛化claim：作者在结论部分声称该方法“可以推广到其他半参数模型”，但没有给出具体证明或例子。这是一个conjecture，而非严格证明的结果。

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

  1. 自适应子抽样：本文使用随机子抽样。能否设计一种自适应子抽样策略，使得子集包含信息量最大的个体（如事件时间极端或协变量异常的个体），从而在更小的m下达到相同效率？这扎根于本文对子集大小m的讨论（Section 2.1：“we randomly select a subset of size m”），作者没有探索非随机抽样。

  2. 时变协变量的高效处理：本文的方法在理论上适用于时变协变量，但有效得分函数的构造依赖于对每个个体在事件时间点的“风险集”的计算。当时变协变量频繁更新时，这一计算可能变得复杂。能否设计一种更高效的近似？这扎根于本文对有效得分函数构造的描述（Section 2.2：“the estimated efficient score function involves integrals over time”），作者没有讨论时变协变量带来的计算挑战。

  3. 高维协变量（p > m）：本文假设p固定且远小于m。当协变量维数p接近或超过子集大小m时，初始MPLE可能不唯一或发散，一步更新也会失效。如何将本文框架推广到高维Cox模型（如使用Lasso或SCAD进行变量选择）？这扎根于本文的假设A3（正定信息矩阵），该假设在高维下不成立。

  4. 与其他半参数模型的连接：作者声称方法可推广，但未给出具体路径。一个自然的问题是：本文的“有效得分函数”构造能否直接迁移到其他半参数模型（如加速失效时间模型、可加风险模型）？ 这扎根于本文结论部分的最后一句（“The proposed method can be extended to other semiparametric models”），这是一个未经验证的conjecture。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论