乳腺癌已成为全球最常见的女性癌病，是最致命的女性癌症之一^[1]，2020年约有68.5万女性死于患乳腺癌^[2]。研究发现部分乳腺癌患者在确诊后十年期间仍存活于复发可能性的阴影之下^[3]，因此针对乳腺癌患者的预后管理显得尤为重要，但其适用性和有效性又受限于临床医生颇具主观性的专业经验。近年来，生存分析^[4]技术被用以辅助像乳腺癌这种发病不明显、病程较长的疾病制定预后措施，研究人员将预后和预测特征即患者病理协变量独立或整合以指导患者的治疗^[5]。鉴于个体化医学的重要性，研究人员结合当下个体化医学研究中流行的用于发现并识别复杂数据之间模式和关系的机器学习技术，对乳腺癌等长病程数据集进行个体生存分析。此类方法基于患者间的病理协变量表达差异实现个性化的生存风险估计，进而绘制患者个体生存曲线，以辅助实现精准医疗。

近年来，个体生存分析研究主要分为基于Cox比例风险假设的生存分析和无参数型离散生存分析。Cox比例风险模型假设独立协变量对生存时间的影响不随时间而变化，以评估多个协变量对生存的影响^[6]。Katzman等 ^[7]基于Cox对数比例风险假设搭建一个深度前馈神经网络DeepSurv，用于预测患者协变量对风险率的影响。Lee等 ^[8]使用一个共享子网络以处理并预测时变协变量，同时设计多个特定原因子网络捕获特定原因风险同协变量之间的关系，构成Dynamic-DeepHit模型用于动态风险预测。Tang等^[9]对删失数据进行最大似然估计并结合常微分方程进行生存分析建模。上述单任务回归模型遵循对数风险假设，导致乳腺癌患者间的风险函数永不相交，虽然可描绘出患者病情在总生存期中一个大致的走势，但难以实现阶段性且个性化的生存风险描述。Yu等^[10]提出多任务线性回归模型(Multi-task Logistic Regression, MTLR)，其采用响应转换的方式绕过比例风险假设以实现生存函数预测，但回归器核心为线性变换，难以捕捉患者生存数据中的非线性关系。Fotso等^[11]通过往MTLR中引入多层感知器，实现灵活的建模，并命名为N-MTLR(Neural Multi-task Logistic Regression)。Hu等^[12]基于时下流行的Transformer网络架构，使用有序回归实现对时变生存概率的优化，但受限于多头注意力机制在病理协变量对风险概率贡献较均衡时难以分配权重。Wang等^[13]提出了多任务带状回归模型(Multi-task Banded Regression, MTBR)，带状校验矩阵的引入为响应转换提供约束，解决风险概率下降缓慢导致的生存概率分配不均的情况。以上的多任务回归模型采用线性变换或多层感知器进行响应转换，使得协变量独立于多个序列之间，难以提取出病理协变量对不同患者生存时间之间影响的异同关系，导致模型一致性指数下降。

为此，本文提出了一种新的多任务回归模型RNN-MTBR(Multi-task Banded Regression via Recurrent Neural Network) 用于乳腺癌个体生存分析。首先，针对以往的多任务模型的回归器无法使序列之间产生前后联系的状况，使用循环神经网络(Recurrent Neural Network, RNN) ^[14]作为多任务回归模型的核心，RNN可通过维护网络过去的记忆，让模型在提取病理协变量同生存时间的关系时能利用到序列的上下文信息。其次，将带状校验矩阵改进并探究其对风险分布的作用效果。最后，在乳腺癌真实数据集上实现了对个体生存分析区分度的优化，成功验证了模型的有效性。

1 乳腺癌生存分析

给定一个含有M名患者的乳腺癌生存数据集，设$\forall i \in \left[ {1,M} \right]\text{，}({X_i},{\delta _i},{t_i})$为一个数据点，其中包含生存时间、与生存时间有关联的p维独立变量${X_i} = \left[ {x_1^i,x_2^i, \cdots ,x_p^i} \right]$和事件类型标签${\delta _i} $(0或1)。生存数据分为完全数据、删失数据。完全数据指从观察起点到发生死亡事件所经历的时间均被观察到的生存数据。删失数据则分为左删失数据、区间删失数据和右删失数据^[15]。左删失表示事件发生在观察时间开始之前。区间删失表示只能得知患者在不连续的2次随访区间中发生终点事件与否，而无法得知终点事件发生的具体时间。右删失表示在随访观察中，已知患者被观察的起始时间，而无法得知患者发生终点事件的时间，其生存时间大于观察时间。本文实验中“删失”数据只包含“右删失”数据。

