Literotica
发布日期：2025-12-17 14:15 点击次数：112💡专注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绘制的 Literotica，就不再叠加了。
本站仅提供存储行状，所有执行均由用户发布，如发现存害或侵权执行，请点击举报。