Random Survival Forest for Censored Functional Data

In a survival study, $N$ subjects are enrolled in a given period $\mathcal{T}=[a,b]$ called a follow-up. In this period, different data are collected for each $i$-th subject, with $i=1,\ldots,N$ . In general, the times $T_{i}$ denote the survival time, and $C_{i}$ denote the censored time, i.e., an event did not occur for the subject. Let $\delta_{i}$ be defined for each subject as follows:

During the follow-up study, $J_{i}$ values are collected for each subject until the event time, denoted by $T_{i}^{*}=\mbox{min}\,(C_{i},T_{i})$, occurs. This information defines the Survival Time Data (STD) given by $(\bm{Y}_{i},\bm{t}_{i},T_{i}^{*},\delta_{i})$, where $\bm{Y}_{i}=(Y_{11},...,Y_{1J_{i}})^{T}$ refers to the observed values vector with $Y_{ij}=Y(t_{ij})$ and $\bm{t}_{i}=(t_{i1},...,t_{iJ_{i}})^{T}$ is the time vector for the $i$-th unit.

In this context, the main challenge of employing the FDA approach and leveraging functional tools is to retrieve data for every subject, even when there are no observations during the entire observation period. For this reason, our approach involves the introduction of CFD, enabling the continuous reconstruction of information for each subject at any given moment. The main issue of this extension is how to replace these data by using different functions according to the number of recording values $J_{i}$. Note that $t_{iJ_{i}}\leq T_{i}^{*}\leq t_{iJ_{i+1}}$ for each $i$-th unit and no longitudinal measurements are available after $T_{i}^{*}$. Thus, the observation of the function $\mathcal{X}_{i}$ consist of $J_{i}$ pairs $(t_{ij},Y_{ij})$ for $i=1,...,N$ and $j=1,\ldots,J_{i}$.

Constructing functional data requires a different approach in survival analysis with censored data. Hence, we aim to develop a methodology that handles censoring and addresses the unique aspects of functional data, such as temporal truncation, to ensure accurate modelling and interpretation of the underlying patterns over time. In our case, we have $J_{i}\neq J$ as observed data, which are $J_{i}$-dimensional. For this reason, the first step is to convert the observed values $\{(Y_{ij},t_{ij}):i=1,...,N\>\text{and}\>j=1,...,J_{i}\}$ into a functional form by considering different recording values. In this connection, the standard approach to estimating the functional form by starting with observed values is the basis approximation. The CFD can be defined as follows:

As in (15), we can express these functions according to numbers of $J_{i}$. For which, using a fixed basis representation and other functional representations, we obtain the CFD as follows:

Moreover, $K_{i}$ is determined individually for $i$-th unit when the number of recorded values exceeds or equals $3$. It is determined using leave-one-out cross-validation. The method minimizes the prediction error of the number of components $K$ choosing the best $\hat{K}_{i}$ for the $i$-th subject as follows:

In addition, considering that in real-world applications the measurement of $Y_{ij}$ may be subject to errors, we consider the random noise $\varepsilon_{ij}$ with $E[\varepsilon_{ij}]=0$ and $var[\varepsilon_{ij}]=\sigma_{ij}^{2}$, where $\varepsilon_{ij}$ are independent across $i$ and $j$ and the errors are considered homoschedastic with $\sigma_{ij}^{2}=\sigma^{2}$. Then, we can express the observed values as $Y_{ij}=X_{i}(t_{ij})+\varepsilon_{ij}$, with $t_{ij}\in[a,T_{i}^{*}]$, max$\{T_{i}^{*}:i=1,...,N\}\leq\tau$, where $\tau$ is the length of the study follow-up. Hence, for different numbers $J_{i}\geq 2$ using the definition (16) we obtain the follow model:

In order to evaluate the coefficients for the (18), we use the SSD criterion as follows:

calculated for all subjects $i=1,...,N$.

