Entropy Balancing for Continuous Treatments

2.1 Identifying Assumptions

First, conditional on observed pre-treatment covariates X, potential outcomes must be independent of the treatment intensity received, i.e.

This assumption is called the conditional independence assumption (CIA, Lechner 2001) also known as the selection-on-observables assumption (Heckman and Robb 1985) and requires that the researcher observes all covariates X that simultaneously determine the selection into different treatment intensities as well as the outcome of interest. The CIA is potentially a very strong assumption. Hence, it needs to be discussed on a case-by-case basis for the application at hand and requires substantive knowledge about the selection mechanism at play. The second assumption requires there to be common support, i.e. the conditional density of treatment needs to be positive over T :

In comparison to the binary treatment case, this condition is potentially much stronger and may be violated, especially in regions of low density in terms of the treatment variable. If this is the case, the sample needs to be trimmed and the DRF is estimated on the subset of observations in order to avoid extrapolation (Crump et al. 2009; Lechner and Strittmatter 2019). Lastly, one needs to assume the so-called stable-unit treatment value assumption (SUTVA, see Rubin 1980), requiring that each individual’s outcome only depends on their own level of treatment intensity. Essentially, this rules out general equilibrium and spill-over effects of treatment (see Imbens and Wooldridge 2009; Manski 2013, for examples).

2.2 Re-weighting Methods Based on the Generalized Propensity Score

Comparing outcomes of individuals with exactly the same set of X but different T to estimate the DRF quickly becomes infeasible with growing dimension of X. To avoid this curse of dimensionality, Hirano and Imbens (2004) show that the conditional independence assumption also holds by conditioning on the generalized propensity score (GPS) R _i = f _T∣X (T _i ∣X _i ), i.e. the conditional density of the treatment intensity evaluated at T _i and X _i . This is because, conditional on the GPS, the covariate distribution is expected to be independent of the treatment intensity received, i.e.

Weighting methods – which will be the main focus of this paper – make use of this fact and re-weight observations based on their GPS. Some of the most prominent methods of this kind will be used for comparison regarding the proposed methodology. The first approach originates from inverse probability weighting (IPW, see Horvitz and Thompson 1952) and is provided by Robins, Hernán, and Brumback (2000) who generalize IPW with weights defined as

Similar to Hirano and Imbens (2004), Robins, Hernán, and Brumback (2000) estimate (un-) conditional densities f _T (T _i ) and f _T∣X (T _i ∣X _i ) based on OLS regressions and the assumption that the treatment intensity follows a normal distribution. However, such distributional assumptions are likely to be violated in practice, which may result in serious bias.^[3] Other procedures, such as the machine-learning based approach by Su, Ura, and Zhang (2019), are more robust but also less intuitive and user-friendly.

The other comparison methods combine the estimation of the GPS with algorithmic optimization in order to minimize the Pearson correlation between covariates and the treatment variable. The methods used are the parametric and the non-parametric covariate balancing generalized propensity score (CBGPS, Fong, Hazlett, and Imai 2018) as well as generalized boosted models (GBM, Zhu, Coffman, and Ghosh 2015). As will become clear in the remainder of the paper, these automated balancing approaches based on the GPS mostly improve upon covariate balance and robustness relative to standard IPW. However, these re-weighting procedures may not achieve satisfactory balance, leaving estimates susceptible to bias due to residual imbalance.

2.3 Estimation of the Dose–Response Function

Once balancing weights are obtained, the DRF can be estimated using parametric or non-parametric regression techniques based on the re-weighted sample. In general, for this to yield consistent estimates, the aforementioned identifying assumptions need to be true and balancing weights have to render the covariate distribution independent of the treatment intensity. If any of these conditions is not fulfilled, resulting estimates will be inconsistent. Moreover, there is a risk of mis-specification bias due to wrong functional form assumptions for the DRF. To avoid this, flexible parametric specifications using (fractional-) polynomials (Sauerbrei and Royston 1999) or (Basis) splines (d. Boor 1978) in the treatment variable or completely non-parametric methods such as local linear regressions (Fan 1992) may be preferred. To provide some guidance to the reader, global polynomial regressions are compared to local linear regressions in the simulations in Section 4.

3.1 Entropy Balancing for Continuous Treatments

