Time dependent survival

I would like to see time-dependent covariates added for survival data. I have many clinical research studies that could use this.

The key idea that makes this possible is that for a Cox model, the first derivative can be written as the matrix product mX, where m is the vector of martingale residuals and X is the n x p matrix of predictors. Here n= number of observations in the data set and p = number of predictors. (The second derivative for a Cox model is not nearly as nice -- this makes gradient boosting with it's reliance on first derivatives a nice match. ) This should fit in well with the sparse matrix forms for X that xgboost has adopted.

The data set for time-dependent outcomes needs to be in counting process form. That is, each subject is represented by multiple rows with time1, time2, status, strata, covariates. The covariate values are those that hold true during the interval from time1 to time2, status = 1 if this interval ended with an event and 0 otherwise. Strata is a special covariate, optional, that divides the subjects into disjoint groups; it opens a lot more doors for special analytical models. This form of data is now very common in software packages such as R, SAS, Stata; so users have a lot of resources and tools for creating one. Each subject can have a unique set of intervals and time points.

A key step in all this is the routine that calculates the martingale residuals, over and over and over again. Code for this is integrated into the R survival package. Input arguments for the underlying C routine are the time1, time2, strata, and status vectors, along with the current vector of predicted values for each row of data, and two order vectors: in R they are order(strata, time1) and order(strata, time2). The routine's run time is O(2n), it has to walk down the two order vectors and keep running sums. I have 20 year's experience tuning and bomb proofing this code, some of the latter perhaps not as necessary with a gbm predictor since it won't have wild outliers for the predicted values. (Users come up with the strangest data sets.)

The naive way to do this computation is O(nd) by the way, where d= number of deaths. This is the reason that the run time for SAS phreg can be >50x longer (or more) than R:coxph for large time dependent data sets. All Cox model routines know how to do it in O(n) when there are no time dependent covariates.

I have not thought through whether further speedups are possible when the linear predictor changes only a little from one computation of m to another. This would be a very interesting think and discussion. I'd be happy to discuss the algorithm and hear ideas to make it better; it's easiest to use a picture though.

What I clearly don't know is how best to set this up in xgboost from a user's point of view. The 'response' is now 3 columns. The strata is neither fish nor fowl; not a response and not quite a predictor. I would need to cooperate with someone who knows xgboost much better.
Preprocessing for the code is fairly simple: a data check that time1 < time2 for all rows, that status = 0 or 1, that strata (if present) does not put every obs into its own group (m will be all zeros in this case, and it is almost certainly a user error), and creation of the two order vectors.

Terry Therneau
author of R survival

The issue names no repository files or tests. Start by reviewing xgboost's current survival/Cox support and the R survival package's martingale-residual routine. Done would require an agreed interface for time1, time2, status, and strata, input validation, efficient residual computation, and tests.

Written by the indexing model from the issue text.