跳转至

Variable selection in functional linear Cox model

作者: Yuanzhen Yue, Stella Self, Yichao Wu, Jiajia Zhang, Rahul Ghosal
来源: Biometrics
主题: 非参数 / 半参数
相关性: 7/10
链接: 期刊页 · arXiv

一、领域脉络与小综述

这个方向是什么

本文解决的根本问题是：在高维功能数据与标量协变量并存的场景下，如何从关联时间-事件结局（如全因死亡率）的多个功能预测变量中进行变量选择。它处于“功能数据分析 (FDA)”与“生存分析”的交汇处，其成熟度适中——功能线性Cox模型已有估计和推断方法，但针对多个功能预测变量的自动变量选择方法此前尚付阙如。

发展脉络

奠基工作（2000s中期-2015）：

  • Tibshirani (1997) 和 Simon et al. (2011) 将Lasso和坐标下降引入Cox模型，建立了标量Cox变量选择的基线。
  • Gellar et al. (2015) 提出了单一功能协变量的功能线性Cox模型，结合惩罚信号回归和混合效应Cox模型，使用AIC/EPIC选择平滑参数。该工作为功能生存分析提供了首个可操作的估计框架。
  • Kong et al. (2018) FLCRM进一步发展了功能线性Cox模型的估计与推断，引入功能主成分分析（FPCA）进行降维，并给出了功能系数斜率函数的收敛速率和检验统计量。

主要进展（2015-2020）：变量选择方法从标量向功能扩展：

  • Zhang (2010) MCP惩罚的提出，解决了Lasso的偏差问题及其变量选择一致性所需的不必要表征条件（strong irrepresentable condition），提供了近乎无偏的高维线性回归变量选择。这是本文组MCP惩罚的理论基础。
  • Breheny & Huang (2015) 将MCP/SCAD推广到分组变量（组MCP、组SCAD），并开发了组下降算法，直接为本文提供了优化框架。
  • Gertheiss, Maity & Staicu (2013) 提出了广义功能线性模型（GFLM）中的变量选择方法，首次将组惩罚应用于功能系数，同时实现了功能系数的稀疏性和平滑性。该方法与本文最为接近，但仅限于广义线性模型（如标量响应），未涉及生存数据。
  • Ghosal et al. (2020) 将组变量选择扩展到功能线性共时回归（functional linear concurrent regression），证明了组MCP/SCAD在高噪声、稀疏观测下依然能准确选择变量。

当前frontier（2020-至今）：处理更高维、更复杂的分布特征与非线性效应：

  • Cui et al. (2021) 提出了可加功能Cox模型，允许功能协变量与log风险之间的非线性关系，进一步提升了灵活性，但其估计和推断仍限于单一功能协变量，且不专注于变量选择。
  • Ghosal et al. (2021, 2023)、Cho et al. (2024) 将功能数据方法应用于身体活动（PA）的日分布特征——使用时变L-moments（分位数、变异性等）而不是简单的均值曲线作为功能预测变量。这一演变提出了一个新的变量选择需求：多个功能协变量（代表不同的分布特征）同时存在时，哪一个/些是真正与结局相关的？
  • Feng et al. (2023) 进一步指出PA的时间窗口（如上午 vs. 下午）也具有不同效应，强化了“滑动窗口-功能特征”方法学的必要性。

本文所处的位置：本文是第一个将同时处理多个功能协变量和标量协变量的组MCP变量选择方法引入功能线性Cox模型的。它填补了Gertheiss et al. (2013)（功能线性模型，非生存）和Gellar et al. (2015)/Kong et al. (2018)（功能Cox模型，无变量选择）之间的空白。