To extend the entropy balancing approach to the case of continuous treatments, a few changes have to made. First, in contrast to the binary case, the estimation of balancing weights w _i is performed for the set of treated units, i.e. units i with T _i > 0. This is because we wish to impute average outcomes of treated units had they received certain doses of the treatment.^[4] Second, the choice of moment conditions need to be altered. Following recent developments in the automated balancing literature, Entropy Balancing for Continuous Treatments (EBCT) focuses on a set of uncorrelatedness conditions based on Pearson correlations.

Analogous to above, define X ̃ i to be the de-meaned version of the (possibly expanded) covariate vector. Furthermore, let T ̃ i p be the de-meaned pth order term of the treatment intensity. Lastly, define the column vector g ( p , X ̃ i , T ̃ i ) = X ̃ i ′ , T ̃ i , … , T ̃ i p , X ̃ i ′ ⋅ T ̃ i , … , X ̃ i ′ ⋅ T ̃ i p ′ . For a given p, the EBCT method aims to solve the following constrained minimization problem:

EBCT minimizes the loss function H(w) subject to the balancing constraints in terms of g ( p , X ̃ i , T ̃ i ) and the normalizing constraints that weights have to sum up to one and be strictly positive. If p = 1 is chosen, resulting weights retain unconditional means of covariates (and their higher-order or cross-moments if included) as well as the treatment intensity and they purge the treatment variable from its correlation with covariates. However, simple uncorrelatedness may not be sufficient to achieve acceptable covariate balance and hence, obtain consistent estimates of the DRF in the next step (Yiu and Su 2018). This is because even if X ̃ i is chosen to be relatively flexible, uncorrelatedness between T ̃ i and X ̃ i does not in general imply independence, i.e. the covariate distribution is not necessarily the same across the entire treatment intensity distribution. This is an issue also touched upon by Vegetabile et al. (2021), who also investigate the performance of EBCT. However, in contrast to this paper, Vegetabile et al. (2021) do not further examine the influence on rendering higher moments of T _i uncorrelated with X ̃ i on the performance of EBCT.^[5]

In order to detect potential residual confounding, it is crucial to flexibly check for remaining dependencies by estimating pseudo-DRFs using covariates (and their higher moments or interactions) as pseudo-outcomes for an additional balancing check. If covariates are truly balanced, these pseudo-DRFs should be completely flat and their derivatives should be zero. If this is not the case, p should to be increased and balance re-assessed. For example, choosing p = 2 means that, in addition to the properties above, the marginal mean of T i 2 is retained and X ̃ i is rendered uncorrelated with T i 2 . It should be clear there is most likely a bias-variance trade-off here: by increasing p, one approximates independence between T _i and X _i more closely and thus, one is likely to reduce bias but this may come at the cost of increased variance of estimates.

3.2 Implementing EBCT

To implement the proposed EBCT approach, this paper also uses the Kullback (1959) entropy metric h(w _i ) = w _i  ln(w _i /q _i ). If no base weights q _i are specified, uniform weights q _i = 1/N ₁ ∀ i are used. This implies that EBCT chooses balancing weights such that they differ as little as possible from baseline weights in terms of the entropy metric while achieving zero correlation conditions in the re-weighted sample. Notice that the loss function attains a minimum at w _i = q _i ∀ i and is undefined for non-positive weights. The latter property allows to drop the positivity constraint on weights, reducing the optimization problem to one with only equality constrains. Using the Lagrange method, the constrained optimization can be re-written as an unconstrained optimization as

where λ and γ are Lagrange-multipliers on the constraints. As ∂ 2 H / ∂ w i 2 > 0 for all w _i > 0 and because the constraints are linear in w _i , the optimization problem (7) has a global minimum if the constraints are consistent (Boyd and Vandenberghe 2004, Chapter 5). In order to reduce the dimensionality of the optimization problem, the implied structure of balancing weights is obtained by re-arranging the first-order condition ∂ L / ∂ w i = 0 and plugging the result into the condition ∂ L / ∂ λ = 0 . This yields the weighting function in terms of the Lagrange-multipliers γ, base weights q _i and the data g ( p , X ̃ i , T ̃ i ) as

