text:
\section{Introduction}
Survival analysis for time-to-event data is used in a variety of disciplines, particularly in medical research. The Cox proportional hazard, accelerated failure time, and proportional odds models and their extensions are popular modeling tools to relate survival times with predictors or covariates in the presence of censorship. Though these are popular models, they are constrained by parametric assumptions such as proportional hazards, constant odds, linearity, etc. Survival trees and forests are popular non-parametric alternatives to these semi-parametric models which have the ability to identify non-linear interactions. One of the main benefits of tree models is their easy interpretation. Moreover, they can cluster subjects automatically according to their survival behavior based on covariates to identify prognostic groups.

Tree-based survival analysis was introduced in the 1980s \citep{leblanc1993survival,bou2011review} primarily through the extension of classification and regression tree (CART) models \citep{breiman1984classification} to censored survival data \citep{ciampi1981approach,gordon1985tree}. The major difference among these papers are the use of different splitting criteria such as log-rank statistics, likelihood ratio statistics, Wicoxon-Gehan statistics, node deviance measures, and different impurity measures. A crucial issue of these tree models is the selection of a single tree of suitable size using forward or backward algorithms.

In contrast to ordinary learning approaches which  construct one learner from training data, ensemble methods construct a set of learners known as base learners and combine them for prediction. Base learners are usually generated from the training data by a standard algorithm like decision tree, neural networks or other kinds of learning algorithms. There are two major kinds of ensemble methods based on the method of construction of the base learners. The first type is sequential ensemble methods where the base learners are generated sequentially and  the dependence among these base learners are used efficiently. The major example of this type  is the popular AdaBoost \citep{freund1997decision}.
The alternative type of ensemble learning  constructs the base learners in a parallel way with bagging \citep{breiman1996bagging}  as an example.

Recently, ensemble trees have become a popular alternative as a predictive survival model building tool \citep{ishwaran2004relative,hothorn2005survival,hothorn2006unbiased,ishwaran2008random,zhou2015rotation} using the concepts of bagging and random forests \citep{breiman1996bagging,breiman2001random}. In terms of application, \citet{garg2011phase} develop an approach to cluster hospital stays using survival trees and \citet{loh2015regression} use regression trees to identify subgroups with differential treatment effects. Indeed, there are a variety of other tree models which have been developed for use in survival and reliability models. \citet{levine2014bayesian} present some methodology for modeling tree survival functions with correlated observations. \citet{barrett2013gaussian} use an empirical Bayes approach by using Gaussian process regression on transformed survival times.  \citet{martino2011approximate,joensuu2012risk} and \citet{fernandez2016gaussian} use Gaussian processes in the survival model framework, and \citet{negassa2005tree}, \citet{su2008interaction} apply tree models for subgroup identification, but all retain the proportional hazards assumption.

Bayesian analysis of survival data is popular due to its flexibility to include prior information, suitability to develop complex models, and ability to construct efficient algorithms for full scale uncertainty quantification \citep{ibrahim2014b,sinha1997semiparametric}.   The Bayesian  classification and regression tree methods were developed by \citet{chipman1998bayesian,denison1998bayesian}.
\citet{clarke2008bayesian} use empirical Bayes priors and ensembles of Weibull tree models to model survival distributions.
The majority of Bayesian tree-based models in survival analysis are through ensembles of trees techniques (in a boosting way) using Bayesian additive regression trees (BART) \citep{chipman2010bart}. \citet{bonato2010bayesian} utilize BART to estimate latent parameters in a gamma process \citep{kalbfleisch1978non}, Weibull, and accelerated failure time models; \citet{henderson2017individualized} use BART in conjunction with a mean-constrained Dirichlet process mixture model to analyze heterogeneous treatment effects, which is reminiscent of \citet{kuo1997bayesian}; \citet{sparapani2016nonparametric} model survival times using a BART probit regression approach.

Even though machine learning algorithms usually focus on the predictive performance, the research in the field of interpretable machine learning \citep{molnar2020interpretable,doshi2017towards} has flourished in recent years due to its direct connection with artificial intelligence (AI) \citep{kuziemsky2019role}. For instance, AI is transforming the industry of healthcare rapidly and telehealth is one of its newest products. Telehealth utilizes digital information and communication technologies to access and manage healthcare services remotely. AI is booming in this industry of telemedicine where an efficient system for doctors is created to analyze, screen, and diagnose different disease conditions. However, a black-box machine-learning based predictive tool is not sufficient  and a more interpretable decision system is needed for this purpose. Hence, most of the ensemble methods (based on boosting or bagging logics) are not suitable in this framework.

To develop interpretable as well as flexible machine learning models for survival data,  we resort to a single tree model. One purpose of the proposed methodology is to provide a Bayesian tree framework which is simple to interpret (i.e.\ obtaining one tree) and flexible, i.e., moves beyond the proportional hazards model \citep{cox1972regression}. The  hazard function at each terminal node of the tree is modeled flexibly using an exponentiated Gaussian process. The posterior space of trees is efficiently searched via Markov chain Monte Carlo (MCMC) by obtaining the marginal likelihood of the observed data via a Laplace approximation \citep{tierney1986accurate}. The combination of interpretability coupled with flexibility makes this model well-suited for many problems, including subgroup identification and identifying prognostic and predictive biomarkers in medical research. In this Bayesian tree framework we sample from the posterior distribution of trees which provides a sample of multiple trees for decision making.

Most existing Bayesian tree survival models have transformed the survival problem into a regression or classification framework. For example,  \citet{bonato2010bayesian} use latent variables and parametric models. Tree ensembles are used primarily on the covariate structure in hierarchical specifications through a regression structure. \citet{sparapani2016nonparametric} transform survival responses into binary outcomes and use BART in this classification setup.  In contrast, we model the survival responses directly using  different hazard functions at each terminal node, and  construct the partition of the covariate space using a single tree structure similar to that of the original Bayesian classification and regression tree (CART) papers \citep{chipman1998bayesian,denison1998bayesian}.   However, our proposed model is different than existing CART models where simple parametric models (like Gaussian or multinomial models) with finite numbers of parameters and conjugate priors have been used at the terminal nodes of the tree. In contrast, we allow an unknown hazard function to model the  survival data within each terminal node of the tree.  An exponentiated Gaussian process is used as a prior for this hazard function which is  flexible but creates a non-conjugate posterior distribution. That way, explicit marginalization is not possible to obtain the marginal likelihood function which is the key step of the existing Bayesian Tree search algorithms for CART as well as BART  models. We propose a Laplace-approximation based method to develop a novel tree search algorithm for non-conjugate Bayesian tree models. In general, this new approach will aid the Bayesian tree models to be applied in more complex modeling scenarios.

Moreover, an important contribution of this paper is to provide theoretical properties of these Bayesian tree models. We  establish  consistency properties for our estimator under the tree partition model. Posterior consistency for survival functions in the Bayesian non-parametric setup  has been   investigated by modeling the  hazard  functions by completely random measures. Some earlier work using Dirichlet process priors can be found in  \citet{ghosh1995consistency}, where weak consistency was shown.  Later, using neutral to right processes such as the Beta and Gamma processes, weak consistency was shown in \citet{kim2003bayesian}. \citet{de2009asymptotics} show similar consistency results  where the hazard rates are a  kernel mixture of completely random measures. Here, we provide a framework for modeling the log hazard by a Gaussian process over partitions, where we can select the smoothness parameter based on observations for each partition. Under flexible conditions we show weak consistency and under slightly stronger conditions we establish the strong or $L^1$ consistency.  Partition modeling has been used for survival functions
\citep{ishwaran2008random} for a discrete nominal variable where $L^1$ consistency for survival function over a compact set  was shown for Kaplan-Meier type estimators \citep{ishwaran2010consistency}. Later,   developments using non i.i.d assumptions can be found in  \citet{cui2017consistency}. In a  different context, theoretical results have been provided for BART mainly in a regression setup (\citet{linero2018,Rock2017,Rock2018}). However, we are the first one to provide consistency theory for a Bayesian tree survival model, while modeling the sudden changes in survival function over the covariate space.

The paper is organized as follows:  Section 2 describes the model, Section 3 shows the theoretical properties of the model, Section 4 discusses the MCMC algorithm, Section 5 shows some simulated examples and an application to a liver survival dataset, Section 6 concludes.
\section{Survival models}
Let $t$ be a continuous nonnegative random variable representing the survival time of the subjects with the density function $f(t)$, the survival function $S(t)$ and the hazard function $h(t)= f(t) / S(t)$. In general, the hazard function depends on the survival time and the vector of covariates $\bx=(x_{1},x_{2},\ldots,x_{p})$. The Cox proportional hazards model \citep{cox1972regression} allow time and covariates to affect the hazard function separately, $h(t \mid \bx)=h_{0}(t){\rm exp}(\bx\bbeta)$, where $h_{0}$ is the baseline hazard and $\bbeta$ is a vector of regression parameters. This model has a strong assumption that the ratio of hazards for two individuals depends on the difference between their linear predictors at any time.

We extend this class of models by developing the Treed Hazards Model (THM) which allows a nonlinear effect of survival time and covariates in a flexible and adaptive way. In the THM model we split the covariate space of $\bx$ with tree partitions and a fit piecewise constant hazard model within each partition element. Let's denote the tree as $T$ with $M$ terminal nodes and the vector of  hazards within tree partition elements (terminal nodes) as $\bh=(h_{1}(\cdot),h_{2}(\cdot),\cdots,h_{M}(\cdot))$ where $h_{i}(\cdot)$ is the hazard
function at the $i$th terminal node.  This is in contrast to specifying a finite number of parameters at each terminal node as in existing Bayesian CART models. The THM model describes the conditional distribution of $t$ given $\bx$ in this tree structure. If $t$ lies in the region corresponding to the $i$th terminal node then $t \mid \bx$ has the distribution $f(t\mid h_{i}(\cdot))$. In the next section we develop the hierarchical tree model for the survival data.

\subsection{The Bayesian hierarchical treed hazards model }
Let $\by = \{y_1,\ldots,y_n\}$ be the survival time for $n$ subjects (realizations of $t$),  $\bdelta = \{\delta_1, \ldots, \delta_n\}$ be the corresponding censoring indicators  where $\delta_j$ is 1 if $y_j$ is observed and 0 if $y_j$ is right censored. Then, we can write the posterior distribution of the  model unknowns $T$ and $\bh$ as $p(T,\bh \mid \by,\bx,\bdelta)\propto p(\by, \bdelta \mid T,\bh)p(\bh \mid T)p(T \mid \bx)$.  In the following sections, we specify the likelihood function $p(\by, \bdelta \mid T,\bh)$ and the prior distributions $p(\bh \mid T)$, $p(T \mid \bx)$ which will be used to obtain the posterior distribution.

\subsubsection{Tree partition and the prior distribution $p(T \mid \bx)$}
First, we describe the tree partition $T$ on the covariate space and the related prior distribution. We assume that the data can be partitioned into $M$ partition elements $\by_1,\ldots,\by_M$ such that $\bigcup_{i=1}^M \by_i = \by$ and $\by_i \cap \by_{i^\prime} = \emptyset$ for all $i \neq i^\prime$.  We further assume that the data within the $i$th partition element are conditionally independent given the hazard function, $h_i(t)$, and assume that the $M$ hazard functions are mutually independent.  The tree partition is constructed through a series of recursive binary splits of the covariate space \citet{chipman1998bayesian}. The binary tree $T$ subdivides the predictor space marginally by choosing one of the $p$ covariates from $\bx$ and then defining a partition rule.

We adopt the same prior as \citet{chipman1998bayesian} which is constructed by splitting the tree nodes with a designated probability and assigning a rule each time a node is split.  Let $S_\ip$ equal 1 if the $\ip$th node is split, otherwise 0, $\ip=1,\ldots,2M-1$ where $M$ is the number of terminal nodes on a tree.  Let $R_\jp = (v_\jp,r_\jp)$, $\jp=1,\ldots,M - 1$ represent the rules of the $M -1$ internal nodes.  Here, $v_\jp \in \{1,\ldots,p\}$ indicates which variable is to be split, and $r_\jp$ is the splitting rule on the $v_\jp$th variable.  The partition rule for continuous variables is formed by choosing a threshold and then grouping observations based on whether or not their covariate value is less than or equal to, or greater than, the splitting threshold.  For categorical variables, the split rule is defined as a subset of the unique categories for that variable, and observations are grouped based on whether their covariate value is in the subset or not.