FPCA is an extension of the classical PCA because the principle is to replace vectors with functions, matrices by linear operators, and, in particular, covariance by auto-covariance operators Ramsay and Silverman (2005). Moreover, scalar products in vector are replaced by scalar products in function space $L^{2}$. Let’s suppose that $X(t)$ is a generic trajectory defined on Hilbert space $L^{2}(\mathcal{T})$ with its mean function $\mu(t)=E[X(t)]$ and its covariance function $G(s,t)=Cov(X(t),X(s))$, with $t,s\in\mathcal{T}$. The covariance function can be expressed by the spectral decomposition as $G(s,t)=\sum_{m=1}^{\infty}\lambda_{m}\xi_{m}(t)\xi_{m}(s)$, where $\{\lambda_{m}\}_{m=1,..,\infty}$, correspond to a set non-increasing eigenvalues such that $\sum_{m}\lambda_{m}<\infty$ and $\{\xi_{m}(t)\}_{m=1,..,\infty}$ correspond to the eigenfunctions. Thus, the trajectory $X(t)$ admits the following Karhunen-Loevè expansion:

where the coefficient $\nu_{m}=\int_{\mathcal{T}}[X(t)-\mu(t)]\xi_{m}(t)\,dt$ corresponds to the $m$-th Functional Principal Component (FPC) score of $X(\cdot)$ and satisfies the condition $E[\nu_{m},\nu_{n}]=\delta_{m\,n}n_{m}$ with $\delta_{m\,n}=1$ if $m=n$ and $0$ otherwise.

For the case of sparse or irregular functional data, a fully non-parametric approach denoted as PACE Yao et al. (2005) is used. Let $Y_{ij^{(h)}}$ be the $j$-hth observation obtained from CFD (16) evaluating the $i$-th curve at time $t_{ij^{(h)}}\in[a,T_{i}^{*}]$, where $h=t_{ij^{(h)}}-t_{i{j-1}^{(h)}}$ is the step size chosen for the time increment by considering $t_{i1}=t_{i1^{(h)}}\leq t_{i2^{(h)}}\leq...\leq t_{iJ_{i-1}^{(h)}}\leq t_{iJ_{i}^{(h)}}=t_{iJ_{i}}$ for $i=1,...,N$. Then, including additional measurement errors $\varepsilon_{ij}$, uncorrelated among them with $E[\varepsilon_{ij^{(h)}}]=0$ and $Var[\varepsilon_{ij^{(h)}}]=\sigma^{2}$, to a model with measurements $Y_{ij^{(h)}}$ calculated at time $t_{ij^{(h)}}$ from (20) we obtain:

where the expansion (21) can be approximated as follows:

Moreover, the mean function, the covariance function and the eigenfunctions are calculated using the local linear smoothers Fan and Gijbels (1992). From (21), we note that the $Cov(Y_{ij^{(h)}},Y_{il^{(h)}}|\,t_{ij^{(h)}},t_{il^{(h)}}=Cov(X(t_{ij^{(h)}},X(t_{il^{(h)}})+\sigma^{2}\delta_{jl}$ for which $G_{i}(t_{ij^{(h)}},t_{il^{(h)}})=(Y_{ij^{(h)}}-\hat{\mu}(t_{ij^{(h)}})\,(Y_{il^{(h)}}-\hat{\mu}(t_{il^{(h)}})$, where $\hat{\mu}(t)$ is a mean function obtained by applying a local linear smoother to the scatterplot $\{(t_{ij^{(h)}},Y_{ij^{(h)}}):1\leq i\leq N,\,1\leq j\leq J_{i}^{(h)}\}$. The estimated covariance surface, $\hat{G}(s,t)$, is calculated considering the raw estimates $G_{i}(t_{ij^{(h)}},t_{il^{(h)}})$ and then applying two-dimensional smoothing to the scatterplot $\{(\,(t_{ij^{(h)}},t_{il^{(h)}})\,;\,G_{i}(t_{ij^{(h)}},t_{il^{(h)}})\,):1\leq i\leq N,\,1\leq j\neq l\leq J_{i}^{(h)}\}$. Nevertheless, the calculation of scores $\nu_{im}$ is done according to the definition of integrals in $L^{2}(\mathcal{T})$, but in this case having that $Y_{ij^{(h)}}$ are available only at discrete random times $t_{ij^{(h)}}$ we can obtain the estimates $\hat{\nu}_{im}$ from (21) as $\hat{\nu}_{im}=\sum_{j=1}^{J_{i}^{h}}(Y_{ij^{(h)}}-\hat{\mu}(t_{ij^{(h)}})\,\hat{\xi}_{m}(t_{ij^{(h)}})\,(t_{ij^{(h)}}-t_{ij-1^{(h)}})$. The goal is to obtain predicted trajectories $\hat{X}_{i}(t)$ for irregular data $Y_{ij^{(h)}}$. In summary, by discretising continuous trajectories with step sizes $h$, we can effectively analyse functional data even when observations are sparse or irregular. This approach allows for the practical implementation of FPCA, capturing essential features of the data while accounting for measurement errors and irregular observation patterns.

An alternative to the Riemann sum given by the previous formula is to assume that in (21), the $\nu_{im}$ and $\varepsilon_{ij^{(h)}}$ are jointly Gaussian. For this reason, defining $\bm{\hat{Y}}_{i}^{(h)}=(Y_{i1^{(h)}},...,Y_{iJ_{i}^{(h)}})^{T}$, $\bm{\mu}_{i}^{(h)}=(\mu(t_{i1^{(h)}}),...,\mu(t_{iJ_{i}^{(h)}}))^{T}$ and $\bm{\xi}_{im}^{(h)}=(\xi_{m}(t_{i1^{(h)}}),...,\xi_{m}(t_{iJ_{i}^{(h)}}))^{T}$ under Gaussian assumptions, we obtain the best prediction for the $m$-th FPC score $\nu_{im}$ for the $i$-th subject as conditional expectation $E[\nu_{im}|\bm{\hat{Y}}_{i}^{(h)}]$ as follows:

where $\bm{\Sigma}_{\bm{\hat{Y}}_{i}^{(h)}}=\text{Cov}(\bm{\hat{Y}}_{i}^{(h)},\bm{\hat{Y}}_{i}^{(h)})=\text{Cov}(\bm{\hat{X}}_{i}^{(h)},\bm{\hat{X}}_{i}^{(h)})+\sigma^{2}\bm{I}$, with $\bm{\hat{X}}_{i}^{(h)}=(\bm{X}(t_{i1^{(h)}}),...,\bm{X}(t_{iJ_{i}^{(h)}}))$ and $\bm{I}\in\mathbb{R}^{J_{i}^{(h)}\times J_{i}^{(h)}}$ identity matrix. Thus, $(\bm{\Sigma}_{\bm{\hat{Y}}_{i}^{(h)}})_{jl}=G(t_{ij^{(h)}},t_{il^{(h)}})+\sigma^{2}\delta_{jl}$, where $\delta_{jl}=1$ if $j=l$ and $0$ otherwise. Equation (23) can be obtained by using the estimates of $\lambda_{m}$, $\bm{\xi}_{im}^{(h)}$, and $\bm{\Sigma}_{\bm{\hat{Y}}_{i}^{(h)}}$ as follows:

Finally, from (24) the prediction for the trajectory $X_{i}(t)$ for $i$th subject can be done considering the first $p$ eigenfunctions as follows:

In supervised learning, Survival Trees (STs) extend the classical Classification Trees (CTs) grown on censored data. Unlike classical decision trees, STs are not directly applicable to censored data; a ST focuses on predicting the time until an event of interest occurs for a subject. These problems introduce a distinctive challenge when dealing with survival data, making certain aspects of implementing ST more intricate than decision trees for classification tasks. In survival studies, we have, for each subject, an observation vector $\bm{Y}_{i}=(Y_{i1},...,Y_{iJ_{i}})^{T}$ that is a data collected in $J_{i}$ values until the event time $T_{i}^{*}=\mbox{min}\,(C_{i},T_{i})$ occurs for the $i$-th unit and the response vector $\bm{Z}=(Z_{1},Z_{2},...,Z_{N})^{T}$ that can have binary values given by $Z_{i}=0$ or $Z_{i}=1$ for $i=1,...,N$. Thus, we observe, for the $i$-th subject, the data given by $(\bm{Y}_{i},\bm{t}_{i},Z_{i},T_{i}^{*},\delta_{i})$, with $\delta_{i}$ censoring indicator.

In traditional CTs, the non-pruned trees lead to identifying homogeneous terminal nodes (leaves) using some splitting criteria based on decreasing an impurity metric, e.g. the Gini or Shannon-Weier indexes. However, in STs, the splitting criterion used in each node $s$ of the tree, denoted by $T$, differs from classical methods when dealing with censored data. Indeed, iteratively, the procedure grows a tree until each node defines, at the very least, a singular distinct event by using survival splitting as a means for maximizing between-node survival differences given by curves defined through parametric (Cox regression) and non-parametric (Kaplan-Meier) approaches Shimokawa et al. (2015).

In a classical decision tree, the features vectors $\bm{X}\in\mathbb{R}^{J}$ have the same size for each unit: on the contrary, in STs the challenge is how to grow a tree having, for each subject, data collected during follow-up study $\mathcal{T}$ given by $Y_{ij}=Y(t_{ij})$ for $j=1,...,J_{i}$ values with $J_{i}\neq J_{k}$ for some $i\neq k$ with $t_{ij}\in\mathcal{T}$. For this reason, starting from these data types and defining the CFD, we deal with new features using the values estimated on CFD and applying the PACE decomposition on them. Therefore, extending Functional Classification Trees Maturo and Verde (2022) to the context of survival analysis, we propose the Functional Survival Classification Tree with Principal Components (FSCTs-FPCs). The latter classifier considers the new data $(\bm{\hat{\nu}}_{i},Z_{i},T_{i}^{*},\delta_{i})$, where $\bm{\hat{\nu}}_{i}=(\hat{\nu}_{i1},\hat{\nu}_{i2},...,\hat{\nu}_{ip})^{T}$ corresponds to the $i$-th score vector obtained from Equation 24. Hence, the features matrix is defined as follows:

Typically, considering the node $s$ of the survival tree $T$, given a score-predictor $\bm{\hat{\nu}}_{m}$, the method finds a value $c$ and selects the best $\bm{\hat{\nu}}_{0}$ such that survival differences between the two inequalities, $\bm{\hat{\nu}}_{0}\leq c$ and $\bm{\hat{\nu}}_{0}>c$, are maximized in order to predict the response $\bm{Z}$. The aim is to define in node $s$ two regions $R_{1}(m,c)=\{\hat{X}(t)\in L^{2}(\mathcal{T})|\bm{\hat{\nu}}_{m}\leq c\}$ and $R_{2}(m,c)=\{\hat{X}(t)\in L^{2}(\mathcal{T})|\bm{\hat{\nu}}_{m}>c\}$ to identify $l=1,...,L$ times for the two groups, with $t_{1}\leq t_{2}\leq...\leq t_{L}$. One of the survival criteria used within node $s$, to select the threshold $c$ and the functional score-prediction $\bm{\hat{\nu}}_{0}$, is the log-rank test as splitting method Ziegler et al. (2007). It can be used to maximize a splitting value $c$ for a predictor variable $\bm{\hat{\nu}}$ as follows:

with $d_{lj}$ number of events in daughter nodes $j\in\{1,2\}$ at time $t_{l}\in\mathcal{T}$ (time of observation), where $d_{l}=d_{l1}+d_{l2}$ is the total number of events within node $s$. Moreover, $r_{lj}$ is the number of individuals who are at risk (alive) in daughter nodes $j\in\{1,2\}$ at time $t_{l}\in\mathcal{T}$, where $r_{l}=r_{l1}+r_{2}$ is the total number of individuals alive within node $s$. $K$ is the total number of times between two daughter nodes. The number $r_{l1}$ corresponds to $T_{l}\geq t_{l}$ for which $\bm{\hat{\nu}}\leq c$, while $r_{l2}$ corresponds to $T_{l}\geq t_{l}$ for which $\bm{\hat{\nu}}>c$ within node $s$. The Equation (27) allows finding the best $\bm{\hat{\nu}}^{*}$ and the best $c^{*}$ which $|L(\bm{\hat{\nu}}^{*},c^{*})|\geq|L(\bm{\hat{\nu}},c)|$ for each $\bm{\hat{\nu}}$ and $c$ chosen randomly. This process is repeated randomly at every node $s$ until the terminal node is reached. The FSCTs-FPCs use different methods for constructing predictive models. Two common alternative approaches include log-rank splitting, which divides data based on survival differences, and the Nelson-Aalen estimator, used for estimating cumulative hazard. The last method helps us to identify significant predictors of survival outcomes in the terminal nodes. Once the survival tree $T$ is built, let $s$ be the terminal node, where $s=1,...,S$, in which are defined the times $t_{1s}\leq t_{2s}\leq...\leq t_{L_{s}s}$ and let $d_{ls}$ and $r_{ls}$ the number of events and individuals at risk at time $t_{ls}$, the survival predictor in $T$ is defined in terms of the predictor within each terminal node $s$. For this, the Cumulative Hazard Function (CHF) and the Survival Function (SF) are obtained by using the Nelson-Aalen estimator and Kaplan-Meier estimator defined as follows:

where $\mathcal{H}_{s}(t)$ and $\hat{S}_{s}(t)$ are the hazard and the survival estimate for a terminal node $s$ obtained considering all distinct event times $t_{ls}\leq t$. However, predicting survival outcomes after the survival tree construction is essential. This evaluation often involves comparing the predicted survival probabilities with the observed outcomes. One standard method for evaluating the importance of a survival tree is using the Concordance Index (C-index), also known as Harrell’s C or the C-statistics, which tells us how accurately the model ranks people’s survival times. A higher C-index indicates better predictive performance of the survival tree model. Other metrics, such as the Brier score, can be used for further evaluation.

The Functional Random Survival Forest (FRSF) extends the classic RSF. In our approach, the method is built for the CFD considering the pairs $\{Z_{i},X_{i}(t)\}$, where $X_{i}(t)$ is a CFD curve for $i$-th unit defined in the metric space $(L^{2}(\mathcal{T}),\|\cdot\|_{2})$, while $Z_{i}$ is a scalar response. Starting from the dataset $D=\{Z_{i},\hat{\bm{\nu}}_{i},T_{i}^{*},\delta_{i}\}_{i=1}^{N}$ used for the realization of the model, where $\hat{\bm{\nu}}_{i}\in\mathbb{R}^{p}$ is the $i$-th score calculated by PACE on CFD, $Z_{i}$ is a binary response, $T_{i}^{*}=\min\,(C_{i},T_{i})$, and $\delta_{i}$ is a censoring indicator at sample $i=1,...,N$. The method generates $B$ bootstrap samples from $D$ by creating the Training data set (In-bag data) and Validation data set (Out-of-bag data) to construct multiple STs.

The following pseudo-code can summarise the FRSF:

The step $8$ of the algorithm can be explained from (28) as follows:

where ${\mathcal{H}}_{b}(t|\hat{\bm{\nu}})$ is the in-bag CHF, while ${\mathcal{H}}_{b}(t|\hat{\bm{\nu}}_{i})$ is the out-of-bag CHF for $b$-th survival tree.

The metrics derived from the ensemble CHF, such as the predictor error and the variable importance (VIMP), provide crucial insights into the model’s goodness-of-fit calculated using only the OOB data. This process evaluates the accuracy of predictors and determines their relative importance in predicting survival outcomes. In FRSF the Breiman-Cutler Variable Importance (VIMP), also known as permutation importance Breiman (2002) is used to quantify the importance of each predictor variable by evaluating its impact on prediction error. The VIMP method permutes the out-of-bag (OOB) values of $\hat{\bm{\nu}}$ within each tree. This permutation process entails randomising the values of $\hat{\bm{\nu}}$, while keeping the values of the other predictor variables unchanged. Moreover, the modified data influences decision paths through the tree structure based on the altered $\hat{\bm{\nu}}$ values. The resulting out-of-bag error, indicative of prediction accuracy with the permuted $\hat{\bm{\nu}}$ is compared to the original error calculated with unaltered $\hat{\bm{\nu}}$ values. The discrepancy between these errors quantifies the importance of $\hat{\bm{\nu}}$ in predicting outcomes. Aggregating these individual importance scores across all trees in the forest yields the permutation importance for variable $\hat{\bm{\nu}}$, providing insights into its contribution to the overall predictive performance of the model.