Introduction to the MTLR Workflow

Data: lung

One can access the lung dataset by loading the survival package.

library(survival)
#Looking at the top 6 rows...
head(lung)
#>   inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
#> 1    3  306      2  74   1       1       90       100     1175      NA
#> 2    3  455      2  68   1       0       90        90     1225      15
#> 3    3 1010      1  56   1       0       90        90       NA      15
#> 4    5  210      2  57   1       1       90        60     1150      11
#> 5    1  883      2  60   1       0      100        90       NA       0
#> 6   12 1022      1  74   1       1       50        80      513       0
#help(lung) #See the basic information of lung.

If you look at the help file for lung you will see the following feature definitions:

• inst: Institution code
• time: Survival time in days
• status: censoring status 1=censored, 2=dead
• age: Age in years
• sex: Male=1 Female=2
• ph.ecog: ECOG performance score (0=good 5=dead)
• ph.karno: Karnofsky performance score (bad=0-good=100) rated by physician
• pat.karno: Karnofsky performance score as rated by patient
• meal.cal: Calories consumed at meals
• wt.loss: Weight loss in last six months

Most importantly you will notice the two features needed for every survival dataset for use of MTLR – an event time (here time), and the indicator identifying if an observation is uncensored/censored (here status). For this example we have status == 1 indicating a right censored individual and status == 2 indicating an uncensored individual. Later on we will be using the Surv function to structure our survival data for MTLR – there are other acceptable formats for the indicator feature (status) – see help(Surv) for more information.

Pre-processing

We will remove inst for this example since this is a categorical feature with 19 unique values and we would like to keep the number of features relatively small.

lung <- lung[,-1]

Before progressing any further we will split our data into a training and testing set. Note that we could stratify our training/testing set by the censor status but for simplicity we skip that for now.

numberTrain <- floor(nrow(lung)*0.8)
set.seed(42)
trInd <- sample(1:nrow(lung), numberTrain)
training <- lung[trInd,]
testing <- lung[-trInd,]

You may also notice that there are some missing values in the data, namely in meal.cal and wt.loss (although ph.ecog, ph.karno, and pat.karno also have missing values). The MTLR package does not handle missing values for users so this must be pre-processed ahead of time. If one passes in data which contains missing values anyway, all rows with missing values will be removed before model training/predictions. To remedy this problem we perform a very basic mean imputation on the dataset. Note that we use the means from the training set to impute the test set.

#Perform imputation
trMeans <- colMeans(training,na.rm=T)
for(i in 1:ncol(training)){
  training[is.na(training[,i]), i] <- trMeans[i]
  testing[is.na(testing[,i]), i] <- trMeans[i]
}

Model Training

Once the dataset has been prepared we can begin to play around with some of the functions found in the MTLR package. Most importantly we will be utilizing the mtlr function to train our model. There are a number of arguments that can be used by mtlr, though only a select few are discussed here. There are only two arguments required to train an mtlr model, formula and data. For formula we must structure our event time feature and censor indicator feature using the Surv function. Since we have time and status as these two features we can create our formula object:

formula <- Surv(time,status)~.

The above says we will be training a model on the survival object created from time and status and using all the other features in our dataset as predictors. If we wanted to select a few features we could do this as well, for example, with age and sex.

formulaSmall <- Surv(time, status)~age+sex

Next, we just need the data argument which in our case is training. We can finally make our first model!

library(MTLR)
fullMod <- mtlr(formula = formula, data = training)
smallMod <- mtlr(formula = formulaSmall, data = training)
#We will print the small model so the output is more compact.
smallMod
#> 
#> Call:  mtlr(formula = formulaSmall, data = training) 
#> 
#> Time points:
#>  [1]  60.6 101.2 155.4 177.0 192.7 210.6 235.4 269.5 291.8 310.6 353.0
#> [12] 386.2 455.1 553.6 688.4
#> 
#> 
#> Weights:
#>           Bias      age     sex
#> 60.62   0.0714  0.05518 -0.0148
#> 101.25  0.0933  0.03607 -0.0275
#> 155.44  0.0863  0.02890 -0.0272
#> 177    -0.0761  0.03339 -0.0465
#> 192.69  0.1011  0.01274 -0.0566
#> 210.62  0.4183  0.01650 -0.0318
#> 235.38 -0.0482  0.00554 -0.0410
#> 269.5  -0.3877 -0.01542 -0.0284
#> 291.81  0.0233  0.01393 -0.0471
#> 310.62 -0.2513  0.02312 -0.0414
#> 353     0.0146  0.01510 -0.0233
#> 386.25 -0.4602 -0.00296 -0.0199
#> 455.12  0.4586  0.01033 -0.0248
#> 553.62 -0.5599  0.01773 -0.0263
#> 688.38 -0.2246  0.02091 -0.0308

There is a lot to take in at first from the output of the mtlr model. The first item is simply the call that was used to build the model. Next is the time points that mtlr used to train the model. If these time points are not specified when constructing the model then mtlr will choose time points based on the quantiles of the event time feature. Additionally, the number of time points is chosen to be the sqrt(N) where N is the number of observations. Since we had 205 training instances and the sqrt(205 = 14.317) mtlr rounded up to 15 time points.

Last, mtlr outputs the weight matrix for the model – these are the weights corresponding to each feature at each time point (additionally notice that we include the bias weights). The row names correspond to the time point for which the feature weight belongs. If you would like to access these weights, they are saved in the model object as weight_matrix so you can access them using smallMod$weight_matrix.

