Research ArticleOpen Access

1]Youfei Yu

2]Jiacong Du

[3]Min Zhang

4]Zhenke Wu

5]Andrew M. Ryan

6]Bhramar Mukherjee

[1]Department of Biostatistics, University of Michigan; E-mail<EMAIL_ADDRESS>
[2]Department of Biostatistics, University of Michigan; E-mail<EMAIL_ADDRESS>
[3]Department of Biostatistics, University of Michigan; E-mail<EMAIL_ADDRESS>
[4]Department of Biostatistics, University of Michigan; E-mail<EMAIL_ADDRESS>
[5]Department of Health Management and Policy, University of Michigan; E-mail<EMAIL_ADDRESS>
[6]Department of Biostatistics, University of Michigan; E-mail<EMAIL_ADDRESS>
Outcome Adaptive Propensity Score Methods for Handling Censoring and High-Dimensionality: Application to Insurance Claims

Outcome-adaptive propensity score methods

Propensity scores are commonly used to reduce the confounding bias in non-randomized observational studies for estimating the average treatment effect. An important assumption underlying this approach is that all confounders that are associated with both the treatment and the outcome of interest are measured and included in the propensity score model. In the absence of strong prior knowledge about potential confounders, researchers may agnostically want to adjust for a high-dimensional set of pre-treatment variables. As such, variable selection procedure is needed for propensity score estimation. In addition, recent studies show that including variables related to treatment only in the propensity score model may inflate the variance of the treatment effect estimates, while including variables that are predictive of only the outcome can improve efficiency.
In this paper, we propose a flexible approach to incorporating outcome-covariate relationship in the propensity score model by including the predicted binary outcome probability (OP) as a covariate. Our approach can be easily adapted to an ensemble of variable selection methods, including regularization methods and modern machine learning tools based on classification and regression trees.
We evaluate our method to estimate the treatment effects on a binary outcome, which is possibly censored, among multiple treatment groups.
Simulation studies indicate that incorporating OP for estimating the propensity scores can improve statistical efficiency and protect against model misspecification. The proposed methods are applied to a cohort of advanced stage prostate cancer patients identified from a private insurance claims database for comparing the adverse effects of four commonly used drugs for treating castration-resistant prostate cancer.

Journal of Causal Inference

§ INTRODUCTION

To obtain an unbiased estimate for the causal average treatment effect (ATE) using data from observational studies, it is important to control for confounding by adjusting for the differences in pre-treatment baseline covariates between treatment groups. A common tool for confounder adjustment is the propensity score, defined as the conditional probability of receiving a specific treatment given a set of covariates [Rosenbaum and Rubin, 1983], which is usually unknown in observational studies and needs to be estimated from the data. The estimated propensity scores can be used in matching [Stuart, 2010, Yang et al., 2016], weighting [Lunceford and Davidian, 2004], regression adjustment [Zhou et al., 2019], among other methods, for estimating the causal treatment effect. Propensity score-based causal estimators rely on robust and accurate estimation for the propensity scores.

