1 Introduction

Consider a typical setting of observational study that two conditions (“treatment” and “control”) are not randomly assigned to the units. Deriving a causal conclusion from such observational data is essentially difficult because the treatment exposure may be related to some covariates that are also related to the outcome. In this case, those covariates may be imbalanced between the treatment groups and the naive mean causal effect estimator can be severely biased.

To adjust for the covariate imbalance, the seminal work of Rosenbaum and Rubin [1] points out the essential role of propensity score, the probability of exposure to treatment conditional on observed covariates. This quantity, rarely known in an observation study, may be estimated from the data. Based on the estimated propensity score, many statistical methods are proposed to estimate the mean causal effect. The most popular approaches are matching [2, 3], stratification [4], and weighting [5, 6]. Theoretically, propensity score weighting is the most attractive among these methods. Hirano et al. [7] show that nonparametric propensity score weighting can achieve the semiparametric efficiency bound for the estimation of mean causal effect derived by Hahn [8]. Another desirable property is double robustness. The pioneering work of Robins et al. [5] augments propensity score weighting by an outcome regression model. The resulting estimator has the so-called double robustness property:

Property 1. If either the propensity score model or the outcome regression model is correctly specified, the mean causal effect estimator is statistically consistent.

In practice, the success of any propensity score method hinges on the quality of the estimated propensity score. The weighting methods are usually more sensitive to model misspecification than matching and stratification, and furthermore, this bias can even be amplified by a doubly robust estimator, brought to attention by Kang and Schafer [9]. In order to avoid model misspecification, applied researchers usually increase the complexity of the propensity score model until a sufficiently balanced solution is found. This cyclical process of modeling propensity score and checking covariate balance is criticized as the “propensity score tautology” by Imai et al. [10] and, moreover, has no guarantee of finding a satisfactory solution eventually.

Recently, there is an emerging interest, particularly among applied researchers, in directly incorporating covariate balance in the estimation procedure, so there is no need to check covariate balance repeatedly [11, 12, 13, 14]. In this paper, we study a method of this kind called Entropy Balancing (hereafter EB) proposed in Hainmueller [15]. In a nutshell, EB solves an (convex) entropy maximization problem under the constraint of exact balance of covariate moments. Due to its easy interpretation and fast computation, EB has already gained some popularity in applied fields [16, 17]. However, little do we known about the theoretical properties of EB. The original proposal in Hainmueller [15] did not give a condition such that EB is guaranteed to give a consistent estimate of the mean causal effect.

In this paper, we shall show EB is indeed a very appealing propensity score weighting method. We find EB simultaneously fits a logistic regression model for the propensity score and a linear regression model for the outcome. The linear predictors of these regression models are the covariate moments being balanced. We shall prove EB is doubly robust (Property 1), in the sense that if at least one of the two models are correctly specified, EB is consistent for the Population Average Treatment effect for the Treated (PATT), a common quantity of interest in causal inference and survey sampling. Moreover, EB is sample bounded [18], meaning the PATT estimator is always within the range of the observed outcomes, and it is semiparametrically efficient if both models are correctly specified. Lastly, The two linear models have an exact correspondence to the primal and dual optimization problem used to solve EB, revealing an interesting connection between doubly robust estimation and convex optimization.

Our discoveries can be summarized in the diagram in Figure 1. Conventionally, the recipe given by Robins and his coauthors is to fit separate models for propensity score and outcome regression and then combine them by a doubly robust estimator (see e. g. [5, 9, 19, 20]). In contrast, Entropy Balancing achieves this goal through enforcing covariate balance. The primal optimization problem of EB amounts to an empirical calibration estimator [21, 22], which is widely popular in survey sampling but perhaps not sufficiently recognized in causal inference [23]. The balancing constraints in this optimization problem result in unbiasedness of the PATT estimator under linear outcome regression model. The dual optimization problem of EB is fitting a logistic propensity score model with a loss function different from the negative binomial likelihood. The Fisher-consistency of this loss function (also called proper scoring rule in statistical decision theory, see e. g. Gneiting and Raftery [24]) ensures the other half of double robustness – consistency under correctly specified propensity score model. Since EB essentially just uses a different loss function, other types of propensity score models, for example the generalized additive models [25] can also easily be fitted. A forthcoming article by Zhao [26] offers more discussion and extension to other weighted average treatment effects.

