2D score-based estimation of heterogeneous treatment effects

2.1 Matching

To tackle the “missing data” problem when estimating treatment effects in randomized experiments in practice, matching serves as a very powerful tool. The main goal of matching is to find matched groups with similar or balanced observed covariate distributions [20]. The exact K -nearest-neighbour matching [2] is one of the most common and easiest to implement and understand methods; ratio matching [21–23], which finds multiple good matches for each treated individual, performs well when there is a large number of control individuals. Rosenbaum [24], Gu and Rosenbaum [25], and Zubizarreta [26,27] developed various optimal matching algorithms to minimize the total sum of distances between treated units and matched controls in a global sense. Abadie and Imbens [28] studied the consistency of covariate matching estimators under large sample assumptions. Instead of greedy matching on the entire covariates, propensity score matching (PSM) by Rubin and Thomas [29] is an alternative algorithm that does not guarantee optimal balance among covariates and reduces dimension sufficiently. Imbens [30] improved propensity score matching with regression adjustment. The additional matching on prognostic factors in propensity score matching was first considered by Rubin and Thomas [22]. Later, Leacy and Stuart [6] demonstrated the superiority of the joint use of propensity and prognostic scores in matching over single score-based matching in low-dimensional settings through extensive simulation studies. Antonelli et al. [7] extended the method to fit to high-dimensional settings and derived asymptotic results for the so-called doubly robust matching estimators (DRME). The sequential works by Aikens et al. [31] and Aikens and Baiocchi [32] pioneered in visualizing the relationship between propensity and prognostic scores as auxiliary in the estimation of treatment effects.

2.2 Subclassification

To understand heterogeneity of treatment effects in the data, subclassification, first used by Cochran [33], is another important research problem. The key idea is to form subgroups over the entire population based on characteristics that are either immutable or observed before randomization. Rosenbaum and Rubin [4,34], and Lunceford and Davidian [35] examined how creating a fixed number of subclasses according to propensity scores removes the bias in the estimated treatment effects, and Yang et al. [36] developed a similar methodology in settings with more than two treatment levels. Full matching [37–39] is a more sophisticated form of subclassification that selects the number of subclasses automatically by creating a series of matched sets. Schou and Marschner [40] presented three measures derived using the theory of order statistics to claim heterogeneity of treatment effect across subgroups. Su et al. [41] pioneered in exploiting standard regression tree methods [42] in subgroup treatment effect analysis. Further, Athey and Imbens [43] derived a recursive partition of the population according to treatment effect heterogeneity. Hill [44] was the first work to advocate for the use of BART [45] for estimating heterogeneous treatment effects, followed by a significant number of research papers focusing on the seminal methodology, including Green and Kern [46], Hill and Su [47], and Hahn et al. [8]. Abadie et al. [12] introduced endogenous stratification to estimate subgroup effects for a fixed number of subgroups based on certain quantiles of the prognostic score. More recently, Padilla et al. [48] combined the fused lasso estimator with score matching methods to lead to a data-adaptive subgroup effects estimator.

2.3 Machine learning for causal inference

For the goal of analyzing treatment effect heterogeneity, supervised machine learning methods play an important role. One of the more common ways for accurate estimation with experimental and observational data is to apply regression [49] or tree-based methods [50]. From a Bayesian perspective, Heckman et al. [51] provided a principled way of adding priors to regression models, and Taddy et al. [52] developed Bayesian nonparametric approaches for both linear regression and tree models. The recent breakthrough work by Wager and Athey [9] proposed the causal forest estimator arising from random forests from Breiman [53]. More recently, Athey et al. [10] took a step forward and enhanced the previous estimator based on generalized random forests. Imai and Ratkovic [54] adapted an estimator from the support vector machine classifier with hinge loss [55]. Bloniarz et al. [56] studied treatment effect estimators with lasso regularization [57] when the number of covariates is large, and Koch et al. [58] applied group lasso for simultaneous covariate selection and robust estimation of causal effects. In the meantime, a series of papers including Qian and Murphy [59], Kü nzel et al. [60], and Syrgkanis et al. [61] focused on developing meta-learners for heterogeneous treatment effects that can take advantage of various machine learning algorithms and data structures.