There has been a large body of work focusing on the estimation and use of propensity scores in a low-dimensional setting, where the sample size is far greater than the number of candidate covariates. In this setting, empirical researchers often estimate the propensity scores by fitting a logistic regression model for binary treatment, or a multinomial logistic regression model for more than two treatment groups. In the era of `big data', large healthcare databases are increasingly being used to conduct comparative effectiveness research. For example, insurance claims data were used to compare the safety of four drugs on the market prescribed for patients with metastatic castration-resistant prostate cancer [Yu et al., 2021]. A massive collection of candidate covariates, such as demographics, socioeconomic status, clinical measurements, and diagnosis codes, can be ascertained from patients' health care or insurance claims data, and the number of covariates may exceed the sample size in each treatment group. Standard parametric regression models may become problematic when the dimensionality increases, and a key challenge is to identify variables to be included in the propensity score model from a high-dimensional set of measured covariates. There has been considerable interest in developing methods that perform variable selection for estimating propensity/balancing scores in the high-dimensional setting. For instance, Schneeweiss et al., 2009 proposed the high-dimensional propensity score (hd-PS) algorithm that selects covariates to facilitate high-dimensional propensity score adjustment using health care claims data. Specifically, they rank each covariate based on its potential for controlling confounding by assessing the covariate's mean and univariate association with the treatment and outcome, and then pick the top $k$ covariates for inclusion in the propensity score modeling.

Parametric models require variable selection and specification of functional forms. A practical difficulty is that misspecification of the model can result in substantial bias of the estimated treatment effect [Kang and Schafer, 2007]. Machine learning methods are well-known for their powerful predictive performance and ability to handle complex and nonlinear relationship between outcomes and covariates. There has been a growing interest in using machine learning techniques for propensity score modeling [Setoguchi et al., 2008, Lee et al., 2010, Ju et al., 2019, Tu, 2019]. Setoguchi et al., 2008 compared several data mining techniques that optimize the prediction of treatment status, including classification and regression trees (CART), pruned CART, and neural networks, in the context of propensity score matching with a continuous outcome. Lee et al., 2010 extended the work of Setoguchi et al., 2008 to the setting with a binary outcome, and evaluated the performance of CART, pruned CART, bootstrap aggregated (bagged) CART, random forests, and boosted CART with regard to propensity score weighting. Both works focused on the low-dimensional setting, and demonstrated that machine learning methods, such as CART and neural networks, are promising alternatives to parametric modeling for the estimation of propensity scores in the presence of nonadditivity and/or nonlinearity in the true treatment assignment mechanism. However, numerical evaluations of these machine learning techniques for high-dimensional covariates and multiple treatment groups with a binary endpoint remain limited.

These aforementioned approaches only consider the treatment-covariate relationship when modeling the propensity/balancing scores, and fail to incorporate the outcome models into the treatment modeling process. Studies illustrate that using covariates that are associated with the treatment but not the outcome will inflate the variance of the estimators of ATE without reducing the bias [De Luna et al., 2011, Shortreed and Ertefaie, 2017]. On the other hand, adding covariates explaining the outcome but not the treatment to the propensity score model can improve the precision of the treatment effect estimates [Brookhart et al., 2006, Shortreed and Ertefaie, 2017]. These findings suggest that a highly predictive model for treatment assignment will not necessarily lead to the most efficient estimators of treatment effects. Therefore, standard variable selection methods designed for prediction, which rely only on the relationship between treatment and covariates, may yield sub-optimal results in the context of causal inference. There is an expanding literature on variable selection methods for causal inference that account for the information about the outcome-covariate relationship. Shortreed and Ertefaie, 2017 proposed the Outcome-Adaptive LASSO (OAL) method, which adopts the adaptive LASSO framework [Zou, 2006] and places smaller adaptive weights on the outcome predictors in the propensity score model to favor the inclusion of outcome-predictive covariates. As a result, heavier penalties are imposed on variables relevant to treatment only than variables predictive of outcome only. Other works that allow the outcome information to contribute to variable selection for propensity score modeling include Outcome Highly Adaptive LASSO proposed by Ju et al., 2020 and Bayesian Adjustment for Confounding proposed by Zigler and Dominici, 2014. However, how the information about the outcome-covariate relationship can be incorporated into tree-based machine learning methods, such as CART, and to what extent can outcome models contribute to the efficiency gain in high-dimensional data settings with multiple treatment groups have been less studied.

In this paper, we propose a flexible approach to incorporating outcome-covariate relationship in the propensity score model by including the predicted outcome probability (OP) as a covariate. Our approach can be easily adapted to an ensemble of variable selection methods, including regularization methods and tree-based machine learning algorithms. We evaluate our approach for estimating the causal average treatment effect with multiple treatment groups and high-dimensional covariates. We use the inverse probability weighting (IPW) estimator, although the estimated propensity scores considered in this study are applicable to any propensity score-based methodology, such as propensity score matching. We also aim to account for censoring at the same time, where the bias due to censoring is controlled for by applying the inverse probability of remaining uncensored as weights to the outcome. We focus on estimating the treatment effects on a binary outcome (that is possibly censored) among multiple treatment groups. To the best of our knowledge, the literature lying within the intersection of high-dimensionality, complexity of the associations between treatment and covariates, multiple treatments, censored outcomes, and outcome-adaptive propensity score modeling is largely absent.

In Section <ref>, we introduce the notation and basic setup of the problem. Section <ref> describes the algorithm for selecting variables to be included in the treatment model. Section <ref> outlines the methods considered for the final treatment model that is used to estimates the propensity scores. Section <ref> presents simulation studies that compare the methods in settings with high-dimensional covariates and potentially complex underlying treatment model (i.e., model with nonlinearity and/or nonadditivity). We also consider a setting where censoring exists. The simulation results illustrate the efficiency gain provided by leveraging the information about the outcome-covariate relationship when estimating the propensity scores. Section <ref> presents an example that compares the risks of possibly right-censored hospitalization and emergency room (ER) visits of four prostate cancer treatments using data from an insurance claims database. The binary outcomes (i.e., whether the patient was hospitalized/admitted to ER within 180 or 360 days of treatment initiation) in this data example were subject to censoring, as the follow-up period tended to terminate early due to drop-out or treatment switch. We close with a brief discussion.

§ DEFINITION OF THE PROBLEM AND NOTATION

§.§ Notation and Assumptions

We consider $n$ independent individuals, indexed by $i$, with $\boldsymbol{X}_i$ being a $p$-dimensional vector of covariates measured prior to receiving the treatment $Z_i$, where $Z_i=j\in\{1,\cdots,J\}$. We let $\mathcal{C}$, $\mathcal{Z}$, $\mathcal{Y}$, and $\mathcal{S}$ be the indices of confounders (i.e., covariates associated with both treatment and outcome), covariates predictive of treatment only, covariates predictive of outcome only, and covariates unrelated to both treatment and outcome (i.e., spurious covariates), respectively. Suppressing the index by $i$, $\boldsymbol{X}_\mathcal{C}$, $\boldsymbol{X}_\mathcal{Z}$, $\boldsymbol{X}_\mathcal{Y}$, and $\boldsymbol{X}_\mathcal{S}$ are mutually exclusive and $\boldsymbol{X}=\boldsymbol{X}_\mathcal{C}\cup\boldsymbol{X}_\mathcal{Z}\cup\boldsymbol{X}_\mathcal{Y}\cup\boldsymbol{X}_\mathcal{S}$. We further let $|\mathcal{C}|$, $|\mathcal{Z}|$, $|\mathcal{Y}|$, and $|\mathcal{S}|$ denote the cardinality of the corresponding set. Re-introducing subject index $i$, we let $T_i$ be the underlying time to the first event of interest for each individual, and $C_i$ be the censoring time. The outcome of interest is whether the event of interest occurs before a prespecified time point $d$, defined by $Y_i=I(T_i<d)$, which results in a possibly censored binary outcome. The information on $Y_i$ may not be completely available due to dropout, study termination, or treatment switch. In the absence of censoring, $T_i$ (and therefore $Y_i$) would be observed for all individuals, and the set of complete data is then $(\boldsymbol{X}_i,Z_i,Y_i)$. When the outcome variable is subject to right-censoring, $Y_i$ is observed only if the individual has not been censored before $d$, and we let $R_i=I\{C_i\ge\min(T_i,d)\}$ be the indicator of observing $Y_i$.

Under the potential outcome framework [Rubin, 1974], each individual is associated with a set of potential outcomes $\{Y^{(1)},\cdots,Y^{(J)}\}$, where $Y^{(j)}$ denotes the potential outcome had the individual received treatment $j$. The causal parameter of interest is the marginal ATE on the outcome between $j$ and $j'$, denoted by $\tau(j,j')=E\{Y^{(j')} \}-E\{Y^{(j)}\}$. The following assumptions are required for statistical identification of causal effects using the observable data [Yu et al., 2022]:

* (Random sampling) The individuals in the study are randomly sampled from the population;

* (Stable Unit Treatment Value Assumption, or SUTVA) For any individual $i$, $i=1,\cdots,n$, if $Z_i=j$, then $Y_i=Y_i^{(j)}$, for all $j=1,\cdots,J$;

* (Unconfoundedness) $\{Y_i^{(1)},\cdots, Y_i^{(J)}\}\condind Z_i|\boldsymbol{X}_i$;

* (Overlap) For all values of $j$ and $\boldsymbol{x}$, $0<\pi_j(\boldsymbol{x})<1$, where $\pi_j(\boldsymbol{x})=pr(Z_i=j|\boldsymbol{x})$;

* (Censoring at random) $ C_i\condind\{T_i^{(1)},\cdots,T_i^{(J)}\}\Big|(Z_i, \boldsymbol{X}_i)$.

§.§ Underlying Models for Outcome, Treatment, and Censoring, and Estimators for Average Treatment Effects

We assume that the true outcome model for $Z_i=j$ is a logistic regression model,

\begin{align*}
\text{logit}\,P(Y_i=1|\boldsymbol{W}_i, Z_i=j)=\boldsymbol{W}_i^T\boldsymbol{\beta}_j,
\end{align*}

where $\boldsymbol{W}_i$ is a $p_w$-dimensional function of $(\boldsymbol{X}_\mathcal{C}^T, \boldsymbol{X}_\mathcal{Y}^T)^T$. The treatment assignment mechanism is governed by a multinomial logistic regression

\begin{align*}
\log\frac{P(Z_i=j|\boldsymbol{V}_i)}{P(Z_i=J|\boldsymbol{V}_i)}=\boldsymbol{V}_i^T\boldsymbol{\alpha}_j,\quad j=1,\cdots,J,
\end{align*}

where $J$ is the reference level and $\boldsymbol{V}_i$ is a $p_v$-dimensional function of $(\boldsymbol{X}_\mathcal{C}^T, \boldsymbol{X}_\mathcal{Z}^T)^T$. Both $\boldsymbol{W}_i$ and $\boldsymbol{V}_i$ may contain nonlinear terms and interactions. We assume that $|\mathcal{C}|+|\mathcal{Y}|\ll p$ and $|\mathcal{C}|+|\mathcal{Z}|\ll p$, where $|\mathcal{C}|$, $|\mathcal{Y}|$, and $|\mathcal{Z}|$ are the numbers of variables in $\boldsymbol{X}_\mathcal{C}$, $\boldsymbol{X}_\mathcal{Z}$, and $\boldsymbol{X}_\mathcal{Y}$, respectively. Note that in practice, $\boldsymbol{X}_\mathcal{C}$, $\boldsymbol{X}_\mathcal{Y}$, and $\boldsymbol{X}_\mathcal{Z}$ are generally unknown and a common practice is to include all candidate covariates in the model. In the case where the ratio of $p$ to $n$ is relatively large, traditional regression models based on maximum likelihood estimation may not converge, and some variable selection procedure is required for model fitting.

With respect to censoring, we assume a proportional hazard model for treatment $j=1,\cdots,J$,

\begin{align*}
\lambda_j(t|\boldsymbol{U}_i)=\lambda_{0j}(t)\exp(\boldsymbol{U}_i^T\boldsymbol{\gamma}_j) ,
\end{align*}

where $\lambda_{0j}(t)$ is the treatment-specific baseline hazard function and $\boldsymbol{U}_i$ is a $p_u$-dimensional function of $\boldsymbol{X}_i$. The censoring model is assumed to have a moderate number of predictors ($p_u\ll n$) that are known a priori.

We estimate the ATE using the IPW estimator,

\begin{align*}
\hat{\tau}(j,j')=\frac{\sum_{i=1}^n\hat{w}_{i}I(Z_i=j')Y_i}{\sum_{i=1}^n\hat{w}_{i}I(Z_i=j')}\,-\, \frac{\sum_{i=1}^n\hat{w}_{i}I(Z_i=j)Y_i}{\sum_{i=1}^n\hat{w}_{i}I(Z_i=j)} ,
\end{align*}

where $\hat{w}_i=\sum_{j=1}^J 1/\hat{\pi}_j(\boldsymbol{X}_i)$ and the propensity scores $\hat{\pi}_j(\boldsymbol{X}_i)$, $j=1\cdots,J$, are required to be estimated from the data in an observational study. In the presence of censoring, the outcome is further weighted by the inverse probability of remaining uncensored at $d$,

\begin{align*}
\hat{\tau}(j,j')=\frac{\sum_{i=1}^n\hat{w}_i^{\ast}R_iI(Z_i=j')Y_i}{\sum_{i=1}^n\hat{w}_i^{\ast}R_iI(Z_i=j')}\,-\, \frac{\sum_{i=1}^n\hat{w}_i^{\ast}R_iI(Z_i=j)Y_i}{\sum_{i=1}^n\hat{w}_i^{\ast}R_iI(Z_i=j)}.
\end{align*}

The weights $\hat{w}_i^{\ast}=\sum_{j=1}^J(\hat{\pi}_j(\boldsymbol{X}_i)\exp{\{\Lambda_{ij}(\min(T_i,C_i,d))\}})^{-1}$ [Yu et al., 2022], where $\Lambda_{ij}(t)$ is the cumulative hazard function of $C_i$ at $t$ for treatment $j$.

§ VARIABLE SELECTION FOR DIMENSIONALITY REDUCTION FOR THE PROPENSITY SCORE MODEL

When the dimension of the covariate vector is high, it tends to be infeasible to fit an unrestricted parametric model, such as a multinomial logistic regression model, for the treatment using all the available covariates. A practical problem for empirical researchers is to identify a subset of covariates to be conditioned on to control for confounding. Typically, in medical research, a list of important covariates will be suggested based on the evidence in the literature and/or expert opinion. However, as the number of available covariates increases, it becomes extremely difficult for human experts to check manually which variables are potential confounders. Alternatively, one can turn to data-driven variable selection approaches, such as LASSO [Tibshirani, 1996], which automatically select important variables for treatment predictions from all the available covariates. Figure <ref> displays a flowchart of several possible routes that can be followed to identify the set of covariates to be included in the final treatment model for estimating the propensity scores. A commonly chosen route is to apply shrinkage methods directly to the original reservoir of covariates (route 0.851 in Figure <ref>), and we label the set of covariates . In this case, variable selection and propensity score estimation are conducted simultaneously in a single step. However, as we noted later in our simulation, following route 0.851 can result in substantially biased effect estimates. For sparse high-dimensional data, a large value of tuning parameter is necessary to select a parsimonious model. However, large penalties at the same time increase the shrinkage of non-zero components, leading to biased estimation [Meinshausen, 2007]. Therefore, reducing the number of spurious covariates entering the final treatment model through variable selection may help improve the performance of the propensity score estimation methods. In addition, such traditional variable selection approaches targeting $\boldsymbol{X}_\mathcal{C}\cup\boldsymbol{X}_\mathcal{Z}$ for the treatment only are sub-optimal. It has been shown that the use of $\boldsymbol{X}_\mathcal{Z}$ for propensity score modeling may inflate the variance of the estimated ATE [Brookhart et al., 2006, Shortreed and Ertefaie, 2017]. On the other hand, including $\boldsymbol{X}_\mathcal{Y}$ predictive of the outcome in the propensity score model can improve the precision of the ATE estimates [Brookhart et al., 2006, Shortreed and Ertefaie, 2017]. Therefore, an optimal propensity model for estimating ATE should include $\boldsymbol{X}_\mathcal{C}$ and $\boldsymbol{X}_\mathcal{Y}$ simultaneously while excluding $\boldsymbol{X}_\mathcal{Z}$.

Flowchart for variable selection and propensity score estimation. Routes 0.851-0.856 correspond to different sets of input variables: 0.851 All, 0.852 Ysel, 0.853 YZsel, 0.854 OP+All, 0.855 OP+Ysel, 0.856 OP+YZsel.

§.§ Using Outcome Model for Variable Selection

To identify $\boldsymbol{X}_\mathcal{C}$, $\boldsymbol{X}_\mathcal{Y}$ and $\boldsymbol{X}_\mathcal{Z}$, we apply a pre-selection procedure to the original set of covariates (), where a regularized model (in this case we choose LASSO) is fitted separately for both the outcome and the treatment. Specifically, the model for a binary outcome is specified as

\begin{align}
\text{logit}P(Y_i=1|\boldsymbol{X}_i)=\boldsymbol{X}_i^T\boldsymbol{\theta},
\label{eq:Chap3_outcome_model}
\end{align}

the coefficients of which are estimated based on the LASSO penalty

\begin{align*}
\hat{\boldsymbol{\theta}} = \argmin_{\boldsymbol{\theta}}\sum_{i=1}^n\left[\{-Y_i\boldsymbol{X}_i^T\boldsymbol{\theta}+\ln(1+e^{\boldsymbol{X}_i^T\boldsymbol{\theta}})\} + \lambda\sum_{k=1}^p |\theta_k|\right],
\end{align*}

where the tuning parameter $\lambda$ is chosen using cross-validation. Here we exclude the treatment variable $Z$ from the outcome model, and only consider the association of $\boldsymbol{X}$ with $Y$ across all treatment groups. For the treatment assignment mechanism, we assume a multinomial logistic regression model,

\begin{align}
P(Z_i=j|\boldsymbol{X}_i) = \frac{\exp(\boldsymbol{X}_i^T\boldsymbol{\psi}_j)}{\sum_{l=1}^J\exp(\boldsymbol{X}_i^T\boldsymbol{\psi}_l)},
\label{eq:Chap3_trt_model}
\end{align}

which is fitted by minimizing the negative penalized log-likelihood

\begin{align*}
N^{-1}\sum_{i=1}^N\log P(Z_i=j|\boldsymbol{X}_i)
+ \lambda \sum_{k=1}^p \|\boldsymbol{\psi}_{\cdot k} \|_2 ,
\end{align*}

where $\boldsymbol{\psi}_{\cdot k}$ is a $J$-dimensional vector containing the $k$-th coefficient for all $J$ treatment groups [Friedman et al., 2010]. This group LASSO penalty function ensures that each variable in $\boldsymbol{X}_i$ is selected or excluded for all $J$ levels, as opposed to each level having its own set of selected variables. We use  to label the set of variables selected by Model (<ref>) and  to label the intersection of the two sets of variables selected by Model (<ref>) and Model (<ref>). In an ideal case,  is identical to $\boldsymbol{X}_\mathcal{C}$ and \ is identical to $\boldsymbol{X}_\mathcal{Y}$. Theoretically, adjusting for $\boldsymbol{X}_\mathcal{C}$ alone in the treatment model is sufficient to remove the confounding bias. Route 0.853 in Figure <ref> corresponds to the case where  is used as predictors for propensity scores estimation.

As we mentioned previously, the ideal propensity score model adjusts for $\boldsymbol{X}_\mathcal{C}\cup\boldsymbol{X}_\mathcal{Y}$ while excluding $\boldsymbol{X}_\mathcal{Z}$. The OAL method proposed by Shortreed and Ertefaie, 2017 intends to achieve this by fitting an adaptive LASSO for the treatment [Zou, 2006], where the adaptive weights are computed based on the coefficients from the outcome model, with smaller coefficients corresponding to larger weights. OAL discourages $\boldsymbol{X}_\mathcal{Z}$ from being selected and encourages the inclusion of $\boldsymbol{X}_\mathcal{Y}$ by imposing heavier penalties on $\boldsymbol{X}_\mathcal{Z}$ than $\boldsymbol{X}_\mathcal{Y}$. There are two limitations of the work of Shortreed and Ertefaie, 2017. The first is that the coefficients of the outcome model were estimated using unpenalized linear regression, which may be fitted for a continuous outcome in the context where the ratio of $p$ to $n$ is relatively large. However, for a binary outcome as considered in our study, it tends to be more difficult for a standard logistic regression to converge. We extend their method by considering a LASSO-fitted outcome model to compute the adaptive weights, in which case the coefficients of covariates not predictive of the outcome can be zero, and therefore the covariates with zero coefficients will be excluded from the treatment model. Another limitation is that the idea of OAL cannot be conveniently extended to machine learning methods that do not rely on regularization.

OAL incorporates the outcome-covariate associations as adaptive weights when selecting variables for propensity score estimation. An alternative is to directly include the covariates predictive of the outcome in the treatment model, namely, . The treatment can then be modeled as a function of  (route 0.852 in Figure <ref>).

One possible problem of using  as input for the final treatment model is that the dimensionality of  can still remain relatively high for sample size $n$. To improve propensity score modeling by utilizing the variables predictive of the outcome only while maintaining a reasonable dimensionality, we propose to incorporate the information contained in $\boldsymbol{X}_\mathcal{Y}$ into the propensity score estimation process in the form of predicted probabilities of the outcome, Outcome Probability (OP) for short, as opposed to directly including those variables in their original form. We estimate OP for each subject by $\hat{p}_i=\exp(\boldsymbol{X}_i^T\hat{\boldsymbol{\theta}})/\{1+\exp(\boldsymbol{X}_i^T\hat{\boldsymbol{\theta}})\}$, and the logit scale of $\hat{p}_i$ is denoted $\hat{p}_i^\ast=\log\{\hat{p}_i/(1-\hat{p}_i)\}$. OP summarizes the information about the outcome-covariate relationship and reduces the dimensionality of the outcome predictors to a one-dimensional vector. After obtaining  which targets $\boldsymbol{X}_\mathcal{C}$, we add OP back to the set of input variables to capture the information in $\boldsymbol{X}_\mathcal{Y}$ (route 0.856 in Figure <ref>, and we denote the combination ). In addition to , as was noted later in the simulation studies, OP can also be combined with all the available covariates and  as input for the final treatment model in Figure <ref> to improve the effect estimates. We denote the combinations  and , which correspond to route 0.854 and route 0.855 in Figure <ref>, respectively. We discuss how OP can be used as predictors to estimate the propensity scores for different treatment modeling techniques in Section <ref>, where we outline several possible choices for the final treatment model.

In summary, both  and  use the outcome model for selecting covariates to account for confounding. Variable selection and propensity score estimation are conducted separately in two steps for methods based on , , , and .

§ POSSIBLE CHOICES FOR THE FINAL TREATMENT MODEL FOR PROPENSITY SCORE ESTIMATION

We consider five different methods (Table <ref>) for estimating the propensity score given a set of candidate predictors. We leave out the OP for now and only consider the alternatives for the final treatment model for routes 0.851-0.853 in Figure <ref>. The first is a multinomial logistic regression model (LOGIS), specified as

\begin{align*}
P(Z_i=j|\tilde{\boldsymbol{X}}_i) = \frac{\exp(\tilde{\boldsymbol{X}}_i^T\boldsymbol{\alpha}_j)}{\sum_{l=1}^J\exp(\tilde{\boldsymbol{X}}_i^T\boldsymbol{\alpha}_l)},
\end{align*}

where $\tilde{\boldsymbol{X}}_i$ is the set of (possibly selected) variables entering the final treatment model. When the dimension of $\tilde{\boldsymbol{X}}_i$ is sufficiently small, the coefficients can be estimated using the maximum likelihood estimator. In the case where dimension reduction is needed, the estimates $\hat{\boldsymbol{\alpha}}$ can be obtained by minimizing the negative log-likelihood with the group LASSO penalty

\begin{align*}
N^{-1}\sum_{i=1}^N\log P(Z_i=j|\tilde{\boldsymbol{X}}_i)+\lambda \sum_{k=1}^p \|\boldsymbol{\alpha}_{\cdot k} \|_2,
\end{align*}

which indicates that each covariate in $\tilde{\boldsymbol{X}}_i$ is associated with all or none of the $J$ levels.

The second is the classification and regression tree (CART) method which performs recursive binary splitting on the feature space in a top-down fashion. At each split, CART agnostically searches for a variable $X$ and a cutpoint such that the response values in each of the resulting nodes lead to the greatest homogeneity [Breiman et al., 1984]. In that sense, CART intrinsically conducts variable selection while growing the tree, as variables that are not predictive of the treatment are less likely to be chosen at each split. We used the Gini index, a measure of the total variance across the $J$ classes, as the metrics for node splitting. Small value of Gini index indicates that observations in this node are predominated by a single class. CART tends to overfit the data. To address the overfitting issue, the common strategy is to first grow a large tree and prune it back in order to retain only part of the tree, as simpler trees tend to be less sensitive to the noises in the data. This method is referred to as pruned CART.

The single tree implementation of both CART and pruned CART, sometimes known as weak learners, may give poor predictions on their own. The ensemble methods, which combine multiple weak learners into one predictive model, have been developed to enhance the predictions. One example is the bootstrap aggregation of the CART algorithm (bagged CART). The bootstrap step of bagged CART involves randomly drawing $n$ observations (i.e., the same size as the original sample) with replacement from the original sample and fitting a CART separately for each bootstrap replicate. The bagging estimates of the probability of subjects being assigned to each class are obtained by averaging the predicted class probabilities from each of the single trees. Another popular ensemble method is the random forests. Similar to bagged CART, random forests build trees based on bootstrap samples of the original observations. What is different from bagged CART is that random forests only considers a random sample of $m$ predictors ($m<p$) at each split, and typically $m\approx\sqrt{p}$ is chosen for classification problems in practice [James et al., 2021]. In our simulation studies and data example, we choose to grow a relatively large number of trees in order to stabilize the out-of-bag error rate.

§.§ Implementation of Incorporating the Predicted Outcome Probabilities

The methods listed in Table <ref> can be applied to , , and  for estimating the propensity scores, which are then used to compute the ATE. The OP is employed for propensity score estimation in different ways for different final treatment models. For the LOGIS method, the propensity scores are estimated by regressing the treatment variable on the union of the input covariates (, , or ) and $\hat{p}_i^\ast$. When the regression model is regularized, no penalty is imposed on $\hat{p}_i^\ast$. In this way the outcome information is guaranteed to be utilized in the propensity score model.

For the tree-based machine learning methods such as CART, there is no straightforward way to force the OP into the tree growing process, where variable selection is intrinsically conducted. We instead conduct a two-step estimation procedure. We first obtain the set of propensity scores $\hat{\boldsymbol{\pi}}(\tilde{\boldsymbol{X}}_i)=\{\hat{\pi}_1(\tilde{\boldsymbol{X}}_i),\cdots,\hat{\pi}_{J-1}(\tilde{\boldsymbol{X}}_i)\}^T$ using the tree-based methods.  Next, we fit a logistic regression model for the treatment as a function of $\hat{p}_i^\ast$ and the propensity scores obtained in route 0.851, 0.852, or 0.853:

\begin{align}
\log\frac{P(Z_i=j|\tilde{\boldsymbol{X}}_i)}{P(Z_i=J|\tilde{\boldsymbol{X}}_i)}=\phi_{0j} + \hat{\boldsymbol{\pi}}(\tilde{\boldsymbol{X}}_i)^T\boldsymbol{\eta}_j+\phi_j\hat{p}_i^\ast.
\label{eq:Chap3_OA_trt_model}
\end{align}

The coefficients $(\boldsymbol{\phi}_0,\boldsymbol{\phi}, \boldsymbol{\eta}_1^T,\cdots,\boldsymbol{\eta}_{J-1})^T$ can be estimated using maximum likelihood. The final propensity scores that take OP into account are then obtained by calculating the predicted probabilities from model (<ref>).

For  (route 0.855) and  (route 0.856), the associations between the covariates and the outcome are used twice in the entire estimation process, one for variable selection using the outcome model, and the other for propensity score estimation using the OP.

1cc]@c@Methods for Constructing
the Treatment Model     c]@c@Simultaneous Variable
Selection and Estimation     c]@c@Ways to Incorporate OP
into the Estimation Process     R package
c]@l@Logistic regression
(LOGIS, standard or penalized)     c]@c@Yes for penalized
logistic regression     c]@c@Used as a regressor for logistic regression.
No penalty is imposed on OP
for penalized logistic regression.     glmnet*, gcdnet
CART     Yes     c]@c@A multinomial logistic regression for the treatment is fitted
as a function of estimated propensity scores and OP     rpart*
Pruned CART     Yes     Same as above     rpart*
Bagged CART     Yes     Same as above     ipred*
Random forests     Yes     Same as above     randomForest*, ranger

Possible choices for the final treatment model.
R packages used to implement the methods in this paper.

§ SIMULATION STUDIES

§.§ Implementation of Methods under Comparison

The comparative methods were implemented in R with default parameters unless otherwise specified. The LASSO algorithm was implemented using the R package glmnet. The tuning parameter $\lambda$ was determined using 10-fold cross-validation with the  criterion for selecting variables in the pre-selection step, as the goal was to reduce the number of covariates entering the final treatment model. For LOGIS based on , which does not involve the pre-selection step, the  criterion was used. CART was implemented using the rpart package with the complexity parameter (cp) being 0.001, which encourages a large and complex tree structure. For pruned CART, the “cp” that corresponded to the smallest 10-fold cross-validated error was used to determine the best pruned tree. Bagged CART was implemented using the ipred package with 200 bootstrap replicates. Random forests were implemented using the randomForest package with 1000 bootstrap replicates. The minimum size of terminal nodes (the  parameter) was set to be 7 in order to make it consistent with the parameters used for CART. For bagged CART and random forests, propensity scores were estimated based on the out-of-bag predictions.

We also extended the OAL approach to three-treatment comparison. Following Shortreed and Ertefaie, 2017, we considered a set of possible values for the tuning parameter $\lambda_n$, $\{n^{-20}, n^{-15}, n^{-10}, n^{-5}, n^{-3}, n^{-1},\\ n^{-0.75}, n^{-0.5}, n^{-0.25}, n^{0.25}, n^{0.49}\}$, and $\lambda_n$ was selected by minimizing a weighted absolute mean difference between treatment groups, a quantity that combines the weighted difference in covariates and the absolute values of the coefficients corresponding to the covariates in the outcome model. Since the adaptive weights for the covariates excluded by the outcome model got inflated to infinity, these covariates were not used to fit the adaptive LASSO for variable selection and propensity score estimation. In that sense, OAL in the high-dimensional setting is equivalent to a LASSO based on  with additional weights imposed on the covariates for variable selection.

§.§ Simulation Setup

For each simulated dataset, $J=3$ treatment groups were compared and $p=100$ covariates were considered, with $|\mathcal{C}|$ confounders, $|\mathcal{Z}|$ related to treatment only, $|\mathcal{Y}|$ related to outcome only, and $|\mathcal{S}|$ spurious predictors. Covariates $\boldsymbol{X}_i$ were generated as follows unless otherwise specified. Half of the covariates (rounded down) in $\boldsymbol{X}_\mathcal{C}$, $\boldsymbol{X}_\mathcal{Z}$, $\boldsymbol{X}_\mathcal{Y}$, and $\boldsymbol{X}_\mathcal{S}$ were generated from a binomial distribution with a probability of 0.3, and the other half were generated from a multivariate normal distribution with a $p'$-dimensional vector of means $\boldsymbol{0}_{p'}$ and a covariance matrix $\Sigma$, where $p'$ is the number of continuous covariates in each subset ($\boldsymbol{X}_\mathcal{C}$, $\boldsymbol{X}_\mathcal{Z}$, $\boldsymbol{X}_\mathcal{Y}$, or $\boldsymbol{X}_\mathcal{S}$) and $\Sigma$ is an identity matrix.

We considered three different simulation settings. Our first setting assumes additivity and linearity for the treatment generating model,

\begin{align}
Z_i\sim\text{Multinomial}\{\pi_1(\boldsymbol{V}_i), \pi_2(\boldsymbol{V}_i), \pi_3(\boldsymbol{V}_i)\},
\label{eq:Chap3:trt_generating_model}
\end{align}

$\pi_j(\boldsymbol{V}_i)=\exp(\boldsymbol{V}_i^T\boldsymbol{\alpha}_j)/\sum_{l=1}^3\exp(\boldsymbol{V}_i^T\boldsymbol{\alpha}_l)$ is the probability of receiving treatment $j$, with $\boldsymbol{V}_i=(1, \boldsymbol{X}_{\mathcal{C}}^T, \boldsymbol{X}_{\mathcal{Z}}^T)^T$. We assumed heterogeneous treatment effects on the outcome, and sampled the potential outcome $Y_i^{(j)}$ from a binomial distribution with probability

\begin{align*}
\end{align*}

where $\boldsymbol{W}_i=(\boldsymbol{X}_\mathcal{C}^T, \boldsymbol{X}_\mathcal{Y}^T)^T$, $(\beta_{01}, \beta_{02}, \beta_{03})=(0, 0.6, 0.4)$, and $\boldsymbol{\beta}_j\propto(1,\cdots,1)^T$. We considered three levels of sparsity for the treatment and outcome models, with $|\mathcal{C}|=|\mathcal{Z}|=|\mathcal{Y}|=5$, $10$, and $20$ for the scenarios with sparse, moderately sparse, and dense models, respectively. As a result, the dimension of $\boldsymbol{\alpha}_j$ differed across the scenarios. The parameter $\boldsymbol{\alpha}=(\boldsymbol{\alpha}_1^T,\boldsymbol{\alpha}_2^T, \boldsymbol{\alpha}_3^T)^T$ was scaled such that $\|\boldsymbol{\alpha}\|_2=5$. Similarly, the signal strength of the outcome model was scaled such that $\|\boldsymbol{\beta}_1\|=3$, $\|\boldsymbol{\beta}_2\|=2$, and $\|\boldsymbol{\beta}_3\|=4$. The true values for $E\{Y^(1)\}$, $E\{Y^(2)\}$, and $E\{Y^(3)\}$ were 0.61, 0.71, and 0.68 for the `sparse' scenario, 0.69, 0.77, and 0.76 for the `moderately sparse' scenario, and 0.76, 0.83, and 0.83 for the `dense' scenario. To examine the performance as the sample size increases, we let the sample size $n$ be 500, 1000, and 2000 for the `sparse' scenario.