子线索聚类

  1. 标量生存分析的变量选择（Tibshirani 1997, Simon et al. 2011, Du et al. 2010, Wu 2012）：将Lasso、弹性网、SCAD、自适应Lasso应用于Cox模型。本文的优化算法和参数选择方法部分借鉴了这些思想，但将目标从标量系数扩展到功能系数。

  2. 功能线性模型，变量选择不（Gertheiss et al. 2013, Ghosal et al. 2020）：这一簇回答了“在非生存结局的回归中，如何从多个功能预测变量中选出重要变量”，并产生了组MCP/SCAD这类成熟工具。本文的直接前驱就在这里——它将这一簇的方法逻辑从标量响应（如连续或二值Y）迁移到时间-事件结局（右删失）。

  3. 功能Cox模型估计（Gellar et al. 2015, Kong et al. 2018, Cui et al. 2021）：提供了功能Cox模型中的估计与推断工具但缺少统一的变量选择机制。本文将其与线索2结合，推导出针对Cox模型的部分似然的惩罚函数和组下降算法。

  4. NHANES身体活动的分布特征分析与流行病学（Leroux et al. 2019, Saint-Maurice et al. 2020, Ghosal et al. 2021, Cho et al. 2024）：这一簇指明了应用场景和科学问题——从步数、MVPA等标量摘要 \(\to\) 日间PA曲线分布 \(\to\) 时变L-moments。本文的应用层就站在这一簇的on top。

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

  1. 如何同时处理多个功能协变量与多个标量协变量的变量选择？ 传统的“一次检查一个功能变量”（Gellar et al. 2015; Cui et al. 2021）无法回答“哪一个分布特征（如均值 vs. 变异性）比另一个更重要”这种多重比较问题。

  2. 如何在部分似然框架下实现组级稀疏？ 每个功能协变量的系数曲线用多基展开，其系数自然构成一组。需确保组内全为0或非0。

  3. 如何自动选择平滑参数与稀疏参数，避免用户主观调参？ 两个正则化项（二阶差分平滑 + 组MCP稀疏）的平衡需自动化。

已知瓶颈：现有方法要么无变量选择能力（Gellar, Kong），要么只能处理标量响应（Gertheiss），要么只处理单一功能变量（Cui）。缺乏一个统一的、支持自动参数选择的、可扩展至约几十个功能预测变量的框架。

⚠️ 作者的 framing

这是作者的说法（必须明确）：作者将缺口 frame 为“现有的Cox变量选择方法专注于标量预测变量；现有的多功能预测方法缺乏变量选择”，因此本文是“显然的下一步”，将Gertheiss et al. (2013)中的组MCP惩罚首次引入功能比例风险设定。他们淡化了两点：(1) 已有同时处理多个功能变量且进行变量选择的方法吗？——没有直接的，但Gertheiss是指标GFLM的，他们只是扩展DOMAIN；(2) 算法的理论保证（Oracle性质、收敛速率）是未证明的，只在模拟中展示经验表现。

什么明显该被引/该存在、却没出现在intro里？ 本文是被引论文 [20] FLCRM (Kong et al., 2018) 提供了类似的估计方法，但用了FPCA而非样条。作者只在引言中一句带过，并未深入讨论与FPCA基方法的比较（如计算复杂性、对稀疏/不规则观测的鲁棒性差异）。这可能是值得进一步探索的议题。

张力

未见明显对立引用。各被引工作之间在方法论上基本兼容渐进：标量变量选择 → 功能线性模型变量选择 → 功能Cox模型估计 → 功能Cox模型变量选择。核心参考文献一致地指向组MCP、功能线性模型框架与Cox模型。

二、最小内核 / 最简单的数学问题

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

符号名册：

记号 含义
\(i = 1,\dots,n\) 独立同分布（iid）个体的索引
\(T_i\) 真实事件时间（潜=潜在变量，非完全可观测）
\(C_i\) 删失时间（注目 = unobservable 若 \(T_i < C_i\)）
\(\Delta_i = I(T_i \le C_i)\) 删失指示（可观测；1=事件发生，0=删失）
\(Y_i = \min(T_i, C_i)\) 可观测的随访时间（到事件/删失）
\(Z_i = (Z_{i1},\dots, Z_{ip})^\top \in \mathbb{R}^p\) 标量协变量向量（可观测）
\(\{ X_{ik}(s): s \in [0,1] \}, k = 1,\dots,K\) 第\(k\)个功能协变量的曲线，定义在标准化域 \(s\in[0,1]\)（如24小时）
\(\beta = (\beta_1,\dots,\beta_p)^\top \in \mathbb{R}^p\) 标量系数（要估的对象）
\(\beta_k(s), s\in[0,1], k=1,\dots,K\) 功能系数/斜率函数（要估的对象）
\(h_i(t)\) 个体\(i\)在时间\(t\)的瞬时风险函数
\(h_0(t)\) 基准风险函数（非参、未指定）

