Propensity Score Estimation Using Classification and Regression Trees in the Presence of Missing Covariate Data

Counterfactual outcomes and estimating causal effects

We adopt a perspective of potential or counterfactual outcomes, formal accounts of which are given for example by Neyman et al. (1935), Holland (1986), Rubin (1974), Holland (1988), and Pearl (2009).

Consider a sequence S=(X1,X2,...,Xn) of variables and let F=(fX1,fX2,...,fXn) be a collection of functions fXj that deterministically map a realisation of the predecessors (Xi:i<j) of Xj and of exogenous variable εXj into a realisation of Xj. We may write the random variable Xj as follows:

For any intervention setting Xt=xt for t in a subset T of {1,...,n}\{j}, the counterfactual version of Xj is obtained by evaluating the right-hand side of eq. (1) with Xt replaced by xt for all t∈T. Specifically, let S = (W,A,Y,R) so that W=fW(εW), A=fA(W,εA), Y=fY(W,A,εY), and R=fR(W,A,Y,εR). W may be thought of as a (random vector of) baseline or pre-exposure variable(s), A denotes the binary exposure of interest, Y the outcome, and R a missing indicator vector of W. A subject’s counterfactual outcomes Y0 and Y1, obtained if exposure A were set possibly contrary to fact to 0 and 1, respectively, are defined such that Y0=fY(W,0,εY) and Y1=fY(W,1,εY).

Causal effects are readily defined in terms of counterfactual outcomes. In this article, the focus is on the causal odds ratio (OR) for the marginal effect of exposure A on binary outcome Y among the exposed (A = 1), that is,

Under consistency as defined by Cole and Frangakis (2009), Y1 is equal to the observed outcome Y for subjects in the exposure group. Y0, on the other hand, is not observed for exposed subjects. We may, however, validly estimate the causal OR under a set of conditions, which includes no interference between subjects (or Stable Unit Treatment Value Assumption, Tchetgen and VanderWeele, 2012), consistency, positivity, and conditional exchangeability (Lesko et al. 2017). To simplify arguments and notation, we shall assume that all of these conditions hold, with the exception of conditional exchangeability, unless otherwise indicated. If there exists a (set of) variable(s) Z such that the potential outcomes are conditionally independent of exposure status given Z, we may write

so that the causal OR may be expressed entirely in terms of observables

W satisfies the definition of Z whenever εY⊥⊥εA|W. In practice, validly estimating E[Y|A=0,Z] may be difficult when Z is multidimensional and Y is rare (Albert and Anderson 1984). In this case, it may be desirable to summarise Z in a single balancing score (Rosenbaum and Rubin 1983).

The propensity score

The propensity score e(W), defined as the conditional probability of exposure given covariates W, satisfies a number of balancing properties. First, covariate(s) W and exposure A are conditionally independent given the propensity score, and conditional exchangeability given covariate(s) W implies conditional exchangeability given e(W) (Rosenbaum and Rubin 1983, Theorems 1 and 3). Thus, the causal OR becomes

This formulation has motivated the propensity score matching approach as discussed by Rosenbaum and Rubin (1983).

Balance may also be attained by inverse probability weighting (Appendix A). To simplify arguments and notation, we assume that W and Y take a discrete joint distribution; however, the results extend to continuous or mixed discrete/continuous distributions. To obtain an IPW estimator of the ATT, let

for realisations w of W and a of A, where I denotes the indicator function taken the value 1 if the argument is true and 0 otherwise. Weighting by φ yields independence between covariate(s) W and A; that is, for all w,

Also, conditional exchangeability given W implies exchangeability following weighting by ϕ; that is, if (Y0,Y1)⊥⊥A|W=w for all w, then

for all y0,y1. Thus, the causal OR becomes

In words, this means that the causal odds ratio is equal to the crude odds ratio of the ATT in the (pseudo-)population that is obtained by weighting each observation by φ.

Ensemble CART methods in the absence of missing data