where λ has been cancelled out. Hence, when p = 1, weights implied by EBCT are a log-linear function of a linear index containing covariates (and their transformations), the treatment intensity and their cross-products. When p > 1, the linear index also contains the higher-order terms of T ̃ i and their interactions with X ̃ i . Substituting this expression into the Lagrange function yields the dual L d as

Differentiating L d with respect to γ yields the first-order conditions

where γ* refer to the multiplier values at the optimum.^[6] As Eq. (10) are non-linear in those multipliers, they have to be solved for numerically. This is done using a quasi-Newton optimization approach. Due to the convexity of the optimization problem, the algorithm tends to converge relatively quickly even in large datasets. Once values for γ* are obtained, balancing weights are backed out using (8) for subsequent analysis. As noted by Hainmueller (2012), optimization can be performed iteratively to limit the influence of units with potentially extreme weights. To do so, the researcher estimates EBCT weights and truncates excessive weights beyond some threshold. For the binary literature, Imbens (2004) suggests a value of 5%. For the estimation of causal effects in the continuous treatment case, the threshold should probably be substantially smaller, e.g. 1–2%, in order to avoid an overly large influence of single observations on the estimated shape of the DRF. Having truncated weights to some threshold, they are used as base weights for a second run of the algorithm. Resulting weights should display smaller maximum weights.

4.1 The Data

The 1987 NMES data includes information on respondents smoking behavior, their medical expenditures and some background characteristics. For the analysis, N = 9, 408 current or previous smokers are used.^[8] For them, the number of pack years smoked (=smoking duration in years ⋅ cigarette packs smoked per day) is generated as a measure of cumulated exposure to smoking. Available background characteristics are the continuous (starting) age, a male indicator and categorical variables on race, seatbelt usage, education, marital status, census region of residence and poverty status. The treatment variable is taken in logarithmic terms to reduce the skewness of the distribution and make the normality assumption by IPW and CBGPS more plausible. A histogram of the resulting distribution can be found in Figure 1. In general, the data are characterized by relatively strong correlations between the treatment variable and covariates, especially the (starting) age. As second and third-order terms of (starting) age are strongly related to the treatment intensity, they are also included in the specification for the estimation of balancing weights.

Figure 1:

Histogram of the treatment intensity.

This graph shows a histogram of the treatment intensity t = ln(pack years).

4.2 Estimating Weights and Inspecting Correlations

First, balancing weights are estimated. Computation times differ substantially from below 1 s for EBCT with p = 1 up to 17 min for npCBGPS.^[9] Having estimated weights, the next step is to check whether the (automated) balancing algorithms did what they were designed to do: purge the treatment variable from its correlation with covariates. Table 1 gives an overview of (mean absolute) Pearson correlations before and after weighting for all control variables and their transformations used.

Before weighting there is a substantial absolute correlation between ln(pack years) and the polynomials in age (38–47%). Correlations with the other covariates are smaller in magnitude but often substantive nonetheless. Compared to the unweighted sample, IPW increases the mean absolute correlation between covariates and the smoking intensity from 12 to 16%, leading to higher correlations in magnitude for several variables. Clearly, this is a very unfavorable balancing outcome. The parametric CBGPS reduces correlations for almost all variables or leads to slight increases in correlations for variables that were initially almost completely uncorrelated with the treatment. Its non-parametric counterpart is about as successful in reducing absolute correlations for covariates that were initially heavily related to the treatment. However, it leads to larger increases in imbalance for variables with lower initial correlations. For example, the correlation between the male indicator and the treatment increases from 14 to 20% after weighting using the npCBGPS. GBM is almost as effective in alleviating correlations between covariates and the treatment as the parametric CBGPS. Only the correlation between the male indicator and the treatment of 14% remains above the 0.1 threshold suggested by Zhu, Coffman, and Ghosh (2015). EBCT eradicates all correlations of the smoking intensity with covariates in this application, independent of which polynomial order p is used to obtain the balancing weights. Comparing the smallest maximum weight shares, one can see that EBCT with p = 1 yields the lowest value with 0.4% compared to 0.7% (GBM), 0.8% (EBCT, p = 2), 1.4 (EBCT, p = 3), 6.6% (npCBGPS), IPW (7.1%) and almost 13% (CBGPS). Hence, especially the parametric CBGPS method achieves better balance only by allowing for relatively extreme weights in this setting.

