%% BioMed_Central_Tex_Template_v1.06
%%                                      %
%  bmc_article.tex            ver: 1.06 %
%                                       %

%%IMPORTANT: do not delete the first line of this template
%%It must be present to enable the BMC Submission system to
%%recognise this template!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                     %%
%%  LaTeX template for BioMed Central  %%
%%     journal article submissions     %%
%%                                     %%
%%          <8 June 2012>              %%
%%                                     %%
%%                                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                 %%
%% For instructions on how to fill out this Tex template           %%
%% document please refer to Readme.html and the instructions for   %%
%% authors page on the biomed central website                      %%
%% http://www.biomedcentral.com/info/authors/                      %%
%%                                                                 %%
%% Please do not use \input{...} to include other tex files.       %%
%% Submit your LaTeX manuscript as one .tex document.              %%
%%                                                                 %%
%% All additional figures and files should be attached             %%
%% separately and not embedded in the \TeX\ document itself.       %%
%%                                                                 %%
%% BioMed Central currently use the MikTex distribution of         %%
%% TeX for Windows) of TeX and LaTeX.  This is available from      %%
%% http://www.miktex.org                                           %%
%%                                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% additional documentclass options:
%[doublespacing]
%  [linenumbers]   - put the line numbers on margins

%%% loading packages, author definitions

%\documentclass[twocolumn]{bmcart}% uncomment this for twocolumn layout and comment line below
\documentclass{bmcart}

%%% Load packages
\usepackage{amsthm,amsmath,bm} % packages for theorems and mathematics
\DeclareMathOperator*{\argmin}{\arg\!\min}
\RequirePackage{natbib} % package for the citation style
%\RequirePackage[authoryear]{natbib}% uncomment this for author-year bibliography
\RequirePackage[colorlinks=true, allcolors=blue]{hyperref} % package for the url references
\usepackage[utf8]{inputenc} %unicode support
%\usepackage[applemac]{inputenc} %applemac support if unicode package fails
%\usepackage[latin1]{inputenc} %UNIX support if unicode package fails
\usepackage{graphicx} % to use the graphicx package for the figures

% extra packages
\usepackage{booktabs}
%\usepackage{fontspec}
\usepackage[a4paper,top=3cm,bottom=3cm,left=3.5cm,right=3.5cm,marginparwidth=3cm]{geometry}
\usepackage[margin=0.50cm]{caption} % margin for the caption of the figures
\usepackage{subcaption}
\usepackage{enumerate}
\usepackage{etoolbox}
\appto\UrlBreaks{\do\-}
\usepackage{tabularx}
\usepackage{ragged2e}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                             %%
%%  If you wish to display your graphics for   %%
%%  your own use using includegraphic or       %%
%%  includegraphics, then comment out the      %%
%%  following two lines of code.               %%
%%  NB: These line *must* be included when     %%
%%  submitting to BMC.                         %%
%%  All figure files must be submitted as      %%
%%  separate graphics through the BMC          %%
%%  submission process, not included in the    %%
%%  submitted article.                         %%
%%                                             %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\def\includegraphic{}
\def\includegraphics{}



%%% Put your definitions there:
\startlocaldefs
\endlocaldefs


%%% Begin ...
\begin{document}%[doublespacing]

%%% Start of article front matter
\begin{frontmatter}