模型：比例风险假设：

\[h_i(t) = h_0(t) \exp \left( \beta^\top Z_i + \sum_{k=1}^K \int_0^1 X_{ik}(s) \beta_k(s) \, ds \right)\]
即，标量协变量的效应是 \(\beta^\top Z_i\)（线性，参数）；功能协变量的效应是每个功能曲线与对应系数函数的内积之和（线性、但用于参数基展开再估计）。\(h_0(t)\)完全非参，不出现在估计中（部分似然可以消去）。

可观测数据与想估的对象： - 可观测： - \(\{(Y_i, \Delta_i, Z_i, X_{i1}(s),\dots,X_{iK}(s))\}_{i=1}^n\). - 其中，功能曲线 \(X_{ik}(s)\) 通常在精细栅格上离散观测（如每隔1分钟）。 - 想估但不可直接观测的： - \(\beta\), \(\beta_k(s)\). - 部分似然只用到这些参数；\(h_0(t)\) 是一个讨厌参数（nuisance），不直接用于变量选择。 - 关键区分：“变量选择”不是选择 \(X_{ik}(s)\) 的“子区域”，而是选择整个 \(k\) 是否与结局相关（即 \(\beta_k(s)=0\) 对几乎所有 \(s\)）。因此这是一个组选择问题。

第二步：最小内核

特例：\(K = 1\) 个功能协变量，\(p = 0\) 个标量协变量。

情况退化至单功能线性Cox模型，变量选择问题消失了（只有一个变量）。为了展示整篇论文的核心概念运作，我们需要扩展特例：

最小问题（去掉多余一般性假设的“彩蛋”版本）：设 \(K=2\)，即只有两个功能协变量 \(X_{i1}(s), X_{i2}(s)\)。真实的数据生成模型中，仅 \(X_{i1}(s)\) 对风险有贡献（即 \(\beta_1^*(s) \neq 0\) 对非零测度的 \(s\)），而 \(\beta_2^*(s)=0\) 对所有 \(s\)。我们的目标是：在 \(n\) 不太小（如 \(n=200\)）、功能曲线在 \(M=1000\) 个时间点观测的设定下，设计一个统计程序使——大概率——选 \(X_{i1}\) 为重要、\(X_{i2}\) 为不重要。

记号简化： - 每个功能协变量用 \(L\) 个B样条基函数 \(\{B_1(s),\dots,B_L(s)\}\) 近似展开。记 \(X_{ik}(s) \approx \sum_{l=1}^L x_{ikl} B_l(s)\)，其中 \(x_{ikl}\) 是观测值在基上的系数（实际中是投影）。 - 对第 \(k\) 个功能系数：\(\beta_k(s) \approx \sum_{l=1}^L b_{kl} B_l(s)\). - 于是，\(\int X_{ik}(s) \beta_k(s) ds \approx \mathbf{b}_k^\top \left( \int \mathbf{B}(s) \mathbf{B}(s)^\top ds \right) \mathbf{x}_k\)，可以吸收进设计矩阵的效率：记 \(U_{ik} = \mathbf{b}_k^\top \tilde{\mathbf{x}}_{ik}\)（其中 \(\tilde{\mathbf{x}}_{ik}\) 是基内积后的版本）。 - 组 \(\mathbf{b}_k = (b_{k1},\dots,b_{kL})^\top\) 构成第 \(k\) 个功能预测变量的“组”向量。在 \(X_{i2}\) 无关时，应有 \(\mathbf{b}_2 \approx 0\)（组级稀疏）。

在这个特例下，目标退化成什么？ 目标变为：从两部分似然出发，最小化：

