摘要: 本文主要分析了来自美国癌症数据库SEER的盲肠癌数据。首先通过随机生存森林模型(Random survival forest)进行盲肠癌的独立预后因素初步筛选，筛选出来的变量为：AJCC (American Joint Committee on Cancer)分期、肿瘤大小、年龄、婚姻状况、组织学分级、化疗状况、种族、放疗状况。然后通过筛选出来的变量分别建立了多因素Cox比例风险回归模型和多因素竞争风险模型。结果表明：Cox比例风险回归模型中，化疗治疗、已婚、肿瘤直径大小在1 cm以上的为盲肠癌患者生存预后的保护因素，年龄大于65岁患者、放疗治疗、AJCC分期大于I、组织学等级高于一级、婚姻状况为其它的因素为危险因素；在竞争风险模型中，化疗治疗、肿瘤直径大小在1 cm以上变量为盲肠癌患者生存预后的保护因素，年龄大于65岁患者、AJCC分期大于I、组织学等级高于I级、放疗治疗都为危险因素。在模型的比较中，竞争风险模型更胜一筹，在对于存在竞争事件的生存分析中，选择基于竞争风险构建的预测模型不仅准确度高，而且更具合理性。
Abstract:
This article mainly analyzes the cecum cancer data from the US cancer database SEER. Firstly, the independent prognostic factors of cecum cancer were preliminarily screened by the Random survival forest, and the variables screened out were: AJCC Stage, Tumor Size, Age, Marital status, Grade, Chemotherapy status, Race, and Radiotherapy status. Then, the multi-factor Cox proportional risk regression model and the multi-factor competitive risk model were established by the filtered variables. The results showed that in the Cox proportional risk regression model, chemotherapy treatment, marriage, and tumor diameter size of more than 1 cm were the protective factors for survival and prognosis of patients with cecum cancer, and patients with age greater than 65 years old, radiotherapy treatment, AJCC stage was greater than I, Grade was higher than grade I, and marital status was other factors as risk factors. In the competitive risk model, chemotherapy therapy and tumor diameter size of more than 1 cm were the protective factors for survival prognosis in patients with cecum cancer, and patients older than 65 years old, AJCC stage greater than I, Grade higher than grade I, and radiotherapy therapy were all risk factors. In the comparison of models, the competitive risk model is superior, and in the survival analysis of the existence of competitive events, the selection of a prediction model based on competitive risk is not only more accurate, but also more reasonable.
1. 引言
结直肠癌(Colorectal cancer)是全球第三大常见恶性肿瘤和第四大癌症相关死亡原因 [1]，因此是一种全球关注的癌症 [2]。根据最近的一项研究，被诊断为盲肠癌的患者的预后仍然比上升结肠癌患者更差，被诊断为盲肠癌的患者需要承担更多的疾病负担 [3]。因此，探索影响盲肠癌患者预后的风险因素将有助于临床医生制定对这些患者有利的个性化诊断和治疗方案 [4]。
本文选取来自于SEER数据库中的5240例盲肠癌患者数据进行生存分析研究，考虑到机器学习等大数据分析手段目前充分应用到统计学学科中，它能够适应各类复杂条件，本文首先引用了随机生存森林模型(Random survival forests)对影响盲肠癌生存预后因素进行初步筛选 [5]，其次基于筛选后的因素分别使用Cox比例风险模型 [6] 和竞争风险模型 [7] 对盲肠癌患者的生存预后进行预测研究。
2. 数据说明与预处理
本文纳入分析的变量为9个，具体的变量说明如表1所示。因变量为生存时间和生存结局，特别这里的生存结局不单只有存活和死亡两种情况，还包括了另一类结局事件即由于出现了竞争事件导致死亡结局的出现，相比较为原来的死亡结局的二分类变量，这里扩充为三分类变量。
处理完数据后剩余2551例，将数据处理并重新进行编码，并按照训练集：测试集为5:5划分，其中训练集为1275例，测试集为1276例。对训练集的数据进行相关模型的建立，基于构建的模型于测试集中进一步检验说明，分析预测的准确性和拟合效果。
Table 1. Variable descriptions and assignment representations
表1. 变量说明及赋值表示
3. 模型预测实证分析
本章首先构建随机生存森林预测模型，并用机器学习中的VIMP (variable importance)法和最小深度法结合对影响盲肠癌生存预后因素进行初步筛选，然后分别构建Cox比例风险模型和竞争风险模型对影响盲肠癌患者预后的因素进行分析确定。
3.1. 随机生存森林模型
用所有纳入分析的9个变量构建随机生存森林模型，表2展示出该模型在训练集和测试集上的结果，该模型默认生成500个二元生存树，平均每个生存树有15个终端节点，模型在训练集上的错误率为27.1%，而在测试集上的错误率为27.8%，测试集和训练集相差不大，效果较好。
Table 2. Random forest model training set and test set comparison
表2. 随机森林模型训练集和测试集对比
图1表示随机生存森林模型生成的生存树的数量与模型预测错误率大小的关系图，随着生存树数量的增加，其预测错误率明显降低；当生存树增加到一定数量后，预测错误率曲线趋于平稳(27.1%)，所以选择该随机生存森林模型树的数量选择合适。
随机生存森林模型可以对变量的重要度进行排名，因此利用该模型对影响盲肠癌生存预后因素进行初步筛选。图2表示VIMP法和最小深度法相结合的散点图。其中，蓝色点代表VIMP值大于0，红色则代表VIMP值小于0；主对角线以上的点表示VIMP排名更高，主对角线以下的点表示最小深度法排名更高。根据综合排名，决定去掉性别变量，则初步筛选剩下影响盲肠癌生存预后的因素为AJCC分期、肿瘤大小、年龄、婚姻状况、组织学分级、化疗状况、种族、放疗状况。
Figure 1. Models with different numbers of survival trees predict error rates
图1. 不同数量生存树下的模型预测错误率
Figure 2. The VIMP method and the minimum depth method combined variable selection
图2. VIMP法和最小深度法结合变量筛选
3.2. Cox比例风险模型
3.2.1. K-M生存曲线
本小节考虑基于机器学习方法初步筛选出来的变量构建Cox比例风险模型。图3~6分别为不同年龄、种族、AJCC分期、化疗状况、放疗状况的K-M生存曲线，并给出对数秩检验，由于篇幅限制，文中只给出前4个变量的生存曲线。除了种族变量没有通过对数秩检验，其它7个变量都通过了对数秩检验。如图3，年龄大于65岁的患者整体生存概率要低于年龄小于65岁的患者；图5中，AJCC分期为IV期的生存概率要比其它分期生存概率低的多，在生存时间为15个月时，生存概率就降到50%；值得注意的是图6中，选择放疗治疗的患者要比选择不放疗治疗的患者生存概率低很多，因为放疗治疗的患者很少，只有当盲肠癌病情更加严重的患者会选择进行放疗治疗，因此其生存概率低于不放疗治疗的患者。
3.2.2. 基于Cox比例风险模型的单因素与多因素分析
为了更进一步分析影响盲肠癌患者生存预后影响因素。表3为Cox比例风险模型的单因素与多因素分析表，单因素和多因素分析结果基本一致。在单因素分析表中，只有种族变量对应的p值均大于0.05，这说明这些变量对于生存预后而言不是独立的预后因素。年龄、AJCC分期、放疗状况、化疗状况、组织学分级、肿瘤大小、婚姻状况这些变量都大致通过了显著性检验，这些都是盲肠癌生存预后独立的预后因素。其中，对于年龄来说，HR危险比为1.72，表示明在其他协变量不变的情况下，年龄大于65岁患者死亡风险率为年龄小于65岁患者的1.72倍。
在多因素分析表中，种族对应的p值均大于0.05，这说明该变量对于生存预后而言不是独立的预后因素。组织学分级、肿瘤大小、婚姻状况这些变量大致都通过了显著性检验。其中，对于年龄来说，HR危险比为2.02，表示在其他协变量不变的情况下，年龄大于65岁患者死亡风险率为年龄小于65岁患者的2.04倍。
Figure 3. K-M survival curves for different Age
图3. 不同年龄的K-M生存曲线
Figure 4. K-M survival curves for different Race
图4. 不同种族的K-M生存曲线
Figure 5. K-M survival curves for different AJCC stage
图5. 不同AJCC分期的K-M生存曲线
Figure 6. K-M survival curves for different Radiotherapy conditions
图6. 不同放疗状况的K-M生存曲线
Table 3. Single-factor and multivariate analysis table for Cox proportional risk models
表3. Cox比例风险模型单因素与多因素分析表
3.2.3. 预测模型的构建
1) 比例风险假定检验
在使用Cox比例风险回归模型时，其模型存在一个前提假设即时间与协变量之间不存在交互作用(PH假定) [8]。本文主要采用时间依赖变量的Cox模型对前提假设条件进行检验。检验结果主要如表4所示，不难看出，所有变量和总体的p值均超过了显著性水平，这说明模型整体接受了原假设即认为时间与协变量之间不存在交互作用，满足了构建模型的前提。
2) Cox比例风险回归模型的建立
通过逐步回归的方式，即满足AIC最小的准则，对纳入分析的变量进行筛选和确定，使得最终模型满足该准则。所筛选出来的变量及系数等最终结果可见表5。
Table 5. Table of Cox proportional risk regression covariate coefficients
表5. Cox比例风险回归协变量系数表
从表5中可以看出最终纳入的分析变量有：年龄、组织学分级、肿瘤大小、AJCC分期、婚姻状况、放疗状况、化疗状况。最终确定的模型为：
其中，系数为正表示该因素为危险因素；系数为负表示保护因素。在该多因素比例风险回归模型中，化疗治疗、已婚、肿瘤直径大小在1 cm以上的为盲肠癌患者生存预后的保护因素，年龄大于65岁患者、放疗治疗、AJCC分期大于I、组织学等级高于一级、婚姻状况为其它的因素都为危险因素。根据查询，患者中选择放疗治疗的患者数量极少，只有病情特别严重的患者选择放疗治疗，但效果不佳。
3) 显著性检验
基于所构建的模型进行显著检验，具体检验结果可见表6。三种情况下的检验p值远远小于0.05，这说明，上述Cox多因素回归模型通过了显著性检验，模型总体效果不错并且具有合理性，同时这7个变量能较好的刻画出生存函数总的变化情况。
3.3. 竞争风险模型
3.3.1. Nelson-Aalen累计风险曲线
与Cox比例风险模型不同的是，考虑到导致出现感兴趣终点事件的发生存在着竞争事件。故这里使用通过计算每个结局的累积发生率函数(CIF)绘制累积生存曲线，并使用Fine-Gray’s检验来比较不同变量的组间之间的风险函数是否存在显著的差异。下文基于随机生存森林模型筛选出的8个变量绘制累积生存曲线，由于篇幅限制，本文给出前4个变量的累积生存曲线，图7~10分别代表的是：年龄、种族、AJCC分期、化疗状况。图中的实线表示的是变量中不同类别盲肠癌患者死亡的发生率，而虚线表示的是在死于竞争事件发生率。例如：图7在年龄变量中，年龄大于65岁的患者的死于盲肠癌累计发生率要远高于竞争事件的死亡发生率，年龄小于65岁的患者死于盲肠癌发生率也远高于竞争事件的死亡发生率。
结果表明：在控制竞争事件后，年龄(p < 0.01)、种族(p = 0.2)、组织学分级(p < 0.01)、肿瘤大小(p < 0.01)、AJCC分期(p < 0.01)、婚姻状况(p < 0.01)、化疗状况(p < 0.01)、放疗状况(p < 0.01)，其中除了种族没有通过显著性检验之外，其余变量都通过了显著性检验，同一种死亡原因的不同组之间死亡率存在差别。
3.3.2. 基于竞争风险模型的单因素与多因素分析
为了更加准确确定盲肠癌患者生存预后影响因素以及比较Cox比例风险模型确定的盲肠癌患者生存预后影响因素，本节考虑控制死亡原因为盲肠癌的竞争风险模型进行建模分析。单因素与多因素分析结果也基本一致。
表7为基于竞争风险模型的单因素与多因素分析表。在单因素分析表中，种族、婚姻状况对应的p值均大于0.05，这说明这些变量对于生存预后而言不是独立的预后因素。年龄、AJCC分期、放疗状况、化疗状况、组织学分级、肿瘤大小这些变量均通过了显著性检验。其中，对于年龄来说，HR危险比为1.17，表示明在其他协变量不变的情况下，年龄大于65岁患者死亡风险率为年龄小于65岁患者的1.17倍。
在多因素分析表中，种族、婚姻状况对应的p值均大于0.05，这说明该变量对于生存预后而言不是独立的预后因素。年龄、AJCC分期、放疗状况、化疗状况、组织学分级、肿瘤大小大致通过了显著性检验。其中，对于年龄来说，HR危险比为1.34，表示明在其他协变量不变的情况下，年龄大于65岁患者死亡风险率为年龄小于65岁患者的1.34倍。
Figure 7. CIF cumulative survival curves at different Age
图7. 不同年龄的CIF累积生存曲线
Figure 8. CIF cumulative survival curves at different Race
图8. 不同种族的CIF累积生存曲线
Figure 9. CIF cumulative survival curves at different AJCC stage
图9. 不同AJCC分期的CIF累积生存曲线
Figure10. CIF cumulative survival curves at different Chemotherapy conditions
图10. 不同化疗状况的CIF累积生存曲线
Table 7. Competitive risk model univariate and multivariate analysis table
表7. 竞争风险模型单因素与多因素分析表
3.3.3. 竞争风险预测模型的构建
通过上述单因素与多因素分析结果，选取年龄、AJCC分期、放疗状况、化疗状况、组织学分级、肿瘤大小这六个变量作为自变量，而因变量则是生存时间和生存结局(DSS)，区别于前面的Cox比例风险回归模型，这里的因变量发生了改变，考虑将竞争事件引入。
Table 8. Multi-factor competitive risk model coefficient table
表8. 多因素竞争风险模型系数表
基于逐步回归以及满足AIC最小的准则，构建了竞争风险模型，从表8中可以看出：除了肿瘤大小个别变量之外，其余变量均通过了显著性检验，说明对于生存预后而言是显著相关的。因此构建的竞争风险模型为：
其中，系数为正表示该因素为危险因素；系数为负表示保护因素。在该多因素竞争风险模型中，化疗治疗、肿瘤直径大小在1 cm以上变量为盲肠癌患者生存预后的保护因素，年龄大于65岁患者、AJCC分期大于I、组织学等级高于I级、放疗治疗等都为危险因素。分析结果基本与Cox比例风险模型一致。
4. 模型预测与评估
4.1. 生存预测列线图
图11和图12是根据Cox比例风险模型和竞争风险模型构建的预测患者一年、三年、五年生存概率列线图。根据列线图 [8] 可以诊断病人的在一年、三年、五年时的生存概率。例如在Cox比例风险模型图11预测中：某盲肠癌患者年龄大于65岁则年龄得分记为30；组织学评级为III级，得分为35分；肿瘤大小为小于1 cm，得分为40分；AJCC分期为III，得分为40分；没有放疗治疗，得分为0分，已婚得分为10分；没有化疗治疗，得分为30，则总分为185分，对应的一年生存概率为0.45、三年生存概率为0.1、五年的生存概率基本为0。而在竞争风险模型图12中患者年龄大于65岁则年龄得分记为8；组织学评级为III级，得分为25分；肿瘤大小为小于1 cm，得分为50分；AJCC分期为III，得分为55分；没有放疗治疗，得分为0分；没有化疗治疗，得分为15，则总分为153分，对应的一年生存概率为0.7左右，三年生存概率0.2，五年的生存概率为0.1。两种模型在生存概率预测上表现不同主要是因为竞争风险模型考虑的死亡原因除盲肠癌之外还有其它原因，所以竞争风险模型中表现的生存概率要比Cox比例风险模型生存概率高。
Figure 11. Cox proportional hazards model survival nomogram
图11. Cox比例风险模型生存列线图
Figure 12. Competing risk model survival nomogram
图12. 竞争风险模型生存列线图
4.2. 模型评估
图13和图14分别表示Cox比例风险模型和竞争风险模型在训练集上的ROC曲线 [9]，在训练集中，Cox比例风险模型预测的一年、三年、五年的ROC曲线对应的AUC值分别是0.799、0.790和0.792，而竞争风险模型对应的AUC值为0.848、0.841和0.830。图15和图16分别表示Cox比例风险模型和竞争风险模型在测试集上的ROC曲线，在测试集中，Cox比例风险模型三个预测时间点对应的AUC值分别为 0.809、0.779、0.752，竞争风险模型三个预测时间点对应的AUC值分别为0.850、0.840和0.825。不论是在训练集还是测试集上，竞争风险模型比Cox比例风险模型效果都要好，预测准确率更高。
Figure13. Cox proportional hazards model training set ROC curve
图13. Cox比例风险模型训练集ROC曲线
Figure14. Competitive risk model training set ROC curve
图14. 竞争风险模型训练集ROC曲线
Figure 15. Cox proportional hazards model test set ROC curve
图15. Cox比例风险模型测试集ROC曲线
Figure 16. Competitive risk model test set ROC curve
图16. 竞争风险模型测试集ROC曲线
5. 结论
通过建立随机生存森林模型、Cox比例风险模型、竞争风险模型进一步研究对于盲肠癌患者生存预后影响因素的探讨。结果表明：Cox比例风险回归模型中，化疗治疗、已婚、肿瘤直径大小在1 cm以上的为盲肠癌患者生存预后的保护因素，年龄大于65岁患者、放疗治疗、AJCC分期大于I、组织学等级高于一级、婚姻状况为其它的因素为危险因素；在竞争风险模型中，化疗治疗、肿瘤直径大小在1 cm以上变量为盲肠癌患者生存预后的保护因素，年龄大于65岁患者、AJCC分期大于I、组织学等级高于I级、放疗治疗等都为危险因素。竞争风险模型剔除了婚姻状况这一变量，其余分析结果与Cox比例风险模型基本一致。
随机生存森林模型在训练集和测试集上的错误率为0.271和0.278；Cox比例风险模型三年预测结果在训练集和测试集上的AUC值为0.790和0.779；竞争风险模型三年预测结果在训练集和测试集上的AUC值为0.841和0.840。模型的比较中，竞争风险模型更胜一筹。在对于存在竞争事件的生存分析中，选择基于竞争风险构建的预测模型不仅准确度高，而且更具合理性。