\begin{fmbox}
\dochead{Research}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the title of your article here     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Survival prediction models since liver transplantation - comparisons between Cox models and machine learning techniques}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors here                   %%
%%                                          %%
%% Specify information, if available,       %%
%% in the form:                             %%
%%   <key>={<id1>,<id2>}                    %%
%%   <key>=                                 %%
%% Comment or delete the keys which are     %%
%% not used. Repeat \author command as much %%
%% as required.                             %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\author[
   addressref={aff1,aff2,aff3},                   % id's of addresses, e.g. {aff1,aff2}
   corref={aff1,aff2,aff3},                       % id of corresponding address, if any
   noteref={n1},
   email={G.Kantidakis@lumc.nl, georgios.kantidakis@eortc.org}   % email address
]{\inits{G}\fnm{Georgios} \snm{Kantidakis}}
\author[
   addressref={aff2},
]{\inits{H}\fnm{Hein} \snm{Putter}}
\author[
   addressref={aff1},
]{\inits{C}\fnm{Carlo} \snm{Lancia}}
\author[
   addressref={aff4},
]{\inits{Jde}\fnm{Jacob de} \snm{Boer}}
\author[
   addressref={aff4},
]{\inits{AE}\fnm{Andries E} \snm{Braat}}
\author[
   addressref={aff1,aff2,aff5},
]{\inits{M}\fnm{Marta} \snm{Fiocco}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors' addresses here        %%
%%                                          %%
%% Repeat \address commands as much as      %%
%% required.                                %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\address[id=aff1]{%                           % unique id
  \orgname{Mathematical Institute (MI) Leiden University}, % university, etc
  \street{Niels Bohrweg 1},                     %
  \postcode{2333 CA}                                % post or zip code
  \city{Leiden},                              % city
  \cny{the Netherlands}                                    % country
}
\address[id=aff2]{%                           % unique id
  \orgname{Department of Biomedical Data Sciences, Section Medical Statistics, Leiden University Medical Center (LUMC)}, % university, etc
  \street{Albinusdreef 2},                     %
  \postcode{2333 ZA}                                % post or zip code
  \city{Leiden},                              % city
  \cny{The Netherlands}                                    % country
}
\address[id=aff3]{%                           % unique id
  \orgname{Department of Statistics, European Organisation for Research and Treatment of Cancer (EORTC) Headquarters}, % university, etc
  \street{Ave E. Mounier 83/11},                     %
  \postcode{1200}                                % post or zip code
  \city{Brussels},                              % city
  \cny{Belgium}                                    % country
}

\address[id=aff4]{%                           % unique id
  \orgname{Department of Surgery, Leiden University Medical Center (LUMC)}, % university, etc
  \street{Albinusdreef 2},                     %
  \postcode{2333 ZA}                                % post or zip code
  \city{Leiden},                              % city
  \cny{the Netherlands}                                    % country
}
\address[id=aff5]{%                           % unique id
  \orgname{Trial and Data Center, Princess Máxima Center for pediatric oncology (PMC)}, % university, etc
  \street{Heidelberglaan 25},                     %
  \postcode{3584 CS}                                % post or zip code
  \city{Utrecht},                              % city
  \cny{the Netherlands}                                    % country
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter short notes here                   %%
%%                                          %%
%% Short notes will be after addresses      %%
%% on first page.                           %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{artnotes}
%\note{Sample of title note}     % note to the article
\note[id=n1]{The majority of this work was done at Leiden University} % note, connected to author
\end{artnotes}

\end{fmbox}% comment this for two column layout

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Abstract begins here                 %%
%%                                          %%
%% Please refer to the Instructions for     %%
%% authors on http://www.biomedcentral.com  %%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstractbox}

\begin{abstract} % abstract
\parttitle{Background} Predicting survival of recipients after liver transplantation is regarded as one of the most important challenges in contemporary medicine. Hence, improving on current prediction models is of great interest.

Nowadays, there is a strong discussion in the medical field about machine learning (ML) and whether it has greater potential than traditional regression models when dealing with complex data. Criticism to ML is related to unsuitable performance measures and lack of interpretability which is important for clinicians. 

\parttitle{Methods}
In this paper, ML techniques such as random forests and neural networks are applied to large data of 62294 patients from the United States with 97 predictors selected on clinical/statistical grounds, over more than 600, to predict survival from transplantation. Of particular interest is also the identification of potential risk factors. A comparison is performed between 3 different Cox models (with all variables, backward selection and LASSO) and 3 machine learning techniques: a random survival forest and 2 partial logistic artificial neural networks (PLANNs). For PLANNs, novel extensions to their original specification are tested. Emphasis is given on the advantages and pitfalls of each method and on the interpretability of the ML techniques. 

\parttitle{Results}
Well-established predictive measures are employed from the survival field (C-index, Brier score and Integrated Brier Score) and the strongest prognostic factors are identified for each model. Clinical endpoint is overall graft-survival defined as the time between transplantation and the date of graft-failure or death. The random survival forest shows slightly better predictive performance than Cox models based on the C-index. Neural networks show better performance than both Cox models and random survival forest based on the Integrated Brier Score at 10 years.

\parttitle{Conclusion}
In this work, it is shown that machine learning techniques can be a useful tool for both prediction and interpretation in the survival context. From the ML techniques examined here, PLANN with 1 hidden layer predicts survival probabilities the most accurately, being as calibrated as the Cox model with all variables.      

\parttitle{Trial registration}
Retrospective data were provided by the Scientific Registry of Transplant Recipients under Data Use Agreement number 9477 for analysis of risk factors after liver transplantation. 

\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The keywords begin here                  %%
%%                                          %%
%% Put each keyword in separate \kwd{}.     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{keyword}
\kwd{Random Survival Forests}
\kwd{Neural Networks}
\kwd{Predictive Performance}
\kwd{Risk Factors}
\kwd{Post-transplantation}
\kwd{Survival analysis}
\end{keyword}

% MSC classifications codes, if any
%\begin{keyword}[class=AMS]
%\kwd[Primary ]{}
%\kwd{}
%\kwd[; secondary ]{}
%\end{keyword}

\end{abstractbox}
%
%\end{fmbox}% uncomment this for twcolumn layout

\end{frontmatter}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Main Body begins here                %%
%%                                          %%
%% Please refer to the instructions for     %%
%% authors on:                              %%
%% http://www.biomedcentral.com/info/authors%%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%% See the Results and Discussion section   %%
%% for details on how to create sub-sections%%
%%                                          %%
%% use \cite{...} to cite references        %%
%%  \cite{koon} and                         %%
%%  \cite{oreg,khar,zvai,xjon,schn,pond}    %%
%%  \nocite{smith,marg,hunn,advi,koha,mouse}%%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%% start of article main body
% <put your article body there>

%%%%%%%%%%%%%%%%
%% Background %%
%%
\section*{Background}
Liver transplantation (LT) is the second most common type of transplant surgery in the United States after kidney \cite{Grinyo2013WhyImportant}. Over the last decades, the success of liver transplants has improved survival outcome for a large number of patients suffering from chronic liver disease everywhere on earth \cite{Merion2005TheTransplantation}. Availability of donor organs is a major limitation especially when compared with the growing demand of liver candidates due to the enlargement of age limits. Therefore, improvement on current prediction models for survival since LT is important. 

In the statistical field, there is an open discussion about the value of machine learning (ML) versus statistical models (SM) for medical application. For survival data, the most commonly applied statistical model is the Cox proportional hazards regression model  \cite{Cox1972RegressionLife-Tables}. This model allows a straightforward interpretation, but is at the same time restricted to the proportional hazards assumption. On the other hand, ML techniques are assumption-free and data adaptive which means that they can be effectively employed for modelling complex data. However, there is a danger of over-fitting. 

Biganzoli \textit{et al.}  proposed a partial logistic regression approach of feed forward neural networks (PLANN) for flexible modelling of survival data \cite{Biganzoli1998FeedApproach.}. By using the time interval as an input in a longitudinally transformed feed forward network with logistic activation and entropy error function, they estimated smoothed discrete hazards at each time interval in the output layer. This is a well known approach for modelling survival neural networks \cite{Wang2019MachineSurvey}.  In 2000, Xiang \textit{et al.} \cite{Xiang2000ComparisonData} compared the performance of 3 existing neural network methods for right censored data (the Faraggi-Simon \cite{Faraggi1995AData}, the Liestol-Andersen-Andersen \cite{LiestlK1994SurvivalNets} and a modification of the Buckley-James method \cite{Buckley1979LinearData}) with Cox models in a Monte Carlo simulation study. None of the networks outperformed the Cox models and they only performed as good as Cox for some scenarios. Lisboa \textit{et al.} extended the PLANN approach introducing a Bayesian framework which can perform Automatic Relevance Determination for survival data (PLANN-ARD) \cite{Lisboa2003ACancer}. Several applications of the PLANN and the PLANN-ARD methods can be found in the literature \cite{Biganzoli2002AData,Biglarian2013NonlinearNetwork,Jones2006AnCarcinoma, Taktak2007Double-blindStudy}. They show potential for neural networks in systems with non-linearity and complex interactions between factors. Here extensions of the PLANN approach for big LT data are examined.   

Predicting survival after LT is hard as it depends on many factors and is associated with donor, transplant and recipient characteristics whose importance changes over time and per outcome measure \cite{Blok2018IdentificationTime}. Models that combine donor and recipient characteristics have usually better performance for predicting overall graft-survival and particularly those that include sufficient donor risk factors have better performance for long-term graft survival \cite{deBoer2019PredictiveTransplantation}. Our aim is to identify potential risk factors and to compare their relative importance using 2 ML methods (random survival forests, survival neural networks) and Cox models. Emphasis is given on interpretability and predictive performance of the models as well as their potential for medical application.  

\section*{Methods}
An analysis is presented on survival data after LT based on 62294 patients from the United States. Information was collected from the United Network of Organ Sharing (UNOS)\footnote{UNOS is a non-profit and scientific organisation in the United States which arranges organ donation and transplantation. For more information visit its website \url{https://unos.org}. }. After extensive pre-processing from a set of more than 600 covariates and imputation of missing values, 97 variables were included in the final dataset based on clinical and statistical considerations (see Additional file 1); 52 donor and 45 liver recipient characteristics. As the UNOS data is large in both number of observations and covariates, it is of interest to see how ML algorithms - which are able to capture naturally multi-way interactions between variables and can deal with big datasets - will perform compared to Cox models. The clinical endpoint is overall graft-survival defined as the time between LT and graft-failure or death. The choice for this endpoint was made for two reasons 1) it is of primary interest for clinicians and 2) it is likely the most appropriate outcome measure since LT \cite{deBoer2019PredictiveTransplantation}.

\subsection*{Data collection and imputation technique}
UNOS manages the Organ Procurement and Transplantation Network (OPTN) and together they collect, organise and maintain data of statistical information regarding organ transplants in the Scientific Registry of Transplant Recipients (SRTR) database\footnote{Dictionary for variables details is provided at: \url{https://www.srtr.org/requesting-srtr-data/saf-data-dictionary/}.}. SRTR gathers data from local Organ Procurement Organisations (OPO) and from OPTN (primary source). It includes data from transplantations performed in the United States from 1988 onwards. This information is used to set priorities and seek improvements in the organ donation process. 

The data provided by UNOS included 62294 patients who underwent LT surgery from 2005 to 2015 (project under DUA number 9477). SAF contained 657 variables for both donors and patients (candidates and recipients). Among these, 97 candidate risk factors - 52 donor and 45 patient characteristics - were selected. This resulted in a final dataset with 76 categorical and 21 continuous variables amounting to 2.2\% missing data overall. To reconstruct the missing values the \texttt{missForest} algorithm \cite{Stekhoven2012Missforest-Non-parametricData} was applied. This is a non-parametric imputation method that does not make explicit assumptions about the functional form of the data and builds a random forest model for each variable (500 trees were used). It specifies the model to predict missing values by using information based on the observed values. This is the most exhaustive and accurate of all random forests algorithms used for missing data imputation, because all possible variable combinations are checked as responses.  

\subsection*{Cox proportional hazard regression models}

In survival analysis, the focus is on the time till the occurrence of the event of interest (here graft-failure or death). The Cox proportional hazards model is usually employed to estimate the effect of risk factors on the outcome of interest \cite{Cox1972RegressionLife-Tables}.

Data with sample size $n$ consist of the independent observations from the triple $(T, D, X)$ i.e. $(t_{1}, d_{1}, x_{1}), \cdots, (t_{n}, d_{n}, x_{n})$. For the $i^{th}$ individual, $t_{i}$ is the survival time, $d_{i}$ the indicator ($d_{i} = 1$  if the event occurred and $d_{i} = 0$ if the observation is right censored) and $x_{i}$ is the vector of predictors $(x_{1},\cdots,x_{p})$. The hazard function of the Cox model with time-fixed covariates is as follows:

 \begin{eqnarray} h(t|X)  = h_{0}(t)\exp(X^{T}\bm{\beta}), \end{eqnarray} where $h(t|X)$ is the hazard at time t given predictor values X, $h_{0}(t)$ is an arbitrary baseline hazard and $\bm{\beta} = (\beta_{1},\cdots,\beta_{p})$ is a parameter vector.

The corresponding partial likelihood can be written as:
 
 \begin{eqnarray} \label{equat:cox_log_partial_lik}
     L(\bm{\beta}) = \prod_{i=1}^{D} \frac{\exp{(\sum_{k = 1}^{p} \beta_{k} X_{ik})}}{\sum_{j\in R(t_{i})} \exp{(\sum_{k = 1}^{p} \beta_{k} Z_{jk})}},
 \end{eqnarray}
 where $D$ is the set of failures, and $R(t_{i})$ is the risk set at time $t_{i}$ of all individuals who are still in the study at the time just before time $t_{i}$.  This function is then maximised over $\bm{\beta}$ to estimate the model parameters. Two other Cox models were 1) a Cox model with a backward elimination or 2) a penalised Cox regression with the Least Angle and Selection Operator (LASSO). 

For the first, a numerically stable version of the backward elimination on factors was used \cite{Lawless1978EfficientModels}. This method estimates the full model and computes approximate Wald statistics by computing conditional maximum likelihood estimates - assuming multivariate normality of estimates. Factors that require multiple degrees of freedom are dropped or retained as a group. 

The latter approach uses a combination of selection and regularisation \cite{Tibshirani1997THEMODEL}. Denote the log-partial likelihood by $\ell(\bm{\beta}) = logL(\bm{\beta})$. The vector $\bm{\beta}$ is estimated via the criterion:

\begin{eqnarray} \label{equat:eq_Lasso1}
  \bm{\hat{\beta}} =\argmin[ \ell(\bm{\beta})], \quad
  \textrm{subject} \ \textrm{to} \ \sum_{j = 1}^{p} |\beta_{j}| \leq s
\end{eqnarray}
with $s$ a user specified positive parameter.

Equation (\ref{equat:eq_Lasso1}) can also be rewritten as

\begin{eqnarray} \label{equat:eq_Lasso2}
    \bm{\hat{\beta}} = \argmin_{\beta} \Big( \ell (\beta ) + \lambda_{LASSO} \sum_{j = 1}^{p} |\beta_{j}| \Big).
\end{eqnarray}
The quantity \(\sum_{j = 1}^{p} |\beta_{j}| \) is also known as the $L_{1}$-norm and performs regularisation to the log-partial likelihood. The term $\lambda_{LASSO}$ is a non-negative constant that assigns the amount of penalisation. Larger values for the parameter mean larger penalty to the $\beta_{j}$ coefficients and enlarged shrinkage towards zero. 

The tuning parameter $s$ in equation (\ref{equat:eq_Lasso1}) or equivalently parameter $\lambda_{Lasso}$ in equation (\ref{equat:eq_Lasso2}) is the controlling mechanism for the variance of the model. Higher values reduce further the variance but introduce at the same time more bias (variance-bias trade off). To find a suitable value for this parameter 5-fold cross-validation was performed to minimise the prediction error; here in terms of the cross-validated log-partial likelihood (CVPL) \cite{Verweij1993Cross-validationAnalysis}

\begin{eqnarray} \label{chap3:eq_lasso_cvpl}
    \textit{CVPL}(s) = \sum_{i = 1}^{n} \big( \ell (\hat{\beta}_{(-i)}(s)) - \ell_{{(-i)}}(\hat{\beta}_{(-i)}(s)) \big),
\end{eqnarray}
where $\ell_{(-i)}(\beta)$ is the partial log-likelihood of equation (\ref{equat:cox_log_partial_lik}) when individual $i$ is excluded. Therefore, the term $\ell (\hat{\beta}_{(-i)}) - \ell_{{(-i)}}(\hat{\beta}_{(-i)})$ represents the contribution of observation $i$. The value that maximizes $\ell_{{(-i)}}(\beta_{(-i)})$ is denoted by $\hat{\beta}_{(-i)}$.


\subsection*{Random forests for survival analysis}
Random Survival Forests (RSFs) are an ensemble tree method for survival analysis of right censored data \cite{Ishwaran2008RandomForests} adapted from random forests \cite{Breiman2001RandomForests}. The main idea of random forests is to get a series of decision trees - which can capture complex interactions but are notorious for their high variance - and obtain a collection averaging their characteristics. In this way weak learners (the individual trees) are turned into strong learners (the ensemble) \cite{Hastie2009TheLearning}. 

For RSFs, randomness is introduced in two ways: bootstrapping a number of patients at each tree $\mathcal{B}$ times and selecting a subset of variables for growing each node. During growing each survival tree, a recursive application of binary splitting is performed per region (called node) on a specific predictor in such a way that survival difference between daughter nodes is maximised and difference within them is minimised. Splitting is terminated when a certain criterion is reached (these nodes are called terminal). The most commonly used splitting criteria are the log-rank test by Segal \cite{Segal1988RegressionData} and the log-rank score test by Hothorn and Lausen \cite{Hothorn2003OnStatistics}. Each terminal node should have at least a pre-specified number of unique events. Combining information from the $\mathcal{B}$ trees, an ensemble cumulative hazard estimate can be calculated using the Nelson-Aalen methodology and subsequently the survival probabilities.

The fundamental principle behind each survival tree is the conservation of events. It is used to define ensemble mortality, a new type of predicted outcome for survival data. This principle asserts that the sum of estimated cumulative hazard estimate over time is equal to the total number of deaths, therefore the total number of deaths is conserved within each terminal node $\mathcal{H}$ \cite{Ishwaran2008RandomForests}. It can also be shown that the total number of deaths is also conserved in a tree grown from the original non-bootstrapped data. RSFs can handle both data with large sample size and vast number of predictors. Moreover, they can reach remarkable stability combining the results of many trees. However, combining an ensemble of trees downgrades significantly the intuitive interpretation of a single tree. 

\subsection*{Survival neural networks}
Artificial neural networks (ANNs) are a machine learning method able to model non-linear relationships between prognostic factors with great flexibility. These systems are inspired from biological neural networks that aimed at imitating the human brain activity \cite{vanGerven2017Editorial:Processing}. A ANN has a layered structure and is based on a collection of connected units called nodes or neurons which comprise a layer. The input layer picks up the signals and passes them through transformation functions to the next layer which is called "hidden". A network may have more than one hidden layer that connects with the previous and transmit signals towards the output layer. Connections between artificial neurons are called edges. Artificial neurons and edges have a weight (connection strength) which adjusts as learning proceeds. It increases or decreases the strength of the signal of each connection according to its sign. For the purpose of training, a target is defined, which is the observed outcome. The simplest form of a NN is the single layer feed-forward perceptron with the input layer, one hidden layer and the output layer \cite{Minsky1969PerceptronsGeometry}.

The application of NNs has been extended to survival analysis over the years \cite{LiestlK1994SurvivalNets}. Different approaches have been considered; some model the survival probability $\mathcal{S}(t)$ directly or the unconditional probability of death $\mathcal{F}(t)$ whereas other approaches estimate the conditional hazard $h(t)$ \cite{Wang2019MachineSurvey}. They can be distinguished according to the method used to deal with the censoring mechanism. Some networks have $k$ output nodes \cite{Lapuerta1995UseDisease.} - where $k$ denotes $k$ separate time intervals - while others have a single output node. 

In this research, the method of Biganzoli was applied, which specifies a partial logistic feed-forward artificial neural network (PLANN) with a single output node \cite{Biganzoli1998FeedApproach.}. This method uses as inputs the prognostic factors and the survival times to increase the predictive ability of the model. Data have to be transformed into a longitudinal format with the survival times being divided into a set of $k$ non-overlapping intervals (months or years) $I_{k} = (\tau_{k-1}, \tau_{k}]$, with $0=\tau_{o}<\tau_{1}<\cdots<\tau_{K}$ a set of pre-defined time points. On the training data, each individual is repeated for the number of intervals he/she was observed in the study and on the test data for all time intervals. PLANN provides the discrete conditional probability of dying $\mathcal{P}(T \in I_{k} \mid T>\tau_{k-1})$  using as transformation function of both input and output layers the logistic (sigmoid) function: 

\begin{eqnarray} \label{chap3:eq_logistic}
    f(\eta) &=& \frac{1}{1 + e^{-\eta}},
\end{eqnarray}
where $\eta = \sum_{i = 1}^{p}w_{i} X_{i}$ is the summed linear combination of the weights $w_{i}$ of input-hidden layer and the input variables $X_{i}$ ($i = 1, 2, \cdots, p$). 

The contribution to the log-likelihood is expanded as $\sum_{\eta \in I_{0}}^{I_{k}} \log(f(\eta)) $ over the intervals at which the specific patient is at risk. The output node is one large target vector with 0 if the event did not occur and 1 if the event occurred in a specific time interval. Therefore, such a network first estimates the hazard for each interval $ h_{k} = P(\tau_{k-1} < T \leq \tau_{k} | T > \tau_{k-1})$ and then $S(t) = \prod_{k: t_{k} \leq t} (1 - h_{k})$. 

In this work, novel extensions in the specification of the PLANN are tested. Two new transformation functions were investigated for the input-hidden layer the rectified linear unit (ReLU)

\begin{eqnarray} \label{chap3:eq_rectified_linear_unit}
    f(\eta) = \eta^{+} = \max{(0, \eta)},
\end{eqnarray}
which is the most used activation function for NNs and the hyperbolic tangent (tanh)

\begin{eqnarray} \label{chap3:eq_hyperbolic_tangent}
    f(\eta) &=& \frac{1 - e^{-2\eta}}{1 + e^{-2\eta}}.
\end{eqnarray}
These functions can be seen as different modulators of the degree of non-linearity implied by the input and the hidden layer.

The PLANN was expanded in 2 hidden layers with same node size and identical activation functions for input-hidden 1 and hidden 1 - hidden 2 layers. The $k$ non-overlapping intervals of the survival times were treated as $k$ separate variables. In this way, the contribution of each interval to the predictions of the model using the relative importance method by Garson \cite{Garson1991InterpretingWeights} and its extension for 2 hidden layers can be obtained (see subsection "Interpretability of the models" below and Additional file 1).   

\subsection*{Model training}
The split sample approach was employed; data was split randomly into two complementary parts, a training set (2/3) and a test set (1/3) under the same event/censoring proportions. To tune a model, 5-fold cross validation was performed in the training set for the machine learning techniques (and for Cox LASSO). Training data was divided into 5 folds. Each time 4 folds were used to train a model and the remaining fold was used to validate its performance and the procedure was repeated for all combination of folds. Tuning of the hyper-parameters was done using grid search and performance of final models was assessed on the test set. 

For RSF the \texttt{randomForestSRC} package \cite{KogalurU2019RandomForestSRC:RF-SRC} of the R programming language was used. Parameters tuned were \texttt{ntree} the number of bootstrapped trees grown (range 1-500), \texttt{mtry} the number of candidate variables examined at each split point (range 5-47), \texttt{nsplit} the number of split points at which an $X$-variable is tested using the log-rank splitting rule (range 3-7) \cite{Segal1988RegressionData} and \texttt{nodesize} (range 10-100) the average number of observations in the terminal nodes across the forest. In general, parameters \texttt{ntree} and \texttt{mtry} are the most fundamental for RSF. The parameter \texttt{ntree} modulates the consistency of the forest's performance and \texttt{mtry} controls an important part of randomness during the growth of decision trees. Parameter \texttt{nsplit} with \texttt{nsplit} $>$ 0 can be used to trigger a randomised selection of exactly nsplit points for each of the mtry variables within a node $h$. Last, parameter \texttt{nodesize} plays an important role in the topology of the trees as it controls the average node size of the forest. Large values in \texttt{nodesize} parameter will essentially force the forest to under-grow whereas small values will lead each tree to keep growing on with more and more noisy variables being selected. The best combination of the parameters was determined based on the error of the forests defined as $E = 1 - C$, where C is Harrell's concordance index \cite{Houwelingen2012DynamicAnalysis} based on the ensemble survival mortality of individuals. Figure \ref{app:fig:rsf_3dtuning} shows the cross-validated error $E$ for different values of tuning parameters \texttt{mtry}, \texttt{nodesize} and \texttt{nsplit}. Optimal values were given by \texttt{mtry} = 12, \texttt{nodesize} = 50 and \texttt{nsplit} = 5.    

\begin{figure}[h!]
%\centering
\includegraphics[width=12cm, height =8cm]{figure1.png}
\caption{\label{app:fig:rsf_3dtuning}\csentence{Cross-validation applied on a 3-D space for the RSF.} Cross-validation applied on a 3-D space for the RSF. Parameters tuned are the number of candidate variables (\texttt{mtry}), the average node size across the forest (\texttt{nodesize}) and the number of split points for each x-variable (\texttt{nsplit}). The bar on the right side of the graph shows the \textbf{error E} corresponding to combinations on training data. Dark blue dots correspond to lower prediction error E and represent the best hyper-parameter combinations whereas green and yellow dots correspond to poorly performed combinations. }
\end{figure}

For survival NNs, standard software of implementation was not available in R. Therefore, data transformation to longitudinal format was required. The task was turned into a classification problem where patients were divided into maximum 10 time intervals on the training set and exactly 10 intervals on the test set (new unseen data) and the intervals were added to the other features as covariates to estimate conditional hazard probabilities. This led to a training set of 194635, a validation set of 415300, and a test set of 207640 observations in long format. Variables were presented in dummy coding and all continuous factors were standardized. Model tuning was performed in R with the package \texttt{keras} \cite{JJAllaireKeras} which is an interface for the original state-of-the-art NN library written in Python programming language. \texttt{keras} runs on top of \texttt{tensorflow} \cite{Abadi2016TensorFlow:Learning} which is a symbolic maths library used for machine learning. Two of the main advantages of the package are that it allows the use of distributed training of deep learning models on clusters of graphic processing units and the specification of many building blocks such as layers, objectives, activation functions, optimisers. 

To narrow down the grid of point combinations, search for training data was performed on a 5-D space of some of the most fundamental hyper-parameters. Those are \texttt{nodesize} the size of nodes in the hidden layer(s) (range 10-130), \texttt{dropout rate} that randomly selects the amount of nodes to be dropped-out with a given probability (0.1, 0.2 or 0.3), \texttt{learning rate}  which is the step size of weight iteration (0.01, 0.1 or 0.2), \texttt{momentum} which helps to accelerate gradient vectors (0.8 or 0.9) and \texttt{weak class weight} that defines the weight of minority class (1, 2.236 or 10). \texttt{nodesize} defines the number of weights of the network and consequently the amount of its complexity. Having 129 inputs in total (119 potentially prognostic variable levels in dummy coding + 10 time intervals), the optimum node size parameter will be somewhere in the range 10 - 130. Regarding the rest of the parameters, \texttt{dropout rate} is a technique that can reduce over-fitting \cite{Srivastava2014Dropout:Overfitting}, \texttt{learning rate} adjusts how fast the stochastic gradient descent iterative method uses stochastic approximation. \texttt{Momentum} can accelerate the stochastic gradient vectors in the right directions and \texttt{weak class weight} can be used for re-weighting unbalanced classes. In this dataset, 30.9 \% of the patients experienced the event of interest so it was investigated if the performance can be improved by re-weighting the minority class with 2.236 (= 69.1 / 30.9). Moreover, in the long data format each patient was replicated for a maximum of 10 intervals on the training data and for exactly 10 intervals on the validation data, it was investigated if the performance can be improved by re-weighting 10 times the minority class.

For survival NNs there is no well-established measure of cross-validation performance as they have rarely been used in practice. Thus, the best combinations of hyper-parameters were selected based on the Integrated Brier Score (IBS) -a measure of that summarises the time-dependent Brier score estimated at different time points until 10 years in one value. The choice of the hyper-parameters and their range (grid search) for both Random Survival Forests and Neural Networks was decided based on expert user recommendations from the world-wide web \cite{RandomClassification, HotExchange}. Afterwards, the search was repeated on a narrower domain centred around highly performing combinations. 

\subsection*{Assessing predictive performance on test data}
To assess the final predictive performance of the models the concordance index, the Brier score, and the Integrated Brier Score (IBS) were applied.

The most popular measure of model performance in a survival context is the concordance index \cite{Harrell1996MultivariableErrors} which computes the proportion of pairs of observations for which the survival times and model predictions order are concordant taking into account censoring. It takes values typically in the range 0.5 - 1 with higher values denoting higher ability of the model to discriminate and 0.5 indicating no discrimination. The C-index cannot be defined for neural network models since it relies on the ordering of individuals according to prognosis and there is no unique ordering between the subjects. At one year individual i may have better survival probability than individual j, but this could be reversed for a different time point.

The C-index provides a rank statistic between the observations that is not time-dependent. Following van Houwelingen and le Cessie \cite{VanHouwelingen1990PredictiveModels} a time-dependent prediction error is defined as 

\begin{eqnarray} \label{chap4:eq_brier_score}
\textit{Brier}(y, \hat{S}(t_{0}|x)) &=& (y - \hat{S}(t_{0}|x))^2,
\end{eqnarray} 
where $\hat{S}(t_{0}|x)$ is the model-based probabilistic prediction for the survival of an individual beyond $t_{0}$ given the predictor $x$, and $y = 1\{ t > t_{0} \}$ is the actual observation ignoring censoring. The expected value with respect to a new observation $Y_{new}$ under the true model $S(t_{0}|x)$ can be written as: 

\begin{equation}
    E[\textit{Brier}(Y_{new}, \hat{S}(t_{0}|x))] = S(t_{0}|x) (1 - S(t_{0}|x)) + (S(t_{0}|x) - \hat{S}(t_{0}|x))^{2}.
\end{equation}
The Brier Score consists of two components: the "true variation" $ S(t_{0}|x) (1 - S(t_{0}|x))$ and the error due to the model $(S(t_{0}|x) - \hat{S}(t_{0}|x))^{2}$. A perfect prediction is only possible if $S(t_{0}|x) = 0$ or $S(t_{0}|x) = 1$. In practice the two components cannot be separated since the true $S(t_{0}|x)$ is unknown. 

To assess the performance of a prediction rule in actual data, censored observations before time $t_{0}$ must be considered. To calculate Brier Score when censored observations are present, Graf proposed the use of inverse probability of censoring weighting \cite{Graf1999AssessmentData.}. Then an estimate of the average prediction error of the prediction model $\hat{S}(t|x)$ at time $t = t_{0}$ is 

\begin{equation} \label{eq:brier_score}
    Err_{Score}(\hat{S}, t_{0}) = \frac{1}{n} \sum_{i} 1\{d_{i} = 1 \vee t_{i} > t_{0}\} \frac{Score(1\{t_{i} > t_{0}\}, \hat{S}(t_{0}|x_{i}))}{\hat{C}(\min(t_{i}-, t_{0}) | x_{i})} 
\end{equation}

In (\ref{eq:brier_score}) $Score$ is the Brier Score for the prediction model. It ranges typically from 0 to 0.25 with a lower value meaning smaller prediction error.

As Brier score is calculated at different time-points, an overall measure of prediction error the Integrated Brier Score (IBS) can be used to summarise the prediction error curves over the whole range up the horizon $\int_{0}^{t_{hor}}Err_{Score}(\hat{S}, t_{0})dt_{0}$ (here $t_{hor}$ = 10 years) \cite{Houwelingen2012DynamicAnalysis}. IBS takes values in the same range as the Brier score.


\subsection*{Interpretability of the models}
Interpretation of models is of great importance for the medical community. It is well known that Cox models offer a straightforward interpretation through hazard ratios.

For neural networks with one hidden layer the connection weights algorithm by Garson  \cite{Garson1991InterpretingWeights} -- later modified by Goh \cite{Goh1995Back-propagationSystems} -- can provide information about the mechanism of the weights. The idea behind this algorithm is that inputs with larger connection weights produce greater intensities of signal transfer. As a result, these inputs will be more important for the model. Garson's algorithm can be used to determine relative importance of each input variable, partitioning the weights in the network. Their absolute values are used to specify percentage of importance. Note that the algorithm does not provide the direction of relationships, so it remains uncertain whether the relative importance indicates a positive or a negative effect. For details about the algorithm see \cite{Olden2002IlluminatingNetworks}. During this work, the algorithm was extended for 2 hidden layers to obtain the relative importance of each variable (for the implementation see algorithm 1 on Additional file 1).   

Random survival forests relies on two methods which can provide interpretability: variable importance (VIMP) and minimal depth \cite{Ishwaran2010High-dimensionalData}. The former is associated with the prediction error before and after the permutation of a prognostic factor. Large importance values indicate variables with strong predictive ability. The latter is related to the forest topology as it assesses the predictive value of a variable by computing its depth compared to the root node of a tree. VIMP is more frequently reported than minimal depth in the literature  \cite{Ishwaran2019StandardSurvival}. For both methods interpretation is available only for variable entities and not for each variable level.


\section*{Results}
Administrative censoring was applied to the UNOS data at 10 years. Median follow-up is equal to 5.36 years (95\% CI: 5.19 - 5.59 years) and it was estimated with reverse Kaplan-Meier \cite{Schemper1996ATime}. Clinical endpoint is overall graft-survival (OGS). From the total number of patients,  69.1\% was alive/censored and 30.9\% experienced the event of interest (graft-failure or death). 3 models were used from the Cox family to predict survival outcome: a) a model with all 97 prognostic factors, b) a model with backward selection and c) a model based on the LASSO method for variable selection. Furthermore, 3 machine learning methods were employed: a) a random survival forest, b) a NN with one hidden layer and c) a NN with two hidden layers. 

