5
$\begingroup$

The xgboost package enables survival modeling using parameter arguments: objective = "survival:cox" and eval_metric = "cox-nloglik".

The predict method for the resulting model only outputs risk scores (same as type = "risk" in the survival::coxph function in r).

How do I use xgboost to predict entire survival curves?

$\endgroup$
0

2 Answers 2

11
$\begingroup$

The proportional hazard model assumes hazard rates of the form: $h(t|X) = h_0(t) \cdot risk(X)$ where usually $risk(X) = exp(X\beta)$. The xgboost predict method returns $risk(X)$ only. What we can do is use the survival::basehaz function to find $h_0(t)$.

Problem is it's not "calibrated" to the actual baseline hazard rate computed in xgboost. What we can do is find some constant $C$ that minimizes the ibrier score between the sample observed death/censorship times and $h_0(t) \cdot risk(X) \cdot C$.

I've implemented this approach in a tiny R package I've written.

$\endgroup$
4
  • $\begingroup$ How to calculate C-index in this setting? Can your package accomplish this? Thanks! $\endgroup$
    – Tommy
    Mar 31, 2022 at 14:07
  • $\begingroup$ My package does not include that functionality. I would recommend using either the pec or riskRegression packages $\endgroup$
    – Iyar Lin
    Apr 1, 2022 at 17:34
  • $\begingroup$ How do censorship times contribute to the ibrier score? $\endgroup$
    – 42-
    Jun 27, 2022 at 22:46
  • $\begingroup$ Not sure I understand the question. The ibrier score is a function of the estimated and observed survival functions, in which censorship times also take a role. $\endgroup$
    – Iyar Lin
    Jun 29, 2022 at 10:13
3
$\begingroup$

The solution to use survival::basehaz() with a coxph model and estimate a constant C, as implemented by survXgboost should be used with caution. When you have binary predictors, coxph coefficients explode, leading to really overestimated baseline hazard, the constant C will not do much and the performance of xgboost will look much worse than what it really is.

The gbm package has a function gbm::basehaz which skips the model, avoiding the compatibility problem that you have in survival::basehaz(), and uses the predict() results to estimate the baseline hazard. It is more reliable and the (cumulative) baseline hazard is as expected.

$\endgroup$
2
  • 1
    $\begingroup$ Note that gbm::basehaz returns the Breslow estimator. This can be done using survival::survfit and specifying stype = 2 as well, which skips the dependency on gbm (as most R users working with survival data will use the survival package). $\endgroup$
    – Erik A
    Jan 31 at 11:51
  • $\begingroup$ I've used the advice from this thread and it's now incorporated to the code. Thanks @ErikA for your contribution. $\endgroup$
    – Iyar Lin
    Feb 1 at 18:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.