Figure 1:

The role of covariate balance in doubly robust estimation. Dashed arrows: conventional procedure to achieve double robustness. Solid arrows: double robustness of Entropy Balancing via covariate balance.

2 Setting

First, we fix some notations for the causal inference problem considered in this paper. We follow the potential outcome language of Neyman [27] and Rubin [28]. In this causal model, each unit i is associated with a pair of potential outcomes: the response Yi(1) that is realized if Ti=1 (treated), and another response Yi(0) realized if Ti=0 (control). We assume the observational units are independent and identically distributed samples from a population, for which we wish to infer the treatment’s effect. The main obstacle is that only one potential outcome is observed: Yi=TiYi(1)−(1−Ti)Yi(0), which is commonly known as the “fundamental problem of causal inference” [29].

In this paper we focus on estimating the Population Average Treatment effect on the Treated (PATT):

The counterfactual mean μ(0|1)=E[Y(0)|T=1] also naturally occurs in survey sampling with missing data [21, 22] by viewing Y(0) as the only outcome of interest (so T=1 stands for non-response).

Along with the treatment exposure Ti and outcome Yi, each unit i is usually associated with a set of covariates denoted by Xi measured prior to the treatment assignment. In a typical observational study, both treatment assignment and outcome may be related to the covariates, which can cause serious confounding bias. The seminal work by Rosenbaum and Rubin [1] suggest that it is possible to correct the confounding bias under the following two assumptions:

Intuitively, the first assumption says that the observed covariates contain all the information that may cause the selection bias, i. e. there is no unmeasured confounding variable, and the second assumption ensures that the bias-correction information is available across the entire domain of X.

Since the covariates X contain all the information of confounding bias, it is important to understand the relationship between T,Y and X. Under Assumption 1 (strong ignorability), the joint distribution of (X,Y,T) is determined by the marginal distribution of X and two conditional distributions given X. The first conditional distribution e(X)=P(T=1|X) is often called the propensity score and plays a central role in causal inference [1]. The second conditional distribution is the density of Y(0) and Y(1) given X. Since we only consider the mean causal effect in this paper, it suffices to study the mean regression functions g0(X)=E[Y(0)|X] and g1(X)=E[Y(1)|X].

To estimate the PATT defined in 1, a conventional weighting estimator based on the propensity score is the inverse probability weighting (IPW) defined as

Here ∑Ti=t is a short-hand notation of summation over all units i such that Ti=t. This will be repeatedly used throughout this paper. In 2, the control units are weighted proportionally to eˆ(Xi)(1−eˆ(Xi))−1 to resemble the full population. The most popular choice of propensity score model is the logistic regression, where logit(e(x))=log[e(x)/(1−e(x))] is modeled by ∑j=1pθjcj(x) and cj(x) are functions of the covariates.

3 Entropy balancing

Entropy Balancing (EB) is an alternative weighting method proposed by Hainmueller [15] to estimate PATT. EB operates by maximizing the entropy of the weights under some pre-specified balancing constraints:

Hainmueller [15] proposes to use the weighted average ∑Ti=0wiEBYi to estimate the counterfactual mean E[Y(0)|T=1]. This gives the Entropy Balancing estimator of PATT

The functions {cj(⋅)}j=1p in 3 are called moment functions of the covariates. They can be any transformation of X, not necessarily polynomial functions. We use c(X) and cˉ(1) to stand for the vector of cj(X) and cˉj(1), j=1,…,p. We shall see the functions {cj(⋅)}j=1p indeed serve as the linear predictors in the propensity score model and the outcome regression model, although at this point it is not even clear that EB attempts to fit any model.

First, we give some heuristics that allows us to view EB as a propensity score weighting method. Since EB seeks to empirically match the control and treatment covariate distributions, we connect EB with density estimation. Let m(x) be the density function of the covariates X for the control population. The minimum relative entropy principle estimates the density of the treatment population by

where H(m˜∥m)=Em˜[log(m˜(X)/m(X))] is the relative entropy between m˜ and m. As an estimate of the distribution of the treatment group, the optimal m˜ of 5 is the “closest” to the control distribution among all distributions satisfying the moment constraints. Now let w(x)=[P(T=1)⋅m˜(x)]/[P(T=0)⋅m(x)] be the population version of the inverse probability weights in 2. Applying a change of measure, we can rewrite 5 as an optimization problem over w(x):