\subsection*{Comparisons between models}
In this section a direct comparison of the 6 models is illustrated in terms of variable importance on the training set and predictive performance on the test set. Specification of the variables with dummy coding included 119 variable levels from the 97 potentially prognostic factors. For NNs - to apply and extend the methodology of Biganzoli - follow-up time was divided into 10 time intervals $(0,1], (1, 2], \cdots , $ $(9, 10] $ denoting years since transplantation. For Cox models and RSF exact time points were used.  All analyses were performed in R programming language version 3.5.3 \cite{RCoreTeamComputing}.

The backward and the LASSO methods selected 28 and 45 predictors respectively. For the Cox model with all predictors the proportional hazards assumption was violated for 17 out of 97 variables.
5-fold cross-validation in the training data resulted in the following optimal hyper-parameters combinations for the machine learning techniques: 

\begin{itemize}
    \item For the Random Survival Forest \texttt{nodesize} = 50, \texttt{mtry} = 12, \texttt{nsplit} = 5 and \texttt{ntree} = 300. Stratified bootstrap sub-sampling of half the patients was used per tree (due to the large training time required).
    \item For the neural network with 1 hidden layer \texttt{activation function} = "sigmoid" (for the input-hidden layer), \texttt{node size} = 85, \texttt{dropout rate} = 0.2, \texttt{learning rate} = 0.2, \texttt{momentum} = 0.9 and \texttt{weak class weight} = 1. 
    \item  For the neural network with 2 hidden layers \texttt{activation function} = "sigmoid" (for the input-hidden 1 and the hidden 1-hidden 2 layers), \texttt{node size} = 110, \texttt{dropout rate} = 0.1, \texttt{learning rate} = 0.2, \texttt{momentum} = 0.9 and \texttt{weak class weight} = 1.
