Mixture survival trees for cancer risk classification

Beilin Jia; Donglin Zeng; Jason J. Z. Liao; Guanghan F. Liu; Xianming Tan; Guoqing Diao; Joseph G. Ibrahim

doi:10.1007/s10985-022-09552-w

寿命数据分析。作者手稿；PMC 2023 8月4日提供。

以最终编辑形式发布为：

寿命数据分析。2022年7月；28(3): 356–379.

2022年4月29日在线发布。数字对象标识：10.1007/s10985-022-09552-w

预防性维修识别码：项目经理10402927

美国国立卫生研究院：美国国立卫生研究院1918361

PMID：35486260

用于癌症风险分类的混合生存树

贾碑林,¹ 曾东林,¹ Jason J.Z.廖,² 广汉F.Liu,^三谭显明,¹ 郭庆雕,⁴和约瑟夫·易卜拉欣¹

作者信息版权和许可信息 PMC免责声明

摘要

在肿瘤学研究中，重要的是了解和表征患者之间的疾病异质性，以便将患者分为不同的风险组，并在正确的时间识别高危患者。然后，可以使用这些信息来确定更同质的患者群体，以开发精确医学。本文提出了一种直接风险分类的混合生存树方法。我们假设患者可以被分为预先指定的风险组，其中每个组都有不同的生存情况。我们提出的基于树的方法设计用于使用EM算法估计潜在的组成员。在递归划分中，使用观测数据的对数似然函数作为分割准则。通过广泛的模拟研究评估了有限样本的性能，并通过乳腺癌的案例研究说明了所提出的方法。

关键词：审查、潜在模型、混合分布、风险分类、基于树的方法

1 介绍

在肿瘤学研究中，由于癌症患者之间存在广泛的异质性，以及进行肿瘤研究的巨大成本，确定癌症复发或死亡的高危或易感人群对于药物开发至关重要。准确的风险分类，即根据患者的临床状态、遗传标记和环境暴露评估患者的生存率，对于在精确医学时代开发靶向癌症治疗也至关重要。

参数或半参数生存模型通常用于评估个体风险。例如，可以拟合Cox比例风险模型(考克斯1972)或其他转换模型(Zeng等人2016)获得作为个体协变量线性函数的风险评分。这种风险评分的一个例子来自弗雷明汉心脏研究。分析(Kannel和McGee 1979年;Kannel等人，1979年;Wilson等人，1998年)研究表明，冠心病风险与年龄、糖尿病、血压、胆固醇水平和吸烟状况有关。另一个例子与晚期口腔癌有关，其中Tseng等人（2020）通过使用Cox比例风险模型和预先指定的风险分层阈值，根据综合临床病理学和遗传数据将患者分为高风险组和低风险组。

还开发了非参数和机器学习方法来评估患者的风险。决策树由于其简单性和可解释性，已成为文献中处理时间到事件数据的流行方法。具体来说，存活树(Ciampi等人，1981年;Marubini等人，1983年;Gordon和Olshen 1985年)为了扩展现有的基于树的连续结果方法来处理生存数据，构建了分裂规则来优化节点内同质性和节点间异质性。最近，人们对各种分裂标准进行了大量讨论，例如似然比检验(Ciampi等人，1987年)，指数对数似然损失(Davis and Anderson 1989年)，完全似然偏差(LeBlanc和Crowley 1992年)，两个子节点生存函数之间的积分绝对差(Moradian等人，2017年)以及评估两个子节点危险性差异的综合协调措施(Sun等人2019年). 将患者分为不同的风险组、生存功能(Ibrahim等人，2009年;周和麦卡德尔2015)或危险功能(Vergara等人，2018年)在每个终末淋巴结处获得并与预先指定的阈值进行比较，以确定每个患者的高风险或低风险组。

半参数和机器学习方法在生存风险分类方面都有一些局限性。首先，参数或半参数模型旨在研究协变量和结果之间的关联，因此实际上并没有开发用于风险分类。例如，由于考克斯模型风险评分的单调性而导致的模型错误指定可能导致风险组的严重错误分类。其次，尽管生存树等机器学习方法对模型错误指定更为稳健，但它们是为生存预测而设计的，而不是直接用于风险分类。此外，递归划分期间的节点分裂通常基于比较嵌套节点的生存函数，以便在终端节点产生许多不同的生存函数。因此，风险分类的决策选择很困难。相反，这些方法必须依赖于生存功能的一些粗略汇总统计数据，例如给定时间点的中位生存率或生存概率，以对患者进行分类，而这些患者很可能会错过个体生存状况的整体情况。最后，所有这些方法都依赖于选择阈值进行分类，这可能是主观的，可能没有临床意义。

相反，一种更直接的风险分类方法是将患者的风险组标签视为缺失数据，以便可以使用观察到的数据推断这些标签。例如，廖和刘（2019）指出，对于某些特定的癌症，如黑色素瘤，患者群体可以近似分解为两个或三个潜在群体，每个潜在亚群具有独特的生存特征。因此，可以考虑将多个生存分布与未知参数混合，以描述每一潜在患者组的生存情况。由此得出的群体成员估计值是协变量的参数函数，直接为每个个体提供风险分类。有限混合模型被广泛用于研究生存数据的异质性(Larson和Dinse 1985;再见1982).沈和（2015）构建了一个结构化的logistic-normal混合模型来识别治疗效果增强的亚组。在建模之前，他们进一步进行了验证性统计测试，以检查亚组的存在。一种准牛顿EM算法(Bussy等人2019年)基于离散生存数据探讨患者风险组。然而，所有这些方法都是参数化的，假设协变量和风险类别之间存在限制性关系，因此存在模型指定错误。当涉及大量协变量时，计算也具有挑战性。

在这项工作中，我们提出了一种用于风险分类的混合生存树方法。具体来说，我们假设患者可以被划分为预先指定的风险组，其中每个组都有不同的生存曲线，这些曲线是使用威布尔分布的一般族建模的。我们根据患者的协变量对各组成员进行非参数建模。对于估计，我们在EM算法中采用基于树的方法来估计组成员。在每次迭代中，观察到的对数似然函数被用作分裂标准，以优化患者生存行为方面的节点内同质性和节点间异质性。由于在迭代中只使用二进制分裂，因此计算速度很快，并且我们表明似然函数随着迭代而增加。更重要的是，我们提出的方法为未来患者提供了直接的风险分类，因为可以使用估计树显式地导出风险组成员。本文的其余部分组织如下。第2节阐述了提出的基于树的生存风险组发现方法的动机。第3节提供了有关提议的基于树的方法的详细信息。第4节通过大量的仿真研究证明了该方法的性能。实际数据应用程序如所示教派。5.

2 IBCSG乳腺癌试验

国际乳腺癌研究小组（IBCSG）进行的一项大型临床试验IBCSG试验VI推动了拟议混合生存树的开发。IBCSG试验VI研究了绝经前淋巴结阳性乳腺癌患者辅助化疗的持续时间和延迟化疗的再次引入(Colleoni等人，2002年;Chi和Ibrahim 2006). 患者按照二乘二析因设计随机接受以下治疗：（a）环磷酰胺、甲氨蝶呤和氟尿嘧啶（CMF）连续六个周期（CMF6）；（B） CMF6加上三个单循环再引入CMF（CMF6+3循环）；（C）连续三个循环的CMF（CMF3）；和（D）CMF3加上三个单独的再引入CMF循环（CMF3+3循环）。患者的生活质量（QOL）也在基线时进行测量，并假设其包含预后信息和反映乳腺癌进展。通过自我评估生活质量问卷评估生活质量的四个方面，包括身体健康、情绪、食欲和感知应对，并在入组后每三个月评估一次患者的生活质量。除了治疗效果和患者的生活质量外，还考虑了无病生存期（分析中重新标度为[0,1]）、事件状态、基线年龄、雌激素受体（ER）状态（1=阳性，0=阴性）和肿瘤阳性淋巴结数（即淋巴结组，1=5个或更多阳性淋巴结，0=其他）。在原始数据中，肿瘤阳性淋巴结的数量被二分法(Chi和Ibrahim 2006). 962名患者在基线检查时有完整的观察结果。无病生存期的中位随访时间为7.47年，事件发生率约为45%。有关试验的更多详细信息，请参见国际乳腺癌（1996）.

