A Framework for Covariate Balance using Bregman Distances

TL;DR

This paper introduces a new framework for covariate balancing in observational studies using Bregman distances, providing a flexible and theoretically grounded approach that aligns with existing methods.

Contribution

It proposes a convex optimization framework for covariate balancing weights based on Bregman distances, unifying and generalizing existing covariate balancing techniques.

Findings

The framework produces weights identical to existing covariate balancing methods.

Numerical studies demonstrate the equivalence and effectiveness of the proposed approach.

Abstract

A common goal in observational research is to estimate marginal causal effects in the presence of confounding variables. One solution to this problem is to use the covariate distribution to weight the outcomes such that the data appear randomized. The propensity score is a natural quantity that arises in this setting. Propensity score weights have desirable asymptotic properties, but they often fail to adequately balance covariate data in finite samples. Empirical covariate balancing methods pose as an appealing alternative by exactly balancing the sample moments of the covariate distribution. With this objective in mind, we propose a framework for estimating balancing weights by solving a constrained convex program where the criterion function to be optimized is a Bregman distance. We then show that the different distances in this class render identical weights to those of other covariate balancing methods. A series of numerical studies are presented to demonstrate these similarities.

Tables

Table 1. Table 1: Average estimate, Monte Carlo standard error, residual mean squared error, and empirical bias of the constant conditional ATE using the four methods for estimating balancing weights described in Section 5.1 . IPW uses inverse probability of treatment weights estimated from a generalized linear model, CBPS uses the covariate balance propensity score weights, SENT uses the constrained optimal solution of the shifted relative entropy, and BENT uses the constrained optimal solution of the binary relative entropy.

	Outcome Scenario	Treatment Assignment Scenario	IPW	CBPS	SENT	BENT
Avg. Estimate (MC Std. Err.)	a	a	19.60 (3.28)	20.02 (0.55)	20.02 (0.53)	20.02 (0.51)
	a	b	20.40 (2.70)	20.00 (0.50)	20.01 (0.49)	20.00 (0.49)
	b	a	19.36 (5.27)	19.67 (3.81)	19.75 (2.90)	19.90 (2.21)
	b	b	15.33 (3.57)	15.26 (3.28)	15.83 (2.73)	16.47 (2.31)
Mean Squared Error (Bias)	a	a	10.91 (-0.40)	0.30 (0.02)	0.28 (0.02)	0.26 (0.02)
	a	b	7.42 (0.40)	0.25 (0.00)	0.24 (0.01)	0.24 (0.00)
	b	a	28.21 (-0.64)	14.63 (-0.33)	8.49 (-0.25)	4.87 (-0.10)
	b	b	34.54 (-4.67)	33.14 (-4.74)	24.79 (-4.17)	17.83 (-3.53)

Table 2. Table 2: Average estimate, Monte Carlo standard error, residual mean squared error, and empirical bias of the linear conditional ATE using the five methods for estimating balancing weights described in Section 5.2 . AIPW uses an augmented inverse probability of treatment approach where the propensity scores are estimated from a generalized linear model, CAL uses the calibration estimated weights, iCBPS uses the covariate balance propensity score weights, hdCBPS is an augmented version of CBPS, and SENT uses the constrained optimal solution of the shifted relative entropy.

	Outcome Scenario	Treatment Assignment Scenario	AIPW	CAL	iCBPS	hdCBPS	SENT
Avg. Estimate (MC Std. Err.)	a	a	19.93 (1.51)	19.93 (1.49)	19.06 (3.94)	19.93 (1.51)	19.93 (1.49)
	a	b	20.09 (1.48)	20.09 (1.47)	18.05 (3.50)	20.09 (1.47)	20.09 (1.47)
	b	a	19.80 (7.35)	20.43 (2.76)	19.72 (3.63)	19.59 (3.02)	20.25 (2.78)
	b	b	15.04 (4.52)	16.91 (2.55)	14.21 (4.48)	16.15 (2.76)	16.77 (2.54)
Mean Squared Error (Bias)	a	a	2.28 (-0.07)	2.22 (-0.07)	16.41 (-0.94)	2.28 (-0.07)	2.23 (-0.07)
	a	b	2.18 (0.09)	2.16 (0.09)	16.08 (-1.95)	2.17 (0.09)	2.16 (0.09)
	b	a	53.98 (-0.20)	7.77 (0.43)	13.26 (-0.28)	9.28 (-0.41)	7.78 (0.25)
	b	b	45.04 (-4.96)	16.01 (-3.09)	53.58 (-5.79)	22.45 (-3.85)	16.84 (-3.23)

Table 3. Table 3: ATT estimates for 30-day unplanned readmission in thoracoscopic versus open lung resection patients. UN denotes the unadjusted results, PSM denotes the propensity score matched results, IPW denotes the inverse probability weighted results, EB denotes the entropy balancing results.

Balancing Method	Risk Difference	Std. Error	95% Confidence Interval
UN	-0.021	0.006	(-0.032, -0.010)
PSM	-0.010	0.006	(-0.021, 0.001)
IPW	-0.015	0.006	(-0.027, -0.004)
EB	-0.016	0.006	(-0.027, -0.004)

Equations

$${\hat{\tau }}_{\text{ATE}}=\frac{1}{n}\sum _{i=1}^n\left[\frac{Z_i{Y}_i}{\pi ({\mathbf{X}}_i)}-\frac{(1-Z_i)Y_i}{1-\pi ({\mathbf{X}}_i)}\right]$$

$${\hat{\tau }}_{\text{ATT}}=\frac{1}{n_1}\sum _{i=1}^n\left[Z_i{Y}_i-\frac{\pi ({\mathbf{X}}_i)(1-Z_i)Y_i}{1-\pi ({\mathbf{X}}_i)}\right]$$

$$\hat{\tau }=\sum _{i=1}^n\frac{(2Z_i-1)p({\mathbf{X}}_i)Y_i}{{\sum }_{i=1}^{n}p({\mathbf{X}}_i)Z_i},$$

$$D_f(\mathbf{p}\Vert \mathbf{q})\equiv f(\mathbf{p})-f(\mathbf{q})-[\nabla f(\mathbf{q}){]}^T(\mathbf{p}-\mathbf{q}).$$

$$D_f(\mathbf{p}\Vert \mathbf{q})=\sum _{i=1}^n\left[p_i\log \left(\frac{p_i}{q_i}\right)-p_i+q_i\right].$$

$$D_f(\mathbf{p}\Vert \mathbf{q})=\sum _{i=1}^n\frac{(p_i-q_i{)}^2}{2}.$$

$$\begin{array}{rl}\text{minimize}& \sum _{i=1}^n{D}_f(p_i\Vert q_i)\\ \text{subject to}& {\mathbf{A}}^T\mathbf{p}=\mathbf{b}\end{array}$$

$$\hat{\mathbf{p}}\equiv \underset{\mathbf{p}\in \Omega \cap {\bar{\Delta }}^n}{\mathrm{argmin}}{D}_f(\mathbf{p}\Vert \mathbf{q}),$$

$$L(\mathbf{p},\mathbf{q},\mathbf{\lambda })\equiv D_f(\mathbf{p}\Vert \mathbf{q})+({\mathbf{A}}^T\mathbf{p}-\mathbf{b}{)}^T\mathbf{\lambda }.$$

$$P_f\left(\mathbf{q},\mathbf{A}\hat{\mathbf{\lambda }}\right)\equiv (\nabla f{)}^{-1}\left(\nabla f(\mathbf{q})-\mathbf{A}\hat{\mathbf{\lambda }}\right)$$

$$\begin{array}{rl}\text{maximize}& L\left[P_f(\mathbf{q},\mathbf{A}\mathbf{\lambda }),\mathbf{q},\mathbf{\lambda }\right]\\ \text{subject to}& \mathbf{\lambda }\in {\Re }^m.\end{array}$$

$$\begin{array}{rl}\text{minimize}& \sum _{i=1}^n{D}_f(p_i\Vert q_i)\\ \text{subject to}& \sum _{i=1}^n{p}_i(2Z_i-1)c_j({\mathbf{X}}_i)=0\text{and}\\ & \sum _{i=1}^n{p}_i{Z}_i{c}_j({\mathbf{X}}_i)=\sum _{i=1}^n{q}_i{c}_j({\mathbf{X}}_i)\text{for all}j=1,2,\dots ,m.\end{array}$$

$$\begin{array}{rl}{L}_{\text{ATE}}(\mathbf{p},\mathbf{q},\mathbf{\lambda })& \equiv \sum _{i=1}^n{D}_f(p_i\Vert q_i)+\sum _{j=1}^m{\lambda }_{j0}\left[\sum _{i=1}^n{p}_i(2Z_i-1)c_j({\mathbf{X}}_i)\right]\\ & +\sum _{j=1}^m{\lambda }_{j1}\left[\sum _{i=1}^n{p}_i{Z}_i{c}_j({\mathbf{X}}_i)-\sum _{i=1}^n{q}_i{c}_j({\mathbf{X}}_i)\right]\end{array}$$

$$\hat{\mathbf{\lambda }}=\underset{\mathbf{\lambda }\in {\Re }^{2m}}{\mathrm{argmax}}\sum _{i=1}^n{L}_{\text{ATE}}\left\{P_f\left[q_i,\sum _{j=1}^m(2Z_i-1)c_j({\mathbf{X}}_i){\lambda }_{j0}+\sum _{j=1}^m{Z}_i{c}_j({\mathbf{X}}_i){\lambda }_{j1}\right],q_i,\mathbf{\lambda }\right\}.$$

$$\hat{p}({\mathbf{X}}_i)=P_f\left[q_i,\sum _{j=1}^m(2Z_i-1)c_j({\mathbf{X}}_i){\hat{\lambda }}_{j0}+\sum _{j=1}^m{Z}_i{c}_j({\mathbf{X}}_i){\hat{\lambda }}_{j1}\right],i=1,2,\dots ,n.$$

$$\begin{array}{rl}\text{minimize}& \sum _{i=1}^n{D}_f(p_i\Vert q_i)\\ \text{subject to}& \sum _{i=1}^n{p}_i(1-Z_i)c_j({\mathbf{X}}_i)=\sum _{i=1}^n{q}_i{Z}_i{c}_j({\mathbf{X}}_i)\text{for all}j=1,2,\dots ,m.\end{array}$$

$$P_f\left[q_i,\sum _{j=1}^m{c}_j({\mathbf{X}}_i){\lambda }_j\right]=\frac{\pi ({\mathbf{X}}_i)}{1-\pi ({\mathbf{X}}_i)}$$

$$L_{\text{ATT}}(\mathbf{p},\mathbf{q},\mathbf{\lambda })\equiv \sum _{i=1}^n{D}_f(p_i\Vert q_i)+\sum _{j=1}^m{\lambda }_j\left[\sum _{i=1}^n{p}_i(1-Z_i)c_j({\mathbf{X}}_i)-\sum _{i=1}^n{q}_i{Z}_i{c}_j({\mathbf{X}}_i)\right].$$

$$\begin{array}{rl}\text{minimize}& \sum _{i=1}^n{D}_f(p_i\Vert q_i)\\ \text{subject to}& \sum _{i=1}^n{p}_i(2Z_i-1)c_j({\mathbf{X}}_i)=0\text{for all}j=1,2,\dots ,m\text{and}\\ & p_i=q_i\text{for all}i\in \{i:Z_i=1\}.\end{array}$$

$$\hat{p}({\mathbf{X}}_i)=P_f\left[q_i,(1-Z_i)\sum _{j=1}^m{c}_j({\mathbf{X}}_i){\hat{\lambda }}_j\right],i=1,2,\dots ,n,$$