\end{itemize}

\subsubsection*{Global performance measures}
The global performance measures on test data are provided in Table \ref{tab:global_measures}. Examining the Integrated Brier Score (IBS), the NNs with 1 and with 2 hidden layers have the lowest (IBS = 0.180) followed by the RSF (IBS = 0.182). Cox models have a comparable performance (IBS = 0.183). The best model in terms of C-index is the Random Survival Forest (0.622) while the Cox models with all variables has slightly worse performance.  

\begin{table}[h!]
\centering
\normalsize % small, normalsize, large, Large
%\scalebox{1.00}{
\begin{tabular}[t]{lccccccc}
\toprule
  & IBS & C-index\\
\midrule
Cox all variables & 0.183 & 0.620\\
Cox backward & 0.183 & 0.615\\
Cox LASSO & 0.183 & 0.614\\
RSF & 0.182 & \textbf{0.622}\\
Neural Network 1h & \textbf{0.180} & -\\
Neural Network 2h & \textbf{0.180} & -\\
\bottomrule
\end{tabular} 
\caption{\label{tab:global_measures} Integrated Brier Score (IBS) and C-index on the test data. Neural network 1h and 2h refer to a neural network with one and two hidden layers respectively.}
%}
%\\[10pt]
%\caption*{Caption without a number}
\end{table}  

