引用本文： 陈浩然, 刘夏阳, 王敏, 杨林, 王嘉阳, 孙海霞, 段永恒, 吴旭生, 尚丽, 钱庆, 和晓峰, 李姣. 机器学习模型在非比例风险生存资料中的应用及案例实践. 中国循证医学杂志, 2024, 24(9): 1108-1116. doi: 10.7507/1672-2531.202401190 复制
既包含结局事件又包含生存时间的数据被称为生存资料，又称为时间-事件结局数据。它因同时具有定量和定性的双重性质，可提供比单纯结局事件更加丰富的数据信息。由于生存资料的这种特殊性，针对它的独特统计分析方法被称为生存分析。生存分析在帮助了解疾病的自然史、评估治疗方法的效果、预测患者的生存期，以及指导临床决策等方面有着重要的意义[1]。
目前的生存分析方法大致包括三个方面：① 生存过程的描述（如Kaplan-Meier曲线）；② 生存过程的比较（如log-rank检验）；③ 影响因素分析及生存预测（如Cox比例风险回归模型，简称Cox回归）[2]。log-rank检验和Cox回归模型是生存资料统计分析中最常用的分析策略，但是它们应用有着重要的前提条件—等比例风险假设，而这一要求在既往应用过程中常被忽视：严若华等[3]发现在PubMed检出的使用Cox回归模型的文章中，仅有0.34%（413/121 342）提及了比例风险假定。美国心脏协会发布的心血管医学统计报告中明确建议需要测试并报告Cox回归模型是否违反了等比例危险假定[4]。但是在实际数据中由于大量噪声，很难满足这个假设条件。鉴于此，一些针对不满足等比例风险假定生存资料的统计分析方法已经提出：① 基于log-rank检验或者Kaplan-Meier检验用于生存过程比较的衍生方法[5]；② 基于传统的Cox回归模型用于生存预后影响因素分析的改良方法或其他创新性方法[2]（表1）。然而，这些针对非比例风险生存数据提出的改良分析方法，从根本上说是为小规模、低维度的临床研究数据服务的。尽管它们大多有着可解释性强和易于实施等优点，但是面对大样本数据、高维特征数据或需要探索更加复杂关联形式的分析要求时，它们使用的范围就相对有限[6]。
表1
传统的非比例风险生存分析方法
随着精准医学时代的到来，可以获得的医学数据的规模和特征的维度急剧增加，因此对分析方法提出了更高的要求[7]。作为人工智能重要分支之一的机器学习模型（包括人工神经网络）正越来越渗透于医药健康研究中[8]。针对生存分析的问题，许多研究者也开发了相关机器学习模型来分析生存预后问题[9]。机器学习模型对于数据分布的限制极少，且有更大的假设空间和拟合效果，另外这些机器学习模型在处理高维特征生存大数据或者存在竞争风险生存数据方面也有着独特的优势[6]。许多机器学习生存模型是对Cox回归模型的改进（如DeepSurv[10]和Modified-DeepSurv[11]），尽管它们在解决Cox回归模型中的非线性和交互作用等挑战方面表现出了有效性，但对于存在非比例风险的生存数据的适用性却受到限制。
本文在解读NPH等相关概念之后，重点概述了可以处理不满足等比例风险假定生存资料的机器学习模型，并进行了基于机器学习模型的脑卒中患者死亡风险的案例研究，以期推动人工智能在非比例风险生存资料分析中的应用。
1 相关概念
1.1 非比例风险
Cox回归模型使用广泛的主要原因在于其没有对生存时间和基础风险函数的分布形式做任何的要求，且可以方便的求得风险比（hazard ratios，HRs）以衡量协变量的影响效应。但是Cox回归模型需要满足等比例风险的假设条件，在该条件下求得HRs才会在不同的生存时间保持一致。Cox回归模型是一个半参数模型，在等比例假设下，HRs计算公式如下：
 |