为了探讨患者的异质性，Kaplan-Meier曲线(卡普兰和迈耶1958)根据40岁以下患者和40岁或40岁以上患者的无病生存率，用以证明患者生存情况的差异。40岁年龄的截止值是从Colleoni等人（2002年）年轻患者可能比老年患者的生存经验更差，因为年龄与闭经的可能性有关。40岁以下患者和40岁以上患者的癌症进展可能不同。图1显示40岁以下的患者对不同治疗的反应不同。例如，对于年龄小于40岁的患者，接受治疗B（CMF6+3周期）和C（CMF3）的患者的生存曲线在6年后趋于平稳，但接受治疗A（CMF6）和D（CMF3+3周期）的患者没有观察到这种现象。此外，在治疗组中，40岁以下患者和40岁以上患者的生存曲线形状不同。对于治疗C（CMF3）和D（CMF3+3周期），40岁以下患者的生存风险高于40岁或40岁以上患者。这些发现表明，患者在治疗方面存在潜在的异质性。由于数据集中还包含其他基线患者特征，为了探讨基线协变量如何决定患者的潜在异质性，我们提出了一种基于树的方法来发现具有不同生存特征的患者的潜在异构性，并将患者直接分为不同的生存风险组。

保存图片、插图等的外部文件。对象名为nihms-1918361-f0001.jpg

在单独的窗口中打开

图1

根据治疗，40岁以下患者的无病生存率的Kaplan-Meier曲线（左图）和40岁或以上患者的无病生存率的Kaplan-Meier曲线（右图）治疗A:CMF6。治疗B：CMF6+3个周期。治疗C CMF3。治疗D CMF3+3个周期。

三。方法

3.1. 混合生存模型

我们假设整个人口包括 $K（K）$ 不同的生存风险组。每组患者将遵循特定的生存情况。更具体地说，我们假设 $k个第个$ 组具有生存分布 $S公司 (t吨, η_{k个})$ ，具有未知参数的参数形式 $η_{k个}$ ，用于 $k个 = 1, \dots, K（K）$ 在本文中，我们假设每个潜在风险组的生存结果遵循威布尔分布，这是生存分析中常用的分布，因为它具有灵活性和可靠性(廖和刘2019). Weibull分布 $k个第个$ 潜在风险组具有以下形式 $S公司 (t吨, η_{k个}) = 经验 {- (\frac{t吨}{λ_{k个}})^{κ_{k个}}}$ ，其中 $η_{k个} = (κ_{k个}, λ_{k个})^{T型}$ , $κ_{k个}$ 是形状参数 $λ_{k个}$ 是比例参数。

让 $T型$ 表示事件发生的时间 $X（X）$ 表示所有基线协变量，可能是高维的。使用基线协变量将每个患者分为一个存活组 $X（X）$ ，我们引入了一个潜在的组成员身份 $B类$ 并假设

P（P） (T型 > t吨 Ş===================================================================================================== B类 = k个, X（X）) = S公司 (t吨, η_{k个})

和

P（P） (B类 = k个 Ş===================================================================================================== X（X）) = \frac{经验 {克_{k个} (X（X）)}}{\sum_{k个 = 1}^{K（K）} 经验 {克_{k个} (X（X）)}} = {第页}_{k个} (X（X）),

对于 $k个 = 1, \dots, K（K）$ ，其中 $克_{k个} (X（X）)$ 是一个非参数函数，并且 $克_{1} (X（X）) = 0$ （我们将子组1设置为参考组）。因此，潜在的群体成员关系决定了患者应该属于哪个群体。这种成员关系取决于通过非参数分布的基线协变量。显然，所提出的模型意味着 $T型$ 采用混合形式：

P（P） (T型 > t吨 Ş===================================================================================================== X（X）) = \sum_{k个 = 1}^{K（K）} S公司 (t吨, η_{k个}) {第页}_{k个} (X（X）),

以及患者属于风险组的概率 $k个$ 考虑到他或她的基线特征， ${第页}_{k个} (X（X）)$ ，是我们用于分配潜在生存风险组成员的风险分数。使用基线协变量对任何新患者进行未来试验 $X（X） = x个$ ，然后我们将该患者分为风险组 $k个$ 具有最大值 ${第页}_{k个} (X（X）)$ 即最可能的组成员。

3.2. 基于树的模型拟合算法

我们提出了一种基于树的算法来估计组分类函数， ${第页}_{k个} (x个)$ , $k个 = 1, \dots, K（K）$ 不同于传统分类树中的类标签，在我们提出的基于树的方法中，生存风险组的标签， $B类$ ，在数据中未知。因此，我们将其视为缺失数据，并使用EM算法来解决风险组成员资格的延迟问题。我们提出的算法的基本思想是递归地划分协变量空间以优化观测数据的对数似然函数，并将相同的生存风险组成员分配给协变量空间中相同子区域的患者。当只有少数患者包含在子节点中时，递归分区停止。

更具体地说，假设有 $K（K）$ 数据中包含的子组。我们有来自 $n个$ i.i.d.患者，表示为

{{Y（Y）}_{我} = {T型}_{我} \land {C类}_{我}, Δ_{我} = 我 ({T型}_{我} \leq {C类}_{我}), {X（X）}_{我}, 我 = 1, \dots, n个},

哪里 ${C类}_{我}$ 是审查时间。假设审查时间独立于 ${T型}_{我}$ 鉴于 ${X（X）}_{我}$ ，观测数据的对数似然函数由下式给出

我 (η, 克_{2}, \dots, 克_{K（K）}; X（X）, Y（Y）, Δ) = \sum_{我 = 1}^{n个} [Δ_{我} 日志 (\sum_{k个 = 1}^{K（K）} （f） ({Y（Y）}_{我}, η_{k个}) {第页}_{k个} ({X（X）}_{我})) + (1 - Δ_{我}) 日志 (\sum_{k个 = 1}^{K（K）} S公司 ({Y（Y）}_{我}, η_{k个}) {第页}_{k个} ({X（X）}_{我}))]

完整的数据log-likelihood函数如下所示

我_{c（c）} (η, 克_{2}, \dots, 克_{K（K）}; X（X）, Y（Y）, Δ, {B类}_{我}) = \sum_{我 = 1}^{n个} \sum_{k个 = 1}^{K（K）} 我 ({B类}_{我} = k个) [Δ_{我} 日志 {（f） ({Y（Y）}_{我}; η_{k个}) {第页}_{k个} ({X（X）}_{我})} + (1 - Δ_{我}) 日志 {S公司 ({Y（Y）}_{我}; η_{k个}) {第页}_{k个} ({X（X）}_{我})}],