$$\hat{\mathbf{\lambda }}=\underset{\mathbf{\lambda }\in {\Re }^m}{\mathrm{argmax}}\sum _{i=1}^n{L}_{\text{ATT}}\left\{P_f\left[q_i,(1-Z_i)\sum _{j=1}^m{c}_j({\mathbf{X}}_i){\lambda }_j\right],q_i,\mathbf{\lambda }\right\}.$$

$$\begin{array}{rl}\text{minimize}& \sum _{i=1}^n\left[p_i\log \left(\frac{p_i}{q_i}\right)-p_i+q_i\right]\\ \text{subject to}& \sum _{i=1}^n{p}_i(1-Z_i)c_j({\mathbf{X}}_i)=\sum _{i=1}^n{q}_i{Z}_i{c}_j({\mathbf{X}}_i)\text{for all}j=1,2,\dots ,m.\end{array}$$

$$\hat{\mathbf{\lambda }}=\underset{\mathbf{\lambda }\in {\Re }^m}{\mathrm{argmax}}\sum _{i=1}^n\left\{-q_i\exp \left[-(1-Z_i)\sum _{j=1}^m{c}_j({\mathbf{X}}_i){\lambda }_j\right]-q_i{Z}_i\sum _{j=1}^m{c}_j({\mathbf{X}}_i){\lambda }_j\right\}.$$

$$\hat{p}({\mathbf{X}}_i)=q_i\exp \left[-(1-Z_i)\sum _{j=1}^m{c}_j({\mathbf{X}}_i){\hat{\lambda }}_j\right],i=1,2,\dots ,n.$$

$$\pi ({\mathbf{X}}_i)=\frac{\exp \left[{\sum }_{j=1}^m{c}_j({\mathbf{X}}_i){\lambda }_j\right]}{1+\exp \left[{\sum }_{j=1}^m{c}_j({\mathbf{X}}_i){\lambda }_j\right]},i=1,2,\dots ,n,$$

$$\sum _{i=1}^n\left[\frac{Z_i{c}_j({\mathbf{X}}_i)}{\pi ({\mathbf{X}}_i)}-\frac{(1-Z_i)c_j({\mathbf{X}}_i)}{1-\pi ({\mathbf{X}}_i)}\right]=0\text{for all}j=1,2,\dots ,m.$$

$$\sum _{i=1}^n\left\{1+\exp \left[-(2Z_i-1)\sum _{j=1}^m{c}_j({\mathbf{X}}_i){\lambda }_j\right]\right\}(2Z_i-1)c_j({\mathbf{X}}_i)=0\text{for all}j=1,2,\dots ,m.$$

$$\begin{array}{rl}\text{minimize}& \sum _{i=1}^n\left[(p_i-1)\log \left(\frac{p_i-1}{q_i-1}\right)-p_i+q_i\right]\\ \text{subject to}& \sum _{i=1}^n{p}_i(2Z_i-1)c_j({\mathbf{X}}_i)=0\text{for all}j=1,2,\dots ,m.\end{array}$$

$$\hat{\mathbf{\lambda }}=\underset{\mathbf{\lambda }\in {\Re }^m}{\mathrm{argmax}}\sum _{i=1}^n\left\{(2Z_i-1)\sum _{j=1}^m{c}_j({\mathbf{X}}_i){\lambda }_j-\exp \left[-(2Z_i-1)\sum _{j=1}^m{c}_j({\mathbf{X}}_i){\lambda }_j\right]\right\}.$$

$$\hat{p}({\mathbf{X}}_i)=1+\exp \left[-(2Z_i-1)\sum _{j=1}^m{c}_j({\mathbf{X}}_i){\hat{\lambda }}_j\right],i=1,2,\dots ,n.$$

Full text

Keywords – Average Treatment Effect, Bregman Divergence, Causal Inference, Covariate Balance

1 Introduction

Causal inference has been a topic of intense interest in the statistical literature. The focus of causal inference methodology deals with the issue of how to properly evaluate treatment effects in a non-randomized setting. In many medical and scientific studies, randomization cannot be performed due to logistic, economic, and/or ethical limitations. Under these circumstances, the emergent challenge is the consistent evaluation of treatment effects in the presence of confounding. For example, in Bhagat et al., (2017), a cohort of patients undergoing lung resection surgery are examined to compare the rate of unplanned readmission following thoracoscopic versus open anatomic lung resections. The mode of surgery is not randomized and there are several pre-operative characteristics which inform the method of surgery, that in turn affect the readmission rate. Methods for causal inference seek to minimize any bias induced by these confounding variables.

Two important concepts in causal inference are the potential outcomes model (Splawa-Neyman et al.,, 1990; Rubin,, 1974) and the propensity score (Rosenbaum and Rubin,, 1983). The potential outcomes approach provides a powerful tool for conceptualizing, estimating, and performing inference regarding causal effects. An overview for implementing the potential outcomes model can be found in Imbens and Rubin, (2015). They demonstrate that a natural quantity which regularly arises when balancing potential confounders between experimental groups in observational studies is the propensity score (Rosenbaum and Rubin,, 1983). The propensity score is defined as the probability of receiving treatment given a set of measured covariates. Based on the assumptions underlying the potential outcomes model and the propensity score, causal inference proceeds in the following stages: (a) a propensity score model is fit using the observed data; (b) diagnostics for covariate balance using the propensity score are evaluated; and (c) estimates of the causal effect are produced by conditioning on the propensity score. Iterating between steps (a) and (b) is often necessary to ensure the homogeneity of the propensity score adjusted covariate distributions.

A key goal for the propensity score model is to achieve covariate balance, which means that the distribution of confounders between the treated and control groups are equivalent. From Rosenbaum and Rubin, (1983), the assumptions of strongly ignorable treatment assignment (defined in Section 2.2), in conjunction with the definition of the propensity score, imply that adjustment on the propensity score alone will theoretically achieve balance. However, this result is based on the population propensity score and does not necessarily hold in finite samples. There have been numerous approaches that address the issue of balancing empirical covariate distributions using weighting estimators. We refer to the weights produced by these methods as balancing weights. One popular method is to construct propensity scores with covariate balance built into the estimation procedure. Imai and Ratkovic, (2013) and Fan et al., (2016) introduced the covariate balance propensity score (CBPS) and its subsequent improvement (iCBPS), both of which use generalized methods of moments to fit a logit model with covariate balance serving as an auxiliary condition. Any resulting estimate of the propensity score will automatically achieve balance by construction.

In the political science literature, Hainmueller, (2012) uses maximum entropy density estimation to find balancing weights to estimate the average treatment effect of the treated. The algorithm, termed entropy balancing, finds the vector of balancing weights that minimize the normalized relative entropy from a vector of sampling weights subject to a set of linear equality constraints about the moments of the covariate distribution. Recent work by Zhao and Percival, (2017) shows how this algorithm enjoys a double-robustness property. The general idea of double-robust estimation is to combine covariate information about the treatment assignment and the outcome model into the weighting estimator (Bang and Robins,, 2005; Kang and Schafer,, 2007). If at least one model is correctly specified, then the resulting causal effect estimate is consistent. When both the outcome and treatment models are correctly specified, then the estimate achieves the semiparametric efficiency bound described by Hahn, (1998). Entropy balancing is limited to finding balancing weights to estimate the average treatment effect of the treated, leaving gaps in the procedure for developing doubly-robust estimators of other estimands. This issue is related to the choice of Hainmueller, (2012) to optimize the normalized relative entropy. By changing the criterion distance function in a convex optimization problem, similar to the one presented by Hainmueller, (2012), we can draw parallels to other covariate balancing methods. A similar idea is proposed by Zhao, (2019), who shows that CBPS and entropy balancing can be generalized by modifying the score function derived from the respective covariate balance problem. Calibration estimators (Deville and Sarndal,, 1992) also produce balancing weights using constrained convex optimization techniques. The proposed methods in Chan et al., (2015) implicitly extends entropy balancing to include other distance functions. However, they restrict their attention to a nonparametric setting, characterizing their methodology as a departure from the propensity score literature.

Our aim is to extend the work of Hainmueller, (2012), Imai and Ratkovic, (2013), and Fan et al., (2016) for finding balancing weights that facilitate causal effect estimation when the treatment assignment is not determined by a logit model. We do so by demonstrating how balancing weights can be computed from Bregman distances (Bregman,, 1967). Bregman distances have multiple geometric properties that allow for easy estimation of the balancing weights. This geometric interpretation of balancing weights complements the implicit geometry found in classic semiparametric inference. Using the results of our framework, we prove that CBPS (Imai and Ratkovic,, 2013) and iCBPS (Fan et al.,, 2016) are doubly-robust estimators of the average treatment effect while assuming the propensity scores follow a logit model. As an extension to CBPS, we propose an estimator for balancing weights akin to the overlap weights discussed by Li et al., (2018). We also show how our framework is consistent with the calibration estimator approach of Chan et al., (2015), thereby bridging the empirical covariate balancing methods of entropy balancing, CBPS, iCBPS, and calibration estimators. We are interested in these methods in particular as they do not incorporate a model of the outcome process into their designs in the spirit of Rubin, (2008).

The outline of this article is as follows. Section 2 defines the general notation and assumptions that will be applied throughout the manuscript. Section 3 describes the methods for finding balancing weights by solving a constrained optimization problem using Bregman distances as the criterion function. Section 4 describes the similarities between our method and other covariate balancing methods. Section 5 summarizes results from two simulation studies comparing different covariate balancing methods. This section also contains the results for a replication study of Bhagat et al., (2017) using a variety of different covariate balancing methods. The real data set illustrates the importance of selecting appropriate covariate balancing methods. Finally, Section 6 concludes with a discussion of the framework and future work.

2 Background and Preliminaries

2.1 Notation and Definitions

Parameters will be denoted using Greek letters, whereas random variables will be denoted with Roman letters. Boldface letters will denote vectors and matrices while non-boldface letters represent scalars. For a matrix $\mathbf{A}$ the transpose is written as $\mathbf{A}^{T}$. The symbol $\nabla f$ denotes the gradient of a function $f$. Let $\mathbf{1}_{n}$ denote the $(n\times 1)$ vector with each entry equal to one. Similarly, let $\mathbf{0}_{n}$ denote the $(n\times 1)$ vector with each entry equal to zero.

Let $\mathbf{X}$ denote a vector of real-valued covariate measurements, $Z$ denote the random treatment assignment with support $\{0,1\}$, and $Y$ denote the real-valued outcome variable. The independent sampling units will be indexed by $i=1,2,\ldots,n$. The $(n\times 1)$ vector of balancing weights will be written as $\mathbf{p}\equiv(p_{1},p_{2},\ldots,p_{n})^{T}$ while the $(n\times 1)$ vector of sampling weights will be written as $\mathbf{q}\equiv(q_{1},q_{2},\ldots,q_{n})^{T}$. We will often write $p_{i}=p(\mathbf{X}_{i})$ for the $i$th subject to emphasize the fact that the balancing weights are conditioned on the covariates. Define $\{c_{j}(\mathbf{X});j=1,2,\ldots,m\}$, as a set of functions that generate linearly independent features to be balanced between treatment groups. We will refer to these quantities as balance functions.

2.2 Potential Outcomes Model