The EB optimization problem 3 is the finite sample version of 6, where the population distribution m is replaced by the empirical distribution of the control units.

Using the Lagrangian multipliers, one can show the solution of 5 belongs to the family of exponential titled distributions of m [30]:

Here, ψ(θ) is the moment generating function of this exponential family. Consequently, the solution of the population EB 6 is

where α=log(P(T=1)/P(T=0)). This is exactly the logistic regression model with predictors c(x).

Notice that EB is different from the maximum likelihood fit of the logistic regression. The dual optimization problem of 3 is

whereas the maximum likelihood solves

It is apparent from 7 and 8 that EB and maximum likelihood use different loss functions. As a remark, the estimating equations defined by 7 are used to augment the estimating equations defined by 8 in the covariate balancing propensity score (CBPS) approach of Imai and Ratkovic [13]. We will compare the empirical performance of these methods in Section 5.

The optimization problem 7 is strictly convex and the unique solution θˆEB can be efficiently computed by Newton method. The EB weights (solution to the primal problem 3) are given by the Karush-Kuhn-Tucker (KKT) conditions: for any i such that Ti=0,

As a final remark, Entropy Balancing bridges two existing approaches of estimating the mean causal effect:

1. The calibration estimator that is very popular in survey sampling [21, 22, 23];

2. The empirical likelihood approach that significantly advances the theory of doubly robust estimation in observation study [18, 31, 32, 33].

EB is a special case of these two approaches. The main distinction is that it uses the Shannon entropy ∑i=1nwilogwi as the discrepancy function, resulting in an easy-to-solve convex optimization. Due to its easy interpretation, Entropy Balancing has already gained some ground in practice (e. g. [16, 17]).

4 Properties of entropy balancing

We give some theoretical guarantees of Entropy Balancing to justify its usage in real applications. The following is the main theorem of this paper, which shows EB is doubly robust even though its original form 3 does not contain a propensity score model or a outcome regression model.

We give two proofs of the first claim in Theorem 1. The first proof reveals an interesting connection between the primal-dual optimization problems 3 and 7 and the statistical property, double robustness, which motivates the interpretation in Figure 1. The second proof uses a stabilization trick in Robins et al. [8].

First proof (sketch) The consistency under the linear model of logit(P(T=1|X)) is a consequence of the dual optimization problem 7. See Section 3 for a heuristic justification via the minimum relative entropy principle and Appendix A for a rigorous proof by using the M-estimation theory.

The consistency under the linear model of Y(0) can be proved by expanding E[Y(0)|X] and ∑Ti=0wiYi. Here we provide an indirect proof by showing that augmenting EB with a linear outcome regression does not change the estimator. Given an estimated propensity score model eˆ(x), the corresponding weights eˆ(x)/(1−eˆ(x)) for the control units, and an estimated outcome regression model gˆ0(x), a doubly robust estimator of PATT is given by

This estimator satisfies Property 1, i. e. if eˆ(x)→e(x) or gˆ0(x)→g(x), then γˆDR is statistically consistent for γ. To see this, in the case that gˆ0(x)→g0(x), the first sum in 10 is consistent for γ and the second sum in 10 has mean going to 0 as n→∞. In the case where gˆ0(x)g0(x) but eˆ(x)→e(x), the second sum in 10 is consistent for the bias of the first sum (as an estimator of γ).

When the estimated propensity score model eˆ(x) is obtained by the EB dual problem 7 and the estimated outcome regression model is gˆ0(x)=∑j=1pβˆjcj(x), we have

Therefore by enforcing covariate balancing constraints, EB implicitly fits a linear outcome regression model and is consistent for γ under this model.□

Second proof. This proof is pointed out by an anonymous reviewer. In a discussion of Kang and Schafer [9], Robins et al. [34] show that one can stabilize the standard doubly robust estimator in a number of ways. Specifically, one trick suggested by Robins et al. [34] is to estimate the propensity score, say e˜(x), by the following estimating equation