哪里 $η = (η_{1}, \dots, η_{K（K）})^{T型}$ , $（f） (t吨, η_{k个}) = - {S公司}^{^{'}} (t吨, η_{k个})$ .

实现研究非参数函数的数值算法 ${第页}_{k个} (X（X）)$ 为了生成决策树，我们首先为 $η_{k个}$ 和 $克_{k个} (X（X）)$ ，用于 $k个 = 1, \dots, K（K）$ . $η_{k个} ’ 秒$ 选择接近0，但具有不同的值，以确保每个子组的生存分布是可识别的。 $克_{k个} (X（X）)$ 取决于加权logistic回归模型的系数，本节稍后将介绍更多详细信息。加权logistic回归模型系数的起始值选择为0。对于每次拆分，基于数据 $ℬ_{第页 + 1}$ 在当前节点中，我们应用EM算法。在E步骤中，根据观测数据计算基于所有参数的预期对数似然，这相当于计算 ${B类}_{我} = k个$ 根据观测数据，

{q个}_{我 k个}^{(第页 + 1)} = {(\frac{（f） ({Y（Y）}_{我}, η_{k个}) {第页}_{k个}^{(第页)} ({X（X）}_{我})}{\sum_{k个 = 1}^{K（K）} （f） ({Y（Y）}_{我}, η_{k个}) {第页}_{k个}^{(第页)} ({X（X）}_{我})})}^{Δ_{我}} {(\frac{S公司 ({Y（Y）}_{我}, η_{k个}) {第页}_{k个}^{(第页)} ({X（X）}_{我})}{\sum_{k个 = 1}^{K（K）} S公司 ({Y（Y）}_{我}, η_{k个}) {第页}_{k个}^{(第页)} ({X（X）}_{我})})}^{1 - Δ_{我}},

哪里

{第页}_{k个}^{(第页)} ({X（X）}_{我}) = \frac{经验 {克_{k个}^{(第页)} ({X（X）}_{我})}}{\sum_{k个 = 1}^{K（K）} 经验 {克_{k个}^{(第页)} ({X（X）}_{我})}},

对于 $我 \in ℬ_{第页 + 1}$ 和 $k个 = 1, \dots, K（K）$ 通过使用 ${{\hat{η}}^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}}$ 从父节点获取。在M步骤中，对于每个特征 $j个$ , $j个 = 1, \dots, 第页$ ，对于每个潜在分割 $x个$ ，我们优化了目标函数

我_{j个} (η, 克_{2}, \dots, 克_{K（K）}; X（X）, x个) {Ş=====================================================================================================}_{η = η^{(第页)}} = \sum_{我 = 1}^{n个} \sum_{k个 = 1}^{K（K）} {q个}_{我 k个}^{(第页 + 1)} [Δ_{我} 日志 (（f） ({Y（Y）}_{我}, η_{k个})) + (1 - Δ_{我}) 日志 (S公司 ({Y（Y）}_{我}; η_{k个})) + 日志 ({第页}_{k个} ({X（X）}_{我 j个}))]

这相当于拟合加权logistic回归模型

{L（左）}_{j个} (θ_{k个}; X（X）, x个) = \sum_{我 \in ℬ_{第页 + 1}} \sum_{k个 = 1}^{K（K）} {q个}_{我 k个}^{(第页 + 1)} 日志 {\frac{经验 {θ_{0 k个} + θ_{1 k个} 我 ({X（X）}_{我 j个} < x个)}}{\sum_{k个 = 1}^{K（K）} 经验 {θ_{0 k个} + θ_{1 k个} 我 ({X（X）}_{我 j个} < x个)}}}

(1)

基于当前数据 $ℬ_{第页 + 1}$ 更具体地说，加权logistic回归模型中的响应变量是每个受试者的组标签 $我$ ，权重是的后验概率 ${B类}_{我} = k个$ 给定观测数据，自变量是相对于特征变量的指示函数 $j个$ .拟合加权logistic回归模型(1)，我们使用一步Newton-Raphson更新。更具体地说， $θ_{0 k个}^{(第页 + 1)}$ 和 $θ_{1 k个}^{(第页 + 1)}$ 由更新

θ_{k个}^{(第页 + 1)} = (θ_{0 k个}^{(第页 + 1)}, θ_{1 k个}^{(第页 + 1)})^{T型} = θ_{k个}^{(第页)} - {(\frac{\partial^{2} {L（左）}_{j个} (θ; X（X）, x个)}{\partial θ_{k个} \partial θ_{k个}^{T型}} {Ş=====================================================================================================}_{θ_{k个} = θ_{k个}^{(第页)}})}^{- 1} \frac{\partial {L（左）}_{j个} (θ; X（X）, x个)}{\partial θ_{k个}} {Ş=====================================================================================================}_{θ_{k个} = θ_{k个}^{(第页)}}

对于 $k个 = 1, \dots, K（K）$ .

下一步是在优化观测数据对数似然的网格点之间搜索分裂变量和对应的分裂值。对特征变量进行彻底搜索 $j个$ , $j个 = 1, \dots, 第页$ 及其相应的潜在分裂值 $x个$ .给， $x个$ 可以是变量的第20、30、…、80百分位 $j个$ 找到拆分器{ ${j个}^{'}$ , ${x个}^{'}$ }使目标函数最大化(1)。然后我们计算相应的{ $θ_{0 k个}^{^{'}}$ , $θ_{1 k个}^{^{'}}$ }和更新 $克_{k个}^{(第页 + 1)} ({X（X）}_{我})$ 在每次迭代时，使用{ $θ_{0 k个}^{^{'}}$ , $θ_{1 k个}^{^{'}}$ }和分裂{ ${j个}^{'}$ , ${x个}^{'}$ }，

克_{k个}^{(第页 + 1)} ({X（X）}_{我}) = 克_{k个}^{(第页)} ({X（X）}_{我}) \times {θ_{0 k个}^{^{'}} + θ_{1 k个}^{^{'}} 我 ({X（X）}_{我 {j个}^{'}} < {x个}^{'}} .

非参数函数 $克_{k个} (\cdot)$ , $k个 = 1, \dots, K（K）$ ，演示了如何在迭代过程中划分协变量空间。为了估计每个风险组的生存分布，我们考虑了整个数据集 $ℬ$ 也就是说，对于树的每个级别，我们将合并每个子节点中的所有数据点。然后使用Newton-Raphson算法更新未知生存参数 $η_{k个}$ , $k个 = 1, \dots, K（K）$ ，通过计算目标函数的一阶和二阶导数 $我 (, 克_{2}, \dots, 克_{K（K）}; {X（X）}_{我}, X（X）)$ 关于 $η$ 以获得 ${\hat{η}}^{(第页 + 1)}$ 更具体地说，在每个级别上，我们只迭代一次来更新生存参数 $η_{k个}$ , $k个 = 1, \dots, K（K）$ ，基于整个数据集。当每个叶包含的内容不超过 $\sqrt{n个}$ 向后执行修剪过程以获得最佳大小的子树。

在完全生长一棵树之后，我们通过计算一系列子树的BIC来获得最佳大小的子树。BIC标准可以通过以下公式计算

B类 我 C类 = - 2 我 (η, 克_{2}, \dots, 克_{K（K）}; X（X）, Y（Y）, Δ) {Ş=====================================================================================================}_{η = \hat{η}, 克_{2} = {\hat{小时}}_{2}, \dots, 克_{K（K）} = {\hat{小时}}_{K（K）}} + 日志 (n个) (我 e（电子） 一 v（v） e（电子） 秒 + 2 K（K）),

哪里 $我 (η, 克_{2}, \dots, 克_{K（K）}; X（X）, Y（Y）, Δ) {Ş=====================================================================================================}_{η = \hat{η}, 克_{2} = {\hat{小时}}_{2}, \dots, 克_{K（K）} = {\hat{小时}}_{K（K）}}$ 是当前树的观测数据对数似然， $n个$ 是中的观察数 $ℬ$ ,叶子是当前树的叶数，以及 $K（K）$ 是数据中潜在组的数量。当找到与最小BIC对应的树时，修剪停止。

3.3. 选择潜在组的数量

如果有关于潜在群体数量的科学或临床知识，我们可以继续利用这些信息。然而，对此类信息进行验证通常非常困难。因此，为了选择数据中包含的潜在风险组的数量，我们假设数据中存在多个潜在风险组数量选择，并基于每个选择生成一棵树。计算每个树的BIC值。然后选择产生最小BIC的树作为最佳树。我们在这里提到，BIC惩罚修剪树的叶子数，而对级别数的惩罚用于BIC计算，以选择最佳数量的潜在组。BIC计算的这种惩罚选择在教派。4.

3.4. 算法的单调似然性质

在本节中，我们证明了基于中所述算法的观测数据对数似然在迭代过程中单调增加教派。3.2.在不损失通用性的情况下，我们选择任意节点， $𝒩^{(第页)}$ ，英寸 $第页第个$ 迭代，并且对应于该节点的数据集由 $ℬ^{(第页)}$ 假设在 $(第页 + 1)$ 迭代，两个数据集， $ℬ_{1}^{(第页 + 1)}$ 和 $ℬ_{2}^{(第页 + 1)}$ ，基于二进制拆分从 $ℬ^{(第页)}$ ，其中 $ℬ^{(第页)} = ℬ_{1}^{(第页 + 1)} \cup ℬ_{2}^{(第页 + 1)}$ 。略带滥用符号，我们将观察到的数据写成 ${Y（Y）}_{光突发事件} = {{Y（Y）}_{我}, Δ_{我}, {X（X）}_{我}, 我 = 1, 2, \dots, n个}$ ，完整数据为 ${Y（Y）}_{c（c）} = {{Y（Y）}_{我}, Δ_{我}, {X（X）}_{我}, {B类}_{我}, 我 = 1, 2, \dots, n个}$ ，其中 ${B类}_{我}$ 是患者的潜在群体成员 $我$ ，缺少的数据为 ${Y（Y）}_{管理信息系统} = {{B类}_{我}, 我 = 1, 2, \dots, n个}$ 。我们得出以下定理。

定理1 对于第th次迭代中的任何节点及其在(第页+ 1)第次迭代，观测数据的对数似然值在两次连续迭代之间增加。那就是，

我 (η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}; {Y（Y）}_{o个 b条 秒} \in ℬ^{(第页)}) \leq 我 (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{o个 b条 秒} \in ℬ_{1}^{(第页 + 1)}) + 我 (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{o个 b条 秒} \in ℬ_{2}^{(第页 + 1)}),

哪里 $我 (η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}; {Y（Y）}_{o个 b条秒} \in ℬ^{(第页)})$ 表示第th次迭代中的观测数据对数似然性，以及 $我 (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{o个 b条秒} \in ℬ_{1}^{(第页 + 1)})$ ; $我 (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{o个 b条秒} \in ℬ_{2}^{(第页 + 1)})$ 表示在两个相应的子节点中观察到的数据对数似然.

证明定理1在中给出附录.

4 模拟研究

4.1、。模拟设置

我们进行了广泛的仿真研究，以评估我们提出的基于树的方法的性能。我们在第一个模拟研究中考虑了两个潜在风险组，并假设第二个模拟研究的数据包含三个潜在风险群。在每个模拟研究中，评估了三种不同的场景，并将结果与参数混合模型、带有风险分层的Cox-PH模型和生存树进行了比较。我们通过以下方式应用参数混合模型Bussy等人（2019年）.混合形式， $（f） (t吨 Ş===================================================================================================== X（X） = x个) = \sum_{k个 = 1}^{K（K）} π_{k个} (x个) {（f）}_{k个} (t吨; α_{k个})$ ，被认为是生存时间条件概率的模型 $T型$ 给定基线协变量，其中他们假设 $（f） (t吨 Ş===================================================================================================== X（X） = x个)$ 是几何分布或威布尔分布。重量， $π_{k个} (x个)$ ，对于每个混合成分，取决于患者基线特征，并具有逻辑回归函数 $π_{k个} (x个) = \frac{经验 {{x个}^{T型} β_{k个}}}{\sum_{k个 = 1}^{K（K）} 经验 {{x个}^{T型} β_{k个}}}$ EM算法用于处理生存风险组标签的延迟。通过Newton-Raphson算法估计协变量的回归系数，然后将患者分为不同的生存风险组。请注意，此方法可以应用于任何参数分布。带有风险分层的Cox PH模型是Tseng等人（2020年）其中，将包含所有协变量的Cox PH模型应用于训练集，并将生存风险组分类的阈值设置为训练集的中值生存时间（观察样本中生存时间的三分位数用于三种潜在组场景）。该阈值将用于对测试集中的生存数据进行二分，并相应地获得生存风险组。生存树是使用R包实现的，零件，其中，将包含所有协变量的生存树拟合到训练集上，并将训练集的中值生存时间（生存时间的三分位数用于三个潜在组场景）设置为生存风险组分类的阈值。生存风险组是通过使用阈值对测试集中的生存数据进行二分得到的。我们不会采用任何变量选择程序来比较不同的方法。我们提出的混合生存树考虑了生存剖面中未观察到的异质性，并使用EM算法对未观察到的组标签进行建模。我们提出的混合生存树的分裂准则是观测数据的对数似然函数，这与R中的生存树方法不同。因此，我们提出的方法中的递归划分可以更好地适应数据中生存分布的混合。

为了模拟两个潜在亚组，十个基线协变量， ${X（X）}_{1}$ , ${X（X）}_{2}, \dots, {X（X）}_{10}$ 模型中考虑了，其中只有少数具有非零效应。 ${X（X）}_{1}$ , ${X（X）}_{2}, \dots, {X（X）}_{10}$ 由标准正态分布独立生成。每个潜在子组的时间到事件数据均由威布尔分布生成。刻度参数的真实值（即。， $λ_{1}$ 和 $λ_{2}$ )威布尔分布的真值分别被设置为1和4.5。， $κ_{1}$ 和 $κ_{2}$ )分别为1和3。2年生存率高风险组为13.5%，低风险组为91.5%。审查时间由指数分布生成，其中平均值由预先指定的30%审查率校准。高风险组约有40%的个体，其余60%的受试者属于低风险组。亚组特定生存曲线如左图所示图2考虑不同的场景来模拟实际应用中的各种情况。具体来说，我们考虑了 $克 (\cdot)$ 。对于所有情况，我们都将高风险群体作为参考，因此 $克_{1} (X（X）)$ 设置为0。对于场景I.1， $克_{2} (X（X）)$ 设置为线性函数，

