R语言临床预测模型系列，从入门到进阶，目前已更新20 篇优秀推文，所有数据和代码均可免费获取，持续更新中，欢迎订阅。
区分度评价：C-statistic的计算
净重新分类指数NRI计算
综合判别改善指数IDI的计算
二分类资料校准曲线的绘制
生存资料校准曲线的绘制
tidymodels不能画校准曲线？
mlr3的校准曲线也是一样画！
......
今天学习C-index的比较。
二分类资料C-index的比较
二分类资料的AUC和C-index是一样的，可以参考AUC的比较：ROC（AUC）曲线的显著性检验。
生存资料C-index的比较
可以使用compareC包，专门用来比较生存资料的C-index。
rm(list = ls())
library(compareC)
还是用之前推文的数据，获取示例数据可以查看历史推文。
代码语言：javascript复制load(file = "../000files/timeROC.RData")
str(df2)
## 'data.frame': 297 obs. of 8 variables:
## $ event : num 0 0 1 0 0 1 0 0 0 0 ...
## $ age : int 59 63 65 73 59 66 56 42 61 48 ...
## $ riskScore: num -0.249 -0.511 -0.211 -0.427 0.279 ...
## $ futime : num 3.03 1.16 1.82 1.52 1.34 ...
## $ gender : num 2 2 2 1 2 2 1 2 2 2 ...
## $ t : num 4 4 4 3 3 3 5 3 NA 4 ...
## $ n : num 1 5 1 1 1 1 3 1 NA 1 ...
## $ m : num 1 1 1 1 1 3 1 1 3 3 ...
只要提供4个参数：time，status，第一个指标，第二个指标，即可。
代码语言：javascript复制compareC(df2$futime,
df2$event,
df2$riskScore,
df2$age
)
## $est.c
## Cxy Cxz
## 0.3383690 0.3894508
##
## $est.diff_c
## [1] -0.05108181
##
## $est.vardiff_c
## [1] 0.002124315
##
## $est.varCxy
## [1] 0.001046384
##
## $est.varCxz
## [1] 0.00112578
##
## $est.cov
## [1] 2.392402e-05
##
## $zscore
## [1] -1.108299
##
## $pval
## [1] 0.2677329
第1行给出了2个C指数；第2行是2个C指数相减的差值；第3行是2个C指数方差相减的差值；第4、5行是方差；第6行是协方差；第7、8行是z值和p值。
两个cox模型的比较
下面再多说一点两个cox模型的比较，也是很简单的，方差分析即可，使用anova()函数。
我们用lung数据集进行演示。
library(survival)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
df1 <- lung %>%
mutate(status=ifelse(status == 1,1,0))
建立两个cox模型：
代码语言：javascript复制cox_fit1 <- coxph(Surv(time, status) ~ age sex ph.ecog ph.karno pat.karno,
data = lung,x = T, y = T)
cox_fit2 <- coxph(Surv(time, status) ~ ph.ecog ph.karno pat.karno,
data = lung,x = T, y = T)
直接使用anova()即可：
anova(cox_fit1,cox_fit2)
## Analysis of Deviance Table
## Cox model: response is Surv(time, status)
## Model 1: ~ age sex ph.ecog ph.karno pat.karno
## Model 2: ~ ph.ecog ph.karno pat.karno
## loglik Chisq Df P(>|Chi|)
## 1 -706.48
## 2 -712.83 12.703 2 0.001745 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rms包也带了一个lrtest()函数，也是直接构建两个cph模型即可进行比较，简单演示一下。
suppressMessages(library(rms))
cox_fit1 <- cph(Surv(time, status) ~ age sex ph.ecog ph.karno pat.karno,
data = lung,x = T, y = T)
cox_fit2 <- cph(Surv(time, status) ~ ph.ecog ph.karno pat.karno,
data = lung,x = T, y = T)
lrtest(cox_fit1,cox_fit2)
##
## Model 1: Surv(time, status) ~ age sex ph.ecog ph.karno pat.karno
## Model 2: Surv(time, status) ~ ph.ecog ph.karno pat.karno
##
## L.R. Chisq d.f. P
## 12.702485990 2.000000000 0.001744577
可以看到P值是一样的~
