Summary
本协议描述了R中的代码，用于评估竞争风险模型的辨别和校准能力，以及用于内部和外部验证的代码。
Abstract
Cox比例风险模型广泛用于临床环境中的生存分析，但它无法应对多种生存结果。与传统的Cox比例风险模型不同，竞争风险模型考虑竞争事件的存在及其与列线图（一种图形计算设备）的组合，列线图是临床医生进行精确预后预测的有用工具。在这项研究中，我们报告了一种建立竞争风险列线图的方法，即评估其辨别力（即相关指数和曲线下面积）和校准（即校准曲线）能力，以及净收益（即决策曲线分析）。此外，还使用原始数据集的自举重采样进行内部验证，并使用已建立的竞争风险列线图的外部数据集进行外部验证，以证明其外推能力。竞争风险列线图应作为临床医生在考虑竞争风险的情况下预测预后的有用工具。
Introduction
近年来，随着精准医学的发展，新兴的预后因素被确定出来，结合分子和临床病理因素的预后模型在临床环境中越来越受到关注。然而，非图形模型，如Cox比例风险模型，具有系数值的结果，临床医生很难理解1。相比之下，列线图是回归模型（包括 Cox 回归模型、竞争风险模型等）的可视化工具，回归模型是专为数学函数2 的近似图形计算而设计的二维图。它能够评估临床模型中不同级别的变量，并计算风险评分（RS）以预测预后。
模型评估在模型构建中是必不可少的，并且通常接受两个特征进行评估：判别和校准。在临床模型中，鉴别是指模型将发生事件的个体与未发生事件的个体（例如死亡患者与存活患者）区分开来的能力，并且一致性指数（C指数）或受试者工作特征曲线下面积（AUC）通常用于表征它3,4.校准是将模型的预测概率与实际概率进行比较的过程，校准曲线已被广泛用于表示它。此外，模型验证（内部和外部验证）是模型构建中的重要步骤，只有经过验证的模型才能进一步外推5。
Cox比例风险模型是医学研究中使用的回归模型，用于调查预后因素与生存状态之间的关联。然而，Cox比例风险模型只考虑结果的两种状态[Y（0,1）]，而研究对象通常面临两种以上的状态，并且出现竞争风险[Y（0,1,2）]1。总生存期（OS）定义为从起源日期（例如治疗）到因任何原因死亡的时间，是生存分析中最重要的终点。然而，OS未能区分癌症特异性死亡和非癌症特异性死亡（例如，心血管事件和其他不相关原因），因此忽略了竞争风险6。在这些情况下，竞争风险模型是预测生存状态的首选，同时考虑竞争风险7。构建和验证Cox比例风险模型的方法已经确立，而关于验证竞争风险模型的报告很少。
在我们之前的研究中，建立了特定的竞争风险列线图，列线图和竞争风险模型的组合以及基于竞争风险模型的风险评分估计8。本研究旨在提出评估和验证已建立的竞争风险列线图的不同方法，该方法应作为临床医生在考虑竞争风险的情况下预测预后的有用工具。
Subscription Required. Please recommend JoVE to your librarian.
Protocol
监测、流行病学和最终结果 （SEER） 数据库是一个开放访问的癌症数据库，仅包含去识别的患者数据 （SEER ID： 12296-Nov2018）。因此，本研究免于浙江大学医学院附属金华医院评审委员会的批准。
1. 数据准备和R包准备
- 准备并导入数据。
>数据集<-读取.csv（“.../乳腺癌数据.xlsx”） #Import 数据。
注意：数据在 补充文件 1 中上传。 - 安装并加载 R 包。
>包<- c（“rms”，“cmprsk”，“mstate”，“survival”，“riskRegression”，”
普罗德利姆”）
> req.pcg <- function（pcg）{
新<- pcg[！（pcg %in% installed.packages（）[， “Package”]）]
if （length（new）） install.packages（new， dependencies = T）
sapply（pcg， require， ch = T）
}
> req.pcg（packages）
注意：使用 rms、 cmprsk、 mstate、 survival、 riskRegression 和 prodlim （http://www.r-projectrg/） 软件包基于 R 软件（版本 3.6.2）执行以下过程。
2. 以两种不同的方法建立竞争风险列线图
- 以直接方法建立竞争风险列线图。
> mod_cph <- cph（生存月数，状态） ~ 因子1+ 因子2+...，
x=T， y=T， surv=T， data=Dataset）
> nom <- 列线图（mod_cph， fun=list（function（x） 1-surv_cph（36， x）...），
funlabel=c（“3-year event1 Prob.”...）， lp=F） #Take 第 36 个月为例。
> mod_crr <- crr（生存月数，状态，失败代码=1，cov1=cov）
> score <- log（log（（1-real.3y），（1-cif.min36）））/（maxbeta/100）
>剧情（标称） - 以加权方法建立竞争风险列线图。
> df.w <- crprep（“生存月”， fstatus“，
data=Dataset， trans=c（1,2）， cens=0，
keep=c（“factor1”，“ factor2”...））
> mod.w <- cph（Surv（Tstart， Tstop， status==1） ~ factor1+factor2+...，
data=df.w， weight=weight.cens， subset=failcode==1， surv=T）
> nom.w <- 列线图（mod.w）
3. 竞争风险列线图的判别能力
- 歧视的C指数
- 将矩阵 cov 拟合到竞争风险模型 mod_crr ，并获得预测的矩阵 SUV。
> SUV <- predict.crr（mod_crr， cov） - 从 suv 获取某个月份的累积发病率，并使用函数 rcorr.cens 计算 C 指数。
> CIF36 <- SUV[其中（SUV[，1]==36），][-1]
> rcorr <- rcorr.cens（1-cif36，Surv（Dataset$Survivalmonths，Dataset$tumordeath））
> cindex <- rcorr[1]
- 将矩阵 cov 拟合到竞争风险模型 mod_crr ，并获得预测的矩阵 SUV。
- 歧视的AUC
- 使用函数 Score （风险回归包）对竞争风险模型的预测性能进行评分。
> fgr.w <- FGR（Hist（生存月，状态）~ 因子1+ 因子2+...，数据=数据集，原因=1）
> 分数 <- 分数（列表（“细灰” = fgr.w）， - 从“分数”中提取 AUC。
>分数$AUC
- 使用函数 Score （风险回归包）对竞争风险模型的预测性能进行评分。
4. 竞争风险模型的校准能力
- 竞争风险模型置信区间为 95% 的校准曲线
- 获取包含每个人在特定故障时间内的累积发生率的数据框。
> CIF36 <- data.frame（CIF36） #Take 第 36 个月为例。
> 同名（cif36.36_o）<-c（“36m”） - 根据估计的累积发病率将队列划分为五个子组，并计算每个子组的平均预测累积发病率。
>组36 <- 切割（CIF36$'36m'，
分位数（CIF36$'36M'， seq（0， 1， 0.2）），
包含。最低 = TRUE，标签 = 1：5）
> mean36 <- as.vector（by（cif36 $'36m'， group36， mean）） - 使用函数 cuminc 计算观测到的累积发生率，即实际累积发生率，然后在某个失效时间内以 95% 置信区间获得观测到的累积发生率。
> cum36 <- cuminc（Dataset$Survivalmonths，Dataset$fstatus，group36）
> obs36 <- 时间点（暨36，数据集$生存月）$est[c（1：5），36]
> obs36var <- timepoints（cum36，Dataset$Survivalmonths）$var[c（1：5），36]
> df <- data.frame（mean36， obs36， obs36var） - 使用 ggplot 函数绘制校准曲线，将预测的累积发生率作为 x 轴，将观测的累积发生率绘制为 y 轴。
> ggplot（df）+ geom_point（aes（x=mean36，y=obs36），col=“red”）+
geom_point（AES（x=mean36，y=obs36），col=“red”，pch=4）+
geom_line（col=“red”，aes（x=mean36，y=obs36））+
geom_errorbar（col=“red”，aes（x=mean36，y=obs36+1.96
*sqrt（obs36var）），
ymin =obs36-1.96*sqrt（obs36var）， ymax = obs36+1.96
*sqrt（obs36var））
geom_abline（lty=3，lwd=2，col=c（rgb（0,118,192，
最大颜色值=255）））
- 获取包含每个人在特定故障时间内的累积发生率的数据框。
- 具有竞争风险模型风险评分的校准曲线
- 评估所有变量的每个水平并获得总 RS。
> Dataset$factor1[Dataset$factor1==1] <- factor1.scale[“Factor1_level1”]
> ... #For 例子，Dataset$histology[Dataset$histology==1]<-histology.scale[“Histology1”]
> Dataset$rs <- Dataset$factor1+Dataset$Factor2+Dataset$Factor3+...
注意：通过对每个变量的点求和来获得每位患者的总 RS。 - 计算频率并计算观察到的不同总风险评分的累积发生率。
> rs.freq <- as.data.frame（table（Dataset$rs））
> obs.36 <- vector（mode=“numeric”， length=nrow（rs.freq））
> for （i in 1： nrow（rs.freq）） {
数据集 <- 子集（Dataset，Dataset$rs== rs.freq [i，1]）
cif.dataset <- cuminc（dataset$Survivalmonths，dataset$death3）
CIF36.数据集<-时间点（CIF.dataset，36）
obs.36[i] <- cif36.dataset$est[1]} - 设置 x 轴的范围并计算总风险评分的预测累积发生率。
> RS <- 范围（nom$总点）
> x.36 <- 序列（最小值），最大值（RS），0.01）
> pre.36 <- 1-（1-cif.min36）^exp（x.36*maxbeta/100） - 使用风险评分绘制校准曲线。
> 图（x.36， pre.36， type='l'...
> par（new=TRUE）
> plot（as.vector（rs.freq[，1]）， obs.36... ）
- 评估所有变量的每个水平并获得总 RS。
5. 竞争风险模型的决策曲线分析
- 获取 stdca 函数以执行决策曲线分析。
>来源（“标准。R") - 从列线图中提取多项式方程以计算生存概率。
> 列线图Ex（nomo = nom）
> 数据集$预测因子 <- A * （Dataset$rs ^3） + B * （Dataset$rs ^2） + C * Dataset$rs + D #predictors 是由已建立的列线图计算的癌症特异性死亡的预测概率 - 执行决策曲线分析。
> stdca（data = 数据集， 结果 = “状态”， ttresult = “生存月份”， 时间点 = 36，
预测变量 = “预测变量”，cmprsk = 真，平滑 = 假，概率 = 假）
注意：为了在存在竞争风险的情况下评估结果，应为 cmprsk 选择 TRUE。
6. 使用引导方法进行内部验证
- 使用自举方法获取平均预测的累积发生率。
- 使用替换对原始数据集（数据集）进行重采样以生成引导数据集 （Dataset_in）。使用引导数据集建立竞争风险模型 （mod.in_crr）。使用函数 predict.crr 预测mod.in_crr和循环 b 时间以生成 suvall.in。
B=b
suvall.in <- 列表（）
for（j in 1：B）{
Dataset_in <- 数据集[sample（c（1：nrow（Dataset）），nrow（Dataset），
替换 = 真），]
附（Dataset_）
科弗。in <- model.matrix（~factor1+ factor2+...）[,-1]
国防部。in _crr <- crr（Survivalmonths， fstatus， failcode=1， cov1=cov.in）
分离（数据集.
SUV.在<- predict.crr（mod. in _crr， cov）
suvall.in[[j]] <- suv.in} - 获取某个月的平均预测累积发生率。
CIF36全部。 内部 <- 向量（mode=“numeric”， length=nrow（Dataset））
对于 （1：B 中的 k） {
CIF36全部。内部<- CIF36全部。内+苏瓦尔。inner[[k]][which（suvall. inner[[k]][，1]==36），][-1]
}
cif36.in <- cif36all.in/B
- 使用替换对原始数据集（数据集）进行重采样以生成引导数据集 （Dataset_in）。使用引导数据集建立竞争风险模型 （mod.in_crr）。使用函数 predict.crr 预测mod.in_crr和循环 b 时间以生成 suvall.in。
- 使用函数 rcorr.cens 的内部交叉验证计算 C 指数。
雷科尔。内<- rcorr.cens（1-cif36.in，Surv（Dataset$Survivalmonths，Dataset$tumordeath））
索引。内<- rcorr.内[1] - 使用交叉内部验证进行校准。
注意：具有内部验证的竞争风险模型的校准曲线代码类似于第 4 节中的代码，而 suv 被 suv.in 取代。
7. 竞争风险模型的外部验证
- 使用外部数据获取预测的累积发生率。使用外部数据变量矩阵 （cov.ex） 获取预测的累积发生率。
suv.ex <- predict.crr（mod_crr，cov.ex）
cif36.ex <- suv.ex [which（suv.ex $time==“36”），][-1] - 使用外部验证计算 C 指数。
rcorr.ex <- rcorr.cens（1-cif36.ex，Surv（Dataset.ex$Survivalmonths，Dataset.ex$tumordeath））
cindex.ex <- rcorr.ex[1] - 使用外部验证进行校准。
注意：具有内部验证的竞争风险模型的校准曲线代码类似于第 4 节中的代码，而 suv 由 suv.ex 替换。
Subscription Required. Please recommend JoVE to your librarian.
Representative Results
在这项研究中，从SEER数据库中检索了乳腺癌患者的数据，并作为示例数据。SEER数据库提供了占美国人口约34.6%的癌症数据，并获得了访问该数据库的许可（参考编号12296-Nov2018）。
采用直接法和加权法分别建立了两个列线图（图1），包括组织学类型、分化分级、T分期和N期。每个变量水平的点数和总点对应的概率几乎相同，但观察到一些细微的差异。Zhang等人引入了一种“加权”方法来建立竞争风险列线图，该方法首先将原始数据转换为加权数据（使用函数crprep），然后使用加权数据构建Cox回归模型（使用函数coxph），最后使用Cox回归模型（使用函数列线图）建立竞争风险列线图9.相比之下，本研究中的“直接”方法与“加权”方法完全不同。简而言之，从竞争风险模型（使用函数crr）生成的参数取代了Cox回归模型中的参数（使用函数coxph），该模型最终用于建立竞争风险列线图（使用函数列线图）。在比较通过“加权”法和“直接”法确定的列线图时，两个列线图总体上相似，但可以观察到一些细微的差异。研究中的“直接”方法更精确，因为它直接获得构建列线图的参数（在第2节中使用公式“score=log（log（（1-real.3y），（1-cif.min36）））/（maxbeta/100）”）。
在 rcorr.cens（X， Surv） 中， X 是具有任何时间点累积发生率的数字预测因子，而 Surv 是包含生存月份和状态的生存对象。当 Surv 将经历竞争事件的患者定义为删失，然后生成生命表时，可评估的有序对是相同的。未经验证的竞争风险模型的C指数为0.7978（95%CI = 0.7650-0.8305），表明该模型具有中等判别能力。自举分析重复500次，然后将500个结果平均以产生用于计算C指数的单一估计。内部验证中的C指数为0.7978（95%CI = 0.7651-0.8305），与原始数据集中的C指数相似。将外部数据集拟合到竞争风险模型中，外部验证中的C指数为0.5071（95%CI = 0.4637-0.5505）。竞争风险模型的AUC是根据研究中的原始数据集计算的。第36个月的AUC为0.8461（95%CI = 0.8232-0.8691），证明了该模型的辨别能力。
如图2A所示，校准曲线上的点接近等效线，每组观测频率的95%置信区间落入等效线，表明模型的精确校准能力。使用内部和外部验证的校准曲线分别如图2B和图2C所示，表明构建的模型在内部验证中具有良好的校准能力，但在外部验证中具有较差的校准能力。
如 补充图1所示，表示观测到的累积发生率的点分布在表示预测累积发生率的线周围，观测到的发生率和预测的发生率之间没有观察到显著差异。决策曲线分析的结果如图 3所示，其中显示了净收益随阈值概率的增加而变化。
图1：用两种方法建立竞争风险列线图。 （A）使用直接方法建立列线图。（B） 使用加权法建立的列线图。组织学：1、浸润性导管癌;2、浸润性小叶癌;3、浸润性导管癌+浸润性小叶癌。等级：1，分化良好;2、中度分化;3、分化差。T级：1、T1级;2、T2级;3、T3阶段;4、T4级。N级：0，N0级;1、N1级;2、N2级;3、N3级。缩写：CSD = 癌症特异性死亡。请点击此处查看此图的大图。
图 2：竞争风险列线图的校准曲线。 （A）具有已建立竞争风险模型置信区间的校准曲线。（B）内部验证中竞争风险模型的校准曲线。（C）外部验证中竞争风险模型的校准曲线。请点击此处查看此图的大图。
图 3：竞争风险列线图的决策曲线分析。净收益与阈值概率作图。“全部”线通过考虑所有遭受癌症特异性死亡的患者来显示净收益，“无”线是通过考虑所有未遭受癌症特异性死亡的患者来显示净收益。请点击此处查看此图的大图。
补充文件1：乳腺癌数据。 职称定义：等级修改，差异化等级;组织学，组织学类型;T期，肿瘤T期;N期，肿瘤N期;生存月，从治疗日期到因任何原因或审查员死亡的时间;死亡、死亡（包括癌症特异性死亡和非癌症特异性死亡）或审查员;死亡3、癌症特异性死亡、非癌症特异性死亡或审查员。 请点击此处下载此文件。
补充图1：已建立的竞争风险模型的风险评分校准曲线。 表示观测到的累积发生率的点分布在表示预测累积发生率的线周围。 请点击此处下载此文件。
Subscription Required. Please recommend JoVE to your librarian.
Discussion
这项研究比较了通过两种不同方法建立的竞争风险列线图，并对已建立的列线图进行了评估和验证。具体而言，本研究提供了一个基于直接方法建立列线图的分步教程，以及计算C指数和绘制校准曲线。
R软件中的均方根包广泛用于Cox比例风险模型的构建和评估，但不适用于竞争风险模型。对于具有多个结果的模型，Zhang等人报告了使用R软件中的riskRegression包对竞争风险模型的验证，该包计算AUC和Brier分数以评估辨别能力，并绘制校准曲线以验证校准能力9,10。
但是，riskRegression 包有一些缺点，例如无法计算 C 指数。C 指数是一致可评估有序对在所有可评估有序对中的比例（C 指数 = 一致可评估有序对/所有可评估有序对）3,11。在传统的双状态结局生存分析中，只有两个患者都死亡（Di = 1，Dj = 1）或死亡患者的生存月数比失败的患者短的配对才能成为可评估的有序对（Di = 1，Dj = 0，Ti < Tj）（i = 患者 i，j = 患者 j，T = 时间， D = 状态，1 = 死亡，0 = 失败）。在竞争风险的背景下，从竞争风险中失败的患者仍然处于风险集中，因为它只能推断经历竞争事件的患者的生存月数比观察到的要长。此外，只有两个患者都遭受死亡（Di = 1，Dj = 1）或遭受死亡的患者的生存月数比失败或经历竞争事件的患者短的配对才能成为可评估的有序对（Di = 1，Dj = 0|2，Ti < Tj）（i = 患者 i，j = 患者 j，T = 时间， D = 状态，2 = 竞争事件，1 = 死亡，0 = 失败）。因此，本研究使用函数rcorr.cens来计算C指数。
在竞争风险模型的内部验证中，本研究应用了自举法，并证明了构建模型12的良好性能。随机拆分为训练数据集和测试数据集可能会有问题，因为测试数据集可能特别容易（或难以）预测13。K-fold方法也可用于交叉验证，而在模型验证5,14中应用的频率较低。这项研究还对竞争风险模型进行了外部验证，但在外部环境中表现不佳。这可能是因为我们用于外部验证的数据只是原始数据的重新采样。
然而，这项研究有几个局限性。首先，我们的方法基于R软件;因此，用户需要一定的编程知识，这可能会限制其目标受众。此外，还有数百行代码，有些代码需要针对不同的数据进行更改;我们希望在未来的研究中开发一个“多合一”的R包，可以应用于各种数据。这项研究没有其他乳腺癌数据来进行外部验证，别无选择，只能从原始数据中重新采样，但外部验证的方法和代码是相同的。最重要的是，研究中假设的协变量和结果之间的线性在现实世界的研究中可能不成立，应考虑相互作用和非线性，对此集成建模可能会有所帮助15。
总之，本研究以“直接”方法建立了竞争风险列线图，并评估了其在原始、内部和外部数据集中的辨别和校准能力。希望竞争风险列线图将作为R中 riskRegression 包的补充，并为处理临床竞争风险事件提供帮助。
Subscription Required. Please recommend JoVE to your librarian.
Disclosures
提交人声明他们没有竞争利益。
Acknowledgments
该研究得到了浙江省医学科技计划项目（批准号2013KYA212）、浙江省自然科学基金一般计划（批准号Y19H160126）和金华市科技局重点计划（批准号2016-3-005,2018-3-001d和2019-3-013）的资助。
Materials
| Name | Company | Catalog Number | Comments |
| R software | None | Not Applicable | Version 3.6.2 or higher |
| Computer system | Microsoft | Windows 10 | Windows 10 or higher |
References
- Andersen, P. K., Gill, R. D. Cox's regression model for counting processes: A large sample study. The Annals of Statistics. 10 (4), 1100-1120 (1982).
- Lubsen, J., Pool, J., vander Does, E. A practical device for the application of a diagnostic or prognostic function. Methods of Information in Medicine. 17 (2), 127-129 (1978).
- Harrell, F. E., Lee, K. L., Mark, D. B. Multivariable prognostic models: Issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics In Medicine. 15 (4), 361-387 (1996).
- Hung, H., Chiang, C. -T. Estimation methods for time-dependent AUC models with survival data. The Canadian Journal of Statistics / La Revue Canadienne de Statistique. 38 (1), 8-26 (2010).
- Moons, K. G. M., et al. Risk prediction models: I. Development, internal validation, and assessing the incremental value of a new (bio)marker. Heart. 98 (9), 683-690 (2012).
- Fu, J., et al. Real-world impact of non-breast cancer-specific death on overall survival in resectable breast cancer. Cancer. 123 (13), 2432-2443 (2017).
- Fine, J. P., Gray, R. J. A proportional hazards model for the subdistribution of a competing risk. Journal of the American Statistical Association. 94 (446), 496-509 (1999).
- Wu, L., et al. Establishing a competing risk regression nomogram model for survival data. Journal of Visualized Experiments. (164), e60684 (2020).
- Zhang, Z., Geskus, R. B., Kattan, M. W., Zhang, H., Liu, T. Nomogram for survival analysis in the presence of competing risks. Annals of Translational Medicine. 5 (20), 403 (2017).
- Zhang, Z. H., et al. Overview of model validation for survival regression model with competing risks using melanoma study data. Annals Of Translational Medicine. 6 (16), 325 (2018).
- Newson, R. Confidence intervals for rank statistics: Somers' D and extensions. Stata Journal. 6 (3), 309-334 (2006).
- Davison, A. C., Hinkley, D. V., Schechtman, E.
- Roecker, E. B. Prediction error and its estimation for subset-selected models. Technometrics. 33 (4), 459-468 (1991).
- Steyerberg, E. W., Harrell, F. E. Prediction models need appropriate internal, internal-external, and external validation. Journal of Clinical Epidemiology. 69, 245-247 (2016).
- Zhang, Z., Chen, L., Xu, P., Hong, Y. Predictive analytics with ensemble modeling in laparoscopic surgery: A technical note. Laparoscopic, Endoscopic and Robotic Surgery. 5 (1), 25-34 (2022).
