一文处分：COX模子的校准弧线（2种法式+进修集+外部考据集）

发布日期：2025-12-24 10:13    点击次数：146

💡专注R话语在🩺生物医学中的使用

设为“星标”，精彩可以过

前边咱们照旧讲过逻辑转头模子的校准弧线的画法，此次咱们学习cox转头的校准弧线画法。

准备数据

使用自带数据lung数据集进行演示。

这个是对于肺癌的生涯数据，一共有228行，10列，其中time是生涯本领，单元是天，status是生涯状态，1是删失，2是物化。其余变量是自变量，道理如下：

inst：机构代码，对于咱们此次建模没啥用age：年纪sex：1是男性，2是女性ph.ecog：ECOG评分。0=无症状，1=有症状但完好意思可以来往，2=每天<50%的本领在床上，3=在床上>50%的本领但莫得卧床，4=卧床不起ph.karno：医师评的KPS评分，领域是0-100，得分越高，健康情景越好，越能隐忍调养给身体带来的反作用。pat.karno：患者我方评的KPS评分meal.cal：用餐时摧毁的卡路里wt.loss：曩昔6个月的体重减小数，单元是磅

library(survival)rm(list = ls())dim(lung)

[1] 228  10

str(lung)

'data.frame':   228 obs. of  10 variables: $ inst     : num  3 3 3 5 1 12 7 11 1 7 ... $ time     : num  306 455 1010 210 883 ... $ status   : num  2 2 1 2 2 1 2 2 2 2 ... $ age      : num  74 68 56 57 60 74 68 71 53 61 ... $ sex      : num  1 1 1 1 1 1 2 2 1 1 ... $ ph.ecog  : num  1 0 0 1 0 1 2 2 1 2 ... $ ph.karno : num  90 90 90 90 100 50 70 60 70 70 ... $ pat.karno: num  100 90 90 60 90 80 60 80 80 70 ... $ meal.cal : num  1175 1225 NA 1150 NA ... $ wt.loss  : num  NA 15 15 11 0 0 10 1 16 34 ...

大多量情况下王人是使用1代表物化，0代表删失，这个数据集用2代表物化，但有的R包会报错，需要详确！

咱们这里给它悔改来：

library(dplyr)library(tidyr)lung <- lung %>%   mutate(status=ifelse(status == 2,1,0))str(lung)

'data.frame':   228 obs. of  10 variables: $ inst     : num  3 3 3 5 1 12 7 11 1 7 ... $ time     : num  306 455 1010 210 883 ... $ status   : num  1 1 0 1 1 0 1 1 1 1 ... $ age      : num  74 68 56 57 60 74 68 71 53 61 ... $ sex      : num  1 1 1 1 1 1 2 2 1 1 ... $ ph.ecog  : num  1 0 0 1 0 1 2 2 1 2 ... $ ph.karno : num  90 90 90 90 100 50 70 60 70 70 ... $ pat.karno: num  100 90 90 60 90 80 60 80 80 70 ... $ meal.cal : num  1175 1225 NA 1150 NA ... $ wt.loss  : num  NA 15 15 11 0 0 10 1 16 34 ...

然后把数据分辩为进修集、测试集，底下这个分辩法式是错的哈，这个示例数据样本量太少了，我思让进修集和测试集样本量王人尽量多一丝，是以用了一个失实的法式进行演示：

set.seed(123)ind1 <- sample(1:nrow(lung),nrow(lung)*0.7)train_df <- lung[ind1,]set.seed(563435)ind2 <- sample(1:nrow(lung),nrow(lung)*0.7)test_df <- lung[ind2, ]dim(train_df)

[1] 159  10

dim(test_df)

[1] 159  10

法式1：rms

library(rms)# 必须先打包数据dd <- datadist(train_df)options(datadist = "dd")units(train_df$time) <- "day" # 单元成立为：天

构建cox比例风险模子。rms可以使用里面重抽样的法式绘制校准弧线，可以选拔bootstrap法粗略交叉考据法，底下咱们选拔500次bootstrap的里面考据法式，计较本领点为第100天的校准弧线：