对于在观察期间内能够完整观察到事件发生的患者个体实例，其生存时间记为实际事件时间$ {t}_{{\rm{o}}}^{i} $。对于在观察期间内失去对患者的观察，或者患者事件发生时间超出观察时间的实例，其生存时间记为删失时间$ {t}_{{\rm{c}}}^{i} $，$t_{\rm{o}}^i $为未删失时间。患者个体生存时间$ {t}_{i}=\mathrm{m}\mathrm{i}\mathrm{n}({t}_{{\rm{o}}}^{i},{t}_{{\rm{c}}}^{i}) $。生存数据如图1所示。患者1为完全数据，因乳腺癌死亡于观察区间之外的患者2 为右删失数据，发生了其他事件死亡的患者3 也为右删失数据。

乳腺癌个体生存分析利用患者病理独立协变量$ {X}_{i} $与生存时间$ {t}_{i} $建立一个乳腺癌风险随时间变化的函数，常见为风险函数h(t) 、风险概率分布函数f(t)或概率累计函数F(t) 。

2 多任务循环神经网络带状回归 2.1 基本原理及流程

多任务方法通过在不同时间间隔上建立回归模型，来计算乳腺癌患者死亡在每个时间间隔中发生的概率，使不同患者间的结果具有更为明显的差异性。多任务回归模型以响应变量的形式来分别代表患者在每一时间间隔的状态(生存或死亡)。假设总生存期分为N个生存子区间，则在第j个生存子区间$ {\tau }_{j} $(左开右闭)内患者的响应变量为

$ {y}_{i}=\left\{\begin{array}{l}1,{\tau }_{j}为患者死亡事件发生所在生存子区间\\ 0,其余生存子区间\end{array}\right. $

(1)

且患者发生死亡事件后的所有时间间隔响应变量值将恒为1。因而患者的生存向量为

$ {\boldsymbol{Y}}={[{y}_{1}=0, \cdots ,{y}_{j-1}=0,{y}_{j}=1,{\cdots ,y}_{N}=1]}^{\mathrm{T}} $

(2)

多任务循环神经网络带状回归模型的处理流程如图2所示。将患者的病理协变量X输入循环神经网络回归器中并以生存向量${\boldsymbol{Y}} $作为监督向量，多个回归器输出的生存响应$ {\boldsymbol{\psi}} \left({\boldsymbol{X}}\right) $再经过带状响应转换，则得到患者在每个生存子区间的风险概率，即可拟合出患者的生存曲线。

2.2 带状校验矩阵的优化

当乳腺癌患者的生存时间分为N个生存子区间，则多任务回归模型的可信度矩阵^[16]为

$ \boldsymbol{A}=\left[\begin{array}{ccc}{\boldsymbol{A}}_{{\boldsymbol{1}}}& & \\ & \begin{array}{ccc}\ddots & & \\ & {\boldsymbol{A}}_{\boldsymbol{i}}& \\ & & \ddots \end{array}& \\ & & {\boldsymbol{A}}_{\boldsymbol{N}}\end{array}\right] $

(3)

式中：${\boldsymbol{A}}_{\boldsymbol{i}}$为组分i的系数矩阵。

MTBR中将各组分系数矩阵设计为正二对角块非奇矩阵，定义为${\boldsymbol{A}}_{\boldsymbol{i}}=\left[\begin{array}{c}1\\ -1\end{array}\right]$。此种系数矩阵设计造成回归器生成每一子区间的生存响应只同上一子区间之间产生联系。为了在生成这一子区间的生存响应时对下一子区间的生存响应有所约束，本文将带状校验系数矩阵优化为只含有主对角线元素和次对角线元素的Toeplitz矩阵^[17]，故带状校验矩阵$ \boldsymbol{B} $为

${\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} {{b_0}}&{{b_1}}&0& \cdots &0\\ {{b_{ - 1}}}&{{b_0}}&{{b_1}}& \cdots &0\\ 0&{{b_{ - 1}}}&{{b_0}}& \cdots &0\\ \vdots & \vdots & \vdots & &{{b_1}}\\ 0&0&0&{{b_{ - 1}}}&{{b_0}} \end{array}} \right] $