克_{2} (X（X）) = 0.4 + 0.2 {X（X）}_{1} - 0.6 {X（X）}_{2} - 0.3 {X（X）}_{三} .

保存图片、插图等的外部文件。对象名为nihms-1918361-f0002.jpg

在单独的窗口中打开

图2

两项模拟研究的真实生存曲线。左图显示了两个潜在群体的真实生存曲线。右图显示了三个潜在群体的真实生存曲线。

在这种情况下，由于参数混合模型是真实的模型，因此预计它将具有更好的性能。对于场景I.2，我们模拟协变量空间上的递归划分，通过考虑指标和交互项来分配风险组成员。因此， $克_{2} (X（X）)$ 设置为

克_{2} (X（X）) = 0.4 + 0.2 我 ({X（X）}_{1} < 0.4, {X（X）}_{三} > 0.2) - 0.6 {X（X）}_{1}^{2} 我 ({X（X）}_{2} > 0.1) + 0.2 我 ({X（X）}_{三} > 0.4 或 {X（X）}_{2} < - 0.2) + 0.2 我 ({X（X）}_{2} < 0.4) .

我们预计我们提出的基于树的方法在这种情况下会表现更好，因为模型指定错误可能是参数混合模型和Cox PH模型中的一个问题。对于场景I.3，我们考虑基于单位圆的生存风险组成员分配的更复杂的决策边界

克_{2} (X（X）) = {X（X）}_{1}^{2} + {X（X）}_{2}^{2} - 1 .

这些设置旨在说明我们提出的基于树的方法的灵活性。

我们考虑使用类似的程序生成三个潜在组的模拟研究的基线协变量和生存数据。从−2到2的截断标准正态分布独立生成十个协变量。注意，截断正态分布仅用于稳定数值计算。审查率定为30%。高、中、低风险组的威布尔分布具有真值1，3.3比例参数分别为7.4，形状参数分别为1、2.5和4。高、中、低风险组的2年生存率分别为13.5%、75.5%和95.5%。亚组特定生存曲线如右图所示图2。在三个潜在群体的模拟研究中，考虑了三个类似的场景。对于情景II.1，约17%的受试者被分配到高风险组，约46%的受试者属于中风险组，其余的受试者属于低风险组。在情景二.2中，高、中、低风险组的受试者比例分别为28%、18%和54%。在场景II.3中，高、中、低风险组的受试者比例分别为29%、35%和36%。非参数函数 $克_{2} (X（X）)$ 和 $克_{三} (X（X）)$ ，分别是情景II.1下的基线协变量的线性组合，

克_{2} (X（X）) = 0.9 + {X（X）}_{1} - 0.8 * {X（X）}_{2} - {X（X）}_{三}, 克_{三} (X（X）) = 0.6 - 0.6 * {X（X）}_{1} - {X（X）}_{2} + 0.7 * {X（X）}_{三} .

对于场景II.2，我们通过设置