Potential outcomes provide a convenient framework for conceptualizing causal effects. This framework was first introduced by Splawa-Neyman et al., (1990) for randomized experiments. The concepts and assumptions necessary to extend this framework to observational data were later formalized by Rubin, (1974). The potential outcomes are denoted with a vector $\left[Y(0),Y(1)\right]^{T}$ with $Y(0)$ and $Y(1)$ corresponding to the counterfactual outcome when $Z=0$ and $Z=1$, respectively. The conditional expectations for the potential outcomes are denoted with $\mu_{0}(\mathbf{X})\equiv\mathbb{E}[Y(0)|\mathbf{X}]$ and $\mu_{1}(\mathbf{X})\equiv\mathbb{E}[Y(1)|\mathbf{X}]$. The random outcome is defined by the transformation $Y\equiv ZY(1)+(1-Z)Y(0)$. Some common causal estimands are the population average treatment effect (ATE), $\tau_{\text{ATE}}\equiv\mathbb{E}\left[Y(1)-Y(0)\right]$, and the population average treatment effect of the treated (ATT), $\tau_{\text{ATT}}\equiv\mathbb{E}\left[Y(1)-Y(0)|Z=1\right]$. In any case, the causal effects are non-identifiable as one of the two required potential outcomes is always missing. This simple observation is the fundamental problem of causal inference. We adopt the setting proposed by Rosenbaum and Rubin, (1983) who describe a set of assumptions that will allow us to find consistent estimates of the treatment effect in observational studies. This includes the following assumptions about the data.

Assumption 1** (Strong Ignorability):**

[Y(0),Y(1)]^{T}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z|\mathbf{X}.

The strong ignorability assumption requires that the vector of potential outcomes be independent of the treatment assignment when we condition on the covariates. This assumption further implies that there is no unmeasured confounding. The implication of Assumption 1 along with the definition of the propensity score as a balance criterion allows us to conclude that [Y(0),Y(1)]^{T}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z|\pi(\mathbf{X}), where $\pi(\mathbf{X})\equiv\Pr\{Z=1|\mathbf{X}\}$ denotes the propensity score (Rosenbaum and Rubin,, 1983).

Assumption 2** (Positivity):**

$0<\Pr\{Z=1|\mathbf{X}\}<1$* for all $\mathbf{X}$.*

The treatment positivity assumption requires the probability that a subject is assigned to the treatment group as opposed to the control group be bounded away from zero and one. Since $\Pr\{Z=1|\mathbf{X}\}$ must be estimated from the covariates, then Assumption 2 equivalently amounts to requiring sufficient overlap between the covariate distributions of the two treatment groups. The feasibility of the convex optimization problems that we will introduce later on are deeply intertwined with Assumption 2. Without sufficient overlap, the estimated balancing weights will either not exist, or be unstable and produce biased estimates of the causal effect.

2.3 Horvitz-Thompson Estimator

The Horvitz-Thompson class of estimators (Horvitz and Thompson,, 1952; Hirano et al.,, 2003) frequently appears in the causal inference literature. For example, the Horvitz-Thompson estimator for the average treatment effect is

[TABLE]

while the Horvitz-Thompson estimator for $\tau_{\text{ATT}}$ is

[TABLE]

where $n_{1}=\sum_{i=1}^{n}Z_{i}$. Hahn, (1998) was able to show that (1) and (2) have optimal asymptotic properties for estimating $\tau_{\text{ATE}}$ and $\tau_{\text{ATT}}$. Even when we substitute a consistent estimator of the propensity score into (1), the estimator for $\tau_{\text{ATE}}$ remains consistent and achieves the semiparametric efficiency bound.

A more general form for causal effect estimation is

[TABLE]

which accommodates several different estimands through the choice of $p(\mathbf{X})$. For example, we will see in Section 3.4 that the estimator for $\tau_{\text{ATT}}$ is similar to estimators of $\tau_{\text{ATE}}$ with additional constraints placed on the balancing weights so that $p(\mathbf{X})=q$ whenever $Z=1$. Equations (1) and (2) provide direction for identifying $p(\mathbf{X})$ within (3) in order to estimate $\tau_{\text{ATE}}$ and $\tau_{\text{ATT}}$, respectively. If the propensity score is known, $\tau_{\text{ATE}}$ can be estimated by setting $p(\mathbf{X})=\pi(\mathbf{X})^{-1}$ when $Z=1$ and $p(\mathbf{X})=[1-\pi(\mathbf{X})]^{-1}$ when $Z=0$. We can also find an estimator for $\tau_{\text{ATT}}$ by setting $p(\mathbf{X})=\pi(\mathbf{X})[1-\pi(\mathbf{X})]^{-1}$ when $Z=0$ and $p(\mathbf{X})=1$ when $Z=1$. When the propensity score is unknown, finding an estimator for $p(\mathbf{X})$ that produces consistent estimates of $\tau_{\text{ATE}}$ and $\tau_{\text{ATT}}$ is relatively straightforward. Estimating balancing weights that also preserve the efficiency of $\hat{\tau}_{\text{ATE}}$ and $\hat{\tau}_{\text{ATT}}$ is a more challenging proposition.

3 Bregman Distances

3.1 Definition

Let $\Delta^{n}\subseteq\Re^{n}$ be a non-empty, convex, and open set with closure $\bar{\Delta}^{n}$. Define $f:\bar{\Delta}^{n}\rightarrow\Re$ to be a continuously differentiable, strictly convex function. The Bregman distance generated by the function $f$ is the difference between $f$ evaluated at $\mathbf{p}\in\bar{\Delta}^{n}$ and the first-order Taylor series approximation of $f$ about $\mathbf{q}\in\Delta^{n}$, evaluated at $\mathbf{p}$. In other words, a Bregman distance $D_{f}:\bar{\Delta}^{n}\times\Delta^{n}\rightarrow\Re$ may be defined as

[TABLE]

Bregman distances are often used to measure the convexity associated with $f$. Since $f$ is strictly convex over $\bar{\Delta}^{n}$, it follows that for $\mathbf{p}\in\bar{\Delta}^{n}$ and $\mathbf{q}\in\Delta^{n}$, $D_{f}(\mathbf{p}\;\|\;\mathbf{q})\geq 0$ with equality holding if and only if $\mathbf{p}=\mathbf{q}$. This implies that $D_{f}(\mathbf{p}\;\|\;\mathbf{q})$ is also strictly convex. A more complete definition of Bregman distances can be found in Chapter 2 of Censor and Zenios, (1998), which includes additional properties that $D_{f}$ must satisfy which are not mentioned here. A visual representation of a Bregman distance can be found in Figure 1.

One of the most common examples of a Bregman distance is the unnormalized relative entropy. Let $f(\mathbf{p})=\sum_{i=1}^{n}p_{i}\log\left(p_{i}\right)$ for $\mathbf{p}\in[0,\infty)^{n}$. We assume $0\log(0)=0$ so that the domain of $f$ includes the boundary points contained within the closure of $\Delta^{n}$. The resulting Bregman distance is written as

[TABLE]

The Euclidean distance is another example of a Bregman distance. By selecting $f(\mathbf{p})=\sum_{i=1}^{n}p_{i}^{2}/2$ for $\mathbf{p}\in\Re^{n}$ we get

[TABLE]

In order to simplify the presentation of the methods, we will only consider Bregman distances that are separable. This means $D_{f}(\mathbf{p}\;\|\;\mathbf{q})=\sum_{i=1}^{n}D_{f}(p_{i}\;\|\;q_{i})$. Note that both the unnormalized relative entropy and the Euclidean distance are separable. Since we require positive weights, we also restrict our focus to convex functions where $\Delta^{n}\subseteq[0,\infty)^{n}$ in order to avoid setting additional constraints for $\mathbf{p}\geq\mathbf{0}_{n}$. Notice that the domain of the unnormalized relative entropy satisfies this condition while the domain of the Euclidean distance does not. In addition, we will assume throughout that the sampling weights $\mathbf{q}\in\Delta^{n}$ are fixed by design and known.

3.2 Constrained Optimization and Duality

For some $\mathbf{q}\in\Delta^{n}$, the value $\hat{\mathbf{p}}\in\bar{\Delta}^{n}$ that minimizes $D_{f}(\mathbf{p}\;\|\;\mathbf{q})$ in an unconstrained setting is $\hat{\mathbf{p}}=\mathbf{q}$. In covariate balance problems, we specify a set of linear constraints that the optimal solution must satisfy. Consider the constrained convex optimization problem to

[TABLE]

where $\mathbf{A}$ is a linearly independent $(n\times m)$ matrix that forms the basis of a linear subspace that defines the constraints of the program and $\mathbf{b}$ is an $(m\times 1)$ vector denoting the margins of those constraints. The entries of $\mathbf{A}$ and $\mathbf{b}$ are denoted with $a_{ij}\in\Re$ and $b_{j}\in\Re$ ($i=1,2,\ldots,n$ and $j=1,2,\ldots,m$), respectively. Equation (4) is often referred to as the primal problem and the corresponding solution is referred to as the primal solution. We denote the set of feasible primal solutions that satisfy the linear constraints in (4) as $\Omega\equiv\left\{\mathbf{p}:\mathbf{A}^{T}\mathbf{p}=\mathbf{b}\right\}$.

Geometrically, the solution to the primal problem is the point

[TABLE]

which is the generalized projection of $\mathbf{q}\in\Delta^{n}$ into $\Omega$. Note that $\Omega\cap\bar{\Delta}^{n}$ is sometimes empty. One solution to avoid this issue is to choose $\mathbf{b}=\mathbf{A}^{T}\tilde{\mathbf{p}}$ where $\tilde{\mathbf{p}}\in\bar{\Delta}^{n}$ so that $\tilde{\mathbf{p}}\in\Omega\cap\bar{\Delta}^{n}$. For example when using the unnormalized relative entropy, we might set $\mathbf{b}=\mathbf{0}_{m}$ with $\tilde{\mathbf{p}}=\mathbf{0}_{n}$. However, this condition is not so obvious for some of the covariate balancing problems that we will encounter later on. Instead, we will assume $\Omega\cap\bar{\Delta}^{n}\neq\emptyset$ throughout. Lemma 1 proves that the generalized projection is unique. The proof appears in Section 2.1 of Censor and Zenios, (1998). For the sake of completeness, a version of this proof is also found in the online supplement.

Lemma 1**.**

Suppose $\Omega\cap\bar{\Delta}^{n}\neq\emptyset$. Then the generalized projection of $\mathbf{q}$ into $\Omega$, defined in (5), is unique.

When $\hat{\mathbf{p}}\in\Omega\cap\Delta^{n}$, the primal problem can be solved by introducing a vector of Lagrangian multipliers. With Lagrangian multipliers, we can formulate the Lagrangian $L:\bar{\Delta}^{n}\times\Delta^{n}\times\Re^{m}\rightarrow\Re$ for any constrained optimization problem in the form of (4) as

[TABLE]

Optimizing the Lagrangian with respect to $\mathbf{p}\in\Delta^{n}$ and $\boldsymbol{\lambda}\in\Re^{m}$ is an unconstrained problem equivalent to (4). Instead of finding the point $\hat{\mathbf{p}}\in\Omega\cap\Delta^{n}$ that minimizes $D_{f}(\mathbf{p}\;\|\;\mathbf{q})$, we find the vector $\hat{\mathbf{p}}$ that minimizes the Lagrangian with respect to $\mathbf{p}\in\Delta^{n}$ and $\hat{\boldsymbol{\lambda}}$ that maximizes the Lagrangian with respect to $\boldsymbol{\lambda}\in\Re^{m}$. In other words, the optimal solution forms a saddle point on (6) over the space $\Delta^{n}\times\Re^{m}$ for a fixed $\mathbf{q}\in\Delta^{n}$.

The following Propositions are used in tandem to obtain balancing weights for treatment effect estimation. A necessary condition for these Propositions is that the function $f$ which generates the Bregman distance $D_{f}$ be zone consistent with respect to $\Omega$. This means that for any $\mathbf{q}\in\Delta^{n}$, the Bregman distance produced by $f$ has its generalized projection of $\mathbf{q}$ into $\Omega$ contained within the open set $\Delta^{n}$.

Proposition 1**.**