We will now briefly describe how CART can be applied to estimate the propensity score. Detailed information can be found elsewhere (Breiman 1996; Ridgeway 1999; Breiman 2001; McCaffrey et al. 2004; Elith et al. 2008; Hastie et al. 2009). CART is a type of supervised learning task that entails finding a set of rules, subject to constraints, that partition the data into regions based on the input data (covariates) such that within regions, target values (e.g., exposure levels) meet a desirable level of homogeneity. Typically, a tree is built in a recursive manner by splitting the dataset into increasingly homogeneous subsets and choosing the splitting rule at each step or node that best splits the data further, with ‘best’ referring to the greatest improvement in terms of some homogeneity metric, such as the Gini index (Therneau and Atkinson 2017). Ensemble techniques by definition fit more than one tree to the data and combine them to form a single predictor of ‘the outcome’ (in the case of propensity score, the assigned exposure). The aim of ensemble techniques is to enhance performance and reduce issues of overfitting by a single tree (Elith et al. 2008; 2008; Hastie et al. 2009). We focus here on two popular CART ensemble methods, namely boostrap aggregated (bagged) CART and boosted CART.

Ignorable missing data

Suppose W=(W1,W2,...,Wp) and R=(R1,R2,...,Rp) are random vectors of length p such that for j = 1,2,...,p, Rj=0 if Wj is missing and Rj=1 if Wj is observed. Following Rubin (1976), define the extended random vector V=(V1,V2,...,Vp) with range to include the special value ∗ to indicate a missing datum: Vj=Wj if Rj=1 and Vj=∗ if Rj=0. Let v be a particular sample realisation of V, so that each vj is either a known quantity or ∗. These values imply a realisation for the random variable R, denoted r. For notational convenience, we write W=(Wobs,Wmis) and V=(Vobs,Vmis) to indicate that each may be partitioned into two vectors corresponding to all j such that rj=1 for observed data and rj=0 for missing data. It is important to note that these partitions are defined with respect to r, the observed pattern of missing data. Given a realisation r of R, and provided that A,Y are observed, covariate data are said to be missing at random (MAR) if Pr(R=r|Wobs,Wmis=u,A,Y) and Pr(R=r|Wobs,Wmis=u′,A,Y) are the same for all u,u^ʹ and at each possible value of the parameter vector φ that fully characterises the missing data mechanism (Rubin 1976). If in addition to MAR, the parameter φ is distinct, in the sense of Rubin (1976), from the vector θ that parameterises the distribution of the data that we would have based inference on had there been no missingness, then missing data is said to be ignorable and it is not necessary to consider the missing data or the missing data mechanism in making inferences about θ (Rubin 1976, 1987; Schafer 1997). Thus, if the missing data mechanism is ignorable, one may validly model the complete data to create imputations for the missing data (Rubin 1987; Van Buuren 2012).

The generalised propensity score

The generalised propensity score e∗(V) is defined as the conditional exposure probability given the extended covariate vector V (D’Agostino and Rubin 2000). That is,

Using the same argumentation to establish the balancing properties of the usual propensity score, it can be shown that the generalised propensity score has the same balancing properties with respect to V as the usual propensity score has with respect to W. Thus, the observed covariate data and missingness information and exposure A are conditionally independent given the generalised propensity score, and conditional exchangeability given the extended covariate(s) V implies conditional exchangeability given the generalised propensity score e∗(V).

To obtain an IPW estimator of the ATT, let

for realisations v of V and a of A, Then, weighting by γ renders V independent of A; that is, for all v,

Also, conditional exchangeability given V implies conditional exchangeability following weighting by γ; that is, if (Y0,Y1)⊥⊥A|V, then

for all y0,y1. Importantly, the propensity score e(W) need not equal the generalised propensity score e∗(V). That is, given the observed covariate data, the unobserved covariate data need not provide the same information about exposure allocation as does the missing data pattern. In addition, neither covariate balance given the propensity score (W⊥⊥A|e(W)) nor balance of the observed data and missingness information given the generalised propensity score (V⊥⊥A|e∗(V)) generally implies covariate balance given the generalised propensity score (W⊥⊥A|e∗(V)).

More crucially perhaps, conditional exchangeability given the generalised propensity score is not guaranteed even if conditional exchangeability given the usual propensity score holds or the generalised propensity score balances both observed and unobserved covariate data (i.e., neither (Y0,Y1)⊥⊥A|e(W) nor W⊥⊥A|e∗(V) nor both imply that (Y0,Y1)⊥⊥A|e∗(V); see Appendix B for an example).