Our second setting assumed that models were `sparse' ($|\mathcal{C}|=|\mathcal{Z}|=|\mathcal{Y}|=5$) and considered a set of association equations for the treatment assignment that varied in degrees of nonlinearity and nonadditivity. In this case, $\boldsymbol{V}_i$ in (<ref>) may contain some transformation of covariates in $(\boldsymbol{X}_\mathcal{C}, \boldsymbol{X}_\mathcal{Y})$ and/or interaction effects. The structure of $\boldsymbol{V}_i$ are shown in Table A.1. Specifically, we considered scenarios with nonlinear main effects and no interactions (NL), linear main and interaction effects (L-L), nonlinear main effects and linear interactions (NL-L), and nonlinear main and interaction effects (NL-NL). All confounders were continuous in this setting, and the outcomes were generated using the same model as was used in the first setting except that $(\beta_{01}, \beta_{02}, \beta_{03})=(0, -0.6, 0.4)$.

Our third setting also assumed that the models were sparse and let censoring come into play. Instead of sampling binary outcomes directly from a binomial distribution, we first generated time to event $T_i$ from a logistic distribution with mean function

\begin{align*}
\beta_{01}I(Z_i=1)+\beta_{02}I(Z_i=2)+\beta_{03}I(Z_i=3)+\boldsymbol{W}_i^T\boldsymbol{\beta}
\end{align*}

and scale parameter $s=6$, where $\boldsymbol{W}_i=(\boldsymbol{X}_\mathcal{C}^T, \boldsymbol{X}_\mathcal{Y}^T)^T$, $(\beta_{01}, \beta_{02}, \beta_{03})=(120,100,115)$ and $\boldsymbol{\beta}=(5,\cdots,5)^T$. Then we obtained the outcome such that $Y_i=I\{T_i<130\}$. We generated the censoring time $C$ using inverse transform sampling [Bender et al., 2005] as a function of $\boldsymbol{U}_i=(1,\boldsymbol{X}_\mathcal{C})^T$. Specifically, we assumed a Cox proportional hazards model with the baseline hazard following a Weibull distribution,

\begin{align*}
C_i^{(j)} = \{\lambda^{-1}\exp(\boldsymbol{U}_i^T\boldsymbol{\gamma}_j)^{-1}\log u\}^{1/\nu},
\end{align*}

with scale parameter $\lambda=0.01$ and shape parameter $\nu=7$, where $u$ was sampled from a $\text{Uniform}(0,1)$ distribution. Note that $\boldsymbol{U}_i$ was assumed to be known in our simulation settings, which resembles our data example where censoring is believed to only depend a low-dimensional set of covariates that can be identified by human experts. The proportion of subject not being observed for $d=130$ was around 22%.

We considered the six routes displayed in Figure <ref> for propensity score estimation. We also present the results for three sets of covariates that are usually unknown in practice: confounders only ($\boldsymbol{X}_\mathcal{C}$), treatment predictors ($\boldsymbol{X}_\mathcal{C}\cup\boldsymbol{X}_\mathcal{Z}$), and outcome predictors ($\boldsymbol{X}_\mathcal{C}\cup\boldsymbol{X}_\mathcal{Y}$). These results were used to illustrate the impact of different groups of covariates on the treatment effect estimates. The ATE were estimated using the IPW method.

§.§ Construction of Confidence Intervals

The confidence intervals (CI) were constructed using bootstrap standard errors based on 200 bootstrap replicates. For the settings without censoring, we considered two possible bootstrap procedures. The first applies variable selection to each bootstrap replicate and refits the models for the treatment and outcome using variables selected in each replicate, and we refer to it as usual bootstrap. The second ignores the variability due to selection for covariates. Instead, for LOGIS based on  and  and all tree-based methods, OP and propensity scores are directly bootstrapped from the OP and propensity scores obtained in the original sample, without refitting the model. For LOGIS based on , , , and , propensity scores are obtained by refitting the final treatment model using variables selected in the original data set. As a result, the usual bootstrap is much more computationally intensive than the modified bootstrap. A trial simulation study under the scenario of linear sparse models for sample sizes of 500 (Table <ref>) and 1000 (Table A.2) showed that usual bootstrap tended to overestimate the standard errors and produce overly conservative CIs for the tree-based methods. For example, most of the coverage rates for the methods that involved OP were above 98%. Over-coverage was also observed for (unpenalized) LOGIS for usual bootstrap. On the other hand, LOGIS based on  had close to nominal coverage using usual bootstrap at the expense of high computational burden, and slight over-coverage using modified bootstrap. For the tree-based methods, standard errors based on modified bootstrap were close to their corresponding Monte Carlo standard deviation. In this case, modified bootstrap remedied the overestimation of the usual bootstrap by dropping the variability due to variable selection and estimated the true variability well across all methods. Therefore, we proceed with the modified bootstrap technique for our simulation studies, which directly samples the estimated propensity scores and OP from the original simulated data set for each bootstrap replicate. For the third setting where censoring existed, we again used modified bootstrap, but with the censoring model refitted for each bootstrap sample. Metrics that were used to compare the various propensity score estimation methods included bias, Monte Carlo standard deviations, standard errors, root mean squared error (RMSE), and coverage rate of 95% CIs. True values were determined using $5\times10^5$ replicates.

1l     3cEmpirical SD     3cSE (usual)     3cCoverage (%; usual)     3cSE (modified)     3cCoverage (%; modified)
(lr)2-4 (lr)5-7 (lr)8-10 (lr)11-13 (lr)14-16
1lEstimators     1l1 vs 2     1 vs 3     2 vs 3     1 vs 2     1 vs 3     2 vs 3     1 vs 2     1 vs 3     2 vs 3     1 vs 2     1 vs 3     2 vs 3     1 vs 2     1 vs 3     2 vs 3
1lOAL     54     48     51     57     51     53     96.1     96.0     95.7     55     50     52     94.1     95.6     94.2

All     52     51     50     56     54     53     59.4     73.2     93.4     55     55     52     59.8     75.9     93.3

Ysel     57     49     54     94     82     85     100     99.7     99.9     58     51     55     95.4     96.0     94.8

YZsel     58     52     54     80     70     74     99.7     99.1     99.4     58     53     55     95.0     95.5     95.0

OP + All     51     46     50     57     52     54     95.8     96.9     95.5     58     56     56     96.3     97.8     96.9

OP + Ysel     57     49     54     94     82     85     100     99.7     99.9     58     51     55     95.4     96.0     94.8

OP + YZsel     55     48     53     80     69     73     99.8     99.4     99.3     55     49     53     95.2     96.0     94.3

All     79     81     76     95     96     89     78.1     87.4     97.0     78     80     74     63.2     75.9     92.0

Ysel     72     70     67     87     86     82     90.9     94.4     97.5     72     72     69     81.6     88.1     94.0

YZsel     69     66     64     83     83     79     93.1     96.4     99.0     70     68     67     86.2     90.8     94.8

OP + All     80     74     79     101     93     102     98.8     98.6     98.9     79     74     79     93.4     93.8     94.2

OP + Ysel     67     61     65     87     81     88     99.4     99.0     98.7     66     63     66     94.2     94.4     94.8

OP + YZsel     62     58     62     82     77     84     99.3     99.5     99.4     63     59     62     94.3     94.7     94.8

1lPruned CART

All     62     63     56     92     92     86     75.3     90.4     99.1     61     61     56     40.1     60.4     90.4

Ysel     63     60     55     85     84     79     87.0     93.2     99.0     60     60     56     61.0     75.5     92.4

YZsel     64     61     55     82     81     76     87.2     93.5     99.3     61     60     57     68.6     80.7     93.0

OP + All     58     53     56     95     87     96     99.8     99.5     99.9     58     52     57     94.3     94.5     94.6

OP + Ysel     56     49     54     83     77     84     100     99.5     99.7     55     50     54     94.2     94.8     94.6

OP + YZsel     54     50     53     78     73     80     99.8     99.6     99.6     55     50     54     95.7     94.9     94.7

1lBagged CART

All     54     53     51     45     45     42     28.9     46.9     85.3     56     57     53     44.1     65.5     91.2

Ysel     71     63     63     54     51     49     87.3     88.7     86.2     71     66     67     93.4     95.6     95.4

YZsel     91     82     81     63     59     57     80.0     82.1     80.9     86     79     81     91.6     92.5     94.8

OP + All     53     46     52     86     84     90     100     99.9     99.9     52     47     51     94.2     95.4     94.4

OP + Ysel     50     45     49     80     80     85     100     99.9     100     50     45     50     94.6     95.2     95.0

OP + YZsel     50     45     49     77     77     82     100     99.8     99.9     50     46     50     94.6     95.0     94.6

1lRandom Forests

All     49     50     48     40     40     38     13.4     30.9     80.4     52     53     50     27.1     51.2     90.0

Ysel     57     52     51     41     41     38     68.4     76.2     84.7     59     57     57     84.9     91.2     95.9

YZsel     94     78     79     46     44     41     66.8     75.7     72.0     79     70     75     91.0     91.2     94.5

OP + All     55     47     53     61     56     59     97.5     98.4     97.2     53     48     52     94.8     94.9     94.0

OP + Ysel     51     45     50     70     66     71     99.9     99.5     99.4     51     46     50     94.5     94.8     94.4

OP + YZsel     50     46     50     79     76     82     99.9     99.7     99.9     51     46     50     94.8     95.0     94.7

[Standard errors (SE) and coverage of 95% confidence intervals estimated by usual bootstrap and modified bootstrap for sample size of 500.]Standard errors (SE) and coverage of 95% confidence intervals estimated by usual bootstrap and modified bootstrap for sample size of 500. The scenario with sparse treatment models was considered. Results were obtained based on 1000 simulated datasets. For each dataset, 200 bootstrap samples were generated.

1l                         3cc]@c@Bias
$\times$1000     3cc]@c@Empirical SD
$\times$1000     3cc]@c@SE (usual)
$\times$1000     3cc]@c@Coverage (%)
(usual)         3cc]@c@SE (modified)
$\times$1000     3cc]@c@Coverage (%)
(lr)2-4 (lr)5-7 (lr)8-10 (lr)11-13 (lr)14-16 (lr)17-19
1lEstimators                1 vs 2                       1 vs 3                      2 vs 3                      1 vs 2                         1 vs 3                         2 vs 3                         1 vs 2                         1 vs 3                        2 vs 3                        1 vs 2                          1 vs 3                         2 vs 3                         1 vs 2                          1 vs 3                         2 vs 3                         1 vs 2                          1 vs 3                         2 vs 3

1lOAL             101                          77                          -25                         54                             48                             51                             57                             51                            53                            96.1                            96.0                           95.7                           55                              50                             52                             94.1                            95.6                           94.2

All                                          96                           75                          -21                         52                             48                             49                             56                             54                            53                            59.4                            73.2                           93.4                           55                              55                             52                             59.8                            75.9                           93.3

Ysel                                         -6                           -2                          5                           56                             49                             55                             94                             82                            85                            100.0                           99.7                           99.9                           58                              51                             55                             95.4                            96.0                           94.8

YZsel                                        -1                           3                           4                           57                             50                             55                             80                             70                            74                            99.7                            99.1                           99.4                           58                              53                             55                             95.0                            95.5                           95.0

OP+All                                       15                           15                          0                           50                             45                             49                             57                             52                            54                            95.8                            96.9                           95.5                           58                              56                             56                             96.3                            97.8                           96.9

OP+Ysel                                      -6                           -2                          5                           56                             49                             55                             94                             82                            85                            100.0                           99.7                           99.9                           58                              51                             55                             95.4                            96.0                           94.8

OP+YZsel                                     -5                           -1                          5                           54                             47                             54                             80                             69                            73                            99.8                            99.4                           99.3                           55                              49                             53                             95.2                            96.0                           94.3

All                                          125                          95                          -30                         79                             78                             74                             95                             96                            89                            78.1                            87.4                           97.0                           78                              80                             74                             63.2                            75.9                           92.0

Ysel                                         74                           57                          -17                         68                             68                             66                             87                             86                            82                            90.9                            94.4                           97.5                           72                              72                             69                             81.6                            88.1                           94.0

YZsel                                        59                           45                          -14                         66                             62                             62                             83                             83                            79                            93.1                            96.4                           99.0                           70                              68                             67                             86.2                            90.8                           94.8

OP+All                                       10                           10                          0                           77                             73                             78                             101                            93                            102                           98.8                            98.6                           98.9                           79                              74                             79                             93.4                            93.8                           94.2

OP+Ysel                                      14                           12                          -3                          63                             61                             64                             87                             81                            88                            99.4                            99.0                           98.7                           66                              63                             66                             94.2                            94.4                           94.8

OP+YZsel                                     14                           10                          -4                          60                             56                             61                             82                             77                            84                            99.3                            99.5                           99.4                           63                              59                             62                             94.3                            94.7                           94.8

1lPruned CART

All                                          134                          102                         -32                         63                             61                             55                             92                             92                            86                            75.3                            90.4                           99.1                           61                              61                             56                             40.1                            60.4                           90.4

Ysel                                         97                           74                          -23                         61                             59                             54                             85                             84                            79                            87.0                            93.2                           99.0                           60                              60                             56                             61.0                            75.5                           92.4

YZsel                                        81                           63                          -19                         64                             60                             54                             82                             81                            76                            87.2                            93.5                           99.3                           61                              60                             57                             68.6                            80.7                           93.0

OP+All                                       -3                           1                           4                           56                             51                             55                             95                             87                            96                            99.8                            99.5                           99.9                           58                              52                             57                             94.3                            94.5                           94.6

OP+Ysel                                      1                            4                           2                           52                             49                             54                             83                             77                            84                            100.0                           99.5                           99.7                           55                              50                             54                             94.2                            94.8                           94.6

OP+YZsel                                     2                            3                           1                           53                             48                             53                             78                             73                            80                            99.8                            99.6                           99.6                           55                              50                             54                             95.7                            94.9                           94.7

1lBagged CART

All                                          120                          91                          -29                         54                             51                             49                             45                             45                            42                            28.9                            46.9                           85.3                           56                              57                             53                             44.1                            65.5                           91.2

Ysel                                         17                           17                          0                           67                             61                             64                             54                             51                            49                            87.3                            88.7                           86.2                           71                              66                             67                             93.4                            95.6                           95.4

YZsel                                        -35                          -26                         8                           90                             80                             85                             63                             59                            57                            80.0                            82.1                           80.9                           86                              79                             81                             91.6                            92.5                           94.8

OP+All                                       -5                           1                           6                           51                             46                             50                             86                             84                            90                            100.0                           99.9                           99.9                           52                              47                             51                             94.2                            95.4                           94.4

OP+Ysel                                      -7                           0                           7                           47                             44                             48                             80                             80                            85                            100.0                           99.9                           100.0                          50                              45                             50                             94.6                            95.2                           95.0

OP+YZsel                                     -7                           -1                          6                           47                             44                             49                             77                             77                            82                            100.0                           99.8                           99.9                           50                              46                             50                             94.6                            95.0                           94.6

2lRandom   Forests

All                                          135                          103                         -32                         50                             48                             47                             40                             40                            38                            13.4                            30.9                           80.4                           52                              53                             50                             27.1                            51.2                           90.0

Ysel                                         53                           43                          -10                         56                             50                             51                             41                             41                            38                            68.4                            76.2                           84.7                           59                              57                             57                             84.9                            91.2                           95.9

YZsel                                        -36                          -22                         14                          91                             72                             78                             46                             44                            41                            66.8                            75.7                           72.0                           79                              70                             75                             91.0                            91.2                           94.5

OP+All                                       -5                           1                           6                           53                             46                             51                             61                             56                            59                            97.5                            98.4                           97.2                           53                              48                             52                             94.8                            94.9                           94.0

OP+Ysel                                      -9                           -1                          8                           48                             44                             49                             70                             66                            71                            99.9                            99.5                           99.4                           51                              46                             50                             94.5                            94.8                           94.4

OP+YZsel                                     -8                           -1                          7                           48                             45                             49                             79                             76                            82                            99.9                            99.7                           99.9                           51                              46                             50                             94.8                            95.0                           94.7

[Standard errors (SE) and coverage of 95% confidence intervals estimated by usual bootstrap and modified bootstrap for sample size of 500.]Standard errors (SE) and coverage of 95% confidence intervals estimated by usual bootstrap and modified bootstrap for sample size of 500. The scenario with sparse treatment models was considered. Results were obtained based on 1000 simulated datasets. For each dataset, 200 bootstrap samples were generated.

§.§ Simulation Results

We present the box plots of bias for the IPW estimates across the simulation settings in Figures <ref>-<ref>. The numerical results for all evaluation metrics are reported in the Supplementary Materials.

§.§.§ Bias

The numerical results for the setting of linear treatment models are presented in Tables A.3-A.9 in the Supplementary Materials. LOGIS that adjusted for confounders ($\boldsymbol{X}_\mathcal{C}$), treatment predictors ($\boldsymbol{X}_\mathcal{C}\cup\boldsymbol{X}_\mathcal{Z}$), and outcome predictors ($\boldsymbol{X}_\mathcal{C}\cup\boldsymbol{X}_\mathcal{Y}$) had close to zero empirical bias as expected (e.g., Table A.3). For the CART family methods, using confounders only generally yielded smaller bias than using treatment predictors, outcome predictors, or . Methods based on pre-selected predictors ( and ) tended to yield smaller empirical bias than methods using all the available predictors as input (Figure <ref>), which indicates that excluding noise variables before fitting the final treatment model can help remove the bias, though the absolute bias was still greater than zero for the tree-based methods. One exception was that random forests based on  resulted in substantial bias in the `sparse' scenario, possibly because there were too few true predictors at each split for random forests to choose from. The OAL method had similar performance in terms of bias to LOGIS based on  and  in the `sparse' scenario (Tables A.3-A.5), while the latter outperformed the former in the `moderately sparse' (Table A.6) and `dense' scenarios (Tables A.8 and A.9).