Let $T = (S_1,\ldots,S_{2M -1},R_1,\ldots,R_{M-1})$ be the collection of split indicators and rules for a particular tree. Then the prior for a tree, $T$, can be evaluated as
\begin{eqnarray} \label{eq:prior}
\pi(T) = \prod_{\ip=1}^{2M-1} \pi(S_\ip) \prod_{\jp=1}^{M-1}\pi(R_\jp).
\end{eqnarray}

Following \citet{chipman1998bayesian}, we define $\pi(S_\ip = 1) = 1 - Pr(S_\ip = 0) = \gamma(1 + d_i)^{-\theta}$ where $d_\ip$ represents the depth of the $\ip$th node (i.e.\ the root node has a depth of 0, child nodes of the root node have depth 1, etc.).  Here, $\gamma$ and $\theta$ are hyperparameters chosen by the modeller and influence the prior for the general size and shape of the tree. Here, $\pi(R_\jp) = 1/(m_\jp m^\star_\jp)$ where $m_\jp$ represents the number of variables on the corresponding node that are available for splitting and $m^\star_\jp$ is the number of splitting rules available on the selected variable for the node to which the rule belongs.  The splitting rules at a given node are selected from the observed data which reach that node from its ancestry.

\subsubsection{The Likelihood function $p(\by, \bdelta \mid T,\bh)$}
 The partition model assumes that given a tree, $T$, the hazard functions at each terminal node are independent and the data within a terminal node, given the hazard function at the node, are independent and identically distributed.  It is important to note that covariates are used to construct the tree, $T$, but are not used in the hazard functions.  Consequently, the likelihood function, given a tree, can be expressed as a simple product of likelihoods
\begin{equation}
p(\by, \bdelta \mid T, \bh) = \prod_{i=1}^{M} p(\by_i, \bdelta_i \mid h_i(\cdot)) \nonumber
\end{equation}
where $\by_i,\ \bdelta_i$ represent the $n_i$ responses and censoring indicators which belong to the $i$th terminal node, $p(\by_i, \bdelta_i \mid h_i(\cdot))$ represents the likelihood of $\by_i,\ \bdelta_i$ given the hazard function in the $i$th terminal node and $M$ is the total number of terminal nodes.

We write the density of an uncensored observation at the $i$th terminal node as $p_i(t)=h_i(t)e^{-\int_0^th_i(s)ds}$ for hazard function $h_i(\cdot)$. Let $p_c(\cdot)$ be the density of the censoring distribution, which has the same support as the true  densities $p_{0,i}(\cdot),\ i = 1, \ldots, m^*$ corresponding to true partition of the covariate space.

We place a Gaussian process prior on the log-hazard rate function on $[0,T_n]$, where $T_n=\max\{y_i\}$, $\log h_{i}(\cdot)=\beta_i+w_i(\cdot)$ and $w_i(t)$ follows Gaussian process with covariance  function $\tau_i^2K_{l_i}(s,t)=\tau_i^2\exp\{-|s-t|/l_i\}$, and let $\pi(\tau_i^2,l_i)=\pi(\tau_i^2)\pi(l_i)$ be the prior on the hyperparameters, and $\pi(\beta_i)$ be the prior on $\beta_i$.

\subsubsection{The hazards model $p(\bh \mid T)$ and the approximate likelihood}

In order to run a tractable Markov chain, a couple of alterations to the above model are used in practice.  First, an approximate likelihood and prior are used  where the domain is decomposed into small-width bins. Also, we normalize the domain between $[0,1]$ and use a mean zero Gaussian process under the proposed covariance kernel.

The approximation assumes the hazard functions are piecewise constant within each bin of the form $h_i(t) = \lambda_{ik}$ for $t \in (s_{k-1}, s_k]$ with $z_k$ denoting the midpoint of the interval.  In all, there are $K$ bins with endpoints $(s_{k-1},s_k]$ for $k=1,\ldots,K$ where $s_0 =0$ and $s_K = 1$.  For simplicity, we use the same grid for each hazard function and assume that the bins all have equal width, inducing a regular grid for $(z_1,\ldots,z_K)$.   The likelihood for data at a single terminal node under these assumptions can be expressed as
\begin{eqnarray*}
p(\by_i, \bdelta_i \mid \lambda_{i1},\ldots,\lambda_{iK}) &=& \prod_{j=1}^{n_i} \lambda_{ik_{ij}}^{\delta_{ij}} \exp\left\{ -\left[(y_{ij} - s_{k_{ij} - 1}) \lambda_{ik_{ij}} + I(k_{ij} > 1)\sum_{g=1}^{k_{ij} - 1} (s_g - s_{g-1}) \lambda_{ig}\right]\right\} \\
&=& \exp\left\{ \sum_{j=1}^{n_i} \delta_{ij} \log(\lambda_{ik_j}) - \left[\sum_{j=1}^{n_i} (y_{ij} - s_{k_{ij}-1})\lambda_{k_{ij}} + \sum_{k=1}^K (s_k - s_{k-1}) \lambda_{ik} m_{ik} \right] \right\} \\
&=& \exp\left\{ \sum_{k=1}^K \left[ n_{ik} \log(\lambda_{ik}) -  \lambda_{ik}\left( \sum_{j: k_{ij} = k} (y_{ij} - s_{k - 1}) + m_{ik}(s_k - s_{k-1}) \right)\right] \right\}
\end{eqnarray*}
where $y_{ij}$ is the $j$th response from the $i$th terminal node, $k_{ij}$ denotes the interval $y_{ij}$ belongs to, $m_{ik}$ is the number of observations of $\by_i$ which are greater than $s_k$, $n_{ik}$ is the number of uncensored observations of $\by_i$ which fall in interval $k$, and $\delta_{ij}$ is the censoring indicator for $y_{ij}$.

We model the hazard function $h_i(t) = \lambda_{ik} =  \exp(\omega_i(z_k))$ for $t \in (s_{k-1},s_k]$ where $\omega_i(t)$ is a zero-mean Gaussian process with covariance function $C_i(t,t^\prime) = \tau_i^2 \exp\{-|t - t^\prime|/l_i \}$.  This covariance kernel is chosen primarily for its computational advantages; see the appendix.  Given a set of $K$ locations, $\Bf_i = (\omega_i(z_1),\ldots,\omega_i(z_K))^T$  follows a multivariate-normal distribution with mean 0 and covariance matrix $\Sigma_i$ with the $i^\prime,j^\prime$th element of $\Sigma_i$ being $C_i(z_{i^\prime},z_{j^\prime})$. The posterior distribution of $\Bf_i$, assuming the covariance function parameters are known, is
\[
 p(\Bf_i \mid \by_i, \bdelta_i) \propto p(\by_i, \bdelta_i \mid \Bf_i) \pi(\Bf_i) = (2\pi)^{-K/2}|\Sigma_i|^{-.5} \exp\left\{g_0(\Bf_i) \right\}
\]
where
\[
g_0(\Bf_i) = \sum_{k=1}^K \left[ n_{ik} w_{ik} -  \exp(w_{ik})\left( \sum_{j: k_{ij} = k} (y_{ij} - s_{k - 1}) + m_k(s_k - s_{k-1}) \right)\right] -\frac{1}{2} \Bf_i ^T \Sigma^{-1}_i  \Bf_i
\]
and $w_{ik}$ is the $k$th element of $\Bf_i$.  Note that $g_0(\Bf_i)$ implicitly depends on the data, bins, and covariance function parameters.  Due to the assumption of independence between partitions, the full joint posterior distribution of a tree, $T$, with $M$ terminal nodes and $\Bf_1,\ldots,\Bf_M$ can be expressed as
\[
p(T,\Bf_1,\ldots,\Bf_M \mid \by, \bdelta) \propto \pi(T) \prod_{i=1}^M p(\by_i, \bdelta_i \mid \Bf_i) \pi(\Bf_i)
\]
where $\pi(T)$ is the prior on the tree as in \eqref{eq:prior}.

The goal in many partition models is to find a suitable partition structure.  In a Bayesian framework, this corresponds to searching the posterior of the tree, i.e.,
\begin{eqnarray}
p(T \mid \by) &=& \int \ldots \int p(T,\Bf_1,\ldots,\Bf_M \mid \by) d\Bf_1 \ldots d\Bf_M \nonumber \\
&\propto & \pi(T) \prod_{i=1}^M \int p(\by_i \mid \Bf_i) \pi(\Bf_i) d\Bf_i. \label{eq:marginal}
\end{eqnarray}
To do this, often a model is chosen in which all node-specific parameters can be integrated out analytically.  When the integrals have an analytic form, often an efficient reversible jump MCMC chain \citep{green1995reversible} can be constructed to search the posterior partition space. Unfortunately there is no analytical solution for the integrals in \eqref{eq:marginal} for the present model.  To overcome this issue, a Laplace approximation is used to approximate each integral. Laplace approximations have been used previously in various survival model settings. \citet{levine2014bayesian} utilize a Laplace approximation to search trees in a Bayesian frailty model and \citet{martino2011approximate} use integrated nested Laplace approximations via data augmentation to fit a Cox proportional hazards model to avoid MCMC. In our method, the Laplace approximation is similar to \citet{levine2014bayesian}, but we do not have random effects and we also optimize the covariance function hyperparameters via empirical Bayes.  \citet{payne2020conditional} utilize a Laplace approximation in a partition model framework for non-survival conditional density estimation using a Voronoi partition.

In usual form, the Laplace approximation approximates the integral with a normal distribution derived from a Taylor expansion of $\log[p(\by_i \mid \Bf_i) p(\Bf_i)]$.  The Laplace approximation has a closed form given by
\begin{equation} \label{eq:laplace}
\int p(\by_i, \bdelta_i \mid \Bf_i) \pi(\Bf_i) d\Bf_i \approx |D_i + \Sigma_i^{-1}|^{\frac{1}{2}}|\Sigma_i|^{-\frac{1}{2}} \exp\left\{g_0(\hat{\Bf}_i)\right\}
\end{equation}
where $D_i$ is a diagonal matrix with $j^\prime$th element $\exp(\hat{w}_{ij^\prime})\left[\sum_{j : k_{ij} = j^\prime} (y_{ij} - s_{j^\prime-1}) + (s_{j^\prime} - s_{j^\prime-1})m_{j^\prime}\right]$ and $\hat{\Bf}_i = \arg\max_{\Bf_i} p(\by_i \mid \Bf_i) \pi(\Bf_i)$, which can be obtained using Newton's method.

Thus far we have assumed that the covariance function parameters are fixed and known.  In practice, we employ empirical Bayes to set these values in each terminal node.  Specifically, we choose $\tau_i$ and $l_i$ to maximize $p(\by_i \mid \Bf_i) \pi(\Bf_i \mid l_i, \tau_i) \pi(l_i) \pi(\tau_i)$ where $\pi(l_i)$, $\pi(\tau_i)$ are t-densities with one degree of freedom and $\pi(\tau_i) = t(\tau_i / 10; 1)$ where $t(\cdot; 1)$ is the density of a t-distribution with one degree of freedom.

At this point, we have a way to approximate $p(T \mid \by)$, up to a normalizing constant, for a single tree utilizing Laplace approximations.  We now need to determine a way to search the approximate posterior.  As is common in many tree models such as \citet{chipman1998bayesian} and \citet{gramacy}, we employ Markov chain Monte Carlo methods.  Due to the varying number of parameters in the prior of the tree structure, a reversible jump Markov chain Monte Carlo algorithm \citep{green1995reversible} will be outlined in Section 4.

\section{Consistency results}
In this section, we show the density estimation consistency  when the true data generating density follows a partition model induced by a tree. If the survival function for the true data generating density has a tree partition form, then first we show that as $n$, the number of observations goes to infinity, the posterior distribution concentrates on  any weak neighborhood of the true data generating density with probability one. Then under slightly stronger  conditions, we establish the $L^1$ consistency.

