小程序开发价格 xgboost生计分析和加快失效模子

发布日期：2024-09-28 13:29    点击次数：70

💡专注R说念话在🩺生物医学中的使用小程序开发价格

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

对于进步算法咱们也曾先容过了GBM模子、Xgboost、lightGBM（Catboost暂不可用）3种，这几种次序皆是守旧生计分析的，咱们之前的推文皆是以分类梗概纪念为例进行先容的。现在对于这几种算法进行生计分析的官方文档也很少（皆是基于Python的…）。

是以咱们今天先容下如安在R说念话中使用xgboost完了生计分析和加快失效模子。

加载R包和数据

library(xgboost)library(survival)

用1示意牺牲，0示意删失：

app

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

展望变量的准备莫得什么格外的：

2024年有五项世界大赛开战，再加上上半年进行决赛的梦百合杯，本赛季的六项世界大赛，已经有三项有了决赛人选。梦百合杯李轩豪胜党毅飞，衢州烂柯杯辜梓豪对垒申真谞，小程序开发专业公司应氏杯谢科迎战一力辽。中国棋手占据了其中四位，中国围棋的“厚度”优势依旧。世界大赛四强八强的人数和人次也能佐证这一点。

西班牙vs法国的半决赛中，上半场第9分钟，姆巴佩吸引吸收后传中，穆阿尼后点包抄头球破门，法国队取得本届杯赛的第一个运动战进球。第21分钟，亚马尔一记漂亮的世界波帮助西班牙扳平比分，他以16岁362天的年龄，成为欧洲杯历史上最年轻的进球者。4分钟后，奥尔默在禁区内大力抽射破门，西班牙2-1逆转，并且将比分保持到终场，挺进决赛。

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、简书同名账号哔哩哔哩：阿越即是我 本站仅提供存储劳动，统共本色均由用户发布，如发现存害或侵权本色，请点击举报。

• 小程序开发价格
• 分析
• 加快
• 生计
• XGBoost