2.4 Applied work

On the application side, the estimation of heterogeneous treatment effects is particularly an intriguing topic in causal inference with broad applications in scientific research. Gaines and Kuklinski [62] estimated heterogeneous treatment effects in randomized experiments in the context of political science. Dehejia and Wahba [63] explored the use of propensity score matching for nonexperimental causal studies with application in economics. Dahabreh et al. [64] investigated heterogeneous treatment effects to provide the evidence base for precision medicine and patient-centred care. Zhang et al. [65] proposed the survival causal tree method to discover patient subgroups with heterogeneous treatment effects from censored observational data. Rekkas et al. [66] examined three classes of approaches to identify heterogeneity of treatment effect within a randomized clinical trial, and Tanniou et al. [67] rendered a subgroup treatment estimate for drug trials.

4.1 Step 1

We first estimate propensity and prognostic scores for all observations in the sample. For propensity score, we apply a logistic regression on the entire covariate X and the treatment indicator Z by solving the optimization problem

where σ ( x ) = 1 1 + exp ( − x ) is the logistic function. With the coefficient vector α ˆ , we compute the estimated propensity scores e ˆ by

For prognostic score, we restrict to the controlled group and regress the outcome variable Y on the covariate X through ordinary least squares: we solve

and we estimate prognostic scores as

4.2 Step 2

Next, we perform a nearest-neighbour matching based on the two estimated scores from the previous step. We adapt the notation from Abadie and Imbens [28], and use the Mahalanobis distance norm as the distance metric in the matching algorithm. Mahalanobis distance matching, first proposed by Hansen [69], has been widely used in matching algorithms for causal inference [6,31]. It has been found to perform well when the dimension of covariates are relatively low [20], and thus, it is more suitable than using standard Euclidean distance in our case when we have only propensity scores and prognostic scores as covariates.

Formally, for the units i and j with estimated propensity scores e ˆ i , e ˆ j and propensity scores p ˆ i , p ˆ j , we define the score-based Mahalanobis distance between i and j by

where Σ denotes the variance-covariance matrix of ( e ˆ , p ˆ ) ⊤ .

Let j k ( i ) be the index j ∈ { 1 , 2 , … , n } that solves Z j = 1 − Z i and

where 1 { ⋅ } is the indicator function. This is the index of the unit that is the k th closest to unit i in terms of the distance between two scores, among the units with the treatment opposite to that of unit i . We can now construct the K -nearest-neighbour set for unit i by the set of indices for the first K matches for unit i ,

We then compute

Intuitively, the construction of Y ˜ gives a proxy of the ITE on each unit. We find K matches for each unit in the opposite treatment group based on the similarity of their propensity and prognostics scores, and the mean of the K matches is used to estimate the unobserved potential outcome for each unit.

4.3 Step 3

The last step involves denoising of the point estimates of the individual treatment effects Y ˜ obtained from Step 2. The goal is to partition all units into subgroups such that the estimated treatment effects would be constant over some 2d intervals of propensity and prognostic scores (see the left of Figure 1).

Figure 1

Left: a hypothetical partition over the 2d space of propensity and prognostic scores with the true values of piecewise constant treatment effects; Right: a sample regression tree T constructed in Step 3.

To perform such stratification, we grow a regression tree on Y ˜ , denoted by T , and the regressors are the estimated propensity scores e ˆ and the estimated prognostic scores p ˆ from Step 1. We follow the very general rule of binary recursive partitioning to build the tree T : allocate the data into the first two branches, using every possible binary split on every covariate; select the split that minimizes Gini impurity, and continue the optimal splits over each branch along the covariate’s values until the minimum node size is reached. To avoid overfitting, we set the minimum node size as 20 in our model. Choosing other criteria such as information gain instead of Gini impurity is another option for splitting criteria. A 10-fold cross validation is also performed at meantime to prune the large tree T for deciding the value of cost complexity. Cost complexity is the minimum improvement in the model needed at each node. The pruning rule follows that if one split does not improve the overall error of the model by the chosen cost complexity, then that split is decreed to be not worth pursuing (see more details in Section 9.2 in the study by Hastie et al. [70]).

