Cox, AFT, and Piecewise-Exponential Objectives

library(fastgbm)
library(survival)

Cox objective (objective = "cox", the default)

Minimizes the negative Cox partial log-likelihood (Breslow handling of tied event times). predict(..., type = "link") returns the log relative-risk score; higher values mean higher risk (shorter expected survival). predict(..., type = "survival") combines the fitted risk score with a Breslow baseline cumulative hazard estimated after training to return survival probabilities at specific horizons: \(\hat S(t \mid x) = \hat S_0(t)^{\exp(f(x))}\).

lung_dat <- na.omit(lung[, c("time", "status", "age", "sex", "ph.ecog")])
x <- as.matrix(lung_dat[, c("age", "sex", "ph.ecog")])

cox_fit <- fastgbm(x, time = lung_dat$time, status = lung_dat$status,
                   objective = "cox", ntrees = 100L, max_depth = 3L, verbose = FALSE)

risk <- predict(cox_fit, x, type = "link")
surv <- predict(cox_fit, x[1:5, ], type = "survival", times = c(90, 180, 365))
surv
#>           [,1]      [,2]      [,3]
#> [1,] 0.9522407 0.8707755 0.6382920
#> [2,] 0.9072084 0.7593050 0.4092590
#> [3,] 0.9582282 0.8863460 0.6760728
#> [4,] 0.9059485 0.7563272 0.4040743
#> [5,] 0.9547701 0.8773314 0.6540165

AFT objective (objective = "aft")

Fits a normal location-scale accelerated failure time model: \(\log T = f(x) + \sigma Z\), \(Z \sim N(0, 1)\). Right-censored observations contribute via the survival function of the normal distribution. Unlike Cox, the linear predictor estimates \(\log(\text{time})\) directly (higher = longer predicted survival), the opposite direction from a Cox risk score – metrics() accounts for this automatically when computing concordance, but if you use the raw linear predictor for anything else, remember to negate it for a risk-like ordering.

set.seed(1)
n <- 200
x_aft <- matrix(rnorm(n * 2), ncol = 2, dimnames = list(NULL, c("x1", "x2")))
time <- exp(1 + 0.6 * x_aft[, 1] - 0.3 * x_aft[, 2] + rnorm(n, sd = 0.3))
status <- rbinom(n, 1, 0.75)

aft_fit <- fastgbm(x_aft, time = time, status = status, objective = "aft",
                   ntrees = 100L, max_depth = 3L, verbose = FALSE)
predict(aft_fit, x_aft[1:5, ], type = "survival", times = c(1, 3, 5))
#>           [,1]      [,2]      [,3]
#> [1,] 0.6969007 0.3565512 0.8543226
#> [2,] 0.2799113 0.8553721 0.7068923
#> [3,] 0.1369927 0.4845023 0.9835060
#> [4,] 0.8928347 0.2912686 0.8493422
#> [5,] 0.5569074 0.9843708 0.6994400

Piecewise-exponential objective (objective = "pexp")

Unlike Cox and AFT, pexp doesn’t fit a single scalar score per subject; it models the hazard jointly over covariates and time. Training data is expanded into person-time rows via the standard “Poisson trick” (one row per subject per time interval they were at risk in), and the ensemble predicts the log hazard rate for each row – an ordinary Poisson-with-offset objective, no baseline assumption needed. This means predict(..., type = "link") isn’t a single “risk score” the way Cox’s is; it defaults to the cumulative hazard at the model’s full fitted time horizon, a fixed and well-defined (if coarser) risk score for ranking. predict(..., type = "survival") evaluates the fitted hazard-over-time surface directly at the requested times.

pexp_fit <- fastgbm(x, time = lung_dat$time, status = lung_dat$status,
                    objective = "pexp", ntrees = 100L, max_depth = 3L, verbose = FALSE)
predict(pexp_fit, x[1:5, ], type = "survival", times = c(90, 180, 365))
#>           [,1]      [,2]      [,3]
#> [1,] 0.9689721 0.8977238 0.7209522
#> [2,] 0.9608723 0.8953829 0.7045686
#> [3,] 0.9681084 0.9153831 0.7564132
#> [4,] 0.9608207 0.8684158 0.6675388
#> [5,] 0.9692716 0.9183865 0.7641165
metrics(pexp_fit, y = Surv(lung_dat$time, lung_dat$status))
#> $objective
#> [1] "pexp"
#> 
#> $metric
#> [1] "cindex"
#> 
#> $value
#> [1] 0.6226105

Concordance

metrics(cox_fit, y = Surv(lung_dat$time, lung_dat$status))
#> $objective
#> [1] "cox"
#> 
#> $metric
#> [1] "cindex"
#> 
#> $value
#> [1] 0.7033984
metrics(aft_fit, y = Surv(time, status))
#> $objective
#> [1] "aft"
#> 
#> $metric
#> [1] "cindex"
#> 
#> $value
#> [1] 0.9107155

Missing values

Missing predictor values are routed natively: at every candidate split, both “missing goes left” and “missing goes right” are evaluated, and the direction that reduces loss more is stored on the node. No imputation is performed unless you do it yourself.

x_missing <- x
x_missing[1, "ph.ecog"] <- NA
predict(cox_fit, x_missing[1:3, ], type = "link")
#> [1] -0.0681736  0.2065150 -0.6186579