Stability of the networks was investigated by rerunning the same models on the test data, and showed that the NN with 1 hidden layer had stable predictive performance and variable importance. In contrast, the NN with 2 hidden layers was quite unstable regarding variable importance. This behavior might be related to the vast amount of weights that had to be trained for this model which can lead to overfitting (in total 26621 connection weights were estimated for a sample size of 41,530 patients in long format; whereas for the NN with 1 hidden layer 11136 connection weights). For the RSF, model obtained remarkable stability in terms of performance error after a particular number of trees (\texttt{ntree} = 300 was selected).

\subsection*{Prediction error curves}
Figure \ref{fig:brier_and_probs} shows the average prediction Brier error over time for all models. Small differences can be observed between Cox models and RSF. The NNs with 1 hidden and with 2 hidden layers have almost identical evolution over time achieving better performance than the Cox models and the RSF. 

\begin{figure}[h!]
\centering
\includegraphics[width=0.97\textwidth]{figure2.png}
\caption{\label{fig:brier_and_probs} \csentence{Prediction error curves for all models.}}
\end{figure}


\subsection*{Variable importance}
\begin{table}[h!]
\centering
\small
%\scalebox{0.97}{
\begin{tabular}[t]{lccc}
\toprule
  & Cox all variables & Cox backward & Cox LASSO\\
  & HR (95\% CI) &  HR (95\% CI) &  HR \\
\midrule
Re-transplantation & 1.602 (1.491-1.721) & 1.608 (1.501-1.722) & 1.558\\
Donor age & 1.010 (1.008-1.011) & 1.011 (1.009-1.012) & 1.009\\
Donor type DCD$^{(a)}$ & 1.483 (1.362-1.616) & 1.443 (1.338-1.556) &  1.298\\
log(Total cold ischemic time) & 1.258 (1.192-1.327) & 1.285 (1.221-1.353) & 1.191\\
Diabetes & 1.173 (1.125-1.225) & 1.176 (1.128-1.226) & 1.136\\
Race Black$^{(b)}$ & 1.240 (1.171, 1.314) & 1.261 (1.193-1.332) & 1.186\\
Life support & 1.343 (1.240-1.454) & 1.375 (1.272-1.487)  & 1.304\\
Recipient age  & 1.007 (1.005-1.009) & 1.008 (1.006-1.010) & 1.006\\
Incidental tumour & 1.314 (1.202, 1.437) & 1.315 (1.203-1.437) & 1.203\\
Hypertensive bleeding & 1.296 (1.185, 1.418) & 1.301 (1.190-1.423) & 1.214\\
HCV$^{(c)}$ serology status &  1.147 (1.091-1.206) & 1.148 (1.094-1.205) & 1.166\\
Pre-treatment status ICU$^{(d)}$ & 1.240 (1.143, 1.346) & 1.253 (1.160-1.354) & 1.164\\
\bottomrule
\end{tabular}
%} 
\\ \caption*{ \small{(a): Donor type DCD (Donor Circulatory Dead) vs DBD (Donor after Brain-Dead),  $\quad$ (b): Race Black vs White,
(c): Chronic hepatitis C virus,  $\quad$ (d): Intense Care Unit vs Non-hospitalised/Hospitalised}}