When the outcome information was not taken into account, LOGIS in general produced less biased estimates than the tree-based methods for each of the routes 0.851-0.853 in Figure <ref>. Bias for , , and  were getting closer as the model became `denser' (Figure <ref>). As sample size increased from 500 to 2000 for the `sparse' scenario (Figure <ref>), the nonzero empirical bias persisted across the methods considered, which illustrates the slow convergence rate of LASSO. The inclusion of OP greatly reduced the bias of the estimates in the case where methods based on , , or  resulted in larger than zero absolute bias across the simulation settings, which highlights the robustness of the treatment effect estimator provided by incorporating OP into the estimation process.

With nonlinearity and nonadditivity in the treatment model, the performance of LOGIS was not inferior to the tree-based methods in terms of bias (Tables A.10-A.17), possibly because in this case, LOGIS approximated nonlinear functions reasonably well. The good approximation of linear methods to nonlinear functions was also observed in Tu, 2019, where multivariable linear regression yielded smaller bias than bagged CART and random forests in some cases with nonlinear and nonadditive associations in the treatment models.

When censoring existed, LOGIS based on , , , and  produced much lower bias than their CART family counterparts. One exception was -based bagged CART, which had close to zero bias (Table A.18).

§.§.§ Statistical Efficiency and RMSE