The final tree T (see the right plot of Figure 1) contains a few terminal nodes, and these are the predicted treatment effects for all units in the data. The values exactly represent a piecewise constant stratification over the 2d space of propensity and prognostic scores.

4.4 Estimation on a new unit

After we obtain the regression tree model T in Step 3, we can now estimate the value of the individual treatment effect corresponding to a new unit with covariate x new .

We first compute the estimated propensity and prognostic scores for the new observation by

where α ˆ and θ ˆ are the solutions to equations (5) and (6), respectively. Then with the estimated propensity score e ˆ new and prognostic score p ˆ new , we can obtain an estimate of the treatment effect for this unit following the binary predictive rules in the tree T .

4.5 High-dimensional estimator

In a high-dimensional setting where the covariate dimension d is much larger than the sample size n , we can estimate the propensity and prognostic scores by adding a lasso ( l 1 -based) penalty [57] instead. This strategy was first proposed and named as DRME by Antonelli et al. [7]. The corresponding optimization problems for the two scores can be written as follows:

The selection of the tuning parameters λ 1 and λ 2 can be determined by any information criteria (AIC, BIC, etc.). In practice, we use 10-fold cross-validation (CV) to select the value of λ . Then we perform the K -nearest-neighbour matching based on propensity and prognostic scores to obtain the estimates of individual treatment effects using equation (7).

We extend the aforementioned estimator with our proposal of applying a regression tree on the estimated propensity and prognostic scores. The procedure for estimation of subgroup heterogeneous treatment effects and estimation on a new unit remains the same as Step 3 for low-dimensional setups.

4.6 Computational complexity

Our method is composed of three steps as introduced before. We first need to implement a logistic regression for estimating propensity score for a sample with size n and ambient dimension d . To solve a logistic regression optimization problem involves using the proximal Newton method, and its computational complexity is of # of iterations ⋅ O ( n d ) [74]. In general, the algorithm converges after a small number of iterations, and we can safely assume it is of a constant order. As a result, the overall time complexity of solving a logistic regression is O ( n d ) . The estimation of prognostic scores also requires a computational cost of O ( n d ) , and for high-dimensional settings, the cost of solving the solution via coordinate descent remains at O ( n d ) [75]. The complexity of a K -nearest-neighbour matching based on the two estimated scores in the second step is of O ( K n ) [76], and the selection of K ≈ log ( n ) leads to a complexity of O ( n log n ) . In the third step, we grow a regression tree based on two estimated scores, and it requires a computational complexity of O ( n log n ) .

The overall computational complexity of our method depends on the comparison between the order of d and log ( n ) . Our method attains a computational complexity of O ( n d ) if the order of d is greater than that of log ( n ) . Otherwise, the complexity becomes O ( n log n ) .

5.1 Simulated data

We first examine on the following simulated data sets under six different data generation mechanisms. We obtain insights from the simulation study by Leacy and Stuart [6] for the models considered in Scenarios 1–3. The propensity score and outcome (prognosis) models in Scenarios 1 and 3 are characterized by additivity and linearity (main effects only), but with different piecewise constant structures in the true treatment effects over a 2d grid of the two scores. We add nonadditivity and nonlinear terms to both propensity and prognosis models in Scenarios 2. In other words, both propensity and prognostic scores are expected to be misspecified in these two models if we apply generalized linear models directly in estimation. Scenario 4 comes from Abadie et al. [12], with a constant treatment effect over all observations. Scenario 5 is modified from the study by Wager and Athey [9] (see Equation 27 there), in which the propensity model follows a continuous distribution instead of a linear structure. In addition, the true treatment effect is not a function of propensity and prognostic scores. In Scenario 6, we break the mold assumption of nicely delineated subgroups in treatment effects over 2d space of the two scores, and use a smooth function instead. A high-dimensional setting ( d ≫ n ) is examined in Scenario 7, where the generative model inherits from Scenario 1. We also include a table summarizing all simulation scenarios in Appendix B to help the audience interpret the result of the simulation study more readily.

We first introduce some notations used in the experiments: the sample size n , the ambient dimension d , as well as the following functions:

Throughout all the models we consider, we maintain the unconfoundedness assumption discussed in Section 3, generate the covariate X following a certain distribution, and entail homoscedastic Gaussian noise ε .

We evaluate the accuracy of an estimator τ ( X ) by the mean-squared error for estimating τ * ( X ) at a random example X , with m = n ∕ 10 , defined by

We record the averaged mean squared error (MSE) over 1,000 Monte Carlo trials for each scenario. For each Monte Carlo trial, we use train-test split to evaluate model performance. In detail, we first obtain n training sample and train models on this set. We then generate a separate testing sample under the same data generation mechanism with sample size n ∕ 10 . We make predictions on the test data using the trained models and compute the MSEs from different models accordingly. In terms of uncertainty quantification, we measure the coverage probability of τ ( X ) with a target coverage rate of 0.95. For endogenous stratification and our proposed method, we use nonparametric bootstrap with replacement to construct the empirical quantiles for each unit. The details on the implementation of nonparametric bootstrap methods are presented in Appendix C. For causal forest, we construct 95% confidence intervals by estimating the standard errors of estimation using the “grf” package. We do not include the results for Scenario 7, a high-dimensional setting, because the asymptotic normality guarantees of nonparametric bootstrap and random forests are not valid in high dimensions [9,79].

Scenario 1. With d ∈ { 2 , 10 , 50 } , n ∈ { 1,000 , 5,000 } , for i = 1 , … , n , we generate the data as follows:

where β e and β p are randomly generated from U [ − 1 , 1 ] d in each iteration.

Scenario 2. We now add some interaction and nonlinear terms to the propensity and prognostic models in Scenario 1, while keeping the setups of the covariate X , the response Y , the true treatment effect τ * , and the error term ε unchanged. We set d = 10 and n = 5,000 in this case.

Again for each simulation, β e and β p are randomly generated from U [ − 1 , 1 ] d .

Scenario 3. In this case, we define the true treatment effect with a more complicated piecewise constant structure over the 2d grid, under the same model used in Scenario 1, with d = 10 and n = 5,000 :

Scenario 4. Setting d = 10 and n = 5,000 , the data are generated as follows:

where β = ( 1 , … , 1 ) ⊤ ∈ R d . Moreover, the treatment indicators for the simulations are such that ∑ i Z i = ⌈ n ∕ 2 ⌉ . By construction, the vector of treatment effects satisfies τ * = 0 .

Scenario 5. The data satisfy

where e 1 = ( 1 , 0 , … , 0 ) ⊤ . We compare the results of different methods under the setting: d = 10 , n = 5,000 .

Scenario 6. For this scenario, we examine the performance when the treatment effect takes on a smooth interval of propensity and prognostic scores. We use the same model in Scenario 1 for generating the two scores with d = 10 and n = 5,000 , and we inherit and modify the generative function (28) from the study by Wager and Athey [9] for treatment effects as follows:

Scenario 7. In the last case, we study the performance of different estimators on a high-dimensional data. The data model follows

We select n = 3,000 and d = 5,000 for analysis.