4.3 Balance Checks via Pseudo-DRFs

As mentioned earlier, uncorrelatedness does not necessarily imply the equality of covariate distributions across different levels of the treatment variable. Hence, it is crucial to check for potentially remaining imbalances before estimating the DRF to avoid bias in the resulting estimates. To inspect balance further, the following analysis estimates pseudo-DRFs using linear regressions of each covariate on a fourth-order polynomial in the treatment variable. This is similar in spirit to the procedure suggested by Smith and Todd (2005) in the binary treatment setting. Figure 2 shows box-plots of the regression R ² for all covariates used in the estimation of balancing weights as well as all sensible two-way interactions. If covariates are truly balanced, the pseudo-DRFs should be completely flat and thus, R ² should be close to zero after weighting.

Figure 2:

Balancing quality – regression R ².

This graph shows boxplots of R ² from a regression of each covariate and all sensible two-way interactions on a fourth-order polynomial of the treatment intensity t = ln(pack years) before and after weighting. Re-weighting approaches employed are inverse probability weighting estimated via OLS (IPW, see Robins, Hernán, and Brumback 2000), (non-) parametric covariate balancing generalized propensity scores (np-/CBGPS, see Fong, Hazlett, and Imai 2018), generalized boosted modeling (GBM, see Zhu, Coffman, and Ghosh 2015) as well as the novel entropy balancing for continuous treatments (EBCT). The parameter p shows up to which power the treatment variable has been rendered uncorrelated with covariates.

As Figure 2 reveals, this is the case for the majority of variables. However, all methods that focus on simple uncorrelatedness as their balancing target display unwanted residual imbalances due to lingering non-linear relationships between covariates and the treatment intensity.^[10] This is despite the relatively flexible specification used for the estimation of balancing weights and shows the importance of checking for covariate balance beyond simple correlations before analyzing outcomes. Increasing the polynomial degree p for EBCT drastically improves balance, for p = 3 only negligible regression R ² are found after balancing.

4.4 Bias and Mean Squared Error

This paragraph provides evidence on the expected performance of effect estimates in terms of (absolute) bias and root mean squared error using the different balancing approaches considered. Similar to the approach by Frölich, Huber, and Wiesenfarth (2017), empirical Monte-Carlo simulations are performed based on real data.^[11] For this analysis, the observed medical expenditures from the N = 9, 408 current or previous smokers are replaced by matched outcomes from the group of the N = 9, 804 non-smokers also available in the dataset. More specifically, the smokers and non-smokers are matched using a non-parametric nearest-neighbor matching with replacement based on the Mahalanobis distance metric. The list of variables to match on includes the same set of variables as described above, except for starting age. The great advantage of this approach is that due to its non-parametric nature, the resulting dataset features realistic (and possibly complex) functional relationships between covariates and the matched outcome. As the matched medical expenditures stems from non-smokers, they are unaffected by the smoking intensity such that the CIA holds and estimates should display completely flat DRFs. This allows to inspect the performance of estimators in a realistic setting. For the simulations, R = 200 samples of size N = 9, 408 are drawn with replacement from the original data. For each of the samples, balancing weights are re-estimated. Outcome regressions are implemented via weighted least squares regression using a cubic polynomial in ln(pack years) as well as through weighted local linear kernel regressions based on the Epanechnikov Kernel (Fan 1992). For the kernel regression, the bandwidth is chosen in a data-driven manner via cross-validation as implemented by the locpol package (Cabrera 2018). From these estimates, the dose–response functions E[Y _i (t)] are obtained as predicted values for an evenly-spaced grid of values of the treatment variable between its first and 99th percentile. Table 2 shows estimates of absolute bias and root mean squared error (RMSE), averaged over all grid values.

The simulation results show that IPW and the parametric CBGPS barely lead to any gains in terms of bias and yield an increase in RMSE relative to the unweighted scenario. When using the polynomial regression, they even increase bias and RMSE compared to the unweighted regression. The non-parametric CBGPS on the other hand reduces bias substantially but leads to an increase in RMSE, implying relatively large variance of estimates. GBM and EBCT with p = 1 perform similarly but still display non-negligible bias. Only when EBCT also renders higher moments of the treatment variable uncorrelated with covariates, bias and also RMSE drop substantially. Moreover, the results show that bias and RMSE is reduced by switching from parametric to non-parametric estimation of the DRF in most cases. When using EBCT with p > 1, DRF estimates using local linear regressions display slightly larger bias and RMSE than using a global polynomial regression. Thus, it appears that non-parametric estimation of the DRF is able to mitigate the effects of residual confounding on resulting estimates at least to some degree.

