竞争风险模型的基于 R 的景观验证

摘要

本协议描述了R中的代码，用于评估竞争风险模型的辨别和校准能力，以及用于内部和外部验证的代码。

摘要

Cox比例风险模型广泛用于临床环境中的生存分析，但它无法应对多种生存结果。与传统的Cox比例风险模型不同，竞争风险模型考虑竞争事件的存在及其与列线图（一种图形计算设备）的组合，列线图是临床医生进行精确预后预测的有用工具。在这项研究中，我们报告了一种建立竞争风险列线图的方法，即评估其辨别力（即相关指数和曲线下面积）和校准（即校准曲线）能力，以及净收益（即决策曲线分析）。此外，还使用原始数据集的自举重采样进行内部验证，并使用已建立的竞争风险列线图的外部数据集进行外部验证，以证明其外推能力。竞争风险列线图应作为临床医生在考虑竞争风险的情况下预测预后的有用工具。

引言

近年来，随着精准医学的发展，新兴的预后因素被确定出来，结合分子和临床病理因素的预后模型在临床环境中越来越受到关注。然而，非图形模型，如Cox比例风险模型，具有系数值的结果，临床医生很难理解¹。相比之下，列线图是回归模型（包括 Cox 回归模型、竞争风险模型等）的可视化工具，回归模型是专为数学函数² 的近似图形计算而设计的二维图。它能够评估临床模型中不同级别的变量，并计算风险评分（RS）以预测预后。

模型评估在模型构建中是必不可少的，并且通常接受两个特征进行评估:判别和校准。在临床模型中，鉴别是指模型将发生事件的个体与未发生事件的个体（例如死亡患者与存活患者）区分开来的能力，并且一致性指数（C指数）或受试者工作特征曲线下面积（AUC）通常用于表征它^3,4.校准是将模型的预测概率与实际概率进行比较的过程，校准曲线已被广泛用于表示它。此外，模型验证（内部和外部验证）是模型构建中的重要步骤，只有经过验证的模型才能进一步外推⁵。

Cox比例风险模型是医学研究中使用的回归模型，用于调查预后因素与生存状态之间的关联。然而，Cox比例风险模型只考虑结果的两种状态[Y（0,1）]，而研究对象通常面临两种以上的状态，并且出现竞争风险[Y（0,1,2）]¹。总生存期（OS）定义为从起源日期（例如治疗）到因任何原因死亡的时间，是生存分析中最重要的终点。然而，OS未能区分癌症特异性死亡和非癌症特异性死亡（例如，心血管事件和其他不相关原因），因此忽略了竞争风险⁶。在这些情况下，竞争风险模型是预测生存状态的首选，同时考虑竞争风险⁷。构建和验证Cox比例风险模型的方法已经确立，而关于验证竞争风险模型的报告很少。

在我们之前的研究中，建立了特定的竞争风险列线图，列线图和竞争风险模型的组合以及基于竞争风险模型的风险评分估计⁸。本研究旨在提出评估和验证已建立的竞争风险列线图的不同方法，该方法应作为临床医生在考虑竞争风险的情况下预测预后的有用工具。

研究方案

监测、流行病学和最终结果 （SEER） 数据库是一个开放访问的癌症数据库，仅包含去识别的患者数据 （SEER ID: 12296-Nov2018）。因此，本研究免于浙江大学医学院附属金华医院评审委员会的批准。

1. 数据准备和R包准备

1. 准备并导入数据。
   >数据集<-读取.csv（".../乳腺癌数据.xlsx"） #Import 数据。
   注意:数据在 补充文件 1 中上传。
2. 安装并加载 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. 以两种不同的方法建立竞争风险列线图

1. 以直接方法建立竞争风险列线图。
   > 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）
   >剧情（标称）
2. 以加权方法建立竞争风险列线图。
   > 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. 竞争风险列线图的判别能力

1. 歧视的C指数
   1. 将矩阵 cov 拟合到竞争风险模型 mod_crr ，并获得预测的矩阵 SUV。
      > SUV <- predict.crr（mod_crr， cov）
   2. 从 suv 获取某个月份的累积发病率，并使用函数 rcorr.cens 计算 C 指数。
      > CIF36 <- SUV[其中（SUV[，1]==36），][-1]
      > rcorr <- rcorr.cens（1-cif36，Surv（Dataset$Survivalmonths，Dataset$tumordeath））
      > cindex <- rcorr[1]
2. 歧视的AUC
   1. 使用函数 Score （风险回归包）对竞争风险模型的预测性能进行评分。
      > fgr.w <- FGR（Hist（生存月，状态）~ 因子1+ 因子2+...，数据=数据集，原因=1）
      > 分数 <- 分数（列表（"细灰" = fgr.w），
   2. 从"分数"中提取 AUC。
      >分数$AUC