Assume that $f$ is zone consistent with respect to $\Omega$. Let $\hat{\mathbf{p}}\in\Omega\cap\Delta^{n}$ be the generalized projection of $\mathbf{q}$ into $\Omega$, as defined in (5). Then $\hat{\mathbf{p}}$ is uniquely determined by

[TABLE]

where $\hat{\boldsymbol{\lambda}}\in\Re^{m}$ is also unique.

The proof of Proposition 1 can be found in the online supplement. Given the result of Proposition 1, the convex optimization problem can be solved by estimating $\hat{\boldsymbol{\lambda}}\in\Re^{m}$ with the dual problem, which is to

[TABLE]

Propositions 1 and 2 imply that the primal solution to (4) can be obtained by plugging the solution to the dual problem into the generalized projection. The proof of Proposition 2 is adapted from Section 3.4 of Bertsekas, (1999) and appears in the online supplement.

Proposition 2**.**

Assume $f$ is zone consistent with respect to $\Omega$. If the primal problem (4) has an optimal solution, then the dual problem (7) also has an optimal solution and the two optimal values are equal.

3.3 Balancing Weights for the ATE

In this section we outline the general strategy and guidelines for obtaining balancing weights for estimating $\tau_{\text{ATE}}$. The proposed method requires solving the primal problem to

[TABLE]

As mentioned in the previous section, (8) can be solved by optimizing the corresponding Lagrangian, which is defined as

[TABLE]

where $\boldsymbol{\lambda}_{0}\equiv(\lambda_{10},\lambda_{20},\ldots,\lambda_{m0})^{T}$ and $\boldsymbol{\lambda}_{1}\equiv(\lambda_{11},\lambda_{21},\ldots,\lambda_{m1})^{T}$ with $\boldsymbol{\lambda}\equiv(\boldsymbol{\lambda}^{T}_{0},\boldsymbol{\lambda}^{T}_{1})^{T}$. The criterion distance function should be selected so that $P_{f}\left[q_{i},\sum_{j=1}^{m}c_{j}(\mathbf{X}_{i})\lambda_{j0}\right]=\pi(\mathbf{X}_{i})^{-1}$.

We can also frame this problem using the notation from Section 3.2. Let $\mathbf{A}_{0}$ be an $(n\times m)$ matrix whose elements consist of $a_{ij0}=(2Z_{i}-1)c_{j}(\mathbf{X}_{i})$, $\mathbf{A}_{1}$ be an $(n\times m)$ with entries $a_{ij1}=Z_{i}c_{j}(\mathbf{X}_{i})$, $\mathbf{b}_{0}=\mathbf{0}_{m}$, and $\mathbf{b}_{1}$ be an $(m\times 1)$ vector with entries $b_{j1}=\sum_{i=1}^{n}q_{i}c_{j}(\mathbf{X}_{i})$. We then combine $\mathbf{A}_{0}$ and $\mathbf{A}_{1}$ to construct $\mathbf{A}=\left[\mathbf{A}_{0},\mathbf{A}_{1}\right]$ while $\mathbf{b}_{0}$ and $\mathbf{b}_{1}$ are concatenated into $\mathbf{b}=\left[\mathbf{b}^{T}_{0},\mathbf{b}^{T}_{1}\right]^{T}$.

After differentiating (9) with respect to $p_{i}$ for some $i=1,2,\ldots,n$ and setting the resulting derivative to zero, we arrive at the generalized projection evaluated at $\boldsymbol{\lambda}\in\Re^{2m}$. The dual objective function is obtained by substituting the generalized projection for $p_{i}$ in the Lagrangian. The dual solution solves the dual problem,

[TABLE]

As a result of Propositions 1 and 2, the balancing weights are uniquely determined by the generalized projection evaluated at $\hat{\boldsymbol{\lambda}}\in\Re^{2m}$,

[TABLE]

It is recommended, and often necessary (see the proofs to Theorems 1 and 2), that one of the balance functions be an intercept - we will assume throughout that $c_{1}(\mathbf{X})=1$ for all $\mathbf{X}$. This constraint implies $\sum_{\{i:Z_{i}=1\}}\hat{p}_{i}=\sum_{\{i:Z_{i}=0\}}\hat{p}_{i}$.

3.4 Balancing Weights for the ATT

Next we consider the problem of finding balancing weights to estimate $\tau_{\text{ATT}}$. This requires solving the primal problem to

[TABLE]

The criterion Bregman distance should be chosen so that the generalized projection resembles the functional form for the odds of treatment. That is,

[TABLE]

where $\boldsymbol{\lambda}\equiv(\lambda_{1},\lambda_{2},\ldots,\lambda_{m})^{T}$. In terms of the notation presented in Section 3.2, (12) is equivalent to (4) by setting $a_{ij}=(1-Z_{i})c_{j}(\mathbf{X}_{i})$, $i=1,2,\ldots,n$, and $b_{j}=\sum_{i=1}^{n}q_{i}Z_{i}c_{j}(\mathbf{X}_{i})$, $j=1,2,\ldots,m$. Similar to the balancing weights for estimating $\tau_{\text{ATE}}$, we set $c_{1}(\mathbf{X})=1$ for all $\mathbf{X}$ so that $\sum_{\{i:Z_{i}=1\}}q_{i}=\sum_{\{i:Z_{i}=0\}}\hat{p}_{i}$. The associated Lagrangian for (12) can then be expressed as

[TABLE]

Notice that when $Z=1$ and $q\in\Delta$, then $P_{f}\left[q,(1-Z)\sum_{j=1}^{m}c_{j}(\mathbf{X})\lambda_{j}\right]=q$. Moreover, recall from Section 3.1 that $D_{f}(q\;\|\;q)=0$ for some $q\in\Delta$. Therefore, (12) can be reconstructed into the equivalent primal problem to

[TABLE]

According to Propositions 1 and 2, the balancing weights are evaluated with the resulting generalized projection,

[TABLE]

where the dual vector is estimated by solving for

[TABLE]

4 Relationship with other Covariate Balance Methods

4.1 Entropy Balancing

Entropy balancing (Hainmueller,, 2012) is a special case of a constrained Bregman distance optimization problem. By setting $f(\mathbf{p})=\sum_{i=1}^{n}p_{i}\log(p_{i})$ and $q_{i}\in(0,\infty)$ for all $i=1,2,\ldots,n$, we can identify the entropy balancing primal problem, which is to

[TABLE]

According to Proposition 2, optimizing (13) is equivalent to maximizing the dual objective function,

[TABLE]

The vector of balancing weights is obtained by evaluating the generalized projection with the solution to the dual problem, which yields

[TABLE]

In Hainmueller, (2012), (13) is written using the normalized relative entropy instead of the unnormalized relative entropy. However, optimizing the normalized relative entropy is simply achieved by (14) and (15) with minor alterations. Let $q^{\prime}_{i}=q_{i}/\sum_{i=1}^{n}q_{i}Z_{i}$ for all $i\in\{i:Z_{i}=1\}$ and constrain the intercept so that $\sum_{i=1}^{n}\hat{p}_{i}(1-Z_{i})=\sum_{i=1}^{n}q^{\prime}_{i}Z_{i}=1$. As previously suggested in Section 3.4, we recommend setting $c_{1}(\mathbf{X})=1$. In doing so, the resulting balancing weights for the control group will sum to one while still satisfying the constraints of the primal problem.

Using the resulting estimating equations for $\boldsymbol{\lambda}$ and $\tau_{\text{ATT}}$ in concordance with results from M-estimation theory (Stefanski and Boos,, 2002), Zhao and Percival, (2017) show that entropy balancing weights produce doubly-robust estimates of $\tau_{\text{ATT}}$. This means if either $\text{logit}[\pi(\mathbf{X})]\in\text{span}\{c_{j}(\mathbf{X}):j=1,2,\ldots,m\}$ or $\mu_{0}(\mathbf{X})\in\text{span}\{c_{j}(\mathbf{X}):j=1,2,\ldots,m\}$, then the balancing weights of (15) applied to (3) is consistent for $\tau_{\text{ATT}}$. If both conditions are satisfied, then the estimator achieves the semiparametric efficiency bound derived by Hahn, (1998) for estimators of $\tau_{\text{ATT}}$. The Horvitz-Thompson estimator for $\tau_{\text{ATT}}$ that substitutes a consistent estimate of the propensity score for $\pi(\mathbf{X})$ in (2), on the other hand, does not achieve the semiparametric efficiency bound.

4.2 Covariate Balance Propensity Scores

Another method for covariate balance, developed by Imai and Ratkovic, (2013), proposes fitting a logit model for the propensity score,

[TABLE]

subject to

[TABLE]

They opted to solve for $\hat{\boldsymbol{\lambda}}\in\Re^{m}$ using generalized method of moments (GMM) while at the same time satisfying (17). The estimated propensity scores can be transformed into balancing weights for estimating $\tau_{\text{ATE}}$ with the inverse probability of treatment weighting estimator. We will refer to the model where the balance functions that appear in (17) are identical to the balance functions within the linear predictor of (16) as the exactly-specified CBPS model.

The weights obtained with an exactly-specified CBPS model can be expressed in an equivalent manner to a constrained optimization problem following our framework. First, notice that the fixed effect coefficients of the logit model can double as a vector of dual variables. Next, observe that (17) can be rewritten as

[TABLE]

The CBPS primal problem can then be constructed in order to

[TABLE]

We call the criterion distance function in (18) the shifted relative entropy which is generated by setting $f(\mathbf{p})=\sum_{i=1}^{n}(p_{i}-1)\log(p_{i}-1)$, $\mathbf{p}\in[1,\infty)^{n}$. Notice that (18) also contains fewer constraints than (8). We assume $q_{i}=2$ for all $i=1,2,\ldots,m$. Assuming uniform sampling weights follows the prevailing philosophy of the causal inference literature in which observational data are typically randomly sampled from the population of interest. The solution to the dual problem for (18) finds

[TABLE]

The principal reason for selecting the shifted relative entropy as the criterion distance function are the resulting balancing weights which resemble the inverse probability of treatment weights,

[TABLE]

A similar derivation of CBPS using the dual function setup was also described by Zhao, (2019).

Fan et al., (2016) identifies a condition that the balance functions must satisfy in order for CBPS to be doubly-robust for estimating $\tau_{\text{ATE}}$. This condition is not obvious from a data analytic context. However, the condition is satisfied if we assume a constant conditional average treatment effect. Under this assumption, we can prove that CBPS is doubly-robust using the balancing weights produced by (20).

Assumption 3** (Constant Conditional ATE):**

For all $\mathbf{X}$, $\mu_{1}(\mathbf{X})-\mu_{0}(\mathbf{X})=\tau$.

Theorem 1**.**

Let Assumptions 1 and 2 be given. Suppose $\mathbb{E}[Y(0)]$, $\mathbb{E}[Y(1)]$, and $\mathbb{E}[c_{j}(\mathbf{X})]$ exist for all $j=1,2,\ldots,m$. Furthermore, assume $\mathbb{V}[Y(0)]<\infty$ and $\mathbb{V}[Y(1)]<\infty$. Then the balancing weights determined by (19) and (20) applied to (3) is doubly-robust in the sense that:

If $\text{logit}[\pi(\mathbf{X})]=\sum_{j=1}^{m}c_{j}(\mathbf{X})\lambda_{j}$ for some $\lambda_{j}\in\Re$, $j=1,2,\ldots,m$, then $\hat{\tau}$ is consistent for $\tau_{\text{ATE}}$; 2. 2.

Under Assumption 3 and if $\mu_{0}(\mathbf{X})=\sum_{j=1}^{m}c_{j}(\mathbf{X})\beta_{j}$ for some $\beta_{j}\in\Re$, $j=1,2,\ldots,m$, then $\hat{\tau}$ is consistent for $\tau_{\text{ATE}}$; 3. 3.