克_{2} (X（X）) = - 0.5 + 0.6 * 我 ({X（X）}_{1} < - 1 或 {X（X）}_{2} > 1) - 0.8 * 我 ({X（X）}_{1} > 0 和 {X（X）}_{三} > 0.2), 克_{三} (X（X）) = 0.5 + 0.4 * 我 ({X（X）}_{1} < - 1 或 {X（X）}_{2} > 1) + 0.4 * 我 ({X（X）}_{1} > 0 和 {X（X）}_{三} > 0.2) .

场景II.3的决策边界基于单位圆，我们在这里设置

克_{2} (X（X）) = ({X（X）}_{1}^{2} + {X（X）}_{2}^{2} - 1) * (我 ({X（X）}_{1} > 0.5) + 1), 克_{三} (X（X）) = ({X（X）}_{1}^{2} + {X（X）}_{2}^{2} - 1) * (我 ({X（X）}_{1} \leq 0.5) + 1) .

我们还考虑了数据中没有群体效应的情况，并进行了模拟研究，以在数据只有一个群体的情况下评估我们提出的方法。结果表明，BIC标准可以准确识别数据中只有一个子组。有关此场景的模拟设置和结果的更多详细信息，请参阅补充信息考虑到高风险组和低风险组之间中等差异的两个潜在亚组的另一个场景也用于评估我们提出的方法。关于非参数函数的三种不同设置 $克 (X（X）)$ 已考虑。有关此模拟研究的更多详细信息，请参阅补充信息.

4.2. 仿真结果

我们基于树的算法的第一步是确定数据中存在的潜在组的数量。为此，我们假设数据中潜在组的数量不同，并为每个选择的潜在组数量生成一棵树。对每棵树计算对树的层数有惩罚的BIC，然后使用它来确定最佳的潜在组数。我们报告了基于1000个重复的模拟结果，并考虑了300、500、1000和2000个样本大小。对于两个潜在群体的场景，潜在群体数量的候选选择为1、2和3。当样本大小等于300时，在70%以上的模拟数据集中选择了正确的潜在组数。随着样本量增加到1000，几乎90%的复制品都能正确选择潜在亚组的数量。对于三个潜在群体的场景，潜在群体数量的候选选择范围为1到4。当样本量为300时，BIC可以从大约70%的模拟数据集中选择合适数量的潜在组。随着样本量增加到1000，比例增加到80%。样本量为2000，超过94%的模拟数据集选择了正确数量的潜在群体。

表1和和22总结中位数预测精度以及中位数绝对偏差（MAD），用于两个和三个潜在组的模拟研究。中值绝对偏差定义为与数据中值的绝对差值的中值，即。， $M（M） A类 D类 = 米 e（电子） d日我一 n个 (Ş===================================================================================================== {X（X）}_{我} - 米 e（电子） d日我一 n个 (X（X）) Ş=====================================================================================================)$ ，其中 $X（X）$ 是定量数据的单变量样本。我们提出的方法的结果是基于BIC修剪的树得到的。对于数值稳定性，对于两个潜在群体场景，威布尔分布估计值限制为（0，7]，并且限制为当实施Newton-Raphson算法时，对于三个潜在组场景（0，15]）。我们还根据平均深度和平均叶数报告了树结构。通过将从训练集获得的决策规则应用于独立生成的样本大小为10的验证集，计算预测准确度，000.最佳准确率计算为1–贝叶斯错误率，其中贝叶斯误差率由以下公式计算 $1 - E类 (\underset{k个}{最大值} P（P） (B类 = k个 Ş===================================================================================================== X（X）))$ 一般来说，预测精度随着样本量的增加而增加，并产生中值绝对偏差的下降趋势。对于大样本量，预测精度也接近最佳准确率。在检查树结构时，即使对于一些复杂的场景，基于树的方法也会生成相对简单的树，这有助于可视化和解释基线协变量如何确定生存风险组分类。对于场景I.1和II.1，我们假设非参数函数中存在线性关系g（·），也就是说，正确指定了参数混合模型。因此，参数混合模型的预测精度优于基于树的方法。具有风险分层的Cox PH模型的结果略好于生存树的结果，这可能也是由于在设置 $克 (\cdot)$ 然而，对于参数混合模型不是真实模型的其他场景，基于树的方法优于参数混合模型。与具有风险分层和生存树的Cox PH模型相比，我们的方法也获得了更好的结果。一般来说，生存树比具有风险分层的Cox PH模型表现得更好，这是合理的，因为生存树可以更好地捕捉非参数关系 $克 (\cdot)$ .

表1

两种潜在群体情景的模拟研究结果

N个	预测准确性				树状结构
N个	我们的方法药物累积。（摩洛哥迪拉姆）	混合物模型药物累积。（摩洛哥迪拉姆）	Cox PH模型药物累积。（摩洛哥迪拉姆）	生存树药物累积。（摩洛哥迪拉姆）	平均深度	树叶平均数量
场景I.1。最佳准确率为0.649。
300	0.590 (0.012)	0.614 (0.013)	0.592（0.012）	0.537 (0.032)	2.79	8.03
500	0.600 (0.014)	0.627 (0.008)	0.606 (0.009)	0.588 (0.023)	3.30	10.36
1000	0.613 (0.009)	0.637 (0.006)	0.612 (0.006)	0.611 (0.011)	3.95	15.12
2000	0.623 (0.007)	0.642 (0.004)	0.618 (0.005)	0.615（0.009）	4.44	20.90
场景I.2。最佳准确率为0.668。
300	0.602 (0.012)	0.593 (0.017)	0.535 (0.015)	0.527（0.043）	2.55	6.97
500	0.606 (0.011)	0.602 (0.011)	0.544 (0.013)	0.567 (0.040)	3	8.88
1000	0.620 (0.015)	0.607 (0.006)	0.553 (0.010)	0.604 (0.009)	3.76	13.68
2000	0.642 (0.012)	0.610 (0.004)	0.562 (0.007)	0.606（0.006）	4.39	20.14
场景I.3。最佳准确率为0.721。
300	0.623 (0.019)	0.561 (0.031)	0.502 (0.005)	0.586 (0.018)	3.51	10.61
500	0.664 (0.038)	0.589 (0.021)	0.506 (0.005)	0.620 (0.012)	4.28	15.35
1000	0.702 (0.010)	0.614（0.008）	0.507 (0.005)	0.622 (0.006)	5.01	23.31
2000	0.709 (0.006)	0.619 (0.004)	0.507 (0.004)	0.627（0.004）	5.32	29.71

在单独的窗口中打开

混合模型：参数混合模型Cox PH模型：具有风险分层的Cox-PH模型Med accu.：中值精度MAD：中值绝对偏差平均深度：平均深度平均无叶数：平均叶数

表2

三种潜在群体情景的模拟研究结果

N个	预测准确性				树状结构
N个	我们的方法药物累积。（摩洛哥迪拉姆）	混合物模型药物累积。（摩洛哥迪拉姆）	Cox PH模型药物累积。（摩洛哥迪拉姆）	生存树药物累积。（摩洛哥迪拉姆）	平均深度	树叶平均数量
场景II.1。最佳准确率为0.663。
300	0.407 (0.044)	0.434 (0.098)	0.479 (0.010)	0.388 (0.088)	2.13	5.17
500	0.421 (0.034)	0.489 (0.041)	0.485 (0.007)	0.413 (0.103)	4.11	15.58
1000	0.434 (0.045)	0.550 (0.010)	0.486 (0.007)	0.469 (0.044)	4.70	23.03
2000	0.467 (0.035)	0.592 (0.006)	0.488 (0.006)	0.522 (0.025)	5.10	30.38
场景II.2。最佳准确率为0.551。
300	0.484 (0.035)	0.328 (0.026)	0.177 (0.009)	0.177 (0.059)	3.81	12.61
500	0.502 (0.026)	0.340 (0.028)	0.177 (0.004)	0.191 (0.038)	4.02	15.01
1000	0.515 (0.020)	0.368 (0.020)	0.179 (0.002)	0.258 (0.018)	4.45	20.40
2000	0.523 (0.016)	0.428 (0.014)	0.187 (0.003)	0.324 (0.008)	4.92	28.10
场景II.3。最佳准确率为0.611。
300	0.390 (0.035)	0.365 (0.071)	0.358 (0.004)	0.401 (0.021)	3.96	13.01
500	0.427 (0.037)	0.404 (0.050)	0.359 (0.004)	0.407 (0.022)	4.33	16.18
1000	0.473 (0.026)	0.445 (0.016)	0.359 (0.003)	0.407 (0.020)	5.15	25.20
2000	0.498 (0.015)	0.449（0.016）	0.359 (0.004)	0.408 (0.013)	5.96	38.47