The boxplots that depict the distribution of MSEs obtained under all scenarios are presented in Figure 2. We can see that for Scenario 1, our proposed estimator achieves better accuracy when the sample size n is large, and it is the best among the three estimators in these cases. The good performance of our method is consistent when we assume a more complex partition on the defined 2d grid as in Scenario 3. In addition, variation in accuracy, measured by the difference between the upper quartiles and the lower quartiles (also referred as the interquartile ranges) in each boxplot becomes smaller, accompanied with the increase in n . The reason for the enhanced performance when applying our method on a large-size data is the nature of nearest neighbour algorithm, which benefits more in learning the boundaries of clusters accurately with a larger sample size. Moreover, causal random forest suffers more in accuracy especially as d gets larger [9], while our method does not impact significantly because it focuses on a low-dimensional space of the two scores instead of the overall covariate space. In Scenario 2, we introduce nonadditivity and nonlinear terms into the data model. Although linear assumptions are violated for both propensity and prognostic models, our method performs better compared to the other two methods regarding accuracy and variability. For a potential outcome model with randomized assignment of treatment and constant treatment effects, as in Scenario 4, our method has the best accuracy among all candidates, even though large noise is added to the true signal. In Scenario 5, our estimator performs slightly worse than causal forests in accuracy since the true treatment function is not directly dependent on propensity and prognostic scores, but it still outperforms endogenous stratification. However, for the previous two cases, our method shows larger variation in accuracy than the others. One of the reasons for this phenomenon is due to the misspecification in estimating the true scores and the large variance in constructing tree models. When the assumption of a sharp stratification of treatment effects in the 2d space of the two scores is invalid, as shown in Scenario 6, our method maintains its superiority over the benchmarks. With subgroups identified over the 2d interval of the scores, the mean squared errors are reduced through our estimator comparably better than the other methods. In a high-dimensional setting such as Scenario 7, we consider modified methods with lasso-regularized regressions for both our methodology and endogenous stratification, and our method remains its high performance as in the low-dimensional setups since the piecewise constant structure in treatment effects in the setting does not change.

Figure 2

Comparison of mean squared errors over 1,000 Monte Carlo simulations under different generative models for our proposed method (PP), endogenous stratification (ST), and causal forest (CF).

We now take a careful look at the visualization comparison between the true treatment effect and the predictions obtained from our method for Scenarios 1, 2, 3, 6, and 7. We confine both the true signal and the predictive model in a 2d grid scaled on the true propensity and prognostic scores, as shown in Figure 3. It is not surprising that our proposed estimators provide a descent recovery of the piecewise constant partition in the true treatment effects over the 2d grid as in Scenarios 1, 2, 3, and 7, with only a small difference in the magnitude of treatment effects. In the setting when treatment effects are smoothly distributed across the grid of the two scores as in Scenario 6, our estimation does not exactly capture the soft boundary of the shape, but still provides us with a simplified delineation that help us understand general heterogeneity in treatment effects over the 2d space.

Figure 3

One instance comparison between the true treatment effects and the estimates from the score-based method for Scenarios 1, 2, 3, 6, and 7 (from the top to the bottom). Plots on the left column depict the true treatment effect on the 2d grid of propensity and prognostic scores, and the right plots demonstrate the estimated treatment effects from our proposed method over the same grid.

With regard to uncertainty quantification, we examine coverage rates with a target confidence level of 0.95 for each method under different scenarios, and the corresponding results, along with the within iteration standard error estimates, are recorded in Table 1. It is quite clear that our proposed method achieves nominal coverage over the other two methods in almost all scenarios. It is worth to notice that although our estimator achieves a better coverage rate than the benchmarks, it suffers from large variation due to the nature of nearest-neighbour matching, especially in the case with small sample size. In general, our method that incorporates both propensity and prognostic scores in identifying subgroups and estimating treatment effects show its advantage in reducing bias, but the resulting large variance is worrisome and it remains space for future study to solve this issue.

We also compare the computational cost of our proposed scored-based method and causal forest, with simulated data from Scenario 1. Endogenous stratification is not considered in this comparison because the method exerts an arbitrary subclassification into a predetermined number of subgroups, and it does not involve a cross-validation process. We record the averaged time consumed to obtain the estimate over 1,000 simulations for both methods. Figure 4 demonstrates that our proposed method (PP) is comparably time-efficient over its competitor, and CF can be very expensive in operational time for large-size problems.

Figure 4

A plot of time per simulation in seconds for our proposed method (PP) and CF against problem size n (50 values from 1,000 up to 10,000). For each method, the time to compute the estimate for one simulated data is averaged over 1,000 Monte Carlo simulations.

These results highlight the promise of our method for accurate estimation of subgroups in heterogeneous treatment effects, all while emphasizing avenues for further work. An immediate challenge is to control the bias in estimation of propensity and prognostic scores. Using more powerful models instead of simple linear regression in estimation is a good way to reduce bias by enabling the trees to focus more closely on the coordinates with the greatest signal. The study of splitting rules for trees designed to estimate causal effects is still in its infancy and improvements may be possible.

5.2 Real data analysis