\caption{\label{tab:importance_cox} Hazard ratios along with their 95\% confidence intervals for the 12 most influential variables for the Cox models. Variables are presented in decreasing order according to the absolute z-score values (12.90 to 5.16) for the Cox model with all variables.}
\end{table}

In this section, the models are compared based on the most prognostic variables identified from the set of 97 predictors - 52 donor and 45 recipient characteristics. The hazard ratios for the 3 Cox models are shown in Table \ref{tab:importance_cox}. The strongest predictor is \textit{re-transplantation}. Having been transplanted before increases the hazard of graft-failure or death by more than 55\%.  The other most detrimental variables are \textit{donor age} and \textit{donor type circulatory dead}. One unit increase for donor age rises the hazard by around 1\% while having received the graft from a donor circulatory versus brain-dead increases the hazard by more than 29\% for all models. The rest of the factors which have an adverse effect are: \textit{cold ischemic time}, \textit{diabetes}, \textit{race}, \textit{life-support}, \textit{recipient age}, \textit{incidental tumour}, \textit{spontaneous hypertensive bleeding}, \textit{serology status of HCV} and \textit{intense care unit before the operation}.

\begin{table}[ht]
\centering
%\tiny
\scalebox{0.87}{
\begin{tabular}[t]{cccccc}
\toprule
 Neural network 1h & Rel-Imp & Neural network 2h & Rel-Imp & RSF & VIMP\\
\midrule
Re-transplantation & 0.035  & Re-transplantation & 0.028 & Donor age & 0.010\\
Life-support & 0.025 & HCV$^{(d)}$ serology status & 0.025 & Re-transplantation & 0.009\\
Pre-treatment status ICU$^{(a)}$ & 0.023 & Life-support & 0.024 & Life support & 0.007\\
Donor type DCD$^{(b)}$ & 0.023 & Donor age & 0.023 & HCV$^{(d)}$ serology status & 0.007\\
Race Black$^{(c)}$ & 0.022 & Diabetes & 0.021 & Pre-treatment status & 0.006\\
HCV$^{(d)}$ serology status & 0.022 & Pre-treatment status ICU$^{(a)}$  & 0.020 & Recipient age & 0.004\\
Diabetes & 0.020 & Working income & 0.020 & Aetiology & 0.003\\
Donor age & 0.020 & Race Black$^{(c)}$ & 0.019 & log(Last serum creatinine)  & 0.003\\
Working income & 0.018 & Previous abdominal surgery & 0.015 & Functional status & 0.002\\
Functional status Total assistance$^{(e)}$ & 0.017 & Donor pre-recovery diuretics & 0.015 & log(Total cold ischemic time) & 0.002\\
Aetiology HCV  & 0.017 & Aetiology Cholestatic & 0.011 & Race & 0.002\\
Hypertensive bleeding & 0.017 & Functional status Total assistance$^{(e)}$ & 0.015 & Diabetes & 0.002\\
\bottomrule
\end{tabular}
}
\\ \caption*{ \footnotesize{(a): Intense Care Unit vs Non-hospitalised/Hospitalised $\quad$ (b): Donor type DCD (Donor Circulatory Dead) vs DBD (Donor after Brain-Dead), (c): Race Black vs White,  $\quad$ (d): Chronic hepatitis C virus, $\quad$ (e): Total assistance vs No assistance}}