This suggests that it is not generally desirable to distribute across exposure groups both the observed data and the missingness information by adjusting for the generalised propensity score. However, there are situations conceivable in which it is appropriate to base inference on the generalised rather than the usual propensity score. Until now, we have assumed an ordering of the variables in which the outcome Y precedes R, the missingness pattern of W. Consequently, Y was defined as a function fY of W, A and exogenous variable εY and not of R. Consider now a setting where S = (W,R,A,Y) so that R forms a predecessor of A and Y (and, therefore, a potential common cause of A and Y). Then, if exchangeability can be attained by conditioning on e∗(V), conditional exchangeability given e(W) need not hold (see Appendix C for an example).

Thus, the choice between adjustment for the generalised versus the usual propensity score should ideally rest on the relative extent to which conditional exchangeability holds given the generalised versus the usual propensity score. In practice, it is not possible to estimate directly the true propensity score when covariate data are missing (Rosenbaum and Rubin 1983; D’Agostino and Rubin 2000; Cham and West 2016). However, under ignorability of missing data, one may ‘recover’ the unobserved data, e.g., via multiple imputation (Rubin 1987; Van Buuren 2012), prior to estimating propensity scores. Henceforth, we assume that exchangeability can be attained by conditioning on the complete covariate data or, therefore, the usual propensity score, if data were not missing. We also assume that missing data is ignorable.

Applying CART to incomplete data

Simulation structure

We performed a series of Monte-Carlo simulation experiments based on the simulation structure described in Setoguchi et al. (2008) with modifications so as to allow for missing data. For n = 2000 subjects, we independently generated 10 covariates Wi (four confounders, three predictors of the exposure only, and three predictors of the outcome only), a binary exposure variable A, and a binary outcome Y (Figure 1). Missing data were introduced into one or two covariates. A number of CART-based approaches were used to estimate propensity scores, before and after the introduction of missing data, and in turn the log odds ratio for the exposure-outcome effect among the treated. For comparison, we also estimated propensity scores in imputed datasets using a correctly specified propensity score model, and using a logistic model with main effects only. The process was repeated 5000 times for each of eight simulation scenarios that varied primarily by missing data mechanism. All simulations were conducted with R-3.2.2 on a Windows 7 (64-bit) platform (R Core Team 2016).

Figure 1:

Complete data structure for simulation experiments. Dashed arcs without arrowheads connecting variables indicate non-zero entries for the corresponding variables in the covariance matrix of the joint distribution of all Wi, i = 1,2,...,10.

Data generation

Data were generated by sequentially going through the following steps. First, the covariates were generated by sampling from a multivariate normal distribution with zero means and unit variances; correlations were set to zero except for the correlations between W1 and W5, W2 and W6, W3 and W8, and W4 and W9, which were set to 0.2, 0.9, 0.2, and 0.9, respectively. Second, covariates W1, W3, W5, W6, W8, and W9 were dichotomised, setting any value to 1 if greater than 0 and to 0 otherwise.

Following Setoguchi et al. (2008), the binary exposure variable A was related to the covariate vector following the propensity score model Pr(A=1|W)=expit{0.8W1−0.25W2+0.6W3−0.4W4−0.8W5−0.5W6+0.7W7−0.25W22−0.4W42+0.7W72+0.4W1W3−0.175W2W4+0.3W3W5−0.28W4W6−0.4W5W7+0.4W1W6−0.175W2W3+0.3W3W4−0.2W4W5−0.4W5W6}. Realisations a for A were generated by drawing a pseudo-value from the uniform(0,1) distribution and setting a to 1 if this number was less than the true propensity score and to 0 otherwise. Consequently, A can be thought of as a exposure that is generated by a non-linear and non-additive propensity score model. This model assigns approximately half of subjects to the exposure group.

Outcomes were generated following the mechanism described by Setoguchi et al. (2008) with slight modifications to increase the outcome fraction (from approximately 2% to 20% or 40%). Specifically, the binary outcome, Y, was modelled as a Bernoulli random variable given A and W: an independent random number, εY, was drawn from the uniform distribution; Y was set to 1 if this number was less than the inverse logit (expit) of a linear transformation η(A,W)=−1+0.3W1−0.36W2−0.73W3−0.2W4+0.71W8−0.19W9+0.26W10+γA of A and W and to 0 otherwise. The true conditional log odds ratio for the exposure-outcome effect was set to 1 or −1 depending on the scenario. The outcome incidence was roughly 40% for scenarios with γ = 1; and 20% for scenarios with γ=−1. The counterfactual outcomes Y0 and Y1 for any subject with realisations w of W and u of εY are found by computing I(u < expit{η(0,w)}) and I(u < expit{η(1,w)}), I denoting the indicator function. With knowledge of the counterfactual outcomes, it can be inferred that with γ = 1, the marginal log odds ratio for the true exposure-outcome effect among the exposed (or treated; ATT) is approximately 0.906; with γ=−1, the marginal log odds ratio is approximately −0.926 (Hernán and Robins 2017). Note that these are different from the conditional causal odds ratios as a result of the non-collapsibility property of the odds ratio.