To illustrate the behaviour of our estimator, we apply our method on the two real-world data sets, one from a clinical study and the other from a complex social survey. Propensity score-based methods are frequently used for confounding adjustment in observational studies, where baseline characteristics can affect the outcome of policy interventions. Therefore, the results from our method are expected to provide meaningful implications for these real data sets. However, one potential limitation of our proposed tree-based approach is that it imposes arbitrary dichotomization with sharp thresholds and thus may obscure potential continuous gradient over the space. For real data analysis in general, our method provides a reasonable simplification that is easy to interpret for the audience and gives some framework for thinking about when such a simplification would be appropriate.

5.2.1 Right heart catheterization (RHC) analysis

While randomized control trials are widely encouraged as the ideal methodology for causal inference in clinical and medical research, the lack of randomized data due to high costs and potential high risks leads to the studies based on observational data. In this section, we are interested in examining the association between the use of RHC during the first 24 h of care in the intensive care unit (ICU) and the short-term survival conditions of the patients. RHC is a procedure for directly measuring how well the heart is pumping blood to the lungs. RHC is often applied to critically ill patients for directing immediate and subsequent treatment. However, RHC imposes a small risk of causing serious complications when administering the procedure. Therefore, the use of RHC is controversial among practitioners, and scientists want to statistically validate the causal effects of RHC treatments. The causal study using observational data can be dated back to Connors et al. [80], where the authors implemented propensity score matching and concluded that RHC treatment lead to lower survival than not performing the treatment. Later, Hirano and Imbens [81] proposed a more efficient propensity-score-based method, and the recent study by Loh and Vansteelandt [82] using a modified propensity score model suggested that RHC significantly affected mortality rate in a short-term period.

A dataset for analysis was first used in the study by Connors et al. [80], and it is suitable for the purpose of applying our method because of its extremely well-balanced distribution of confounders across levels of the treatment [83]. The treatment variable Z in the data indicates whether a patient received a RHC within 24 hours of admission. The binary outcome Y is defined based on whether a patient died at any time up to 180 days since admission. The original data consisted of 5,735 participants with 73 covariates. We preprocess the full data in the way suggested in the Hirano and Imbens [81] and Loh and Vansteelandt [82], by removing all observations that contain null values in covariates, dropping the singular covariate in the reduced data, and encoding categorical variables into dummy variables. The resulted data contains 2,707 observations and 72 covariates, with 1,103 in the treated group ( Z = 1 ) and 1,604 in the control group ( Z = 0 ). Among the 72 observed covariates, there are 21 continuous, 25 binary, and 26 dummy variables transformed from the original six categorical variables.

The result of the prediction model from our proposed method is reported in Figure 5. We observe that the sign of estimated treatment effects varies depending on the value of the propensity score and prognostic score. This particular pattern implies that RHC procedures indeed offer both benefits and risks in affecting patients’ short-term survival conditions. Specifically, we are interested in the occurrence of large positive treatment effects (increase in chance of death) from the estimation. An estimated treatment effect of 0.17 (with a 95% confidence interval between − 0.25 and 1.00) is observed on the group of patients with propensity scores less than 0.64 and prognostic scores less than 0.74, and this group accounts for 56% of the entire sample. For the rest of the sample, we observe relatively small treatment effects in magnitude. Under the scenario of RHC data, a smaller propensity score means that the patient is less likely to receive RHC procedures after admitting to the ICU, and it is related to the availability of RHC procedures at the hospital to which the patient is admitted. A smaller prognostic score tells that the patient has lower underlying chance of death. One possible explanation for this significant positive treatment on this certain group is that drastic change in treatment procedures that were applied to patients who do not actually need the aggressive style of care largely undermine patients’ health conditions after admission and increase the mortality rate. In summary, our findings generally agree with the results and explanations presented in the study by Connors et al. [80], and they offer some insights for practitioners to decide whether they should apply RHC procedures to patients.

Figure 5

Prediction model for RHC data. Left: a tree diagram of predictions on each split based on the two scores, with 95% confidence intervals provided on the bottom. Right: a 2d plot of prediction stratification on the intervals of propensity and prognostic scores, with a colour scale on magnitude of treatment effects.

