Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
166,509 characters · 25 sections · 139 citation commands
Seemingly unrelated Bayesian additive regression trees for cost-effectiveness analyses in healthcare
Many research questions in health economics are concerned with trading off the costs and benefits of a medical intervention. The most prominent examples are in cost-effectiveness analysis (CEA), where we wish to decide whether a new innovative treatment is worth the associated increase in costs compared to usual care. We therefore need to estimate the average treatment effects on both costs and health. However, in order to get coherent measures of uncertainty, the two treatment effects must be estimated jointly in order to account for the correlation between them baio2012bayesian. Furthermore, if the CEA is performed with observational data, where the treatment assignment is not randomised, we additionally have to adjust for confounding bias in the analysis. These points are elaborated in Section (ref).
In this paper, we aim to estimate the cost-effectiveness of a novel treatment for physical trauma rehabilitation, called the transmural trauma care model (TTCM), using data gathered under a study by wiertsema2019cost which we will henceforth refer to as the TTCM data. The treatment assignment is not randomised and the number of potential confounders is large relative to the sample size. Of the multiple cost and effectiveness outcomes the authors investigated, we focus on healthcare-related costs and health-related quality of life. It is of interest to jointly estimate both outcomes, and reasonable to assume both that the outcomes are non-linearly related to the available predictors and that each outcome may depend on different subsets of predictors, which may in turn interact in complex ways. The challenge of CEA under these circumstances motivated the development of our novel methodology, though we also anticipate its use in other CEA studies and broader healthcare settings.
In the causal inference literature, there is wide agreement that flexible nonparametric methods, which do not impose strong parametric assumptions on the regression functions, are the best tools for estimating treatment effects with observational data dorie2019automated,rudolph2023all. It is not straightforward, however, to apply this knowledge in the context of CEAs. Seemingly unrelated regression models zellner1962efficient, the most recommended statistical method for CEAs willan2004regression,el2022scoping, impose strong linearity assumptions, which can bias the inferences if there are strong non-linear relationships between the variables of interest. On the other hand, there is a distinct lack of nonparametric regression methods which can handle multivariate outcomes. We adopt a Bayesian perspective and fill this gap by developing a nonparametric version of SUR. The idea is to replace the linear predictors in the SUR model by sums of regression trees. In the univariate case, this regression method has become known as Bayesian additive regression trees (BART). BART has already demonstrated competitive performance for univariate responses dorie2019automated,rudolph2023all and it seems plausible that this efficacy will extend to situations with multiple outcomes of interest. However, by embedding BART in the SUR framework, we also seek to overcome some limitations of existing multivariate BART extensions.
chipman2010bart introduced BART as an ensemble method, where each learner is a tree following the Bayesian CART approach previously proposed by the same authors chipman1998bayesian. Considering a univariate response vector $\mathbf{y} \in \mathbb{R}^{n}$ and a set of predictors $\mathbf{X} \in \mathbb{R}^{n\times p}$, which may be of mixed type, the main objective of regression modelling is to estimate the conditional expectation $\mathbb{E}\lbrack y_i\:\vert\:\mathbf{x}_i\rbrack$, also called the regression function. As a nonparametric model, BART offers high flexibility in estimating this conditional expectation. The tendency of decision trees to overfit is mitigated by the Bayesian underpinnings of the method, allowing informative prior distributions to impose regularisation in a principled and transparent fashion, as well as the additive nature of the ensemble. These features facilitate better generalisation, in a similar vein to the gradient boosting approach of friedman2001greedy. The success of BART is widely reported in the literature across a broad spectrum of applications janizadeh2021novel,sarti2023bayesian,yee2023evaluating. In addition, theoretical work has demonstrated the frequentist optimality of BART under certain conditions linero2018bayesian,rockova2020semi,rockova2020posterior.
The BART model was originally developed for univariate responses but subsequent formulations emerged to adapt BART to multivariate responses. Examples include the Bayesian additive vector autoregressive tree florian2022vabart for multivariate time-series and formulations by peruzzi2022spatial for multivariate spatial data. Other applied work involves an extension of BART to forecast the tails of multivariate responses clark2023tail. um2023bayesian adapted BART to cover not only multivariate responses but also the assumption of a skew-normal distribution; throughout this paper, we refer to the variant without skewness as mvBART. Furthermore, mcjames2024bayesian proposed an extension of Bayesian causal forests hahn2020bayesian for multivariate responses. All aforementioned approaches share the limitation that the tree structure must be identical for each component of the outcome vector. It is easy to envisage situations for which this is inappropriate: a specific covariate may be strongly associated with one outcome but independent of another; e.g., factors which govern costs may be unrelated to quality of life, and vice versa.
In this study, we propose a novel variant of BART termed seemingly unrelated BART (suBART) which is designed to handle multivariate continuous responses and address this key limitation. chipman2010bart previously drew parallels to SUR models zellner1962efficient and alluded to the potential for extending BART in this direction. Our framework differs from the aforementioned multivariate BART extensions, which assume a single set of trees with correlated multivariate terminal node parameters. Instead, we jointly fit individual ensembles of trees for each outcome and model their interdependence through correlated error terms. Thus, suBART also differs from merely applying univariate BART models to each outcome independently. Motivated by our investigation into the cost-effectiveness of the TTCM intervention, we further extend suBART to incorporate propensity scores, in the spirit of hahn2020bayesian, as befits causal analyses. Beyond CEA settings, we also develop probit suBART, an extension of suBART to accommodate multivariate binary outcomes. While not the main focus of this paper, we envision this version of the model being useful in marketing applications with correlated binary outcomes; see the examples in rossi2005bayesian. Our approach is similar to that of chakraborty2016bayesian, who offers a version of seemingly unrelated BART for exclusively continuous outcomes with an adaptive number of trees. However, it is not specifically tailored to causal inference objectives typical of CEA and lacks an available open-source software implementation. Moreover, the non-informative inverse-Wishart prior adopted by chakraborty2016bayesian on the covariance matrix of the error terms fails to properly calibrate the prior on the residual variance parameters, which is crucial for regularising the contributions of the trees in the ensemble. We address these gaps by providing a comprehensive framework, which covers either continuous or binary outcomes, with appropriately specified priors, and a practical implementation through an R package named subart\footnote{Implemented using Rcpp and available at: \url{https://github.com/MateusMaiaDS/subart}.}.
This article proceeds as follows: Section (ref) provides background theory on CEA and motivates the development of the suBART model in the context of the application to the TTCM data from wiertsema2019cost. Section (ref) then elaborates on the theoretical underpinnings of the suBART methodology and Section (ref) discusses the posterior inference for both the multivariate continuous and multivariate binary outcome settings. Section (ref) considers different simulation scenarios to compare the predictive performance of suBART against standard competitors. The empirical findings of our application of suBART to the TTCM data are presented in Section (ref), along with an additional simulation experiment based on these data which assesses suBART's estimation of treatment effects in the CEA setting. Finally, Section (ref) summarises the proposed methodologies, highlighting both their limitations and potential for further extension, and presents conclusions regarding the health economic application. Additional results, comparisons, and findings are included in the appendices.
Our main motivation in developing suBART is its potential applicability in cost-effectiveness analyses of healthcare treatments, particularly in cases where the treatment assignment is not randomised. When estimating treatment effects from such observational data, simple parametric models often lead to severe bias and flexible nonparametric models are preferable hernan2024causal. Naturally, such concerns extend to CEA settings, where we want to estimate two treatment effects simultaneously. However, there is a lack of statistical methods fit for such purposes. Given that BART has proven to be very useful for univariate causal inference tasks hill2011bayesian,dorie2019automated,rudolph2023all, we consider it a promising basis for multivariate CEA.
We continue by sketching some relevant ideas from health economics and causal inference in Section (ref), drawing on the detailed treatments found in gabrio2019bayesian and li2023bayesian. Subsequently, we give a detailed presentation of the TTCM data in Section (ref).
The fundamental problem of cost-effectiveness analyses in health economics is to determine which of two competing healthcare treatments --- usually, but not always, for the same disease --- should be implemented. Often one is both more effective and more expensive than the other, which raises the question whether the increase in health is worth the added expenses. Henceforth, we let $c_i$ and $q_i$ respectively denote the healthcare costs and the health-related quality of life associated with a patient. We also suppose that there are two different treatments of interest, $t = 0$ and $t = 1$. Using the usual potential outcomes notation, we let $c_i(t)$ denote the costs associated with patient $i$, had they received treatment $t$. We do likewise for $q_i(t)$. We furthermore suppose that there is some vector of baseline characteristics $\mathbf{x}_i$ which have an effect on both the outcomes $c_i$ and $q_i$, as well as the treatment indicator $t_i$.
Given a sample of $n$ observations, we wish to estimate the mixed average treatment effect (MATE)\footnote{The MATE is closely related to the population average treatment effect (PATE), although the terminology for treatment effects is not consistent across the literature. We elect to use the terminology of li2023bayesian, according to whom the PATE requires a generative model for the covariates $\mathbf{x}_i$ and thus necessitates additional assumptions. We work with the MATE in Equation (ref) to simplify the analyses. Note that this common approach actually corresponds to what is called the PATE by imbens2015pate and other authors.} on the costs in this sample, which we define as
We now make the assumption of ignorability, which means that conditional on the baseline covariates $\mathbf{x}_i$, the treatment $t_i$ is independent of the potential outcomes $c_i(0)$ and $c_i(1)$. Under this assumption, we may rewrite our treatment effect as \[\Delta_c = \frac{1}{n} \sum^n_{i = 1} \mathbb{E}\left\lbrack c_i\:\vert\:t_i = 1, \mathbf{x}_i\right\rbrack - \mathbb{E}\left\lbrack c_i\:\vert\:t_i = 0, \mathbf{x}_i\right\rbrack.\] Thus, $\Delta_c$ is completely specified by the conditional expectations $\mathbb{E}\lbrack c_i\:\vert\:t_i, \mathbf{x}_i\rbrack$. We proceed in the same manner for $\Delta_q$, the MATE on the patient's quality of life. It is then customary to combine the two treatment effects into a utility function, the incremental net benefit (INB): \[\text{INB}_{\lambda} \coloneqq \lambda \Delta_q - \Delta_c.\] The scalar parameter $\lambda$ is called the willingness-to-pay. Roughly speaking, $\lambda$ quantifies how much cost (in the given currency) a decision-maker is willing to trade for a one-unit increase in healthcare-related quality of life for one patient. The decision rule is then simple: if the INB is at most zero, we say that treatment $1$ is not cost-effective, and treatment $0$ should be implemented. If the INB is larger than zero, we consider treatment $1$ to be cost-effective and worthy of being implemented. Throughout this paper, we give $\lambda$ values in units of €1,000.
To illustrate the importance of modelling the two outcomes $c$ and $q$ jointly, let us assume for simplicity that the joint distribution of $\Delta_c$ and $\Delta_q$ is bivariate normal. Then
Consequently, the probability of cost-effectiveness depends on the covariance of $\Delta_c$ and $\Delta_q$. Without the normality assumption, this probability can usually not be found explicitly, but the same principle applies nonetheless: the probability of cost-effectiveness depends on the joint distribution of $\Delta_c$ and $\Delta_q$ lothgren2000definition, gabrio2019bayesian. It follows that we must model $c_i$ and $q_i$ jointly, as modelling them separately would enforce the unrealistic assumption that the treatment effects $\Delta_c$ and $\Delta_q$ are independent. This belief is seldom appropriate, since empirical cost and health data are often strongly correlated willan2004regression. We expound on these points through simulation experiments in Appendix (ref).
We hence use the suBART model developed below to jointly estimate the conditional expectations $\mathbb{E}\lbrack c\:\vert\:t, \mathbf{x}_i\rbrack$ and $\mathbb{E}\lbrack q\:\vert\:t, \mathbf{x}_i\rbrack$. The treatment effects and INB can then be obtained as functions of these estimates. Our approach mirrors that of hahn2020bayesian: we first estimate propensity scores (using probit BART), and then condition the suBART model on all covariates, the treatment indicator, and the estimated propensity scores.
Our data of interest, pertaining to $n=140$ patients suffering from traumatic injuries, were originally collected and analysed by wiertsema2019cost. The two treatment options are usual care and the novel Transmural Trauma Care Model (TTCM), denoted by $t=0$ and $t=1$, respectively. The treatment assignment was not randomised. The outcomes we use, for $c_i$ and $q_i$ respectively, are the costs from the healthcare perspective and generic healthcare-related quality of life. As is usual in CEAs, the former is an aggregate measure: $c_i$ comprises the total costs acquired from hospital records as well as several questionnaires conducted over the course of nine months following treatment, in which patients were surveyed on their use of various healthcare resources. The responses --- which relate to issues such as hospital stays, medication use, and surgeries --- were then converted to costs. Conversely, the effectiveness outcome $q_i$ was calculated from one single survey administered nine months after treatment using the EQ-5D-3L instrument lamers2006dutch. The data also include $p=11$ baseline covariates, with respective sample sizes of $83$ and $57$ in the two treatment groups. The ratio of covariates to observations is thus reasonably large. We reproduce the table of baseline variables in Table (ref).
The original dataset had some missing observations --- for survey items related to the outcome variables only --- which wiertsema2019cost dealt with through multiple imputation. $17\%$ of patients did not complete any follow-up questionnaires, and hence were missing all information on $q_i$ and some survey information on $c_i$ (though hospital records were available for all patients). Additionally, $39\%$ and $7\%$ of respondents were missing some (but not all) survey items related to $c$ and $q$, respectively. As missing data is not the subject of this paper, we avoid this complication by simply using one imputed dataset, obtained through predictive mean matching vink2014predictive, and treating that as complete data. Specifically, the imputation is applied to the missing survey items prior to the calculation of $c$ and $q$. The imputation was done using the mice package in R van2011mice, with default settings for everything other than setting the imputation method to predictive mean matching. It follows that our analyses of these data are not directly comparable to the original one in wiertsema2019cost and we do not claim that ours are more valid in this regard. We discuss this issue further in Section (ref).
The statistical analysis in wiertsema2019cost is based on simple linear models\footnote{More specifically: the SUR model, as presented in Section (ref).}, without any interaction terms or other transformations of the covariates. This is imposes rather strong assumptions on the regression functions, which we find difficult to justify --- and as alluded to in the introduction, violations of these assumptions can lead to severely biased results. We hence found it necessary to develop the more flexible suBART model, which allows us to jointly estimate both treatment effects without strong functional assumptions.
We start by reviewing the original BART model in the univariate setting in Section (ref), in order to provide context for what is to follow. We then present our novel extensions to the multivariate continuous outcome setting in Section (ref) and the multivariate binary outcome setting in Section (ref). Specific details regarding posterior inference for suBART models are deferred to Section (ref).
BART was designed to solve the classic regression problem of the form $y_i = \mathbb{E}\left\lbrack y_i\:\vert\:\mathbf{x}_i\right\rbrack + \varepsilon_i$, where $y_i$ is a univariate response variable for observation $i=1,\ldots,n$, $\mathbf{x}_i$ is a $p-$dimensional predictor, and $\varepsilon_i \sim \text{N}(0,\sigma^2)$. The idea is to find a flexible approximation for the conditional expectation $\mathbb{E}\lbrack y_i\:\vert\:\mathbf{x}_i\rbrack$ by expressing it as a sum of regression trees. A regression tree consists of two components:
We now define a function $\mathbf{x}_i,\mathcal{T},\mathcal{M} \mapsto g(\mathbf{x}_i,\mathcal{T},\mathcal{M})$ as follows: $\mathbf{x}_i \in \mathcal{A}_\ell$ for exactly one $\ell$; then $g(\mathbf{x}_i,\mathcal{T},\mathcal{M}) \coloneqq \mu_\ell$. In the case of a single regression tree, we may then define the regression function $\mathbf{x}_i \mapsto \mathbb{E}\lbrack y_i\:\vert\:\mathbf{x}_i\rbrack$ by $\mathbb{E}\lbrack y_i\:\vert\:\mathbf{x}_i\rbrack \coloneqq g(\mathbf{x}_i,\mathcal{T},\mathcal{M})$. Figure (ref) shows an illustration of a simple regression tree, including the regression function it implies, for the case $p = 1$.
To extend this idea to additive regression trees, we consider not just one tree, but multiple trees $\mathcal{T}_1,\ldots,\mathcal{T}_m$, each with their own corresponding partitions and leaf nodes. Then we let $\mathbb{E}\lbrack y_i\:\vert\:\mathbf{x}_i\rbrack \coloneqq \sum^m_{t=1} g(\mathbf{x}_i,\mathcal{T}_t,\mathcal{M}_t)$. The definition of the statistical model now becomes \[y_i = \sum^m_{t=1} g\left(\mathbf{x}_i,\mathcal{T}_t,\mathcal{M}_t\right) + \varepsilon_i,\] which is the basic BART model presented in chipman2010bart. The model is typically not identified, since different sets of trees can lead to the same regression function. However, this is not a problem, since the individual trees are rarely of direct interest.
The sum of trees framework can also be used to model the conditional expectation of a binary response $\mathbf{y}$, which takes values in $\{0,1\}$. This is the probit BART model, again proposed originally in chipman2010bart. The model is easier to present and analyse when it is cast in terms of a continuous latent variable. Suppose $\mathbf{z}$ is such that
As before, we model $z_i$ as $z_i = \sum^m_{t=1} g\left(\mathbf{x}_i,\mathcal{T}_t,\mathcal{M}_t\right) + \varepsilon_i$, with $\varepsilon_i \sim \text{N}(0,1)$. This then implies
where $\Phi(\cdot)$ is the cumulative distribution function of the standard normal distribution.
For both types of outcome, given a sample of size $n$ and a univariate outcome vector $\mathbf{y} \in \mathbb{R}^{n}$ associated with a predictor matrix $\mathbf{X} \in \mathbb{R}^{n\times p}$, we are interested in sampling from the joint posterior distribution $\pi(\boldsymbol{\mathcal{T},\mathcal{M}} \:\vert\:\mathbf{y}, \mathbf{X})$ where $\boldsymbol{\mathcal{T}} = (\mathcal{T}_{1}, \ldots, \mathcal{T}_{m})$ denotes the collection of all trees and $\boldsymbol{\mathcal{{M}}} = (\mathcal{M}_{1}, \ldots, \mathcal{M}_{m})$ denotes their corresponding mean parameters. To obtain the posterior, it is necessary to define priors for both the trees and the terminal node parameters. Assuming the independence of leaf parameters conditional on the tree structures, chipman2010bart define the joint prior distribution as
To achieve conjugacy, it is typically assumed that the residual variance parameter $\sigma^{2}$ follows an inverse-gamma distribution and $\mu_{t\ell} \sim N(0,\sigma_{\mu}^{2})$, with $\sigma_{\mu}^{2} \coloneqq 1/(4\kappa^{2}m)$ being proportional to the number of trees in order to regularise the contribution of each tree. The definition of $\pi(T_{t})$ includes specifying the probability of a non-terminal node as $\alpha (1+\gamma_{t\ell})^{-\beta}$ where $\gamma_{t\ell}$ denotes the depth of that node. The hyperparameters $\alpha$ and $\beta$ take the default values suggested in chipman2010bart of $0.95$ and $2$, respectively, to favour shallow trees.
Once the prior is defined, a sampler for the conditional posterior of the trees and leaves, given $\sigma^2$, can be obtained. Referring to $\boldsymbol{\mathcal{T}}_{(-t)} \coloneqq \boldsymbol{\mathcal{T}} \setminus \{\mathcal{T}_{t}\}$ and $\boldsymbol{\mathcal{M}}_{(-t)} \coloneqq \boldsymbol{\mathcal{M}} \setminus \{\mathcal{M}_{t}\}$, an MCMC sampler is built by sequentially sampling from $\pi(\mathcal{T}_{t} \:\vert\:\boldsymbol{\mathcal{T}}_{(-t)},\boldsymbol{\mathcal{M}},\mathbf{y},\mathbf{X},\sigma^{2})$ and $\pi(\mathcal{M}_{t} \:\vert\:\boldsymbol{\mathcal{T}}, \boldsymbol{\mathcal{M}}_{(-t)}, \mathbf{y}, \mathbf{X}, \sigma^{2})$. It can be shown that $\mathcal{T}_{t}$ and $\mathcal{M}_{t}$ depend on $(\boldsymbol{\mathcal{T}}_{(-t)},\mathbf{y})$ only through the partial residuals $\mathbf{r}_{t} \coloneqq \mathbf{y} - \sum_{j\neq t}^{m}g(\mathbf{X},\mathcal{T}_{j},\mathcal{M}_{j})$. This can be used to construct a Bayesian back-fitting algorithm hastie2000bayesian. The successive draws become
where new trees and splitting rules are sampled through a Metropolis-Hastings step calculated using the integrated-likelihood for the tree $\mathcal{T}_{t}$ over the leaf parameters $\mathcal{M}_{t}$. See chipman2010bart, kapelner2016bartmachine, and tan2019bayesian for further details of the model and additional information the algorithmic implementation on which that of the suBART models is based.
We consider a regression problem of the form
where $\mathbf{y}^{(j)}$ represents the $j$-th component of a $d$-variate outcome, and $(\varepsilon_i^{(1)},\ldots,\varepsilon_i^{(d)})^\top \sim \text{MVN}_d(\mathbf{0}_d, \boldsymbol{\Sigma})$. In principle, different models can be used for each conditional expectation in Equation (ref). If the conditional expectations are all assumed to be linear in $\mathbf{x}_i$, we obtain the classic SUR model zellner1962efficient\footnote{A more modern presentation from a Bayesian perspective, closer in spirit to the present work, may be found in greenberg2012introduction, Chapter 9.} . We instead want to allow the possibility that the conditional expectations are non-linear. Given that the BART model has been shown to be a viable model in the one-dimensional setting, it seems reasonable to expect this viability to extend to the multivariate case. We therefore proceed by assigning an ensemble of regression trees to each $\mathbb{E}\lbrack y_i^{(j)}\:\vert\:\mathbf{x}_i\rbrack$ as follows:
where
The model setup is straightforward, and indeed similar to the original BART model in chipman2010bart. Our model is comprised of $d$ univariate BART models, which are linked through correlated error terms $\boldsymbol{\varepsilon}^{(j)}$. Although allowing an outcome-specific number of trees $m_j$ is possible, we assume that the number of trees $m$ is common across all $d$ outcomes for simplicity. However, we stress that each outcome is assigned its own separate collection of $m$ trees. Thus, there are $m\times d$ trees in total. This is a key difference compared to other multivariate BART models mcjames2024bayesian, um2023bayesian, which typically assume that all outcomes share the same $m$ trees and that the leaf parameters associated with each tree are vectors $\boldsymbol{\mu}_{t\ell}$ which follow a $d$-dimensional multivariate normal distribution. Although using $m \times d$ trees in suBART may be computationally prohibitive when $d$ is large, this is rarely the case in causal inference settings, and indeed we have only $d=2$ in the TTCM application.
A consequence of this assumption in existing multivariate BART versions is the implication that the $p$-dimensional predictor $\mathbf{x}_i$ is the same for all $d$ outcomes. However, in systems of equations such as (ref) and (ref), the set of predictors $\mathbf{x}_i$ need not be of common dimension $p$ for each outcome. Indeed, our subart software implementation allows different subsets of the predictors to be used for each constituent univariate BART model. Though we will henceforth assume that the predictors are the same for all outcomes, for simplicity, it is important to note that having the same set of predictors $\mathbf{x}_i$ available as candidates to form splitting rules in each set of trees does not imply that the trees for different outcomes will split on the same components of $\mathbf{x}_i$ at the same cutoff values. Unlike the standard SUR, imposing restrictions on the $\mathbf{x}_i$ for different outcomes is not required to ensure that different outcomes depend on different covariates under suBART, as different sets of trees will tend to form different partitions on different subsets of the feature space $\mathbf{X}$ anyway, by the inherent nature of the tree-generating process in BART. Although practitioners can impose such restrictions nonetheless, to strictly guarantee that the trees for a given outcome do not depend on certain predictors, this is an appealing property, in the sense that suBART minimises the need to pre-specify the parametric form of the models for each conditional expectation.
We assume that the $m$ regression trees for each outcome are all independent of each other and of the covariance matrix a priori, as per chipman2010bart; i.e,
We also assume that the leaves of a tree are conditionally independent given the tree structure, i.e.,
With this setup, prior distributions for $\mathcal{T}^{(j)}_{t}, \mu^{(j)}_{t\ell}$, and $\boldsymbol{\Sigma}$ are sufficient to specify the joint prior distribution of all model parameters. The tree structure $\mathcal{T}^{(j)}_{t}$ is assigned the same prior as in the original work by chipman2010bart, with default hyperparameters $\alpha = 0.95$ and $\beta = 2$ to favour shallow trees and avoid overfitting.
The shrinkage prior used for the leaf node parameters is also aligned with the approach of standard BART. As noted earlier, this prior is formulated conditionally on the tree $\mathcal{T}^{(j)}_t$. We assume that each outcome component $\mathbf{y}^{(j)}$ is re-scaled such that $y^{(j)}_i \in \lbrack -0.5, 0.5\rbrack$, using its respective observed maximum and minimum values, as outlined by chipman2010bart. This enables us to specify the prior probability that $\mathbb{E}\lbrack y^{(j)}_i\:\vert\:\mathbf{x}_i\rbrack$ lies within the rescaled interval. Consequently, the prior on the leaf node parameters is then
where $\sigma^{(j)}_{\mu} = 1/(2\kappa\sqrt{m})$. We suggest $\kappa = 2$ as a default choice, which assigns a prior probability of $0.95$ to the event $\{\mathbb{E}\lbrack y_i^{(j)}\:\vert\:\mathbf{x}_i\rbrack\in \lbrack-0.5, 0.5\rbrack\}$. As per the standard BART, having the prior variance of the leaf nodes be inversely proportional to $m$ helps to prevent overfitting. Crucially, the regularisation of each set of $m$ trees per outcome does not depend on the total number of outcomes $d$. Hence, $d$ being large does not incur any additional risk of overfitting.
In the univariate BART model, the error variance is typically assigned an inverse-gamma prior, which is conditionally conjugate and can thus be easily incorporated into a Gibbs sampler. Furthermore, by choosing the hyperparameters accordingly, it is possible to put an informative prior on the error variance. The recommended approach of chipman2010bart involves obtaining a `data-based overestimate' $\hat{\sigma}^2$ of the error variance $\sigma^2$ by fitting a linear model via least squares and calculating the standard deviation of the sample residuals\footnote{If the number of observations is smaller than the dimension of $\mathbf{X}$, we instead estimate $\hat{\sigma}$ using the LASSO in the glmnet R package glmnetR.}. Presumably, the true value of $\sigma^2$ is smaller than $\hat{\sigma}^2$, since more of the variation in $\mathbf{y}$ can be explained by the flexible BART model than the linear model. Therefore, a large prior probability $\alpha_\sigma$ is assigned to the event $\{\sigma^2 < \hat{\sigma}^2\}$ (for example, $\alpha_\sigma=0.95$).
We now wish to generalise this idea to the multivariate setting: for all components $j$ of the outcome vector $\left(\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(d)}\right)$, we obtain an overestimate $\hat{\sigma}^2_{j}$ of $\sigma^{2}_{j}$ in a similar fashion and want to assign the large prior probability $\alpha_\sigma$ to the event $\{\sigma^2_j < \hat{\sigma}^2_j\}$. Additionally, we also would like to control the prior on the correlations $\rho_{jk}$, where $j \neq k$. Since there is usually not strong prior information about these correlations, it is preferable for the prior to not be too informative in any direction. Finally, the chosen solution should be computationally tractable and readily incorporated into MCMC samplers. One such candidate is the inverse-Wishart distribution, which is often used as a prior for covariance matrices. As it is a straightforward multivariate generalisation of the inverse-gamma distribution, it facilitates easy computations. However, although it has the desirable property that the prior on $\Sigma$ is invariant to the ordering of the outcome vector $\mathbf{y}_i$, the inverse-Wishart prior is unfortunately not well-suited to meet our other aforementioned goals; among other problems, it imposes a strong prior dependency between the variances and the correlations and, when the true variances are small, overestimates variances and induces downward bias towards zero on the correlations alvarez2014bayesian.
In order to have more control over the variance parameters, we thus instead adapt an approach by huang2013simple and parameterise the covariance matrix $\boldsymbol{\Sigma}$ in a hierarchical fashion as follows:
As this prior is based on the inverse-Wishart, it shares the property of invariance to the ordering of the outcome vector. A similar prior was proposed by giannone2015prior in the context of vector autoregressions. Under our prior, the implied prior distribution for the correlations is given by
which matches the prior implied under the standard inverse-Wishart prior barnard2000modeling. Crucially, it does not depend on $A_1,\ldots,A_d$. It is uniform if and only if $\nu = 2$. For higher values, the prior increasingly concentrates around zero. Moreover, $\nu=1$ leads to a strongly informative and rarely reasonable prior on the correlations, which favours extreme values near $\pm 1$. Thus, we consider $\nu = 2$ a reasonable default choice, since we usually do not have any strong prior information on the correlations. It is also worth recalling from Section (ref) that the response vector is typically bivariate in the CEA setting which motivated the development of suBART, such that there is only one such correlation parameter.
From the previous definitions, the induced prior for the standard deviations is $\sigma_{j} \sim \text{Half-$t$}(\nu, A_{j})$, which in turn implies a scaled $F$-distribution as the prior for $\sigma_j^2$. Conversely, an inverse-Wishart prior would imply that the variances have inverse-gamma distributions. As half-$t$ priors on standard deviations are recommended over inverse-gamma priors on variances in univariate settings by gelman2006prior and others, the construction due to huang2013simple is preferable to the inverse-Wishart prior in multivariate settings for similar reasons. Using this construction, we can specify an uninformative prior on each $\sigma_{j}$, by setting a very large value for the associated hyperparameter $A_j$. However, in the spirit of the standard BART approach, we prefer to specify informative, data-dependent priors, by choosing each $A_j$ in such a way that $\Pr(\sigma_j < \hat{\sigma}_j)=\alpha_\sigma$ is ensured a priori.
As well as ensuring that Equation (ref) is uniform, another appealing property of specifying $\nu=2$ is the heaviness of the tails of the induced half-$t$ prior on $\sigma_j$, which makes the prior more robust to the choice of the scale $A_j$ and thereby reduces sensitivity to the data-based overestimate $\hat{\sigma}_j$ used for prior calibration. Although gelman2006prior highlights the special case of the half-Cauchy prior with $\nu=1$ and even heavier tails, we encountered numerical difficulties in our simulation experiments and the TTCM application with $\nu=1$, particularly in situations where the sample size is small, the dimension of $\mathbf{X}$ is large, and the variability of the multivariate outcome is almost entirely explained by $\mathbf{X}$.
We calibrate the prior for each $\sigma_j$ by setting up the following equation for each outcome
and solving it for $A_{j}$, which can be done through numerical root-finding. This expression is the cumulative distribution function of a $\text{Half-$t$}$-distributed random variable with degrees of freedom $\nu$ and scale parameter $A_j$, evaluated at $\hat{\sigma}_j$. By rewriting Equation (ref) in terms of regularised incomplete beta functions, it can be shown that this expression is continuous as well as strictly decreasing in $A_j$ and approaches $1$ and $0$ as $A_j$ approaches $0$ and $\infty$, respectively. Thus, the solution for $A_j$ exists and is unique.
While the previous model was designed to jointly model multiple continuous outcome variables, we now turn our attention to binary outcomes. We will present a generalisation of the linear multivariate probit model chib1998analysis, where the linear predictors are replaced by sums of regression trees. Alternatively, it can also be seen as a multivariate generalisation of the probit BART model.
Suppose that we have observed some predictor variables $\mathbf{x}_i$ and an associated collection of binary outcome vectors $\mathbf{y}_i \in \{0,1\}^d$ whose dependence on $\mathbf{x}_i$ we want to model. In particular, we do not want to assume that the components of $\mathbf{y}_i$ are conditionally independent, given $\mathbf{x}_i$; there may be some leftover correlation which is not explained by $\mathbf{x}_i$. As in the basic probit BART model, we cast the multivariate version in terms of latent variables $\mathbf{z}_i = (z_i^{(1)},\ldots,z_i^{(d)})$; the construction is exactly as per Equation (ref) for each outcome $j$. Then, the probit suBART model is given by
where
Conditional on the latent variables $\mathbf{z}_i$, the model is essentially the same as the suBART model presented earlier. The only difference is that for each error term $\boldsymbol{\varepsilon}^{(j)}$, we fix the variance at $1$. Without doing this, the model would be unidentified. This issue is not specific to probit suBART, but also arises in the linear multivariate probit model. See chib1998analysis for details. It follows that $\boldsymbol{\Sigma}$, the covariance matrix of the error terms, is equal to $1$ in each diagonal entry, and hence must be a correlation matrix.
The dependence structure of the trees, leaves, and covariance matrix is the same as for the suBART model. The priors for the trees are also exactly the same. Although the other priors are broadly similar, there are some differences that are worth briefly highlighting. The prior for the terminal node parameters $\mu_{t\ell}^{(j)}$ again follows Equation (ref), though the calibration of the variance hyperparameter requires more care. Taking ${\sigma^{(j)}_{\mu}} = q_z/(\kappa\sqrt{m})$, with $\kappa = 2$ as a default choice, we assign a prior probability of $0.95$ to the event $\{\mathbb{E}\lbrack z_i^{(j)}\:\vert\:\mathbf{x}_i\rbrack \in \lbrack -q_z,q_z\rbrack\}$. On the probability scale, this means that $\{\Pr(y_i^{(j)} = 1 \:\vert\:\mathbf{x}_i) \in \lbrack\Phi(-q_z), \Phi(q_z)\rbrack\}$. For example, when taking $q_z = 3$ as per chipman2010bart, we assign a prior probability of $0.95$ to the event $\{\Pr(y_i^{(j)} = 1\:\vert\:\mathbf{x}_i) \in \lbrack0.0013,0.9987\rbrack\}$. This is reasonable for many applications, since extremely small or large probabilities are uncommon.
In the probit setting, $\boldsymbol{\Sigma}$ is a correlation matrix and hence must be positive definite and have all diagonal entries equal to $1$. These restrictions make it difficult to choose a prior for $\boldsymbol{\Sigma}$ which has desirable properties and facilitates easy sampling. chib1998analysis present a prior (and related sampling strategy) which we found to be extremely inefficient in our application. We thus instead adapt an approach by zhang2020parameter (see also barnard2000modeling, where most of the following results originate, and zhang2006sampling). We introduce an auxiliary parameter $\mathbf{D}$, which is a $d \times d$ diagonal matrix. We then define $\mathbf{W} \coloneqq \mathbf{D}^{\frac{1}{2}} \boldsymbol{\Sigma} \mathbf{D}^{\frac{1}{2}}$, and assume the prior $\mathbf{W} \sim \text{Inv-Wishart}_d (\nu + d - 1, \mathbf{I}_d)$, where $\mathbf{I}_d$ is the $d$-dimensional identity matrix. As before, the induced prior on $\boldsymbol{\Sigma}$ is invariant to the ordering of the outcome vector.
Given $\mathbf{W}$, we can recover $\mathbf{D}$ and $\boldsymbol{\Sigma}$ thanks to the identities $\mathbf{D} = \text{diag}(\mathbf{W})$ and $\boldsymbol{\Sigma}= \mathbf{D}^{-\frac{1}{2}} \mathbf{W} \mathbf{D}^{-\frac{1}{2}}$. The induced marginal prior density of $\boldsymbol{\Sigma}$ is
$$\pi\left(\boldsymbol{\Sigma}\right) \propto \left(\det \boldsymbol{\Sigma}\right)^{\frac{1}{2}(\nu+d-1)(d-1)-1} \left(\prod_{j=1}^{d} \det \left\lbrack\boldsymbol{\Sigma}\right\rbrack_{jj} \right)^{-\frac{\nu + d - 1}{2}},$$ where $\lbrack\boldsymbol{\Sigma}\rbrack_{jj}$ is the $j$-th principle submatrix of $\boldsymbol{\Sigma}$. It can be shown that the marginal prior density for the correlations is again given by Equation (ref), as per the suBART model for continuous outcomes above, despite the different priors assumed for $\boldsymbol{\Sigma}$. The hyperparameter $\nu$ plays a similar role as before and we again consider $\nu=2$ to be a reasonable default choice.
This section describes strategies and algorithmic details for conducting posterior inference under suBART and probit suBART for multivariate continuous and multivariate binary outcomes, respectively. For both frameworks, sampling is performed using a Metropolis-within-Gibbs sampler based on their respective priors and model specifications. Given an observed sample $\mathbf{y}_{1}, \ldots, \mathbf{y}_n$, with $\mathbf{y}_{i} \in \mathbb{R}^{d}$, all computations are carried out conditionally on the covariates $\mathbf{X} \in \mathbb{R}^{n\times p}$. Thus, for the sake of readability, we will consider $\mathbf{X}$ fixed and not condition on it explicitly. Additionally, the expression $\pi(\theta\:\vert\:\boldsymbol{\Theta}\setminus\{\ldots\})$ denotes the conditional distribution of $\theta$ with respect to all parameters except $\theta$ itself and the ones listed within the braces. For example, $\pi(\mathcal{T}^{(j)}_t\:\vert\:\boldsymbol{\Theta}\setminus\{\mathcal{M}^{(j)}_t\})$ refers to the distribution of tree $\mathcal{T}^{(j)}_t$ given all parameters except $\mathcal{T}^{(j)}_t$ and $\mathcal{M}^{(j)}_t$. This slightly unusual notation reduces the complexity of the expressions which follow. Given that we routinely condition on all but two parameters, we find it clearer to highlight what is not being conditioned on.
For brevity, we define the estimate for $y_{i}^{(j)}$ as $\hat{y}^{(j)}_{i} \coloneqq \linebreak{}\sum_{t=1}^{m}g(\mathbf{x}_i, \mathcal{T}^{(j)}_{t}, \mathcal{M}^{(j)}_{t})$ and $\mathbf{\hat{y}}_{i}=(\hat{y}_{i}^{(1)},\ldots,\hat{y}^{(d)}_{i})^{\top}$, since there are now multiple components of $\mathbf{y}_{i}$. Analogously, we define the residuals from a tree $t$ associated with the $j$-th component as $\mathbf{r}^{(j)}_t \coloneqq \{r^{(j)}_{t1}, \ldots, r^{(j)}_{tn}\}$, where $r^{(j)}_{ti} \coloneqq y^{(j)}_{i} - \sum_{k \neq t}^{m}g(\mathbf{x}_{i}, \mathcal{T}^{(j)}_{k}, \mathcal{M}^{(j)}_{k})$. As per the standard BART, we are interested in sampling from the posterior distribution $$\pi(\boldsymbol{\Theta}\:\vert\:\mathbf{Y}) = \pi\left(\left(\boldsymbol{\mathcal{T}}^{(1)},\boldsymbol{\mathcal{M}}^{(1)}\right),\ldots,\left(\boldsymbol{\mathcal{T}}^{(d)},\boldsymbol{\mathcal{M}}^{(d)}\right),\boldsymbol{\Sigma}, a_1,\ldots,a_d\:\vert\:\mathbf{Y}\right),$$ which, due to the back-fitting algorithm hastie2000bayesian and properties of the multivariate normal distribution, can be obtained from sequential draws from a collection of conditional distributions. In the multivariate continuous outcomes setting, we have that $\mathbf{y}_i \sim \text{MVN}_d \left(\hat{\mathbf{y}}_i, \boldsymbol{\Sigma}\right) $. For the following, we will also need the conditional distribution of any component $y^{(j)}_i$ given all other components $\mathbf{y}_{i}^{(-j)}$. Using a well-known result (see e.g., bierens2004introduction, Section 5.3), this can be found in closed form:
where $\boldsymbol{\Sigma}_{(-j)(-j)}$ is the submatrix obtained by excluding the $j$-th row and column from $\boldsymbol{\Sigma}$. Analogously, $\boldsymbol{\Sigma}_{j(-j)}$ is the vector obtained by selecting the $j$-th row and excluding the $j$-th column. Using the conditional distribution in Equation (ref), Metropolis-Hastings draws from the posterior distribution of trees $\pi(\mathcal{T}^{(j)}_{t}\:\vert\:\mathbf{r}_{t}^{(j)}, \boldsymbol{\Theta}\setminus\{\mathcal{M}_{t}^{(j)}\})$ can be obtained in essentially the same way as for the univariate BART model; see Appendix A of kapelner2016bartmachine for details. Ultimately, as described in Section (ref), the sampler for the joint posterior distribution of the trees $\mathcal{T}^{(j)}_{t}$ and their parameters $\mathcal{M}_{t}^{(j)}$, for the $j$-th component of $\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(d)}$, is given by successively drawing from $\mathcal{T}^{(j)}_{t} \:\vert\:\mathbf{r}^{(j)}_{t},\boldsymbol{\Theta}\setminus\{\mathcal{M}_{t}^{(j)}\}$ and $\mathcal{M}^{(j)}_{t} \:\vert\:\mathbf{r}^{(j)}_{t},\boldsymbol{\Theta}$.
Notably, the algorithm reduces to that of the standard BART approach when $d=1$. Indeed, each draw above can also be viewed as univariate BART --- albeit with distinct mean and variance parameters --- since it is conditioned on the values of all other components in $\mathbf{Y}$, as illustrated by Equation (ref). The full structure of the suBART sampler is given in Algorithm (ref), but we first describe the remaining required posterior conditional distributions.
Firstly, the posterior distribution for $\mu^{(j)}_{t\ell}$ is given by
where $u_{i}^{(j)} \coloneqq \boldsymbol{\Sigma}_{j(-j)} \boldsymbol{\Sigma}_{(-j)(-j)}^{-1} (\mathbf{y}^{(-j)}_{i} - \hat{\mathbf{y}}_{i}^{(-j)})$, $v^{(j)}\coloneqq \Sigma_{jj} - \boldsymbol{\Sigma}_{j(-j)} \boldsymbol{\Sigma}_{(-j)(-j)}^{-1} \boldsymbol{\Sigma}_{(-j)j}$, and $\smash{n_{t\ell}^{(j)}}$ denotes the number of observations in the given terminal node. The derivation again closely follows that for the univariate BART model; see Appendix B.1 of tan2019bayesian for details. Indeed, when $d=1$, these definitions imply that the $u_i^{(j)}$ term vanishes and $v^{(j)}=\Sigma_{jj}=\sigma^{2}$, resulting in an expression identical to the original BART formulation. Secondly, the conditional posteriors of $\boldsymbol{\Sigma}$ and the auxiliary parameters $a_1,\ldots,a_d$ take simple forms, due to the conditional conjugacy in the construction of huang2013simple. We have
where $(\boldsymbol{\Sigma}^{-1})_{jj}$ denotes the $j$-th entry along the diagonal of $\boldsymbol{\Sigma}^{-1}$, and
where $\mathbf{S} = \sum_{i=1}^{n} (\mathbf{y}_i - \hat{\mathbf{y}}_i) (\mathbf{y}_i - \hat{\mathbf{y}}_i)^\top$ and $\mathbf{S}_0$ is as given in Equation (ref).
In the multivariate binary setting, the goal is to draw samples from the similar posterior $$\pi(\boldsymbol{\Theta}\:\vert\:\mathbf{Y}) = \pi\left(\left(\boldsymbol{\mathcal{T}}^{(1)},\boldsymbol{\mathcal{M}}^{(1)}\right),\ldots,\left(\boldsymbol{\mathcal{T}}^{(d)},\boldsymbol{\mathcal{M}}^{(d)}\right),\boldsymbol{\Sigma},\mathbf{D}\:\vert\:\mathbf{Y}\right).$$ Note that $\mathbf{W}$ is a deterministic function of $\boldsymbol{\Sigma}$ and $\mathbf{D}$, and is hence omitted from the above distribution. However, it will be more convenient to work with the joint posterior of the parameters and the latent variables $$\pi\left(\boldsymbol{\Theta},\mathbf{Z}\:\vert\:\mathbf{Y}\right) = \pi\left(\left(\boldsymbol{\mathcal{T}}^{(1)},\boldsymbol{\mathcal{M}}^{(1)}\right),\ldots,\left(\boldsymbol{\mathcal{T}}^{(d)},\boldsymbol{\mathcal{M}}^{(d)}\right),\boldsymbol{\Sigma},\mathbf{D},\mathbf{Z}\:\vert\:\mathbf{Y}\right),$$ for which the sampling algorithm is very similar to the previously presented Algorithm (ref) in the continuous setting. For brevity, we discuss only the required modifications to Algorithm (ref) without presenting a new algorithm in full.
The updates for the trees and leaf nodes stay essentially the same, with the one difference being that the latent variables $\mathbf{Z}$ replace the data $\mathbf{Y}$. An important additional step is that the latent variables for each component $j$ should be updated after line 10 in Algorithm (ref). In a similar manner to Equation (ref), the marginal distribution of $z_i^{(j)}$ can be obtained via
where we exploit the fact that $\Sigma_{jj}=1$, since $\boldsymbol{\Sigma}$ is a correlation matrix in the probit setting. However, sampling the latent variables also requires conditioning on $\mathbf{Y}$. In doing so, we find that $\pi(z^{(j)}_i \:\vert\:\mathbf{z}^{(-j)}_i, \mathbf{Y}, \boldsymbol{\Theta})$ follows a truncated normal distribution with the same location and scale parameters as Equation (ref). However, we encounter a case distinction for the support of this conditional posterior distribution based on the values of the associated response. If $y_i^{(j)}=0$, which implies that $z_i^{(j)} \le 0$, the support is truncated to $(-\infty, 0\rbrack$. Conversely, if $y_i^{(j)}=1$, which implies that $z_i^{(j)} > 0$, the support is truncated to $(0, \infty)$. In each case, we draw the sample through the method proposed by robert1995simulation.
The other major difference for the probit suBART sampler concerns the update of $\boldsymbol{\Sigma}$. Since the $a_j$ parameters only apply to $\boldsymbol{\Sigma}$ in the continuous setting, their updates are of course dropped. In the probit setting, accounting also for the auxiliary parameter $\mathbf{D}$, we can write the associated conditional posterior distribution as
where $\pi(\boldsymbol{\Sigma},\mathbf{D})=\pi(\mathbf{W})$ is the aforementioned inverse-Wishart prior on $\mathbf{W}$. The term $(\det \mathbf{D})^{\frac{d - 1}{2}}$ is the Jacobian determinant which arises due to the change of variables $\mathbf{W} \mapsto \boldsymbol{\Sigma},\mathbf{D}$. This distribution is not of known form and hence cannot be sampled from directly. We instead proceed using the parameter-expanded Metropolis-Hastings (PX-MH) algorithm of zhang2020parameter. This defines a proposal for $\mathbf{W}^{(k+1)} \:\vert\:\mathbf{W}^{(k)} \sim \text{Inv-Wishart}_{d}\left(\nu_{\text{prop}}, \nu_{\text{prop}} \times \mathbf{W}^{(k)} \right)$ where $k$ denotes the current MCMC iteration and $\nu_{\text{prop}}$ is a tuning parameter which can be fine-tuned to balance the trade-off between the the acceptance rate of the PX-MH algorithm and the autocorrelation in the resulting chains. Such tuning should account for the sample size and the dimension $d$, as outlined in prior studies zhang2006sampling. We set $\nu_{\text{prop}}$ so to as to ensure an acceptance rate of 20-30%.
In this section, we evaluate the predictive performance of the proposed models through experiments with simulated data. Section (ref) and Section (ref) are devoted to simulation designs in which the responses $\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(d)}$ are all continuous and each outcome $\mathbf{y}^{(j)} \in \{0,1\}$ is binary, respectively. We undertake a comparative analysis of suBART against benchmark models, including the standard BART model applied independently to each response, the multivariate BART (mvBART) model, and a Bayesian linear seemingly unrelated regression (SUR) model. We also consider probit versions of each model, where available. This comparative study aims to explore various aspects of the models, including predictive performance and their ability to accommodate assumptions regarding correlation among responses and/or assumptions of linearity. Furthermore, our simulations aim to elucidate the primary distinctions between suBART and mvBART. For example, the splits generated by trees under the mvBART framework entail a splitting rule in all components of $\mathbf{Y}$, potentially deviating from an accurate representation of the true function $f(\mathbf{X})$ for some scenarios. Consequently, each response variable in our experiments is generated using a different subset of covariates. Additionally, a significant improvement in predictive performance capacity is anticipated when compared with the linear SUR model as, for the most part, the responses in our experiments are almost all assumed to be non-linearly related to the covariates.
These assumptions were tested over $100$ replications of each simulation scenario, using different sample sizes of $ n_{\text{train}} = n_{\text{test}} = \{250, 500, 1000 \}$ for training and test samples respectively. The metrics employed to assess differences in model performance included the root mean squared error (RMSE), the continuous ranked probability score raftery2007crps, and the prediction interval (PI) coverage for the multivariate regression cases, while the logarithmic loss, the accuracy (ACC), and the credible interval (CI) coverage of the probabilities from $\Phi(z_{i}^{(j)})$ were used for the multivariate probit scenarios. The posterior means for $\hat{y}_{i}^{(j)}$, $\sigma_{j}$, and $\rho_{jk}$ are obtained by averaging the corresponding posterior samples, while the 50% PIs in the case of continuous outcomes are estimated by drawing from the predictive distribution of $y_{i}^{(j)}$ and computing the 25-th and 75-th percentiles.
Our choice for the inverse-Wishart hyperparameter is $\nu = 2$, reflecting our lack of prior information about the correlation structure huang2013simple. For the probit version of suBART, our selection of the proposal degrees of freedom $\nu_{\text{prop}}$ varies according to the dimension of the response and the sample size. We adopt the default values of $\nu_{\text{prop}}=n_{\text{train}}/10$ when $d=2$ and $\nu_{\text{prop}}=n_{\text{train}}/2$ when $d=3$. These settings appear to ensure sufficiently well-behaved sampling in our simulation experiments. For the MCMC settings for all methods under our experiments with binary responses, we set a total of $N_{\text{MCMC}} = 10000$ iterations, of which $N_{\text{burn-in}} = 2000$ draws are discarded as burn-in. Otherwise, for our continuous response experiments, $N_{\text{MCMC}}=5000$ and $N_{\text{burn-in}}=1000$ suffice for satisfactory convergence for all methods. The higher number of iterations for the probit version of suBART is required since the sampler for $\boldsymbol{\Sigma}$ is less efficient in this case. See zhang2006sampling for details on these aspects of the sampler and Appendix (ref) for a discussion of the satisfactory convergence of the continuous and probit versions of suBART for these simulation experiments.
Lastly, we note the software implementations for each model included in the comparisons. The suBART models are fitted using our own subart implementation and the BART models are fitted using the dbarts package dorie2024package, while different software is used to fit the linear Bayesian SUR models (henceforth BayesSUR) depending on the nature of the outcome variables. We use the R package surbayes alt2020surbayes for the continuous response experiments in Section (ref), which is based on the direct Monte Carlo scheme and default priors from zellner2010direct. For the binary response experiments in Section (ref), we use the probabilistic programming language Stan stanusersguide, through the rstan package rstanref which provides an \textsf{R} interface for this library. The mvBART model was evaluated using the \texttt{skewBART} implementation provided by um2023bayesian, specifically by setting the argument \texttt{do_skew=FALSE} of the main \texttt{MultiskewBART()} function. None of our subsequent analyses use the sparsity-inducing Dirichlet prior on the splitting probabilities linero2018highdim which is available as an option in the \texttt{skewBART} package, as the TTCM application and our simulated scenarios each have only a moderate number of covariates. We have not yet incorporated this prior into our \texttt{subart} software implementation for the same reason. Furthermore, it is worth noting that the \texttt{skewBART} package is currently limited to continuous scenarios with two dimensions, thereby results were constrained to such cases. The default hyperparameter settings were retained for all competing models with the exception of the number of trees for the tree-based models. We use $m=100$ trees \emph{per outcome} for suBART and independent BART and --- for the sake of fairness --- configure mvBART with $m \times d$ trees.
The simulation scenario described by the system of equations below was created to accommodate different types of complexity. In this experiment, the values of the response are non-linear functions modified from examples given in friedman1991multivariate and breiman1996bagging for a multivariate response scenario. The third response is exceptional in the sense that the generating function is purely linear. Note that correlated noise $(\varepsilon_i^{(1)}, \ldots, \varepsilon_i^{(d)})^\top \sim \text{MVN}_{d}(\mathbf{0}_d, \boldsymbol{\Sigma})$ is subsequently added to the $d$-dimensional response.
\noindentFriedman \#1:
The $p=10$ predictors are generated from a uniform distribution. It is notable that not all predictors are used to build the responses. However, all predictors are used for model fitting. It is anticipated that the tree-based models will be able to identify the uninformative noise variables. It is also essential to emphasise that each component of the outcome vector $\mathbf{Y}$ is derived from a distinct set of predictors. However, no restrictions are imposed on which predictors are associated with which response during model fitting. As previously mentioned, it is anticipated that mvBART may encounter challenges in accurately approximating the true generating functions in such cases. For each tree, the partitioning of the covariate space is shared across all responses, which may be an unrealistically restrictive assumption.
We varied the dimension of the covariance error matrix within $d = \{2,3\}$ and defined $\boldsymbol{\Sigma}$ accordingly, with specific values assigned to each parameter $\sigma_{j}$ and each correlation parameter $\rho_{jk}$, for all $j\neq k$. The error covariance parameters were set as detailed in Table (ref) and it is the first two responses $\smash{y_i^{(1)}}$ and $\smash{y_i^{(2)}}$ which comprise the $d=2$ setting. The restriction to $d=2$ enables consideration of the mvBART model in the comparisons, owing to the aforementioned limitation of the skewBART software to bivariate outcome settings.
A comparison of results is depicted in the boxplots in Figure (ref) and Figure (ref) which confirm previous assumptions about suBART performance. These figures illustrate the results for Friedman \#1 with $d=2$ and $d=3$, respectively, with $n_{\text{train}} = n_{\text{test}} = 1000$ in each case. In general, suBART exhibits either slightly superior or competitive predictive performance compared to BART and mvBART, as evidenced by small average RMSE and CRPS values over the test samples. Furthermore, all tree-based methods exhibit a clear superiority over BayesSUR in estimating the non-linear responses. The primary discrepancy occurs when $j=3$, where BayesSUR has the best performance owing to the linearity of this response.
In terms of uncertainty estimation, Figure (ref) shows that all methods exhibit reasonable coverage rates for the prediction intervals when $d=2$, except for the first outcome where the coverage for both mvBART and BayesSUR is notably high. Figure (ref) also shows satisfactory coverage for both suBART and mvBART, while the coverage under BayesSUR is only close to the nominal rate when $j = 3$, which is again due to the linearity of this response.
For proper uncertainty quantification, it is essential to correctly estimate the elements of the covariance matrix $\boldsymbol{\Sigma}$. Table (ref) displays the RMSE and coverage rates of 50% credible intervals for the correlation parameters $\rho_{jk}$ for all $j \neq k$ provided by suBART, mvBART, and BayesSUR. Notably, correlation values for the BART model are not provided as it assumes independence (i.e., $\rho_{jk}= 0\:\forall\:j\ne k$), and the mvBART estimates are provided only when $d=2$, due to the aforementioned limitations of the skewBART package. Standard deviation parameters $\sigma_j$ are also given where applicable. From the results, it is clear that suBART outperforms BART, mvBART, and BayesSUR in terms of coverage, demonstrating its superior ability to estimate correlation structures. The coverage values of zero for BayesSUR are particularly notable and suggest an inability to accurately estimate correlations when assuming linear regressions for non-linear responses. In terms of RMSE, suBART is superior to BART and BayesSUR and comparable to mvBART, albeit only in the $d=2$ setting.
Equivalent figures and tables summarising the results for the remaining scenarios, with $d=3$ and/or different sample sizes, yield the same conclusions as above and have been deferred to Appendix (ref) for brevity.
The experiments for binary responses are aligned with those from Section (ref), wherein different sample sizes of $ n_{\text{train}} = n_{\text{test}} = \{250, 500, 1000 \}$ are used for training and test samples respectively. The simulation of the latent variables $z^{(j)}$ is described by the system of equations below. Other than when $j=3$, the values of the latent variables are non-linear functions. Recall that correlated noise $(\varepsilon_i^{(1)}, \ldots, \varepsilon_i^{(d)})^\top \sim \text{MVN}_{d}(\mathbf{0}_d, \boldsymbol{\Sigma})$ is subsequently added to the $d$-dimensional latent variable and that the generating process for each binary response follows Equation (ref) thereafter. As the generative model is invariant to the values of the variances chib1998analysis, we set them to $1$ for convenience, but the true correlations are set to the same values as in Table (ref).
\noindentFriedman \#2:
The models evaluated for these experiments are the probit suBART and probit extensions of the standard BART and BayesSUR. The mvBART model is excluded, since the skewBART software does not offer a probit implementation of the model, for any dimensionality. We persist in evaluating settings with varying dimension $d=\{2,3\}$ with the $d=2$ setting again comprising the first two responses. The results for $d=3$, with $n_{\text{train}}=n_{\text{test}}=1000$, are summarised in Figure (ref) and are consistent with the findings from Section (ref). When logarithmic loss and ACC are considered as metrics for evaluating predictive performance, suBART either exhibits superior results or comparable averages. Both tree-based models outperform Bayesian SUR, with the exception of the linear third response in the $d=3$ setting, as per the continuous-outcome simulation experiments in Section (ref).
Regarding calibration, it is evident that suBART and BART are satisfactory across all scenarios, showing coverage rates close to the nominal value. On the other hand, due to the inherent linearity of BayesSUR, its calibration performance is notably poorer, though the third linear response in Figure (ref) is again an exception in this regard, as per Section (ref). As per the Friedman \#1 simulations in Section (ref), the results for remaining sample sizes are deferred to Appendix (ref), as they lead to similar conclusions.
The results regarding the estimation of the correlation parameters are presented in Table (ref), which includes the RMSE and CI coverage for the correlation parameters $\rho_{jk}$ associated with binary responses when $n_{\text{train}}=1000$. As BART assumes independent errors, it is excluded from these results. Overall, these results are consistent with Table (ref) in clearly demonstrating superior performance for suBART in terms of estimating correlation structures.
We now apply the continuous suBART model to the TTCM data from wiertsema2019cost described previously in Section (ref). We do this in two different ways. We first construct a simulation experiment in Section (ref) using the baseline covariates contained in the TTCM data. We then apply six different models (detailed below) to 1000 datasets created through simulation. As the true treatment effects are known, this will allow us to evaluate how well each model estimates them. Subsequently, we apply the same six models to the full dataset from the TTCM study in Section (ref), and compare the results.
In each case, we compare three methods for estimating the treatment effects: (1) suBART, (2) mvBART (using the aforementioned R package skewBART), and (3) BayesSUR (using the aforementioned R package surbayes). As before, we use $m=100$ trees for each outcome under suBART and configure mvBART with $m\times d$ trees. Additionally, we do not impose any restrictions on the sets of covariates associated with each response, for any of these methods. This means that all sets of trees for suBART and the single set of multivariate trees for mvBART are \emph{allowed} to form splitting rules using all covariates in Table (ref) and all linear regressions for BayesSUR also share all covariates. Following wiertsema2019cost, we do not specify any interaction effects or non-linear terms in the linear predictors for BayesSUR, owing to the difficulty of pre-specifying appropriate functional forms in the presence of a large amount of covariates and the associated challenges in terms of model selection. In any case, dorie2019automated found that linear models perform poorly in causal settings even when also including interactions, polynomial terms, and regularisation to avoid overfitting. Conversely, BART-based methods are well-equipped to automatically capture low-order interactions and non-linearities rockova2020semi,rockova2020posterior.
In addition, we evaluate each method again with the set of predictors augmented using point estimates of the propensity scores obtained via probit BART. This procedure is inspired by the univariate ps-BART method proposed by hahn2020bayesian, which induces a covariate-dependent prior on the regression function. This step can substantially reduce bias due to regularisation-induced confounding, thereby making the models less prone to attributing the effect of confounders to the treatment variable, and yield results which are more robust to model misspecification li2023bayesian\footnote{There is some conceptual disagreement as to whether one should merely include a point estimate of the propensity score among the covariates, or condition on the propensity score; the latter implies a joint model for the treatment and outcomes. hahn2020bayesian argue convincingly in favour of the former; see in particular Section 8.3. For an opposing viewpoint arguing for joint modelling, see for example li2023bayesian, Section 5. Our adoption of the point estimate approach is also well-supported by extensive simulations conducted by dorie2019automated which show that ps-BART performs well from a frequentist perspective.}. We first estimate propensity scores for each patient through probit BART, which give the probability that a patient received treatment $1$, conditional on their baseline characteristics, and then add the posterior mean propensity score estimates to the set of predictors $\mathbf{x}_i$ used to estimate the conditional expectations $\mathbb{E}\lbrack c\:\vert\:t, \mathbf{x}_i\rbrack$ and $\mathbb{E}\lbrack q\:\vert\:t, \mathbf{x}_i\rbrack$ via the chosen model. Thus, we expand the comparison to include what we refer to as ps-suBART, ps-mvBART, and ps-BayesSUR, which are straightforward adaptations of the univariate ps-BART method. Although the use of BART-based propensity scores in conjunction with linear SUR deviates from usual practice in the applied health economics literature, we nonetheless use the same set of propensity scores estimated via probit BART for each method, in order to ensure the comparison is fair in this regard. Expanding the comparison to include versions of each method with and without propensity scores will help to establish the extent to which differences in results are attributable to differences in model specification or due to the inclusion of propensity scores.
Given that there are several unordered categorical covariates with more than two levels in the TTCM dataset --- Table (ref) reports Medical history, Trauma type, and Fracture region --- such variables require particular attention. With probit BART (when estimating the propensity scores) and suBART, we treat such predictors by adapting the method of ordering categories proposed by breiman1984classification. More specifically, whenever such a variable is split on, we order the levels according to the mean values of the partial residuals and then treat the reordered variable in the same way as a continuous variable. This approach makes the tree-growing process more effective by allowing splitting rules which group multiple levels and can reduce the risk of overfitting by allowing for shallower trees. By contrast, both mvBART and BayesSUR accommodate these predictors as dummy variables.
In this section, we present the results of a simulation experiment based on the TTCM data. We keep the baseline covariates fixed and define generative functions for the potential outcomes and propensity scores. We then simulate treatments and outcomes for each patient and use all six previously discussed models to estimate the treatment effects and $\text{INB}_\lambda$ at representative values of $\lambda=20\mathrm{k}$ and $\lambda=50\mathrm{k}$, where $\lambda=20\mathrm{k}$ corresponds to €20,000. By repeating this procedure $1000$ times, we get an idea of how well each method performs at recovering the true parameters used in the simulation. We summarise the performance in terms of bias, standard deviation (SD), and RMSE, all evaluated using the posterior mean as a point estimate. We further include the coverage rates of 50% credible intervals (CI coverage), and the average width of 50% credible intervals (CI width). We stress that the chosen $\lambda$ values are not indicative of the cost of the TTCM intervention, which averages €272 per patient wiertsema2019cost. Rather, $\lambda=20\mathrm{k}$, for example, is indicative of a situation in which decision-makers would be willing to pay €20,000 per additional unit of healthcare-related quality of life, which is much greater as the treatment effect of $\Delta_q$ is near-zero. In any case, both €20,000 and €50,000 are commonly used thresholds for determining whether an intervention represents value for money in the CEA literature drummond2015methods,gabrio2019bayesian.
We first present the generative models used. The ideas are inspired by similar simulation experiments in the literature --- hill2011bayesian and hahn2020bayesian in particular --- but there are also important novelties. There is first the setting with two correlated outcomes, which is characteristic of the cost-effectiveness context for which the suBART model was developed. We also use somewhat larger variances for the error terms in the simulation (not necessarily in absolute terms, but relative to the variation in $\mathbb{E}\lbrack y_i\:\vert\:\mathbf{x}_i\rbrack$). In our experience, this is more reflective of real cost-effectiveness datasets.
Our simulation includes all covariates from Table (ref). We refer to the covariates in the same order as listed in the table (so $ x^{(1)}_i$ refers to age, $x^{(2)}_i$ refers to gender, and so on). Here, all continuous covariates are normalised prior to generating the outcomes. For the sake of clarity, we introduce the additional notation $\mu_c(\mathbf{x}_i) \coloneqq \mathbb{E}\left\lbrack c_i(0)\:\vert\:\mathbf{x}_i\right\rbrack$ and $\tau_c(\mathbf{x}_i) \coloneqq \mathbb{E}\left\lbrack c_i(1)\:\vert\:\mathbf{x}_i\right\rbrack - \mathbb{E}\left\lbrack c_i(0)\:\vert\:\mathbf{x}_i\right\rbrack$, adopt analogous definitions for $q$, and further let $\overline{\mu_q(\mathbf{x}_i)}$ and $\text{sd}(\mu_q(\mathbf{x}_i))$ denote the corresponding sample mean and standard deviation, respectively.
\noindentData-generating process for outcomes:
$$
\sim MVN_2 \left(
,
\right).$$ The generating functions for $c_i$ and $q_i$ each involve subsets of the available covariates which only partially overlap. The coefficients are chosen such that the simulated outcomes are reflective of values typically observed in CEA data. We emphasise that the data-generating process is linear for $c_i$ and non-linear for $q_i$. In particular, the treatment effect $\tau_q(\mathbf{x}_i)$ is heterogeneous, given its dependence on $\smash{x_i^{(11)}}$, while $\tau_c(\mathbf{x}_i)$ is constant. These settings imply true values for $\Delta_c, \Delta_q$, $INB_{20\mathrm{k}}$, and $INB_{50\mathrm{k}}$ of $500$, $0.0399$, $299.24$, and $1489.10$, respectively. We specify a negative correlation of $-0.25$ for the errors to model a common occurrence in CEA studies: patients who report lower well-being also have higher associated costs \citep{glick2014economic}. The results of additional experiments using correlation values of $0$ and $-0.5$ --- considering only ps-suBART and independent BART models --- are deferred to Appendix (ref).
We further define the propensity scores $\pi_i \coloneqq \Pr(t_i = 1 \:\vert\:\mathbf{x}_i)$ as follows:
By design, the propensity scores are strongly dependent on $\mu_q(\mathbf{x}_i)$, which is the expected health outcome under treatment $0$. This is an example of targeted selection hahn2020bayesian, which induces strong confounding with respect to the outcome $q_i$, and reflects how treatments are assigned in reality: patients expected to have a bad outcome under treatment $0$ are more likely to be given treatment $1$ by medical professionals. Further confounding with respect to $c_i$ arises by including $\smash{x^{(10)}_i}$, on which $\mu_c(\mathbf{x}_i)$ also depends, in Equation (ref). The coefficients are chosen to ensure roughly balanced sample sizes for the two treatment groups.
The performance metrics for $\Delta_c$, $\Delta_q$, and $\text{INB}_\lambda$ (under all six models, across $1000$ simulation replications) are presented in Table (ref), Table (ref), and Table (ref), respectively. We see that including propensity scores generally leads to smaller bias and improved coverage rates, while slightly increasing the SD of the estimates. Including propensity scores also lowers the RMSE in all but one case while increasing the average CI width for all methods. Across all estimands, the intervals from BayesSUR are wider than those from mvBART, which in turn are wider than those from suBART, with or without propensity scores. For $\Delta_c$, the BayesSUR models have the smallest bias, which is expected given the linearity of $\mu_c(\mathbf{x}_i)$ in the data-generating process. For $\Delta_q$, ps-suBART exhibits the lowest bias by far, and also the lowest RMSE by a reasonable margin. BayesSUR exhibits substantial bias for $\Delta_q$, although this is mitigated through the inclusion of propensity scores. For $\text{INB}_\lambda$, ps-suBART's superiority in terms of bias, RMSE, and coverage notably holds at both $\lambda=20\mathrm{k}$ and $\lambda=50\mathrm{k}$.
The results also show that the non-linear BART models perform only slightly worse than linear BayesSUR when the conditional expectation is actually linear ($\Delta_c$ in Table (ref)), and much better when it is non-linear ($\Delta_q$ in Table (ref)). The subpar performance of SUR in the second case can be improved significantly by the inclusion of propensity scores, which seem to absorb much of the non-linearity, although we emphasise that this is a somewhat unnatural model which is unlikely to be used by anyone in practice; if we estimate propensity scores through a non-linear model, it would seem appropriate to also use a non-linear outcome model. Comparing ps-suBART and ps-mvBART specifically, we see a better performance from ps-suBART across all estimands and metrics. Overall, ps-suBART exhibits the smallest bias and RMSE, while also being competitive in terms of interval width and coverage.
Figure (ref) shows the distribution of the estimated propensity scores. Although there is some overlap, it is evident that the treatment groups are quite imbalanced with respect to their baseline characteristics and that some form of covariate adjustment is necessary to avoid biased results. As already described, we reevaluate each model with estimated propensity scores included as an additional covariate. In Figure (ref), we show highest density regions of kernel density estimates of the posterior distributions of $\Delta_c$ and $\Delta_q$ --- obtained through suBART, mvBART, and BayesSUR, as well as their counterpart models which also incorporate the estimated propensity scores as additional predictors --- in the form of a cost-effectiveness plane (CEP; see gabrio2019bayesian for more details). In general, for both $\Delta_c$ and $\Delta_q$, the centers of the distributions are shifted further away from $0$ and posterior uncertainty is greater when the propensity scores are incorporated. These differences are most pronounced between ps-BayesSUR and BayesSUR. A notable distinction between suBART and BayesSUR is that the latter appears to be more uncertain about $\Delta_q$. The same applies to the comparison between ps-suBART and ps-BayesSUR.
In Table (ref), we provide summary statistics for $\Delta_c$, $\Delta_q$, and the INB at the same representative values of $\lambda = 20\mathrm{k}$ and $\lambda=50\mathrm{k}$ adopted in Section (ref). We see some considerable differences between ps-suBART and the other methods. This suggests that there may be strong non-linear functional relationships between the covariates and the outcomes (such that suBART benefits from relaxing the linearity assumption of BayesSUR, without requiring pre-specification of the functional forms) and that those relationships may differ for the two outcomes $c_i$ and $q_i$ (such that suBART benefits from relaxing the mvBART assumption of a common tree structure). With the exception of ps-BayesSUR, the INB values from all models are much lower than those from ps-suBART. For all models, the estimated treatment effects are consistently larger in absolute value for each method when propensity scores are included. The same is true for $\text{INB}_\lambda$ at both $\lambda$ values. Surprisingly, for both $\Delta_q$ and $\text{INB}_\lambda$, ps-BayesSUR produces far higher estimates than all other models, including ps-suBART. We recall in this context that in our simulation study, the standard error of the $\Delta_q$ estimates was the highest by far for ps-BayesSUR. Hence this estimate tends to be rather volatile, and we are inclined to put more faith in the estimates from the other models.
Rather than relying on a single $\lambda$ value, we also show the probability of cost-effectiveness as a function of the willingness-to-pay $\lambda$ in Figure (ref). This plot, called a cost-effectiveness acceptability curve lothgren2000definition, is a highly-important tool in guiding the decision of which medical intervention to implement. The probabilities are simply estimated by counting all posterior draws for which $\text{INB}_{\lambda} > 0$ and dividing this count by the total number of posterior draws.
For suBART and BayesSUR, inclusion of the propensity scores leads to a uniformly higher probability of cost-effectiveness. Hence, these analyses give greater certainty to the conclusion that the treatment is cost-effective. The greatest difference in this regard is observed for BayesSUR. Although mvBART does not exhibit the same effect, all three models involving propensity scores give similar results: if $\lambda$ exceeds a value of around $10000$, then the probability of cost-effectiveness exceeds $0.95$, a commonly used benchmark in CEAs. However, we stress that the superiority of any one method cannot be established from these estimated CEACs, as the true $\text{INB}_{\lambda}$ is unknown for all values of $\lambda$. In light of the simulation study based on the TTCM data presented in Section (ref), which demonstrated the superiority of ps-suBART in recovering the true INB, we tend to put most trust in the estimates from this model.
In Appendix (ref), we sketch two additional analyses of the TTCM data, which also showcase some additional features of our computational implementation.
In this paper, we introduce the suBART model for multivariate outcomes both as means of accounting for non-linearities and interactions in the seemingly unrelated regression framework and as a means of addressing the key limitation of existing multivariate BART approaches which assume a single set of trees, such that the entire response vector is partitioned in the same way by the splitting rules in the ensemble. By modelling each component of the outcome via univariate BART, suBART captures non-linearities in the relationship between the response and predictors while allowing for and detecting the different subsets of covariates --- along with the interactions between them --- associated with each response without enforcing common tree structures. Another defining feature of suBART is that it allows each individual outcome to be associated with a distinct ensemble of trees while handling dependencies between the outcomes through correlated error terms. In this respect, inference on $\boldsymbol{\Sigma}$ is key. Running independent BART models implicitly assumes that $\boldsymbol{\Sigma}$ is diagonal, which has been shown to harm performance in our analyses of the TTCM data (see Appendix (ref)).
We further develop a probit version of the suBART model to handle multivariate binary outcomes. The effectiveness of suBART and probit suBART is demonstrated through extensive simulation experiments. In the simulation experiments in Section (ref), it is shown that the model adequately captures non-linear responses, accurately estimates the covariance structure for multivariate responses, and generally outperforms its main competitors --- including other tree-based alternatives and the Bayesian linear SUR --- from the points of view of both predictive accuracy and uncertainty quantification. Notably, suBART is consistently superior to the strategy of applying the standard BART model independently to each response, which demonstrates the benefits of modelling the covariance between multivariate response under our suBART framework. Furthermore, the results show that the model exhibits enhanced flexibility compared to its direct multivariate counterpart, mvBART. This flexibility stems from suBART permitting variation in splitting rules across each response, by allowing the trees for each outcome component to differ rather than imposing common tree structures. This enables a more accurate representation, especially when different outcome components depend on distinct sets of predictors. As a cautionary note, we recall that BayesSUR exhibited competitive performance when estimating conditional expectations which merely depend linearly on the predictors. In situations where linearity can be reasonably assumed for every outcome, BayesSUR would remain a viable alternative to suBART.
The main focus of this paper is the application of suBART within the context of cost-effectiveness analysis, a setting in healthcare research where it is of interest to jointly estimate the healthcare costs and the health-related quality of life associated with two or more treatment options. Throughout Section (ref), we find some remarkable differences, depending on the particular method used for analysing the TTCM data. This applies both to the additional simulation experiment in Section (ref) and the main application in Section (ref). It is of course expected to see large differences between linear SUR and the two BART-based models, given the very different model assumptions. The significant differences between suBART and mvBART are arguably even more interesting. To reiterate, suBART assigns each outcome its own tree ensemble, while the trees for all outcomes are the same for mvBART. The large differences in results suggest to us that the assumption of a shared tree structure is a significant one, which can strongly impact the results. Unlike mvBART, suBART can accommodate data where the dependence on $\mathbf{X}$ is very different for different components of the outcome vector. In the CEA context, this applies particularly in situations where factors which govern the costs are unrelated to the quality of life, and vice versa. Such situations do occur in practise: when investigating the effect of total knee replacement, dakin2012rationing found that the patients' sex was a strong predictor of healthcare costs yet had no measurable relationship with quality of life, while the exact opposite was true for the patients' age. On the other hand, a shared tree structure may lead to estimates which are more precise, in the sense of having a smaller variance, since there are more data available to inform the tree structure. um2023bayesian claim this as an advantage of their method. Our results are somewhat inconclusive in this respect: when estimating the MATE on the costs, mvBART exhibits slightly smaller standard deviation than suBART, while both methods are similar in this respect for the quality of life and the INB.
According to our simulation experiment in Section (ref), the ps-suBART model demonstrates superior performance over competing methodologies. Recall that this extension includes propensity scores estimated via BART as an additional predictor in suBART. As suBART in particular benefits greatly from the incorporation of propensity scores, and given that this experiment is explicitly based on the TTCM data, we are inclined to believe that the ps-suBART results are the most reliable for the real TTCM application. Overall, we consider the ps-suBART model to be a natural adaptation of the univariate BART approach and expect it to be a very useful tool in the analysis of cost-effectiveness data.
Despite the extensive array of comparisons and scenarios and the additional insights gleaned by suBART in the CEA setting, there remains ample opportunity for further exploration of various extensions to suBART. We delineate some of these possibilities below in light of limitations identified in the simulation experiments and real data application.
We hope to incorporate these advancements into our forthcoming research plans and anticipate suBART's adoption in other CEA settings to inspire further developments.
\interlinepenalty=1000
\setcounter{subsection}{0}
\setcounter{figure}{0} \setcounter{table}{0}
In Section (ref), we discussed the importance of jointly modelling $c_i$ and $q_i$, so as not to enforce the unrealistic assumption that the treatment effects $\Delta_c$ and $\Delta_q$ are independent. This raises two potential questions; namely, what do we lose if we use suBART but the errors are independent, and what do we lose if the errors are not independent, but we assume that they are?
We address both questions by performing additional simulation experiments. We repeat the experiment based on the TTCM data from Section (ref), using two models to analyse the simulated datasets, for simplicity: ps-suBART, and independent BART models equipped with propensity scores as an additional covariate, with analogous priors as we do for the suBART model (i.e. independent half-$t$ priors on the standard deviations of the error terms). We refer to the latter model as ps-indBART. For the sake of fairness of the comparison, we use our own subart implementation for both models, with the code modified to accommodate the independent BART models as appropriate.
We use three different choices for the correlation coefficient $\rho\colon 0, -0.25$, and $-0.5$, where the first value corresponds to the case of independent errors and the second value of $-0.25$ has already been considered for ps-suBART in Section (ref). From Table (ref), we observe that the two models produce essentially the same results when $\rho = 0$. Hence, we do not lose anything from using suBART when the true correlation is zero. On the other hand, when the correlation is non-zero, as in Table (ref) and Table (ref), ps-suBART produces smaller bias and standard deviations, as well as improved coverage rates for all estimands. The differences in the results become more pronounced in each successive table, as the strength of the correlation increases.
In light of these results, we note that the presence of the covariance term in Equation (ref), which gives the formula for $\Pr\left(\text{INB}_{\lambda} > 0\right)$, shows how the independence assumption is likely to yield an inappropriate posterior value of $\Pr\left(\text{INB}_{\lambda} > 0\right)$. In the same way, separate BART models are likely to give inaccurate posterior quantiles for $\text{INB}_{\lambda}$ itself. This might be a cause of the lower coverage rates from ps-indBART in Tables (ref) and (ref). In this context, we draw attention to the average width of the credible intervals for $\text{INB}_{\mathbf{20k}}$ and $\text{INB}_{\mathbf{50k}}$: when $\rho$ is negative, the intervals from suBART are consistently wider than those from indBART. This is indicative of a larger posterior variance, and theoretically consistent with our specification of negative correlations. To see this, observe that $\text{Var}( \text{INB}_{\lambda}) = \lambda^2 \text{Var}(\Delta_q) + \text{Var}(\Delta_c) - 2 \lambda \text{Cov}\lbrack \Delta_q, \Delta_c\rbrack$. A negative residual correlation between the observed $q_i$ and $c_i$ values, conditional on the covariates, is likely to cause the posterior covariance $\text{Cov}\lbrack \Delta_q, \Delta_c\rbrack$ to be negative too. Hence, the posterior variance of $\text{INB}_{\lambda}$ will increase, along with the width of the credible intervals.
\setcounter{figure}{0} \setcounter{table}{0}
Throughout this paper, we adopted the consistent setting of 5000 MCMC iterations, with 1000 iterations discarded as burn-in, for our comparative experiments and application with continuous responses. Convergence appears satisfactory for all methods. However, due to the inefficiencies of the PX-MH algorithm for updating $\boldsymbol{\Sigma}$ under the probit version of suBART, we require longer chains and careful specification of the proposal tuning parameter $\nu_{\text{prop}}$ for these experiments. In particular, higher values $\nu_{\text{prop}}$ increase the acceptance rate at the expense of inducing significant autocorrelation in the chains. See zhang2006sampling for a discussion of the difficulties of optimising this trade-off. In order to maintain an acceptance rate in the range of 20-30% as the dimension and sample size varied, we found $\nu_{\text{prop}}=n_{\text{train}}/10$ (when $d=2$), $\nu_{\text{prop}}=n_{\text{train}}/2$ (when $d=3$), $N_{\text{MCMC}}=10000$, and $N_{\text{burn-in}}=2000$ to strike an appropriate balance between these two concerns and the runtime of the experiments. We hence adopted these MCMC settings for all applications of the probit version of our model and, for fairness, all other competing methods in the presence of binary responses.
In order to demonstrate the satisfactory convergence of the probit version of our model under these tuning parameter and MCMC settings, the MCMC samples for entries of $\boldsymbol{\Sigma}$ under one randomly chosen replicate for the Friedman \#2 experiment are shown in Figure (ref). Specifically, these traceplots depict the posterior samples of $\rho_{12}$ under the $d=2$ (left) and $d=3$ (right) simulation settings, with $n_{\text{train}}=1000$. In each case, the samples concentrate around their respective true values of $0.75$ and $0.8$ after the burn-in period elapses. The chosen values of $\nu_{\text{prop}}=n_{\text{train}}/10$ and $\nu_{\text{prop}}=n_{\text{train}}/2$ for the $d=2$ and $d=3$ cases achieve acceptance rates of $0.28$ and $0.25$, respectively. Similar convergence behaviour and acceptance rates are achieved for the other entries of $\boldsymbol{\Sigma}$, across all evaluated sample sizes and replicates. For brevity, we omit the corresponding traceplots.
For completeness, we also show the convergence behaviour of the continuous version of our model for the same $\rho_{12}$ under a randomly chosen replicate of the Friedman \#1 scenario, with $d=\{2,3\}$ and $n_{\text{train}}=1000$, in Figure (ref). There are evidently no convergence issues in this scenario. Indeed, further examination of traceplots across all relevant simulation experiments and configurations of sample size and dimensionality suggest that our sampler is adequate in the case of continuous responses.
\setcounter{figure}{0} \setcounter{table}{0}
In order to evaluate the performance of suBART against its competitors, we conducted experiments outlined in Section (ref), where we varied $n_{\text{train}} = n_{\text{test}} = \{250,500,1000\}$ and $d = \{2,3\}$. This appendix summarises the remaining results omitted from the main paper, for different configurations of sample size and dimensionality of the outcome. Recall that mvBART is excluded wherever $d=3$, due to the limitations of the skewBART package. Overall, the findings illustrated in the boxplots below align with the conclusions drawn in Section (ref). In particular, suBART exhibits reasonable performance metrics for both continuous (Friedman \#1) and binary (Friedman \#2) outcome scenarios, either matching or surpassing its tree-based model counterparts. Notably, suBART outperforms the linear BayesSUR across all scenarios, with the exception of the responses whose dependence on the predictors is linear. Additionally, we provide a summary of the estimated correlation parameters through tables displaying the RMSE and 50% CI coverage for all $\rho_{kj}$ where $j\neq k$. Across all scenarios, suBART consistently approaches the true values and provides credible intervals with proper coverage rates.
\setcounter{figure}{0} \setcounter{table}{0}
In Appendices (ref) and (ref), we present some additional insights gleaned from the TTCM application by our flagship ps-suBART model. In the former, we graphically explore the heterogeneity in the estimated causal effects, showing how they vary for different values of the covariates and propensity scores. In the latter, we investigate variable importance scores using the ideas of chipman2010bart. However, it is important to emphasise that we do so with reference to each individual outcome, given that the covariates (or propensity scores) may be strongly related to one outcome but not necessarily to another. This is not possible under mvBART, which enforces a common tree structure and does not permit variation in splitting rules across each response.
In our analyses of the TTCM data, we focussed on the MATE, which measures the expected treatment effect averaged over the observed sample. We did this because the interest in cost-effectiveness analysis lies primarily in average effects, not individual ones. Nonetheless, we may sometimes also want to analyse individual treatment effects, and in particular how they depend on the patients' baseline characteristics. We thus present here some ideas on how to do this, using the TTCM data. For brevity, we restrict attention to the results obtained from the flagship ps-suBART model.
We first recall the following definition, which was made in Section (ref): $$ \tau_c(\mathbf{x}_i) = \mathbb{E}\left\lbrack c_i\:\vert\:t_i = 1, \mathbf{x}_i\right\rbrack - \mathbb{E}\left\lbrack c_i\:\vert\:t_i = 0, \mathbf{x}_i\right\rbrack. $$ This quantity is called the conditional average treatment effect (CATE), because it depends on $\mathbf{x}_i$, unlike the MATE $\Delta_c$ defined in Equation (ref). We also define $\tau_q(\mathbf{x}_i)$ in the same way. In general, $\tau_c(\mathbf{x}_i)$ is different for different $i$, unless the dependence of $\mathbb{E}\left\lbrack c_i\:\vert\:t_i, \mathbf{x}_i\right\rbrack$ on $t_i$ and $\mathbf{x}_i$ is linear. See li2023bayesian for details. We now further define the conditional incremental net benefit (CINB)\footnote{There is no standard terminology for this estimand in the literature. This acronym is our own invention.} by
Again, the CINB is the counterpart of the INB defined earlier, where the former depends on $\mathbf{x}_i$ and the latter does not.
In Figure (ref) we plot the estimated propensity scores against the posterior means of $\tau_c(\mathbf{x}_i)$, $\tau_q(\mathbf{x}_i)$, and $\text{CINB}_{20}(\mathbf{x}_i)$, respectively. The plots indicate some considerable heterogeneity. In particular, we see that, while the effect on costs is negative for all patients, there seems to be a further negative relationship with the propensity scores: patients who are more likely to receive the treatment also benefit more from it. For the effect on $q_i$ there is no clear relationship with the propensity scores. Finally, the relationship between the CINB and the propensity scores is positive, which makes sense in light of the definition in Equation (ref).
In addition to the propensity scores, such plots can also be produced for any individual covariate. We show similar plots against the continuous `TTO' variable and the categorical `surgery' variable from the TTCM dataset in Figure (ref) and Figure (ref), where strong relationships and similar levels of heterogeneity are evident. Interestingly, Figure (ref) shows that, for patients who receive surgery, the treatment is both more effective and more cost-saving.
Finally, it is also worth noting that the observed heterogeneity in the above plots of $\tau_c(\mathbf{x}_i)$ against various predictors is indicative of non-linearity in $\mathbb{E}\left\lbrack c_i\:\vert\:t_i, \mathbf{x}_i\right\rbrack$. The same applies to $\tau_q(\mathbf{x}_i)$ and $\text{CINB}_{20}(\mathbf{x}_i)$. Thus, there is evidence to support the flexibility of ps-suBART in capturing non-linear relationships.
Figure (ref) shows the percentage of times, averaged across all trees in the outcome-specific ensemble, that a given predictor was used to form a splitting rule for a given outcome in the TTCM application. Similar plots were shown to demonstrate the quantification of variable importance in the original univariate BART setting by chipman2010bart. Here, the variable importance scores in Figure (ref) illustrate the practical implications of a defining feature of suBART, namely that the relative importance of the predictors differs in relation to each outcome. Crucially, this behaviour is achieved without having to pre-specify the functional form of the relationships between the cost and quality outcomes and the available predictors, as a result of BART's intrinsic variable-selection properties. Some interesting findings are worth briefly commenting on. Firstly, we have variables which are important to both outcomes, such as the treatment indicator. Secondly, we have variables which are relatively unimportant for both outcomes, with the hospital admission indicator being the second-least and least used variable, respectively. Thirdly, we have variables which differ in the extent of their importance for each outcome: trauma type is the most used variable for the $c_i$ outcome, while being number six in the ranking for $q_i$. Finally, we note that tailoring the suBART framework to causal inference objectives by supplying propensity scores as a predictor has indeed led to the propensity scores being frequently used to form splitting rules for both outcomes.