在单独的窗口中打开

混合模型：参数混合模型Cox PH模型：具有风险分层的Cox PH模型Med accu:中值精度MAD:中值绝对偏差Ave Depth：平均深度Ave no leaves:平均叶片数

5 IBCSG数据

我们将建议的方法应用于中讨论的IBCSG数据教派。2研究患者生存结果的潜在异质性，并研究与这种异质性相关的重要变量。

为了探讨不同治疗程序下患者的异质性，我们将我们的方法应用于每种治疗下患者的数据集。根据BIC，在治疗B（CMF6+3周期）和治疗C（CMF3）的患者中检测到两个潜在亚群，在治疗A（CMF6）和治疗D（CMF3+3周期）的患者之间未发现任何潜在亚群。分组选择结果总结如下表3

表3

基于BIC准则的群组选择结果

治疗组	潜在群体数量
治疗组	1	2	三
A类	257.5	278	279
B类	258.2	203.5	234.3
C类	225	201.7	219.5
D类	250	250.5	262

在单独的窗口中打开

接下来，我们通过假设数据中存在两个潜在组来分析接受治疗B（CMF6+3周期）和C（CMF3）的患者的数据。图3显示了治疗B（CMF6+3周期）和治疗C（CMF3）下患者的拟合树。我们还绘制了每片叶子中患者的Kaplan-Meier曲线，该曲线位于每片叶子的下方。基于图3我们发现，对于接受治疗B（CMF6+3周期）的患者，潜在风险组成员分配与阳性淋巴结数量、身体健康、食欲和情绪相关。当观察接受治疗C（CMF3）的患者时，年龄、阳性淋巴结数量、身体健康状况和情绪是潜在风险组成员分配的重要因素。对于接受C类治疗（CMF3）的患者，在Colleoni等人（2002年）：老年患者的生存经验明显更好。绘制的决策树图3基于BIC准则进行修剪，有助于防止过拟合。

保存图片、插图等的外部文件。对象名称为nihms-1918361-f0003.jpg

在单独的窗口中打开

图3

IBCSG乳腺癌试验中接受治疗B（CMF6+3周期）和治疗C（CMF3）的患者的混合生存树。注：。上部面板：接受治疗B（CMF6+3周期）的患者的决策树。下面板：接受治疗C（CMF3）的患者的决策树曲线。

按照中所述的相同程序，将结果与参数混合模型和具有风险分层的Cox PH模型进行比较教派。4根据我们提出的方法，对于接受B组治疗（CMF6+3周期）的患者，高风险组和低风险组的Weibull分布估计值分别为shape=1.86，scale=0.33，shape=1.45，scale=2.48。对于接受C组（CMF3）治疗的患者，高风险组和低风险组的Weibull分布估计值分别为shape=1.82 scale=0.35，shape=1.28，scale=2.67。治疗B（CMF6+3周期）和治疗C（CMF3）患者的Kaplan-Meier曲线如图4。我们还根据从我们的方法中获得的估计Weibull分布参数绘制了高风险组和低风险组的生存分布图4。对于两种治疗中的高危组，我们的方法得出了与威布尔分布相似的估计值，但尾部只有少数受试者。由于相对稀疏的观察结果，低风险组的差异似乎很大。还使用logrank检验比较了两个潜在群体的生存分布。请注意，logrank测试在这个实际数据应用程序中作为一个特别的度量，是我们试图比较两个潜在群体的生存分布，我们不会将logrank检验视为一个正式测试。logrank检验的结果表明，对于接受治疗B（CMF6+3周期）和治疗C（CMF3）的患者，两个潜在组的生存分布显著不同（p值均小于0.0001）。

保存图片、插图等的外部文件。对象名为nihms-1918361-f0004.jpg

在单独的窗口中打开

图4

接受治疗B（CMF6+3周期）和治疗C（CMF3）的患者的Kaplan-Meier曲线

为了与参数混合模型、带风险分层的Cox-PH模型和生存树方法进行比较，我们将这些方法应用于同一数据集，并报告了生存风险组分类的一致性。对于具有风险分层的Cox PH模型和生存树方法，我们使用从观察样本中获得的生存时间的第一个四分位作为阈值，将患者分为高风险组和低风险组，因为我们无法从观察样本计算中位生存时间，并且使用生存时间的第四分位进行分类，在预测方面取得了最高的一致性与我们的方法相比，子组成员身份。对于治疗B（CMF6+3周期），大约94.4%的患者具有参数混合模型和基于树的方法预测的相同风险组成员。在比较Cox PH模型与风险分层和基于树的方法时，48.8%的患者被分配到相同的生存风险组。当比较生存树方法和我们的方法时，所有患者的预测生存风险组成员相同。对于接受C类治疗（CMF3）的患者，在应用参数混合模型和基于树的方法时，有61.5%的患者具有相同的预测生存风险组成员。比较分层Cox PH模型和我们的方法，46.8%的患者被划分为同一风险组。与生存树方法相比，61.9%的患者被划分为同一风险组。

对于治疗B（CMF6+3周期）的模型，基于参数混合模型和具有风险分层的Cox PH模型的结果，正节点数和情绪对潜在风险组成员分配有显著影响。根据混合生存树，这两个协变量也很重要图3在混合生存树中，第一个分裂是大于4个阳性节点的数量，将181名肿瘤中阳性节点少于4个的患者和其余67名肿瘤中有4个以上阳性节点的患者分开。根据44.5的心境值进行第二次分割，这意味着对于肿瘤中有4个以上阳性淋巴结的患者，其中17个心境小于44.5的患者包括一组，其余50名患者的生存情况被认为相似。

对于具有处理C（CMF3）的模型，我们发现在参数混合模型和具有风险分层的Cox PH模型中，年龄和正节点数与潜在群成员分配显著相关，这也与从混合生存树中看到的结果一致图3。混合生存树中的第一次分裂是在43岁时，第二级分裂是根据年龄和肿瘤中阳性节点的数量进行的。

在生存树中，接受治疗B（CMF6+3周期）和C（CMF3）的患者的树非常相似。根据肿瘤中小于4的阳性节点数，从两棵树中进行第一次分割。对于肿瘤中有4个以上阳性淋巴结的患者，第二次分割是基于接受治疗B（CMF6+3周期）的患者的情绪值44.5和接受治疗C（CMF3）的患者情绪值49.5。

研究这些变量如何将患者分为不同的风险组更为有趣。虽然我们发现了决定潜在风险组成员分配的相似协变量集，但基于参数混合模型方法和Cox PH模型，很难检测协变量之间的交互作用。这两种方法只能线性描述协变量和风险组之间的关系。相比之下，混合生存树提供了基线协变量如何确定不同生存风险组的更直接的可视化。例如，在将接受治疗B（CMF6+3周期）的患者分为不同风险组时，阳性淋巴结数量、身体健康和食欲之间存在交互作用，以及年龄、阳性淋巴结数、，为C组（CMF3）患者分配风险组成员时的身体健康和情绪。

生存风险分类的目的是优化风险组之间生存结果的异质性。从这个角度出发，通过对数秩检验来评估由参数混合模型、带有风险分层的Cox-PH模型和生存树方法估计的生存分布之间的差异，以进一步将这些方法与混合生存树进行比较。结果表明，对于治疗B（CMF6+3周期）和治疗C（CMF3）的所有方法，高风险组和低风险组的两种生存分布存在显著差异（p值均<0.05）。对于治疗B（CMF6+3周期）的患者，当测试混合生存树估计的生存分布之间的差异时，logrank检验统计量的值为50.1，对于参数混合模型、Cox PH模型和生存树方法，相应的检验统计量为48.8、4.5和50.1，分别是。这意味着混合生存树更善于区分接受治疗B（CMF6+3周期）的患者之间不同的生存状况。对于接受治疗C（CMF3）的患者，混合生存树、参数混合模型、Cox PH模型和生存树方法的logrank检验统计值分别为11.5、15.9、10.3和34.4。从参数混合模型、具有风险分层的Cox PH模型和生存树方法获得的Kaplan-Meier曲线也绘制在图4为了更好地说明所有方法之间分类结果的差异以及高风险组和低风险组患者生存分布的差异。