We considered ignorable missing data mechanisms for introducing missing data.

MCAR missingness. For all subjects, irrespective of complete data, values of W3 were set to missing with probability p, characterising the MCAR mechanism. The missingness probability of the other variables was set to zero.

MAR missingness. Let M3 be a missing indicator variable that takes the value of one if and only if the value of W3 is missing. Similarly, define M4 to be the missing indicator variable pertaining to W4. Given the full data, W3 and W4 were set to missing independently of one and other and with probability equal to Pr(M3=1|W,A,Y)=p and Pr(M4=1|W,A,Y)=expit{α0+α1W1+α2A+α3Y}. The missingness probability of the other variables was set to zero.

Scenarios

We evaluated the performance of various CART-based methods in eight scenarios (Table 1). The intercepts α0 in scenarios five through eight were chosen so as to yield roughly the same average proportion of missing data points per generated dataset of 24000 data points (2000 records on 10 covariates, one exposure and one outcome variable), namely 3%. The average proportion of missing data points and the fraction of incomplete records were largest in scenario 2 (5% and 60%, respectively; see Table 1). In all of the scenarios considered, data are ‘missing at random’ and it is assumed that there is conditional exchangeability given measured covariates (i.e., (Y0,Y1)⊥⊥A|W).

Note that in scenarios 3 trough 8, conditioning on M may break the independence between A and unobserved outcome predictor εY through what is known as collider stratification (cf. Pearl, 2009). One might therefore expect that that discarding incomplete records in these scenarios would result in bias. In scenarios 3 and 4, however, covariate missingness M is conditionally independent of exposure status and covariate data given the outcome (i.e., M⊥⊥(A,W)|Y). As a result, in these scenarios, the conditional OR for the effect of A on Y given W is equal to the conditional OR given W among the complete cases (Westreich 2012). Bias of complete case estimators in scenarios 3 and 4 therefore cannot be attributed to collider stratification, despite the presence of an unobserved outcome predictor. Instead, it could result from the non-collapsibility of the odds ratio and changes in the covariate distribution brought about by narrowing the focus of inference to the complete cases (Hernán and Robins 2017).

Estimators

Bagged CART was based on 100 bootstrap replicates (Lee et al. 2010). We imposed complexity constraints on the tree fitting algorithm using the rpart package default control settings. For boosted CART, we used 20000 iterations, a shrinkage parameter of 0.0005, and an iteration stopping rule based on the mean Kolmogorov-Smirnov test statistic (Lee et al. 2010; McCaffrey et al. 2004).

The CART methods were combined with several common approaches to handling missing data: leaving missingness information as is (i.e., subjecting incomplete data directly to the CART algorithm); complete case analysis (CCA); and multiple imputation (MI). MI was implemented with the mice package (version 2.46.0) using the logreg and norm options to impute missing binary and continuous variables, respectively, and otherwise default settings (Van Buuren and Groothuis-Oudshoorn 2011). Imputation models included, apart from the variable to be imputed, all other variables, including the outcome, as untransformed main effects only. Propensity score analysis was performed within imputed datasets using the respective sets of estimated propensity scores (Penning de Vries and Groenwold 2016) and results were combined using Rubin’s (1987) rules.

In addition to using CART, as stated, we also estimated propensity scores in imputed datasets using a correctly specified propensity score model (LRc), and using a logistic model with main effects only (LRm).