To inspect the performance of the estimators further, Figure 3 displays estimates of the bias as a function of the treatment variable. To indicate variability of estimates, 95% normal confidence bands are shown based on the Monte-Carlo error (Koehler, Brown, and Haneuse 2009). For the sake of brevity, only the results for non-parametric DRF estimates are shown.

Figure 3:

Empirical Monte-Carlo simulation: bias in detail.

The graph shows further results from the empirical Monte-Carlo simulation based on 200 replications. Displayed are the estimated bias as well as 95% normal confidence intervals based on the Monte-Carlo error a function of the treatment variable t =ln(pack years). For the simulation, observed outcomes of smokers are replaced by matched outcomes from non-smokers. Estimates shown are obtained via weighted local linear kernel regression based on the Epanechnikov kernel with bandwidths chosen via cross-validation. Re-weighting approaches employed are inverse probability weighting estimated via OLS (IPW, see Robins, Hernán, and Brumback 2000), (non-) parametric covariate balancing generalized propensity scores (np-/CBGPS, see Fong, Hazlett, and Imai 2018), generalized boosted modeling (GBM, see Zhu, Coffman, and Ghosh 2015) as well as the novel entropy balancing for continuous treatments (EBCT). The parameter p shows up to which power the treatment variable has been rendered uncorrelated with covariates.

Figure 3 shows that in the unweighted dataset, naive DRF estimates are significantly biased downward in the lower tail and biased upward in the upper tail of the treatment intensity distribution. While IPW and (np-)CBGPS remove some of that bias, they display excessive variance in at least some parts of the distribution of the treatment variable. GBM and EBCT with p = 1 have much lower variance, but are still display a significant upward bias in the upper tail. Only when higher moments of T _i are also rendered uncorrelated using EBCT, bias turns insignificant. Comparing the width of confidence intervals for EBCT with p = 2 and p = 3, one can see an increase in the variance mainly in the tails when increasing the degree of polynomial in T _i to be balanced.

4.5 Estimating the Effects of Smoking

This part estimates the actual DRF of the smoking intensity regarding observed medical expenditures. Due to the relatively poor performance of IPW and (np-)CBGPS in the simulations, estimates based on those weighting procedures are not presented in the main text. All results can be found in Figures A.1 and A.2 in the Appendix. To obtain standard errors, the bootstrap (Efron and Tibshirani 1986; MacKinnon 2006) is used as Vegetabile et al. (2021) find that it performs well in combination with EBCT. Standard errors are estimated using 200 bootstrap replications, including the estimation of balancing weights (and the kernel bandwidth) in each replication.^[12] Resulting estimates of the DRFs along with 95% confidence bands based on the normal approximation are displayed in Figure 4.

Figure 4:

Effects of smoking on medical expenditures.

This graph shows the estimated dose–response function (DRF) between the log(pack years) and medical expenditures based on an either weighted linear regressions using cubic specification or weighted local linear kernel regressions based on the Epanechnikov kernel with bandwidths chosen via cross-validation. Results are only shown for generalized boosted modeling (GBM, see Zhu, Coffman, and Ghosh 2015) as well as the novel entropy balancing for continuous treatments (EBCT). The full set of results can be found in Figures A.1 and A.2. The shaded grey areas represent 95% normal confidence intervals of the DRF based on bootstrapped standard errors obtained using R = 200 replications.

Overall, differences in point estimates of the DRF between global polynomial and local linear regressions are small, although confidence bands are somewhat larger for the non-parametric estimates. All estimates show an increase in medical expenditures with large smoking intensities. While GBM and EBCT with p = 1 show a slight decline in medical expenditures for low treatment intensities, estimates based on EBCT with p > 1 no longer display such a pattern. Hence, also rendering higher orders of the treatment intensity uncorrelated with covariates leads to the most reasonable estimates of the DRF as smoking can be expected to have non-negative effects on medical expenditures.