根据评审员的意见，我们还分析了数据，其中包括那些缺少基线QOL测量但在最近一次就诊时有可用信息的数据。数据共包含1475名患者。数据中的缺失与生活质量测量有关。其中322名患者在任何就诊时都没有记录生活质量信息，因此直接从分析中删除。在其余患者中，962名患者在基线检查时进行了生活质量测量，其中191名患者缺少基线检查时的生活质量测量值，但这些测量值在最近一次就诊时可用。我们在接受基于BIC的治疗B、C和D的患者中检测到两个潜在组。对于接受治疗B和C的患者，预测的组标签与我们在本节中报告的结果高度一致。根据logrank检验，接受治疗D的患者的生存情况不显著。有关更多详细信息，请参阅补充信息.

6 结束语

在本文中，我们提出了一种基于树的算法来探索具有不同生存特征的患者的潜在异质性，并进一步将患者分为不同的生存风险组。利用潜在群成员关系，使用EM算法对数据中的混合成分进行建模。我们在递归划分的框架内提出了一个新的分裂准则。通过优化每次拆分中的观测数据对数似然函数，我们最大限度地提高了节点内的同质性和节点间的异质性。我们提出的基于树的算法能够同时估计未知生存分布和预测潜在子群成员。还提出了一个简单且可解释的树状结构，以描述基线协变量如何确定患者中未观察到的异质性。仿真研究表明，我们提出的方法在各种情况下都能很好地工作，特别是当潜在子组成员分配通过非线性关系依赖于基线协变量时。在我们提出的方法中，我们假设每个亚组中的生存分布为Weibull，并且在模拟研究中，我们将结果与不需要此假设的其他方法（如具有风险分层的Cox PH模型）进行比较。请注意，这些比较可能并不完全公平，因为模拟研究中的场景都符合我们的假设。另一方面，刘和廖（2020）结果表明，当研究人群由具有不同生存模式的混合亚群组成时，具有两个或三个威布尔分布分量的混合PH模型与Cox PH模型表现相似，这表明混合威布尔模型的灵活性。

使用我们的方法，我们可以根据患者的生存情况直接将患者分为不同的风险组（例如，高风险组与低风险组）。由于确定的类别具有不同的生存分布，因此每个类别都具有临床意义，对应于生存轨迹较长或较短的患者。因此，获得的分类可以用于区分患者的亚组，至少在以下方向上是这样的。首先，获得的潜在类别可用于未来临床试验的患者招募。例如，我们可以从高危人群中招募更多患者来授权试验。另一个潜在的应用如实际数据应用所示，我们的方法可用于确定一组患者，与其他患者相比，这些患者可能从一种治疗中受益更多，并基于选定的重要协变量探索该组患者的基线特征。

还可以考虑对所提出的基于树的方法进行一些扩展。在我们提出的算法中，潜在子群的生存分布是通过具有未知参数的威布尔分布来建模的。很容易将我们的算法推广到其他参数分布，如指数分布和对数正态分布。此外，可以进一步削弱参数生存分布假设，也可以考虑Cox PH模型。除了生存分布假设外，识别与生存风险组分类相关的重要协变量也是一个重要问题，特别是当数据集中使用大量协变量时。尽管树结构提供了关于协变量效应的直接解释，但反映变量重要性的汇总统计仍然是可取的。我们未来工作的一个潜在方向是将集成方法（如打包）纳入当前基于树的算法中。随机森林中的可变重要性计算可以很容易地应用到我们的算法中，并且可以提高预测性能。当数据中存在大量协变量时，由于样本量的原因，所提出的基于树的算法可能不稳定，一种策略是遵循以下两步方法廖等（2020）首先使用机器学习技术选择重要的初始变量，然后将选择的变量应用到我们提出的算法中。

此外，可以在我们提出的基于树的方法框架内研究患者的治疗反应异质性。通过假设患者的潜在异质性与异质性生存和治疗反应相关，可以建立一种半参数方法来同时估计患者的基线生存和治疗效果。使用基于树的方法检查这些问题的潜在子组识别目前正在调查中。

附录

证明我们将完整的数据日志写成

我_{c（c）} (η, 克_{2}, \dots, 克_{K（K）}; {Y（Y）}_{c（c）}) = \sum_{我 = 1}^{n个} \sum_{k个 = 1}^{K（K）} 我 ({B类}_{我} = k个) [Δ_{我} 日志 {（f） ({Y（Y）}_{我}; η_{k个}) {第页}_{k个} ({X（X）}_{我})} + (1 - Δ_{我}) 日志 {S公司 ({Y（Y）}_{我}; η_{k个}) {第页}_{k个} ({X（X）}_{我})}] .

为了表明观测数据的对数似然在连续迭代中增加，我们计算了两次连续迭代之间观测数据对数似然的差异，如下所示

我 (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{o个 b条 秒} \in ℬ_{1}^{(第页 + 1)}) + 我 (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{o个 b条 秒} \in ℬ_{2}^{(第页 + 1)}) - 我 (η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}; {Y（Y）}_{o个 b条 秒} \in ℬ^{(第页)}) = {E类 [我_{c（c）} (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{c（c）}) Ş===================================================================================================== {Y（Y）}_{o个 b条 秒} \in ℬ_{1}^{(第页 + 1)}, η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}] + E类 [我_{c（c）} (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{c（c）}) Ş===================================================================================================== {Y（Y）}_{o个 b条 秒} \in ℬ_{2}^{(第页 + 1)}, η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}] - E类 [我_{c（c）} (η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}; {Y（Y）}_{c（c）}) Ş===================================================================================================== {Y（Y）}_{o个 b条 秒} \in ℬ^{(第页)}, η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}]} - {E类 [我_{c（c）} (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{米 我 秒} Ş===================================================================================================== {Y（Y）}_{o个 b条 秒}) Ş===================================================================================================== {Y（Y）}_{o个 b条 秒} \in ℬ_{1}^{(第页 + 1)}, η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}] + E类 [我_{c（c）} (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{米 我 秒} Ş===================================================================================================== {Y（Y）}_{o个 b条 秒}) Ş===================================================================================================== {Y（Y）}_{o个 b条 秒} \in ℬ_{2}^{(第页 + 1)}, η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}] - E类 [我_{c（c）} (η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}; {Y（Y）}_{米 我 秒} Ş===================================================================================================== {Y（Y）}_{o个 b条 秒}) Ş===================================================================================================== {Y（Y）}_{o个 b条 秒} \in ℬ^{(第页)}, η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}]} \equiv (我) - (我 我) .

数量( $我$ )根据中所述算法为非负教派。三主要论文的 ${η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}}$ 满足

E类 [我_{c（c）} (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{c（c）}) Ş===================================================================================================== {Y（Y）}_{o个 b条 秒} \in ℬ_{1}^{(第页 + 1)}, η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}] + E类 [我_{c（c）} (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{c（c）}) Ş===================================================================================================== {Y（Y）}_{o个 b条 秒} \in ℬ_{2}^{(第页 + 1)}, η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}] \geq E类 [我_{c（c）} (η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}; {Y（Y）}_{c（c）}) Ş===================================================================================================== {Y（Y）}_{o个 b条 秒} \in ℬ^{(第页)}, η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}] .

对于( $我我$ )，我们有