Then one can estimate PATT by the IPW estimator 2 by replacing eˆ(Xi) with e˜(Xi). This estimator is sample bounded (the estimator is always within the range of observed values of Y) and doubly robust with respect to the parametric specifications of e˜(x)=e˜(x;θ) and gˆ0(x)=gˆ0(x;β). The only problem with 11 is it may not have a unique solution. However, when logit(e(x)) and g0(x) are assumed linear in c(x), 11 corresponds to the first order condition of the EB dual problem 7. Since 7 is strictly convex, it has an unique solution and e˜(X;θ) is the same as the EB estimate eˆ(X;θ). As a consequence, γˆEB is also doubly robust.

To prove the second claim in Theorem 1, we compute the asymptotic variance of γˆEB using the M-estimation theory. To state our results, we need to introduce four differently weighted covariance-like functions for two random vectors a1 and a2 of length p:

It is obvious that H≥K and usually G≥H. To make the notation more concise, c(X) will be abbreviated as c and Y(0) as 0 in subscripts. For example, Hc,0=Hc(X),Y(0), Gc,1=Gc(X),Y(1) and Kc=Kc(X),c(X).

The proof of Theorem 2 is given in Appendix A. The H, G and K matrices in Theorem 2 can be estimated from the observed data, yielding approximate sampling variances for γˆEB and γˆIPW. Alternatively, variance estimates may be obtained via the empirical sandwich method (e. g. [35]). In practice (particularly in simulations where we compare to a known truth), we find that the empirical sandwich method is more stable than the plug-in method, which is consistent with the suggestion in Lunceford and Davidian [20] for PATE estimators.

To complete the proof of the second claim in Theorem 1, we compare these variances with the semiparametric variance bound of γ with unknown e(X) derived by Hahn ([8], Theorem 1):

After some algebra, one can express V∗ in terms of H⋅,⋅ and G⋅,⋅ defined above:

Now assume logit(P(T=1|X))=θTc(X) and E[Y(t)|X]=β(t)Tc(X),t=0,1, it is easy to verify that

Similarly, Gc,t=Gcβ(t),t=0,1. From here it is easy to check VEB and V∗ are the same. Since Entropy Balancing reaches the efficiency bound in this case, obviously VEB<VIPW when both models are linear.

If logit(P(T=1|X))=θTc(X) is true but E[Y(t)|X]=β(t)Tc(X) is not true for some t=0,1, there is no guarantee that EB has the smaller asymptotic variance. In practice, the features c(X) in the outcome regression models are almost always correlated with Y. This correlation compensates the slight efficiency loss of not maximizing the likelihood function in logistic regression. As a consequence, the variance VEB in 12 is usually smaller than VIPW in 13. This efficiency advantage of EB over IPW is verified in the next section using simulations.

5 Simulations

5.1 Kang-Schafer example

We use the simulation example in Kang and Schafer [9] to compare EB weighting with IPW (after maximum likelihood logistic regression) and the over-identified Covariate Balancing Propensity Score (CBPS) proposed by Imai and Ratkovic [13]. The simulated data consist of {Xi,Zi,Ti,Yi},i=1,…,n}. Xi and Ti are always observed, Yi is observed only if Ti=1, and Zi is never observed. To generate this data set, Xi is distributed as N(0,I4), Zi is computed by first applying the following transformation:

Zi1=exp(Xi1/2),Zi2=Xi2/(1+exp(Xi1))+10,Zi3=(Xi1Xi3+0.6)3,Zi4=(Xi2+Xi4+20)2.

Next we normalize each column such that Zi has mean 0 and standard deviation 1.

In one setting, Yi is generated by Yi=210+27.4Xi1+13.7Xi2+13.7Xi3+13.7Xi4+εi, εi∼N(0,1) and the true propensity scores are ei=expit(−Xi1+0.5Xi2−0.25Xi3−0.1Xi4). In this case, both Y and T can be correctly modeled by (generalized) linear model of the observed covariates X.

In the other settings, at least one of the propensity score model and the outcome regression model is incorrect. In order to achieve this, the data generating process described above is altered such that Y or T (or both) is linear in the unobserved Z instead of the observed X, though the parameters are kept the same.

For each setting (4 in total), we generated 1000 simulated data sets of size n = 200 and 1000, then apply various methods discussed earlier including

1. IPW, CBPS: the IPW estimator in 2 with propensity score estimated by logistic regression or CBPS (since the estimand is overall mean, we use the CBPS weights tailored for estimating PATE);