等比例风险假定即HRs与时间无关，其不随时间变化而变化。比如在针对某种药物开展的随机对照试验中，试验组和对照组的HRs在任何的随访时间都是相等的，即药物的治疗效应对生存率的影响不随时间的改变而改变。而当HRs随时间变化，则为NPH。NPH下计算的HRs仅代表随访期间估计的HRs的平均水平，但是这种平均估计是不准确的，尤其对于存在交叉治疗效应的生存数据[12]。
1.2 非比例风险检验
NPH检验是生存资料统计分析过程中的一个重要步骤，既往研究已经揭示了不同的方法以探索是否满足等比例风险假定，本研究在既往研究的基础上建议依据变量类型的不同而采取不同的方法进行检验[2,3,12,13]：① 连续型变量：Schoenfeld残差图法、Martingale残差图法、Schoenfeld残差和时间线性相关分析法、时依协变量假设检验法；② 分类型变量：Kaplan-Meier生存曲线图、log（-log（生存概率））与log（生存时间）作图法（log-log作图法）、Schoenfeld残差图法、Schoenfeld残差和时间线性相关分析法。此外，建议针对每个协变量至少使用两种以上的策略以判断其是否满足等比例风险假设。
2 基于机器学习算法的非比例风险模型
传统的生存分析统计方法在数据规模较小、特征维度较低及算力有限的情况下尤为适用，且由于它们能提供协变量的参数估计而可解释性更强[14]。以针对非比例风险生存资料衍生的Cox回归模型为例（表1），它们依旧可以提供HRs，从而很容易理解在某个时刻下个体暴露于不同危险因素状态下的发病风险；其他创新性的方法也各有其参数估计的方法（如加速失效时间模型的time-ratio）[12]。随着大数据时代的到来，传统的统计分析方法由于其自身理论基础（如假设检验和抽样估计）受到了使用的限制[14]。因此许多研究者开发了基于机器学习模型的生存统计分析方法，以期利用更先进算法的优势分析生存预后问题。本文在检索PubMed和CNKI数据库的基础上，重点汇总和概述了目前已发表文献使用较为普遍的、适用于非比例风险生存资料的机器学习模型（表2）。
表2
基于机器学习的非比例风险生存分析方法
2.1 一般的机器学习模型
2.1.1 随机生存森林（random survival forests，RSF）
RSF是由Ishwaran等[15]于2008年提出的一种将随机森林和传统生存分析方法相结合的树集成学习方法，是随机森林基于右删失生存数据的扩展。该模型采用自助法的采样方式随机抽取样本形成多个二元独立生存树，进而将所有的树集合形成RSF。RSF计算每课树的累计风险函数和生存函数，并对其进行集和得到整个森林的预测值，因此可不受等比例风险假设的限制[16]。RSF算法为机器学习模型在生存资料分析中应用最广的模型之一。李淼等[17]基于RSF模型建立了肺癌患者死亡风险的预测模型，并揭示了RSF模型不仅可用于右删失生存数据的分析，还能发现重要的影响因素及变量之间的交互作用。
2.1.2 极端随机生存树（extremely randomized survival trees，ERST）
ERST模型是RSF模型的一个变种，两者生存树的节点分裂标准均为log-rank分数。但是与RSF相比，ERST的随机性在计算分割的方式上更进一步。与RSF一样，使用的是候选特征的随机子集，但不是寻找最具区分力的特征值，而是随机抽取每个候选特征值，并从这些随机生成的特征值中挑选最佳的阈值来划分生存树。一般来说，ERST模型的方差相对于RSF进一步减少，但是偏差相对于RSF进一步增大。在某些时候，ERST模型的泛化能力比RSF更好。
2.1.3 生存支持向量机（survival support vector machines，SSVM）
在解决小到中等规模样本、非线性及高维模式识别中表现出许多特有优势的支持向量机模型，同样能和生存资料结合，形成了SSVM模型[18]。SSVM模型丰富的核函数选择为它带来了更大的灵活性，可有效处理生存问题中的非线性、非可加性和交互作用等问题[19]。由于支持向量机模型并不要求估计参数的潜在分布形式，因此可以解决生存资料中的非比例风险假设问题。目前使用支持向量机解决生存问题主要包括三种方法：① 回归方法：回归方法基于支持向量回归的思想，旨在找到一个函数，该函数使用协变量xi估计作为连续结果值的生存时间yi[20]；② 排序方法：排序方法将基于支持向量机的生存分析视为具有有序目标变量的分类问题，它的目的是预测个体之间的风险等级，而不是估计生存时间[20]；③ 两者都有。Van Belle等[19]认为当需处理高维度的生存资料时，建议使用基于回归的SSVM模型。
2.1.4 梯度提升生存分析（gradient boosting survival，GBS）
GBS模型是基于Boosting的用于右删失生存资料的回归树集成模型。尽管部分开发的GBS模型仍需遵循等比例风险这一假设，但是已经呈现出适用于高维数据、可以处理协变量间复杂非线性关联的优势[1,21]。Hothorn等[22]将基于加速失效时间模型的逆概率删失加权方法和梯度提升算法相结合，开发了GBS模型以解决等比例风险假设的问题。Chen等[23]基于GBS模型构造了性能良好的阴茎鳞状细胞癌患者的死亡风险预测模型。
2.1.5 多任务逻辑回归（multi-task logistic regression，MTLR）
不像Cox回归模型针对风险函数建模，MTLR通过联合多个局部的逻辑回归模型直接对生存函数建模，这使得它可以处理右删失和存在时间依赖效应的协变量（即不满足等比例风险假定）的生存资料[24]：MTLR建立了一系列不同生存时间点的逻辑回归模型，并通过包含特定于不同时间点的参数函数将这些逻辑回归模型联合起来，从而捕捉时间依赖的协变量的效应。Gu等[25]发现在预测癌症患者生存风险方面，MTLR的性能优于Cox回归模型。
2.2 深度学习模型
2.2.1 Cox-Time
Cox-time是一种基于神经网络的相对风险模型，它规避了传统Cox回归模型的比例风险假设限制。Cox-time将生存时间（t）视为常规协变量，以模拟时间变量与其他协变量之间的相互作用g（t，X），这类似于时依协变量Cox模型采用的方法，其中协变量x的非比例风险效应可以通过包括时间相关协变量［如x*t和x*log（t）］来建模：
 |
