R语言使用BOOT重抽样获取cox回归方程C-index(C指数)可信区间(2)

bootstrap自采样目前广泛应用与统计学中,其原理很简单就是通过自身原始数据抽取一定量的样本(也就是取子集),通过对抽取的样本进行统计学分析,然后继续重新抽取样本进行分析,不断的重复这一过程N(大于500次以上)次,然后得到N个统计结果,然后进行区间分析,得到最终结果。
在这里插入图片描述
上一章我们简单介绍了BOOT重抽样获取回归方程系数95%可信区间,可能大家对BOOT重抽样的用处感觉还不是很明显。BOOT重抽样在我们统计中处理数据还是很有用的,本期我们来介绍一下怎么使用BOOT重抽样获取cox回归方程C-index(C指数)可信区间,这也是一个粉丝向我问的问题,我觉得蛮有典型性和实用性的,因此就拿出来讲讲。首先我们看看什么是C-index(C指数),C-index,C指数即一致性指数(concordance index),用来评价模型的预测能力。c指数是指所有病人对子中预测结果与实际结果一致的对子所占的比例。我们在既往的文章《手把手教你使用R语言建立COX回归并画出列线图(Nomogram)》中已经介绍了怎么计算C指数,我们继续以原来文章的数据和方法为例进行演示。
我们先导入数据和R包,这里使用的是survival包的肺癌数据,为什么我们要导入rms包呢,因为我们等会需要Hmisc包的rcorrcens函数来计算C指数

library(survival)
library(foreign)
library(rms)
bc<-cancer
bc <- na.omit(bc)

在这里插入图片描述
我们来看看数据,inst: 机构代码,time: 以天为单位的生存时间,status: 状态:审查状态 1=审查,2=死亡,age: 年龄,sex: 男=1 女=2,ph.ecog:由医师评定的 ECOG 表现评分。ph.karno:由医师评定的 Karnofsky 表现评分(差=0-好=100),pat.karno:由患者评定的 Karnofsky 性能评分,meal.cal:用餐时消耗的卡路里,wt.loss:过去六个月的体重减轻。
我们随便挑几个变量组成COX回归方程

f <- cph(Surv(time, status) ~ age + sex + ph.ecog + pat.karno +wt.loss, 
         x=T, y=T, surv=T, data=bc)

计算C指数

rcorrcens(Surv(time, status) ~ predict(f), data = bc)

在这里插入图片描述
1-0.344=0.656,即为C-index
虽然我们算出C指数了,但是rcorrcens函数并没有提供它的可信区间,如果您的论文需要提供这个数据,我们可以通过BOOT重抽样获取。
BOOT重抽样其实就是对函数的反复抽样,因此关键在于设计好抽样函数,这是后台粉丝发给我的代码,说计算不出来结果

c_index <- function(formula, data, indices) {
tran.data <- data[indices,]
vali.data <- data[-indices,]
fit <- coxph(formula, data=tran.data)
result<-survConcordance(Surv(vali.data$time,vali.data$death)~predict(fit,vali.data))
index<-as.numeric(result$concordance)
return(index)
}

我觉得问题主要在分成验证集和建模集这一句,

tran.data <- data[indices,]
vali.data <- data[-indices,]

BOOT重抽样本来就是取子集,没有必要再分成两个数据集,如果你一定要这样取也是可以的,换成另外的方法,我后期会说到。
现在我对粉丝的函数稍稍修改,改成如下:

c_index <- function(data,indices){
  dat <- data[indices,]
  vames<-c("age", "sex", "ph.ecog","pat.karno","wt.loss")
  FML <- as.formula(paste('Surv(time, status)~',paste(vames, collapse = "+")))
  fit<- coxph(FML,data =dat )
  pr1<-predict(fit,newdata=dat)
  Cindex=rcorrcens(Surv(time, status) ~ pr1, data =dat)[1]
  Cindex=1-Cindex
  Cindex 
}

大家注意一下,这个函数有两个重要的地方,dat <- data[indices,]这句话是一定要有的,因为boot函数就是靠它来重抽样的,第二就是在function中cox回归不能直接写出来,大家看看我是怎么表达出来的。写好函数我们来调试一下,给它取个子集,也就是重抽样1次

c_index(bc,1:100)

在这里插入图片描述
调试成功后就可以重抽样了,

results <- boot(data=bc, statistic=c_index, R=500)

在这里插入图片描述
导出结果,后面的是标准误

print(results) 

在这里插入图片描述
如果你想查看每个抽样的结果

results$t

在这里插入图片描述
可以查看抽样分布,可以看出,我们的抽样分布还是不错的

plot(results)

在这里插入图片描述
计算可信区间

boot.ci(results,conf = 0.95)

在这里插入图片描述
OK,介绍完毕,觉得有用的话多多分享哟。

  • 10
    点赞
  • 79
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 23
    评论
R语言中,我们可以使用`survcomp`包中的`pairwise_survdiff()`函数对两两之间的C-INDEX值进行比较。 首先,我们需要用`coxph()`函数拟合每个模型,并使用`survC1()`函数计算每个模型的C-INDEX值。 ``` r library(survcomp) library(survival) # 导入数据 data <- read.csv("D:/5放射诊断/R生存分析/nafld.csv") # 将CACSgrades和CADRADS转换为因子型变量 data$CACSgrades <- as.factor(data$CACSgrades) data$CADRADS <- as.factor(data$CADRADS) # 拟合模型 model1 <- coxph(Surv(time, MACE) ~ sex + age + NAFLD + Hypertension + Diabetes + Dyslipidemia + smoking, data = data) model2 <- coxph(Surv(time, MACE) ~ sex + age + NAFLD + Hypertension + Diabetes + Dyslipidemia + smoking + CACSgrades, data = data) model3 <- coxph(Surv(time, MACE) ~ sex + age + NAFLD + Hypertension + Diabetes + Dyslipidemia + smoking + CACSgrades + CADRADS, data = data) # 计算C-INDEX值 cindex1 <- survC1(model1, data, time = "time", event = "MACE") cindex2 <- survC1(model2, data, time = "time", event = "MACE") cindex3 <- survC1(model3, data, time = "time", event = "MACE") ``` 接下来,我们可以使用`pairwise_survdiff()`函数对两两之间的C-INDEX值进行比较,以检查它们之间是否有显著差异。 ``` r # 对两两之间的C-INDEX值进行比较 pairwise_survdiff(list("Model 1 vs. Model 2" = cindex1$cindex, "Model 1 vs. Model 3" = cindex2$cindex, "Model 2 vs. Model 3" = cindex3$cindex)) ``` 这将输出一个表格,显示每对模型之间的p值和显著性水平。如果p值小于0.05,则意味着这两个模型的C-INDEX值之间有显著差异。
评论 23
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

天桥下的卖艺者

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值