(4)

且必须保持矩阵主对角线占优，故矩阵需满足各行和(对称矩阵中同各列和)为1，且$ {b}_{0} > {b}_{-1}+{b}_{1} $。后续实验将研究$ {b}_{-1} $、$ {b}_{1} $变化对患者生存概率分布的影响。

2.3 生存函数计算及损失函数设计

本文中乳腺癌患者的风险分布函数为

$ {\boldsymbol{f}}\left({\tau }_{s},{\boldsymbol{X}}\right) =\frac{\mathrm{e}\mathrm{x}\mathrm{p}\left({\boldsymbol{\psi}} \left({\boldsymbol{X}}\right) \varsigma {\boldsymbol{\varDelta}} {\boldsymbol{Y}}\right) }{{\boldsymbol{Z}}\left({\boldsymbol{\psi}} \left({\boldsymbol{X}}\right) {\boldsymbol{\varsigma}} \right) } $

(5)

式中：$ {\boldsymbol{\varsigma}} $为带状校验矩阵$ \boldsymbol{B} $的逆矩阵，$ {\boldsymbol{\varDelta}} $为下三角矩阵，${\boldsymbol{Z}}\left(\cdot \right)$为归一化操作。乳腺癌患者生存函数为

$ S\left({\tau }_{s},{\boldsymbol{X}}\right) =1-\sum _{k=1}^{s}\frac{\mathrm{e}\mathrm{x}\mathrm{p}\left({\boldsymbol{\psi}} \left({\boldsymbol{X}}\right) {\boldsymbol{\varsigma}} {\boldsymbol{\varDelta}} \widehat{y}\right) }{{\boldsymbol{Z}}\left({\boldsymbol{\psi}} \left({\boldsymbol{X}}\right) {\boldsymbol{\varsigma}} \right) } $

(6)

模型的损失函数由两部分组成，第1部分为累计损失函数，公式为

$ \begin{split}{L}_{1}=&\sum _{i=1}^{K}\bigg(\delta \left({k}_{i}\ne \mathrm{\varnothing }\right) \sum _{j=0}^{N}|\widehat{F}({\tau }_{j}|{\boldsymbol{X}}) -{y}_{j}|+\\ & \delta \left({k}_{i}=\mathrm{\varnothing }\right) \sum _{j=0}^{{\tau }_{{\rm{e}}}}|\widehat{F}({\tau }_{j}|{\boldsymbol{X}}) -{y}_{j}|\bigg) \end{split}$

式中：$\widehat{F}\left(\cdot \right)$为在离散生存子区间上的生存估计，K为患者总数，$ k $为样本事件情况，$ \tau $为生存子区间，$ {\tau }_{{\rm{e}}} $为事件发生或观察到删失时所在生存子区间。

第2部分为排序损失函数，公式为

$ {L}_{2}=\sum _{k=1}^{K}\sum _{i\ne j}\delta ({k}_{i}={k}_{j},{y}_{i} < {y}_{j}) \mathrm{e}\mathrm{x}\mathrm{p}(\widehat{F}({\tau }_{k}|{{\boldsymbol{X}}}_{j}) -\widehat{F}({\tau }_{k}|{{\boldsymbol{X}}}_{i}) ) $

(8)

3 实验与讨论 3.1 数据集与实验预设置

本文采用METABRIC^[18]和GBSG^[19]2种常用的乳腺癌真实数据集来评估模型的性能。METABRIC数据集包含1 980例乳腺癌患者的9项特征记录，其中约58%的患者死于患乳腺癌，中位生存时间为116个月。GBSG数据集包含2232例乳腺癌患者的记录，其中1546例为淋巴结阳性乳腺癌患者，以及686例随机临床试验的患者完整数据，删失数据约占57%。

本文采用生存分析研究中常用的一致性指数(Concordance index, C-index) ^[20]作为评估指标，C-index常被用于评价预后模型的预测精度，判断出模型的区分能力，公式为

$ \mathrm{C}{\text{-}}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}=\frac{\displaystyle\sum_{i,j}\varphi ({t}_{i} < {t}_{j}) \varphi ({\eta }_{i} < {\eta }_{j}) {\delta }_{i}}{\displaystyle\sum_{i,j}\varphi ({t}_{i} < {t}_{j}) {\delta }_{i}} $