Cox-time和Cox回归模型均为基于连续时间变量的模型，因此具有无需进行生存时间离散化的优势（表1）[26]。Yang等[27]和Adeoye等[28]利用Cox-time算法开发了结直肠癌或口腔癌患者的死亡风险预测模型。
2.2.2 DeepHit
DeepHit是一个使用深度神经网络直接学习生存时间分布的离散时间模型。它的网络架构是由单个共享子网和多个特定子网组成，其损失函数同时利用了生存时间和相对风险。DeepHit模型的独特优势是可以处理竞争风险事件：在医学方面，竞争风险极为普遍，例如当研究的结局是心血管疾病死亡时，因其他原因导致的死亡会阻碍该结局的发生，这两种死亡结局互为竞争风险事件。此时若依旧运用单死亡结局终点的分析方法，将会由于竞争风险事件的存在而导致对终点事件概率的估计产生偏差。目前的竞争风险模型（如原因别风险模型和Fine-Gray模型）均基于传统统计学方法，因此基于深度学习的DeepHit模型将有更大的应用价值。
2.2.3 N-MTLR（neural multi-task logistic regression）
N-MTLR是一种基于神经网络的模型，它在不同的时间间隔上构建不同的神经网络来估计每个时间间隔内结局事件发生的概率[29]。N-MTLR与MTLR不同的是引入了深度学习架构，从而解决了MTLR无法对生存资料中非线性关联建模的问题。Fotso等[29]发现，当分析存在复杂非线性关联的生存数据时，N-MTLR的性能显著高于Cox回归模型和NTLR模型。Kiessling等[30]使用N-MTLR来估计选择性腹主动脉瘤修复后的死亡风险：该研究发现N-MTLR比Kaplan-Meier方法显示出更精确的估计，并且能够发现Kaplan-Meier无法检测到的患者亚组之间的生存差异。
3 案例研究
3.1 数据获得与清洗
案例研究数据来自美国重症监护医学数据库（Medical Information Mart for Intensive Care-Ⅳ，MIMIC-IV）。MIMIC-IV包含了2008年至2019年间美国三级学术医疗中心住院患者的临床电子病历的数据信息。我们使用SQL、PostgreSQL和Navicat Premium软件提取了2 565名中风患者的数据（图1），包括临床特征、实验室测量、合并症、生命体征和疾病严重程度评估等68个变量（表3）。我们分别对于连续型缺失变量和分类型缺失变量进行了中位数填充和众数填充，此外分别对分类变量和连续型变量进行了哑变量化和标准化。本文机器学习算法使用Python（version 3.8）编程语言实现。其中，两种集成学习模型使用了sksurv软件包（version 0.17.2），而两种深度学习模型采用了pycox（version 0.2.3）、torch （version 1.8.0+cpu）和torchtuples（version 0.2.2）软件包，其他的统计分析使用SAS（version 9.4）软件实现。
图1
研究对象纳入流程图
表3
部分变量定义及其缺失比例
3.2 非比例风险检验
针对分类型变量和连续型变量均采用了两种方法检验等比例风险假设（结果见表4、图2和图3）：表4中Schoenfeld残差和生存时间线性相关法检验显示两个分类型变量（DIAGNOSIS_TYPE、marital_status）和两个连续型变量（max_sodium、min_glucose）均不满足等比例风假定（P<0.001）；图2 log（−log（生存概率））与log（生存时间）图形显示不同的卒中诊断类型曲线交叉，不同婚姻状态的生存曲线同样显示交叉，这表明它们的等比例风险假设不成立；图3 Martingale残差图中sodium_max的实际路径在模拟路径之外，glucose_min的schoenfeld残差图在生存时间少于15天内的线性斜率不等于0，这表明不满足等比例风险假设。以上分析结果表明案例卒中生存资料不满足等比例风险假设，因此我们采用适用于非比例风险假设的机器学习和深度学习生存分析算法进行脑卒中死亡风险预测。
表4
基于Schoenfeld残差和生存时间线性相关法检验比例风险假设
图2
分类型变量比例风险假定检验
a：“DIAGNOSIS_TYPE”变量log-log作图法；b：“marital_status”变量生存曲线作图法。
图3
连续型变量比例风险假定检验
a：“sodium_max”变量Martingale残差图法；b：“glucose_min”变量Schoenfeld残差图法。
3.3 模型构建与评价
对于处理后的数据分别使用两种集成学习模型（随机生存森林和极端生存树）和两种深度学习模型（Cox-time和DeepHit）构建脑卒中死亡风险预测模型。模型的训练采用留出法（训练集∶测试集=7∶3），即在训练集上训练模型，在测试集上评价模型的性能。本研究使用Harrell’s 一致性指数（concordance index，C-index）对比不同模型的区分度[31]；使用积分布里尔评分（integrated brier score，IBS）对比不同模型的校准度[32]。本研究没有对各模型的超参数进行调整，即使用它们的默认超参数。除此之外，针对表现最好的随机生存森林模型采用排列重要性（permutation importance）的方法识别影响脑卒中患者死亡风险的重要特征[33]。
表5显示集成机器学习模型的区分度性能和校准度性能均更好：其中随机生存森林的区分度性能最高（C-index=0.773），且其校准度也更好（IBS=0.151）。基于神经网络的算法Cox-time和DeepHit表现相对较差，这可能是由于本案例研究数据规模相对较小且没有调节深度学习生存模型超参数。基于排列重要性算法得出影响脑卒中ICU住院患者死亡风险重要的特征包括年龄和日最大体温等（图4），这也与既往研究的发现一致[34, 35]。
表5
脑卒中患者死亡风险预测的4种机器学习生存算法的性能对比
图4
基于排列重要性算法获得的影响脑卒中患者死亡风险的前10个重要变量
4 讨论
Cox比例风险模型是生存分析中最常用的模型之一，它基于等比例风险假定，可以有效地分析生存数据并得出相关结论。如果生存数据不符合等比例风险假定，那么使用Cox比例风险模型进行分析就是不合适的，研究人员需要考虑其他的生存分析方法。但是传统的非比例风险生存分析方法存在无法对非线性关联建模和不适用于高维数据的固有缺陷[9]。鉴于人工智能在处理医学大数据方面的巨大作用，许多研究者将机器学习应用到生存分析领域并开发出了性能优良算法。因此本文系统概述了适用于非比例风险生存资料的机器学习生存模型并比较了各种模型的优缺点。并以脑卒中患者死亡风险研究作为研究实例比较了多种机器学习模型的预测性能，从而强调了这些机器算法在生存分析研究中应有更多的应用机会。
生存资料的NPH检验问题近些年得到越来越多的关注，主要原因在于统计学家发现肿瘤免疫靶向药物的临床试验的生存资料大多并不满足等比例风险假定（呈现延迟治疗效应）[36]，因此早在2018年就于杜克大学举办的“存在非比例危害的肿瘤学临床试验”会议中提出要开发新的统计方法[37]。但是长久以来开发的传统的非比例风险生存分析方法可能只适合于临床试验的统计分析，主要在于临床试验的数据规模和变量数目都小得多。而在如今的精准医学时代，数据的规模、维度和复杂程度都呈现急剧增加，尤其是一些组学数据还呈现出特征维度远远高于数据量大小的特点，这给既往的非比例生存分析方法带来了巨大的挑战[6]。
Huang等[9]首次开展了将机器学习的方法应用于生存预后分析方面文献的综述，这篇研究对各种机器学习生存分析模型给出了详细的介绍，但是依旧没有关注到这些模型是否适用于非比例风险生存资料。而本研究弥补了这一疏漏之处—重点汇总了适用于非比例风险生存资料的机器学习模型。除此之外，本研究还基于真实世界数据对等比例风险假定的检验方法和机器学习生存模型做了实例探索，以期为研究人员开展类似研究提供更多的参考价值。本文案例研究的目的并不是澄清最优的机器学习生存分析模型，因此并未对各个模型的超参数进行调整。相反，为生存分析选择最合适的机器学习模型应该基于自身的研究问题、数据集规模、特征维度或数据集的平衡程度等方面。
尽管基于机器学习的模型因其具备更强的处理高维数据能力和数据灵活性等优势，在预测复杂生存数据方面可能表现更佳，但相较于传统的非比例风险生存分析方法，我们认为它们在以下几个方面存在明显的局限性：1）模型的可解释性：由于传统的非比例风险生存分析方法主要在于因果关系的解释，而基于机器学习的非比例风险生存分析方法在于结局预测，因此它们的结果可解释性远不如HRs易于理解和接受[38]。现有的机器学习可解释性方法也较局限[39]，其解释结果依然不够清晰（如图4），将来的机器学习生存算法应在提高临床决策可解释性方面得到进一步发展。2）更大的样本量需求[40, 41]：在将基于机器学习的分析方法应用于复杂的生存资料时尤为显著，主要由以下几个内外部因素引起：① 机器学习模型的复杂性增加（例如，深度学习模型具有大量参数和非线性结构），因此为了确保模型能够有效地泛化和保持稳定性，通常需要更大的样本量；② 复杂的生存资料涉及更高维度的特征空间，因此需要更多的样本以充分覆盖特征空间，以避免过度拟合或欠拟合；③ 机器学习模型通常需要更多的数据进行有效的训练和调优，这包括捕获数据中的模式和规律，以及进行参数调整以提高性能；④ 机器学习模型对数据中的噪声和复杂性更为敏感，更大的样本量能够提供更多信息，帮助模型准确区分真实模式和噪声，从而提高鲁棒性和泛化能力。3）过拟合风险高[42]：由于机器学习模型具有足够的灵活性和参数数量，可以几乎完美地拟合训练数据，甚至是其中的噪声，这种过度拟合会导致模型在训练数据上表现良好，但在新数据上表现不佳。相比之下，传统的生存分析模型在一定程度上更加稳定，过拟合的风险较低。关于机器学习的模型更强调采取如交叉验证、正则化等方式以减少过拟合的风险[43]。除了上述局限性之外，基于机器学习的生存分析模型也存在参数调优困难[44]、计算资源需求高[45]和稳定性差[46]等缺陷。
既包含结局事件又包含生存时间的数据被称为生存资料，又称为时间-事件结局数据。它因同时具有定量和定性的双重性质，可提供比单纯结局事件更加丰富的数据信息。由于生存资料的这种特殊性，针对它的独特统计分析方法被称为生存分析。生存分析在帮助了解疾病的自然史、评估治疗方法的效果、预测患者的生存期，以及指导临床决策等方面有着重要的意义[1]。
目前的生存分析方法大致包括三个方面：① 生存过程的描述（如Kaplan-Meier曲线）；② 生存过程的比较（如log-rank检验）；③ 影响因素分析及生存预测（如Cox比例风险回归模型，简称Cox回归）[2]。log-rank检验和Cox回归模型是生存资料统计分析中最常用的分析策略，但是它们应用有着重要的前提条件—等比例风险假设，而这一要求在既往应用过程中常被忽视：严若华等[3]发现在PubMed检出的使用Cox回归模型的文章中，仅有0.34%（413/121 342）提及了比例风险假定。美国心脏协会发布的心血管医学统计报告中明确建议需要测试并报告Cox回归模型是否违反了等比例危险假定[4]。但是在实际数据中由于大量噪声，很难满足这个假设条件。鉴于此，一些针对不满足等比例风险假定生存资料的统计分析方法已经提出：① 基于log-rank检验或者Kaplan-Meier检验用于生存过程比较的衍生方法[5]；② 基于传统的Cox回归模型用于生存预后影响因素分析的改良方法或其他创新性方法[2]（表1）。然而，这些针对非比例风险生存数据提出的改良分析方法，从根本上说是为小规模、低维度的临床研究数据服务的。尽管它们大多有着可解释性强和易于实施等优点，但是面对大样本数据、高维特征数据或需要探索更加复杂关联形式的分析要求时，它们使用的范围就相对有限[6]。
表1
传统的非比例风险生存分析方法
随着精准医学时代的到来，可以获得的医学数据的规模和特征的维度急剧增加，因此对分析方法提出了更高的要求[7]。作为人工智能重要分支之一的机器学习模型（包括人工神经网络）正越来越渗透于医药健康研究中[8]。针对生存分析的问题，许多研究者也开发了相关机器学习模型来分析生存预后问题[9]。机器学习模型对于数据分布的限制极少，且有更大的假设空间和拟合效果，另外这些机器学习模型在处理高维特征生存大数据或者存在竞争风险生存数据方面也有着独特的优势[6]。许多机器学习生存模型是对Cox回归模型的改进（如DeepSurv[10]和Modified-DeepSurv[11]），尽管它们在解决Cox回归模型中的非线性和交互作用等挑战方面表现出了有效性，但对于存在非比例风险的生存数据的适用性却受到限制。
本文在解读NPH等相关概念之后，重点概述了可以处理不满足等比例风险假定生存资料的机器学习模型，并进行了基于机器学习模型的脑卒中患者死亡风险的案例研究，以期推动人工智能在非比例风险生存资料分析中的应用。
1 相关概念
1.1 非比例风险
Cox回归模型使用广泛的主要原因在于其没有对生存时间和基础风险函数的分布形式做任何的要求，且可以方便的求得风险比（hazard ratios，HRs）以衡量协变量的影响效应。但是Cox回归模型需要满足等比例风险的假设条件，在该条件下求得HRs才会在不同的生存时间保持一致。Cox回归模型是一个半参数模型，在等比例假设下，HRs计算公式如下：
|
等比例风险假定即HRs与时间无关，其不随时间变化而变化。比如在针对某种药物开展的随机对照试验中，试验组和对照组的HRs在任何的随访时间都是相等的，即药物的治疗效应对生存率的影响不随时间的改变而改变。而当HRs随时间变化，则为NPH。NPH下计算的HRs仅代表随访期间估计的HRs的平均水平，但是这种平均估计是不准确的，尤其对于存在交叉治疗效应的生存数据[12]。
1.2 非比例风险检验
NPH检验是生存资料统计分析过程中的一个重要步骤，既往研究已经揭示了不同的方法以探索是否满足等比例风险假定，本研究在既往研究的基础上建议依据变量类型的不同而采取不同的方法进行检验[2,3,12,13]：① 连续型变量：Schoenfeld残差图法、Martingale残差图法、Schoenfeld残差和时间线性相关分析法、时依协变量假设检验法；② 分类型变量：Kaplan-Meier生存曲线图、log（-log（生存概率））与log（生存时间）作图法（log-log作图法）、Schoenfeld残差图法、Schoenfeld残差和时间线性相关分析法。此外，建议针对每个协变量至少使用两种以上的策略以判断其是否满足等比例风险假设。
2 基于机器学习算法的非比例风险模型
传统的生存分析统计方法在数据规模较小、特征维度较低及算力有限的情况下尤为适用，且由于它们能提供协变量的参数估计而可解释性更强[14]。以针对非比例风险生存资料衍生的Cox回归模型为例（表1），它们依旧可以提供HRs，从而很容易理解在某个时刻下个体暴露于不同危险因素状态下的发病风险；其他创新性的方法也各有其参数估计的方法（如加速失效时间模型的time-ratio）[12]。随着大数据时代的到来，传统的统计分析方法由于其自身理论基础（如假设检验和抽样估计）受到了使用的限制[14]。因此许多研究者开发了基于机器学习模型的生存统计分析方法，以期利用更先进算法的优势分析生存预后问题。本文在检索PubMed和CNKI数据库的基础上，重点汇总和概述了目前已发表文献使用较为普遍的、适用于非比例风险生存资料的机器学习模型（表2）。
表2
基于机器学习的非比例风险生存分析方法
2.1 一般的机器学习模型
2.1.1 随机生存森林（random survival forests，RSF）
RSF是由Ishwaran等[15]于2008年提出的一种将随机森林和传统生存分析方法相结合的树集成学习方法，是随机森林基于右删失生存数据的扩展。该模型采用自助法的采样方式随机抽取样本形成多个二元独立生存树，进而将所有的树集合形成RSF。RSF计算每课树的累计风险函数和生存函数，并对其进行集和得到整个森林的预测值，因此可不受等比例风险假设的限制[16]。RSF算法为机器学习模型在生存资料分析中应用最广的模型之一。李淼等[17]基于RSF模型建立了肺癌患者死亡风险的预测模型，并揭示了RSF模型不仅可用于右删失生存数据的分析，还能发现重要的影响因素及变量之间的交互作用。
2.1.2 极端随机生存树（extremely randomized survival trees，ERST）
ERST模型是RSF模型的一个变种，两者生存树的节点分裂标准均为log-rank分数。但是与RSF相比，ERST的随机性在计算分割的方式上更进一步。与RSF一样，使用的是候选特征的随机子集，但不是寻找最具区分力的特征值，而是随机抽取每个候选特征值，并从这些随机生成的特征值中挑选最佳的阈值来划分生存树。一般来说，ERST模型的方差相对于RSF进一步减少，但是偏差相对于RSF进一步增大。在某些时候，ERST模型的泛化能力比RSF更好。
2.1.3 生存支持向量机（survival support vector machines，SSVM）
在解决小到中等规模样本、非线性及高维模式识别中表现出许多特有优势的支持向量机模型，同样能和生存资料结合，形成了SSVM模型[18]。SSVM模型丰富的核函数选择为它带来了更大的灵活性，可有效处理生存问题中的非线性、非可加性和交互作用等问题[19]。由于支持向量机模型并不要求估计参数的潜在分布形式，因此可以解决生存资料中的非比例风险假设问题。目前使用支持向量机解决生存问题主要包括三种方法：① 回归方法：回归方法基于支持向量回归的思想，旨在找到一个函数，该函数使用协变量xi估计作为连续结果值的生存时间yi[20]；② 排序方法：排序方法将基于支持向量机的生存分析视为具有有序目标变量的分类问题，它的目的是预测个体之间的风险等级，而不是估计生存时间[20]；③ 两者都有。Van Belle等[19]认为当需处理高维度的生存资料时，建议使用基于回归的SSVM模型。
2.1.4 梯度提升生存分析（gradient boosting survival，GBS）
GBS模型是基于Boosting的用于右删失生存资料的回归树集成模型。尽管部分开发的GBS模型仍需遵循等比例风险这一假设，但是已经呈现出适用于高维数据、可以处理协变量间复杂非线性关联的优势[1,21]。Hothorn等[22]将基于加速失效时间模型的逆概率删失加权方法和梯度提升算法相结合，开发了GBS模型以解决等比例风险假设的问题。Chen等[23]基于GBS模型构造了性能良好的阴茎鳞状细胞癌患者的死亡风险预测模型。
2.1.5 多任务逻辑回归（multi-task logistic regression，MTLR）
不像Cox回归模型针对风险函数建模，MTLR通过联合多个局部的逻辑回归模型直接对生存函数建模，这使得它可以处理右删失和存在时间依赖效应的协变量（即不满足等比例风险假定）的生存资料[24]：MTLR建立了一系列不同生存时间点的逻辑回归模型，并通过包含特定于不同时间点的参数函数将这些逻辑回归模型联合起来，从而捕捉时间依赖的协变量的效应。Gu等[25]发现在预测癌症患者生存风险方面，MTLR的性能优于Cox回归模型。
2.2 深度学习模型
2.2.1 Cox-Time
Cox-time是一种基于神经网络的相对风险模型，它规避了传统Cox回归模型的比例风险假设限制。Cox-time将生存时间（t）视为常规协变量，以模拟时间变量与其他协变量之间的相互作用g（t，X），这类似于时依协变量Cox模型采用的方法，其中协变量x的非比例风险效应可以通过包括时间相关协变量［如x*t和x*log（t）］来建模：
|
Cox-time和Cox回归模型均为基于连续时间变量的模型，因此具有无需进行生存时间离散化的优势（表1）[26]。Yang等[27]和Adeoye等[28]利用Cox-time算法开发了结直肠癌或口腔癌患者的死亡风险预测模型。
2.2.2 DeepHit
DeepHit是一个使用深度神经网络直接学习生存时间分布的离散时间模型。它的网络架构是由单个共享子网和多个特定子网组成，其损失函数同时利用了生存时间和相对风险。DeepHit模型的独特优势是可以处理竞争风险事件：在医学方面，竞争风险极为普遍，例如当研究的结局是心血管疾病死亡时，因其他原因导致的死亡会阻碍该结局的发生，这两种死亡结局互为竞争风险事件。此时若依旧运用单死亡结局终点的分析方法，将会由于竞争风险事件的存在而导致对终点事件概率的估计产生偏差。目前的竞争风险模型（如原因别风险模型和Fine-Gray模型）均基于传统统计学方法，因此基于深度学习的DeepHit模型将有更大的应用价值。
2.2.3 N-MTLR（neural multi-task logistic regression）
N-MTLR是一种基于神经网络的模型，它在不同的时间间隔上构建不同的神经网络来估计每个时间间隔内结局事件发生的概率[29]。N-MTLR与MTLR不同的是引入了深度学习架构，从而解决了MTLR无法对生存资料中非线性关联建模的问题。Fotso等[29]发现，当分析存在复杂非线性关联的生存数据时，N-MTLR的性能显著高于Cox回归模型和NTLR模型。Kiessling等[30]使用N-MTLR来估计选择性腹主动脉瘤修复后的死亡风险：该研究发现N-MTLR比Kaplan-Meier方法显示出更精确的估计，并且能够发现Kaplan-Meier无法检测到的患者亚组之间的生存差异。
3 案例研究
3.1 数据获得与清洗
案例研究数据来自美国重症监护医学数据库（Medical Information Mart for Intensive Care-Ⅳ，MIMIC-IV）。MIMIC-IV包含了2008年至2019年间美国三级学术医疗中心住院患者的临床电子病历的数据信息。我们使用SQL、PostgreSQL和Navicat Premium软件提取了2 565名中风患者的数据（图1），包括临床特征、实验室测量、合并症、生命体征和疾病严重程度评估等68个变量（表3）。我们分别对于连续型缺失变量和分类型缺失变量进行了中位数填充和众数填充，此外分别对分类变量和连续型变量进行了哑变量化和标准化。本文机器学习算法使用Python（version 3.8）编程语言实现。其中，两种集成学习模型使用了sksurv软件包（version 0.17.2），而两种深度学习模型采用了pycox（version 0.2.3）、torch （version 1.8.0+cpu）和torchtuples（version 0.2.2）软件包，其他的统计分析使用SAS（version 9.4）软件实现。
图1
研究对象纳入流程图
表3
部分变量定义及其缺失比例
3.2 非比例风险检验
针对分类型变量和连续型变量均采用了两种方法检验等比例风险假设（结果见表4、图2和图3）：表4中Schoenfeld残差和生存时间线性相关法检验显示两个分类型变量（DIAGNOSIS_TYPE、marital_status）和两个连续型变量（max_sodium、min_glucose）均不满足等比例风假定（P<0.001）；图2 log（−log（生存概率））与log（生存时间）图形显示不同的卒中诊断类型曲线交叉，不同婚姻状态的生存曲线同样显示交叉，这表明它们的等比例风险假设不成立；图3 Martingale残差图中sodium_max的实际路径在模拟路径之外，glucose_min的schoenfeld残差图在生存时间少于15天内的线性斜率不等于0，这表明不满足等比例风险假设。以上分析结果表明案例卒中生存资料不满足等比例风险假设，因此我们采用适用于非比例风险假设的机器学习和深度学习生存分析算法进行脑卒中死亡风险预测。
表4
基于Schoenfeld残差和生存时间线性相关法检验比例风险假设
图2
分类型变量比例风险假定检验
a：“DIAGNOSIS_TYPE”变量log-log作图法；b：“marital_status”变量生存曲线作图法。
图3
连续型变量比例风险假定检验
a：“sodium_max”变量Martingale残差图法；b：“glucose_min”变量Schoenfeld残差图法。
3.3 模型构建与评价
对于处理后的数据分别使用两种集成学习模型（随机生存森林和极端生存树）和两种深度学习模型（Cox-time和DeepHit）构建脑卒中死亡风险预测模型。模型的训练采用留出法（训练集∶测试集=7∶3），即在训练集上训练模型，在测试集上评价模型的性能。本研究使用Harrell’s 一致性指数（concordance index，C-index）对比不同模型的区分度[31]；使用积分布里尔评分（integrated brier score，IBS）对比不同模型的校准度[32]。本研究没有对各模型的超参数进行调整，即使用它们的默认超参数。除此之外，针对表现最好的随机生存森林模型采用排列重要性（permutation importance）的方法识别影响脑卒中患者死亡风险的重要特征[33]。
表5显示集成机器学习模型的区分度性能和校准度性能均更好：其中随机生存森林的区分度性能最高（C-index=0.773），且其校准度也更好（IBS=0.151）。基于神经网络的算法Cox-time和DeepHit表现相对较差，这可能是由于本案例研究数据规模相对较小且没有调节深度学习生存模型超参数。基于排列重要性算法得出影响脑卒中ICU住院患者死亡风险重要的特征包括年龄和日最大体温等（图4），这也与既往研究的发现一致[34, 35]。
表5
脑卒中患者死亡风险预测的4种机器学习生存算法的性能对比
图4
基于排列重要性算法获得的影响脑卒中患者死亡风险的前10个重要变量
4 讨论
Cox比例风险模型是生存分析中最常用的模型之一，它基于等比例风险假定，可以有效地分析生存数据并得出相关结论。如果生存数据不符合等比例风险假定，那么使用Cox比例风险模型进行分析就是不合适的，研究人员需要考虑其他的生存分析方法。但是传统的非比例风险生存分析方法存在无法对非线性关联建模和不适用于高维数据的固有缺陷[9]。鉴于人工智能在处理医学大数据方面的巨大作用，许多研究者将机器学习应用到生存分析领域并开发出了性能优良算法。因此本文系统概述了适用于非比例风险生存资料的机器学习生存模型并比较了各种模型的优缺点。并以脑卒中患者死亡风险研究作为研究实例比较了多种机器学习模型的预测性能，从而强调了这些机器算法在生存分析研究中应有更多的应用机会。
生存资料的NPH检验问题近些年得到越来越多的关注，主要原因在于统计学家发现肿瘤免疫靶向药物的临床试验的生存资料大多并不满足等比例风险假定（呈现延迟治疗效应）[36]，因此早在2018年就于杜克大学举办的“存在非比例危害的肿瘤学临床试验”会议中提出要开发新的统计方法[37]。但是长久以来开发的传统的非比例风险生存分析方法可能只适合于临床试验的统计分析，主要在于临床试验的数据规模和变量数目都小得多。而在如今的精准医学时代，数据的规模、维度和复杂程度都呈现急剧增加，尤其是一些组学数据还呈现出特征维度远远高于数据量大小的特点，这给既往的非比例生存分析方法带来了巨大的挑战[6]。
Huang等[9]首次开展了将机器学习的方法应用于生存预后分析方面文献的综述，这篇研究对各种机器学习生存分析模型给出了详细的介绍，但是依旧没有关注到这些模型是否适用于非比例风险生存资料。而本研究弥补了这一疏漏之处—重点汇总了适用于非比例风险生存资料的机器学习模型。除此之外，本研究还基于真实世界数据对等比例风险假定的检验方法和机器学习生存模型做了实例探索，以期为研究人员开展类似研究提供更多的参考价值。本文案例研究的目的并不是澄清最优的机器学习生存分析模型，因此并未对各个模型的超参数进行调整。相反，为生存分析选择最合适的机器学习模型应该基于自身的研究问题、数据集规模、特征维度或数据集的平衡程度等方面。
尽管基于机器学习的模型因其具备更强的处理高维数据能力和数据灵活性等优势，在预测复杂生存数据方面可能表现更佳，但相较于传统的非比例风险生存分析方法，我们认为它们在以下几个方面存在明显的局限性：1）模型的可解释性：由于传统的非比例风险生存分析方法主要在于因果关系的解释，而基于机器学习的非比例风险生存分析方法在于结局预测，因此它们的结果可解释性远不如HRs易于理解和接受[38]。现有的机器学习可解释性方法也较局限[39]，其解释结果依然不够清晰（如图4），将来的机器学习生存算法应在提高临床决策可解释性方面得到进一步发展。2）更大的样本量需求[40, 41]：在将基于机器学习的分析方法应用于复杂的生存资料时尤为显著，主要由以下几个内外部因素引起：① 机器学习模型的复杂性增加（例如，深度学习模型具有大量参数和非线性结构），因此为了确保模型能够有效地泛化和保持稳定性，通常需要更大的样本量；② 复杂的生存资料涉及更高维度的特征空间，因此需要更多的样本以充分覆盖特征空间，以避免过度拟合或欠拟合；③ 机器学习模型通常需要更多的数据进行有效的训练和调优，这包括捕获数据中的模式和规律，以及进行参数调整以提高性能；④ 机器学习模型对数据中的噪声和复杂性更为敏感，更大的样本量能够提供更多信息，帮助模型准确区分真实模式和噪声，从而提高鲁棒性和泛化能力。3）过拟合风险高[42]：由于机器学习模型具有足够的灵活性和参数数量，可以几乎完美地拟合训练数据，甚至是其中的噪声，这种过度拟合会导致模型在训练数据上表现良好，但在新数据上表现不佳。相比之下，传统的生存分析模型在一定程度上更加稳定，过拟合的风险较低。关于机器学习的模型更强调采取如交叉验证、正则化等方式以减少过拟合的风险[43]。除了上述局限性之外，基于机器学习的生存分析模型也存在参数调优困难[44]、计算资源需求高[45]和稳定性差[46]等缺陷。
