This document demonstrates how to use the efficient implementation of Survival Support Vector Machines as proposed in
Pölsterl, S., Navab, N., and Katouzian, A., Fast Training of Support Vector Machines for Survival Analysis, Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2015, Porto, Portugal, Lecture Notes in Computer Science, vol. 9285, pp. 243-259 (2015)
The source code and installation instructions are available at https://github.com/tum-camp/survival-support-vector-machine.
The main class of interest is survival.svm.FastSurvivalSVM
, which implements the different optimizers for training
a Survival Support Vector Machine. Training data consists of $n$ triplets $(\mathbf{x}_i, y_i, \delta_i)$, where
$\mathbf{x}_i$ is a $d$-dimensional feature vector, $y_i > 0$ the survival time or time of censoring, and $\delta_i \in \{0,1\}$ the binary event indicator. Using the training data, the objective is to minimize the following function:
The hyper-parameter $\alpha > 0$ determines the amount of regularization to apply: a smaller value increases the amount of regularization and a higher value reduces the amount of regularization. The hyper-parameter $r \in [0; 1]$ determines the trade-off between the ranking objective and the regresson objective. If $r = 1$ it reduces to the ranking objective, and if $r = 0$ to the regression objective. If the regression objective is used, it is advised to log-transform the survival/censoring time first.
In this example, I'm going to use the ranking objective ($r = 1$) and grid search to determine the best setting for the hyper-parameter $\alpha$.
The class survival.svm.FastSurvivalSVM
adheres to interfaces used in scikit-learn and thus it is possible to combine it with auxiliary classes and functions from scikit-learn. Here, I'm going to use GridSearchCV to determine which set hyper-parameters performs best for the Veteran's Lung Cancer data. Since, we require an event indicator $\delta_i$, which is boolean, and the survival/censoring time $y_i$ for training, we have to create a structured array that contains both information.
But first, we have to import the classes we are going to use.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pandas
import seaborn as sns
from sklearn.cross_validation import ShuffleSplit
from sklearn.grid_search import GridSearchCV
from survival.column import categorical_to_numeric, standardize
from survival.io import loadarff
from survival.metrics import concordance_index_censored
from survival.svm import FastSurvivalSVM
/home/sebp/miniconda3/envs/ssvm/lib/python3.5/site-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`. "`IPython.html.widgets` has moved to `ipywidgets`.", ShimWarning)
Next, load data of the Veteran's Administration Lung Cancer Trial from disk and convert it to numeric values. The data consists of 137 patients and 6 features. The primary outcome measure was death (Status
, Survival_in_days
).
The original data can be retrieved from http://lib.stat.cmu.edu/datasets/veteran.
Note that it does not matter how you name the fields corresponding to the event indicator and time, as long as the event indicator comes first.
data = loadarff("../data/veteran.arff")
data_x = data.drop(["Status", "Survival_in_days"], axis=1)
x = categorical_to_numeric(standardize(data_x))
y = np.empty(dtype=[("event", bool), ("time", float)], shape=data.shape[0])
y["event"] = (data["Status"] == "'dead'").values
y["time"] = data["Survival_in_days"].values
Now, we are essentially ready to start training, but before let's determine what the amount of censoring for this data is and plot the survival/censoring times.
n_censored = y.shape[0] - y["event"].sum()
print("%.1f%% of records are censored" % (n_censored / y.shape[0] * 100))
6.6% of records are censored
val, bins, patches = plt.hist((y["time"][y["event"]],
y["time"][-y["event"]]),
bins=30, stacked=True)
plt.legend(patches, ["Time of Death", "Time of Censoring"])
/home/sebp/miniconda3/envs/ssvm/lib/python3.5/site-packages/ipykernel/__main__.py:2: DeprecationWarning: numpy boolean negative, the `-` operator, is deprecated, use the `~` operator or the logical_not function instead. from ipykernel import kernelapp as app
<matplotlib.legend.Legend at 0x7f6229580390>
First, we need to create an initial model with default parameters that is subsequently used in the grid search. We are going to use a Red-Black tree to speed up optimization.
estimator = FastSurvivalSVM(optimizer="rbtree", max_iter=1000, tol=1e-6, random_state=0)
Next, we define a function for evaluating the performance of models during grid search. We use Harrell's concordance index.
def score_survival_model(model, X, y):
prediction = model.predict(X)
result = concordance_index_censored(y['event'], y['time'], prediction)
return result[0]
The last part of the setup specifies the set of parameters we want to try and how many repetitions of training and testing we want to perform for each parameter setting. In the end, the parameters that on average peformed best across all test sets (200 in this case) are selected. GridSearchCV can leverage multiple cores by evaluating multiple parameter settings concurrently (I use 4 jobs in this example).
param_grid = {'alpha': 2. ** np.arange(-12, 13, 2)}
cv = ShuffleSplit(x.shape[0], n_iter=200, train_size=0.5, random_state=0)
gcv = GridSearchCV(estimator, param_grid, score_survival_model, cv=cv,
n_jobs=4, iid=False, refit=False)
Finally, start the hyper-parameter search. This can take a while since a total of 13 * 200 = 2600
fits have to be evaluated.
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
gcv = gcv.fit(x, y)
Let's check what is the best average performance across 200 random train/test splits we got and the corresponding hyper-parameters.
gcv.best_score_, gcv.best_params_
(0.6932455298137199, {'alpha': 0.0009765625})
Finally, we retrieve all 200 test scores for each parameter setting and visualize their distribution by box plots.
cv_scores = {}
for i in range(len(gcv.grid_scores_)):
v = gcv.grid_scores_[i]
name = "%.5f" % v.parameters["alpha"]
cv_scores[name] = v.cv_validation_scores
sns.boxplot(pandas.DataFrame(cv_scores))
_, xtext = plt.xticks()
for t in xtext:
t.set_rotation("vertical")
This section demonstrates how to use the efficient implementation of Kernel Survival Support Vector Machines as proposed in
Pölsterl, S., Navab, N., and Katouzian, A., An Efficient Training Algorithm for Kernel Survival Support Vector Machines 4th Workshop on Machine Learning in Life Sciences, 23 September 2016, Riva del Garda, Italy
As kernel we are going to use the clinical kernel, because it distinguishes between continuous, ordinal, and nominal attributes.
from survival.svm import FastKernelSurvivalSVM
from survival.kernels import clinical_kernel
To use GridSearchCV
with a custom kernel, we need to pre-compute the squared kernel matrix and pass it to GridSearchCV.fit
later. It would also be possible to construct FastKernelSurvivalSVM
with kernel="rbf"
(or any other built-in kernel), which does not require pre-computing the kernel matrix.
kernel_matrix = clinical_kernel(data_x)
kssvm = FastKernelSurvivalSVM(optimizer="rbtree", kernel="precomputed", random_state=0)
kgcv = GridSearchCV(kssvm, param_grid, score_survival_model, cv=cv,
n_jobs=4, iid=False, refit=False)
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
kgcv = kgcv.fit(kernel_matrix, y)
Now, print the best average concordance index the corresponding parameters.
kgcv.best_score_, kgcv.best_params_
(0.70387140748796451, {'alpha': 0.015625})
Finally, we visualize the distribution of test scores obtained via cross-validation.
kcv_scores = {}
for i in range(len(kgcv.grid_scores_)):
v = kgcv.grid_scores_[i]
name = "%.5f" % v.parameters["alpha"]
kcv_scores[name] = v.cv_validation_scores
sns.boxplot(pandas.DataFrame(kcv_scores))
_, xtext = plt.xticks()
for t in xtext:
t.set_rotation("vertical")