4. 竞争风险模型的校准能力

1. 竞争风险模型置信区间为 95% 的校准曲线
   1. 获取包含每个人在特定故障时间内的累积发生率的数据框。
      > CIF36 <- data.frame（CIF36） #Take 第 36 个月为例。
      > 同名（cif36.36_o）<-c（"36m"）
   2. 根据估计的累积发病率将队列划分为五个子组，并计算每个子组的平均预测累积发病率。
      >组36 <- 切割（CIF36$'36m'，
      分位数（CIF36$'36M'， seq（0， 1， 0.2）），
      包含。最低 = TRUE，标签 = 1:5）
      > mean36 <- as.vector（by（cif36 $'36m'， group36， mean））
   3. 使用函数 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）
   4. 使用 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）））
2. 具有竞争风险模型风险评分的校准曲线
   1. 评估所有变量的每个水平并获得总 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。
   2. 计算频率并计算观察到的不同总风险评分的累积发生率。
      > 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]}
   3. 设置 x 轴的范围并计算总风险评分的预测累积发生率。
      > RS <- 范围（nom$总点）
      > x.36 <- 序列（最小值），最大值（RS），0.01）
      > pre.36 <- 1-（1-cif.min36）^exp（x.36*maxbeta/100）
   4. 使用风险评分绘制校准曲线。
      > 图（x.36， pre.36， type='l'...
      > par（new=TRUE）
      > plot（as.vector（rs.freq[，1]）， obs.36... ）

5. 竞争风险模型的决策曲线分析

1. 获取 stdca 函数以执行决策曲线分析。
   >来源（"标准。R")
2. 从列线图中提取多项式方程以计算生存概率。
   > 列线图Ex（nomo = nom）
   > 数据集$预测因子 <- A * （Dataset$rs ^3） + B * （Dataset$rs ^2） + C * Dataset$rs + D #predictors 是由已建立的列线图计算的癌症特异性死亡的预测概率
3. 执行决策曲线分析。
   > stdca（data = 数据集， 结果 = "状态"， ttresult = "生存月份"， 时间点 = 36，
   预测变量 = "预测变量"，cmprsk = 真，平滑 = 假，概率 = 假）
   注意:为了在存在竞争风险的情况下评估结果，应为 cmprsk 选择 TRUE。

6. 使用引导方法进行内部验证

1. 使用自举方法获取平均预测的累积发生率。
   1. 使用替换对原始数据集（数据集）进行重采样以生成引导数据集 （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}
   2. 获取某个月的平均预测累积发生率。
      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
2. 使用函数 rcorr.cens 的内部交叉验证计算 C 指数。
   雷科尔。内<- rcorr.cens（1-cif36.in，Surv（Dataset$Survivalmonths，Dataset$tumordeath））
   索引。内<- rcorr.内[1]
3. 使用交叉内部验证进行校准。
   注意:具有内部验证的竞争风险模型的校准曲线代码类似于第 4 节中的代码，而 suv 被 suv.in 取代。

7. 竞争风险模型的外部验证

1. 使用外部数据获取预测的累积发生率。使用外部数据变量矩阵 （cov.ex） 获取预测的累积发生率。
   suv.ex <- predict.crr（mod_crr，cov.ex）
   cif36.ex <- suv.ex [which（suv.ex $time=="36"），][-1]
2. 使用外部验证计算 C 指数。
   rcorr.ex <- rcorr.cens（1-cif36.ex，Surv（Dataset.ex$Survivalmonths，Dataset.ex$tumordeath））
   cindex.ex <- rcorr.ex[1]
3. 使用外部验证进行校准。
   注意:具有内部验证的竞争风险模型的校准曲线代码类似于第 4 节中的代码，而 suv 由 suv.ex 替换。

结果

在这项研究中，从SEER数据库中检索了乳腺癌患者的数据，并作为示例数据。SEER数据库提供了占美国人口约34.6%的癌症数据，并获得了访问该数据库的许可（参考编号12296-Nov2018）。

采用直接法和加权法分别建立了两个列线图（图1），包括组织学类型、分化分级、T分期和N期。每个变量水平的点数和总点对应的概率几乎相同，但观察到一些细微的差异。Zha...

讨论

这项研究比较了通过两种不同方法建立的竞争风险列线图，并对已建立的列线图进行了评估和验证。具体而言，本研究提供了一个基于直接方法建立列线图的分步教程，以及计算C指数和绘制校准曲线。

R软件中的均方根包广泛用于Cox比例风险模型的构建和评估，但不适用于竞争风险模型。对于具有多个结果的模型，Zhang等人报告了使用R软件中的riskRegression包对竞?...

材料

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