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

设为“星标”，精彩不错过

关于提升算法我们已经介绍过了GBM模型、Xgboost、lightGBM（Catboost暂不可用）3种，这几种方法都是支持生存分析的，我们之前的推文都是以分类或者回归为例进行介绍的。目前关于这几种算法进行生存分析的官方文档也很少（都是基于Python的…）。

所以我们今天介绍下如何在R语言中使用xgboost实现生存分析和加速失效模型。

加载R包和数据

library(xgboost)library(survival)

用1表示死亡，0表示删失：

lung$status <- ifelse(lung$status == 2,1,0)lung <- na.omit(lung) # 去掉NA

简单的划分下数据，按照7：3划分训练集和测试集：

set.seed(123)ind <- sample(1:nrow(lung),nrow(lung)*0.7)train_df <- lung[ind,]test_df <- lung[- ind, ]dim(train_df)## [1] 116  10dim(test_df)## [1] 51 10

生存分析

下面就是设置下xgboost的参数，注意生存分析objective = "survival:cox"，以及eval_metric = "cox-nloglik"。

数据格式要设置为xgboost的专用格式，生存分析要注意因变量的设置。直接用Surv是不行的，xgboost中会把生存时间为负值的当做是删失，所以我们设置y时要如下设置：

# 选择参数的值param <- list(objective = "survival:cox",              booster = "gbtree",              eval_metric = "cox-nloglik",              eta = 0.03,              max_depth = 3,              subsample = 1,              colsample_bytree = 1,              gamma = 0.5)# 准备预测变量和结果变量x <- as.matrix(train_df[, c("age","sex","ph.ecog","ph.karno","pat.karno")])# 设置y，把删失的改成负值y <- ifelse(train_df$status == 1, train_df$time, -train_df$time)# 放进专用的格式中train.mat <- xgb.DMatrix(data = x, label = y)train.mat## xgb.DMatrix  dim: 116 x 5  info: label  colnames: yes

然后就可以进行生存分析了：

set.seed(1)xgb.fit <- xgb.train(params = param,                      data = train.mat,                      nrounds = 1000,                     watchlist = list(val2 = train.mat),                     early_stopping_rounds = 50)## [1]  val2-cox-nloglik:3.869080 ## Will train until val2_cox_nloglik hasn't improved in 50 rounds.## ## [2]  val2-cox-nloglik:3.856530 ## [3]  val2-cox-nloglik:3.844888 ## [4]  val2-cox-nloglik:3.834070 ## 省略## [222]    val2-cox-nloglik:3.367639 ## [223]    val2-cox-nloglik:3.367639 ## [224]    val2-cox-nloglik:3.367639 ## Stopping. Best iteration:## [174]    val2-cox-nloglik:3.367639xgb.fit## ##### xgb.Booster## raw: 238.9 Kb ## call:##   xgb.train(params = param, data = train.mat, nrounds = 1000, watchlist = list(val2 = train.mat), ##     early_stopping_rounds = 50)## params (as set within xgb.train):##   objective = "survival:cox", booster = "gbtree", eval_metric = "cox-nloglik", eta = "0.03", max_depth = "3", subsample = "1", colsample_bytree = "1", gamma = "0.5", validate_parameters = "TRUE"## xgb.attributes:##   best_iteration, best_msg, best_ntreelimit, best_score, niter## callbacks:##   cb.print.evaluation(period = print_every_n)##   cb.evaluation.log()##   cb.early.stop(stopping_rounds = early_stopping_rounds, maximize = maximize, ##     verbose = verbose)## # of features: 5 ## niter: 224## best_iteration : 174 ## best_ntreelimit : 174 ## best_score : 3.367639 ## best_msg : [174] val2-cox-nloglik:3.367639 ## nfeatures : 5 ## evaluation_log:##     iter val2_cox_nloglik##        1         3.869080##        2         3.856530## ---                      ##      223         3.367639##      224         3.367639

用在测试集上，计算风险分数：

test_x <- as.matrix(test_df[, c("age","sex","ph.ecog","ph.karno","pat.karno")])risk_score <- predict(xgb.fit, newdata = test_x) # newdata如果是训练集可以获取训练集的风险分数risk_score##  [1] 0.4721478 1.8668407 3.4487274 0.5234923 0.5299142 0.2525677 0.3887208##  [8] 0.8716910 5.8021851 0.2647004 0.6711149 2.4503310 0.4688006 1.4702020## [15] 0.7774669 0.3813822 0.9095489 0.6189218 0.3667074 0.3612074 0.4357747## [22] 0.9482961 0.4830345 0.6624518 0.4522163 0.3241239 0.5598004 1.4702020## [29] 2.3615212 2.0817354 0.6228051 1.8675944 0.5460874 1.4702020 0.5198651## [36] 0.2530192 9.6702890 0.4650362 0.5829992 0.9305520 0.2432446 0.4080229## [43] 0.4403237 0.2779934 0.7739328 0.2030164 0.2839084 0.3376622 0.4515978## [50] 0.3893496 0.2779934hist(risk_score)