If conditions 1 and 2 are both satisfied, then

[TABLE]

where

[TABLE]

As an extension to CBPS, consider the primal problem using the Bregman distance generated by setting $f(\mathbf{p})=\sum_{i=1}^{n}p_{i}\log(p_{i})+(1-p_{i})\log(1-p_{i})$, $\mathbf{p}\in[0,1]^{n}$:

[TABLE]

If we assume $q_{i}=1/2$ for all $i=1,2,\ldots,n$, then according to Propositions 1 and 2 the solution to (21) is

[TABLE]

where the dual solution is obtained by solving for

[TABLE]

The Bregman distance in this case is referred to as the binary relative entropy. This distance is useful for finding balancing weights that produce estimates for a special case of the weighted average treatment effect called the optimally weighted average treatment effect (OWATE) (Crump et al.,, 2006),

[TABLE]

A consistent estimator for $\tau_{\text{OWATE}}$ is also consistent for $\tau_{\text{ATE}}$, with the smallest variance, when we are given Assumption 3 and the potential outcomes have equal variance. Li et al., (2018) further motivates the use of estimators for $\tau_{\text{OWATE}}$ when there is poor overlap between the treated and control groups. Equations (22) and (23) provide a dual interpretation of the covariate balance scoring rule for estimating $\tau_{\text{OWATE}}$ considered by Zhao, (2019). By replacing $p(\mathbf{X})$ with (22) in (3), and using arguments similar to the proof of Theorem 1, we derive a doubly-robust estimator for $\tau_{\text{OWATE}}$ with the usual asymptotic properties.

Corollary 1**.**

Under the same assumptions and conditions as Theorem 1, the balancing weights determined by (22) and (23) applied to (3) is doubly-robust for estimating $\tau_{\text{OWATE}}$ with asymptotic variance

[TABLE]

4.3 Improved Covariate Balance Propensity Scores

The iCBPS approach (Fan et al.,, 2016) improves upon the CBPS method described in Section 4.2 to better accommodate heterogeneous treatment effects. The objective of this method is to fit a logit model of the propensity score subject to the constraints

[TABLE]

Fan et al., (2016) uses GMM to estimate $\hat{\boldsymbol{\lambda}}\in\Re^{m}$ in (16) subject to (25). This modified approach can be adapted to fit into our proposed framework, with the balancing weights being estimated using dual optimization techniques instead of GMM.

Using the same criterion Bregman distance as the one used in (18), we can obtain balancing weights that satisfy (25) as follows. Assume $q_{i}=2$ for all $i=1,2,\ldots,m$. Define the primal problem for iCBPS as

[TABLE]

As opposed to (18), the iCBPS primal problem follows our guidelines in Section 3.3 more closely. The resulting dual solution solves for

[TABLE]

The covariate balancing weights differ slightly from (20) due to the additional constraints in (26) with

[TABLE]

Note that with the GMM approach, $\boldsymbol{\lambda}\in\Re^{m}$ whereas with our method, $\boldsymbol{\lambda}\in\Re^{2m}$. This implies that exact balance between covariates is not necessarily achieved with the proposed methods of Fan et al., (2016). A notable deviation from our own recommendations exists within how $\mathbf{b}$ is specified. Even though $q_{i}=2$ for all $i=1,2,\ldots,n$, we set $b_{j1}=\sum_{i=1}^{n}c_{j}(\mathbf{X}_{i})$ for all $j=1,2,\ldots,m$. If we were to follow the setup in Section 3.3, we would set $b_{j1}=\sum_{i=1}^{n}2c_{j}(\mathbf{X}_{i})$. However, this distinction is minor in the context of uniform sampling weights and should produce similar results.

We now show that the weights produced by (27) and (28) applied to (3) is doubly-robust given a linear conditional average treatment effect, defined in Assumption 4. Note that Assumption 4 is less stringent than Assumption 3 which was necessary to prove Theorem 1.

Assumption 4** (Linear Conditional ATE):**

For all $\mathbf{X}$, $\mu_{1}(\mathbf{X})-\mu_{0}(\mathbf{X})=\sum_{j=1}^{m}c_{j}(\mathbf{X})\alpha_{j}$ where $\alpha_{j}\in\Re$ for all $j=1,2,\ldots,m$.

Theorem 2**.**

Let Assumptions 1 and 2 be given. Suppose $\mathbb{E}[Y_{i}(0)]$, $\mathbb{E}[Y(1)]$ and $\mathbb{E}[c_{j}(\mathbf{X})]$ exist for all $j=1,2,\ldots,m$. Assume $\mathbb{V}[Y(0)]<\infty$ and $\mathbb{V}[Y(1)]<\infty$. Then the balancing weights determined by (27) and (28) applied to (3) is doubly-robust in the sense that:

If $\text{logit}[\pi(\mathbf{X})]=\sum_{j=1}^{m}c_{j}(\mathbf{X})\lambda_{j0}$ for some $\lambda_{j0}\in\Re$, $j=1,2,\ldots,m$, then $\hat{\tau}$ is consistent for $\tau_{\text{ATE}}$; 2. 2.

Under Assumption 4 and if $\mu_{0}(\mathbf{X})=\sum_{j=1}^{m}c_{j}(\mathbf{X})\beta_{j}$ for some $\beta_{j}\in\Re$, $j=1,2,\ldots,m$, then $\hat{\tau}$ is consistent for $\tau_{\text{ATE}}$; 3. 3.

If conditions 1 and 2 are both satisfied, then

[TABLE]

where

[TABLE]

4.4 Calibration Estimators

Chan et al., (2015) describes a class of estimators originally introduced by Deville and Sarndal, (1992) for survey sampling called calibration estimators. One of the contributions from Chan et al., (2015) shows how calibration estimators can be applied to covariate balance problems. For some separable generalized distance function $G:\Delta^{n}\rightarrow\Re$, calibration estimators find balancing weights that solve the primal problem to

[TABLE]

Chan et al., (2015) assume uniform sampling weights. Equation (29) is solved by defining the functions $h(p)\equiv G(1-p)$ and $g(v)\equiv h\left[(\nabla h)^{-1}(v)\right]+v-v(\nabla h)^{-1}(v)$ with $p\in\Delta$ and $v\in\Re$ to write the dual objective functions, which we use to solve for

[TABLE]

The resulting balancing weights are obtained by evaluating the first derivative of $g$ at $\hat{\boldsymbol{\lambda}}_{0}$ and $\hat{\boldsymbol{\lambda}}_{1}$,

[TABLE]

Similar to our solution for finding balancing weights to estimate $\tau_{\text{ATE}}$, the dual variable $\boldsymbol{\lambda}\equiv(\boldsymbol{\lambda}_{0},\boldsymbol{\lambda}_{1})^{T}$ has $2m$ entries. Tseng and Bertsekas, (1987) and Chan et al., (2015) show that (30) can be solved using any strictly concave $g(v)$, $v\in\Re$, assuming that a feasible solution for (29) exists. Therefore, calibration estimators are not necessarily restricted to Bregman distances. However, if $G$ is a monotone increasing transformation of $D_{f}$ with respect to $\mathbf{p}\in\Delta^{n}$, then (29) can be constructed using Bregman distances so that the primal solutions are equivalent. The only difference from the methods we present in Section 3.3 is with the construction of $\mathbf{A}$ and $\mathbf{b}$. In Theorem 3 we identify the conditions for which our proposed method is equivalent to the calibration estimator approach of Chan et al., (2015).

Theorem 3**.**

Suppose we have a generalized distance $G(\mathbf{p})$ which is a monotone increasing transformation of some Bregman distance $D_{f}(\mathbf{p}\;\|\;\mathbf{q})$ with respect to $\mathbf{p}\in\Delta^{n}$ and $\mathbf{q}\in\Delta^{n}$ is uniform. Then for

[TABLE]

$\tilde{\mathbf{p}}=\hat{\mathbf{p}}$* where $\tilde{\mathbf{p}}=\operatorname*{arg\,min}_{\mathbf{p}\in\Omega_{0}\cap\Delta^{n}}G(\mathbf{p})$ and $\hat{\mathbf{p}}=\operatorname*{arg\,min}_{\mathbf{p}\in\Omega_{1}\cap\Delta^{n}}D_{f}(\mathbf{p}\;\|\;\mathbf{q})$.*

Theorem 1 of Chan et al., (2015) shows that calibration estimators can be used to produce consistent estimates for $\tau_{\text{ATE}}$ while also attaining the semiparametric efficiency bound described by Hahn, (1998). This is accomplished using a nonparametric setup where the balance functions represent a basis for uniformly approximating $\mu_{0}(\mathbf{X})$, $\mu_{1}(\mathbf{X})$, and $\pi(\mathbf{X})$. Given this result and Theorem 3 implies that for a sufficiently rich set of balance functions, the Bregman distance weights in conjunction with (3) can produce consistent and efficient estimates of $\tau_{\text{ATE}}$. This result is quite useful when the balance functions that determine either the outcome or the treatment assignment are unknown.

Without further defining the distance to be optimized in the primal problem, (30) and (31) are less flexible when considering non-uniform sampling weights. This is especially important when developing iterative estimation algorithms or dealing with more complex balance designs where the data are not sampled uniformly from the population of interest. Furthermore, calibration estimators, as they are described in Chan et al., (2015), achieve a three-way balance between the treated, the controls, and the combined treatment groups for estimating $\tau_{\text{ATE}}$. As shown in Theorem 1, this condition is not required when Assumption 3 holds. Zhao, (2019) also noted that this condition is not required to achieve global efficiency using covariate balance scoring rules.

5 Numerical Studies

5.1 Homogeneous Treatment Effect Simulation

In this section, we demonstrate the utility of the proposed methodology using simulated data that assumes a constant conditional average treatment effect (Assumption 3). We generate $1000$ replications of several datasets determined by one of $72$ experimental scenarios. For each dataset, we find balancing weights from four different covariate balancing methods to estimate $\tau_{\text{ATE}}$. They are:

(IPW) Inverse probability of treatment weights where the propensity score follows a logit model fit using maximum likelihood estimation; 2. 2.

(CBPS) Inverse probability of treatment weights where the propensity score is fit to an exactly-specified logit model subject to (17). The propensity scores are fit using generalized method of moments as implemented in the CBPS package (Fong et al.,, 2019); 3. 3.

(SENT) Balancing weights that are estimated by minimizing the shifted relative entropy following the results of (27) and (28). Using these balancing weights instead of (19) and (20) allows us to test the effect of over-specifying the linear constraints when we know Assumption 3 is satisfied; 4. 4.

(BENT) Balancing weights that are estimated by minimizing the binary relative entropy subject to the constraints in (21) via the dual and primal solutions of (22) and (23).

We consider an extensive set of experimental scenarios adapted from those examined by Kang and Schafer, (2007). These scenarios vary the sample size $n\in\{200,1000\}$, the error variance $\sigma^{2}\in\{2,5,10\}$, the generative process that determines the treatment assignment (indexed by $\{a,b\}$), the outcome process (indexed by $\{a,b\}$), and the correlation between the potential outcomes, $\rho\in\{-0.3,0,0.5\}$. The covariates to be balanced (i.e the balance functions) are distributed as $X_{1},X_{2},X_{3},X_{4}\sim\mathcal{N}(0,1)$. Define the transformations $U_{1}=\exp(X_{1}/2)$, $U_{2}=X_{2}/[1+\exp(X_{1})]+10$, $U_{3}=(X_{1}X_{3}/25+0.6)^{3}$ and $U_{4}=(X_{2}+X_{4}+20)^{2}$. The vector $(U_{1},U_{2},U_{3},U_{4})^{T}$ is subsequently standardized to have a mean of zero and marginal variances of one.

