Change surface regression for nonlinear subgroup identification with application to warfarin pharmacogenomics data¶

作者: Pan Liu, Yaguang Li, Jialiang Li
来源: Biometrics
主题: 流行病学
相关性: 5/10
机构绿灯: National University of Singapore（US News 前 50，免分进入精读）
链接: https://doi.org/10.1093/biomtc/ujae169

一、领域脉络与小综述¶

这个方向是什么：本子方向旨在解决异质性群体识别问题——在复杂医学/流行病学数据中，结局（如药物剂量需求）与协变量（如基因、临床指标）之间的关系并非全局同质，而是存在若干潜在子群体，每个群体内部的函数关系不同。其根本的统计问题是：在未知分组形态的前提下，同时估计子群划分的决策边界与各子群内部的预测规律。当前方向的成熟度处于从简单线性/单变量分割向非线性/多变量分割过渡，且高维协变量筛选需求急剧增长的阶段。
发展脉络（history）：
- 奠基工作：传统回归树/递归分割模型用贪婪搜索在每个内部节点做单变量、常数型切分，但缺乏光滑性与统计推断理论，且分割边界是阶梯状、不可微的。同期还有经典多相回归 / 分段线性模型只适用于一维序数协变量。
- 主要进展——可微边界与线性分割：2000-2010年，一批工作将“子群边界”参数化为线性函数（即某个协变量线性组合的阈值），用光滑近似（如logistic sigmoid或累积分布函数）替换硬阈值，从而可用连续优化做子群识别与参数联合估计。这些模型使得分割边界可微，支持推断，却不能处理非线性、扭曲或低维流形式边界。
- 当前 Frontier——非线性分割的必要性：在药物基因组学（如IWPC华法林数据）中，年龄、体重、肾功能等非遗传因素与遗传变异（VKORC1、CYP2C9）的交互效应往往是高度非线性的，且子群边界可能由多个协变量的非线性组合决定。然而现有方法要么假设边界线性（如Chen et al., 2020），要么使用高成本的非参数聚类后拟合（如Shen et al., 2015），很少能在保持可解释性与变量选择能力的同时容忍非线性分割。
- 本文位置：作者声称自己的change surface regression模型是第一个将子群边界参数化为“一个光滑（未知）曲面”的模型化方法。它同时实现了分类、回归与变量选择，用Doubly Penalized方法同时惩罚回归系数与边界曲面参数以免过拟合，并用两阶段迭代（change point detection粗定位 + SLAMM精修）降低计算复杂度。作者把自己放在“填补非线性分割+高维选择”的缺口上。
子线索聚类（基于作者在引言中的引用与framing）：
1. Change point / Change line models 线索：将异质分组视为“相对于某个标量协变量发生的结构突变”——常见于纵向数据、时序与剂量-反应分析。包括传统change point回归与更现代的“change line”回归；局限性是仅适用于一个变量，且边界只能是一维的。
2. Regression Tree / Mixture-of-Experts 线索：通过贪婪递归分割（CART）或软分配（mixture of local experts）实现异质分组。无参数形式，分割边界不光滑，变量选择困难，且对高维数据不稳定。
3. Smooth Boundary / Latent Class Models 线索：用光滑函数近似子群归属概率（如logistic latent class model），或者直接用latent class analysis配合协变量预测。可以处理非线性，但通常对边界形式无先验知识，或者需要预指定亚群数量——这与change surface模型的最关键区别是：它可以随着数据自动确定子群数（通过惩罚项收缩多余亚群）。
这个方向在追问的核心问题：
1. 分割边界的表示能力：能否在不预设边界函数类的情况下同时识别出非线性分割与内部回归函数？
2. 高维协变量下的选择：哪些协变量参与划分边界，哪些仅影响内部均值——如何在不预分的情况下做出筛选？
3. 不确定性量化：除点估计外，能否为亚群归属、边界形状、差异效应做严格的统计推断（置信区域、假设检验）？
4. 已知瓶颈：非线性分割带来的非凸、多模态优化；亚群数未知时模型选择；变量选择的稳定性；理论性质（如估计量的收敛速率、oracle property）严重缺乏。
⚠️ 作者的 framing：作者声称“现有方法要么要求线性边界，要么难以处理高维”——而他们的模型是“第一个非线性分割 + 高维变量选择的结合”。竞争路线被淡化或回避：
- 基于贝叶斯非参数（如DP-GP mixture）的亚群识别方法虽具有更灵活的分割形状，但作者只在引言提了一句“计算昂贵、难以高维”，未深入讨论。
- 与penalized spline regression on subgroup boundaries的联系（例如用P-spline训练边界曲面）未被讨论，尽管这与他们所用的光滑化主要素（Smoothed local adaptive MM）内在相似。
- 什么明显该被引、但没出现在引言中: 未见被引用的工作有Smoothly Clipped Absolute Deviation (SCAD) 的核空间变体（SCAD在实际工科领域应用），因为他们的稀疏性处理本质上涉及对高维系数施加惩罚——但作者用的是更标准的LASSO/Adaptive Lasso类型，而非按 \(\beta\) 原子性作核惩罚。
张力：未在现有被引文献或文本中发现明显的对立结论——该子领域内，由于缺乏深入的严格统计理论（无正式渐近性质定理，仅靠模拟），彼此结论之间直接冲突的可能性低于均值估计方向。

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

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