图片

image-20240204181610457

简单做个生存分析看看，根据风险分数高低分组，然后做K-M生存曲线：

groups <- ifelse(risk_score>median(risk_score),"high","low")library(survminer)# 根据rs_train_group建立生存函数f_train <- survfit(Surv(test_df$time, test_df$status) ~ groups)f_train## Call: survfit(formula = Surv(test_df$time, test_df$status) ~ groups)## ##              n events median 0.95LCL 0.95UCL## groups=high 25     20    291     223     460## groups=low  26     20    426     348     643ggsurvplot(f_train,           data = test_df,           surv.median.line = "hv",            #legend.title = "Risk Group",           #legend.labs = c("Low Risk", "High Risk"),           pval = TRUE,            ggtheme = theme_bw()           )

图片

image-20240204181624495

这个结果是没有意义的，但是可以看到在中间的部分，高风险组的死亡率是明显高于低风险组的。

计算C-index，借助Hmisc包实现的：

library(Hmisc)rcorr.cens(as.numeric(risk_score),           Surv(test_df$time, test_df$status))##        C Index            Dxy           S.D.              n        missing ##      0.4178606     -0.1642789      0.1153034     51.0000000      0.0000000 ##     uncensored Relevant Pairs     Concordant      Uncertain ##     40.0000000   2094.0000000    875.0000000    450.0000000

brier score怎么计算？直接改用mlr3吧，满足你对机器学习生存分析的所有幻想→2024年了，高阶生存分析我强烈推荐你用mlr3。

但是要注意，mlr3只是一个框架，它并不能解决你所有问题，很多时候你可能还是要回到这种基础R包中来。

加速失效模型

在R语言中通过xgboost实现加速失效模型挺复杂的，主要是要把数据准备成它需要的格式。

其中的label（也就是因变量，生存时间）需要用一个区间表示，具体需要的形式如下所示：

图片

rm(list = ls())library(xgboost)lung$status <- ifelse(lung$status == 2,1,0)lung <- na.omit(lung) # 去掉NA

预测变量的准备没有什么特别的：

X <- as.matrix(lung[, c("age","sex","ph.ecog","ph.karno","pat.karno")])dtrain <- xgb.DMatrix(X)

结果变量的准备很复杂，主要是根据上面的表格。如果不是删失数据，那么生存时间就是闭区间[a,a]，如果是右删失就是闭区间[a,正无穷]，如果是左删失就是[0,b]，区间删失就是[a,b]。

下面我展示下未删失数据和右删失数据的设置方法，还是用lung这个数据集。如果status是1，就是未删失数据，那么它的上下限就都是生存时间；如果status是0，那就是右删失数据，那么它的上下限就是[生存时间,正无穷]。

# 设置未删失和右删失数据的上下区间：y_lower_bound <- ifelse(lung$status == 1, lung$time, Inf)y_upper_bound <- ifelse(lung$status == 1, lung$time, Inf)

设置好之后使用setinfo函数把这个信息添加到dtrain这个对象中：

setinfo(dtrain, 'label_lower_bound', y_lower_bound)## [1] TRUEsetinfo(dtrain, 'label_upper_bound', y_upper_bound)## [1] TRUE

这样最难的一步就准备好了。接下来就是设置下参数，就可以训练了。

params <- list(objective='survival:aft', # 加速失效模型               eval_metric='aft-nloglik', # 评价指标，不要乱选               aft_loss_distribution='normal',               aft_loss_distribution_scale=1.20,               tree_method='hist',               learning_rate=0.05,               max_depth=2)watchlist <- list(train = dtrain)bst <- xgb.train(params, dtrain, nrounds=30, watchlist)## [1]  train-aft-nloglik:20.251825 ## [2]  train-aft-nloglik:19.003819 ## [3]  train-aft-nloglik:17.959822 ## [4]  train-aft-nloglik:17.048453 ## 省略## [29] train-aft-nloglik:14.974870 ## [30] train-aft-nloglik:15.246215

可以看到这样就可以了~后面的分析就又都是重复的了，我就不详细介绍了。

参考资料https://github.com/dmlc/xgboost/issues/9979https://xgboost.readthedocs.io/en/stable/tutorials/aft_survival_analysis.html

联系我们，关注我们

免费QQ交流群1：613637742免费QQ交流群2：608720452公众号消息界面关于作者获取联系方式知乎、CSDN、简书同名账号哔哩哔哩：阿越就是我 本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。