\[-\ell_n(\beta_1(\cdot),\beta_2(\cdot)) + \lambda_1 \sum_{k=1}^2 \mathcal{P}_{\text{smooth}}(\beta_k) + \lambda_2 \sum_{k=1}^2 \rho( \| \mathbf{b}_k \| ; \text{group MCP})\]
其中： - \(-\ell_n\) 是Cox部分对数似然的负值（仅使用观测的 \((Y_i,\Delta_i)\)）。 - \(\mathcal{P}_{\text{smooth}}\) 是二阶差分惩罚（使 \(\beta_k(s)\) 平滑）。 - \(\rho(x; \gamma) = \gamma \int_0^x (1 - t/(\gamma \lambda_2))_+ dt\) 是MCP惩罚在组上的形式。当组范数 \(\| \mathbf{b}_k \|\) 小于 \(\gamma \lambda_2\) 时，惩罚项线性降低；当超过 \(\gamma \lambda_2\) 时，惩罚项变为常数，不再收缩该组（保留大系数）。

为什么成立？ MCP有一个特性：它对大系数几乎没有惩罚（退化为常数），从而避免了Lasso对大系数的系统性偏差。在非凸性允许的范围内，组MCP会驱动 \(\mathbf{b}_2\) 向零收缩；一旦在一轮迭代中 \(\mathbf{b}_2\) 的范数跨过阈值 \(\gamma \lambda_2\)，它便不再被惩罚（从而数值上不会消失），但若实际效应为零，它会稳定于0。Breheny & Huang (2015) 的组下降算法允许高效求解该非凸目标的局部最小值（可能在主分支）。

因此，本项工作的核心数学任务是：定义一个检视组（每个功能系数向量 \(\mathbf{b}_k\) 为一组），构造一个惩罚，在Cox部分似然下迫使无关组的全部L个系数同时为0。这是Gertheiss (2013)（针对标量响应GFLM）在Cox生存设定下的直接对应。

三、这篇论文做了什么

三句话

  1. 研究了什么问题：在功能线性Cox模型下，同时存在多个功能协变量和标量协变量时，进行变量选择（即识别哪些功能预测变量的系数曲线 \(\beta_k(s)\) 恒为 0）。这是一个组稀疏的惩罚部分似然估计问题。

  2. 核心工具/方法：用B样条基展开近似功能系数，然后在部分似然中添加两个惩罚项：(i) 二阶差分惩罚 (smoothness)；(ii) 组MCP惩罚 (group selection)——从而实现平滑性与组稀疏性的联合。使用Breheny & Huang (2015)的组下降算法优化，并通过BIC/GCV自动选择 \(\lambda_1\) (平滑) 和 \(\lambda_2\) (稀疏) 两个调参，外加MCP形参数 \(\gamma\)。

  3. 主要结论：模拟表明方法在不同信噪比、不同组规模下具有准确的变量选择能力（真阳性率 >0.9 在多数配置下，假阳性率 <0.1）和较好的系数估计精度（MISE低）。NHANES应用识别出多个日间PA分布特征（如时变均值、变异性L-moments）以及年龄、性别、BMI、吸烟等标量协变量为全因死亡率的显著预测因子。

关键设定与假设

  • 记号（已在第二节交代）：增加向量、矩阵记号—— \(\mathbf{X}_{ik} \in \mathbb{R}^{M_k}\) 是功能变量在栅格观测点上的离散化值；\(\mathbf{B}_k \in \mathbb{R}^{M_k \times L_k}\) 是基矩阵；\(\mathbf{\Phi}_k = \mathbf{B}_k^\top \mathbf{B}_k \approx \int \mathbf{B}_k(s) \mathbf{B}_k^\top(s) ds\)。
  • 假设：
  • (A1) 比例风险：\(h_i(t) = h_0(t) \exp(\eta_i)\)，其中 \(\eta_i = \beta^\top Z_i + \sum_{k=1}^K \int X_{ik}(s) \beta_k(s) ds\).
  • (A2) 删失独立于事件时间给定 \(\mathbf{Z},\mathbf{X}\)（条件独立删失） ：\(T_i \perp C_i \mid (Z_i, X_i(\cdot))\).
  • (A3) 功能协变量的支撑在规范化域 \([0,1]\)；功能系数 \(\beta_k(s)\) 光滑（\(s\) 的二阶导数存在且有限）。
  • (A4) 基函数用 \(L_k\) 个等间距CubicB样条（统一长度）；这一选择使得L用在中等范围（模拟中 \(L=12\) 左右）。没有可去基/FPCA的自适应选择。
  • 有哪些设定是相比Gertheiss (2013)更紧/更弱？
  • Gertheiss针对广义线性模型（如正常、二值响应），这里是部分似然（生存）。这使得梯度/海森矩阵计算复杂得多（二阶导涉及时间序统计量）。
  • Gertheiss使用类似自适应LASSO的权重——本文在核心方法中同样使用了类似的自适应权重 \(w_k = 1/\|\hat{\beta}^{(0)}_k\|\)（基于未惩罚估计），但这是针对组之间的变量进行二次加权（addressing 组间尺度差异）。
  • 本文没有提供Oracle性质的理论证明（只提及以Gertheiss(2013)和Zhao(2010)的结果为参考）。