We can also plot the weights for a mtlr model. Before we printed the small model but here we will look at the weights for the complete model.

plot(fullMod)

By default, plot will only look at the 5 features which had the largest sum of absolute values across time (the most influence). You can alter these specifications by playing with the arguments in plot.

Model Predictions

Now that we have trained a MTLR model we should make some predictions! This is where our testing set and the predict function will come into play. Note that there are a number of predictions we may be interested in acquiring. First, we may want to view the survival curves of our test observations.

survCurves <- predict(fullMod, testing, type = "survivalcurve")
#survCurves is pretty large so we will look at the first 5 rows/columns.
survCurves[1:5,1:5]
#>       time         1         2         3         4
#> 1   0.0000 1.0000000 1.0000000 1.0000000 1.0000000
#> 2  60.6250 0.9197721 0.9515804 0.8685207 0.9201473
#> 3 101.2500 0.8427539 0.9092633 0.7588228 0.8429653
#> 4 155.4375 0.7720628 0.8716582 0.6670640 0.7726787
#> 5 177.0000 0.7080682 0.8373624 0.5913061 0.7086728

When we use the predict function for survival curves we will be returned a matrix where the first column (time) is the list of time points that the model evaluated the survival probability for each observation (these will be the time points used by mtlr and an additional 0 point). Every following column will correspond to the row number of the data passed in, e.g. column 2 (named 1) corresponds to row 1 of testing. Each row of this matrix gives the probabilities of survival at the corresponding time point (given by the time column). For example, testing observation 1 has a survival probability of 0.919 at time 60.625.

Since these curves may be hard to digest by observing a matrix of survival probabilities we can also choose to plot them.

plotcurves(survCurves, 1:10)

Here we have specified that we want to observe the survival curves for the first 10 observations (corresponding to the first 10 rows of testing). You will notice that these curves have been smoothed whereas before we only had probabilities for certain time points. We have performed a monotonic spline fit to those survival probabilities to produce the curves you see here.

Additionally, you may have specific plot specifications you want to make. plotcurves is simply returning a ggplot2 object so specifications can be made like you would make to any other ggplot2 graphic. For example, plotcurves(survCurves, 1:10) + ggplot2::xlab("Days") would change the x-axis label to “Days” instead of “Time”.

#Mean
meanSurv <- predict(fullMod, testing, type = "mean_time")
head(meanSurv)
#> [1] 319.0312 407.2113 262.0245 317.3220 327.1613 336.3179
#Median
medianSurv <- predict(fullMod, testing, type = "median_time")
head(medianSurv)
#> [1] 276.5360 378.9466 197.3210 274.5491 279.7247 296.5696

Survival Probability at Event Time

The last prediction type supported is acquiring the observations survival probability at the respective event time. However, in order to use this prediction, the event time (whether censored or uncensored) must be included in the features passed into the predict function.

survivalProbs <- predict(fullMod, testing, type = "prob_event")
head(survivalProbs)
#> [1] 0.58742054 0.15956481 0.06579507 0.56822951 0.19219566 0.80750612
#To see what times these probabilities correspond to:
head(testing$time)
#> [1]  210  883 1022  218  567  144

You will notice that some of these survival probabilities correspond to 0 (usually those with very large event times). We again have drawn the linear extension for the survival time if the event time could not be mapped onto the survival curve.

Miscellaneous Commands

mtlr_cv(formula,training, C1_vec = c(0.01,0.1,1,10,100))
#> $best_C1
#> [1] 1
#> 
#> $avg_loss
#>     0.01      0.1        1       10      100 
#> 2.334196 2.173089 2.127982 2.152218 2.183193

create_folds

As we mentioned, mtlr_cv uses an internal k-fold cross validation to evaluate the loss. We also export the function (create_folds) used to create these cross-validation folds as it is creating folds in a unique way.

These folds can be deterministic, semi-deterministic, or totally random. The deterministic folds arise by stratifying folds by censor status and attempting to create equal ranges in the event times within each fold. This is done by first stratifying the survival dataset into a censored and uncensored portion and then sorting each portion by the event time. These portions are then numbered off into k different folds (see figure below). This option corresponds to “fullstrat”.

The “fullstrat” method of cross-fold validation. Here “delta” refers to the feature used to identify censored/uncensored observations.

The semi-deterministic method stratifies the folds by censor status but ignores event time. In the figure above this would be the same process but skipping the sorting of time. This option is called “censorstrat”. The completely random method ignores the censor status and the event time and randomly assigns observations to folds – the “random” option.

We include the create_folds function in the event one would want to perform k-fold cross-validation on their dataset instead of using a training/testing split. For example, one could set up for 5-fold cross-validation like so:

#Recall we are using the lung dataset.
testInd <- create_folds(lung$time, lung$status, nfolds = 5, foldtype = "fullstrat")

Now testInd is a list of length 5 where each of the 5 items are the indices of the test set to be used for each fold. One could now use a for loop and iterate over 5-folds and index lung by the test indices (lung[testInd,]) for the test set and all but the test indices (lung[-testInd,]) for the training set and train these 5 models.

Evaluation

The MTLR package does not directly supply functions on how to evaluate the survival curves, but instead the methods to produce these curves. For a detailed discussion on a variety of evaluation metrics and access to some R-scripts containing evaluation code please see the paper Effective Ways to Build and Evaluate Individual Survival Distributions.