# 成立cox转头模子，本领点选拔1年coxfit1 <- cph(Surv(time, status) ~ sex + ph.ecog + ph.karno + wt.loss,              data = train_df, x = T, y = T, surv = T,              time.inc = 100 # 100天              )# m=40暗示以40个样本为1组，一般取4-6组，咱们这个数据样本量太少了# u=100和上头的time.inc对应cal1 <- calibrate(coxfit1, cmethod = "KM", method = "boot",                  u = 100, m = 40, B = 500)     # 由于样本量太少出现了教学：    Using Cox survival estimates at  100 days        Warning in groupkm(cox, Surv(y[, 1], y[, 2]), u = u, cuts = orig.cuts): one    interval had < 2 observations        Warning in groupkm(cox, Surv(y[, 1], y[, 2]), u = u, cuts = orig.cuts): one    interval had < 2 observations        Warning in groupkm(cox, Surv(y[, 1], y[, 2]), u = u, cuts = orig.cuts): one    interval had < 2 observations        Warning in groupkm(cox, Surv(y[, 1], y[, 2]), u = u, cuts = orig.cuts): one    interval had < 2 observations

进修集

接下来等于绘画，可以径直使用plot()函数：

plot(cal1)

图片

也可以提真金不怕火数据，我方画，以竣事更多的细节适度：

plot(cal1,     lwd = 2, # 短处线粗细     lty = 1, # 短处线类型，可选0-6     errbar.col = c("#2166AC"), # 短处线阵势     xlim = c(0.7,1),ylim= c(0.7,1), # 坐标轴领域     xlab = "Prediced OS",ylab = "Observed OS",     cex.lab=1.2, cex.axis=1, cex.main=1.2, cex.sub=0.6) # 字体大小lines(cal1[,c('mean.predicted',"KM")],       type = 'b', # 连线的类型，可以是"p","b","o"      lwd = 3, # 连线的粗细      pch = 16, # 点的体式，可以是0-20      col = "tomato") # 连线的阵势box(lwd = 2) # 边框粗细abline(0,1,lty = 3, # 对角线为虚线       lwd = 2, # 对角线的粗细       col = "grey70" # 对角线的阵势       ) 

图片

测试集

然后是外部考据集的校准弧线：

# 赢得测试集的揣测的生涯概率，这一步有莫得王人行estimates <- survest(coxfit1,newdata=test_df,times=100)$survvs <- val.surv(coxfit1, newdata = test_df,               S = Surv(test_df$time,test_df$status),               est.surv = estimates,# 这一步有莫得王人行               u = 100 # 本领点，亦然选100天               )plot(vs)

图片

法式2：riskRegression

这个R包也极端好用，然而这种法式弗成有缺失值，是以我先把缺失值去掉，然后再分辩进修集和测试集，然而由于样本量太少，这里的分辩法式是不正确的哈。

# 删除缺失值df2 <- na.omit(lung)# 分辩数据set.seed(123)ind1 <- sample(1:nrow(df2),nrow(df2)*0.9)train_df <- df2[ind1,]set.seed(563435)ind2 <- sample(1:nrow(df2),nrow(df2)*0.9)test_df <- df2[ind2, ]dim(train_df)

[1] 150  10

dim(test_df)

[1] 150  10

# 构建模子cox_fit1 <- coxph(Surv(time, status) ~ sex + ph.ecog + ph.karno,                  data = train_df,x = T, y = T)

进修集

# 绘画library(riskRegression)set.seed(1)cox_fit_s <- Score(list("fit1" = cox_fit1),               formula = Surv(time, status) ~ 1,               data = train_df,               plots = "calibration",               conf.int = T,               B = 500,               M = 50, # 每组的东谈主数               times=c(100) # 本领点选100天               )plotCalibration(cox_fit_s,cens.method="local",# 减少输出日记                xlab = "Predicted Risk",                ylab = "Observerd RISK")

图片

诚然亦然可以用ggplot2绘画的。

# 赢得数据data_all <- plotCalibration(cox_fit_s,plot = F,cens.method="local")# 数据调理plot_df <- data_all$plotFrames$fit1# 绘画library(ggplot2)ggplot(plot_df, aes(Pred,Obs))+  geom_point()+  geom_line(linewidth=1.2)+  scale_x_continuous(limits = c(0,0.5),name = "Predicted Risk")+  scale_y_continuous(limits = c(0,0.5),name = "Observerd Risk")+  geom_abline(slope = 1,intercept = 0,lty=2)+  theme_bw()

图片

测试集

使用起来完好意思相似，只需要提供测试集即可：

set.seed(1)cox_fit_s <- Score(list("fit1" = cox_fit1),               formula = Surv(time, status) ~ 1,               data = test_df, # 测试集               plots = "calibration",               B = 500,               M = 50,               times=c(100) # 本领点               )plotCalibration(cox_fit_s,cens.method="local",# 减少输出日记                xlab = "Predicted Risk",                ylab = "Observerd Risk")

图片

这个效果亦然可以用ggplot2绘制的，就不再重叠了。

本站仅提供存储就业，悉数试验均由用户发布，如发现存害或侵权试验，请点击举报。