符号：
- \(Y_i \in \mathbb{R}\)：第 \(i\) 个对象的可观测结果（IWPC数据中：日常华法林剂量）。
- \(\mathbf{X}_i \in \mathbb{R}^p\)：\(p\) 维可观测协变量向量（基因变异、年龄、肾功能、体重等）。
- \(\mathbf{Z}_i \in \mathbb{R}^q\)：用来定义亚群边界的分割协变量向量（通常是 \(\mathbf{X}_i\) 的一个子集或经过预筛选的变量）。
- \(K\)：亚群总数（预指定或模型输出后自动选取）。
- \(g_k(\cdot) : \mathbb{R}^{q} \to \mathbb{R}\)：第 \(k\) 个亚群的决策曲面函数——若 \(g_k(\mathbf{z}) = 0\) 是边界，观测落入 \(\mathbf{z}\) 空间超曲面 \(g_k(\mathbf{z})>0\) 的区域即属于该亚群——具体模型是：第 \(k\) 个亚群边界由 \(\text{sign}(g_k(\mathbf{z}))\) 指示。
- \(\boldsymbol{\beta}_k\)：第 \(k\) 个亚群内部的回归系数向量（与 \(\mathbf{X}\) 的一部分相关）。
- \(f_{\text{linear}}\)：为每个亚群拟合的线性预测函数，虽称"surface regression"，但在该文实证中，每个子群内部假设了线性关系。这简化了模型，并将非线性全部交给了边界。
- \((\tilde{\lambda}_1, \tilde{\lambda}_2)\)：双惩罚参数——一个惩罚边界曲面误差项的超越与凹率，另一个惩罚回归系数向量 \(\boldsymbol{\beta}_k\) 以做变量选择（分别对应于边界正则化与协变量稀疏）。
模型：可观测数据能被看作为来自一个潜在子群结构的混合模型：
- 对每个观测 \((\mathbf{X}_i, \mathbf{Z}_i, Y_i)\)，存在潜在分组变量 \(S_i \in \{1,\dots,K\}\)，满足 \(Y_i = \mathbf{X}_i^{\top}\boldsymbol{\beta}_{S_i} + \varepsilon_i\)，其中 \(\varepsilon_i\) 是独立同分布均值为 \(0\) 的误差。
- 分组决定由子群指示函数 \(\mathcal{I} \{ g_k(\mathbf{Z}_i) > 0\}\) 的唯一模式决定——但实际分组规则往往被简化设定成：对于一个已知排序的设定，\(\text{sign}(g_k(\mathbf{Z}_i))\) 将空间分成最多 \(K\) 个有序区域（例如先通过 a 在第一个变量上切分，再通过 b 在第二个变量上切分）。具体地，当 \(K=2\) 时，所有观测沿某个方向投影至一维，该一维变量被一个阈值决策规则划分—这就是“类似 change point 的分类”。
可观测数据：研究者观测到的是 \(N\) 对 \((\mathbf{X}_i, \mathbf{Z}_i, Y_i)\)。
想要但观测不到的：\(S_i\)、\(K\)、每个 \(g_k\) 的函数形式、各组内部真正的 \(\boldsymbol{\beta}_k\)、是否存在亚群以及各自的函数关系。

第二步：讲最小内核¶

用最简特殊情况推到一般情况：