Let  $P_1^*,\dots,P_{m^*}^*$ be the true partition of the covariate space denoted as  $\mathscr{P}^*=\{P^*_1,\dots,P_{m^*}^*\}$, and  let $h_{0,i}(t)>0,\ t>0$ for $i=1,\dots,m^*$,  be the true hazard function for partition $P_i^*$ and the corresponding true uncensored data generating  density for the $i$th partition is given by  $p_{0,i}(t)=h_{0,i}(t)e^{-H_{0,i}(t)}$ where $H_{0,i}(t)=\int_0^th_{0,i}(s)ds$. As mentioned earlier, $p_c$ is the density of the censoring distribution. Let $F_c$ and and $F_{0,i}$ be corresponding cdfs and $1-F_{0,i}(t)=e^{-H_{0,i}(t)}$. For a generic partition $\mathscr{P}=\{P_1,\dots, P_{k'}\}$ of the covariate space and the partition densities $p_{i}(\cdot)$,  for observation $y,\ {\bf x}$, and the censoring indicator function $\delta$, the data generating density $p(\cdot)$  is given by
\begin{eqnarray}
p(y,\delta \mid {\bf x},\mathscr{P})=\sum_{i=1}^{k'}\left[p_{i}(y)\int_y^\infty p_c(s)ds{\bf I}_{\delta=1}+p_c(y)\int_y^\infty p_{i}(s)ds{\bf I}_{\delta=0}\right]{\bf I}_{{\bf x}\in P_i}.
\label{likelihood1}
\end{eqnarray}

We show that under a Gaussian process tree partition model, for any bounded continuous function $g_1(\cdot)$, the posterior probability of the weak neighborhood around true data generating density $p^*$, $U^{\epsilon}_{g_1}=\{p:|\int g_1p-\int g_1p^*|<\epsilon\} $, where $\epsilon>0$,  approaches one almost surely, as $n$ approaches infinity.

We observe  $(y_1,\delta_1),\ldots,(y_n,\delta_n)$ at observed covariates ${\bf x}_1,\ldots,{\bf x}_n$, where ${\bf x}_j$ are i.i.d. for $j=1,\dots,n$. Let  $\#(T)$ denote the number of partitions/terminal nodes for a tree $T$. Denoting the generic prior distribution by $\pi(\cdot)$, for $T_n=\max_i\{y_i\}$,  from Section 2 we have the following prior:
\begin{eqnarray}
P_1,\dots,P_{\#(T)}&\sim&\pi(T),\nonumber\\
\log h_i(t)=\beta_i+w_i(t) &\sim & GP(\beta_i,\mathscr{K}), \nonumber \\
\mathscr{K}(s,t)&=&\tau_i^2 K\left(\frac{s}{l_i},\frac{t}{l_i}\right),\ t\in [0,T_n];\ i=1,\cdots,\#(T),\nonumber\\
\pi(\tau^2_i,l_i)&=&\pi(\tau^2_i)\pi(l_i),\nonumber\\
\beta_i&\sim&\pi(\beta_i).
\label{prior}
\end{eqnarray}

We assume the following:
 \begin{itemize}
 \item[A1] $\log (h_{0,i}(t))=\int_0^\infty \tilde{g}_i(s)K(s/l_i^*,t/l_i^*)ds+\beta_{0,i}$, for some bounded, Lipschitz-continuous and integrable functions $\tilde{g}_i$'s, with $ w_{0.i}(t)=\int_0^\infty \tilde{g}_i(s)K(s/l_i^*,t/l_i^*)ds.$
 \item[A2] $E_{p_{0,i}}[t],\ E_{p_c}[t]<\infty$.
 \item[A3] $\max_t|\int_0^{T'}\tilde{g}_i(s)K(s/l_i^*,t/l_i^*)ds-w_{0,i}(t)|=o(1/T')$, as $T'\uparrow \infty$,  where $w_{0,i}(t)$'s are defined in [A1].
 \item[A4]  The covariates are assumed to be bounded. Each continuous covariate is partitioned into $m_{n}$ equi-spaced  grids for the tree splitting rule, such that $m_n\rightarrow \infty$ as $n\uparrow \infty$ and   $-\log \pi(T)= o(n)$ for any  tree $T$ such that $\#(T)<\infty$.
 \item[A5] $t^\alpha\preceq max\{-\log (1-F_c(t)),-\log (1-F_{0,i}(t))\}; t>0$, for some $\alpha>0$ and for all $i=1,\ldots,m^*$. Also, $p_c$ and $p_{0,i}$'s have the same support.  Here $\preceq$ denotes less than equal to a constant multiple.
 \item [A6]Priors $\pi(\tau^2_i)$ and  $\pi(l_i)$ have strictly positive continuous densities on the positive real line, and $\pi(\beta_i)$ has a strictly positive  continuous density on the real line.
 \end{itemize}

 Conditions [A1] and [A3] ensure that a constant shift of the  true log-hazard function remains in the Gaussian process support. Conditions [A1], [A4], [A5], and [A6] guarantee enough prior probability in a small neighborhood around  the true data generating density under the partition-tree model. Condition [A2] is needed for prior Kullback-Leibler support.

 Suppose the  true partition $\mathscr{P}^*$  is achieved for  a tree $T^*$, and let $\hat{\mathscr{P}}^*=\hat{\mathscr{P}}_T^*=\{\hat{P}_l^*:l=1,\cdots,m^*\}$ corresponding to $T=\hat{T}^*$ be an approximating partition, with $\hat{T}^*$ having the same number of terminal nodes as in $T^*$. We define the following distance between these two partitions  by
\[d(\mathscr{P}^*,\hat{\mathscr{P}}^*)=\text{min}_{c(1),\dots,c(m^*)} \sum_{i=1}^{m^*}\mu_x(P_i^*\Delta\hat{P}_{c(i)}^*)  \] where $c(1),\ldots,c(m^*)$ is a permutation of $1,\ldots,m^*$, and  $\mu_x$ is a measure associated with the covariate space which can be an absolutely continuous Lebesgue measure, a counting measure, or a product of both. Let $\mathscr{P}^*\cap \hat{\mathscr{P}}^*=\cup_{i=1}^{m^*}P_i^*\cap \hat{P}^*_{\hat{c}(i)}$ and  $\mathscr{P}^*\Delta \hat{\mathscr{P}}^*=\cup_{i=1}^{m^*}P_i^*\Delta \hat{P}^*_{\hat{c}(i)}$, where $\hat{c}(\cdot)$ is the minimizing permutation in the definition of $d(\mathscr{P}^*,\hat{\mathscr{P}}^*)$. For notational simplicity,  later $P_i^*\cap \hat{P}^*_{\hat{c}(i)}$ will be denoted by $P_i^*\cap \hat{P}^*_i$, that is $\hat{P}^*_i$ approximates $P^*_i$.

Now, we can state the following results. The proofs for the results are relegated to the appendix.

 \begin{proposition}
Under [A5], $T_n=o(n^\beta)$ almost surely,  for any $\beta>0$, where $T_n=\max_i\{y_i\}$.
\label{max_sup}
\end{proposition}
 \begin{proposition}
For partition $\mathscr{P}^*$, given $\xi>0$ there exists a set $S_T^{\xi}=\{T: d( \hat{\mathscr{P}}_T^*,\mathscr{P}^*)<\xi , \#(T)=\#(T^*)\}$ such that  under [A4], $-\log \pi(S_T^{\xi})=o(n)$.
\label{tree_support1}
\end{proposition}
Next, using Proposition \ref{max_sup} and the results from \citet{van2008reproducing},  \citet{aurzada2008small}  on Gaussian process small ball probability, we show the following.
\begin{lemma}

Let $\tilde{{w}}_{0,i}^n(t)=\int_0^{T_n}\tilde{g}_i(s)K(s/l_i^*,t/l_i^*)ds$ for $0\leq t\leq T_n$.  For any $\xi>0$, $S^i_\xi=\{w_i(\cdot):|w_i(t)-\tilde{w}^n_{0,i}(t)|<\xi/T_n, 0\leq t\leq T_n\}$. Then under  the Gaussian process prior on $w_i(\cdot)$ from \eqref{prior}, under A5, A6, $-\log\pi(S^i_\xi)=o(n)$.
\label{gp_support}
\end{lemma}

 Proposition \ref{tree_support1} and Lemma \ref{gp_support} establish a large enough prior probability around the truth, and Proposition \ref{max_sup} bounds $T_n$ and is used in Lemma \ref{gp_support}.  Now, we state the result determining the Kullback-Leibler neighborhood around the truth, which combined with Propositions \ref{max_sup}, \ref{tree_support1}, and Lemma \ref{gp_support} will guarantee that we have sufficient prior probability around a small Kullback-Leibler neighborhood around the truth. The support for the $p_{0,i},\ p_c$ are assumed to be on the positive real line, under which $T_n$ increases to infinity almost surely,  as $n$ increases to infinity. If the supports are the same for $p_{0,i},\ p_c$ and bounded then  the integral in assumption [A1] is replaced by the integral over the support of $p_{0,i}$'s and the  following Theorem \ref{thm1} will hold.

 \begin{lemma}
Let $p=\hat{p}$ be an estimator of $p^*(t,\delta,{\bf x})$ with hazard function $h_i(t)=e^{w_i(t)+\beta_i}$  for $0\leq t\leq T_n$ in partition $\hat{P}_i$ and $h_i(t)=h_i(T_n)$ for $t>T_n$, where $\hat{\mathscr{P}}^*=\{\hat{P_i}^*\}_{i=1,\dots,m^*}$ a partition of the covariate space,  and $d(\mathscr{P},\hat{\mathscr{P}^*})<\epsilon$.
Then, given $\epsilon>0$, there exists $\xi>0$ such that  for  $|\beta_i-\beta_{0,i}|<\xi$ and  $|w_i(t)-w_{0,i}(t)|<\xi$, for $0\leq t\leq T_n$, in $\hat{\mathscr{P}^*}\cap {\mathscr{P}^*}$, we have under [A1]--[A3] ,  $KL(p^*,p)=\int p^* \log\frac{p^*}{p}<k\epsilon$ for large $n$, for some universal constant $k>0$ .
\label{kl_sup}
\end{lemma}

 Finally, we can state the following consistency theorem.

 \begin{theorem}

For any bounded continuous function $g$, and for $\epsilon>0$, let $U^{\epsilon}_g=\{p:|\int g p^*(\cdot)-\int g p(\cdot)|<\epsilon\}$, be a weak neighborhood of the true joint density $p^*(\cdot)$. Under the model in equation \ref{likelihood1}, prior in equation \ref{prior}, under assumptions [A1]--[A6], and under the true data generating distribution $p^*$, $\Pi_n(U^{\epsilon}_g\mid \cdot)\rightarrow 1$ almost surely as $n\rightarrow \infty$, where $\Pi_n(\cdot)$ denotes the posterior distribution based on $n$ observations.
\label{thm1}
\end{theorem}

Theorem \ref{thm1} establishes consistency in  weak topology with $\Pi_n(\mid \cdot)$ being weakly consistent at $\{ {p_{0,i}}_{i=1,\cdots,m^*},p_c,\mathscr{P}^* \}$. Even though  we state  Theorem \ref{thm1} under the support condition given in [A3], the consistency result will hold true for more  general support  under additional mild conditions, given in the following.

\begin{remark}
Theorem \ref{thm1} holds for piecewise constant hazard functions which eventually converges to a constant that is  $\log h_{0,i}(t)=c_{0,i}$ for $t\geq T_{0,i}\geq 0$ for $i=1,\ldots,m^*$, if we assume $c_0l_i^{\alpha_{0,i}}\leq\pi(l_i)$ for some $\alpha_{0,i}>0,c_0>0$ for $|l_i|\leq c_l, c_l>0$.
\label{rmrk1}

\end{remark}

\begin{remark}
Theorem \ref{thm1} holds for bounded Lipschitz continuous $\log h_{0,i}(t)$'s  if we assume $c_0l_i^{\alpha_{0,i}}\leq\pi(l_i)$ for some $\alpha_{0,i}>0,c_0>0$ for $|l_i|\leq c_l, c_l>0$.
\label{rmrk2}
\end{remark}
We can also slightly relax the tree support condition as follows.
\begin{remark}

If the true partition is not induced by a tree, but given any $\epsilon>0$, can be approximated by a partition ${\mathscr{P}}_1^*$, induced by a  tree ${T}_1^*$ (may depend upon $\epsilon$), such that  $d(\mathscr{P}^*,{\mathscr{P}}_1^*)<\epsilon$, then Theorem \ref{thm1} holds.
\label{rmrk3}

\end{remark}

Let $D_n$ be the maximum number of terminal nodes allowed for $n$ observations.  Note that  $D_n$ will induce a total  less than  $D_n{(mm_n)}^{D_n}$ many possible tree formations where $m_n$ is the number of grids for a covariate and we have $m$ covariates and they are continuous. Establishing $L^1$ or strong consistency will require controlling or bounding the size of this model space. Let $\{M_n\}_{n\geq 1}$ be a  positive sequence increasing to infinity.
For establishing strong consistency, we will need to assume the following conditions.
For  positive sequences $\{a_n\}_{n\geq 1},\{b_n\}_{n\geq 1}$,  $a_n\prec b_n$, implies $a_n/b_n\rightarrow 0$, as $n\uparrow \infty$.

\begin{itemize}

\item[B1]  The covariates are assumed to be bounded. Each continuous covariate is partitioned into $m_{n}$ equi-spaced  grids for the tree splitting rule, such that $m_n\rightarrow \infty$ as $n\uparrow \infty$, and $\log m_n\sim O(\log n)$, and   $-\log \pi(T)= o(M_n^2)$ for any  tree $T$ such that $\#(T)<\infty$.
\item[B2]  $D^{1+\delta'}_nM_n^2  \log m_n=o(n)$, for some $\delta'>0$.
\item[B3] $M_n^2\prec n$ and $n/M_n^2= o(n^{\delta})$ for any $\delta>0$.

 \item[B4]  $\log (h_{0,i}(t))=\int_0^\infty \tilde{g}_i(s)K(s/l_i^*,t/l_i^*)ds+\beta_{0,i}$, for some bounded, Lipschitz-continuous, functions $\tilde{g}_i$'s, that are supported on a compact set. Let $ w_{0.i}(t)=\int_0^\infty \tilde{g}_i(s)K(s/l_i^*,t/l_i^*)ds.$

 \item[B5] $\pi(l_{i'})$ is supported on equi-spaced $M_n^{d'}$, $d'>0$, many grids in  a closed interval  not containing zero, that is a subset of the positive real line,   and that contains  $l_i^*$ as an interior point, for any $i',i$, and  $-\log \pi(l)=o(M_n^2)$, for $l$ in the support.

 Also,
 $\pi(\tau^2_i)$  and $\pi(\beta_i)$ have strictly positive  continuous densities on their supports, positive real line and real line, respectively, and $min\{-\log P(\tau_i^2>t),-\log P(|\beta_i|>t)\}\succeq t^{\alpha'}$, $\alpha'>0,t>0$.

\end{itemize}

Condition $B2$ bounds the $L^1$ covering number of the prior space. Conditions $B3,B4$ are needed for   controlling  the ratio of the joint density  for the parameter value around the truth and the true joint density, almost surely, for a small neighborhood around the truth and to guarantee exponential decaying prior probability outside a compact set. Condition  $B1$ is needed to approximate the true data generating tree, and  establish sufficient prior probability in that small neighborhood.

We can now state the following strong consistency result. For the strong consistency we consider estimated densities that are supported on $[0,T_n]$ and show that they will remain in a small $L^1$ neighborhood of the true density with probability one.

\begin{theorem}

Let $\epsilon>0$ and $U^{\epsilon}=\{p:\int | p^*(\cdot) -p(\cdot)|<\epsilon\}$, be a $L^1$ neighborhood of the true joint density $p^*(\cdot)$. Under the model in equation \ref{likelihood1}, prior in equation \ref{prior}, under assumptions [A2],[A4]--[A5],[B1]--[B5] and under the true data generating distribution $p^*$, $\Pi_n(U^{\epsilon}\mid \cdot)\rightarrow 1$ almost surely as $n\rightarrow \infty$, where $\Pi_n(\cdot)$ denotes the posterior distribution based on $n$ observations.
\label{thm2}
\end{theorem}

Technical details and related proofs of the results are given in  Appendix ii.
\section{MCMC}
Since the dimension of the tree parameter space is variable, we use a reversible jump MCMC algorithm to search the posterior parameter space  \citep{green1995reversible}. The algorithm is inspired by those of \citet{chipman1998bayesian}  and \citet{gramacy} with a few modifications. The primary goal in constructing the algorithm is to obtain draws from $\tilde{p}(T \mid \by)$, the approximate posterior obtained using the Laplace approximations.  This is accomplished in the reversible jump MCMC algorithm by modifying an existing tree with one of four possible moves: grow, prune, change, and swap.

\begin{itemize}
\item {\bf Grow:} Randomly select a terminal node and add two child nodes based on a split rule.  The split rule is selected by first randomly selecting a variable and then randomly selecting an available rule of that variable (i.e.\ the rule is drawn from the prior).  This increases the number of terminal nodes of the tree by one, thereby ``growing" the tree.

\item {\bf Prune:} The opposite of a grow step.  Randomly select a node whose children are both terminal nodes and delete the split rule and the children nodes.  This reduces the number of terminal nodes of the tree by one, thereby ``pruning" the tree.

\item {\bf Change:} Randomly select an interior node and change the existing split rule.  If the current rule is based on a categorical variable, a new rule is drawn randomly from the prior by first randomly selecting a variable, and then randomly selecting a rule from that variable.  If the current rule is based on a numeric variable, there are two possible changes that can be made with pre-specified probability:  1) Move the current rule to the next largest or smallest rule available for the current variable, 2) Draw a new rule from the prior.

\item {\bf Swap:} Randomly choose a parent and child pair which are both internal nodes and swap their rules.  If, however, the parent and child's rules are based on the same continuous variable, a rotate step is employed \citep{gramacy}.  If both children of the parent node have the same rule, both children swap their rules with the parent node.
\end{itemize}

The MCMC algorithm proceeds as follows:

\begin{enumerate}
	\item Initialize a starting tree, $T$.
	\item Propose a new tree, $T^1$, from the current tree using a grow, prune, change, or swap step (as appropriate based on current tree).
	\item Accept $T^1$ with probability
	\begin{eqnarray} \label{eq:acceptance}
	\alpha = \min\left\{1,\frac{q(T \mid T^1) \tilde{p}(T^1 \mid \by)}{q(T^1 \mid T) \tilde{p}(T \mid \by)}  \right\}
	\end{eqnarray}

	where $q(\cdot \mid \cdot)$ is the proposal distribution and $\tilde{p}(T \mid \by)$ is the approximate posterior of $T$ obtained from the Laplace approximation; otherwise keep the current tree, $T$.
	\item Repeat steps 2-3 until the stationary distribution is reached and the desired number of samples is obtained.
\end{enumerate}

In practice we replace the approximate posterior of the tree in \eqref{eq:acceptance} with the likelihood-prior product in \eqref{eq:marginal} since the normalizing constants cancel.

From their inception, Bayesian tree models have been notorious for their difficulty in mixing when using MCMC procedures.   Various methods have been proposed to overcome this issue including restarting the MCMC chain multiple times \citep{chipman1998bayesian}, parallel tempering \citep{gramacy2010categorical}, and multiset sampling \citep{leman2009}.  We implement parallel tempering \citep{geyer1991markov} modeled after \citet{gramacy2010categorical}.  Additional technical details regarding the MCMC, are provided in the appendix.

\section{Simulations \& Applications}
\subsection{Preliminaries}
For the following analyses, we use the same set of tuning parameters for the MCMC chain.  We set the probability of grow, prune, change, and swap steps to be .25.  For the change step for a continuous variable, we choose to move the rule to an adjacent rule with probability .75, otherwise a rule is drawn from the prior.  We run each analysis with 8 parallel chains with the minimum inverse-temperature of .1, and propose a trade between the parallel chains every 10 MCMC iterations. Trees are required to have a minimum of 25 observations at a terminal node.
\subsection{Tree partition with non-proportional hazards} \label{sec:tp}
Our first simulation demonstrates the flexibility and interpretability of the THM.  One thousand data points were simulated from the following model:
\begin{equation} \nonumber
\begin{aligned}
X_1  & \sim  \text{Uniform}(0,10), \\
X_2 & \sim  \text{Multinomial}(\text{A,B,C,D};\ .25,.25,.25,.25), \\
Y & \sim  \begin{cases}
\text{Weibull}(5,2) & \text{ if } X_1 > 5 \text{ and } X_2 \in \{\text{A,B} \} \\
\text{Weibull}(1,5) & \text{ if } X_1 \leq 5 \text{ and } X_2 \in \{\text{A} \} \\
\text{Weibull}(.5,.9) & \text{ if } X_1 \leq 5 \text{ and } X_2 \in \{\text{B} \} \\
\text{Weibull}(5,5) & \text{ if } X_1 \leq 3 \text{ and } X_2 \in \{\text{C,D} \} \\
\text{Weibull}(.5,.5) & \text{ if } 3 < X_1 \leq 7 \text{ and } X_2 \in \{\text{C,D} \} \\
G & \text{ if } X_1 > 7 \text{ and } X_2 \in \{\text{C,D} \} \\
\end{cases}
\end{aligned}
\end{equation}
where $G$ is a random variable drawn from an arbitrary piecewise linear survival function.  The data were also right censored with a 5\% censoring rate. The Markov chain Monte Carlo algorithm ran for $10^4$ iterations.  The tree with the highest posterior probability from the chain is shown in Figure \ref{fig:sim9tree} along with the posterior and true survival curves.  Note that since we are using observed data values to choose the rules at each node, the rules for numeric variables will not be exactly the same as those used to generate the tree.  This example highlights the strength of the THM, simultaneously producing easy interpretation and flexibility.  Interpretation is simple since only one tree is needed, but the the various shapes of the survival curves can still be modeled flexibly at each terminal node.

\begin{figure}
\begin{center}
  \subfigure[]{\includegraphics[trim=40 30 40 10, clip,width = .45\textwidth]{figs/sim9tree.eps}}
  \subfigure[]{\includegraphics[trim=0 10 0 10, clip,width = .45\textwidth]{figs/sim9survival.eps}}
  \caption{(a) The tree with the highest posterior probability from the model in Section \ref{sec:tp}.  (b) The posterior mean (dotted), 95\% credible intervals (dashed), and the true survival curve (solid) in each partition element.}
\label{fig:sim9tree}
\end{center}
\end{figure}

\subsection{Cox proportional hazards} \label{sec:lph}
Although the present partition model assumes the data are generated from a tree structure, it still performs well when this assumption does not hold, e.g., when dealing with the classic linear proportional hazards model \citep{cox1972regression}.  One thousand observations were simulated from an exponential proportional hazards model with an 8\% censoring rate.  The model is
\begin{equation}
\begin{aligned}
h(t \mid \mathbf{X} = \mathbf{x}) =\ & \frac{1}{2} \exp\left\{\mathbf{x} \mathbf{\beta} \right\} \\
\mathbf{X} =\ & (X_1,\ldots,X_6) \\
X_1,\ldots,X_5  \sim\ & \text{Uniform(0,1)} \\
X_6  \sim\ & \text{Bernoulli(.5)} \\
\beta =\ & (-1, 1, 2, 0, 0, -2)^T.
\end{aligned}
\end{equation}
where $h(\cdot \mid \cdot)$ represents the hazard function.

Note that we have two spurious covariates ($X_4,X_5$) which do not influence the hazard function.  The MCMC algorithm ran for $10^4$ iterations.  In the tree with the highest posterior probability (8 terminal nodes), neither of the two spurious covariates were included.  To give a measure of performance, Figure \ref{fig:sim7survival} plots the posterior survival curve and the exponential survival curve using the median true hazard value of the data that reach the terminal node.  This comparison serves as an indicator of where the survival distribution should lie in that region of the covariate space, and indicates the model captures the important features of the curves well.

\begin{figure}[t]
\begin{center}
\includegraphics[scale=.55]{figs/sim7survival.eps}
\caption{The posterior mean (dotted) and 95\% credible intervals (dashed).  The survival curve corresponding to the median true hazard function of the data in each terminal node is also plotted (solid).}
\label{fig:sim7survival}
\end{center}
\end{figure}

\subsection{Prognostic \& predictive biomarker, comparison with BART, survival random forests, and Cox}
Identifying predictive and prognostic biomarker effects is important in drug development.  Sometimes there is a hypothesis that a certain sub-population, based on a predictive biomarker, will respond better to a treatment.   Other times the value of a prognostic biomarker is related to a patient's prognosis.  It is possible a biomarker may be both prognostic and predictive.

A single biomarker, $b$, and a treatment assignment indicator, $a$, are simulated for 1,000 subjects.  The biomarker is both prognostic and predictive, and there are two treatment arms. The time to event, $y$, is simulated from a Weibull model with random right-censoring in 7.8\% of subjects,
\begin{equation}
\begin{aligned}
    a &\sim Binomial(1, .5), \\
    b &\sim Uniform(0, 1), \\
    y &\sim Weibull(1 + 2ab, 1 + 5b),
\end{aligned}
\end{equation}
where $Weibull(\lambda, k)$ is a Weibull distribution with scale $\lambda$ and shape $k$.  For small biomarker values $b$ nearer to 0, there is a small treatment effect, whereas for values of $b$ closer to 1, there is a large treatment effect creating a predictive biomarker.  For a fixed treatment assignment $a$, the shape of the survival curve is a function of $b$, creating a prognostic effect for the biomarker.

The tree partition model was fit to the simulated data, and the tree with the highest posterior probability is plotted in Figure \ref{fig:pred_prog}. The tree and the estimated survival curves at each terminal node capture the predictive and prognostic effects.  The split rule of the root node indicates that there is a difference in survival functions between the two treatment arms.  Descending either right or left from the root node, all remaining splits are on the biomarker, indicating a prognostic biomarker effect.  For two patients with similar biomarker values, but different treatment assignments, their survival curves can be quite different.  For instance, a patient with a biomarker value of .8 has a much better survival curve with treatment assignment 1 verses 0; this indicates that the biomarker is also predictive.

The simulated data were fit to three additional models: a frequentist Cox proportional-hazard model with a non-parametric baseline hazard an an interaction between $a$ and $b$; Bayesian additive regression trees (BART) \citep{sparapani2016nonparametric}; and survival random forests \citep{ishwaran2008random}.  The survival curves, the 95\% confidence/credible intervals, and the true treatment effect for a patient with $a = 1$ and $b = .10$ are plotted in Figure \ref{fig:pred_prog_comparison}; for the survival random forest method, the 95\% intervals are calculated from the predictions of all trees.  The changing shape of the survival curve is not captured well in the Cox model since the proportional hazard assumption is too restrictive.  This is problematic; even if the Cox model identifies a significant treatment by biomarker interaction, which is often used to determine a predictive biomarker effect, the estimates of the survival curve may be inaccurate.  On the other hand, BART and random forests estimate the survival curves well, but often custom post-processing of the results in a potentially large-volume covariate space are required to better understand which predictors are prognostic, predictive, or both.  The THM automatically provides both clear interpretation and flexible estimation of the survival curves.
\begin{figure}[t]
    \centering
    \includegraphics[scale = .55]{figs/pred_prog.eps}
    \caption{Tree with the highest posterior probability with the posterior mean (dotted) and 95\% credible intervals (dashed) at each terminal node.}
    \label{fig:pred_prog}
\end{figure}
\begin{figure}[t]
    \centering
    \includegraphics[scale = .55]{figs/pred_prog_comparison2.eps}
    \caption{The true survival curve (solid) for a patient with $a = 1$ and $b = .10$, estimated mean (dotted) and 95\% confidence/credible intervals (dashed) for the frequentist Cox, BART, and THM.}
    \label{fig:pred_prog_comparison}
\end{figure}
\subsection{Analysis of Primary Biliary Cirrhosis Data}

Primary biliary cirrhosis survival data collected from 1974-1984 by the Mayo Clinic and obtained via the R survival package \citep{survival-package} were fit using the THM.  Covariates included treatment group (D-penicillamine or placebo), age, sex, edema status (0 no edema, 0.5 untreated or successfully treated, 1 edema despite diuretic therapy), serum bilirunbin (mg/dl), serum albumin (g/dl), and standardised blood clotting time.  Only subjects who were randomized to treatment and had no missing values were included in the analysis for a total of 276 subjects.  The outcome of interest was time to death, and patients were censored at the end of study or if they received a transplant, yielding a censoring rate of 60\%.

Ten-thousand MCMC samples were drawn and the tree with the highest posterior probability is shown in Figure \ref{fig:pbc} with the survival curves plotted at each terminal node.  The black curves represent the posterior mean and 95\% credible intervals from the tree model, and the red lines provide mean and 95\% confidence intervals from the Kaplan-Meier curves.  The plot indicates that the THM is very flexible, and close to the Kaplan-Meier estimates.

The tree indicates that patients with bilirubin levels less than 1.9 have the highest survival rates.  For subjects with no edema and bilirubin levels greater than 1.9 and less than or equal to 5.7, patients under the age of 45 had better survival.  Patients with bilirubin levels greater than 1.9 and edema had the worst survival.  These conclusions are consistent with the findings of \citet{dickson1989prognosis}, who used a Cox proportional hazards model on primary biliary cirrhosis data, however, the tree method provides greater flexibility in estimating the survival curves, especially in this case where the hazards are not proportional.  Last, we note that the tree did not split on treatment status, which is consistent with the findings of \citet{gong2004d} who performed a systematic review of the effect of D-penicillamine on patients with primary biliary cirrhosis and concluded that it does not decrease the risk of mortality.

\begin{figure}
    \centering
    \includegraphics[scale=.55]{figs/pbc_km.eps}
    \caption{The tree with the highest posterior probability and the posterior mean (black, dotted), 95\% credible intervals (black, dashed), and the Kaplan-Meier estimates (red) at each terminal node.}
    \label{fig:pbc}
\end{figure}

\section{Conclusion}
The main advantage the tree hazards model (THM) is its ease in interpretation coupled with its flexibility.  Inference is easily obtained from inspecting posterior trees and the flexibly estimated survival curves at the terminal nodes.  In both simulations and applications, the model performs well in identifying interpretable partitions and flexibly estimating survival curves.  Furthermore, we have established weak consistency for the partition model under the proposed kernel, where we can accommodate flexible forms of the hazard function under moderate conditions.

\bibliographystyle{apalike}

\section{Introduction}
Survival analysis for time-to-event data is used in a variety of disciplines, particularly in medical research. The Cox proportional hazard, accelerated failure time, and proportional odds models and their extensions are popular modeling tools to relate survival times with predictors or covariates in the presence of censorship.
Bayesian analysis of survival data is popular due to its flexibility to include prior information, suitability to develop complex models, and ability to construct efficient algorithms for full scale uncertainty quantification \citep{ibrahim2014b,sinha1997semiparametric}. Bayesian semi-parametric regression models for survival data have been developed using the Gamma process, the Dirichlet process, and Polya tree priors and their mixtures \citep{kalbfleisch1978non,christensen1988modelling,walker1999bayesian,hanson2002modeling,kuo1997bayesian,gelfand1995bayesian,gelfand2003bayesian}. Though these are popular models, they are constrained by parametric assumptions such as proportional hazards, constant odds, linearity, etc. Recently, several flexible Bayesian models have been developed  to accommodate data complexities using dependent processes and its extensions \citep{de2009bayesian,jara2011class,nieto2014bayesian,nipoti2018bayesian,riva2021survival}.

An alternative approach is to utilize survival trees and forest-based models which are popular non-parametric alternatives to semi-parametric models and have the ability to identify non-linear interactions among variables of interest. Another benefit of tree models is their easy interpretation. Moreover, they can cluster subjects automatically according to their survival behavior based on covariates to identify prognostic groups.  Tree-based survival analysis was introduced in the 1980s \citep{leblanc1993survival,bou2011review} primarily through the extension of classification and regression tree (CART) models \citep{breiman1984classification} to censored survival data \citep{ciampi1981approach,gordon1985tree}. The major difference among these papers is the use of different splitting criteria such as log-rank statistics, likelihood ratio statistics, Wilcoxon-Gehan statistics, node deviance measures, and different impurity measures. A crucial issue of these tree models is the selection of a single tree of suitable size using forward or backward algorithms.

In contrast to ordinary learning approaches which  construct one learner from training data, ensemble methods construct a set of learners known as base learners and combine them for prediction. Base learners are usually generated from the training data by a standard algorithm like decision tree, neural networks or other kinds of learning algorithms. There are two major kinds of ensemble methods based on the method of construction of the base learners. The first type is sequential ensemble methods where the base learners are generated sequentially and  the dependence among these base learners are used efficiently. The major example of this type  is the popular AdaBoost \citep{freund1997decision}.
The alternative type of ensemble learning  constructs the base learners in a parallel way with bagging \citep{breiman1996bagging}  as an example. Recently, ensemble trees have become a popular alternative as a predictive survival model building tool \citep{ishwaran2004relative,hothorn2005survival,hothorn2006unbiased,ishwaran2008random,zhou2015rotation} using the concepts of bagging and random forests \citep{breiman1996bagging,breiman2001random}.

  The Bayesian  classification and regression tree methods were developed by \citet{chipman1998bayesian,denison1998bayesian}.
\citet{clarke2008bayesian} use empirical Bayes priors and ensembles of Weibull tree models to model survival distributions.
The majority of Bayesian tree-based models in survival analysis are through ensembles of trees techniques (in a boosting way) using Bayesian additive regression trees (BART) \citep{chipman2010bart}. \citet{bonato2010bayesian} utilize BART to estimate latent parameters in a gamma process \citep{kalbfleisch1978non}, Weibull, and accelerated failure time models; \citet{henderson2017individualized} use BART in conjunction with a mean-constrained Dirichlet process mixture model to analyze heterogeneous treatment effects, which is reminiscent of \citet{kuo1997bayesian}; \citet{sparapani2016nonparametric} model survival times using a BART probit regression approach.

Even though machine learning algorithms usually focus on the predictive performance, the research in the field of interpretable machine learning \citep{molnar2020interpretable,doshi2017towards} has flourished in recent years due to its direct connection with artificial intelligence (AI) \citep{kuziemsky2019role}. For instance, AI is transforming the industry of healthcare rapidly and telehealth is one of its newest products. Telehealth utilizes digital information and communication technologies to access and manage healthcare services remotely. AI is booming in this industry of telemedicine where an efficient system for doctors is created to analyze, screen, and diagnose different disease conditions. However, a black-box machine-learning based predictive tool is not sufficient  and a more interpretable decision system is needed for this purpose. Hence, most of the ensemble methods (based on boosting or bagging logics) are not suitable in this framework.

To develop interpretable as well as flexible machine learning models for survival data,  we resort to a single tree model. One purpose of the proposed methodology is to provide a Bayesian tree framework which is simple to interpret (i.e.\ obtaining one tree) and flexible, i.e., moves beyond the proportional hazards model \citep{cox1972regression}. The  hazard function at each terminal node of the tree is modeled flexibly using an exponentiated Gaussian process. The posterior space of trees is efficiently searched via Markov chain Monte Carlo (MCMC) by obtaining the marginal likelihood of the observed data via a Laplace approximation \citep{tierney1986accurate}. The combination of interpretability coupled with flexibility makes this model well-suited for many problems, including subgroup identification and identifying prognostic and predictive biomarkers in medical research. In this Bayesian tree framework we sample from the posterior distribution of trees which provides a sample of multiple trees for decision making.

Most existing Bayesian tree survival models have transformed the survival problem into a regression or classification framework. For example,  \citet{bonato2010bayesian} use latent variables and parametric models. Tree ensembles are used primarily on the covariate structure in hierarchical specifications through a regression structure. \citet{sparapani2016nonparametric} transform survival responses into binary outcomes and use BART in this classification setup.  In contrast, we model the survival responses directly using  different hazard functions at each terminal node, and  construct the partition of the covariate space using a single tree structure similar to that of the original Bayesian classification and regression tree (CART) papers \citep{chipman1998bayesian,denison1998bayesian}.   However, our proposed model is different than existing CART models where simple parametric models (like Gaussian or multinomial models) with finite numbers of parameters and conjugate priors have been used at the terminal nodes of the tree. In contrast, we allow an unknown hazard function to model the  survival data within each terminal node of the tree.  An exponentiated Gaussian process is used as a prior for this hazard function which is  flexible but creates a non-conjugate posterior distribution. That way, explicit marginalization is not possible to obtain the marginal likelihood function which is the key step of the existing Bayesian Tree search algorithms for CART as well as BART  models. We propose a Laplace-approximation based method to develop a novel tree search algorithm for non-conjugate Bayesian tree models. In general, this new approach will aid the Bayesian tree models to be applied in more complex modeling scenarios.  Moreover, an important contribution of this paper is to   establish posterior consistency properties for our estimator under the tree partition model.

The paper is organized as follows:  Section 2 describes the model, Section 3 shows the theoretical properties of the model, Section 4 discusses the MCMC algorithm, Section 5 shows some simulated examples and an application to a liver survival dataset, Section 6 concludes.
\section{Survival models}
Let $t$ be a continuous nonnegative random variable representing the survival time of the subjects with the density function $f(t)$, the survival function $S(t)$ and the hazard function $h(t)= f(t) / S(t)$. In general, the hazard function depends on the survival time and the vector of covariates $\bx=(x_{1},x_{2},\ldots,x_{p})$. The Cox proportional hazards model \citep{cox1972regression} allow time and covariates to affect the hazard function separately, $h(t \mid \bx)=h_{0}(t){\rm exp}(\bx\bbeta)$, where $h_{0}$ is the baseline hazard and $\bbeta$ is a vector of regression parameters. This model has a strong assumption that the ratio of hazards for two individuals depends on the difference between their linear predictors at any time.

We extend this class of models by developing the Treed Hazards Model (THM) which allows a nonlinear effect of survival time and covariates in a flexible and adaptive way. In the THM model we split the covariate space of $\bx$ with tree partitions and a fit piecewise constant hazard model within each partition element. Let's denote the tree as $T$ with $M$ terminal nodes and the vector of  hazards within tree partition elements (terminal nodes) as $\bh=(h_{1}(\cdot),h_{2}(\cdot),\cdots,h_{M}(\cdot))$ where $h_{i}(\cdot)$ is the hazard
function at the $i$th terminal node.  This is in contrast to specifying a finite number of parameters at each terminal node as in existing Bayesian CART models. The THM model describes the conditional distribution of $t$ given $\bx$ in this tree structure. If $t$ lies in the region corresponding to the $i$th terminal node then $t \mid \bx$ has the distribution $f(t\mid h_{i}(\cdot))$. In the next section we develop the hierarchical tree model for the survival data.

\subsection{The Bayesian hierarchical treed hazards model }
Let $\by = \{y_1,\ldots,y_n\}$ be the survival time for $n$ subjects (realizations of $t$),  $\bdelta = \{\delta_1, \ldots, \delta_n\}$ be the corresponding censoring indicators  where $\delta_j$ is 1 if $y_j$ is observed and 0 if $y_j$ is right censored. Then, we can write the posterior distribution of the  model unknowns $T$ and $\bh$ as $p(T,\bh \mid \by,\bx,\bdelta)\propto p(\by, \bdelta \mid T,\bh)p(\bh \mid T)p(T \mid \bx)$.  In the following sections, we specify the likelihood function $p(\by, \bdelta \mid T,\bh)$ and the prior distributions $p(\bh \mid T)$, $p(T \mid \bx)$ which will be used to obtain the posterior distribution.

\subsubsection{Tree partition and the prior distribution $p(T \mid \bx)$}
First, we describe the tree partition $T$ on the covariate space and the related prior distribution. We assume that the data can be partitioned into $M$ partition elements $\by_1,\ldots,\by_M$ such that $\bigcup_{i=1}^M \by_i = \by$ and $\by_i \cap \by_{i^\prime} = \emptyset$ for all $i \neq i^\prime$.  We further assume that the data within the $i$th partition element are conditionally independent given the hazard function, $h_i(t)$, and assume that the $M$ hazard functions are mutually independent.  The tree partition is constructed through a series of recursive binary splits of the covariate space \citet{chipman1998bayesian}. The binary tree $T$ subdivides the predictor space marginally by choosing one of the $p$ covariates from $\bx$ and then defining a partition rule.

We adopt the same prior as \citet{chipman1998bayesian} which is constructed by splitting the tree nodes with a designated probability and assigning a rule each time a node is split.  Let $S_\ip$ equal 1 if the $\ip$th node is split, otherwise 0, $\ip=1,\ldots,2M-1$ where $M$ is the number of terminal nodes on a tree.  Let $R_\jp = (v_\jp,r_\jp)$, $\jp=1,\ldots,M - 1$ represent the rules of the $M -1$ internal nodes.  Here, $v_\jp \in \{1,\ldots,p\}$ indicates which variable is to be split, and $r_\jp$ is the splitting rule on the $v_\jp$th variable.  The partition rule for continuous variables is formed by choosing a threshold and then grouping observations based on whether or not their covariate value is less than or equal to, or greater than, the splitting threshold.  For categorical variables, the split rule is defined as a subset of the unique categories for that variable, and observations are grouped based on whether their covariate value is in the subset or not.

Let $T = (S_1,\ldots,S_{2M -1},R_1,\ldots,R_{M-1})$ be the collection of split indicators and rules for a particular tree. Then the prior for a tree, $T$, can be evaluated as
\begin{eqnarray} \label{eq:prior}
\pi(T) = \prod_{\ip=1}^{2M-1} \pi(S_\ip) \prod_{\jp=1}^{M-1}\pi(R_\jp).
\end{eqnarray}

Following \citet{chipman1998bayesian}, we define $\pi(S_\ip = 1) = 1 - Pr(S_\ip = 0) = \gamma(1 + d_i)^{-\theta}$ where $d_\ip$ represents the depth of the $\ip$th node (i.e.\ the root node has a depth of 0, child nodes of the root node have depth 1, etc.).  Here, $\gamma$ and $\theta$ are hyperparameters chosen by the modeller and influence the prior for the general size and shape of the tree. Here, $\pi(R_\jp) = 1/(m_\jp m^\star_\jp)$ where $m_\jp$ represents the number of variables on the corresponding node that are available for splitting and $m^\star_\jp$ is the number of splitting rules available on the selected variable for the node to which the rule belongs.  The splitting rules at a given node are selected from the observed data which reach that node from its ancestry.

\subsubsection{The Likelihood function $p(\by, \bdelta \mid T,\bh)$}
 The partition model assumes that given a tree, $T$, the hazard functions at each terminal node are independent and the data within a terminal node, given the hazard function at the node, are independent and identically distributed.  It is important to note that covariates are used to construct the tree, $T$, but are not used in the hazard functions.  Consequently, the likelihood function, given a tree, can be expressed as a simple product of likelihoods
\begin{equation}
p(\by, \bdelta \mid T, \bh) = \prod_{i=1}^{M} p(\by_i, \bdelta_i \mid h_i(\cdot)) \nonumber
\end{equation}
where $\by_i,\ \bdelta_i$ represent the $n_i$ responses and censoring indicators which belong to the $i$th terminal node, $p(\by_i, \bdelta_i \mid h_i(\cdot))$ represents the likelihood of $\by_i,\ \bdelta_i$ given the hazard function in the $i$th terminal node and $M$ is the total number of terminal nodes.

We write the density of an uncensored observation at the $i$th terminal node as $p_i(t)=h_i(t)e^{-\int_0^th_i(s)ds}$ for hazard function $h_i(\cdot)$. Let $p_c(\cdot)$ be the density of the censoring distribution, which has the same support as the true  densities $p_{0,i}(\cdot),\ i = 1, \ldots, m^*$ corresponding to true partition of the covariate space.

We place a Gaussian process prior on the log-hazard rate function on $[0,T_n]$, where $T_n=\max\{y_i\}$, $\log h_{i}(\cdot)=\beta_i+w_i(\cdot)$ and $w_i(t)$ follows Gaussian process with covariance  function $\tau_i^2K_{l_i}(s,t)=\tau_i^2\exp\{-|s-t|/l_i\}$, and let $\pi(\tau_i^2,l_i)=\pi(\tau_i^2)\pi(l_i)$ be the prior on the hyperparameters, and $\pi(\beta_i)$ be the prior on $\beta_i$.

\subsubsection{The hazards model $p(\bh \mid T)$ and the approximate likelihood}

In order to run a tractable Markov chain, a couple of alterations to the above model are used in practice.  First, an approximate likelihood and prior are used  where the domain is decomposed into small-width bins. Also, we normalize the domain between $[0,1]$ and use a mean zero Gaussian process under the proposed covariance kernel.

The approximation assumes the hazard functions are piecewise constant within each bin of the form $h_i(t) = \lambda_{ik}$ for $t \in (s_{k-1}, s_k]$ with $z_k$ denoting the midpoint of the interval.  In all, there are $K$ bins with endpoints $(s_{k-1},s_k]$ for $k=1,\ldots,K$ where $s_0 =0$ and $s_K = 1$.  For simplicity, we use the same grid for each hazard function and assume that the bins all have equal width, inducing a regular grid for $(z_1,\ldots,z_K)$.   The likelihood for data at a single terminal node under these assumptions can be expressed as
\begin{eqnarray*}
p(\by_i, \bdelta_i \mid \lambda_{i1},\ldots,\lambda_{iK}) &=& \prod_{j=1}^{n_i} \lambda_{ik_{ij}}^{\delta_{ij}} \exp\left\{ -\left[(y_{ij} - s_{k_{ij} - 1}) \lambda_{ik_{ij}} + I(k_{ij} > 1)\sum_{g=1}^{k_{ij} - 1} (s_g - s_{g-1}) \lambda_{ig}\right]\right\} \\
&=& \exp\left\{ \sum_{j=1}^{n_i} \delta_{ij} \log(\lambda_{ik_j}) - \left[\sum_{j=1}^{n_i} (y_{ij} - s_{k_{ij}-1})\lambda_{k_{ij}} + \sum_{k=1}^K (s_k - s_{k-1}) \lambda_{ik} m_{ik} \right] \right\} \\
&=& \exp\left\{ \sum_{k=1}^K \left[ n_{ik} \log(\lambda_{ik}) -  \lambda_{ik}\left( \sum_{j: k_{ij} = k} (y_{ij} - s_{k - 1}) + m_{ik}(s_k - s_{k-1}) \right)\right] \right\}
\end{eqnarray*}
where $y_{ij}$ is the $j$th response from the $i$th terminal node, $k_{ij}$ denotes the interval $y_{ij}$ belongs to, $m_{ik}$ is the number of observations of $\by_i$ which are greater than $s_k$, $n_{ik}$ is the number of uncensored observations of $\by_i$ which fall in interval $k$, and $\delta_{ij}$ is the censoring indicator for $y_{ij}$.

We model the hazard function $h_i(t) = \lambda_{ik} =  \exp(\omega_i(z_k))$ for $t \in (s_{k-1},s_k]$ where $\omega_i(t)$ is a zero-mean Gaussian process with covariance function $C_i(t,t^\prime) = \tau_i^2 \exp\{-|t - t^\prime|/l_i \}$.  This covariance kernel is chosen primarily for its computational advantages; see the appendix.  Given a set of $K$ locations, $\Bf_i = (\omega_i(z_1),\ldots,\omega_i(z_K))^T$  follows a multivariate-normal distribution with mean 0 and covariance matrix $\Sigma_i$ with the $i^\prime,j^\prime$th element of $\Sigma_i$ being $C_i(z_{i^\prime},z_{j^\prime})$. The posterior distribution of $\Bf_i$, assuming the covariance function parameters are known, is
\[
 p(\Bf_i \mid \by_i, \bdelta_i) \propto p(\by_i, \bdelta_i \mid \Bf_i) \pi(\Bf_i) = (2\pi)^{-K/2}|\Sigma_i|^{-.5} \exp\left\{g_0(\Bf_i) \right\}
\]
where
\[
g_0(\Bf_i) = \sum_{k=1}^K \left[ n_{ik} w_{ik} -  \exp(w_{ik})\left( \sum_{j: k_{ij} = k} (y_{ij} - s_{k - 1}) + m_k(s_k - s_{k-1}) \right)\right] -\frac{1}{2} \Bf_i ^T \Sigma^{-1}_i  \Bf_i
\]
and $w_{ik}$ is the $k$th element of $\Bf_i$.  Note that $g_0(\Bf_i)$ implicitly depends on the data, bins, and covariance function parameters.  Due to the assumption of independence between partitions, the full joint posterior distribution of a tree, $T$, with $M$ terminal nodes and $\Bf_1,\ldots,\Bf_M$ can be expressed as
\[
p(T,\Bf_1,\ldots,\Bf_M \mid \by, \bdelta) \propto \pi(T) \prod_{i=1}^M p(\by_i, \bdelta_i \mid \Bf_i) \pi(\Bf_i)
\]
where $\pi(T)$ is the prior on the tree as in \eqref{eq:prior}.

The goal in many partition models is to find a suitable partition structure.  In a Bayesian framework, this corresponds to searching the posterior of the tree, i.e.,
\begin{eqnarray}
p(T \mid \by) &=& \int \ldots \int p(T,\Bf_1,\ldots,\Bf_M \mid \by) d\Bf_1 \ldots d\Bf_M \nonumber \\
&\propto & \pi(T) \prod_{i=1}^M \int p(\by_i \mid \Bf_i) \pi(\Bf_i) d\Bf_i. \label{eq:marginal}
\end{eqnarray}
To do this, often a model is chosen in which all node-specific parameters can be integrated out analytically.  When the integrals have an analytic form, often an efficient reversible jump MCMC chain \citep{green1995reversible} can be constructed to search the posterior partition space. Unfortunately there is no analytical solution for the integrals in \eqref{eq:marginal} for the present model.  To overcome this issue, a Laplace approximation is used to approximate each integral. Laplace approximations have been used previously in various survival model settings. \citet{levine2014bayesian} utilize a Laplace approximation to search trees in a Bayesian frailty model and \citet{martino2011approximate} use integrated nested Laplace approximations via data augmentation to fit a Cox proportional hazards model to avoid MCMC. In our method, the Laplace approximation is similar to \citet{levine2014bayesian}, but we do not have random effects and we also optimize the covariance function hyperparameters via empirical Bayes.  \citet{payne2020conditional} utilize a Laplace approximation in a partition model framework for non-survival conditional density estimation using a Voronoi partition.

In usual form, the Laplace approximation approximates the integral with a normal distribution derived from a Taylor expansion of $\log[p(\by_i \mid \Bf_i) p(\Bf_i)]$.  The Laplace approximation has a closed form given by
\begin{equation} \label{eq:laplace}
\int p(\by_i, \bdelta_i \mid \Bf_i) \pi(\Bf_i) d\Bf_i \approx |D_i + \Sigma_i^{-1}|^{\frac{1}{2}}|\Sigma_i|^{-\frac{1}{2}} \exp\left\{g_0(\hat{\Bf}_i)\right\}
\end{equation}
where $D_i$ is a diagonal matrix with $j^\prime$th element $\exp(\hat{w}_{ij^\prime})\left[\sum_{j : k_{ij} = j^\prime} (y_{ij} - s_{j^\prime-1}) + (s_{j^\prime} - s_{j^\prime-1})m_{j^\prime}\right]$ and $\hat{\Bf}_i = \arg\max_{\Bf_i} p(\by_i \mid \Bf_i) \pi(\Bf_i)$, which can be obtained using Newton's method.

Thus far we have assumed that the covariance function parameters are fixed and known.  In practice, we employ empirical Bayes to set these values in each terminal node.  Specifically, we choose $\tau_i$ and $l_i$ to maximize $p(\by_i \mid \Bf_i) \pi(\Bf_i \mid l_i, \tau_i) \pi(l_i) \pi(\tau_i)$ where $\pi(l_i)$, $\pi(\tau_i)$ are t-densities with one degree of freedom and $\pi(\tau_i) = t(\tau_i / 10; 1)$ where $t(\cdot; 1)$ is the density of a t-distribution with one degree of freedom.

At this point, we have a way to approximate $p(T \mid \by)$, up to a normalizing constant, for a single tree utilizing Laplace approximations.  We now need to determine a way to search the approximate posterior.  As is common in many tree models such as \citet{chipman1998bayesian} and \citet{gramacy}, we employ Markov chain Monte Carlo methods.  Due to the varying number of parameters in the prior of the tree structure, a reversible jump Markov chain Monte Carlo algorithm \citep{green1995reversible} will be outlined in Section 4.

\section{Consistency results}
In this section, we show the density estimation consistency  when the true data generating density follows a partition model induced by a tree. If the survival function for the true data generating density has a tree partition form, then first we show that as $n$, the number of observations goes to infinity, the posterior distribution concentrates on  any weak neighborhood of the true data generating density with probability one. Then under slightly stronger  conditions, we establish the $L^1$ consistency.

Let  $P_1^*,\dots,P_{m^*}^*$ be the true partition of the covariate space denoted as  $\mathscr{P}^*=\{P^*_1,\dots,P_{m^*}^*\}$, and  let $h_{0,i}(t)>0,\ t>0$ for $i=1,\dots,m^*$,  be the true hazard function for partition $P_i^*$ and the corresponding true uncensored data generating  density for the $i$th partition is given by  $p_{0,i}(t)=h_{0,i}(t)e^{-H_{0,i}(t)}$ where $H_{0,i}(t)=\int_0^th_{0,i}(s)ds$. As mentioned earlier, $p_c$ is the density of the censoring distribution. Let $F_c$ and and $F_{0,i}$ be corresponding cdfs and $1-F_{0,i}(t)=e^{-H_{0,i}(t)}$. For a generic partition $\mathscr{P}=\{P_1,\dots, P_{k'}\}$ of the covariate space and the partition densities $p_{i}(\cdot)$,  for observation $y,\ {\bf x}$, and the censoring indicator function $\delta$, the data generating density $p(\cdot)$  is given by
\begin{eqnarray}
p(y,\delta \mid {\bf x},\mathscr{P})=\sum_{i=1}^{k'}\left[p_{i}(y)\int_y^\infty p_c(s)ds{\bf I}_{\delta=1}+p_c(y)\int_y^\infty p_{i}(s)ds{\bf I}_{\delta=0}\right]{\bf I}_{{\bf x}\in P_i}.
\label{likelihood1}
\end{eqnarray}
We observe  $(y_1,\delta_1),\ldots,(y_n,\delta_n)$ at observed covariates ${\bf x}_1,\ldots,{\bf x}_n$, where ${\bf x}_j$ are i.i.d. for $j=1,\dots,n$. Let  $\#(T)$ denote the number of partitions/terminal nodes for a tree $T$. Denoting the generic prior distribution by $\pi(\cdot)$, for $T_n=\max_i\{y_i\}$,  from Section 2 we have the following prior: $P_1,\dots,P_{\#(T)}\sim\pi(T)$; $\log h_i(t)=\beta_i+w_i(t) \sim  GP(\beta_i,\mathscr{K})$;    $mathscr{K}(s,t)=\tau_i^2 K\left(\frac{s}{l_i},\frac{t}{l_i}\right),\ t\in [0,T_n];\ i=1,\cdots,\#(T)$;
$\pi(\tau^2_i,l_i)=\pi(\tau^2_i)\pi(l_i)$; $\beta_i\sim\pi(\beta_i)$.
\begin{theorem}
For any bounded continuous function $g$, and for $\epsilon>0$, let $U^{\epsilon}_g=\{p:|\int g p^*(\cdot)-\int g p(\cdot)|<\epsilon\}$, be a weak neighborhood of the true joint density $p^*(\cdot)$. Under the model in equation \ref{likelihood1}, prior in equation \ref{prior}, under assumptions [A1]--[A6] in supporting materials, and under the true data generating distribution $p^*$, $\Pi_n(U^{\epsilon}_g\mid \cdot)\rightarrow 1$ almost surely as $n\rightarrow \infty$, where $\Pi_n(\cdot)$ denotes the posterior distribution based on $n$ observations.
\label{thm1}
\end{theorem}

Theorem \ref{thm1} establishes consistency in  weak topology with $\Pi_n(\mid \cdot)$ being weakly consistent at $\{ {p_{0,i}}_{i=1,\cdots,m^*},p_c,\mathscr{P}^* \}$. Even though  we state  Theorem \ref{thm1} under the support condition given in [A3], the consistency result will hold true for more  general support  under additional mild conditions, given in the following. We can now state the following strong consistency result. For the strong consistency we consider estimated densities that are supported on $[0,T_n]$ and show that they will remain in a small $L^1$ neighborhood of the true density with probability one.

\begin{theorem}

Let $\epsilon>0$ and $U^{\epsilon}=\{p:\int | p^*(\cdot) -p(\cdot)|<\epsilon\}$, be a $L^1$ neighborhood of the true joint density $p^*(\cdot)$. Under the model in equation \ref{likelihood1} and related prior distributions, under assumptions [A2],[A4]--[A5],[B1]--[B5] in supporting materials and under the true data generating distribution $p^*$, $\Pi_n(U^{\epsilon}\mid \cdot)\rightarrow 1$ almost surely as $n\rightarrow \infty$, where $\Pi_n(\cdot)$ denotes the posterior distribution based on $n$ observations.
\label{thm2}
\end{theorem}
Technical details and related proofs of the results are given in  Appendix ii.

\section{MCMC}
Since the dimension of the tree parameter space is variable, we use a reversible jump MCMC algorithm to search the posterior parameter space  \citep{green1995reversible}. The algorithm is inspired by those of \citet{chipman1998bayesian}  and \citet{gramacy} with a few modifications. The primary goal in constructing the algorithm is to obtain draws from $\tilde{p}(T \mid \by)$, the approximate posterior obtained using the Laplace approximations.  This is accomplished in the reversible jump MCMC algorithm by modifying an existing tree with one of four possible moves: grow, prune, change, and swap.

\begin{itemize}
\item {\bf Grow:} Randomly select a terminal node and add two child nodes based on a split rule.  The split rule is selected by first randomly selecting a variable and then randomly selecting an available rule of that variable (i.e.\ the rule is drawn from the prior).  This increases the number of terminal nodes of the tree by one, thereby ``growing" the tree.

\item {\bf Prune:} The opposite of a grow step.  Randomly select a node whose children are both terminal nodes and delete the split rule and the children nodes.  This reduces the number of terminal nodes of the tree by one, thereby ``pruning" the tree.

\item {\bf Change:} Randomly select an interior node and change the existing split rule.  If the current rule is based on a categorical variable, a new rule is drawn randomly from the prior by first randomly selecting a variable, and then randomly selecting a rule from that variable.  If the current rule is based on a numeric variable, there are two possible changes that can be made with pre-specified probability:  1) Move the current rule to the next largest or smallest rule available for the current variable, 2) Draw a new rule from the prior.

\item {\bf Swap:} Randomly choose a parent and child pair which are both internal nodes and swap their rules.  If, however, the parent and child's rules are based on the same continuous variable, a rotate step is employed \citep{gramacy}.  If both children of the parent node have the same rule, both children swap their rules with the parent node.
\end{itemize}

The MCMC algorithm proceeds as follows: 1. Initialize a starting tree, $T$. 2. Propose a new tree, $T^1$, from the current tree using a grow, prune, change, or swap step (as appropriate based on current tree). 3. Accept $T^1$ with probability $\alpha = \{1,\frac{q(T \mid T^1) \tilde{p}(T^1 \mid \by)}{q(T^1 \mid T) \tilde{p}(T \mid \by)}\}$
where $q(\cdot \mid \cdot)$ is the proposal distribution and $\tilde{p}(T \mid \by)$ is the approximate posterior of $T$ obtained from the Laplace approximation; otherwise keep the current tree, $T$. 4. Repeat steps 2-3 until the stationary distribution is reached and the desired number of samples is obtained.

In practice we replace the approximate posterior of the tree in $\alpha$ with the likelihood-prior product in \eqref{eq:marginal} since the normalizing constants cancel.

From their inception, Bayesian tree models have been notorious for their difficulty in mixing when using MCMC procedures.   Various methods have been proposed to overcome this issue including restarting the MCMC chain multiple times \citep{chipman1998bayesian}, parallel tempering \citep{gramacy2010categorical}, and multiset sampling \citep{leman2009}.  We implement parallel tempering \citep{geyer1991markov} modeled after \citet{gramacy2010categorical}.  Additional technical details regarding the MCMC, are provided in the appendix.

\section{Simulations \& Applications}
\subsection{Preliminaries}
For the following analyses, we use the same set of tuning parameters for the MCMC chain.  We set the probability of grow, prune, change, and swap steps to be .25.  For the change step for a continuous variable, we choose to move the rule to an adjacent rule with probability .75, otherwise a rule is drawn from the prior.  We run each analysis with 8 parallel chains with the minimum inverse-temperature of .1, and propose a trade between the parallel chains every 10 MCMC iterations. Trees are required to have a minimum of 25 observations at a terminal node.
\subsection{Tree partition with non-proportional hazards} \label{sec:tp}
Our first simulation demonstrates the flexibility and interpretability of the THM.  One thousand data points were simulated from the following model:

$X_1   \sim  \text{Uniform}(0,10)$; $X_2 \sim  \text{Multinomial}(\text{A,B,C,D}; .25,.25,.25,.25); $Y  \sim  \text{Weibull}(5,2)  \text{ if } X_1 > 5 \text{ and } X_2 \in \{\text{A,B}\}$;
$Y  \sim \text{Weibull}(1,5)  \text{ if } X_1 \leq 5 \text{ and } X_2 \in \{\text{A}\} $; Y \sim $\text{Weibull}(.5,.9)  \text{ if } X_1 \leq 5 \text{ and } X_2 \in \{\text{B}\}$;
$ Y \sim $\text{Weibull}(5,5)  \text{ if } X_1 \leq 3 \text{ and } X_2 \in \{\text{C,D}\}$;
$ Y\sim $\text{Weibull}(.5,.5) \text{ if } 3 \leq X_1 \leq 7 \text{ and } X_2 \in \{\text{C,D}\};$$
$Y\sim G  \text{ if } X_1 > 7 \text{ and } X_2 \in \{\text{C,D} \}$

where $G$ is a random variable drawn from an arbitrary piecewise linear survival function.  The data were also right censored with a 5\% censoring rate. The Markov chain Monte Carlo algorithm ran for $10^4$ iterations.  The tree with the highest posterior probability from the chain is shown in Figure \ref{fig:sim9tree} along with the posterior and true survival curves.  Note that since we are using observed data values to choose the rules at each node, the rules for numeric variables will not be exactly the same as those used to generate the tree.  This example highlights the strength of the THM, simultaneously producing easy interpretation and flexibility.  Interpretation is simple since only one tree is needed, but the the various shapes of the survival curves can still be modeled flexibly at each terminal node.

\begin{figure}
\begin{center}
  \subfigure[]{\includegraphics[trim=40 30 40 10, clip,width = .45\textwidth]{figs/sim9tree.eps}}
  \subfigure[]{\includegraphics[trim=0 10 0 10, clip,width = .45\textwidth]{figs/sim9survival.eps}}
  \caption{(a) The tree with the highest posterior probability from the model in Section \ref{sec:tp}.  (b) The posterior mean (dotted), 95\% credible intervals (dashed), and the true survival curve (solid) in each partition element.}
\label{fig:sim9tree}
\end{center}
\end{figure}

\subsection{Cox proportional hazards} \label{sec:lph}
Although the present partition model assumes the data are generated from a tree structure, it still performs well when this assumption does not hold, e.g., when dealing with the classic linear proportional hazards model \citep{cox1972regression}.  One thousand observations were simulated from an exponential proportional hazards model with an 8\% censoring rate.  The model is $h(t \mid \mathbf{X} = \mathbf{x}) = \frac{1}{2} \exp\left\{\mathbf{x} \mathbf{\beta} \right\}$,       \mathbf{X} = (X_1,\ldots,X_6) $
$X_1,\ldots,X_5  \sim\text{Uniform(0,1)}
,$X_6  \sim\ & \text{Bernoulli(.5)} $
$\beta =\ & (-1, 1, 2, 0, 0, -2)^T$ where $h(\cdot \mid \cdot)$ represents the hazard function.

Note that we have two spurious covariates ($X_4,X_5$) which do not influence the hazard function.  The MCMC algorithm ran for $10^4$ iterations.  In the tree with the highest posterior probability (8 terminal nodes), neither of the two spurious covariates were included.  To give a measure of performance, Figure \ref{fig:sim7survival} plots the posterior survival curve and the exponential survival curve using the median true hazard value of the data that reach the terminal node.  This comparison serves as an indicator of where the survival distribution should lie in that region of the covariate space, and indicates the model captures the important features of the curves well.

\begin{figure}[t]
\begin{center}
\includegraphics[scale=.55]{figs/sim7survival.eps}
\caption{The posterior mean (dotted) and 95\% credible intervals (dashed).  The survival curve corresponding to the median true hazard function of the data in each terminal node is also plotted (solid).}
\label{fig:sim7survival}
\end{center}
\end{figure}

\subsection{Prognostic \& predictive biomarker, comparison with BART, survival random forests, and Cox}
Identifying predictive and prognostic biomarker effects is important in drug development.  Sometimes there is a hypothesis that a certain sub-population, based on a predictive biomarker, will respond better to a treatment.   Other times the value of a prognostic biomarker is related to a patient's prognosis.  It is possible a biomarker may be both prognostic and predictive.

A single biomarker, $b$, and a treatment assignment indicator, $a$, are simulated for 1,000 subjects.  The biomarker is both prognostic and predictive, and there are two treatment arms. The time to event, $y$, is simulated from a Weibull model with random right-censoring in 7.8\% of subjects,   $a \sim Binomial(1, .5)$, $b \sim Uniform(0, 1)$
    $y \sim Weibull(1 + 2ab, 1 + 5b)$
where $Weibull(\lambda, k)$ is a Weibull distribution with scale $\lambda$ and shape $k$.  For small biomarker values $b$ nearer to 0, there is a small treatment effect, whereas for values of $b$ closer to 1, there is a large treatment effect creating a predictive biomarker.  For a fixed treatment assignment $a$, the shape of the survival curve is a function of $b$, creating a prognostic effect for the biomarker.

The tree partition model was fit to the simulated data, and the tree with the highest posterior probability is plotted in Figure \ref{fig:pred_prog}. The tree and the estimated survival curves at each terminal node capture the predictive and prognostic effects.  The split rule of the root node indicates that there is a difference in survival functions between the two treatment arms.  Descending either right or left from the root node, all remaining splits are on the biomarker, indicating a prognostic biomarker effect.  For two patients with similar biomarker values, but different treatment assignments, their survival curves can be quite different.  For instance, a patient with a biomarker value of .8 has a much better survival curve with treatment assignment 1 verses 0; this indicates that the biomarker is also predictive.

The simulated data were fit to three additional models: a frequentist Cox proportional-hazard model with a non-parametric baseline hazard an an interaction between $a$ and $b$; Bayesian additive regression trees (BART) \citep{sparapani2016nonparametric}; and survival random forests \citep{ishwaran2008random}.  The survival curves, the 95\% confidence/credible intervals, and the true treatment effect for a patient with $a = 1$ and $b = .10$ are plotted in Figure \ref{fig:pred_prog_comparison}; for the survival random forest method, the 95\% intervals are calculated from the predictions of all trees.  The changing shape of the survival curve is not captured well in the Cox model since the proportional hazard assumption is too restrictive.  This is problematic; even if the Cox model identifies a significant treatment by biomarker interaction, which is often used to determine a predictive biomarker effect, the estimates of the survival curve may be inaccurate.  On the other hand, BART and random forests estimate the survival curves well, but often custom post-processing of the results in a potentially large-volume covariate space are required to better understand which predictors are prognostic, predictive, or both.  The THM automatically provides both clear interpretation and flexible estimation of the survival curves.
\begin{figure}[t]
    \centering
    \includegraphics[scale = .55]{figs/pred_prog.eps}
    \caption{Tree with the highest posterior probability with the posterior mean (dotted) and 95\% credible intervals (dashed) at each terminal node.}
    \label{fig:pred_prog}
\end{figure}
\begin{figure}[t]
    \centering
    \includegraphics[scale = .55]{figs/pred_prog_comparison2.eps}
    \caption{The true survival curve (solid) for a patient with $a = 1$ and $b = .10$, estimated mean (dotted) and 95\% confidence/credible intervals (dashed) for the frequentist Cox, BART, and THM.}
    \label{fig:pred_prog_comparison}
\end{figure}
\subsection{Analysis of Primary Biliary Cirrhosis Data}

Primary biliary cirrhosis survival data collected from 1974-1984 by the Mayo Clinic and obtained via the R survival package \citep{survival-package} were fit using the THM.  Covariates included treatment group (D-penicillamine or placebo), age, sex, edema status (0 no edema, 0.5 untreated or successfully treated, 1 edema despite diuretic therapy), serum bilirunbin (mg/dl), serum albumin (g/dl), and standardised blood clotting time.  Only subjects who were randomized to treatment and had no missing values were included in the analysis for a total of 276 subjects.  The outcome of interest was time to death, and patients were censored at the end of study or if they received a transplant, yielding a censoring rate of 60\%.

Ten-thousand MCMC samples were drawn and the tree with the highest posterior probability is shown in Figure \ref{fig:pbc} with the survival curves plotted at each terminal node.  The black curves represent the posterior mean and 95\% credible intervals from the tree model, and the red lines provide mean and 95\% confidence intervals from the Kaplan-Meier curves.  The plot indicates that the THM is very flexible, and close to the Kaplan-Meier estimates.

The tree indicates that patients with bilirubin levels less than 1.9 have the highest survival rates.  For subjects with no edema and bilirubin levels greater than 1.9 and less than or equal to 5.7, patients under the age of 45 had better survival.  Patients with bilirubin levels greater than 1.9 and edema had the worst survival.  These conclusions are consistent with the findings of \citet{dickson1989prognosis}, who used a Cox proportional hazards model on primary biliary cirrhosis data, however, the tree method provides greater flexibility in estimating the survival curves, especially in this case where the hazards are not proportional.  Last, we note that the tree did not split on treatment status, which is consistent with the findings of \citet{gong2004d} who performed a systematic review of the effect of D-penicillamine on patients with primary biliary cirrhosis and concluded that it does not decrease the risk of mortality.

\begin{figure}
    \centering
    \includegraphics[scale=.55]{figs/pbc_km.eps}
    \caption{The tree with the highest posterior probability and the posterior mean (black, dotted), 95\% credible intervals (black, dashed), and the Kaplan-Meier estimates (red) at each terminal node.}
    \label{fig:pbc}
\end{figure}

\section{Conclusion}
The main advantage the tree hazards model (THM) is its ease in interpretation coupled with its flexibility.  Inference is easily obtained from inspecting posterior trees and the flexibly estimated survival curves at the terminal nodes.  In both simulations and applications, the model performs well in identifying interpretable partitions and flexibly estimating survival curves.  Furthermore, we have established weak consistency for the partition model under the proposed kernel, where we can accommodate flexible forms of the hazard function under moderate conditions.

\bibliographystyle{apa}

redpajama_set_name: RedPajamaArXiv


Coming to the American Ballet Theater's production of Onegin from a music and opera background, and a literary as well, the advantages of dance are immediately apparent.
Physical expression of emotion is clear. Although the choreographic vocabulary is not extensive in the classic form, the opportunities for imaginative creation of scenes which exist only in the character's heads are abundant.
David Hallberg, returning from a serious injury, dances Onegin. Hee Seo is a fragile but firm Tatiana.
The letter scene, missing the delicious Tchaikovsky music here, is breathtaking. Onegin is imagined by Tatiana behind a mirror. All the feeling that she has poured into her letter is now danced to the actual presence of the absent Onegin. The orchestra plays Tchaikovsky, but there are only glimpses of the Onegin music. Instead choreographer John Cranko has set the dance to a pastiche of the composer's other music.
Always the dramatists who take the poem on are compelled to make Onegin at the least insensitive. If we haven't read the poem, we don't know that he has freed his serfs and is a well-read man whose library is filled with books of poetry, history and art.
The corps de ballet is exciting. Out dancing in the woods at one moment we see them fly by, leaping across the stage and back again, a whirl of legs and arms and flying skirts. At a ball, they capture Pushkin: laughter, curtsy, bow, mazurka, gallop, and waltz.
After the ballroom scene, in which Onegin reappears in Tatiana's life, the push and pull of desire is danced in the most powerful moments of the evening. Tatiana and Onegin seem bound by an elastic chord, pulling apart from each other and thrust back again, in multiple positions, always returning and then pushing back. We may know the outcome, but the rush of feelings in this dance leaves us proud of a woman who will do her duty and keep her vows and yet is saddened by the loss of love.
The sets are marvelous. Yes, we have the white birches which every designer has decided signify Russia. Tatiana's country home, indicated by full columns stage right and a suggestive wall on the left, is comfortable and pleasing. A wonderfully domed ceiling filled with stained glass is used for ballrooms. A smaller dome of vines and leaves suggests a lover's nest. Magnificent candelabras descend for the final ballroom scene.
The forest where the duel takes place is created by dark lighting amidst the birch trees. Taking the place of two male 'seconds' during the duel, two female figures implore the duelists to stop their fight. In their touching huddled and then outstretched movements to and from the characters, the urgency of their request is dramatically conveyed.
While you miss the themes of Tchaikovsky's take on the drama, if you put aside the often bland musical choices made, this ballet version of the story is thrilling. One even imagines the last scene being projected during an opera version to give physical form to the final confrontation between Tatiana and Onegin.


Westmont and its affiliated law firm, Fredric Marro and Associates, P.C. provide merger and acquisition services on a national basis. Due to our unique preparation approach and assessment of multi-state requirements, our services are organized, concise and more economical. Acquisition services begin with pre-acquisition assessment and notification, Form A preparation and filing and continue through post-acquisition fulfillment. With Westmont, you will not be paying to train our staff. The average years of staff experience exceeds 20 years. Our team is comprised of lawyers, accountants, actuaries, former regulators and regulatory specialists as required.