In general, estimates produced by routes 0.854, 0.855, and 0.856 had smaller variability than those resulting from routes 0.851, 0.852 and 0.853, respectively, for each method across all the settings, which illustrates the advantage of leveraging the information about the outcome model in terms of statistical efficiency.

Methods based on , , and  had smaller variabilities and RMSE of the effect estimates across the simulation settings compared to methods based on , , and , respectively (Figures B.1-B.3). For example, for the scenario with moderately sparse models, the percentage of reduction in RMSE ranged from 1.6% to 20.6% for LOGIS based on  compared to , and from 26.3% to 41.6% for random forests based on  compared to  (Table A.6). The variabilities for , , and  were close to one another, and one was not uniformly smaller than the other two. In general, regardless of whether OP was incorporated, CART had the largest variability and RMSE among the methods for all scenarios considered.

In the presence of censoring, OAL and LOGIS based on  had the smallest RMSE among the methods under comparison (Table A.18). The inclusion of OP reduced the variability of the effect estimates for all method considered. The bias for LOGIS, CART, and pruned CART decreased when the OP was taken into account, with  resulting in the smallest bias for each method. However, we still observed residual bias, especially for the CART and pruned-CART.

§.§.§ Coverage of 95% CI

For the tree-based methods and unpenalized LOGIS, standard errors obtained using the modified bootstrap were close to the corresponding Monte Carlo standard deviations in the scenarios with `sparse' models (regardless of whether the treatment models were linear and/or additive). The modified bootstrap tended to slightly overestimate the variability for LOGIS with LASSO penalty. For example, the ratio of standard errors to Monte Carlo standard deviations ranged from 1.12 to 1.21 for  in the scenario with linear associations and `sparse' representation of the models for sample size of 500. The tree-based methods achieved close to nominal coverage of 95% for , , and  for most of the scenarios considered. In the case where the coverage fell below the nominal level, for example the `dense' scenario (Table A.8), the under-coverage was mainly caused by the empirical bias rather than the underestimation of the standard errors.

§ DATA ANALYSIS

We applied the algorithms considered in Figure <ref> and the propensity score estimation methods listed in Table <ref> to a dataset of patients with metastatic castration-resistant prostate cancer (mCRPC) from the Optum Clinformatics Data Mart [Ross et al., 2021]. Patients who used at least one of the six drugs (docetaxel, abiraterone, enzalutamide, sipuleucel-T, cabazitaxel, and radium-223) approved to treat mCRPC from January 1, 2014, to December 31, 2019 were included in the analytic cohort. We excluded the patients who received cabazitaxel ($n=56$) or radium-223 ($n=28$) as their first-line therapy from our analysis, since there were much fewer samples in these two groups than the other four. The adverse effects of the four drugs for mCRPC were compared, with the outcome being the occurrence of at least one emergency room (ER) visit and hospitalization, respectively, within 180 or 360 days of treatment initiation. In the previous analysis, the working model for the outcome was specified as

\begin{align*}
\log\frac{P(Y_i=1|\boldsymbol{A}_i,\boldsymbol{B}_i)}{P(Y_i=0|\boldsymbol{A}_i,\boldsymbol{B}_i)} = \boldsymbol{A}_i^T\boldsymbol{\beta}_A+\boldsymbol{B}_i^T\boldsymbol{\beta}_B,
\end{align*}

where $\boldsymbol{A}_i$ contained the sociodemographic factors and other relevant covariates, and $\boldsymbol{B}_i$ contained five pre-existing comorbid conditions, including diabetes, hypertension, arrhythmia, congestive heart failure, and osteoporosis. Specifically, $\boldsymbol{A}_i$ included age, race, education level, household income, geographic region, insurance product type, whether the insurance plan is administrative services only (ASO), metastatic status of cancer, year of first prescription, and provider type [Caram et al., 2019]. In this analysis, we increased the granularity of the comorbid conditions and considered a list of phenotype codes (phecodes), which are aggregations of the International Classification of Diseases (ICD) codes that represent clinically meaningful phenotypes [Denny et al., 2010]. We matched the ICD codes in the claims data to the list of phecodes and identified 1042 phecodes as the original reservoir of predictors. These 1042 phecodes represent 16 broad categories of diseases (circulatory system, congenital anomalies, dermatological diseases, endocrine/metabolic diseases, genitourinary diseases, hematopoietic diseases, infectious diseases, injuries and poisonings, mental disorders, musculoskeletal diseases, neoplasms, neurological diseases, respiratory diseases, sense organs, and symptoms). The working model for the outcome became

\begin{align*}
\log\frac{P(Y_i=1|\boldsymbol{A}_i,\boldsymbol{M}_i)}{P(Y_i=0|\boldsymbol{A}_i,\boldsymbol{M}_i)} = \boldsymbol{A}_i^T\boldsymbol{\beta}_A+\boldsymbol{M}_i^T\boldsymbol{\beta}_M,
\end{align*}

where $\boldsymbol{M}_i$ denotes the list of phecodes of dimension 1042. In the pre-selection step, a LASSO was fitted to select the phecodes that are predictive of the outcome, with the coefficients $\boldsymbol{\beta}_A$ unpenalized. Covariates predictive of the treatment were selected in a similar manner using a multinomial logistic regression with LASSO-type penalty. In other words, covariates in $\boldsymbol{A}_i$ (such as age, race, and household income) are always adjusted for in the treatment and outcome models. The six routes in Figure <ref> were followed to obtain six different sets of estimated propensity scores. Note that standard multinomial logistic regression models adjusting for  and  yielded substantial standard errors for the estimates of ATE (results not shown). Instead, we fitted the models using the LASSO-type penalty (with $\boldsymbol{\beta}_A$ not being penalized) to reduce the variability of the estimates. For censoring, we fitted a Cox model adjusting for $\boldsymbol{A}_i$ and $\boldsymbol{B}_i$, which was a low-dimensional set of covariates. All covariates were coded binary, and the covariates with more than two levels were represented by dummy variables. The standard errors and CIs were constructed using the modified bootstrap procedure.

§.§ Data Analysis Results

The sample sizes for the ER visit cohort and the hospitalization cohort were 7678 and 7709, respectively. The descriptive statistics of the data and the proportions of patients being censored for each treatment group are summarized elsewhere [Yu et al., 2022]. Specifically, the overall proportions of patients being censored within 180 days and 360 days were 20.8% and 32.6%, respectively, for ER visits, and 24.6% and 40.6%, respectively, for hospitalization. Sipuleucel-T group had larger percentage of censored patients than the other three groups.

Propensity scores estimated by different routes in Figure <ref> were highly correlated for each treatment level (results not shown). In general, we observed larger correlations among propensity scores estimated by LOGIS than those estimated by the tree-based methods. For example, correlations between LOGIS based on  and  ranged from 0.97 to 1, while the correlations between random forests based on  and  ranged from 0.93 to 0.99.