\caption{\label{tab:importance_nn_rsf} The 12 most prognostic factors for the neural networks with 1 and 2 hidden layers (Rel-Imp: relative importance) and for the Random Survival Forest (VIMP: variable importance). Note that the NN utilises time intervals as one of the input variables (check the contribution of time intervals in Table 1 of Additional file 1). For RSF importance is measured for each variable without distinction for each level. }
\end{table}

In Table \ref{tab:importance_nn_rsf} the most prognostic factors for the machine learning techniques are presented. The top predictors are provided in terms of relative importance (Rel-Imp) for the PLANN models and in terms of variable importance (VIMP) for the RSF. For the NNs, the strongest predictor is \textit{re-transplantation} (Rel-Imp 0.035 for 1 hidden and 0.028 for 2 hidden layers), which is the second strongest for the RSF (VIMP 0.009). According to the tuned RSF, the most prognostic factor for the overall graft-survival of the patient is \textit{donor age} (VIMP 0.010).  

Other strong prognostic variables for the NN with 1 hidden layer are \textit{life support} (Rel-Imp 0.025), \textit{intense care unit before the operation} (Rel-Imp 0.023) and \textit{donor type circulatory dead versus brain-dead} (Rel-Imp 0.023). For the NN with 2 hidden layers other very prognostic variables are \textit{serology status for HCV} (Rel-Imp 0.025), \textit{life support} (Rel-Imp 0.024) and donor age (Rel-Imp 0.023).

For the RSF \textit{life support} (VIMP 0.007), \textit{serology status for HCV} (VIMP 0.007) and \textit{intense care unit before the operation} (VIMP 0.006). Note that variable \textit{total cold ischemic time} which was identified as the 4th most prognostic for the Cox model with all variables and the 10th most prognostic for random survival forest is not in the list of the 12 most prognostic for both NNs.

\subsection*{Individual predictions}

\begin{figure}[h!]
\centering
\includegraphics [width=13.75cm, height = 13cm]{figure3.png}
\caption{\label{fig:predictions} \textbf{a)} Predicted survival probabilities for 3 new hypothetical patients using the Cox model with all variables (solid lines), the tuned RSF (short dashed lines) and the tuned NN with 1 hidden layer (long dashed lines). The green lines correspond to a reference patient with the median values for the continuous and the mode value for categorical variables. The patient in the orange line has diabetes (the other covariates as in reference patient). The patient in the red line has been transplanted before and has diabetes simultaneously (the other covariates as in reference patient). Values for 10 prognostic variables for the reference patient are provided in Table 2 of Additional file 1. \textbf{b)} Predicted survival probabilities for 3 patients selected from the test data based on the Cox model with all variables (solid lines), the tuned RSF (short dashed lines) and the tuned NN with 1 hidden layer (long dashed lines). Green lines correspond to a patient censored at 1.12 years. Patient in the orange line was censored at 6.86 years. Patient in the red line died at 0.12 years. Values for 10 prognostic variables for the patients are provided in Tables 3-5 of Additional file 1. 
}
\end{figure}

In this section, the predicted survival probabilities are compared for 3 new hypothetical patients and 3 patients from the test data. 

In Figure \ref{fig:predictions}a) the patient with reference characteristics shows the best survival. The highest probabilities are predicted by the RSF and the lowest by the Cox model. The same pattern occurs for the patient that suffers from diabetes (orange lines). The patient with diabetes who has been transplanted before has the worst survival predictions. In this case the NN predicts the highest survival probabilities and the Cox model built using all the prognostic factors the lowest. 

In Figure \ref{fig:predictions}b) the estimated survival probabilities are showed by the Cox model with all variables, the tuned RSF and the tuned PLANN with 1 hidden layer for 3 patients from the test set. The first patient shows the highest survival predictions by the 3 models. The RSF provides the highest survival probabilities and the NN the lowest. The second patient experiences lower survival probabilities (orange lines) whereas the third patient shows the lowest survival probabilities overall. For the second patient the NN predicts the lowest survival probabilities over time and for the third the Cox model.   

In general, the random survival forest provides the most optimistic survival probabilities whereas the most pessimistic survival probabilities are predicted by either the Cox model or the NN (more often by the Cox model).

\subsection*{Calibration}
Here 4 methods are compared: Cox model with all variables, RSF, PLANN 1 hidden and 2 hidden layers based on the calibration on the test data. For each method, the predicted survival probabilities at each year are estimated and the patient data are split into 10 equally sized groups based on the deciles of the probabilities. Then the survival probabilities along with their 95\% confidence intervals are calculated using the Kaplan-Meier methodology \cite{Kaplan1958NonparametricObservations}.

In figure \ref{fig:calibration_2y} the results are showed at 2 years since LT. The Cox model with all variables and the PLANN with 1 hidden layer are both well calibrated. The RSF and the PLANN with 2 hidden layers tend to overestimate the survival probabilities for the patients at higher risk. Survival neural network with 1 hidden layer seems to be the most reliable for predictions between the ML techniques. Calibration plots at 5 and 10 years can be found in Additional file 2.  

\begin{figure}[h!]
\centering
\includegraphics [width=13.5cm, height = 10cm]{figure4.png}
\caption{\label{fig:calibration_2y} Calibration plots at 2 years on the test data: \textbf{a)} Cox model with all variables, \textbf{b)} Random Survival Forest, \textbf{c)} Partial Logistic Artificial Neural Network with 1 hidden layer, \textbf{d)} Partial Logistic Artificial Neural Network with 2 hidden layers.
}
\end{figure}

\section*{Discussion}
With the rise of computational power and technology on the 21$^{st}$ century, more and more data have been collected in the medical field to identify trends and patterns which will allow building better allocation systems for patients, provide more accurate prognosis and diagnosis as well as more accurate identification of risk factors. During the past few years, machine learning (ML) has received increased attention in the medical area. For instance, in the area of LTs graft failure or primary non-function might be predicted at decision time with ML methodology \cite{Lau2017Machine-LearningTransplantation}. Briceño \textit{et al.} \cite{Briceno2014UseStudy} created a NN process for donor-recipient matching specifying a binary classification survival output (recipient or graft survival) to predict 3-month graft mortality.

In this study statistical and ML models were estimated for patients from the US post-transplantation. Random survival forest performed better than Cox models with respect to the C-index. This shows the ability of the model to discriminate between low and high risk groups of patients. The C-index was not estimated for NN because a natural ordering of subjects is not feasible. Therefore, the Brier score was measured each year for all methods. The Cox models / RSF showed results similar to the RSF having slightly smaller total prediction error (in terms of IBS). The NNs performed in general better than the Cox models or the RSF and had very similar performance over time. 

Special emphasis was given on the interpretability of the models. Results showed that Cox models (via hazard ratios) and the NNs with one/two hidden layer(s) (via relative importance) identified similar predictors. Both methods identified \textit{re-transplantation} as the strongest predictor and \textit{donor age}, \textit{diabetes}, \textit{life support} and \textit{race} as relatively strong predictors. According to RSF, the most prognostic variables were \textit{donor age}, \textit{re-transplantation}, \textit{life support} and \textit{serology status of HCV}. \textit{Aetiology} and \textit{last serum creatinine} were selected as the $7^{th}$ and the $8^{th}$ most prognostic. This raises a known concern about the RSF bias towards continuous variables and categorical variables with multiple levels \cite{Loh1997SPLITTREES} (\textit{aetiology} has 9 levels: metabolic, acute, alcoholic, cholestatic, HBV, HCV, malignant, other cirrhosis, other unknown). As continuous and multilevel variables incorporate larger amount of information than categorical, they tend to be favoured by the splitting rule of the forest during binary partitioning. Such bias was reflected in the variable importance results.  

When comparing statistical models with machine learning techniques with respect to interpretability, Cox models offer a straightforward interpretation through the hazard ratios. On the contrary, for both neural networks and random survival forests the sign of the prediction is not provided (if the effect is positive or negative). Additionally, for NNs interpretation is possible for different variable levels (with the method of Garson and its extension), whereas for RSF only the total effect of a variable is shown.