主要结果

  • 模拟评估变量选择：在多种配置（K=5或10个功能变量，只有1-2个真正有效，不同SNR和样本量 n=300, 500, 1000）下，本文方法（VSFLCOX——Variable Selection in Functional Linear Cox model）：
  • 真阳性率 (TPR): 在最高SNR ~0.97-1.0，中等SNR ~0.9-0.96。
  • 假阳性率 (FPR): 保持 <0.1 甚至 <0.05 (n=1000) 在多数配置下。
  • 系数估计的MISE：显著低于未惩罚估计（因为无效组的不必要噪声被剔除）。
  • 比较对象：组Lasso、组SCAD。组MCP表现出优于组Lasso的高TPR和低于组SCAD的FPR（在中等SNR下）。
  • 模拟评估参数选择：BIC调参的自动选择在与人工选择（已知真实模型）的对比中显示出很高的稳健性（Corr>0.9）。
  • 与Gertheiss w/ Adaptive weighting 差异：只在N=100-300高噪声低SNR时，自适应权重才小幅改善性能；在常规设定下，自适应的增加带来的好处有限。

证明路线与技术技巧（本文为应用型论文，但仍有可挖掘的技术细节）

  • 整体路线（估计与选择）：
  • 基与惩罚结构：每个 \(k\) 的系数曲线 \(\beta_k(s)\) 用 \(L\) 个B样条基展开。将\(\ell_2\)范数惩罚用于组收缩，二阶差分惩罚用于控制平滑。
  • 部分似然：风险集序贯结构用于求偏似然的梯度和海森矩阵；使用坐标下降：对每个组\(k\)，固定其他\(j\neq k\)的系数，用线搜索解决单组（L变量子问题）；
  • 组下降算法（Breheny & Huang 2015）的核心：对单组在部分似然下解一个L维二次近似的MCP惩罚回归，应用 必要条件 (univariate thresholding)与次梯度条件。由于MCP的非凸性，无法保证全局最优解，但可以求到“主分支”（main branch in critical point graph: 从原点出发，随着\(\lambda\)递减，路径连续）。
  • 参数选择：对 \(\lambda_2\)（稀疏）使用BIC（或GCV）；对 \(\lambda_1\)（平滑）和外型参数 \(\gamma\)（MCP形参）使用网格搜索+外循环BIC。BIC用模型的非零组数目+所有\(\beta\)参数的B样条系数总数目作为参数数目。
  • 关键跳跃点：在使用组下降时，Cox部分似然的Hessian矩阵因删失的时变风险集而庞大（与事件数目\(O(n^2)\)成正比）。算法只利用对角线近似（等价于对有删失数据的独立Cox风险集做近似），这对非凸MCP是风险点。作者没有从理论上分析这种近似是否会偏向于局部极小较差的解，但模拟会经验性地验证其有效性。
  • 技术技巧点名：
  • B样条基与标准化正交化（Gram矩阵正交化 \(\mathbf{\Phi}_k^{1/2} \mathbf{b}_k\)）来去相关，使可以独立地进行单变量阈值化。这在Breheney (2015)组下降算法中很关键。
  • MCP group penalty：非凸但坏点少，支持连续路径搜索。
  • Cox部分似然的log / exp量 用事件时间排序和链式前向/后向求和加速（Cox蒙特卡洛技巧如风险集的‘提前计算’）。
  • 自适应加权：以\(\frac{1}{\|\mathbf{b}_k^{(0)}\|}\)乘以组惩罚，对较大系数组惩罚小，对小的惩罚大（类似adaptive group lasso）。