The number of phecodes in each disease group selected by the outcome and/or the treatment model in the pre-selection step for each of the two endpoints are reported in Tables A.19-A.22 in the Supplementary Materials. For example, 54 phecodes were selected by the outcome model and 69 were selected by the treatment model, with 12 lying in the intersection for ER visits within 180 days. Among the 12 phecodes predictive of both treatment and outcome, 4 were associated with neoplasm (cancer of prostate, secondary malignancy of respiratory organs, secondary malignant neoplasm, secondary malignant neoplasm of liver), 2 were associated with circulatory system (congestive heart failure NOS and congestive heart failure nonhypertensive), 2 were associated with mental disorders (delirium dementia and amnestic and other cognitive disorders, tobacco use disorder), 1 was associated with endocrine/metabolic system (type 2 diabetes), 1 was associated with hematopoietic system (lymphadenitis), 1 was associated with respiratory system (abnormal findings examination of lungs), and 1 was associated with symptoms (nausea and vomiting). For hospitalization within 180 days, 63 and 79 phecodes were selected by the outcome and treatment model, respectively, and the intersection contained 10 phecodes. Among the 10 phecodes, 7 were overlapped with those identified for ER visits within 180 days (congestive heart failure NOS, congestive heart failure nonhypertensive, lymphadenitis, cancer of prostate, secondary malignant neoplasm, secondary malignant neoplasm of liver, and nausea and vomiting), with 3 additional phecodes associated with neoplasm (cancer of stomach, malignant neoplasm of head, face, and neck, and secondary malignancy of respiratory organs). We also note that a number of phecodes for genitourinary system were identified to be related either to the treatment (e.g., acute cystitis, chronic renal failure (CKD), nephritis, nephrosis, renal sclerosis, other disorders of the kidney and ureters, renal failure, retention of urine, and urinary tract infection) or to the outcome (e.g., chronic kidney disease stage IV, functional disorders of bladder, hyperplasia of prostate, lump or mass in breast, other disorders of prostate, and prostatitis), while the intersection was empty, possibly due to the fine granularity of the phecodes.

Figure <ref> showed results for (penalized) LOGIS, CART, and Bagged CART for the 180-day risks differences in ER visits among the four treatment groups. Docetaxel users exhibited significantly higher 180-day risks of at least one ER visit than the users of the two oral drugs (abiraterone and enzalutamide), a finding that is consistent with previous studies [Yu et al., 2021, Tonyali et al., 2017]. For example, compared to docetaxel users,
the 180-day risk of ER visits for abiraterone users is 11.2% lower (estimated ATE: -0.112 (-0.147,-0.078 )) and for enzalutamide users is 13.5% lower (estimated ATE: -0.135 (-0.177, -0.093)), from the penalized LOGIS method using . Findings from CART and Bagged CART are consistent with the penalized LOGIS method. The 180-day risk differences between abiraterone and enzalutamide users were generally not significant, except that some of the results yielded by random forests indicated significantly lower risk for the enzalutamide group (Figures <ref> and B.4). For the 360-day time window, enzalutamide users showed significantly lower risk of ER visits in most cases (Figures B.5 and B.6).

Similar directions of the comparisons among the four treatment groups were observed for 180-day and 360-day risks in hospitalization. In particular, patients who received enzalutamide as their first-line therapy had significantly lower risk of hospitalization than those who received abiraterone, which is consistent with the findings of a study based on a French insurance system database [Scailteux et al., 2021].

As was observed in the simulation studies, bagged CART and random forests using  as input yielded estimates with large standard errors. In general, incorporating the OP into the treatment model reduced the variability of the estimates and led to narrower 95% CIs. Greater efficiency gains were noted for the CART family methods compared to LOGIS. For example, for the 360-day risk of ER visits, the ratios of CI widths of  over  ranged from 0.98 to 1 and from 0.61 to 0.66 for LOGIS and bagged CART, respectively. When the OP were not included in the treatment model, the estimates yielded by LOGIS tended to have lower variance than those produced by CART family methods. For example, for the 360-day risk of ER visits, the ratios of CI widths of bagged CART over CI widths of LOGIS using  as input ranged from 1.06 to 2.01. Consistent with what was observed in the simulation studies, the point estimates and confidence widths of  and  were very close across the propensity score estimation methods.

§ DISCUSSION

In this paper, we proposed to incorporate the outcome-covaraite relationship by including predicted outcome probability (OP) in the propensity score model. We examined the traditional multinomial logistic regression method and a set of machine learning techniques combined with OP for estimating ATE. Although the idea of using outcome models to improve efficiency has been explored for regularization methods, how to incorporate outcome-covariate relationship into tree-based machine learning methods has been less studied. Simulation studies show that simultaneous variable selection and propensity score estimation (i.e., methods based on ) in a high-dimensional setting led to substantial bias for the LOGIS method, possibly because of the slow convergence rate resulting from the large number of noise variables. Similar pattern was observed for the tree-based methods. We showed that the inclusion of OP can improve the robustness and statistical efficiency of the treatment effect estimators. If the variable selection step had satisfactory performance in terms of identifying the set of important confounders and controlling for the bias, then the benefits of including OP in terms of bias reduction may be minimal. On the other hand, if methods based on , , and  produced biased estimates of ATE, then further adjusting for OP in the treatment model can help reduce the bias. OP alleviates the bias by adding back the information about the confounders that are potentially missed by the variable selection procedure.

The LOGIS method (both standard and penalized depending on the dimension of the covariates) outperformed the tree-based machine learning methods in the majority of scenarios, especially when OP is used. The bias for the LOGIS method could be smaller than the tree-based methods even under conditions of nonlinearity and/or nonadditivity. On the other hand, the nonparametric machine learning methods such as bagged CART and random forests could produce less biased estimates than multinomial logistic regression when only the treatment-covariate relationship was considered. The performance of tree-based methods may be improved by optimizing the tuning parameters, such as minimum size of terminal nodes, maximum number of terminal nodes, and number of trees to grow. In addition, standard cross-validation procedure, which focuses on out-of-sample performance, is often used to optimize tuning parameters for accurate predictive performance in practice and may not have desired characteristics for selecting tuning parameters for causal inference (i.e., unbiased treatment effect estimates) [Häggström and de Luna, 2014]. The selection criteria for the tuning parameters in the context of treatment effect estimation could be a topic of future work.

§ ACKNOWLEDGMENT

This research is supported by NIH (1UG3CA267907), NSF (DMS 1712933) and NIH (R01HG008773).

§ CONFLICT OF INTEREST

Authors state no conflict of interest.

§ DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from OptumInsight but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Data are however available from the authors upon reasonable request and with permission of OptumInsight (https://ihpi.umich.edu/member-resources/data-and-methods/available-datasets).

[Rosenbaum and Rubin, 1983]
P. R. Rosenbaum and D. B. Rubin.
The central role of the propensity score in observational studies for
causal effects.
Biometrika, 700 (1):0 41–55, April 1983.
ISSN 0006-3444.

[Stuart, 2010]
E. A. Stuart.
Matching methods for causal inference: A review and a look
Statistical science : a review journal of the Institute of
Mathematical Statistics, 250 (1):0 1–21, February 2010.
ISSN 0883-4237.

[Yang et al., 2016]
S. Yang, G. W. Imbens, Z. Cui, D. E. Faries, and Z. Kadziola.
Propensity score matching and subclassification in observational
studies with multi-level treatments.
Biometrics, 720 (4):0 1055–1065, 2016.
ISSN 1541-0420.

[Lunceford and Davidian, 2004]
J. K. Lunceford and M. Davidian.
Stratification and weighting via the propensity score in estimation
of causal treatment effects: A comparative study.
Statistics in Medicine, 230 (19):0 2937–2960,
October 2004.
ISSN 1097-0258.

[Zhou et al., 2019]
T. Zhou, M. R. Elliott, and R. J. A. Little.
Penalized Spline of Propensity Methods for Treatment
Journal of the American Statistical Association, 1140
(525):0 1–19, January 2019.
ISSN 0162-1459.

[Yu et al., 2021]
Y. Yu, M. Zhang, X. Shi, M. E. V. Caram, R. J. A. Little, and B. Mukherjee.
A comparison of parametric propensity score-based methods for causal
inference with multiple treatments and a binary outcome.
Statistics in Medicine, January 2021.
ISSN 1097-0258.

[Schneeweiss et al., 2009]
S. Schneeweiss, J. A. Rassen, R. J. Glynn, J. Avorn, H. Mogun, and M. A.
High-dimensional propensity score adjustment in studies of treatment
effects using health care claims data.
Epidemiology (Cambridge, Mass.), 200 (4):0
512–522, July 2009.
ISSN 1531-5487.

[Kang and Schafer, 2007]
J. D. Y. Kang and J. L. Schafer.
Demystifying Double Robustness: A Comparison of Alternative
Strategies for Estimating a Population Mean from Incomplete
Statistical Science, 220 (4):0 523–539,
November 2007.
ISSN 0883-4237, 2168-8745.

[Setoguchi et al., 2008]
S. Setoguchi, S. Schneeweiss, M. A. Brookhart, R. J. Glynn, and E. F. Cook.
Evaluating uses of data mining techniques in propensity score
estimation: A simulation study.
Pharmacoepidemiology and Drug Safety, 170
(6):0 546–555, June 2008.
ISSN 1099-1557.

[Lee et al., 2010]
B. K. Lee, J. Lessler, and E. A. Stuart.
Improving propensity score weighting using machine learning.
Statistics in medicine, 290 (3):0 337–346,
February 2010.
ISSN 0277-6715.

[Ju et al., 2019]
C. Ju, M. Combs, S. D. Lendle, J. M. Franklin, R. Wyss, S. Schneeweiss, and
M. J. van der Laan.
Propensity score prediction for electronic healthcare databases using
super learner and high-dimensional propensity score methods.
Journal of Applied Statistics, 460 (12):0
2216–2236, September 2019.
ISSN 0266-4763.

[Tu, 2019]
C. Tu.
Comparison of various machine learning algorithms for estimating
generalized propensity score.
Journal of Statistical Computation and Simulation, 890
(4):0 708–719, March 2019.
ISSN 0094-9655.

[De Luna et al., 2011]
X. De Luna, I. Waernbaum, and T. S. Richardson.
Covariate selection for the nonparametric estimation of an average
treatment effect.
Biometrika, 980 (4):0 861–875, December 2011.
ISSN 0006-3444.

[Shortreed and Ertefaie, 2017]
S. M. Shortreed and A. Ertefaie.
Outcome-adaptive lasso: Variable selection for causal inference.
Biometrics, 730 (4):0 1111–1122, 2017.
ISSN 1541-0420.

[Brookhart et al., 2006]
M. A. Brookhart, S. Schneeweiss, K. J. Rothman, R. J. Glynn, J. Avorn, and
T. Stürmer.
Variable selection for propensity score models.
American journal of epidemiology, 1630 (12):0
1149–1156, June 2006.
ISSN 0002-9262.

[Zou, 2006]
H. Zou.
The Adaptive Lasso and Its Oracle Properties.
Journal of the American Statistical Association, 1010
(476):0 1418–1429, December 2006.
ISSN 0162-1459.

[Ju et al., 2020]
C. Ju, D. Benkeser, and M. J. van der Laan.
Robust inference on the average treatment effect using the outcome
highly adaptive lasso.
Biometrics, 760 (1):0 109–118, 2020.
ISSN 1541-0420.

[Zigler and Dominici, 2014]
C. M. Zigler and F. Dominici.
Uncertainty in Propensity Score Estimation: Bayesian Methods
for Variable Selection and Model Averaged Causal Effects.
Journal of the American Statistical Association, 1090
(505):0 95–107, January 2014.
ISSN 0162-1459.

[Rubin, 1974]
D. B. Rubin.
Estimating causal effects of treatments in randomized and
nonrandomized studies.
Journal of Educational Psychology, 660 (5):0
688–701, 1974.
ISSN 1939-2176(Electronic),0022-0663(Print).

[Yu et al., 2022]
Y. Yu, M. Zhang, and B. Mukherjee.
An inverse probability weighted regression method that accounts for
right-censoring for causal inference with multiple treatments and a binary
arXiv preprint arXiv:2205.08052, 2022.

[Tibshirani, 1996]
R. Tibshirani.
Regression Shrinkage and Selection via the Lasso.
Journal of the Royal Statistical Society. Series B
(Methodological), 580 (1):0 267–288, 1996.
ISSN 0035-9246.

[Meinshausen, 2007]
N. Meinshausen.
Relaxed Lasso.
Computational Statistics & Data Analysis, 520
(1):0 374–393, September 2007.
ISSN 0167-9473.

[Friedman et al., 2010]
J. Friedman, T. Hastie, and R. Tibshirani.
Regularization Paths for Generalized Linear Models via
Coordinate Descent.
Journal of statistical software, 330 (1):0
1–22, 2010.
ISSN 1548-7660.

[Breiman et al., 1984]
L. Breiman, J. Friedman, C. J. Stone, and R. A. Olshen.
Classification and Regression Trees.
CRC Press, January 1984.
ISBN 978-0-412-04841-8.

[James et al., 2021]
G. James, D. Witten, T. Hastie, and R. Tibshirani.
Tree-Based Methods.
In G. James, D. Witten, T. Hastie, and R. Tibshirani, editors,
An Introduction to Statistical Learning: With Applications
in R, Springer Texts in Statistics, pages 327–365. Springer
US, New York, NY, 2021.
ISBN 978-1-07-161418-1.

[Bender et al., 2005]
R. Bender, T. Augustin, and M. Blettner.
Generating survival times to simulate Cox proportional hazards
Statistics in Medicine, 240 (11):0 1713–1723,
ISSN 1097-0258.

[Ross et al., 2021]
R. D. Ross, X. Shi, M. E. Caram, P. A. Tsao, P. Lin, A. Bohnert, M. Zhang, and
B. Mukherjee.
Veridical causal inference using propensity score methods for
comparative effectiveness research with medical claims.
Health Services and Outcomes Research Methodology, 210
(2):0 206–228, 2021.

[Caram et al., 2019]
M. E. V. Caram, R. Ross, P. Lin, and B. Mukherjee.
Factors Associated With Use of Sipuleucel-T to Treat
Patients With Advanced Prostate Cancer.
JAMA network open, 20 (4):0 e192589, April
ISSN 2574-3805.

[Denny et al., 2010]
J. C. Denny, M. D. Ritchie, M. A. Basford, J. M. Pulley, L. Bastarache,
K. Brown-Gentry, D. Wang, D. R. Masys, D. M. Roden, and D. C. Crawford.
PheWAS: Demonstrating the feasibility of a phenome-wide scan to
discover gene–disease associations.
Bioinformatics, 260 (9):0 1205–1210, May
ISSN 1367-4803.

[Tonyali et al., 2017]
S. Tonyali, H. B. Haberal, and E. Sogutdelen.
Toxicity, Adverse Events, and Quality of Life Associated
with the Treatment of Metastatic Castration-Resistant Prostate
Current Urology, 100 (4):0 169–173, November
ISSN 1661-7649.

[Scailteux et al., 2021]
L.-M. Scailteux, F. Despas, F. Balusson, B. Campillo-Gimenez, R. Mathieu,
S. Vincendeau, A. Happe, E. Nowak, S. Kerbrat, and E. Oger.
Hospitalization for adverse events under abiraterone or enzalutamide
exposure in real-world setting: A French population-based study on
prostate cancer patients.
British Journal of Clinical Pharmacology, July 2021.
ISSN 1365-2125.

[Häggström and de Luna, 2014]
J. Häggström and X. de Luna.
Targeted smoothing parameter selection for estimating average causal
Computational Statistics, 290 (6):0
1727–1748, December 2014.
ISSN 1613-9658.

[Box plots of empirical bias for 2000 inverse probability weighted estimates for the ATE under scenarios with different levels of sparsity.]Box plots of empirical bias for 2000 inverse probability weighted estimates for the ATE under scenarios with different levels of sparsity. The rows represent scenarios and columns represent treatment pairs. Each simulated dataset contained 500 samples.

[Box plots of empirical bias for 2000 inverse probability weighted estimates for the ATE for different sample sizes.]Box plots of empirical bias for 2000 inverse probability weighted estimates for the ATE for different sample sizes. The scenario with sparse treatment and outcome models was considered. The rows represent sample sizes and columns represent treatment pairs.

[Box plots of empirical bias for 2000 inverse probability weighted estimates for the ATE under scenarios with various degrees of nonlinearity and nonadditivity in the treatment generating model.]Box plots of empirical bias for 2000 inverse probability weighted estimates for the ATE under scenarios with various degrees of nonlinearity and nonadditivity in the treatment generating model. The rows represent scenarios and columns represent treatment pairs. Each simulated dataset contained 500 samples.

[Average treatment effects for ER visits within 180 days of treatment initiation for LOGIS, CART, and bagged CART.]Average treatment effects for ER visits within 180 days of treatment initiation for LOGIS, CART, and bagged CART. Data were obtained from Optum Clinformative Data Mart. Total sample size was $N=7678$ ($N_A=2757$, $N_D=2311$, $N_E=2043$, $N_S=567$). Confidence intervals that exclude zero are highlighted in orange. Abbreviations: A, abiraterone; D, docetaxel; E, enzalutamide; S, sipuleucel-T.

[Average treatment effects for hospitalization within 180 days of treatment initiation for LOGIS, CART, and bagged CART.]Average treatment effects for hospitalization within 180 days of treatment initiation for LOGIS, CART, and bagged CART. Data were obtained from Optum Clinformative Data Mart. Total sample size was $N=7709$ ($N_A=2766$, $N_D=2320$, $N_E=2051$, $N_S=572$. Confidence intervals that exclude zero are highlighted in orange. Abbreviations: A, abiraterone; D, docetaxel; E, enzalutamide; S, sipuleucel-T.

§ SUPPLEMENTARY MATERIALS

§.§ Supplementary Tables

Scenario     Predictors for Treatment Generating Model
Linear main effects only (L)$^\ast$     $\begin{aligned}[t]&(X_{\mathcal{C}1}, X_{\mathcal{C}2}, X_{\mathcal{C}3},  X_{\mathcal{C}4}, X_{\mathcal{C}5}, X_{\mathcal{Z}1}, X_{\mathcal{Z}2}, X_{\mathcal{Z}3},  X_{\mathcal{Z}4}, X_{\mathcal{Z}5})\end{aligned}$

Nonlinear main effects only (NL)     c]@c@$\begin{aligned}[t](X_{\mathcal{C}1}(X_{\mathcal{C}1}>0), \exp\{X_{\mathcal{C}2}\}^{-1/2}, |X_{\mathcal{C}3}|,  \log|X_{\mathcal{C}4}+1|, |X_{\mathcal{C}5}|^{-1/2},\\ X_{\mathcal{Z}1}, X_{\mathcal{Z}2}, X_{\mathcal{Z}3},  X_{\mathcal{Z}4}^2, \log|X_{\mathcal{Z}5}|)\end{aligned}$

c]@l@Linear main effects
and linear interactions (L-L)     c]@c@$\begin{aligned}[t](X_{\mathcal{C}1}, X_{\mathcal{C}2}, X_{\mathcal{C}3},  X_{\mathcal{C}4}, X_{\mathcal{C}5},\\ X_{\mathcal{Z}1}, X_{\mathcal{Z}2}, X_{\mathcal{Z}3},  X_{\mathcal{Z}4}, X_{\mathcal{Z}5},\\ X_{\mathcal{Z}1}\times X_{\mathcal{Z}2}, X_{\mathcal{Z}1}\times X_{\mathcal{Z}3}, X_{\mathcal{Z}1}\times X_{\mathcal{Z}4}, X_{\mathcal{Z}1}\times X_{\mathcal{Z}5},\\ X_{\mathcal{Z}2}\times X_{\mathcal{Z}3}, X_{\mathcal{Z}2}\times X_{\mathcal{Z}4}, X_{\mathcal{Z}2}\times X_{\mathcal{Z}5}, X_{\mathcal{Z}3}\times X_{\mathcal{Z}4},\\ X_{\mathcal{Z}3}\times X_{\mathcal{Z}5}, X_{\mathcal{Z}4}\times X_{\mathcal{Z}5}, X_{\mathcal{C}3}\times X_{\mathcal{Z}1})\end{aligned}$

c]@l@Non-linear main effects
and linear interactions (NL-L)     c]@c@$\begin{aligned}[t](X_{\mathcal{C}1}(X_{\mathcal{C}1}>0), \exp\{X_{\mathcal{C}2}\}^{-1/2}, |X_{\mathcal{C}3}|,  \log|X_{\mathcal{C}4}+1|, |X_{\mathcal{C}5}|^{-1/2},\\ X_{\mathcal{Z}1}, X_{\mathcal{Z}2}, X_{\mathcal{Z}3},  X_{\mathcal{Z}4}^2, \log|X_{\mathcal{Z}5}|,\\ X_{\mathcal{Z}1}\times X_{\mathcal{Z}2}, X_{\mathcal{Z}1}\times X_{\mathcal{Z}3}, X_{\mathcal{Z}1}\times X_{\mathcal{Z}4}, X_{\mathcal{Z}1}\times X_{\mathcal{Z}5},\\ X_{\mathcal{Z}2}\times X_{\mathcal{Z}3},  X_{\mathcal{Z}2}\times X_{\mathcal{Z}4}, X_{\mathcal{Z}2}\times X_{\mathcal{Z}5}, X_{\mathcal{Z}3}\times X_{\mathcal{Z}4},\\ X_{\mathcal{Z}3}\times X_{\mathcal{Z}5}, X_{\mathcal{Z}4}\times X_{\mathcal{Z}5}, X_{\mathcal{C}3}\times X_{\mathcal{Z}1})\end{aligned}$

c]@l@Non-linear main effects
and non-linear interactions (NL-NL)     c]@c@$\begin{aligned}[t](X_{\mathcal{C}1}(X_{\mathcal{C}1}>0), \sqrt{\exp\{X_{\mathcal{C}2}\}}, |X_{\mathcal{C}3}|,  \log|X_{\mathcal{C}4}+1|, \sqrt{|X_{\mathcal{C}5}|},\\ X_{\mathcal{Z}1}, X_{\mathcal{Z}2}, X_{\mathcal{Z}3},  X_{\mathcal{Z}4}^2, \log|X_{\mathcal{Z}5}|,\\ X_{\mathcal{C}1}(X_{\mathcal{C}1}>0)\times\sqrt{\exp\{X_{\mathcal{C}2}\}}, (X_{\mathcal{C}1}(X_{\mathcal{C}1}>0)\times |X_{\mathcal{C}3}|, \\ X_{\mathcal{C}1}(X_{\mathcal{C}1}>0)\times\log|X_{\mathcal{C}4}+1|,  X_{\mathcal{C}1}(X_{\mathcal{C}1}>0)\times \sqrt{|X_{\mathcal{C}5}|},\\ \sqrt{\exp\{X_{\mathcal{C}2}\}}\times |X_{\mathcal{C}3}|,  \sqrt{\exp\{X_{\mathcal{C}2}\}}\times \log|X_{\mathcal{C}4}+1|, \\ \sqrt{\exp\{X_{\mathcal{C}2}\}}\times\sqrt{|X_{\mathcal{C}5}|}, |X_{\mathcal{C}3}|\times \log|X_{\mathcal{C}4}+1|, \\ |X_{\mathcal{C}3}|\times \sqrt{|X_{\mathcal{C}5}|}, \log|X_{\mathcal{C}4}+1|\times \sqrt{|X_{\mathcal{C}5}|}\\ |X_{\mathcal{C}3}|\times X_{\mathcal{Z}1}\end{aligned}$
$^\ast$ The baseline scenario.     1l

Design matrix for the treatment generating models of various degrees of nonlinearity and/or nonadditivity.

1l                         3cc]@c@Bias
$\times$1000     3cc]@c@Empirical SD
$\times$1000     3cc]@c@SE (usual)
$\times$1000     3cc]@c@Coverage (%)
(usual)         3cc]@c@SE (modified)
$\times$1000     3cc]@c@Coverage (%)
(lr)2-4 (lr)5-7 (lr)8-10 (lr)11-13 (lr)14-16 (lr)17-19
1lEstimators                1 vs 2                       1 vs 3                      2 vs 3                      1 vs 2                         1 vs 3                         2 vs 3                         1 vs 2                         1 vs 3                        2 vs 3                        1 vs 2                          1 vs 3                         2 vs 3                         1 vs 2                          1 vs 3                         2 vs 3                         1 vs 2                          1 vs 3                         2 vs 3

1lOAL             95                             70                             -25                            37                             34                             36                             39                             35                             37                             96.0                            94.9                           94.5                           38                              34                             36                             94.9                            94.9                           94.9

All                                          72                             56                             -16                            37                             37                             35                             40                             37                             37                             55.5                            65.9                           92.7                           41                              40                             38                             59.8                            74.2                           94.6

Ysel                                         -4                             -2                             2                              38                             35                             36                             48                             41                             44                             98.2                            97.3                           97.3                           39                              35                             37                             95.4                            95.0                           94.9

YZsel                                        -3                             -1                             2                              39                             37                             38                             46                             40                             43                             97.4                            96.0                           95.9                           40                              36                             38                             95.0                            95.3                           95.4