2. EB: the Entropy Balancing estimator (the EB weights are used to estimate the unobserved mean E[Y|T=0]);

3. IPW+DR, CBPS+DR: the doubly robust estimator in 10 with propensity score estimated by logistic regression or CBPS.

The simulation results are presented in Figure 2 and Figure 3. Figure 2 and Figure 3 show the covariate imbalance before adjustment in terms of standardized difference. Figure 2 and Figure 3 show the mean estimates given by the five different methods. First, notice that the doubly robust estimator “IPW+DR” performs poorly when both models are misspecified (bottom-right panel in Figure 2) and Figure 3). In fact, all the three doubly robust methods are worse than just using IPW. Second, the three doubly robust estimators have exactly the same variance if the Y model is correct (top two panels in Figure 2 and Figure 4). It seems that how one fits the propensity score model has no impact on the final estimate. This is related to the observation in Kang and Schafer [9] that, in this example, the plain OLS estimate of Y actually outperforms any method involving the propensity score model. Discussion articles such as Robins et al. [34] and Ridgeway and McCaffrey [36] find this phenomenon very uncommon in practice and is most likely due to the estimated inverse probability weights are highly variable, which is a bad setting for doubly robust estimators.

Figure 2:

Kang-Schafer example: sample size n = 200. Both propensity score model and outcome regression model can be correct or incorrect, so there are four scenarios in total. We generate 1,000 simulations in each scenario. (a) Covariate imbalance before adjustment. (b) Mean estimates. The methods are: Inverse Propensity Weighting (IPW), Covariate Balancing Propensity Score (CBPS), Entropy Balancing (EB), and doubly robust versions of the first two (IPW+DR, CBPS+DR). Target mean is 210 and is marked as a black horizontal line to compare the biases. Numbers printed at Y = 230 are the sample standard deviations to compare efficiency.

Figure 3:

Kang-Schafer example: sample size n=200. Both propensity score model and outcome regression model can be correct or incorrect, so there are four scenarios in total. We generate 1,000 simulations in each scenario. (a) Covariate imbalance before adjustment. (b) Mean estimates. The methods are: Inverse Propensity Weighting (IPW), Covariate Balancing Propensity Score (CBPS), Entropy Balancing (EB), and doubly robust versions of the first two (IPW+DR, CBPS+DR). Target mean is 210 and is marked as a black horizontal line to compare the biases. Numbers printed at Y = 230 are the sample standard deviations to compare efficiency.

Regarding Entropy Balancing (EB), we find that:

1. If both T and Y models are misspecified, EB has smaller bias than the conventional “IPW+DR” or “CBPS+DR”. So EB seems to be less affected by such unfavorable setting.

2. When T model is correct but Y model is wrong (bottom-left panel in Figure 2 and Figure 3), EB has the smallest variance among all estimators. This supports the conclusion of our efficiency comparison of IPW and EB in Section 4.

Finally notice that the same simulation setting is used in Tan [18] to study the performance of a number of doubly robust estimators. The reader can compare the Figure 2 and Figure 3 with the results there. The performance of Entropy Balancing is comparable to the best estimator in Tan [18].

5.2 Lunceford-Davidian example

We provide another simulation example by Lunceford and Davidian [20] to verify claims in Theorems 1 and 2. In this simulation, the data still consist of {(Xi,Zi,Ti,Yi),i=1,…,n}, but all of them are observed. Both Xi and Zi are three dimensional vectors. The propensity score is only related to X through:

logit(P(Ti=1))=β0+∑j=1βjXij.

Note the above does not involve elements of Zi. The response Y is generated according to

Yi=ν0+∑j=13νjXij+ν4Ti+∑j=13ξjZij+εi;εi∼N(0,1).

The parameters here are set to be ν=(0,−1,1,−1,2)T, and β is set as:

βno=(0,0,0,0)T,βmoderate=(0,0.3,−0.3,0.3)T,orβmoderate=(0,0.3,−0.3,0.3)T,or

The choice of β depends on the level of association of T and X. ξ is based on a similar choice on the level of association of Y and Z:

ξno=(0,0,0)T,ξmoderate=(−0.5,0.5,0.5)T,orξstrong=(−1,1,1)T.

The joint distribution of (Xi,Zi) is specified by taking Xi3∼Bernoulli(0.2) and then generate Zi3 as Bernoulli with