The probability that a subject receives the treatment is then determined using the inverse logit link function,

[TABLE]

Scenarios $a$ and $b$ distinguish whether the log odds of the propensity score is either linear or non-linear with

[TABLE]

The treatment indicators are generated by sampling $Z_{i}\sim\text{Bin}\left(1,\pi^{(k)}_{i}\right)$. For the outcome process, we use the bivariate model

[TABLE]

where $\ell\in\{a,b\}$ indexes

[TABLE]

Once the potential outcomes have been generated, the observed outcome is the potential outcome corresponding to the observed treatment assignment. Each of the covariate balancing methods listed above are provided the design matrix with an intercept and the four original covariates; $X_{i1}$, $X_{i2}$, $X_{i3}$, and $X_{i4}$ for $i=1,2,\ldots,n$. The causal effect is then estimated using (3) where we substitute $p(\mathbf{X}_{i})$ with the balancing weights estimated by each method.

We found that the correlation between the potential outcomes did not affect the resulting estimates of $\tau_{\text{ATE}}$. In addition, the effects of altering $\sigma^{2}$ and $n$ had anticipated results. Lower values of $\sigma^{2}$ led to lower standard errors of the causal effect estimate whereas smaller values of $n$ led to larger standard errors. Therefore, we report the results for $\rho=0$, $n=200$, and $\sigma^{2}=10$ in Table 1 and Figure 2. The complete results appear in the online supplement.

For all the methods that we tested, if either the outcome model or the treatment assignment is correctly specified, then the causal effect estimate is unbiased. We see in Table 1 and Figure 2 that the balancing weights obtained with SENT perform as well, or better in some cases, than the exactly-specified CBPS model, even though the balancing weights obtained with SENT have twice as many constraints. The Monte Carlo standard error and bias of the estimates for $\tau_{\text{ATE}}$ are smallest when BENT is used to estimate the balancing weights for every scenario we examined. This is expected since these weights are used for estimating the $\tau_{\text{OWATE}}$ (Crump et al.,, 2006), and because every condition necessary to ensure that an estimator for $\tau_{\text{OWATE}}$ is also an estimator for $\tau_{\text{ATE}}$ are met. The Monte Carlo standard errors and mean square error are also uniformly smaller for the average treatment effect estimates when using balancing weights estimated by CBPS and SENT versus IPW. This result indicates that methods which exactly balance the empirical covariate distributions perform better in finite sample settings. When both the outcome and treatment assignment models are misspecified, the four methods for finding balancing weights all produced biased estimates of $\tau_{\text{ATE}}$. In these completely misspecified scenarios, the balancing weights estimated with BENT produce the least amount of bias and the lowest standard error for estimating $\tau_{\text{ATE}}$.

5.2 Heterogeneous Treatment Effect Simulation

In this section we simulate an additional $72$ scenarios with a linear conditional average treatment effect to test our proposed methods under Assumption 4. We use the same covariate distributions for $(X_{1},X_{2},X_{3},X_{4})^{T}$ and $(U_{1},U_{2},U_{3},U_{4})^{T}$ as in Section 5.1. We also recycle the conditional mean functions $\mu_{i}^{(\ell)}$ , $\ell\in\{a,b\}$, from (33). To generate the linear conditional average treatment effects, define

[TABLE]

For outcome scenarios $a$ and $b$, the bivariate outcome model is defined as

[TABLE]

from which we sample $n\in\{200,1000\}$ entries. Each unit’s treatment assignment is sampled from $\text{Bin}\left(1,\pi^{(k)}_{i}\right)$, where $\pi^{(k)}_{i}$ is determined by $\eta^{(k)}_{i}$, $k\in\{a,b\}$, which are defined in (32). Similar to the simulations conducted in Section 5.1, we also vary $\sigma^{2}\in\{2,5,10\}$ and $\rho\in\{-0.3,0,0.5\}$. For this set of scenarios we examine five different covariate balance methods:

(AIPW) Augmented inverse probability weights which uses the estimator

[TABLE]

$\hat{\mu}_{1}(\mathbf{X})$ is fit using linear regression on the treated group and $\hat{\mu}_{0}(\mathbf{X})$ is fit using linear regression on the controls. $\hat{\pi}(\mathbf{X})$ is fit with logistic regression; 2. 2.

(CAL) Calibration weights that solve (30) and (31) for $g(v)=\exp(-v)$, $v\in\Re$. This is equivalent to minimizing the unnormalized relative entropy subject to the linear constraints in (29). The R package ATE developed by Haris and Chan, (2015) is used to estimate these balancing weights; 3. 3.

(iCBPS) Inverse probability of treatment weights where the propensity scores follow an exactly-specified logit model which must satisfy (25) following the results of Fan et al., (2016). The propensity score is estimated using generalized method of moments as implemented in the CBPS package (Fong et al.,, 2019); 4. 4.

(hdCBPS) An augmented version of CBPS that extends (34) by using regularized regression techniques to debiased estimates of $\mu_{1}(\mathbf{X})$, and $\mu_{0}(\mathbf{X})$. The R package CBPS (Fong et al.,, 2019) is used to implement this method; 5. 5.

(SENT) Balancing weights that minimize the shifted relative entropy conditioned on the constraints in (26) where the estimates are obtained from (27) and (28).

For CAL, iCBPS, and SENT, the resulting balancing weights are substituted for $p(\mathbf{X})$ in (3) to estimate $\tau_{\text{ATE}}$. The augmented approach of AIPW was first proposed by Robins et al., (1994) while hdCBPS uses the augmented estimator proposed by Ning et al., (2018).

As with the previous simulation study, it appears that the correlation between potential outcomes is inconsequential while the Monte Carlo standard errors predictably decrease when either $n$ increases or $\sigma^{2}$ decreases. A representative selection of results from the experiment where $\rho=0$, $n=200$, and $\sigma^{2}=10$ are found in Table 2 and Figure 3. The complete results can be found in the online supplement. This simulation demonstrates that all five methods enjoy the doubly-robust property described in Theorem 2. For each scenario, CAL and SENT had similar levels of bias and variation despite using different criterion distance functions. Even though the constraints of iCBPS and SENT are the same, the differences between the optimization techniques of the two methods becomes quite apparent. The Monte Carlo standard error of the estimates for $\tau_{\text{ATE}}$ using balancing weights obtained with SENT is smaller than the standard error of the estimates using balancing weights found with iCBPS. AIPW, CAL, and SENT performed similarly whenever the outcome was correctly specified. However, when the outcome model is misspecified and the propensity score is correctly specified, the Monte Carlo standard error and mean squared error were greater with AIPW compared to SENT and CAL. This suggests, and is further confirmed by hdCBPS, that methods which exactly-balance covariate distributions can improve the efficiency of a doubly-robust estimator in finite samples. hdCBPS performed about as well as SENT and CAL. However, this method was proposed to alleviate issues encountered with high-dimensional covariate data, rendering many of its benefits redundant in this low-dimensional simulation study.

5.3 Illustrative Example of Unplanned Readmissions after Lung Resection

Next, we investigate the results of a real data set using different weighting and matching methods. In Bhagat et al., (2017), the odds of unplanned, 30-day readmissions are compared between lung cancer patients that receive thoracoscopic versus open lung resections. The study identified 9,510 patients that underwent some form of lung resection from the American College of Surgeons - National Safety and Quality Innovation Program (ACS-NSQIP) database. Of those 9,510 patients, 4,935 (51.9%) received a thoracoscopic resection and 4,575 (48.1%) received an open anatomic resection. The study analysis carried out a greedy one-to-one matching of patients using the estimated propensity score as the criterion matching function (Ho et al.,, 2007). The propensity scores were fit with standard logistic regression. This algorithm matched 3,399 thoracoscopic lung resection patients to 3,399 open anatomic lung resection patients, dropping 2,712 patients (28.5%). In doing so, the “treated” group are assumed to be the patients that receive thoracoscopic lung resections and represent a random sample of the target population. Thus the casual effect being estimated is the average treatment effect of the treated.

We replicated the study conducted in Bhagat et al., (2017) by estimating balancing weights using two different methods. The first method uses entropy balancing (EB) where the estimated balancing weights are obtained with (14) and (15). Recall that these balancing weights applied to (3) is doubly-robust (Zhao and Percival,, 2017). The second method fits a propensity score model using logistic regression (IPW). With the fitted propensity score, we then use (2) substituting $\hat{\pi}(\mathbf{X})$ for $\pi(\mathbf{X})$ to estimate $\tau_{\text{ATT}}$. The causal effect estimates using the propensity score matched (PSM) cohort from the original paper are also reported along with the unadjusted (UN) results in Table 3.

Figure 4 shows the amount of imbalance observed for each of the covariates among those included in the covariate balancing models. We see that across each covariate, entropy balancing perfectly balances the first sample moments of the covariate distribution between the two treatment groups. Logistic regression appears to be less adequate at balancing the covariate moments than matching. However, aside from hospital length of stay, each of the weighted mean differences fell within the conservative 0.05 unit threshold using the inverse probability of treatment weights. The unadjusted differences do not share the same success as their adjusted counterparts, suggesting that some form of balancing should be implemented. After estimating $\tau_{\text{ATT}}$, notice in Table 3 that the estimated risk difference is significant when using either the inverse probability of treatment weights or entropy balancing, but is not significant when using propensity score matching. The difference is likely due to the 2,712 patients that were omitted when matching. This discrepancy illuminates and emphasizes the importance of selecting the most appropriate method for balancing covariate data, even within a large observational study.

6 Discussion

The generalized projection of a Bregman distance from a vector of sampling weights onto a set of intersecting hyperplanes is a powerful and flexible tool for normalizing data. In particular, this process is quite useful for constructing balancing weights for estimating causal effects. Using properties of dual optimization, we identify a doubly-robust estimator for $\tau_{\text{ATE}}$ and the optimally weighted average treatment effect (Crump et al.,, 2006) in Theorem 1 and Corollary 1. We also show that the dual interpretation of improved CBPS (Fan et al.,, 2016) is doubly-robust in Theorem 2. In Theorem 3, we present the conditions for which the balancing weights produced by (11) are the same as the weights produced by (31) suggested by Chan et al., (2015). When the true balance functions are unknown, we can use nonparametric methods similar to those suggested by Hirano et al., (2003) and Chan et al., (2015) within our framework to achieve global efficiency.

In the simulation studies we conducted, we observed that the balancing weights that are typically used to estimate $\tau_{\text{OWATE}}$ had the best performance for estimating $\tau_{\text{ATE}}$ when Assumption 3 is satisfied. We also observed that including additional constraints as in (8) sometimes had better performance than estimators that require fewer constraints, like CBPS. When we assume a linear conditional average treatment effect, our dual interpretation of iCBPS performed better than the analogous GMM estimator. We then apply our framework to a real data set of lung resection patients. Here we demonstrate how the choice of balancing method can have a critical impact on the results of a study.

There are several limitations to our proposed framework. First, each sampling unit’s treatment assignment is assumed to be independent from the treatment assignment of the other sampling units. This assumption is sometimes called the no interference assumption. Health outcomes research is rich in observational data from the emergence of the electronic health record. While numerous in size, these datasets are more complex with patients being clustered within regions, hospitals, clinics, and/or practicing physicians. These are all factors that need to be accounted for in some way. How to extend these methods to clustered data settings is currently under investigation. Second, linear equality constraints are often quite stringent. If a particular covariate is difficult to balance, our proposed framework will sometimes fail to find the appropriate balancing weights. Zubizarreta, (2015) proposes using stable balancing weights which places linear inequality constraints on the weighted sample moments of the covariate distribution while minimizing the Euclidean distance. In more recent work, Wang and Zubizarreta, (2019) have combined this interval constrained optimization approach with calibration estimators. There is also the issue where the balance functions that generate either the outcome or treatment assignment are high-dimensional. This problem is not examined in the presented work. Ning et al., (2018) propose using an augmented approach with the covariate balance propensity scores of Imai and Ratkovic, (2013) and Fan et al., (2016) in the spirit of Robins et al., (1994) and Farrell, (2015). Their proposed methodology boasts compelling results as the dimension of covariate distribution increases. It is possible that our approach could be extended using some of the methods proposed by Ning et al., (2018).