E类 [我_{c（c）} (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{米 我 秒} Ş===================================================================================================== {Y（Y）}_{o个 b条 秒}) Ş===================================================================================================== {Y（Y）}_{o个 b条 秒} \in ℬ_{1}^{(第页 + 1)}, η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}] - E类 [我_{c（c）} (η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}; {Y（Y）}_{米 我 秒} Ş===================================================================================================== {Y（Y）}_{o个 b条 秒}) Ş===================================================================================================== {Y（Y）}_{o个 b条 秒} \in ℬ_{1}^{(第页 + 1)}, η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}] = E类 [日志 {\frac{{L（左）}_{c（c）} (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{米 我 秒} Ş===================================================================================================== {Y（Y）}_{o个 b条 秒})}{{L（左）}_{c（c）} (η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}; {Y（Y）}_{米 我 秒} Ş===================================================================================================== {Y（Y）}_{o个 b条 秒})}} Ş===================================================================================================== {Y（Y）}_{o个 b条 秒} \in ℬ_{1}^{(第页 + 1)}, η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}] \leq 日志 (E类 [\frac{{L（左）}_{c（c）} (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{米 我 秒} Ş===================================================================================================== {Y（Y）}_{o个 b条 秒})}{{L（左）}_{c（c）} (η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}; {Y（Y）}_{米 我 秒} Ş===================================================================================================== {Y（Y）}_{o个 b条 秒})} Ş===================================================================================================== {Y（Y）}_{o个 b条 秒} \in ℬ_{1}^{(第页 + 1)}, η^{(第页)}, 克_{2}^{(第页)}, \dots, 克_{K（K）}^{(第页)}]) = 日志 \int_{𝒴_{米 我 秒} ({Y（Y）}_{o个 b条 秒})} {L（左）}_{c（c）} (η^{(第页 + 1)}, 克_{2}^{(第页 + 1)}, \dots, 克_{K（K）}^{(第页 + 1)}; {Y（Y）}_{米 我 秒} Ş===================================================================================================== {Y（Y）}_{o个 b条 秒}) d日 {Y（Y）}_{米 我 秒} = 0,

哪里 ${L（左）}_{c（c）} (\cdot) = 经验 {我_{c（c）} (\cdot)}$ 是完整的数据似然函数。

根据Jensen不等式，上述不等式成立。数据集的类似参数 $ℬ_{2}^{(第页 + 1)}$ 很容易获得。因此( $我我$ )为非正，这意味着观测数据的似然函数在连续迭代中增加。□

工具书类

Bussy S、Guilloux A、Gaíffas S、Jannot A-S（2019年）C-mix：删失持续时间的高维混合模型，及其在遗传数据中的应用.统计方法医学研究 28(5):1523–1539 [公共医学][谷歌学者]
Chi Y-Y，Ibrahim JG（2006）多变量纵向和多变量生存数据的联合模型.生物识别 62(2):432–445 [公共医学][谷歌学者]
Ciampi A、Bush R、Gospodarowicz M、Till J（1981）非霍奇金淋巴瘤患者生存经验相关预后因素分类方法：基于一系列982例患者：1967年至1975年.癌症 47(三):621–627 [公共医学][谷歌学者]
Ciampi A，Chang C，Hogg S，McKinney S（1987）递归分区：生物统计学中勘探数据分析的通用方法。生物统计学纽约州施普林格，第23-50页[谷歌学者]
Colleoni M、Litman H、Castiglione-Gertsch M、Sauerbrei W、Gelber R、Bonetti M、Coates A、Schumacher M、Bastert G、Rudenstam C等（2002）乳腺癌辅助化疗的持续时间：两项随机试验的联合分析，研究三个疗程与六个疗程的cmf.英国癌症杂志 86(11):1705–1714[PMC免费文章][公共医学][谷歌学者]
考克斯DR（1972）回归模型和生命表.J Roy Stat Soc Ser B（Methodol公司） 34(2):187–202[谷歌学者]
Davis RB，Anderson JR（1989年）指数生存树.统计医学 8(8):947–961 [公共医学][谷歌学者]
《永别了VT》（1982）使用混合模型分析长期幸存者的生存数据.生物识别 38(4):1041–1046 [公共医学][谷歌学者]
Gordon L，Olshen RA（1985）树结构生存分析.癌症治疗代表 69(10):1065–1069 [公共医学][谷歌学者]
易卜拉欣·N，库杜斯·A（2009）多变量生存数据和竞争风险预测分类的决策树.INTECH Open Access Publisher，伦敦[谷歌学者]
国际乳腺癌研究小组（1996年）绝经前淋巴结阳性乳腺癌患者辅助化疗的持续时间和再次应用.临床肿瘤学杂志 14(6)：1885年至1894年[公共医学][谷歌学者]
Kannel W，McGee D（1979）糖尿病和糖耐量是心血管疾病的危险因素：弗雷明翰研究.糖尿病护理 2(2):120–126 [公共医学][谷歌学者]
Kannel WB、Feinleib M、McNamara PM、Garrison RJ、Castelli WP（1979年）家庭冠心病调查：弗雷明翰后代研究.美国传染病学期刊 110(三)：281–290[公共医学][谷歌学者]
Kaplan EL，Meier P（1958）不完全观测的非参数估计.美国统计协会 53(282):457–481[谷歌学者]
Larson MG，Dinse GE（1985年）竞争风险数据回归分析的混合模型.J Roy Stat Soc Ser C（应用统计） 34(三):201–211[谷歌学者]
LeBlanc M，Crowley J（1992年）截尾生存数据的相对风险树.生物识别 48(2):411–425 [公共医学][谷歌学者]
Liao JJ、Farooqui MZ、Marinello P、Hartzel J、Anderson K、Ma J、Gause CK（2020）使用人工智能工具回答重要的临床问题：keynote-183多发性骨髓瘤经验.康普临床试验第106179页。[公共医学][谷歌学者]
廖俊杰，刘广发（2019）临床试验中用于拟合时间-事件数据的灵活参数生存模型.药物统计 18(5):555–567 [公共医学][谷歌学者]
刘广发、廖俊杰（2020）比例风险约束下基于柔性混合模型的时间-事件数据分析.生物医药统计杂志 30(5):783–796 [公共医学][谷歌学者]
Marubini E、Morabito A、Valsecchi M（1983）预测因素和风险组：使用适用于截尾生存数据的算法得出的一些结果.统计医学 2(2):295–303 [公共医学][谷歌学者]
Moradian H、Larocque D、Bellavance F（2017）存活森林中的L1分裂规则.终身数据分析 23(4):671–691 [公共医学][谷歌学者]
沈J，何X（2015）用结构化逻辑正态混合模型进行子群分析的推断.美国统计协会 110(509):303–312[谷歌学者]
孙毅、邱S-H、王M-C（2019）Roc-guided生存树和集合.生物识别 [PMC免费文章][公共医学][谷歌学者]
曾Y-J，王H-Y，林T-W，卢J-J，谢C-H，廖C-T（2020）晚期口腔癌患者生存风险分层的机器学习模型的开发.JAMA Netw开放三(8)：e2011768–e2011768[PMC免费文章][公共医学][谷歌学者]
Vergara P、Tzou WS、Tung R、Brombin C、Nonis A、Vaseghi M、Frankel DS、Di Biase L、Tedrow U、Mathuria N等（2018）识别室性心动过速消融术患者生存和复发风险的预测评分：i-vt评分.循环心律失常电生理 11(12)：e006730[PMC免费文章][公共医学][谷歌学者]
Wilson PW、D'Agostino RB、Levy D、Belanger AM、Silbershatz H、Kannel WB（1998）利用危险因素分类预测冠心病.循环 97(18):1837–1847 [公共医学][谷歌学者]
曾德、毛立、林德（2016）基于区间相关数据的半参数变换模型的极大似然估计.生物特征 103(2):253–271[PMC免费文章][公共医学][谷歌学者]
周毅，McArdle JJ（2015）生存树和生存集合方法的原理和应用.心理测量学 80(三):811–833[PMC免费文章][公共医学][谷歌学者]

用于癌症风险分类的混合生存树

贾碑林

曾东林

Jason J.Z.廖

广汉F.Liu

谭显明

郭庆雕

约瑟夫·易卜拉欣

摘要

1 介绍

2 IBCSG乳腺癌试验

三。 方法

3.1. 混合生存模型

3.2. 基于树的模型拟合算法

3.3. 选择潜在组的数量

3.4. 算法的单调似然性质

4 模拟研究

4.1、。模拟设置

4.2. 仿真结果

表1

表2

5 IBCSG数据

表3

6 结束语

附录

工具书类

三。方法