ML techniques are inherently based on mechanisms introducing randomisation and therefore very small changes are expected between different iterations of the same algorithm. To evaluate stability of performance, ML models were run several times under the same parametrisation. RSF were consistently stable after a certain number of trees (300 were selected). This was not the case for the NNs where instability is a common problem. It is challenging to tune a NN due to many hyper-parameter combinations available and the lack of a consistent global performance measure for survival data. IBS was used to tune the novel NNs, which may be the reason of instability for the NN with 2 hidden layers together with the large number of weights. Note also that the NN with 1 hidden layer is well calibrated whereas the NN with 2 hidden layers is less calibrated on the test data.

This is the first study where ML techniques are applied to transplant data where a comparison with the traditional Cox model was investigated. To construct the survival NN, the original problem had to be converted into a classification problem where exact survival times were transformed into (maximum) 10 time intervals denoting years since transplantation. On the other hand, for the Cox models and the RSF exact time to event was used. Recently, a new feed forward NN has been proposed for omics data which calculates directly a proportional hazards model as part of the output node using exact time information \cite{Ching2018Cox-nnet:Data}. A survival NN with exact times may lead to better predictive performance.  For UNOS data, 69.1\% of the recipients were alive/censored and 30.9\% had the event of interest. Results above were based on these particular percentages for censoring and events (for the NNs the percentages varied because of the reformulation of the problem). 

It might be useful to investigate how the number of variables affects the performance of the models. Here 97 variables were pre-selected supported by clinical and statistical reasons. It might be interesting to repeat the analyses on a smaller group of predictors, implementation time can be drastically reduced as the calculation complexity depends on sample size and predictors multiplicity. Alongside, predictive accuracy might be increased as some noisy factors will be removed from the dataset increasing the signal of potentially prognostic variables.

\section*{Conclusions}
In this work two alternatives to the Cox model from machine learning for medical data with large total sample size (over 60000 patients) and many predictors (97 in total) were discussed. 

RSF showed better performance than the Cox models with respect to C-index so it can be a useful tool for prioritisation of particular high risk patients. NNs showed better prediction performance in terms of Integrated Brier score. However, both ML techniques required a non-trivial implementation time. Cox models are preferable in terms of straightforward interpretation and fast implementation. Our study suggests that some caution is required when ML methods are applied to survival data. Both approaches can be used for exploratory and analysis purposes as long as the advantages and the disadvantages of the methods are presented.

\section*{List of abbreviations}
BS, Brier score;
CVPL, cross-validated log-partial likelihood; 
DCD, Donor Circulatory Dead; 
HBV, Chronic hepatitis B virus; 
HCV, Chronic hepatitis C virus;
IBS, Integrated Brier score; 
LASSO, least angle and selection operator; 
LT, liver transplantation; 
LUMC, Leiden University Medical Center;
ML, machine learning; 
NN(s), artificial neural network(s); 
OGS, overall graft-survival; 
OPO, Organ Procurement Organisations; 
OPTN, Organ Procurement and Transplantation Network; 
PLANN, partial logistic artificial neural network; 
PLANN-ARD, partial logistic artificial neural network - automatic relevance determination; 
Rel-Imp, relative importance; 
RSF, random survival forest; 
SM, statistical model; 
SRTR, Scientific Registry of Transplant Recipients; 
UNOS, United Network of Organ Sharing; 
VIMP, variable importance.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Backmatter begins here                   %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{backmatter}

\section*{Declarations}

\subsection*{Ethics approval and consent to participate} 
The ethics committee of Leiden University Medical Center (LUMC) approved the study. For all patients informed consent was provided to use the data for scientific research. 


\subsection*{Consent for publication}
The study was submitted to a functioning Institutional Review Board (IRB) for review and approval. Consent was provided for publication. 

\subsection*{Availability of data and materials}
The research data for this project is private. Unauthorized use is a violation of the terms of the Data Use Agreement with the U.S. Department of Health and Human Services. More information and instructions for researchers to request UNOS data can be found at \url{https://unos.org/data/}.
The R-code that was used to perform the analysis is available \href{https://github.com/GKantidakis/Survival-prediction-models-since-liver-transplantation/blob/master/analysis_code.R}{here}.


\subsection*{Competing interests}
The authors declare that they have no competing interests. The data reported here have been supplied by the Minneapolis Medical Research Foundation (MMRF) as the contractor for the Scientific Registry of Transplant Recipients (SRTR). The interpretation and reporting of these data are the responsibility of the author(s) and in no way should be seen as an official policy of or interpretation by the SRTR or the U.S. Government.

This study used data from the Scientific Registry of Transplant Recipients (SRTR). The SRTR data system includes data on all donor, wait-listed candidates, and transplant recipients in the US, submitted by the members of the Organ Procurement and Transplantation Network (OPTN). The Health Resources and Services Administration (HRSA), U.S. Department of Health and Human Services provides oversight to the activities of the OPTN and SRTR contractors.



\subsection*{Funding}
Georgios Kantidakis's work as a Fellow at EORTC Headquarters was supported by a grant from the EORTC Soft Tissue and Bone Sarcoma Group and Leiden University as well as from the EORTC Cancer Research Fund (ECRF). 


\subsection*{Authors' contributions}
JDB and AEB requested the data to the Scientific Registry of Transplant Recipients (SRTR) and provided clinical input. GK, HP, CL and MF designed the models. GK carried out the statistical analysis. GK wrote the manuscript and HP, MF critically revised it. All authors read and approved the final version.  

\subsection*{Acknowledgements}
The authors would like to thank the United Network of Organ Sharing (UNOS) and Scientific Registry of Transplant Recipients (SRTR) for providing the data about liver transplantation to Leiden University Medical Center (LUMC) under DUA number 9477.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  The Bibliography                       %%
%%                                                         %%
%%  Bmc_mathpys.bst  will be used to                       %%
%%  create a .BBL file for submission.                     %%
%%  After submission of the .TEX file,                     %%
%%  you will be prompted to submit your .BBL file.         %%
%%                                                         %%
%%                                                         %%
%%  Note that the displayed Bibliography will not          %%
%%  necessarily be rendered by Latex exactly as specified  %%
%%  in the online Instructions for Authors.                %%
%%                                                         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% if your bibliography is in bibtex format, use those commands:
\bibliographystyle{vancouver} % Style BST file (bmc-mathphys, vancouver, spbasic).
\bibliography{bibliography}      % Bibliography file (usually '*.bib' )

% for author-year bibliography (bmc-mathphys or spbasic)
% a) write to bib file (bmc-mathphys only)
% @settings{label, options="nameyear"}
% b) uncomment next line
%\nocite{label}

% or include bibliography directly:
% \begin{thebibliography}
% \bibitem{b1}
% \end{thebibliography}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Figures                       %%
%%                               %%
%% NB: this is for captions and  %%
%% Titles. All graphics must be  %%
%% submitted separately and NOT  %%
%% included in the Tex document  %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%% Do not use \listoffigures as most will included as separate files


%% Use of \listoftables is discouraged.
%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Additional Files              %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section*{Additional Files}
  \subsection*{Additional file 1 --- Supplementary material}
    Additional file 1 includes the Garson's algorithm for 2 hidden layers, a table with the relative importance of the time intervals for the neural networks with 1 and 2 hidden layes, criteria for variable pre-selection, a plot of survival and censoring distributions and 4 plots with individual patient characteristics.
  \subsection*{Additional file 2 --- Calibration plots at 5 and 10 years since LT.}
    Additional file 2 contains calibration plots at 5 and 10 years for a) a Cox model with all prognostic factors, b) a Random Survival Forest with all prognostic factors, c) a Partial Logistic Artificial Neural Network with 1 hidden layer with all prognostic factors and d) a Partial Logistic Artificial Neural Network with 2 hidden layers with all prognostic factors.


\end{backmatter}
\end{document}