In addition to addressing some of the limitations identified in the previous paragraph, in future work we would also like to expand these methods to incorporate multivalued treatment assignments. This would entail modifying the Horvitz-Thompson estimator and also requires extending the constraint matrix $\mathbf{A}$ and target margins $\mathbf{b}$ to facilitate covariate balance between all pairwise combinations of the treatment assignments. Finally, we would like to further investigate methods for generalizing causal effect estimates to a target population, which would involve estimating $\mathbf{q}$ prior to estimating $\mathbf{p}$.

Acknowledgments

Special Thanks: We would like to thank Dr. Robert Meguid and the Adult and Child Consortium for Health Outcomes Research and Delivery Science (ACCORDS) program for making the unplanned readmissions data available to us. We would also like to thank Dr. Peter DeWitt for his help in deploying the associated R package cbal.

Funding information: This research is supported by a pilot grant from the Data Science to Patient Value (D2V) initiative from the University of Colorado. We acknowledge partial support from NSF award number DMS-1420451.

Disclaimer: This manuscript will be submitted to the Department of Biostatistics and Informatics in the Colorado School of Public Health, University of Colorado Denver, in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Biostatistics for Kevin P. Josey.

Appendix A R Package and Simulation Code

The R package used to fit balancing weights as the generalized projection of Bregman distance is still in development with a working version available at https://github.com/kevjosey/cbal. The code used to conduct the simulation study in Section 5 is available at the following URL: https://github.com/kevjosey/cbal-sim. The code for replicating the study by Bhagat et al., (2017) is available from the authors upon request.

Appendix B Technical Proofs

B.1 Proof of Lemma 1

Proof.

For any $\mathbf{w}\in\Omega\cap\bar{\Delta}^{n}$, define

[TABLE]

The set $\mathcal{S}$ is bounded due to the requirement that Bregman distances have bounded partial level sets (Definition 2.1.1 of Censor and Zenios, (1998)) and closed since $D_{f}(\mathbf{p}\;\|\;\mathbf{q})$ is continuous in $\mathbf{p}\in\bar{\Delta}^{n}$. Therefore, the non-empty intersection $\mathcal{T}\equiv\Omega\cap\bar{\Delta}^{n}\cap\mathcal{S}$ is bounded. Since $\Omega\cap\bar{\Delta}^{n}$ and $\mathcal{S}$ are closed, $\mathcal{T}$ is also closed and hence compact. Thus, $D_{f}(\mathbf{p}\;\|\;\mathbf{q})$, takes its infimum subject to $\mathbf{p}\in\Omega\cap\bar{\Delta}^{n}$ at $\hat{\mathbf{p}}\in\mathcal{T}$.

To prove that $\hat{\mathbf{p}}\in\Omega\cap\bar{\Delta}^{n}$ is unique, suppose there are two points $\hat{\mathbf{p}},\tilde{\mathbf{p}}\in\Omega\cap\bar{\Delta}^{n}$ such that $\hat{\mathbf{p}}\neq\tilde{\mathbf{p}}$ and $D_{f}(\hat{\mathbf{p}}\;\|\;\mathbf{q})=D_{f}(\tilde{\mathbf{p}}\;\|\;\mathbf{q})=\operatorname*{arg\,min}_{\mathbf{p}\in\Omega\cap\bar{\Delta}^{n}}D_{f}(\mathbf{p}\;\|\;\mathbf{q})$, for some $\mathbf{q}\in\Delta^{n}$. By the convexity of $\Omega\cap\bar{\Delta}^{n}$, $(\hat{\mathbf{p}}+\tilde{\mathbf{p}})/2\in\Omega\cap\bar{\Delta}^{n}$. By the strict convexity of $f$, we have

[TABLE]

which is a contradiction. ∎

B.2 Proof of Proposition 1

Proof.

Assume $\hat{\mathbf{p}}=\operatorname*{arg\,min}_{\mathbf{p}\in\Omega\cap\bar{\Delta}^{n}}D_{f}(\mathbf{p}\;\|\;\mathbf{q})$ exists. By Lemma 1, $\hat{\mathbf{p}}$ is a global minimum subject to $\mathbf{p}\in\Omega$. Due to the zone consistency assumption, if $\mathbf{q}\in\Delta^{n}$ then the minimum must be contained in $\Delta^{n}$ and not $\bar{\Delta}^{n}$. This assumption allows us to apply the Lagrange multiplier theorem (Bertsekas,, 1999) which proves the existence of a unique $\hat{\boldsymbol{\lambda}}\in\Re^{m}$ that satisfies

[TABLE]

Note that due to the strict convexity of $f$, $\nabla f:\Delta^{n}\rightarrow\Re^{n}$ is strictly increasing and is therefore injective. Thus, the transformation $(\nabla f)^{-1}:\Re^{n}\rightarrow\Delta^{n}$ also has a unique mapping. If we solve (35) for $\hat{\mathbf{p}}$, we get

[TABLE]

∎

B.3 Proof of Proposition 2

Proof.

Let $\boldsymbol{\lambda}\in\Re^{m}$ and $\hat{\mathbf{p}}=\operatorname*{arg\,min}_{\mathbf{p}\in\bar{\Delta}^{n}}L(\mathbf{p},\mathbf{q},\boldsymbol{\lambda})$. Since $f$ is zone consistent, a minimum of the Lagrangian with respect to $\mathbf{p}\in\Delta^{n}$ must satisfy

[TABLE]

By the strict convexity of $f$ (see the proof in Section B.2), we can solve for $\hat{\mathbf{p}}$ to get

[TABLE]

From (36) we know that for all $\mathbf{p}\in\Omega\cap\Delta^{n}$,

[TABLE]

According to the definition of the generalized projection, minimizing the right hand side subject to $\mathbf{p}\in\Omega\cap\Delta^{n}$ yields the solution to the primal problem and

[TABLE]

for all $\boldsymbol{\lambda}\in\Re^{m}$. By Proposition 1, there exists some $\hat{\boldsymbol{\lambda}}\in\Re^{m}$ such that $\hat{\mathbf{p}}=P_{f}\left(\mathbf{q},\mathbf{A}\hat{\boldsymbol{\lambda}}\right)$. Therefore,

[TABLE]

Substituting the result of (38) into (37), we have

[TABLE]

for all $\boldsymbol{\lambda}\in\Re^{m}$ which implies

[TABLE]

∎

B.4 Proof of Theorem 1

Proof.

Denote $\tau^{*}=\mathbb{E}[Y(1)-Y(0)]$ and suppose

[TABLE]

where $\beta^{*}_{j}\in\Re$, $j=1,2,\ldots,m$, denotes the true coefficient values of the outcome model. Let $\hat{p}(\mathbf{X}_{i})$ be determined by (20) where $\hat{\boldsymbol{\lambda}}$ is solved with (19).

Recall that another doubly-robust estimator for the average treatment effect is

[TABLE]

where

[TABLE]

so that $\hat{\mu}_{0}(\mathbf{X}_{i})\rightarrow_{p}\mu^{*}_{0}(\mathbf{X}_{i})$ and $\hat{\mu}_{1}(\mathbf{X}_{i})\rightarrow_{p}\mu^{*}_{1}(\mathbf{X}_{i})$. If we substitute $\hat{p}(\mathbf{X}_{i})$ for $\hat{\pi}(\mathbf{X}_{i})^{-1}$ and $[1-\hat{\pi}(\mathbf{X}_{i})]^{-1}$ and assume, without loss of generality, that $\sum_{i=1}^{n}Z_{i}\hat{p}(\mathbf{X}_{i})=n$, we get

[TABLE]

This shows that the balancing weight method is algebraically equivalent to the augmented regression approach of Robins et al., (1994) under the condition that the outcome model is correctly specified. A similar proof was applied to entropy balancing for estimating $\tau_{\text{ATT}}$ in Zhao and Percival, (2017).

Now, let $\text{logit}[\pi^{*}(\mathbf{X})]=\text{logit}[\pi(\mathbf{X};\boldsymbol{\lambda}^{*})]=\sum_{j=1}^{m}c_{j}(\mathbf{X})\lambda^{*}_{j}$ and suppose $\mu_{1}(\mathbf{X})$ and $\mu_{0}(\mathbf{X})$ are unknown. Denote the $(m\times 1)$ vector of the linearly independent balance functions as

[TABLE]

To show consistency under these conditions, we follow the approach of Tsiatis, (2006) and Kennedy, (2016) by employing influence functions about the parameters. Define the estimating equation for $\boldsymbol{\lambda}$ as

[TABLE]

which is the first order condition of (19). As a result of Propositions 1 and 2, $\hat{\boldsymbol{\lambda}}$ must satisfy

[TABLE]

Next, define the estimating equation for $\tau$ as

[TABLE]

The estimating equation defined in (42) is solved with the Horvitz-Thompson estimator which satisfies

[TABLE]

Define $\boldsymbol{\theta}\equiv(\boldsymbol{\lambda}^{T},\tau)^{T}$ and the the stacked estimating equations as

[TABLE]

Under standard regularity assumptions described in Tsiatis, (2006), it can be shown using Taylor’s theorem that

[TABLE]

It is then straightforward to identify the influence function for $\tau$ from (43), which we denote as

[TABLE]

This allows us to express

[TABLE]

When $\text{logit}[\pi^{*}(\mathbf{X})]=\sum_{j=1}^{m}c_{j}(\mathbf{X})\lambda^{*}_{j}$, then $\mathbb{E}\left[\boldsymbol{\zeta}(\mathbf{X},Z;\boldsymbol{\lambda}^{*})\right]=\mathbf{0}_{m}$. In a similar fashion, and given Assumption 2, it is trivial to show $\mathbb{E}\left[\psi(\mathbf{X},Y,Z;\boldsymbol{\lambda}^{*},\tau^{*})\right]=0$. After applying the weak law of large numbers to (45), we conclude $\hat{\tau}\rightarrow_{p}\tau^{*}$. Notice that Assumption 3 is not required to prove this result.

For simplicity, denote $\phi^{*}(\mathbf{X},Y,Z)=\phi(\mathbf{X},Y,Z;\boldsymbol{\lambda}^{*},\boldsymbol{\tau}^{*})$. Following the classical central limit theorem, the influence function for $\tau$ can also be used to derive the asymptotic distribution which is

[TABLE]

To simplify the algebra further, we introduce a few nuisance parameters and rewrite (41) and (42). Define $\boldsymbol{\gamma}\equiv\mathbb{E}[\mathbf{c}(\mathbf{X})]$ and $\mu_{0}\equiv\mathbb{E}[Y(0)]=\mathbb{E}[\mu_{0}(\mathbf{X})]$. If $\mathbf{c}(\mathbf{X})$ includes an intercept term (i.e. $c_{1}(\mathbf{X})=1$), then by Propositions 1 and 2 we can equivalently express (41) as

[TABLE]

for any $\boldsymbol{\gamma}\in\Re^{m}$. Equation (42) can be rewritten as

[TABLE]