5.2.2 National medical expenditure survey (NMES)

For the next experiment, we analyze a complex social survey data. In many complex surveys, data are not usually well balanced due to potential biased sampling procedure. To incorporate score-based methods with complex survey data requires an appropriate estimation on propensity and prognostic scores. DuGoff et al. [84] suggested that combining a propensity score method and survey weighting is necessary to achieve unbiased treatment effect estimates that are generalizable to the original survey target population. Austin et al. [85] conducted numerical experiments and showed that greater balance in measured baseline covariates and decreased bias is observed when natural retained weights are used in propensity score matching. Therefore, we include sampling weight as an baseline covariate when estimating propensity and prognostic scores in our analysis.

In this study, we aim to answer the research question: how smoking habit affects medical expenditure over lifetime, and we use the same data set as given in the study by Johnson et al. [86], which is originally extracted from the 1987 NMES. The NMES included detailed information about frequency and duration of smoking with a large nationally representative data base of nearly 30,000 adults, and that 1987 medical costs are verified by multiple interviews and additional data from clinicians and hospitals. A large amount of literature focus on applying various statistical methods to analyze the causal effects of smoking on medical expenditures using the NMES data. In the original study by Johnson et al. [86], the authors first estimated the effects of smoking on certain diseases and then examined how much those diseases increased medical costs. In contrast, Rubin [87], Imai and van Dyk [88], and Zhao et al. [89] proposed to directly estimate the effects of smoking on medical expenditures using propensity-score-based matching and subclassification. Hahn et al. [8] applied Bayesian regression tree models to assess heterogeneous treatment effects.

For our analysis, we explore the effects of extensive exposure to cigarettes on medical expenditures, and we use pack-year as a measurement of cigarette measurement, the same as given in the studies by Imai and van Dyk [88] and Hahn et al. [8]. Pack-year is a clinical quantification of cigarette smoking used to measure a person’s exposure to tobacco, defined by

Pack-year = Number of cigarettes per day 20 × Number of years smoked .

Following that, we determine the treatment indicator Z by the question whether the observation has a heavy lifetime smoking habit, which we define to be greater than 17 pack-years, the equivalent of 17 years of pack-a-day smoking.

The subject-level covariates X in our analysis include age at the times of the survey (between 19 and 94), age when the individual started smoking, gender (male, female), race (white, black, other), marriage status (married, widowed, divorced, separated, never married), education level (college graduate, some college, high school graduate, other), census region (Northeast, Midwest, South, West), poverty status (poor, near poor, low income, middle income, high income), seat belt usage (rarely, sometimes, always/almost always), and sample weight. We select the natural logarithm of annual medical expenditures as the outcome variable Y to maintain the assumption of heteroscedasticity in random errors. We preprocess the raw data set by omitting any observations with missing values in the covariates and excluding those who had zero medical expenditure. The resulting restricted data set contains 7,903 individuals, with 4,014 in the treated group ( Z = 1 ) and 3,889 in the controlled group ( Z = 0 ).

The prediction model obtained from our method, as shown in Figure 6, is simple and easy to interpret. We derive a positive treatment effect across the entire sample, and the effect becomes significant when the predicted potential outcome is relatively low (less than 5.7). These results indicate that more reliance on smoking will deteriorate one’s health condition, especially for those who currently do not have a large amount of medical expenditure. Moreover, we observe a significant positive treatment effect of 1.4 (with a 95% confidence interval between 0.56 and 2.79), and in other words, a certain and substantial increase in medical expenditure for the subgroup with propensity score less than 0.17. It is intuitive to assume that a smaller possibility of engaging in excessive tobacco exposure is associated with healthier living styles. This phenomenon is another evidence that individuals who are more likely to stay healthy may suffer more from excessive exposure to tobacco products. In all, these results support policymakers and social activists who advocate for nationwide smoking ban.

Figure 6

Prediction model for NMES data. Left: a tree diagram of predictions on each split based on the two scores, with 95% confidence intervals provided on the bottom. Right: a 2d plot of prediction stratification on the intervals of propensity and prognostic scores, with a colour scale on magnitude of treatment effects.