Within each (multiply imputed) dataset, the ATT was estimated from a logistic model with robust variance estimation using the survey package (Lumley 2014, version 3.31). We used both propensity score matching and inverse probability weighting. Matching was performed using a greedy 1:1 nearest neighbour algorithm, matching exposed (A = 1) to unexposed individuals (A = 0) (Austin 2011a). For any given (imputed) dataset, matching was performed on the logit of the propensity score, using a calliper distance of 20% of the standard deviation of the logit propensity score estimates (Austin 2011b). With the ATT as the estimand, IPW weights were defined as 1 for exposed subjects and PS/(1−PS) for unexposed subjects (PS denoting the estimated propensity score). To avoid undefined weights (1/0) or logit propensity scores (logit(0)), we placed bounds on the estimated propensity scores, truncating all estimates less than 0.001 to 0.001 and setting estimates greater than 0.999 to 0.999. MI-based estimates were pooled using Rubin’s rules to yield for each original dataset a single effect estimate, standard error estimate, and 90% confidence interval (90%CI).

Performance metrics

We evaluated the performance of the various methods through several measures: bias, estimated by the mean deviation of the estimated from the true marginal exposure-outcome effect on the log scale; empirical standard error; mean estimated standard error; mean squared error (MSE); and 90%CI coverage, estimated by the percentage of the 5000 data sets in which the 90%CI included the true exposure-outcome effect. Based on 5000 simulation runs, the Monte Carlo standard error for the true coverage probability of 0.90 is √(0.90(1−0.90)/5000)≈0.0042, implying that the estimated coverage probability is expected to lie with 95% probability between 0.893 and 0.907. Empirical coverage rates outside this interval provide evidence against the true coverage probabilities being equivalent to the nominal level of 0.90. The primary interest, however, was to gauge the effect of missing data on the various effect estimators. Therefore, we also compared, for each scenario, the effect estimates before and after the introduction of missing data.

Bias

Before the introduction of missing data, baCART and bCART showed small to no bias (with absolute values ranging from 0.000 to 0.011 on the log odds ratio scale). MI+LRc performed generally well and to a similar extent as MI+baCART and MI+bCART. MI+LRm consistently underestimated the true effect when inference was based on IPW; this trend was weaker for inference based on propensity score matching (Supplementary Table 1). Among all CART-based missing data approaches considered, multiple imputation yielded the least biased estimators overall (with a maximum absolute value of bias of 0.026 versus 0.221 and 0.138 for CART-only and CCA estimators, respectively), whereas bCART deviated on average the most from the true effect after the introduction of missingness.

As expected, baCART and bCART were biased (with −0.029 and −0.039, respectively, for scenario 1 and −0.064 and −0.072 for scenario 2) under MCAR in the direction of confounding by W3, whereas CCA and MI produced exposure-outcome effect estimates that were on average very close to the true effect. In scenarios 3 and 4, where missingness was outcome-dependent, bCART was biased toward the null after the introduction of missingness (with bias estimates of −0.117 and 0.064 for scenario 3 and 4, respectively, where the causal log odds ratios were approximately 0.906 and −0.926); baCART was downwardly biased in both scenarios (with bias estimates of −0.088 and −0.112). Estimators based on CCA or MI with CART were considerably less biased. In scenarios 5 through 8, CCA estimators systematically underestimated the true effect, particularly when the effect of the exposure or the outcome on the missingness probability was large (scenarios 6 and 7, where bias estimates ranged from −0.116 to −0.138). In these scenarios (5 through 8), baCART produced estimates that were on average close to the true effect, except in scenario 6, where the effect was clearly underestimated (estimated bias −0.050). bCART resulted in estimates that deviated in the same direction and to a similar or greater extent from the true effect as compared with CCA estimators. Again, MI with CART resulted in estimates that were on average close to the true effect. Increasing the effect of covariate W1 on the missingness probability (scenario 8 versus 5) had no evident impact on the results of any of the estimators.

Other performance

As expected, discarding incomplete records (CCA) resulted in relatively large empirical standard errors. Interestingly, MI+LRc had the largest empirical standard error in most scenarios, probably as a consequence of the complexity of the fitted propensity score models. In comparing empirical and mean estimated standard errors, note that multiple imputation produced generally conservative estimates of the standard error. This is consistent with previous observations (Van Buuren 2012). Among the CART-based estimators, the MSE was largest for CCA in nearly all scenarios. MI estimators had consistently small MSE. Overall, the best performance in terms of MSE was attained by MI estimators, followed by baCART and bCART. Multiple imputation with CART resulted in empirical coverage rates close to or slightly higher than the nominal 90% and those of the other estimators.