The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
54,472 characters
Type 2 Tobit Sample Selection Models with Bayesian Additive Regression Trees
\title{Type 2 Tobit Sample Selection Models with Bayesian Additive Regression Trees}
\author{\hspace{-0pt} Eoghan O'Neill$^{a}$\thanks{Email: [email removed] . Acknowledgements: The author gratefully acknowledges helpful comments from Mikhail Zhelonkin, Chen Zhou, Richard Paap, Dennis Fok, Adam Iqbal, F. Javier Rubio, and Emmanuel O. Ogundimu.}
\\
\vspace{-7pt}
\\
{\centering {\small{$^a$University College Dublin}}}}
\vspace{-5pt}
\date{\today}
\maketitle
\vspace{-10pt}
\begin{abstract}
This paper introduces Type 2 Tobit Bayesian Additive Regression Trees (TOBART-2). BART can produce accurate individual-specific treatment effect estimates. However, in practice estimates are often biased by sample selection. We extend the Type 2 Tobit sample selection model to account for nonlinearities and model uncertainty by including sums of trees in both the selection and outcome equations. A Dirichlet Process Mixture distribution for the error terms allows for departure from the assumption of bivariate normally distributed errors. Soft trees and a Dirichlet prior on splitting probabilities improve modeling of smooth and sparse data generating processes. We include a simulation study and an application to the RAND Health Insurance Experiment dataset.
\end{abstract}
\newpage
\section{Introduction}
Treatment effect estimates obtained from observational data are often biased by non-random sample selection. \citet{heckman1974shadow} proposed to jointly model the outcome and the selection of variables to produce unbiased estimates.{\footnote{The Heckman selection model was named the Type 2 Tobit model by {\citet{amemiya1984tobit}} .
}} However, this model assumes a linear functions of covariates in the selection and outcome equations, and joint normality of the error terms.
To address these issues, we introduce a type 2 Tobit model with separate sums of trees in the selection and outcome equations instead of linear combinations of covariates. This model is henceforth referred to as TOBART-2. This builds on the Bayesian formulation of type 2 Tobit models introduced by \cite{omori2007efficient, chib2009estimation}, and \cite{van2011bayesian}. The Markov chain Monte Carlo implementation is most similar to that of \cite{van2011bayesian}. However, Type 2 Tobit models that impose a normality assumption can produce biased estimates when applied to data with small deviations from normality.\footnote{Many other papers have provided alternatives to the assumption of joint normality of the errors in sample selection models. \cite{gallant1987semi} use Hermite series, \cite{marchenko2012heckman} use a t-distribution, \cite{ding2014bayesian} applies a t-distribution in a Bayesian model, and \cite{ogundimu2016sample} consider a skew-normal distribution. \cite{zhelonkin2016robust} propose an outlier-robust estimator.} A more robust method, TOBART-2-NP, contains a Dirichlet Process Mixture (DPM) of bivariate Gaussian distributions for the error term, as described by \cite{van2011bayesian}.\footnote{A similar DPM of normal distributions is used to model the error terms in the instrumental variable BART model of \cite{mcculloch2021causal}, which is based on \cite{conley2008semi}.}
The advantages of TOBART-2 are that it accounts for nonlinearity, variable selection, and model uncertainty, while flexibly modelling the errors to allow for non-normality. In addition to point estimates, the method produces credible intervals for outcome predictions and treatment effect estimates. As noted by \cite{van2011bayesian}, posterior inference can be applied to the covariance of error terms to indicate the extent of sample selection.
Variable importance measures can be obtained for the selection and outcome equations.
Data analysis competitions and other simulation studies have demonstrated the impressive performance of BART \citep{hahn2019atlantic, dorie2019automated, wendling2018comparing} in estimating treatment effects. BART models have been formulated to account for confounding \citep{hill2011bayesian, hahn2020bayesian, kim2022bayesian}, sparsity \citep{linero2018bayesianB}, smoothness \citep{linero2018bayesianA}, and unobserved confounders with instrumental variables \citep{mcculloch2021causal}. To the best of our knowledge, the only other paper that describes a BART model that allows for a form of selection on unobservables is the multivariate missing not at random outcome model of \cite{goh2024joint}, created concurrently with TOBART-2.\footnote{Some papers describe BART-based methods for data with selection on observables, e.g. \cite{wang2023improving, gao2021treatment, elliott2023improving, bisbee2019barp}. } However, \cite{goh2024joint} model missingness as a function of covariates and latent outcomes, and do not model selection through a joint distribution of error terms as in standard Heckman selection model formulations.\footnote{TOBART-2 is arguably a more interpretable model with a more clearly specified form of selection on unobservables, and it allows for direct specification of the prior on selection and allows for inference on the level of selection on unobservables. The method of \cite{goh2024joint} is not generalized to non-normal errors, although such an extension would be straightforward. The general approach to selection modelling of TOBART-2 is more common in the econometrics literature, whereas the approach of \cite{goh2024joint} is more often applied in other disciplines.}
Bayesian Model Averaging of linear Type 2 Tobit models has been considered by \cite{jordan2012tobit} and \cite{ eicher2012robust}. Nonlinear approaches to Type 2 Tobit models include kernel regression \citep{ahn1993semiparametric, newey1990semiparametric, lee1994semiparametric, chen2010semiparametric}, splines \citep{marra2013estimation,das2003nonparametric, newey2009two}, and neural networks \citep{zhang2021deep}.\footnote{However, \cite{zhang2021deep} do not evaluate their methods using simulations with non-zero correlation of the error terms in the selection and outcome equations. i.e. there is no selection on unobservables.} \cite{chib2009estimation} formulate a semiparametric Bayesian model with sample selection and endogeneity by using a second-order Markov Process prior. However, \cite{chib2009estimation} apply the restrictive assumption of jointly normally distributed errors. \cite{kim2019bayesian} formulate a Bayesian sample selection model with a scale mixture of bivariate normal distributions for the error term, and a Berstein polynomial regression model for the outcome. However, the selection equation contains a linear function of covariates. \cite{wiesenfarth2010bayesian} make use of penalized splines in the outcome and selection equations of a Bayesian sample selection model, although the errors are restrictively modelled as being jointly normally distributed. \cite{wiemann2022correcting} introduce a Bayesian method with splines and a bivariate copula.
\cite{fan2007sample} apply trees and model averaging to data with sample selection. \cite{zadrozny2004learning, tran2017selection}, and many others combine machine learning methods with inverse probability weighting approaches that cannot account for selection on unobservables. \cite{alaimo20222} apply gradient boosted trees to tax revenue data with sample selection, and account for model uncertainty by using a bootstrap procedure. However, this approach is also based on the inverse propensity methods of \cite{zadrozny2004learning} and therefore does not account for selection on unobservables.
\cite{brewer2024addressing} describe control function approaches that allow for machine learning methods to be combined with a Heckman-style two-step framework. This method leads to improved predictive accuracy in simulations with moderate to high levels of selection. However, this approach relies on the assumption of joint normality and \cite{brewer2024addressing} do not discuss efficiency or methods for inference. \cite{zadrozny2001learning} similarly apply a control function approach with decision trees. \cite{zhu2017nonasymptotic} apply a LASSO penalty using a control function approach and semiparametric framework. \cite{ogundimu2022lasso, ogundimu2024lasso} applies LASSO and adaptive LASSO in both the selection and outcome equations. \cite{iqbal2023bayesian} apply spike-and-slab priors to the coefficients in the selection and outcome equations to produce a Bayesian sample selection model for high dimensional data. \cite{bia2024double} describe a double machine learning approach to estimation of the average treatment effect for high-dimensional data with sample selection.
In addition to TOBART-2, this paper introduces the following:
\begin{itemize}
\item To allow for data-informed calibration of the prior variance of the outcome equation errors, the prior of \cite{ding2014bayesian} is altered to a bivariate normal prior with an additional data-informed hyperparameter. This is motivated by the fact that good performance of standard BART without hyperparameter tuning is partly attributable to well-calibrated hyperparameters. The sampler of \cite{ding2014bayesian} is also improved by removing unnecessary imputation of the missing outcomes.\footnote{This was noted as an area for future research by \cite{ding2014bayesian}.}
\item We provide new insights pertaining to the \cite{van2011bayesian} and \cite{omori2007efficient} Tobit-2 priors. The CDF of the correlation between selection and outcome errors is derived for both priors. We derive the conditions in which the implied prior on the correlation specified by \cite{van2011bayesian} is unimodal, and when it is concentrated at +1 and -1. The prior CDF of the outcome variance is derived for prior calibration. The prior correlation distribution specified by \cite{omori2007efficient} can take various unimodal and bimodal forms, depending non-trivially on hyperparameters.\footnote{This contrasts with a note by \cite{van2011bayesian}, which perhaps refers to particular hyperparameter values.} These results and prior calibration are discussed in Appendix \ref{imp_details_app}.
\item We consider a variation on the MCMC sampler of \cite{van2005bayesian, van2011bayesian} in which the outcome equation variance and the covariance between the selection and outcome equation errors are jointly sampled. See Appendix \ref{phigamma_norm_ig_sec}.
\item We describe an alternative TOBART-2 implementation that extends the marginalized BART method described by \cite{collins2023improved} to a sample selection model. This implementation, TOBART2marg, marginalizes out both the terminal node parameters and the covariance term determining sample selection when updating trees in the outcome equation, and then makes a joint draw of all marginalized parameters. This is motivated by the observation of \cite{van2005bayesian, van2011bayesian} in the linear model setting that there is high dependence between the outcome equation coefficient draws and the covariance parameter draws. See Appendix \ref{TOBART2Marg_app}.
\end{itemize}
An \texttt{R} package containing implementations of TOBART-2 and the Bayesian linear model described by \cite{van2011bayesian} is available at \url{https://github.com/EoghanONeill/TobitBART}. The remainder of the paper is structured as follows: In section \ref{methods_sec} we describe the TOBART-2 model and MCMC implementation, in section \ref{criteria_sec} we discuss the choice of estimand and evaluation metrics, section \ref{simulation_sec} contains a simulation study, section \ref{application_sec} contains an application to a real world data, and section \ref{conclusion_sec} concludes the paper.
\section{Methods}\label{methods_sec}
\subsection{Review of Bayesian Additive Regression Trees (BART)}
\noindent \textbf{Description of Model and Priors}
\medskip
Suppose there are $n$ observations, and the $n \times p$ matrix of explanatory variables, $X$, has $i^{th}$ row $x_i=[x_{i1},...,x_{ip}]$. Following the notation of \cite{chipman2010bart}, let $T$ be a binary tree consisting of a set of interior node decision rules and a set of terminal
nodes, and let $M = \{ \mu_1 , ..., \mu_b \}$ denote a set of parameter values associated with each of the $b$ terminal nodes of $T$. The interior node decision rules are binary splits of the predictor space into the sets $\{ x_{is} \le c \}$ and $\{ x_{is} > c \}$ for continuous $x_{s}$. Each observation's $x_i$ vector is associated with a single terminal node of $T$, and is assigned the $\mu$ value associated with this terminal node. For a given $T$ and $M$, the function $g(x_i;T,M)$ assigns a $\mu \in M$
to $x_i$. This gives the single tree model $Y \sim g(x_i;T,M) + \varepsilon \ , \ \varepsilon \sim N(0,\sigma^2)$ \citep{chipman1998bayesian}. For the standard BART model, the outcome is determined by a sum of trees,
$$Y_i = f(x_i) +\varepsilon_i = \sum_{j=1}^m g(x_i ; T_j, M_j)+\varepsilon_i \ \ , \ \ \varepsilon_i \overset{i.i.d}{\sim} N(0, \sigma^2)$$
where $g(x_i;T_j,M_j)$ is the output of a decision tree. $T_j$ refers to a decision tree indexed by $j=1,...,m$, where $m$ is the total number of trees in the model. $M_j$ is the set of terminal node parameters of $T_j$.
Prior independence is assumed across trees $T_j$ and across terminal node means $M_j = (\mu_{1j}...\mu_{b_j j})$, where $1,...,b_j$ indexes the terminal nodes of tree $j$. The prior of \cite{chipman2010bart} has the form:
$$p(M_1,...,M_m,T_1,...,T_m,\sigma) \propto \left[ \prod_j \left[ \prod_k p(\mu_{kj}|T_j) \right] p(T_j)\right]p(\sigma) $$
In standard BART, $\mu_{kj} | T_j \overset{i.i.d}{\sim} N(0,\sigma_0^2)$ where $\sigma_0 = \frac{0.5}{e \sqrt{m}}$ and $e$ is a user-specified hyper-parameter.
\cite{chipman2010bart} set a regularization prior on the tree size and shape $p(T_j)$. The probability that a given node within a tree $T_j$ is split into two child nodes is $\alpha (1+d_h)^{-\beta}$, where $d_h$ is the depth of internal node $h$, and the parameters $\alpha$ and $\beta$ determine the size and shape of $T_j$ respectively. \cite{chipman2010bart} use uniform priors on available splitting variables and splitting points.
\cite{chipman2010bart} assume that the model precision $\sigma^{-2}$ has a conjugate prior distribution $\sigma^{-2} \sim Ga(\frac{v}{2}, \frac{v \lambda}{2})$ with degrees of freedom $v$ and scale $\lambda$.
BART predictions are averages of sum-of-tree models. Therefore model uncertainty is taken into account and there are two levels of regularization. Firstly, greater prior probability is placed on models with shallower trees with fewer splitting points. Secondly, over-fitting is further avoided through the prior on the terminal node parameters $\mu_{kj}$, as in standard Bayesian linear regression.
\bigskip
\noindent \textbf{BART Implementation}
\medskip
Samples can be taken from the posterior distribution $p((T_1, M_1),...,(T_m,M_m), \sigma | y)$ by a Bayesian backfitting MCMC algorithm. This algorithm is a Gibbs sampler, involving $m$ successive draws from
$(T_j , M_j )| T_{(j)} , M_{(j)} , \sigma , y $ for $j=1,...,m$, where $T_{(j)} , M_{(j)} $ are the trees and parameters for all trees except the $j^{th}$ tree, followed by a draw of $\sigma $ from the full conditional $\sigma | T_1,...,T_m,M_1,...,M_m,y$.
The $q^{th}$ Gibbs iteration defines the sum of trees function $f_q^*(.)= \sum_{j=1}^m g_q(. \ ; T_{j}^*, M_{j}^*)$. After burn-in, the sequence of $f^*$ draws, $f_1^*,...,f_Q^*$ may be regarded as an approximate, dependent sample of size $Q$ from $p(f|y)$.
To estimate the unknown function $f(x)$,
the expectation $E(f(x)|y)$ is approximated by
$\frac{1}{Q} \sum_{q=1}^Q f_{q}^* (x)$. Posterior credible intervals can be obtained from quantiles of the draws $f_{q}^* (x)$.
\subsection{Type 2 Tobit}
\subsubsection{Type 2 Tobit and TOBART Model}
In the standard (linear) Type 2 Tobit model \citep{omori2007efficient}, separate parameters and variables enter the selection and outcome equations:
$$
y_i =
\begin{cases}
y_i^* & \text{if } z_i^* \ge 0\\
n.a. & \text{otherwise}
\end{cases}
, \ z_i^* = \bm{w}_i'\bm{\theta} + \xi_i , \ y_i^* = \bm{x}_i'\bm{\beta} + \eta_i , \ \begin{pmatrix}
\xi_i \\
\eta_i
\end{pmatrix}
\sim \mathcal{N}(\bm{0}, \Sigma) \text{ where } \Sigma =
\begin{pmatrix}
1 & \gamma \\
\gamma & \phi + \gamma^2
\end{pmatrix} $$
In Type 2 TOBART, the linear combinations $\bm{w}_i'\bm{\theta}$ and $\bm{x}_i'\bm{\beta} $ are replaced by sums-of-trees denoted by $f_z(\bm{w}_i)$ and $f_y(\bm{x}_i)$ respectively.
$$z_i^* = f_z(\bm{w}_i) + \xi_i \ , \ y_i^* = f_y(\bm{x}_i) + \eta_i $$
For linear Type 2 Tobit, the coefficient priors are $ \bm{\theta} \sim \mathcal{N}(\bm{\theta}_0, \bm{\Theta}_0) \ , \ \bm{\beta} \sim \mathcal{N} (\bm{\beta}_0, \bm{B}_0 ) $, and for TOBART the priors are $f_z \sim BART $ and $f_y \sim BART $. The covariance matrix parameters have the priors $ \gamma \sim \mathcal{N} (\gamma_0, G_0) $ and $ \phi \sim IG(\frac{n_0}{2}, \frac{S_0}{2}) $. An alternative prior, $\gamma | \phi \sim \mathcal{N} (\gamma_0, \tau \phi)$, is introduced by \cite{van2011bayesian}, who notes that dependence between $\gamma$ and $\phi$ through $\tau>0$ allows various shapes of the prior on correlation between $\xi$ and $\nu$, and an induced prior correlation with less probability mass near $\pm 1 $.
Details of the Type 2 Tobit and TOBART Gibbs sampler are provided in appendix \ref{app_mcmc_tbart2norm}.
We also consider a prior for the distribution of the error terms similar to the prior introduced by \cite{ding2014bayesian}. \cite{ding2014bayesian} assumes that the errors have a bivariate t-distribution:
$$
\begin{pmatrix}
\xi_i \\
\eta_i
\end{pmatrix}
\sim \ \mathcal{N}(\bm{0}, \alpha \Omega/q_i) \ , \ \text{where} \ q_i \sim \alpha \chi_{\nu}^2/\nu \ , \ i=1,\dots,N
$$
where $ \Omega = \begin{pmatrix}
1 & \rho \tilde{\sigma}_2 \\
\rho \tilde{\sigma}_2 & \tilde{\sigma}_2^2 \\
\end{pmatrix} $
and the prior is $ \alpha \sim b / \xi_{c}^2 $ and
$ \tilde{\Sigma} = \begin{pmatrix}
\tilde{\sigma}_1 & 0 \\
0 & 1 \\
\end{pmatrix} \Omega \begin{pmatrix}
\tilde{\sigma}_1 & 0 \\
0 & 1 \\
\end{pmatrix}
\sim W_2^{-1} (\nu_0, I_2) $,
which is equivalent to $p(\Omega) \sim (1 - \rho^2)^{-3/2} \tilde{\sigma}_2^{- (\nu_0 +3 )} \exp \Big\{ - \frac{1}{2 \tilde{\sigma}_2^2 (1- \rho^2) } \Big\}$ and $\tilde{\sigma}_1^2|\Omega \sim \{ (1 - \rho^2) \chi_{\nu_0}^2 \}^{-1} $ , where $\tilde{\sigma}_1 = \sqrt{\tilde{\Sigma}_{1,1}} $ , and $\tilde{\sigma}_2 = \sqrt{\tilde{\Sigma}_{2,2}} = \sqrt{\Omega_{2,2}} $.
However, instead of specifying a t-distribution, we maintain a normal specification for the error term, and allow for departures from normality through a mixture of normal distributions.
$$
\begin{pmatrix}
\xi_i \\
\eta_i
\end{pmatrix}
\sim \ \mathcal{N}(\bm{0}, \Omega) \ , \ \tilde{\Sigma} = \begin{pmatrix}
\tilde{\sigma}_1 & 0 \\
0 & 1 \\
\end{pmatrix}
\Omega \begin{pmatrix}
\tilde{\sigma}_1 & 0 \\
0 & 1 \\
\end{pmatrix}
\sim W_2^{-1} (\nu_0, I_2)
$$
As noted by \cite{ding2014bayesian}, this prior imposes a $U([-1,1])$ marginal prior on $\rho = \text{corr} (\xi_i, \eta_i) $ when $\nu_0 = 3$. In contrast, the priors of \cite{van2011bayesian} and \cite{omori2007efficient} can be notably non-uniform.\footnote{The prior of \cite{van2011bayesian} can place most mass around zero, or closer to $\pm 1$ depending on the value of $\tau$. See Appendix \ref{imp_details_app} and the appendix of \cite{iqbal2023bayesian} for further discussion.}
\bigskip
However, this prior does not allow for adequate calibration of the (marginal) prior on $\tilde{\sigma}_2 $ relative to the data-informed prior calibration described by \cite{chipman2010bart} and \cite{mcculloch2021causal}. In particular, it can be shown that (see appendix \ref{marginal_prior_app} for derivation)
$$ p(\tilde{\sigma}_2^2 ) \propto \left( \frac{1}{\tilde{\sigma}_2^2 } \right)^{\frac{\nu_0-1}{2} + 1}
\exp \left( - \frac{1}{ 2 \tilde{\sigma}_2 ^2} \right) $$
therefore the marginal prior is $\tilde{\sigma}_2^2 \sim \Gamma^{-1}( \frac{\nu_0-1}{2} , \frac{1}{2})$. This prior does not allow for both a uniform marginal prior on $\rho$ and a data-informed marginal prior on $\tilde{\sigma}_2^2 $. Therefore, we specify the following prior
$$
\begin{pmatrix}
\xi_i \\
\eta_i
\end{pmatrix}
\sim \ \mathcal{N}(\bm{0}, \Omega) \ , \ \tilde{\Sigma} = \begin{pmatrix}
\tilde{\sigma}_1 & 0 \\
0 & 1 \\
\end{pmatrix}
\Omega \begin{pmatrix}
\tilde{\sigma}_1 & 0 \\
0 & 1 \\
\end{pmatrix}
\sim W_2^{-1} (\nu_0, c I_2)
$$
where we have introduced the hyperparameter $c$. Similarly to \cite{ding2014bayesian}, we can re-express this prior as
$p(\Omega) \sim (1 - \rho^2)^{-3/2} \tilde{\sigma}_2^{- (\nu_0 +3 )} \exp \Big\{ - \frac{c}{2 \tilde{\sigma}_2^2 (1- \rho^2) } \Big\}$ and $\tilde{\sigma}_1^2|\Omega \sim \{ \frac{(1 - \rho^2)}{c} \chi_{\nu_0}^2 \}^{-1} $ or $\tilde{\sigma}_1^2|\Omega \sim \Gamma^{-1} \left( \frac{\nu_0}{2}, \frac{c}{2(1-\rho^2)} \right) $. It can be shown (see Appendix \ref{marginal_prior_app}) that this implies the following marginal priors:
$$ p(\rho) \propto (1 - \rho^2)^{\frac{\nu - 3}{2} } \ \text{ and } \ p(\tilde{\sigma}_2^2) \propto \left( \frac{1}{\tilde{\sigma}_2^2} \right)^{\frac{\nu_0-1}{2} + 1} \exp \left( - \frac{ c }{ 2 \tilde{\sigma}_2^2 } \right) $$
This prior allows us to set $\nu_0 =3$ for a uniform marginal $\rho$ distribution and set $c$ such that the $q^{th}$ quantile of the marginal prior on $\tilde{\sigma}_2^2$, i.e. $\Gamma^{-1}( \frac{\nu_0-1}{2}, \frac{c}{2}) $, is equal to an estimated variance of $\eta_i$, denoted by $\hat{\sigma}_2^2$. This is similar to the prior calibration approach of \cite{chipman2010bart} and \cite{mcculloch2021causal}. We set $q = 0.95$. Setting $\hat{\sigma}_2^2$ equal to the sample variance of the observed outcomes might not give a well-calibrated prior if the variance of the observed outcomes is far from the variance of the full set of outcomes, therefore it is preferable to use the estimated outcome variance from a linear Tobit model estimated by maximum likelihood.\footnote{The value of $c$ is $\hat{\sigma}_2^2$ multiplied by the $(1-q)^{th}$ quantile of the $\chi_{\nu_0}^2$ distribution. We also considered a number of data-informed prior calibration methods for the priors of \cite{omori2007efficient} and \cite{van2011bayesian}. See Appendix \ref{imp_details_app}.}
\subsection{Nonparametric Type 2 Tobit}
\subsubsection{Nonparametric Type 2 Tobit and TOBART Models}
The models presented here are based on the semiparametric Bayesian Type 2 Tobit model implementation introduced by \cite{van2011bayesian}.\footnote{The notation is more similar to that of \cite{omori2007efficient}.} In the standard (linear) Type 2 Tobit model, separate parameters and variables enter the selection and outcome equations:
$$
y_i =
\begin{cases}
y_i^* & \text{if } z_i^* \ge 0\\
n.a. & \text{otherwise}
\end{cases} , \ z_i^* = \bm{w}_i'\bm{\theta} + \xi_i , \ y_i^* = \bm{x}_i'\bm{\beta} + \eta_i, \begin{pmatrix}
\xi_i \\
\eta_i
\end{pmatrix}
\overset{i.i.d.}{\sim} \mathcal{N}\Bigg(
\begin{pmatrix}
\mu_{i1} \\
\mu_{i2}
\end{pmatrix}
, \Sigma_i =
\begin{pmatrix}
1 & \gamma_i \\
\gamma_i & \phi_i + \gamma_i^2
\end{pmatrix}\Bigg)
$$
$$\{
\begin{pmatrix}
\mu_{i1} &
\mu_{i2}
\end{pmatrix}'
, \gamma_i, \phi_i \} | H \sim H \ , \ H | \alpha, H_0 \sim \mathcal{DP} (\alpha, H_0)$$
In Type 2 TOBART, the linear combinations $\bm{w}_i'\bm{\theta}$ and $\bm{x}_i'\bm{\beta} $ are replaced by sums-of-trees denoted by $f_z(\bm{w}_i)$ and $f_y(\bm{x}_i)$ respectively.
$$z_i^* = f_z(\bm{w}_i) + \xi_i \ , \ y_i^* = f_y(\bm{x}_i) + \eta_i \ , \
$$
For linear Tobit, the priors are $ \bm{\theta} \sim \mathcal{N}(\bm{\theta}_0, \bm{\Theta}_0) \ , \ \bm{\beta} \sim \mathcal{N} (\bm{\beta}_0, \bm{B}_0 ) $, and for TOBART $f_z \sim BART $ and $f_y \sim BART $. $H$ is a discrete distribution of the parameters $ \vartheta_i = \{
\begin{pmatrix}
\mu_{i1} &
\mu_{i2}
\end{pmatrix}'
, \gamma_i, \phi_i \} $ . The base measure $H_0$ is defined by
$\begin{pmatrix}
\mu_{i1} &
\mu_{i2}
\end{pmatrix}' \sim \mathcal{N} (\bm{0}, \Omega) \ , \ \gamma_i \sim \mathcal{N} (\gamma_0, G_0) \ , \ \phi_i \sim IG(\frac{n_0}{2}, \frac{S_0}{2}) $ and
\cite{van2011bayesian} sets $\Omega = 10 \bm{I}_2$.\footnote{An alternative would be to set it equal to $\frac{\hat{\sigma}_y^2}{0.016} \bm{I}_2$ such that the prior probability that a $\mu_{i}$ value lies in $(-10,10)$ is $0.8$. This is similar to the hyperparameter settings described by \cite{chipman1998bayesian} and \cite{conley2008semi}, albeit without a conjugate prior. The other base distribution hyperparameters are set to the same values as for standard TOBART-2. Alternatively, as suggested by \cite{george2019fully}, scale and shape of the prior variance of the outcome error can be set to values that imply a tighter base distribution, as the mixture distribution creates additional spread for the overall prior.} However, since $ f_z(\bm{w}_i) $ is non-linear, we can set $\mu_{i1}$ to 0 and still have a flexible selection model, $\Pr(z_i^* >0) = \Phi ( f_z(\bm{w}_i))$, with a standard normal marginal distribution of $\xi_i$ for identification of $f_z(\bm{w}_i)$.
The parameter $\alpha$ determines the distribution of the weights given to each element of $H$. As $\alpha$ increases, more elements receive non-negligible weight.
As in the parametric error model, \cite{van2011bayesian} uses the alternative prior $\gamma_i | \phi_i \sim \mathcal{N} (\gamma_0, \tau \phi_i)$.
\cite{george2019fully} and \cite{van2011bayesian} apply different priors to $\alpha$. \cite{george2019fully} choose the maximum and minimum number of components in $G$, denoted by $I_{min}$ and $I_{max}$, then solve for $\alpha_{min}$ such that the mode of the number of elements (denoted by $I$) of $H$ is $I_{min}$. Similarly $\alpha_{max}$ is derived from $I_{max}$. Then
$p(\alpha) \propto (1- \frac{ \alpha - \alpha_{min}}{\alpha_{max} - \alpha_{min}} )^{\psi} $ .
The default values are $I_{min} = 1$, $I_{max} = [(0.1)*n]$, where $[\bullet]$ denotes the floor function, and $\psi = 0.5$.
\cite{mcculloch2021causal} also specify a DPM prior on the errors of an IV-BART model and set $I_{min} = 2$, $I_{max} = [(0.1)*n] + 1$.
\cite{van2011bayesian} applies the prior from \cite{escobar1994estimating}, $\alpha \sim \text{Gamma}(c_1, c_2)$, with $c_1=2$ and $c_2=2$. Details of the sampler introduced by \cite{van2011bayesian} are provided in appendix \ref{app_mcmc_tbart2np}. We apply the sampler of \cite{van2011bayesian} by default and provide an option for the prior of \cite{george2019fully}.
\subsubsection{Type 2 TOBART Gibbs Sampler}
Samples of $y_i^*$ and $z_i^*$ in the Gibbs sampler for TOBART are the same as in the standard Tobit sampler. However, in contrast to the draws of coefficient parameters for linear Tobit, each tree must be drawn separately conditional on all other trees. Detailed outlines of the full TOBART-2 and TOBART-2-NP samplers are included in appendices \ref{app_mcmc_tbart2norm} and \ref{app_mcmc_tbart2np} respectively. An alternative sampler, in which the outcome equation tree draws are made unconditional on the terminal nodes and covariance term $\gamma$, and in which the outcome equation leaf parameters and $\gamma$ are sampled jointly, is described in appendix \ref{TOBART2Marg_app}. This alternative sampler is intended to reduce the dependence between the outcome equation tree samples and $\gamma$ samples. The sampler of \cite{van2011bayesian} is similarly motivated.
\subsubsection{Inference about degree of dependence}
As noted by \cite{van2011bayesian}, the samples of the error terms allow for inference on the degree of dependence between the errors in the selection and outcome equations.
At each iteration $t$ of the TOBART-2-NP Gibbs sampler, we may generate a pseudo-sample $\{u_{i,t} \}_{i=1}^n$ where $u_{i,t} \sim \mathcal{N}(\mu_{i,t}, \Sigma_{i,t})$ and calculate a dependence measure. The set of values of this measure across all iterations of the Gibbs sampler are an approximate sample from its posterior distribution.
\subsection{Note on prior hyperparameters}
\cite{iqbal2023bayesian} apply the \cite{van2011bayesian} prior with $n_0 = 2$ and $S_0= 2$. It is claimed that these values give an induced prior on the outcome variance $ \phi + \gamma^2$ that is similar to an Inverse Gamma. However, the induced prior mean is $ E[ \phi + \gamma^2 ] = (1 + \tau)\frac{S_0}{n_0-2} + g_0^2 $ and therefore is undefined if $n_0=2$, which it would also be for the Inverse Gamma. We note in Appendix \ref{imp_details_app} that the prior distribution of $ \phi + \gamma^2$ is a non-trivial sum of an inverse gamma distributed variable and a (dependent) square of a location and scale shifted t-distributed variable.
Appendix \ref{imp_details_app} contains an additional discussion of some methods for calibration of $n_0$ and $S_0$, and of potential disadvantages of the recommendation of \cite{iqbal2023bayesian} to set $\tau=5$ in general. In particular, large values of $\tau$ increase the prior mean of of $\phi + \gamma^2$. Moreover, we observe that larger values of $\tau$ can result in a more peaked bimodal prior for $\rho$ with modes closer to $\pm1$.\footnote{An alternative would be to apply a hyperprior to $\tau$. We did not find this to be useful.}
\subsection{Note on Estimands and Evaluation criteria}\label{criteria_sec}
\subsubsection{Prediction Estimands}
We assume that the data generating process is as described for the Type 2 TOBART model. To predict the observed outcome, $Y_i$ for selected observations, we may use the following estimand:
$$ E[Y_i | Z_i^* \ge 0 , \bm{x}_i, \bm{w}_i ] =
f_{y}(\bm{x}_i) + E[\eta_i |f_z(\bm{w}_i) + \xi_i \ge 0 ] = f_{y}(\bm{x}_i) + \gamma \frac{\phi(f_z(\bm{w}_i) )}{\Phi(f_z(\bm{w}_i) )} $$
We expect any nonlinear machine learning method trained on the selected observations to estimate this equally as well as TOBART-2.
For the nonparametric TOBART-2 DGP, we instead have
$$ E[Y_i | Z_i^* \ge 0 , \bm{x}_i, \bm{w}_i,
\mu_{i1} ,
\mu_{i2}
, \Sigma_i ] = $$
$$ f_{y}(\bm{x}_i) + E[\eta_i |f_z(\bm{w}_i) + \xi_i \ge 0 ,
\mu_{i1} ,
\mu_{i2}
, \Sigma_i ] =
f_{y}(\bm{x}_i) + \mu_{i2} + \gamma_i \frac{\phi(f_z(\bm{w}_i) + \mu_{i1})}{\Phi(f_z(\bm{w}_i) + \mu_{i1})} $$
Then we predict the outcome by averaging over MCMC draws indexed by $d=1,\dots,D$
$$ E[Y_i | Z_i^* \ge 0, \bm{x}_i, \bm{w}_i ] \approx \frac{1}{D} \sum_{d=1}^{D} E[Y_i | Z_i^* \ge 0 , f_{y}^{(d)}, f_{z}^{(d)},
\mu_{i1}^{(d)} ,
\mu_{i2}^{(d)}
, \Sigma_i^{(d)}, \bm{x}_i, \bm{w}_i ] = $$
$$ \frac{1}{D} \sum_{d=1}^{D} \Bigg\{ f_{y}^{(d)}(\bm{x}_i) + \mu_{i2}^{(d)} + \gamma_i^{(d)} \frac{\phi(f_z^{(d)}(\bm{w}_i) + \mu_{i1}^{(d)})}{\Phi(f_z^{(d)}(\bm{w}_i) + \mu_{i1}^{(d)})} \Bigg\} $$
If the task is to predict the latent outcome regardless of whether the test observation will be selected or not, or if there will be no selection mechanism in the test data, and instead all outcomes will be observed, then we predict $Y_i$ with the following estimand for the TOBART-2 DGP:
$$ E[Y_i^* | \bm{x}_i, \bm{w}_i ] = f_{y}(\bm{x}_i) $$
Therefore, for this DGP, a machine learning method trained on the selected data would target an estimand containing a bias term $ \gamma \frac{\phi(f_z(\bm{w}_i) )}{\Phi(f_z(\bm{w}_i) )} $. In contrast, a Tobit-2-based approach, such as TOBART-2, directly models $f_y$, and we can use the estimate $\frac{1}{D} \sum_{d=1}^{D} f_{y}^{(d)}(\bm{x}_i) $ for standard TOBART-2 and $\frac{1}{D} \sum_{d=1}^{D} \Big\{ f_{y}^{(d)}(\bm{x}_i) + \mu_{i2}^{(d)} \Big\} $ for TOBART-2-NP.
\subsubsection{Treatment Effect Estimands}
Potential outcomes under treatment and control allocation are denoted by $Y_i(1)$ and $Y_i(0)$ respectively. Similarly, $Y_i^*(1)$ and $Y_i^*(0)$ denote potential latent outcomes, and $Z_i^*(1)$ and $Z_i^*(0)$ denote potential selection equation latent variable values. The binary treatment variable is denoted by $T_i$.
Generally, the estimand of interest is an effect on the latent outcome $Y_i^*$, not an effect on the observed outcome $Y_i$. If treatment affects the latent outcome, but has no effect on selection, i.e. we now have $f_y(\bm{x}_i, T_i)$ and $f_z(\bm{w}_i)$, then the estimand of interest for the TOBART-2 DGP is:
$$ E[Y_i^*(1) | \bm{x}_i, \bm{w}_i ] - E[Y_i^*(0) | \bm{x}_i, \bm{w}_i ] $$
Assume that there is unconfoundedness of treatment, possibly conditional on covariates, i.e. $ Y_i^*(1), Y_i^*(0) \perp T_i | \bm{x}_i$. Then TOBART-2 directly models the above estimand as $f_y(\bm{x}_i, 1) - f_y(\bm{x}_i, 0) \approx \frac{1}{D} \sum_{d=1}^{D} f_y^{(d)}(\bm{x}_i, 1) - f_y^{(d)}(\bm{x}_i, 0) $. A nonlinear machine learning method naively trained only on selected observations targets the same estimand because \footnote{The naive method can be an ``S-learner'' with one trained model for all selected observations, with treatment included as a covariate, or a ``T-learner'' with two models separately trained on selected treatment and control observations \citep{kunzel2019metalearners}.}
$$ E[Y_i^*(1) | Z_i^* \ge 0, \bm{x}_i, \bm{w}_i ] - E[Y_i^*(0) | Z_i^* \ge 0, \bm{x}_i, \bm{w}_i ] = f_y (\bm{x}_i, 1) - f_y (\bm{x}_i, 0) + \gamma \Bigg( \frac{\phi(f_z(\bm{w}_i) )}{\Phi(f_z(\bm{w}_i) )} - \frac{\phi(f_z(\bm{w}_i) )}{\Phi(f_z(\bm{w}_i) )} \Bigg) $$
$$ = f_y (\bm{x}_i, 1) - f_y (\bm{x}_i, 0) . $$
Now suppose that treatment enters both the outcome and selection equations, so that we have $f_y(\bm{x}_i, T_i)$ and $f_z(\bm{w}_i, T_i)$, and assume $ Y_i^*(1), Y_i^*(0), Z_i^*(1), Z_i^*(0), \perp T_i | \bm{x}_i,\bm{w}_i$. TOBART-2 still directly models the estimand of interest $f_y (\bm{x}_i, 1) - f_y (\bm{x}_i, 0)$. However, the treatment effect estimates produced by a machine learning method naively trained only on the selected observations would target an estimand containing a selection bias term:
$$ E[Y_i^*(1) | Z_i^* (1) \ge 0, \bm{x}_i, \bm{w}_i ] - E[Y_i^*(0) | Z_i^*(0) \ge 0, \bm{x}_i, \bm{w}_i ] = $$
$$ f_y (\bm{x}_i, 1) - f_y (\bm{x}_i, 0) + \gamma \Bigg( \frac{\phi(f_z(\bm{w}_i,1) )}{\Phi(f_z(\bm{w}_i,1) )} - \frac{\phi(f_z(\bm{w}_i,0) )}{\Phi(f_z(\bm{w}_i,0) )} \Bigg) $$
Therefore TOBART-2 is expected to produce more accurate ITE estimates for DGPs in which treatment affects sample selection.
\bigskip
\section{Simulation Study}\label{simulation_sec}
\subsection{Description of Simulations}
In addition to the simulations described here, we provide a simulation study with the linear DGP of \cite{iqbal2023bayesian} in Appendix \ref{further_sims}, although this DGP does not include an excluded instrument.
\subsubsection{ \cite{brewer2024addressing} Simulations}
We simulate the data generating processes described by \cite{brewer2024addressing}, with less observations and variables for computational feasibility. We also consider non-normally distributed errors. Covariates $X_1,...,X_{10}$ are generated from a multivariate normal distribution with $\text{Cov}(X_{ki},X_{ji}) = 0.3^{|k-j|}$. The excluded instrument, $W_{exc,i}$, that appears in the selection equation, but not the outcome equation is generated by $W_{exc,i} = \sum_{j=1}^{10} 0.05 X_{ji} + e_i \ , \ e_{w,i} \sim \mathcal{N}(0,0.75) $.
We consider the following DGPs, and two more DGPs in Appendix \ref{further_sims}:
\begin{enumerate}
\item Outcome: $Y_i^* = \sum_{j=1}^{10} \frac{0.4}{j^2} X_{ji} + \eta_i$ \\
Selection: $ Z_{i}^* = 1.25 + \sum_{j=1}^{10} \frac{0.1}{(10.5 - j)^2} X_{jo} + W_{exc,i} + \xi_i $. The errors $\eta_i$ and $\xi_i$ are bivariate standard normal with correlation $\rho \in \{ 0,0.45,0.9 \} $.
\item Outcome: $Y_i^* = -0.25 1.25 \sin \left( \frac{\pi}{4} + 0.75 \pi \sum_{j=1}^{10} \frac{0.4}{j^2} X_{ji} \right) + \eta_i$ \\
Selection: $ Z_{i}^* = 1.25 + \sum_{j=1}^{10} \frac{0.1}{(10.5 - j)^2} X_{jo} + W_{exc,i} + \xi_i $. The errors $\eta_i$ and $\xi_i$ are bivariate standard normal with correlation $\rho \in \{ 0,0.45,0.9 \} $.
\end{enumerate}
Appendix \ref{further_sims} includes results for the above DGPs with the following non-normal distributions for the errors. \footnote{TOBART-2-NP results will be added to an updated version of this paper at a later date.}
\begin{itemize}
\item $\eta_i \sim \rho V_i + \rho t_{\nu = 5} $ and $\xi_i \sim t_{\nu = 5} $, where $t_{nu = 5}$ is a standard t-distribution with $5$ degrees of freedom. The correlation is $\rho = \frac{1}{sqrt{2}} $.
\item The errors are generated from a mixture of standard bivariate normals $ 0.3 \mathcal{MVN}_2 ( (0,-2.1)', \Sigma) + 0.7 \mathcal{MVN}_2 ( (0,0.9)', \Sigma) $ where $\Sigma_{12} = 0.85$ , and the implied correlation is $0.5$.
\end{itemize}
The numbers of training and test observations are set to $2500$ and $500$ respective;y, and the number of repetitions for each simulation scenario is set to $10$.
\subsubsection{Heterogeneous Treatment Effect Simulations}
When designing a simulation study to evaluate heterogeneous treatment effect estimation, we consider the following:
\begin{itemize}
\item There is unlikely to be a notable difference between TOBART-2 and a naive application of a nonparametric regression method to sample selected data unless treatment has a non-negligible effect on the inverse Mill's ratio for a substantial number of observations, assuming bivariate normally distributed errors. Furthermore, if the bias of naive application of nonparametric regression methods is to be meaningful, it must be of non-trivial magnitude relative to the true effect on the latent outcome.
\item To identify the true conditional mean of the latent outcome (i.e. without the inverse Mills ratio term) and identify heterogeneous effects, we require the excluded variables to have a large impact on the inverse Mill's ratio for a range of values of both treatment and other covariates.
\item We require sufficiently many treated, untreated, selected, and unselected observations.
\end{itemize}
We generate $p$ independent standard normally distributed covariates $X_{1},\dots, X_p \sim \mathcal{N}(0,1)$. A binary treatment variable is generated with selection on observables. The number of covariates is $p=50$, and $X_1, \dots, X_5$ are excluded from the outcome equation data, i.e. they are in $\bm{w}_i$, but not $\bm{x}_i$. There are 4000 training observations and 4000 test observations. The propensity score function is:
$$\pi(\bm{x}_i) = 0.7 \Phi(0.5 x_{6,i}x_{7,i} + 0.3 x_{8,i}^2) +0.1 $$
The treatment variable is denoted by $T_i \sim \text{Bernoulli}(\pi(\bm{x}_i))$. The selection and outcome equations are:
$$ Y_i^* = 0.5x_{6,i} + 0.5 x_{6,i}x_{7,i} + 0.25 (x_{3,i}-1)^2 T_i + \eta_i $$
$$ Z_{i}^* = 0.5 x_{1,i}x_{6,i} + x_{1,i}x_{6,i}T_i + \xi_i $$
This implies that the heterogeneous effects on the inverse Mills ratio are equal to:
$$ \frac{\phi(f_z(\bm{w}_i,1) )}{\Phi(f_z(\bm{w}_i,1) )} - \frac{\phi(f_z(\bm{w}_i,0) )}{\Phi(f_z(\bm{w}_i,0) )} = \frac{\phi (0.5 x_{1,i}x_{6,i} + x_{1,i}x_{6,i} ) }{ \Phi(0.5 x_{1,i}x_{6,i} + x_{1,i}x_{6,i} ) } - \frac{\phi (0.5 x_{1,i}x_{6,i} ) }{ \Phi(0.5 x_{1,i}x_{6,i} ) } $$
The heterogeneous effects on the latent outcome are equal to:
$$ f_y (\bm{x}_i, 1) - f_y (\bm{x}_i, 0) = 0.25 (x_{3,i}-1)^2 $$
Let the potential outcomes of the selection indicator variable under treatment and control allocation be denoted by $S_i(1)$ and $S_i(0)$ respectively. We use the following measures to evaluate all methods:
\medskip
\noindent $\bullet$ MSE of estimated effect on the probability of censoring
$$ \frac{1}{n_{test}} \sum_{i=1}^{n_{test}} \Big( \big(\hat{\Pr}(S_i(1)=0|X) - \hat{\Pr}(S_i(0)=0|X) \big) - \big({\Pr}(S_i(1)=0|X) - {\Pr}(S_i(0)=0 |X) \big)\Big)^2 $$
$\bullet$ RMSE of estimated treatment effect on the latent outcome:
$$ \sqrt{\frac{1}{n_{test}} \sum_{i=1}^{n_{test}} \Big( \big(\hat{\mathbb{E}}(Y_i(1)=0|X) - \hat{\mathbb{E}}(Y_i(0)=0|X) \big) - \big({\mathbb{E}}(Y_i(1)=0|X) - {\mathbb{E}}(Y_i(0)=0 |X) \big)\Big)^2 }$$
\subsection{Simulation Results}\label{sim_results}
\subsubsection{\cite{brewer2024addressing} Simulations }
Table \ref{Brewer_res_table} contains the results for the simulation study from \cite{brewer2024addressing}.\footnote{TOBART-2-NP results will be added to an updated version of this paper at a later date.} As anticipated, the sparse linear Bayesian model described by \cite{iqbal2023bayesian} produces better results for DGP 1, because a linear model is appropriate for a DGP with linear selection and outcome equations. BART and TOBART do not improve estimates of sample selection probabilities relative to the sparse linear model, because for both DGPs the true model for the latent selection variable is linear.
Latent outcome predictions for the nonlinear DGP2 are notably more accurate for Soft BART and TOBART models than for linear Tobit models or BART. When there is no sample selection, there is no notable difference between the SoftBART and Soft TOBART models. When there is sample selection in unobservables, TOBART models produce much more accurate nonlinear DGP predictions than the other models. Moreover, confidence interval coverage is close to 95\% and posterior means of the correlation between selection and outcome equation errors are very close to the true values.
The best performing TOBART models are marginalized TOBART-2 and Soft TOBART-2. A combination of the two models, i.e. marginalized soft TOBART-2, produces similar estimated to unmarginalized Soft TOBART-2, although it produces more accurate estimates of $\rho$ when $\rho=0.45$.\footnote{Marginalized SoftBART surprisingly produces much less accurate esitmates of selection probabilities. Perhaps this is explainable by differences in the implementations of unmarginalized and marginalized tree samplers.} We do not observe large differences in results for different priors on the distribution of the error terms, and different calibrations of these priors. Perhaps the methods would be more sensitive to the prior if the sample size were much smaller than 2000.
\begin{table}[ht]
\centering
\begin{tabular}{p{0.1cm}p{0.6cm}|p{1.08cm}p{1.08cm}p{1.08cm}p{1.0cm}p{1.0cm}p{1.08cm}p{1.08cm}p{1.08cm}p{1.08cm}p{1.08cm}p{1.08cm}p{1.08cm}}
\hline
D G P & corr & Tobit VH & Tobit Ding & Tobit Omori & BART & Soft BART & Sparse Tobit & TO BART 2 VH & TO BART 2 Ding & TO BART 2 marg & Soft TO BART 2 VH & Soft TO BART 2 marg \\
\hline
\hline
\hline
\multicolumn{12}{c}{$f_y(\bm{x})$ prediction RMSE, relative to BART}\\
\hline
1 & 0.00 & 0.677 & 0.666 & 0.676 & 1.000 & 0.554 & 0.365 & 0.854 & 0.865 & 0.850 & 0.576 & 0.705 \\
2 & 0.00 & 2.350 & 2.349 & 2.350 & 1.000 & 0.746 & 2.556 & 0.997 & 0.998 & 0.986 & 0.760 & 0.772 \\
\hline
1 & 0.45 & 0.396 & 0.440 & 0.395 & 1.000 & 0.792 & 0.212 & 0.496 & 0.557 & 0.487 & 0.361 & 0.396 \\
2 & 0.45 & 1.543 & 1.551 & 1.543 & 1.000 & 0.832 & 1.679 & 0.645 & 0.680 & 0.643 & 0.505 & 0.505 \\
\hline
1 & 0.90 & 0.240 & 0.272 & 0.239 & 1.000 & 0.897 & 0.123 & 0.257 & 0.306 & 0.262 & 0.190 & 0.207 \\
2 & 0.90 & 0.989 & 1.001 & 0.989 & 1.000 & 0.913 & 1.075 & 0.389 & 0.414 & 0.401 & 0.310 & 0.306 \\
\hline
\hline
\hline
\multicolumn{12}{c}{$f_y(\bm{x})$ 95\% prediction interval mean coverage}\\
\hline
1 & 0.00 & 0.871 & 0.868 & 0.870 & 0.985 & 0.982 & 0.956 & 0.953 & 0.949 & 0.972 & 0.978 & 0.994 \\
2 & 0.00 & 0.211 & 0.208 & 0.213 & 0.973 & 0.934 & 0.157 & 0.933 & 0.930 & 0.943 & 0.934 & 0.990 \\
\hline
1 & 0.45 & 0.853 & 0.781 & 0.854 & 0.885 & 0.684 & 0.946 & 0.954 & 0.914 & 0.972 & 0.966 & 0.989 \\
2 & 0.45 & 0.212 & 0.226 & 0.211 & 0.884 & 0.737 & 0.153 & 0.928 & 0.908 & 0.950 & 0.931 & 0.993 \\
\hline
1 & 0.90 & 0.809 & 0.728 & 0.811 & 0.647 & 0.511 & 0.895 & 0.947 & 0.909 & 0.985 & 0.952 & 0.993 \\
2 & 0.90 & 0.197 & 0.225 & 0.196 & 0.683 & 0.550 & 0.125 & 0.897 & 0.896 & 0.932 & 0.910 & 0.991 \\
\hline
\hline
\hline
\multicolumn{12}{c}{Selection Probability MSE, relative to BART}\\
\hline
1 & 0.00 & 2.154 & 2.151 & 2.153 & 1.000 & 0.727 & 0.249 & 1.409 & 1.454 & 1.808 & 0.735 & 1.173 \\
2 & 0.00 & 2.143 & 2.144 & 2.155 & 1.000 & 0.727 & 0.252 & 1.387 & 1.423 & 1.775 & 0.716 & 1.162 \\
\hline
1 & 0.45 & 2.115 & 2.139 & 2.122 & 1.000 & 0.724 & 0.196 & 1.556 & 1.565 & 1.533 & 0.686 & 0.892 \\
2 & 0.45 & 2.128 & 2.143 & 2.134 & 1.000 & 0.724 & 0.213 & 1.539 & 1.589 & 1.667 & 0.687 & 0.908 \\
\hline
1 & 0.90 & 2.158 & 2.258 & 2.155 & 1.000 & 0.625 & 0.232 & 1.370 & 1.422 & 1.684 & 0.566 & 1.643 \\
2 & 0.90 & 2.200 & 2.274 & 2.200 & 1.000 & 0.625 & 0.248 & 1.380 & 1.464 & 1.680 & 0.603 & 0.963 \\
\hline
\hline
\hline
\multicolumn{12}{c}{Mean correlation estimate for selection and outcome equation errors}\\
\hline
1 & 0.00 & -0.001 & 0.000 & -0.001 & 0.000 & 0.000 & -0.018 & -0.040 & -0.024 & -0.004 & -0.016 & 0.002 \\
2 & 0.00 & -0.004 & -0.001 & -0.003 & 0.000 & 0.000 & -0.035 & -0.032 & -0.019 & -0.009 & -0.023 & -0.005 \\
\hline
1 & 0.45 & 0.418 & 0.242 & 0.418 & 0.000 & 0.000 & 0.416 & 0.410 & 0.240 & 0.454 & 0.427 & 0.446 \\
2 & 0.45 & 0.373 & 0.219 & 0.376 & 0.000 & 0.000 & 0.350 & 0.413 & 0.244 & 0.461 & 0.427 & 0.448 \\
\hline
1 & 0.90 & 0.894 & 0.703 & 0.894 & 0.000 & 0.000 & 0.891 & 0.894 & 0.692 & 0.933 & 0.911 & 0.920 \\
2 & 0.90 & 0.751 & 0.542 & 0.751 & 0.000 & 0.000 & 0.709 & 0.884 & 0.682 & 0.923 & 0.900 & 0.913 \\
\hline
\end{tabular}
\caption{Results for \cite{brewer2024addressing} simulation study. RMSE and MSE results are relative to BART.} \label{Brewer_res_table}
\end{table}
\subsubsection{CATE simulation results}
Table \ref{cate_res_table} contains the results for the treatment effect simulation study. Nonlinear methods outperform linear methods for this simple quadratic DGP with interactions. Tree-based methods with soft splitting are more accurate than their counterparts with hard splitting rules. Marginalized TOBART outperforms vanilla implementations of TOBART, suggesting that without marginalization TOBART sampler did not converge or mix well.
Surprisingly, SoftBART with the inverse Mills ratio as a covariate, following the general approach of \cite{brewer2024addressing}, is almost as accurate as Soft TOBART-2 across all measures except PEHE of effects on selection and outcome CATE posterior interval length. Soft marginalized TOBART-2 produces the best estimates of outcome treatment effects, and surprisingly the worst predictions of selection treatment effects. It also produces the smallest posterior intervals, albeit with over-coverage.
\begin{table}[ht]
\centering
\begin{tabular}{p{0.6cm}p{0.4cm}p{0.7cm}|p{0.9cm}p{1cm}p{1.3cm}p{1.3cm}p{1cm}p{1cm}p{1cm}p{1cm}p{1cm}p{1cm}}
\hline
n & p & corr & Tobit VH & Sparse Bayes Tobit & BART control & Soft BART control & TO BART 2 marg & TO BART 2 VH & TO BART 2 Ding & Soft TO BART 2 VH & Soft TO BART 2 Ding & Soft TO BART 2 marg \\
\hline
& & & \multicolumn{10}{l}{ { \blue
\textit{PEHE Selection Test}} } \\
4000 & 50 & 0.50 & 1.000 & 0.997 & 0.703 & 0.646 & 0.544 & 0.655 & 0.626 & 0.593 & 0.409 & 1.738 \\
& & & \multicolumn{10}{|l}{ { \blue
\textit{PEHE Selection Training}} } \\
4000 & 50 & 0.50 & 1.000 & 0.997 & 0.710 & 0.646 & 0.547 & 0.656 & 0.630 & 0.594 & 0.408 & 1.739 \\
& & & \multicolumn{10}{|l}{{ \blue \textit{PEHE Outcome Test}}}\\
4000 & 50 & 0.50 & 1.000 & 1.013 & 0.403 & 0.295 & 0.345 & 0.451 & 0.431 & 0.295 & 0.297 & 0.270 \\
& & & \multicolumn{10}{|l} { { \blue \textit{PEHE Outcome Training}} }\\
4000 & 50 & 0.50 & 1.000 & 1.013 & 0.402 & 0.291 & 0.341 & 0.448 & 0.428 & 0.291 & 0.293 & 0.265 \\
& & & \multicolumn{10}{|l} { {\blue \textit{Length Outcome CATE PI Test}}}\\
4000 & 50 & 0.50 & 0.217 & 0.247 & 0.733 & 0.678 & 0.624 & 0.461 & 0.507 & 0.602 & 0.611 & 0.582 \\
& & & \multicolumn{10}{|l}{ { \blue \textit{Coverage Outcome CATE PI Test}} }\\
4000 & 50 & 0.50 & 0.144 & 0.226 & 0.938 & 0.978 & 0.932 & 0.790 & 0.836 & 0.974 & 0.973 & 0.979 \\
& & & \multicolumn{10}{|l}{ { \blue \textit{Length Outcome CATE PI Training}}}\\
4000 & 50 & 0.50 & 0.217 & 0.247 & 0.733 & 0.682 & 0.623 & 0.462 & 0.509 & 0.604 & 0.613 & 0.582 \\
& & & \multicolumn{10}{|l}{ { \blue \textit{Coverage Outcome CATE PI Training}}}\\
4000 & 50 & 0.50 & 0.145 & 0.228 & 0.938 & 0.980 & 0.932 & 0.788 & 0.839 & 0.977 & 0.976 & 0.981 \\
\hline
\end{tabular}
\caption{Results for Treatment Effect simulation study. PEHE results are relative to linear Bayesian Tobit. PI = 95 \% Posterior Interval} \label{cate_res_table}
\end{table}
\FloatBarrier
\section{Data Application}\label{application_sec}
\subsubsection{Description of RAND Health Insurance Experiment}
We estimate the heterogeneous effects of allocation to a deductible plan for health insurance. The dataset has been studied by \cite{deb2002structure}, \cite{cameron2005microeconometrics}, and others. The outcome variable is the log of medical expenses in dollars. Selected observations are those with non-zero medical expenditures. There are 15 other covariates.\footnote{The other covariates are: the log of one plus the coinsurance rate, the log of the participation incentive payment, physical limitations (binary), number of chronic diseases, good/fair/poor self-rated health (binary variables), log of family income, log of family size, education of household head in years, age in years, female (binary), child (1 if has child below 18 years), interaction of female and child variables, black (binary).} None of the covariates are excluded from the outcome equation. Further applications to female earning data \citep{mroz1987sensitivity} and a job training experiment for disadvantaged youths \citep{schochet2001national} are included in Appendix \ref{further_application_app}.
\subsubsection{Data Application Results}
The data application results are presented below.\footnote{TOBART-2-NP results will be added to an updated version of this paper at a later date.} Table \ref{mroz_dependency_table} shows that the posterior mean and 95\% posterior interval for the correlation parameter $\rho$ are comparable to the estimates produced by \cite{van2011bayesian}. Table \ref{mroz_ate_table} contains estimates of average treatment effects of individual deductible plans on probabilities of zero expenditures and the log of medical expenditures. The average treatment effect estimates are larger in magnitude for the Bayesian methods than for standard Tobit. The average treatment effect on the latent outcome estimate produced by Soft-TOBART-2 is remarkably close to the effect estimated by \cite{van2011bayesian} with a Bayesian linear Tobit model.
Appendix \ref{data_app} contains 5-fold cross validated Brier scores of selection probabilities, MSE of observed outcome predictions, and coverage of observed outcome posterior intervals. All models produce similar results. This is unsurprising because sample selection models are not expected to produce better predictions of the observed outcomes, but rather better predictions of the mean of latent outcomes, which cannot be evaluated using observed sample selected data.
\FloatBarrier
\begin{table}[ht]
\centering
\begin{tabular}{llll}
\hline
\hline
& Mean & 2.5\% Quantile & 97.5\% Quantile \\
\hline
TOBART-2 & 0.679 & 0.519 & 0.763 \\
\cite{van2011bayesian} & 0.712 & 0.634 & 0.784 \\
\hline
\end{tabular}
\caption{Mean and posterior distributions of dependence measures for selection and outcome equation error terms of TOBART model applied to RAND Health Insurance Experiment data, and original results from \cite{van2011bayesian}.}
\label{mroz_dependency_table}
\end{table}
\begin{table}[H]
\centering
\begin{tabular}{lp{1.25cm}p{1.15cm}p{1.15cm}p{1.15cm}p{1.25cm}p{1.5cm}p{1.25cm}p{1.5cm}p{1.5cm}}
\hline
\hline
& Tobit 2-step & Tobit ML & BART & Soft BART & Bayesian Splines & TO BART 2 & Soft \newline TO BART 2 & Bayes Tobit \cite{van2011bayesian} \\
\hline
Prob. Censored & -0.039 & -0.031 & -0.063 & -0.085 & -0.031 & -0.079 & -0.04 & \\
Latent Outcome & -0.042 & -0.108 & & & -0.13 & -0.133 & -0.147 & -0.147 \\
Observed Outcome & -0.025 & -0.053 & -0.13 & -0.011 & & -0.157 & -0.084 & \\
\hline
\end{tabular}
\caption{Estimates of average treatment effects of individual deductible plans on the probability of censoring (i.e. zero expenditure), the the latent outcome, and observed log of medical expenses (conditional on expenses being positive) for RAND Health Insurance Experiment data.}
\label{mroz_ate_table}
\end{table}
\FloatBarrier
\FloatBarrier
\section{Conclusion}\label{conclusion_sec}
The results of this paper suggest that TOBART-2 can produce more accurate outcome predictions and treatment effect estimates than linear sample selection models when the DGP is non-linear.
However, for some DGPs, a naive application of BART or softBART to the selected data, with an estimate of the inverse Mills ratio as a covariate, also provides estimates with comparable accuracy.
This paper also describes alternative MCMC implementations and prior settings. Improvements in predictive accuracy can be obtained by the use of soft splitting rules and marginalization in MCMC draws of trees. Prior sensitivity, particularly for small sample sizes, would be an interesting area for future research. Possible extensions of the model include heteroskedastic type 2 Tobit models \citep{donald1995two, de2022generalized}, with modelling of variances by Bayesian trees \citep{pratola2020heteroscedastic}.
\bibliographystyle{agsm}
\bibliography{tobart_references}