(9)

式中：$\varphi (\cdot )$为患者对的比较，$ \eta $为生存子区间的风险值。

为了验证模型的可靠性，将数据集以6:2:2的比例分为训练集、验证集和测试集，并采用更换随机种子的方式进行10次实验，基于正态分布法计算C-index的95%置信区间。RNN-MTBR的模型采用循环神经网络作为核心，隐藏层维度为32，选用Adam优化器，Dropout概率设为0.1，学习率设为0.01。

3.2 实验结果 3.2.1 带状校验矩阵系数对患者生存概率分布的影响

带状校验矩阵系数变化会对RNN-MTBR模型性能产生影响，本文通过固定带状校验矩阵主对角线系数$ {b}_{0} $，变更次对角线系数$ {b}_{-1} $、$ {b}_{1} $的大小关系进行实验，其中${b}_{0}=\dfrac{5}{9}$，${b}_{-1}\in \left[0,\dfrac{4}{9}\right]$，${b}_{1}\in \left[0,\dfrac{4}{9}\right]$，且$ {b}_{0}+{b}_{-1}+ {b}_{1}=1 $。为了更明显地看出带状校验矩阵系数对乳腺癌患者生存概率的作用效果，实验将得到的患者相邻生存子区间的生存概率进行线性10等分，并随机抽取其中5名患者的生存曲线进行绘制。如图3所示，随着$ {b}_{-1} $的增大，患者生存风险呈现逐步前移的现象，反之，随着$ {b}_{1} $的增大，患者生存风险则出现逐步后移。而且能明显看出，在引入带状矩阵后，患者生存率的骤变情况能够得到明显的缓解。

3.2.2 带状校验矩阵优化对性能的影响

图4为各组分为正二对角块非奇矩阵的带状校验矩阵与优化为只含有主对角线元素和次对角线元素的Toeplitz矩阵生存模型C-index的对比。经比较，优化后带状校验矩阵对模型性能的提升更大。

3.2.3 生存分析模型性能对比

CoxCC^[21]、DeepSurv^[7]、DeepHit^[8]、SODEN^[9]、N-MTLR^[11]、Transformer-Based DSA^[12]和MTBR^[13]共7种生存分析模型作为RNN-MTBR的实验对比模型。表1为各生存分析模型在乳腺癌数据集测试集上的C-index值及95%置信区间。从表1可以看出，本文提出的RNN-MTBR模型效果在变更回归器核心为循环神经网络和优化了带状校验矩阵后，模型在2个数据集上的预测性能相比原先的MTBR模型均有所提升，且同其他主流的生存分析模型相比具有最好的预测性能。

图5给出了几种模型分别在GBSG数据集中拟合相同的5例乳腺癌患者的生存曲线的效果。CoxCC和DeepSurv属于比例风险模型，该种模型默认所有个体拥有相同的基线风险，且实际风险与基线风险成正比，使得计算得到的患者间风险不会出现交叉，这与现实相悖。DeepHit没有考虑对连续时间的尺度离散化，导致患者间生存曲线差异化不明显。N-MTLR将多层感知器引入多任务回归中，能提取患者生存数据中的非线性关系，但存在响应转换在阈值附近波动的现象，导致生存曲线中较多出现生存率下降幅度减小的情况。MTBR通过引入带状校验矩阵有效抑制了响应转换波动在阈值附近出现的现象，使得患者生存曲线能呈现较明显的下降趋势，但患者之间差异性的体现仍不明显。RNN-MTBR将MTBR核心替换为RNN，输入数据则得以在各时间步之间实现上下文关联，并能在不同时间步之间共享权重，并合理调整带状校验矩阵，使得患者生存概率下降幅度随时间变化显得更为合理且患者之间的生存曲线能有较为明显的差异性，从而方便个性化预后管理的施行。

4 结论

本文采用循环神经网络作为多任务回归器，优化设计了带状校验矩阵，提出了一种多任务循环神经网络带状回归模型，并应用于乳腺癌个体生存分析。在2个乳腺癌真实数据集的生存分析中，多任务循环神经网络带状回归模型的C-index较医学界常用的Cox比例风险模型提升了0.02，且有效缩小了95%置信区间。研究还通过可视化的生存曲线验证了新模型对患者间区分效果的增益性。