OP+All                                       14                             13                             -2                             37                             35                             35                             40                             36                             38                             95.3                            94.1                           95.8                           43                              40                             41                             96.6                            97.6                           97.2

OP+Ysel                                      -4                             -2                             2                              38                             35                             36                             48                             41                             44                             98.2                            97.3                           97.3                           39                              35                             37                             95.4                            95.0                           94.9

OP+YZsel                                     -3                             -1                             2                              37                             35                             36                             45                             39                             42                             97.7                            96.8                           97.0                           38                              34                             37                             94.8                            95.0                           94.8

All                                          112                            82                             -30                            58                             58                             55                             69                             69                             65                             65.8                            81.6                           95.7                           57                              58                             55                             51.0                            70.8                           91.9

Ysel                                         61                             43                             -18                            50                             50                             47                             63                             62                             59                             87.6                            93.5                           98.4                           52                              51                             50                             78.1                            86.8                           94.9

YZsel                                        51                             38                             -13                            48                             47                             46                             60                             59                             56                             91.9                            95.3                           97.4                           50                              49                             48                             84.3                            89.4                           95.1

OP+All                                       7                              4                              -3                             55                             52                             56                             70                             65                             71                             98.2                            98.5                           98.3                           56                              52                             56                             94.0                            94.4                           94.4

OP+Ysel                                      13                             7                              -6                             45                             44                             45                             61                             57                             61                             98.8                            98.6                           99.2                           46                              44                             46                             94.0                            95.1                           95.3

OP+YZsel                                     17                             11                             -6                             43                             41                             43                             57                             54                             58                             98.4                            98.4                           99.4                           44                              41                             44                             93.5                            95.0                           94.7

1lPruned CART

All                                          125                            94                             -31                            45                             43                             38                             68                             68                             63                             53.0                            80.4                           99.1                           42                              42                             39                             19.4                            40.2                           88.6

Ysel                                         86                             64                             -22                            44                             43                             38                             62                             61                             57                             78.6                            88.6                           98.5                           42                              42                             39                             48.4                            66.0                           92.8

YZsel                                        75                             56                             -19                            46                             45                             39                             60                             58                             55                             80.6                            90.7                           97.8                           43                              42                             40                             60.0                            73.4                           93.9

OP+All                                       -4                             -1                             3                              39                             36                             38                             68                             62                             68                             99.9                            100.0                          99.9                           40                              36                             39                             94.0                            95.0                           95.9

OP+Ysel                                      -1                             -1                             1                              36                             35                             37                             58                             54                             59                             99.9                            99.3                           99.5                           38                              34                             37                             94.6                            94.9                           95.6

OP+YZsel                                     1                              0                              0                              37                             35                             37                             55                             52                             56                             99.1                            99.5                           99.5                           38                              35                             37                             94.8                            95.4                           95.2

1lBagged CART

All                                          109                            82                             -27                            38                             38                             36                             31                             31                             29                             10.1                            28.2                           79.1                           40                              40                             38                             24.4                            46.9                           91.1

Ysel                                         -1                             0                              1                              52                             46                             48                             39                             36                             34                             84.8                            87.0                           85.0                           52                              48                             49                             94.4                            95.8                           96.0

YZsel                                        -51                            -40                            11                             66                             58                             59                             46                             41                             39                             69.8                            74.8                           81.0                           63                              57                             60                             85.8                            86.9                           95.2

OP+All                                       -4                             1                              4                              37                             34                             36                             62                             62                             68                             99.9                            99.9                           100.0                          38                              33                             37                             94.7                            95.0                           94.7

OP+Ysel                                      -7                             -2                             5                              34                             32                             34                             57                             60                             64                             99.9                            99.9                           100.0                          36                              32                             35                             94.3                            95.6                           94.7

OP+YZsel                                     -7                             -2                             5                              34                             32                             34                             55                             57                             60                             99.3                            100.0                          99.9                           36                              32                             35                             94.5                            95.2                           94.6

2lRandom   Forests

All                                          126                            95                             -30                            34                             35                             34                             28                             28                             26                             2.5                             13.2                           73.1                           37                              38                             35                             7.0                             27.8                           88.2

Ysel                                         32                             25                             -7                             44                             39                             37                             30                             29                             27                             69.3                            77.2                           84.5                           43                              41                             41                             86.2                            91.4                           96.7

YZsel                                        -71                            -49                            22                             86                             65                             63                             34                             32                             29                             48.3                            57.3                           64.3                           62                              52                             58                             73.2                            78.6                           93.5

OP+All                                       -4                             1                              4                              38                             35                             37                             43                             39                             42                             96.6                            96.9                           96.4                           39                              34                             38                             95.3                            95.2                           94.0

OP+Ysel                                      -8                             -2                             6                              35                             32                             35                             50                             47                             52                             98.6                            99.5                           98.7                           36                              32                             36                             94.7                            95.2                           94.7

OP+YZsel                                     -8                             -2                             6                              35                             32                             35                             56                             54                             59                             99.5                            99.8                           99.8                           36                              32                             36                             94.8                            95.6                           94.7

[Standard errors (SE) and coverage of 95% confidence intervals estimated by usual bootstrap and modified bootstrap for sample size of 1000.]Standard errors (SE) and coverage of 95% confidence intervals estimated by usual bootstrap and modified bootstrap for sample size of 1000. The scenario with sparse treatment models was considered. Results were obtained based on 1000 simulated datasets. For each dataset, 200 bootstrap samples were generated.

1l     3cBias$\times 1000$     3cMCSD$\times 1000$     3cRMSE$\times 1000$     3cSE$\times 1000$     3cCoverage (%; modified)
(lr)2-4 (lr)5-7 (lr)8-10 (lr)11-13 (lr)14-16
1lEstimators     1l1 vs 2     1 vs 3     2 vs 3     1 vs 2     1 vs 3     2 vs 3     1 vs 2     1 vs 3     2 vs 3     1 vs 2     1 vs 3     2 vs 3     1 vs 2     1 vs 3     2 vs 3
1lNaive     165     124     -41     50     52     49     172     135     64     50     52     48     9.7     34.1     85.5

1lOAL     5     5     0     56     50     52     57     50     52     55     50     52     94.1     95.6     94.2

Confounder     -1     1     2     58     53     55     58     53     55     57     52     54     94.3     94.3     94.3

Treatment     2     3     1     73     63     65     73     63     65     69     61     62     94.1     94.2     93.4

Outcome     -1     0     2     57     49     53     57     49     53     56     50     53     94.5     96.1     94.1

All     95     73     -22     52     51     50     108     89     55     55     55     52     59.8     75.9     93.9

Ysel     -6     -3     3     57     49     54     57     49     54     58     51     55     95.4     96.0     94.8

YZsel     -6     -3     3     58     51     54     58     52     54     58     53     55     95.0     95.5     95.0

OP + All     16     13     -2     51     46     50     54     48     50     58     56     56     96.3     97.8     96.9

OP + Ysel     -6     -3     3     57     46     54     57     59     54     58     51     55     95.4     96.0     94.8

OP + YZsel     -5     -2     2     55     48     53     56     48     53     55     49     53     95.2     96.0     94.3

Confounder     55     41     -14     70     66     65     89     78     66     69     68     66     86.7     90.6     94.6

Treatment     94     71     -24     72     72     66     119     101     70     73     72     68     73.8     82.7     93.6

Outcome     68     51     -17     71     67     65     99     84     67     71     70     68     81.0     89.1     94.8

All     123     93     -30     79     81     76     146     123     81     78     80     74     63.2     75.9     92.0

Ysel     72     54     -17     72     71     66     102     89     70     72     72     69     81.6     88.1     94.0

YZsel     60     44     -15     69     66     63     91     79     66     70     68     67     86.2     90.8     94.8

OP + All     10     8     -1     80     74     79     80     76     79     79     74     79     93.4     93.8     94.2

OP + Ysel     12     9     -3     67     61     65     68     62     65     66     63     66     94.2     94.4     94.8

OP + YZsel     15     9     -6     62     58     62     64     59     62     63     59     62     94.3     94.7     94.8

1lPruned CART

Confounder     81     61     -20     64     61     55     103     86     59     61     60     57     69.4     81.4     93.3

Treatment     110     83     -27     66     66     58     128     106     64     65     64     59     58.1     72.4     92.3

Outcome     92     69     -22     63     59     54     111     91     59     60     60     56     63.1     77.0     93.0

All     133     100     -32     62     63     56     147     118     66     61     61     56     40.1     60.4     90.4

Ysel     96     72     -24     63     60     55     115     94     60     60     60     56     61.0     75.5     92.4

YZsel     84     63     -20     64     61     55     105     87     59     61     60     57     68.6     80.7     93.0

OP + All     -2     1     2     58     53     56     58     52     56     58     52     57     94.3     94.5     94.6

OP + Ysel     1     2     1     56     49     54     56     49     54     55     50     54     94.2     94.8     94.6

OP + YZsel     2     2     0     54     50     53     54     50     53     55     50     54     95.7     94.9     94.7

1lBagged CART

Confounder     -50     -37     12     93     83     86     105     91     87     88     80     84     90.1     90.0     93.8

Treatment     43     33     -10     81     76     71     92     83     72     79     74     71     89.0     92.0     93.9

Outcome     -2     0     2     76     65     70     76     65     70     75     70     71     94.2     96.4     95.0

All     119     90     -30     54     53     51     131     104     59     56     57     53     44.1     65.5     91.2

Ysel     17     14     -3     71     63     63     73     65     63     71     66     67     93.4     95.6     95.2

YZsel     -34     -27     7     91     82     81     97     85     82     86     79     81     91.6     92.5     94.8

OP + All     -4     -1     4     53     46     52     53     46     52     52     47     51     94.2     95.4     94.4

OP + Ysel     -7     -2     5     51     45     49     51     45     50     50     45     50     94.6     95.2     95.0

OP + YZsel     -7     -2     5     51     45     49     51     45     49     50     46     50     94.6     95.0     94.6

1lRandom Forests

Confounder     -22     -13     9     77     67     71     80     68     71     75     68     72     94.7     95.1     95.0

Treatment     62     48     -14     62     58     57     87     75     58     63     61     59     82.6     89.4     94.8

Outcome     35     28     -7     58     52     54     67     59     54     62     59     59     92.3     95.7     96.2

All     135     102     -33     49     51     48     144     114     58     52     53     50     27.1     51.2     90.0

Ysel     53     41     -11     57     52     51     77     66     53     59     57     57     84.9     91.2     94.5

YZsel     -36     -24     13     94     78     79     101     81     82     79     70     75     91.0     91.2     94.0

OP + All     -5     -1     4     55     47     53     55     47     52     53     48     52     94.8     94.9     94.0

OP + Ysel     -9     -2     6     51     45     50     52     45     50     51     46     50     94.5     94.8     94.4

OP + YZsel     -8     -2     5     50     46     50     51     45     50     51     46     50     94.8     95.0     94.7

[Simulation results for the scenario with sparse models and sample size of 500.]Simulation results for the scenario with sparse models and sample size of 500. Standard errors were estimated based on 200 bootstrap replications. Results were obtained using 2000 simulated datasets.

1l     3cBias$\times 1000$     3cMCSD$\times 1000$     3cRMSE$\times 1000$     3cSE$\times 1000$     3cCoverage (%; modified)
(lr)2-4 (lr)5-7 (lr)8-10 (lr)11-13 (lr)14-16
1cEstimators     1l1 vs 2     1 vs 3     2 vs 3     1 vs 2     1 vs 3     2 vs 3     1 vs 2     1 vs 3     2 vs 3     1 vs 2     1 vs 3     2 vs 3     1 vs 2     1 vs 3     2 vs 3
1lNaive     165     125     -41     36     37     33     169     130     53     36     37     34     0.5     8.3     78.4

1lOAL     -2     0     2     38     33     36     38     33     36     38     34     36     94.9     94.9     94.9

Confounder     -3     0     2     40     36     37     40     36     37     40     36     38     94.5     95.0     95.2

Treatment     -2     0     1     50     43     43     50     43     43     49     43     43     94.2     94.7     94.4

Outcome     -3     0     2     38     33     36     38     33     36     38     34     37     94.8     95.1     94.9

All     71     55     -16     38     36     35     81     66     38     41     40     38     59.8     74.2     94.6

Ysel     -6     -3     3     38     34     36     39     34     36     39     35     37     95.4     95.0     94.9

YZsel     -5     -2     3     40     36     37     40     36     37     40     36     38     95.0     95.3     95.4

OP + All     12     11     -1     37     33     35     39     35     35     43     40     41     96.6     97.6     97.2

OP + Ysel     -6     -3     3     38     34     36     39     34     36     39     35     37     95.4     95.0     94.9

OP + YZsel     -5     -3     3     38     34     36     39     34     36     38     34     37     94.8     95.0     94.8