PZi3=1|Xi3=0.75Xi3+0.25(1−Xi3).

Conditional on Xi3, (Xi1,Zi1,Xi2,Zi2) is then generated as multivariate normal N(aXi3,BXi3), where a1=(1,1,−1,−1)T, a0=(−1,−1,1,1)T and

B0=B1=10.5−0.5−0.50.51−0.5−0.5−0.5−0.510.5−0.5−0.50.51.

Figure 4 shows the covariate imbalance in the three settings before adjustment.

Figure 4:

Lunceford-Davidian example: covariate imbalance before adjustment.

The data generating model implies that the true PATT is γ=2. Since the outcome Y depends on both X and Z, we always fit a full linear model of Y using X and Z, if such model is needed. T only depends on X, so it is not necessary to include Z in propensity score modeling. However, as pointed out by Lunceford and Davidian [20], it is actually beneficial to “overmodel” the propensity score by including Z in the model. Here we will try both possibilities, the “full” modeling of T using both X and Z, and the “partial” modeling of T using only X. Since the estimand is PATT in this case, we use the over-identified CBPS weights tailored for estimating PATT.

We generated 1000 simulated data sets and the results are shown in Figure 4 for “full” propensity score modeling and Figure 6 for “partial” propensity score modeling. We make the following comments about these two plots:

1. IPW and all other estimators are always consistent, no matter what level of association is specified. This is because the propensity score model is always correctly specified.

2. When using the “full” propensity score modeling, all doubly robust estimators (EB, IPW+DR, CBPS+DR and EB+DR) have almost the same sample variance. This is because all of them are asymptotically efficient.

3. CBPS, to our surprise, does not perform very well in this simulation. It has smaller variance than IPW but this comes with the price of some bias. If we use the partial propensity score model (only involve X, Figure 6), this bias is smaller but still not negligible. While it is not clear what causes this bias, one possible reason is that the optimization problem of CBPS is nonconvex, so the local solution which is used to construct γ estimator could be far from the global solution. Another possibility is that CBPS uses GMM or Empirical Likelihood to combine likelihood with imbalance penalty, which is less efficient than maximum likelihood directly. Thus, although the estimator is asymptotically unbiased, the convergence spend to the true γ is quite slower than IPW. CBPS combined with outcome regression (CB+DR) fixes the bias and inefficiency issue occurred in CBPS without outcome regression.

4. EB, in contrast, performs quite well in this simulation. It has relatively small variance, particularly if we use the “full” model in which both X and Z are balanced.

5. The difference between EB and EB+DR is that while EB only balances “partial” or “full” covariates, EB+DR additionally combines a outcome linear regression model on all the covariates. As shown in the first proof of Theorem 1, when the “full” covariates are used, EB is exactly the same as EB+DR. We can observe this from Figure 4. When EB only balances “partial” covariates, the two methods are different and indeed EB+DR is more efficient in Figure 6 since it fits the correct Y model.

6. Using the “full” propensity score model improves the efficiency of pure weighting estimators (IPW, CBPS and EB) a lot, but has very little impact on estimators that involves an outcome regression model (IPW+DR and CBPS+DR) compared to “partial” propensity score modeling. Although EB could be viewed as fitting a outcome model implicitly, the “partial” EB estimator only uses X in the outcome model, that is precisely the reason why it is not efficient. Thus there are both robustness and efficiency reasons that one should include all relevant covariates in EB, even if the covariates affect only one of T and Y.

Figure 5:

Results of the Lunceford-Davidian example (full propensity score modeling). The propensity score model and outcome regression model, if applies, are always correctly specified, but the level of association between T or Y with X or Z could be different, ended up with 9 different scenarios. X are confounding covariates and Z only affects the outcome. We generate 1,000 simulations of 1,000 in each scenario and apply five different estimators. The true PATT is 2 and is marked as a black horizontal line to compare the biases of the methods. Numbers printed at Y = 5 are the sample standard deviation of each method, in order to compare their efficiency.

Figure 6:

Results of the Lunceford-Davidian example (partial propensity score modeling). The settings are exactly the same as Figure 4 except the methods here don’t use Z in their propensity score models.

In summary, EB outperforms IPW in all the simulations, making it an appealing alternative to the conventional propensity score weighting methods.