Outcome Adaptive Propensity Score Methods for Handling Censoring and High-Dimensionality: Application to Insurance Claims

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}\geq\min(T_{i},d)\}$ be the indicator of observing $Y_{i}$.

Under the potential outcome framework [19], 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^{\prime}$, denoted by $\tau(j,j^{\prime})=E\{Y^{(j^{\prime})}\}-E\{Y^{(j)}\}$. The following assumptions are required for statistical identification of causal effects using the observable data [20]:

1. I.

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

2. II.

   (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$;

3. III.

   (Unconfoundedness) $\{Y_{i}^{(1)},\cdots,Y_{i}^{(J)}\}\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}Z_{i}|\boldsymbol{X}_{i}$;

4. IV.

   (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})$;

5. V.

   (Censoring at random) $C_{i}\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}\{T_{i}^{(1)},\cdots,T_{i}^{(J)}\}\Big{|}(Z_{i},\boldsymbol{X}_{i})$.

We assume that the true outcome model for $Z_{i}=j$ is a logistic regression model,

$\displaystyle\text{logit}\,P(Y_{i}=1|\boldsymbol{W}_{i},Z_{i}=j)=\boldsymbol{W}_{i}^{T}\boldsymbol{\beta}_{j},$

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

$\displaystyle\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,$

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$,

$\displaystyle\lambda_{j}(t|\boldsymbol{U}_{i})=\lambda_{0j}(t)\exp(\boldsymbol{U}_{i}^{T}\boldsymbol{\gamma}_{j}),$

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,

$\displaystyle\hat{\tau}(j,j^{\prime})=\frac{\sum_{i=1}^{n}\hat{w}_{i}I(Z_{i}=j^{\prime})Y_{i}}{\sum_{i=1}^{n}\hat{w}_{i}I(Z_{i}=j^{\prime})}\,-\,\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)},$

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$,

$\displaystyle\hat{\tau}(j,j^{\prime})=\frac{\sum_{i=1}^{n}\hat{w}_{i}^{\ast}R_{i}I(Z_{i}=j^{\prime})Y_{i}}{\sum_{i=1}^{n}\hat{w}_{i}^{\ast}R_{i}I(Z_{i}=j^{\prime})}\,-\,\frac{\sum_{i=1}^{n}\hat{w}_{i}^{\ast}R_{i}I(Z_{i}=j)Y_{i}}{\sum_{i=1}^{n}\hat{w}_{i}^{\ast}R_{i}I(Z_{i}=j)}.$

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}$ [20], where $\Lambda_{ij}(t)$ is the cumulative hazard function of $C_{i}$ at $t$ for treatment $j$.

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 (All), 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

$\displaystyle\text{logit}P(Y_{i}=1|\boldsymbol{X}_{i})=\boldsymbol{X}_{i}^{T}\boldsymbol{\theta},$

(1)

the coefficients of which are estimated based on the LASSO penalty

$\displaystyle\hat{\boldsymbol{\theta}}=\operatorname*{arg\,min}_{\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],$

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,

$\displaystyle 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})},$

(2)

which is fitted by minimizing the negative penalized log-likelihood

$\displaystyle N^{-1}\sum_{i=1}^{N}\log P(Z_{i}=j|\boldsymbol{X}_{i})+\lambda\sum_{k=1}^{p}\|\boldsymbol{\psi}_{\cdot k}\|_{2},$

where $\boldsymbol{\psi}_{\cdot k}$ is a $J$-dimensional vector containing the $k$-th coefficient for all $J$ treatment groups [23]. 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 Ysel to label the set of variables selected by Model (1) and YZsel to label the intersection of the two sets of variables selected by Model (1) and Model (2). In an ideal case, YZsel is identical to $\boldsymbol{X}_{\mathcal{C}}$ and Ysel\YZsel 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 \scalebox{0.85}{3}⃝ in Figure 1 corresponds to the case where YZsel 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 [14] intends to achieve this by fitting an adaptive LASSO for the treatment [16], 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 [14]. 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, Ysel. The treatment can then be modeled as a function of Ysel (route \scalebox{0.85}{2}⃝ in Figure 1).

One possible problem of using Ysel as input for the final treatment model is that the dimensionality of Ysel 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 YZsel 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 \scalebox{0.85}{6}⃝ in Figure 1, and we denote the combination OP+YZsel). In addition to YZsel, as was noted later in the simulation studies, OP can also be combined with all the available covariates and Ysel as input for the final treatment model in Figure 1 to improve the effect estimates. We denote the combinations OP+All and OP+Ysel, which correspond to route \scalebox{0.85}{4}⃝ and route \scalebox{0.85}{5}⃝ in Figure 1, respectively. We discuss how OP can be used as predictors to estimate the propensity scores for different treatment modeling techniques in Section 4, where we outline several possible choices for the final treatment model.

In summary, both Ysel and YZsel 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 Ysel, YZsel, OP+Ysel, and OP+YZsel.

The methods listed in Table 1 can be applied to All, Ysel, and YZsel 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 (All, Ysel, or YZsel) 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 \scalebox{0.85}{1}⃝, \scalebox{0.85}{2}⃝, or \scalebox{0.85}{3}⃝:

$\displaystyle\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}.$

(3)

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 (3).

For OP+Ysel (route \scalebox{0.85}{5}⃝) and OP+YZsel (route \scalebox{0.85}{6}⃝), 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.

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 lambda.1se 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 All, which does not involve the pre-selection step, the lambda.min 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 nodesize 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 [14], 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 Ysel with additional weights imposed on the covariates for variable selection.

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^{\prime}$-dimensional vector of means $\boldsymbol{0}_{p^{\prime}}$ and a covariance matrix $\Sigma$, where $p^{\prime}$ 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,

$\displaystyle Z_{i}\sim\text{Multinomial}\{\pi_{1}(\boldsymbol{V}_{i}),\pi_{2}(\boldsymbol{V}_{i}),\pi_{3}(\boldsymbol{V}_{i})\},$

(4)

where $\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

$\displaystyle P\{Y_{i}^{(j)}=1|\boldsymbol{W}_{i}\}=\text{expit}\{\beta_{0j}+\boldsymbol{W}_{i}^{T}\boldsymbol{\beta}_{j}\},$

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 (4) 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

$\displaystyle\beta_{01}I(Z_{i}=1)+\beta_{02}I(Z_{i}=2)+\beta_{03}I(Z_{i}=3)+\boldsymbol{W}_{i}^{T}\boldsymbol{\beta}$

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 [26] 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,

$\displaystyle C_{i}^{(j)}=\{\lambda^{-1}\exp(\boldsymbol{U}_{i}^{T}\boldsymbol{\gamma}_{j})^{-1}\log u\}^{1/\nu},$

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 1 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.

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 All and OP+All 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 Ysel, YZsel, OP+Ysel, and OP+YZsel, 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 2) 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 OP+All 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\times 10^{5}$ replicates.

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

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 [20]. 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 1 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 YZsel and OP+YZsel ranged from 0.97 to 1, while the correlations between random forests based on YZsel and OP+YZsel 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 5 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 [6, 30]. 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 OP+YZsel. 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 5 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 [31].

As was observed in the simulation studies, bagged CART and random forests using YZsel 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 OP+All over All 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 All 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 OP+Ysel and OP+YZsel were very close across the propensity score estimation methods.