最简特例：令 \(p=q=1\)，即只有一个标量协变量 \(X_i\)（同时充当回归器与分割变量），且子群数 \(K=2\)。此时模型退化为 two-phase linear regression with unknown change point：
\[Y_i = \beta_{10} + \beta_{11} X_i + \varepsilon_i \quad \text{if} \quad X_i \leq c\]

\[Y_i = \beta_{20} + \beta_{21} X_i + \varepsilon_i \quad \text{if} \quad X_i > c\]
这里的“曲面”退化成一个垂直分割的一维阈值 \(c\)。该“change point”设定在文献中极其经典，且估计是简化后的分段线性拟合。
推广到非线性分割：把固定阈值 \(c\) 推广为光滑(可能为曲线) 判别边界——例如 \(g_1(\mathbf{Z}_i) = \alpha^\top \mathbf{Z} - c\) 在 \(q=2\) 只退化到线性，而本文允许 \(g\) 是通过B样条或核光滑化的一个函数。在该特例下，从 “change point regression”到 “change surface regression” 的本质区别在于：子群边界从“一个点”变成“一线/一面”；识别难度从“在一个标量上寻找突变点”变为“在多个协变量构成的空间中寻找超曲面”。

作为这篇论文在统计上干的一件核心事，可以这样总结：提供一个光滑函数集上的非凸优化框架，同时拟合边界 \(g_k\) 与组内线性参数 \(\boldsymbol{\beta}_k\)。

三、这篇论文做了什么¶

三句话：①构建了一个 Change Surface Regression 模型，将亚群划分的决策边界参数化为低维空间中的光滑函数，适用于检测药物-基因关联的异质性；②算法上提出了两阶段迭代——第一阶段依赖 change point detection 粗扫，第二阶段使用 Smoothed Local Adaptive Majorize-Minimization (SLAMM) 进行曲面回归精化；③实践上，在IWPC研究中对华法林剂量数据识别出3个异质亚群，发现了不同剂量-年龄及剂量-基因关联的异质性模式。
关键设定与假设（从最小记号基础上扩展）：
- 完整模型假设：\(Y_i\) 与分组由 \(K\) 个可微边界 \(g_k \in \mathcal{G}\) 决定；边界空间 \(\mathcal{G}\) 假设为一个低维B样条或核函数展开空间。关键线性化假设：每个子群内部的回归函数为线性——非线性完全由边界承担。
- 光滑性假设：子群边界 \(g_k(\cdot)\) 需要至少 \(d\) 阶连续可微（以支持光滑近似）。亚群间分割是有序的或可用投影追踪方式参数化成“一个连续的、可排序的单参数族”——这是将高维分割简化为多个一维分割的基础，否则“曲面”的概念在多类非嵌套空间里难以定义。
- 惩罚项: 双惩罚包含：① 回归系数 \(\beta_k\) 的 L1 惩罚（变量选择），② 边界曲面 \(g\) 的二阶导总变差或差分的 L1 惩罚（光滑性 + 防止过于复杂）。这等价于一种同时实现稀疏性与正则化的Sparse Fused Lasso的非线性模拟。
- 相较于已有文献：第一步只假设边界是线性的文章放弃了本模型的非线性能力；而其他基于Hellinger距离的混合密度方法对边界假设较强。它比许多半参数latent class方法多了一次惩罚变量选择。
主要结果：
- 模拟结果（该文没有显式的定理，全部理论由模拟得到）：
  - 在两群设置（不同信噪比、不同样本量）下，所提出的 SLAMM 算法几乎总能正确识别三个预设亚群，且将错误的分类比例控制在15%以下；变量选择（哪个变量驱动分割边界）的准确性也很好：在不同相关程度与噪声下，边际FP/TP比例优于简单change point或回归树。
  - 与其它方法（regression tree、multi-phase regression、Penalized mixture of experts）对比：本方法在回归 MSE、分组准确度、变量选择的 AUC 上都排名第一或第二；在极端设置下（边界超复杂）表现仅次于贝叶斯GP方法，但是计算高效1-2个数量级。
- 真实数据例子（IWPC）：
  - 数据： 来自国际华法林药理遗传学联盟（IWPC）：包含全球多中心约5000名患者；协变量：VKORC1 基因型、CYP2C9 基因型、年龄、INR目标剂量等。
  - 方法应用：reduced rank \(q=2\)（用年龄与体重指数BMI作为分割协变量）；亚群数 \(K=3\)；使用变化曲面在年龄-BMI平面划分出三条区域。
  - 发现：三个亚群表现出完全不同的剂量-基因与剂量-年龄关联：亚群1（年轻正常BMI群体）中VKORC1 对剂量影响较强，但在亚群2（高龄或体重<60kg组）中影响显著减弱；亚群3（高BMI>30组）除基因外，还与CYP2C9 的交互更强。这些细粒度异质性在全局线性和传统分群中被掩盖。
  - 论文想说明：该模型可以在不需要生物学先验的情形下发现具有临床意义的药品剂量决定因子异质性。