真实例子与应用

数据集： - NHANES 2003-2006：年龄>=50岁的3618名参与者。 - 功能预测变量：从加速度计数据（7天，以分钟为单位）提炼出11个“日间分布 L-moments”功能曲线： 1. PA均值 curve (分段)$\tilde= $ 即 \(X_1(s)\) - 当天时间\(s\)的原始平均活动水平。 2-5. 其他L-moments（变异性、不对称性、峰度等 with time-of-day dependency）.这些是“分布数据对象”Ghosal (2021)的利用。 - 标量协变量：年龄、性别、种族、教育、BMI、吸烟状况、晨尿总钙/C反应蛋白等。 - 结局：全因死亡率（NDI至2011年）。

怎么把方法用上去： - 将11个L-moment曲线作为功能协变量（K=11）；按目前文献所建议的自动参数选择流程。选择后得到 ≈4 个选择的显著功能变量 (PA mean X1, 变异 X2, 早上6-10点的特定曲线系数) 以及一众标量变量（Age, Smoking 等)。 - 为进一步证实选择有效性，通过引入10个伪功能协变量（噪声序列），本文方法成功地将它们的组系数归于0（FPR有控制）。

这个例子想说明什么： 1. 实际相关性：日间PA分布特征（不止均值）与全因死亡率显著相关，且这种模式是调整了人口学/生活方式协变量后的独立效应。 2. 方法论好处：在多功能变量的生存分析中，自动变量选择可以避免人为功夫（将许多无用的L-moment曲线手动排除）而只保留重要的分布特征。特别是PA变异性（mid-morning variability）首次被识别为死亡风险的关键因子（而传统均值曲线方法可能忽略这一信号）。 3. 稳健性：通过插入冗余变量未被选中，验证了变量选择的一致性。

🔎 结论是否比证明窄

是，有观察到的“证明 - 语句与结论范围”的缺口： - 摘要说“efficient variable selection method”，但全文没有证明方法的Oracle性质（variable selection consistency / estimation consistency）。对比Zhang (2010)对高维线性回归MCP的严谨Oracle证明，本文仅提供了模拟。 - Introduction中说“We propose a novel variable selection method”；但核心部分（样条+组MCP+组下降）在Breheny & Huang (2015, group descent)和Gertheiss et al. (2013, 功能变量选择)中完整出现过。真正的novelty主要是将这两者插到Cox部分似然中，并用PA L-moment数据做实证。可接受的scope是“方法应用”，不是纯新理论。 - 文中没有指出高维（p或K很大时理论选择一致性的条件假设）。目前的模拟最多只到K=10（10个功能变量），但结论没有讨论延伸到更高维时的计算收敛性或统计一致性极限。

四、开放问题（扎根具体语句）

  1. 后选择推断：“…explore developing uncertainty quantification methods for the estimated functional effects in VSFLCox using post-selection inference approaches (Lee et al., 2016; Taylor and Tibshirani, 2018) and sample-splitting strategies (Wasserman and Roeder, 2009).” ——这是一个直白 gap：本文只做了点估计和变量选择，没有给出选择后的置信区间或假设检验（如Lee 2013的post-selection inference对于Lasso的功能扩展），尤其是针对功能数据特有的组性状函数的推断。

  2. 非线性效应：“The proposed variable selection method could be extended to accommodate additive models (Cui et al., 2021) which would offer greater flexibility by allowing nonlinear additive effects of the functional and scalar predictors.” ——目前的函数效应是线性的（内积形式）。推广到非线性的可加结构会导致基组极大膨胀（\(L × T\) 倍），会使组稀疏惩罚失效。要如何针对这个非线性场景设计新的稀疏结构是有待回答的。

  3. 理论Oracle性质：全文未给出任何估计一致性公式。作者仅在模拟部分展示经验表现，但在引言或结论中完全没有与Oracle性质相关的引理。这是一个明确的信仰缺口——能否为组MCP在Cox部分似然设定下推导变量选择一致性的渐近条件（类似Zhang 2010中强烈不必要条件的缺失）？是否这个设定天然地由于基准风险的存在而更困难？

Maintained by 陈星宇 · Homepage · Source on GitHub

评论