for any $\mu_{0}\in\Re$. Moreover, let $\mathbf{H}_{\mathbf{U},\mathbf{V}}(W)=\mathbb{E}\{W[\mathbf{U}-\mathbb{E}(\mathbf{U})][\mathbf{V}-\mathbb{E}(\mathbf{V})]^{T}\}$, $\mathbf{H}_{\mathbf{U}}(W)=\mathbf{H}_{\mathbf{U},\mathbf{U}}(W)$ and

[TABLE]

Given Assumption 2, the asymptotic variance is expanded into

[TABLE]

If we suppose $\mu_{1}(\mathbf{X})=\mu^{*}_{1}(\mathbf{X})$ and $\mu_{0}(\mathbf{X})=\mu^{*}_{0}(\mathbf{X})$ as in (39), we can write

[TABLE]

This allows us to reduce (46) into the more succinct format of

[TABLE]

where the second equality holds due to the conditional variance identity. ∎

Remark 1**.**

We avoid using sampling weights in Theorem 1 due to the observation that

[TABLE]

for some $q\in(1,\infty)$. To resolve this issue, we define $q(1)$ and $q(0)$ such that $q(1)-1=[q(0)-1]^{-1}$, which can be used to identify the sampling weights as $q\equiv Zq(1)+(1-Z)q(0)$. It is straightforward to show

[TABLE]

We would then assume $\pi^{*}(\mathbf{X})^{-1}=1+[q(1)-1]\exp\left[-\sum_{j=1}^{n}c_{j}(\mathbf{X})\lambda_{j}\right]$ in condition 1 of Theorem 1 and substitute this notation accordingly into the proof.

B.5 Proof of Corollary 1

Proof.

Since the constraints placed on the weighting estimator in Theorem 1 are identical to the constraints in Corollary 1, the consistency of the Horvitz-Thompson estimator under Assumption 3 and the condition that $\mu_{0}(\mathbf{X})$ is correctly specified is proven in a nearly identical manner to Theorem 1 using balancing weights that are instead determined by (22) and (23). The only difference would be to use an augmented regression estimator for $\tau_{\text{OWATE}}$ instead of $\tau_{\text{ATE}}$.

Assume $\text{logit}[\pi^{*}(\mathbf{X})]=\text{logit}[\pi(\mathbf{X};\boldsymbol{\lambda}^{*})]=\sum_{j=1}^{m}c_{j}(\mathbf{X})\lambda^{*}_{j}$ and let

[TABLE]

Our approach to proving consistency under the condition that the treatment assignment is correctly specified is similar to the proof in Section B.4. However, we instead define

[TABLE]

which is the first order condition of (23) including a normalizing constant in the denominator (Hirano et al.,, 2003), and

[TABLE]

From the proof of Theorem 1, the influence function for $\tau$ is identified as

[TABLE]

When $\text{logit}[\pi^{*}(\mathbf{X})]=\sum_{j=1}^{m}c_{j}(\mathbf{X})\lambda^{*}_{j}$ and given Assumption 2, then

[TABLE]

so that $\mathbb{E}\left[\psi^{\prime}(\mathbf{X},Y,Z;\boldsymbol{\lambda}^{*},\tau^{*})\right]=0$. In a similar fashion, we can show $\mathbb{E}\left[\boldsymbol{\zeta}^{\prime}(\mathbf{X},Z;\boldsymbol{\lambda}^{*})\right]=\mathbf{0}_{m}$. We then apply the weak law of large numbers to conclude $\hat{\tau}\rightarrow_{p}\tau^{*}$.

Remark 2**.**

If we set $\tau^{*}=\mathbb{E}[Y(1)-Y(0)]=\mathbb{E}[\mu_{1}(\mathbf{X})-\mu_{0}(\mathbf{X})]$, then we would require Assumption 3 to hold in order for $\mathbb{E}\left[\psi^{\prime}(\mathbf{X},Y,Z;\boldsymbol{\lambda}^{*},\tau^{*})\right]=0$.

Recall from the proof of Theorem 1 that

[TABLE]

Notice that

[TABLE]

where $\psi(\mathbf{X},Y,Z;\boldsymbol{\lambda},\tau)$ is defined in (42) and

[TABLE]

where $\boldsymbol{\zeta}(\mathbf{X},Z;\boldsymbol{\lambda})$ is defined in (41). If $\mu^{*}_{1}(\mathbf{X})$ and $\mu^{*}_{0}(\mathbf{X})$ follow the form of (39), then

[TABLE]

This means

[TABLE]

where $\phi(\mathbf{X},Y,Z;\boldsymbol{\lambda},\tau)$ is defined in (44). According to the proof of Theorem 1, we know

[TABLE]

which allows us to write

[TABLE]

∎

Remark 3**.**

The proof of Theorem 1, and by extension Corollary 1, exploits the identity

[TABLE]

where $\omega^{*}(\mathbf{X})$ is a known weighting function. In the case of Theorem 1, $\omega^{*}(\mathbf{X})=1$, while for Corollary 1, $\omega^{*}(\mathbf{X})=\pi^{*}(\mathbf{X})[1-\pi^{*}(\mathbf{X})]$. In the case where we find estimators for $\tau_{\text{ATT}}$, we require

[TABLE]

The identity in (47) requires careful selection of the generalized projection, which should complement the assumed functional form of the propensity score. To generalize Theorem 1 would require modifying the statement of condition 1 to assume $\Pr\{Z=1|\mathbf{X}\}=P_{f}\left[q,\sum_{j=1}^{n}c_{j}(\mathbf{X})\lambda_{j}\right]^{-1}$ and $\Pr\{Z=0|\mathbf{X}\}=P_{f}\left[q,-\sum_{j=1}^{n}c_{j}(\mathbf{X})\lambda_{j}\right]^{-1}$. If this condition is satisfactory, the proof of Theorem 1 would be the same albeit with changes to the notation. However, we found this condition to be counter-intuitive in instances where we violate (47).

B.6 Proof of Theorem 2

Proof.

First, we will prove that $\hat{\tau}$ is consistent for $\tau_{\text{ATE}}$ when $\mu_{0}(\mathbf{X})$ lies on the span of a linearly independent set of balance functions, $\{c_{j}(\mathbf{X}):j=1,2,\ldots,m\}$. Let $\hat{p}(\mathbf{X}_{i})$ be determined by (28) where $\hat{\boldsymbol{\lambda}}$ is solved by (27). Assume that the conditional means of the potential outcomes are

[TABLE]

where $\alpha^{*}_{j}$ and $\beta^{*}_{j}$ are the true values of the regression coefficients, and $\tau^{*}=\mathbb{E}\left[\mu^{*}_{1}(\mathbf{X})-\mu^{*}_{0}(\mathbf{X})\right]$. Notice that if $c_{1}(\mathbf{X}_{i})=1$ for all $i=1,2,\ldots,n$, then $n=\sum_{i=1}^{n}\hat{p}(\mathbf{X}_{i})Z_{i}$. We will again show the algebraic equivalency of the exact balancing weight estimator with the doubly-robust estimator

[TABLE]

where

[TABLE]

so that $\hat{\mu}_{0}(\mathbf{X}_{i})\rightarrow_{p}\mu^{*}_{0}(\mathbf{X}_{i})$ and $\hat{\mu}_{1}(\mathbf{X}_{i})\rightarrow_{p}\mu^{*}_{1}(\mathbf{X}_{i})$. We then substitute $\hat{p}(\mathbf{X}_{i})$ for $\hat{\pi}(\mathbf{X}_{i})^{-1}$ and $[1-\hat{\pi}(\mathbf{X}_{i})]^{-1}$ depending to whether $Z_{i}=1$ and $Z_{i}=0$, respectively. This gives us

[TABLE]

Now we let $\text{logit}[\pi^{*}(\mathbf{X})]=\text{logit}[\pi(\mathbf{X};\boldsymbol{\lambda}^{*}_{0})]=\sum_{j=1}^{m}c_{j}(\mathbf{X})\lambda^{*}_{j0}$ and assume $\mu_{1}(\mathbf{X})$ and $\mu_{0}(\mathbf{X})$ are unknown. We will reuse much of the notation and definitions presented in the proof of Theorem 1. First, we redefine the estimating equations for $\boldsymbol{\lambda}$ as

[TABLE]

with $\boldsymbol{\zeta}(\mathbf{X},Z;\boldsymbol{\lambda})=\left[\boldsymbol{\zeta}_{0}(\mathbf{X},Z;\boldsymbol{\lambda})^{T},\boldsymbol{\zeta}_{1}(\mathbf{X},Z;\boldsymbol{\lambda})^{T}\right]^{T}$. We also redefine the estimating equation for $\tau$ as

[TABLE]

Recall from the proof of Theorem 1 that under the standard regularity assumptions of Tsiatis, (2006), we have

[TABLE]

Next, observe that

[TABLE]

The previously mentioned regularity assumptions allow us to interchange differentiation and expectation to obtain the first derivative of (50) with respect to $\boldsymbol{\lambda}_{1}$, which yields a negative-definite matrix (given that the balance functions are linearly independent). This implies (50) can only be zero if $\lambda^{*}_{j1}=0$ for all $j=1,2,\ldots,m$. It is then trivial to show $\mathbb{E}[\boldsymbol{\zeta}_{0}(\mathbf{X},Z;\boldsymbol{\lambda}^{*})]=0$ and $\mathbb{E}[\psi(\mathbf{X},Y,Z;\boldsymbol{\lambda}^{*},\tau^{*})]=0$. This means that the expected value of (49) is zero, and we can apply the weak law of large numbers to arrive at our desired result where $\hat{\tau}\rightarrow_{p}\tau^{*}$.

Due to the central limit theorem, we know $\sqrt{n}(\hat{\tau}-\tau^{*})\rightarrow_{d}\mathcal{N}\left\{0,\mathbb{E}\left[\phi^{*}(\mathbf{X},Y,Z)^{2}\right]\right\}$ where $\phi^{*}(\mathbf{X},Y,Z)=\phi(\mathbf{X},Y,Z;\boldsymbol{\lambda}^{*},\tau^{*})$ is defined in (44). Without loss of generality (and by the assumption that $c_{1}(\mathbf{X})=1$) redefine the estimating equations

[TABLE]

for some $\boldsymbol{\gamma}\in\Re^{m}$, and

[TABLE]

for some $\mu_{0}\in\Re$. To lighten the notation further, we will write $\psi^{*}\equiv\psi(\mathbf{X},Y,Z;\boldsymbol{\lambda}^{*},\tau^{*})$ and $\boldsymbol{\zeta}^{*}\equiv\boldsymbol{\zeta}(\mathbf{X},Z;\boldsymbol{\lambda}^{*})$. We then expand the asymptotic variance into more digestible components with

[TABLE]

Suppose $\mu_{1}(\mathbf{X})=\mu^{*}_{1}(\mathbf{X})$ and $\mu_{0}(\mathbf{X})=\mu^{*}_{0}(\mathbf{X})$ which are defined in (48). Then

[TABLE]

and, with some algebra,

[TABLE]

This allows us to further reduce (51) into

[TABLE]

where the last equality holds due to the conditional variance identity and by recognizing

[TABLE]

∎

B.7 Proof of Theorem 3

Proof.

Let $\mathbf{p}\in\Omega_{1}$ which satisfies the constraints

[TABLE]

This implies $\sum_{i=1}^{n}p_{i}(1-Z_{i})c_{j}(\mathbf{X}_{i})=b_{j}$ so that $\mathbf{p}\in\Omega_{0}$. Thus $\Omega_{1}=\Omega_{0}$.

According to Lemma 1, the generalized projection is unique implying

[TABLE]

Because $G$ is a monotone increasing transformation of $D_{f}$, the optimal point is preserved and we conclude

[TABLE]

∎

Tables and Figures