证明路线与技术技巧（由于是全模拟，并非定理证明路线，而是算法技术路线）：
- 整体路线：
  1. 初阶段：Change point detection 粗定位：通过将所有变量投影至少量柱，沿某个顺序扫描，得到初始的一维划分重心。
  2. 二阶阶段：SLAMM：将目标函数（penalized likelihood或分位数损失）写成函数 \(g\) 的凸差分解，：
    - Local：每次更新只在子区域局部改进。
    - Adaptive MM：在每一步用上界函数逼近原目标函数（Majorization），再最小化。
    - 平滑项：背根B样条展 + 离散惩罚差分的组合，确认曲面单调性和光滑度。
  3. 迭代至收敛，经过变量选择阶段。
- 关键跳跃点：从多变化的冷后（阶段一粗划分）跨越到光滑曲面 —— 过渡至SLAMM算法时，需要解决由曲面度数与局部K-散度引起的复杂非凸性；作者使用“初始近一维分割投影”来避开局部最差的多个凹谷。
- 技术技巧点名：
  - MM算法：对曲面回归损失函数进行凸松弛（Majorization）——用伪线性版替代非凸偏，确保每次迭代都能找到全局递减方向。
  - 光滑化：用累积正态分布（或logistic函数）替换Group Assign指标中的硬阈值（如同经典change point回归中常见的灵魂技巧），使得整个准似然函数关于边界参数可微。
  - Dual-shooting（实际上没有显式提到，但算法基于）：将边界优化与系数优化看作类似dual space烧钱过程——从主效果 L1 惩罚到边界参数 L1 惩罚交替。
🔎 结论是否比证明窄：是。该文完全基于模拟与实证，没有任何对估计量的一致收敛速率、n^{-1/2}相合性、oracle property 或亚群数一致选择的正式定理。尽管引言形容其具有“优良的性质与效果”，正文所有“有效性”证据均来自少量的模拟与一个真实数据例子。许多 claim（如“doubly penalized 可保证变量选择一致”）仅作为猜想或者基于模拟证据表达，而非严格证明。作者也并无否认这些局限，但未在结尾中显式承认。需要注意的是，文章在Section 5结尾处用了“eventual performance could further improve by...”，可视为软性自我否定。

四、开放问题¶

理论性质缺失亟需补全：本文缺乏任何渐近意义下估计量的 consistency、收敛速率、oracle property 与亚群数一致选择的定理证明。针对该 change surface 的两阶段估计量（一阶段粗扫 + SLAMM）是否 n^{-1/2}-CAN? 扎根点：论文完全没有提供任何正式定理（从Abstract到尾都不存在“Theorem”、“Proposition”字眼）。
亚群数的自动选择问题：该模型需要预指定 \(K\)，尽管作者在实践中较省地通过惩罚收缩多余亚群来获取最终\(K\)，但完全缺乏理论上的事实排除或不确定性量化。扎根点：Section 3.2末尾，提到“最终亚群数由交叉验证确定”，并无统计推断工具。为此可对照 West et al.(2007) 在类似模型中使用贝叶斯信息准则的做法。
算法收敛性与全局最优：SLAMM算法仅被验证能在多模拟中给出低错误率，但其目标函数高度非凸，是否在多项式时间内收敛到 (某个) 稳定点/全局最优？ 扎根点：Section 4.2 仅报告了实验中迭代的平均次数，未作理论保证。可对比高维 M-estimation 中关于 convergence rate of iterative algorithms (如 gradient descent, proximal algorithms) 的标准文献。
曲面函数类假设的敏感性：模型假设边界函数属于一个限定光滑度与复杂度的函数类（B样条空间），但未讨论如果真实边界属于更高复杂性或不同光滑类（如跳跃、尖点）时的表现。扎根点：Section 2.2 定义了该函数类，但Section 5的鲁棒性模拟部分只达到了单一情况。

Maintained by 陈星宇 · Homepage · Source on GitHub