2.1 Intuition: targeting weights rather than propensity scores

To gain intuition, the left panel of Figure 2 presents the true (the x-axis) and estimated weights (the y-axis) for each observation. The proposed method (triangles) estimates inverse probability weights better than the MLE (circles). Most of the triangles are within or near the shaded narrow corridor around the diagonal line, which indicates that the difference between the true and estimated weights is less than two, whereas many circles lie far from the narrow corridor.

Figure 2. A numerical example: The true and estimated weights and propensity scores with a misspecified propensity score model. Notes: This figure compares the performance of the proposed estimator (triangles) and MLE estimator (circles) with a misspecified propensity score model in terms of inverse probability weights (left panel) and propensity scores (right panel). In the left (right) panel, the x-axis represents the true weights (propensity scores), and the y-axis represents the estimated ones. The shaded area in the left panel indicates that the difference between the true and estimated weights is small (less than two) whereas the shaded area in the right panel indicates the corresponding area where the difference in weights is small (not the difference in propensity scores).

The right panel shows the true (the x-axis) and estimated propensity scores (the y-axis) for each observation. Overall, the MLE (circles) estimates propensity scores better than the proposed method (triangles) but the proposed method restricts the difference more effectively than the MLE when the difference in the resulting inverse probability weights is large shown as the shaded area. This explains why the proposed method approximates the weights better than the MLE.

These results provide intuition when the difference between the two estimators is large: It is when some units have large estimated weights, which is likely to happen when there are unmodeled nonlinear effects of covariates on treatment assignment, there are unmodeled interaction effects of covariates on treatment assignment, or the numbers of the treated and controlled units are unbalanced. In addition, the proposed estimator is especially better than the MLE when the outcomes of these units are extreme (relative to the true average outcome), but not so much when the opposite is true. The difference in the estimated weights between the two estimators is small when no units have large estimated weights.

The root mean squared errors (RMSEs) of the estimated weights are 1.18 for the proposed method and 3.10 for the MLE, and RMSEs of the estimated propensity scores are 0.198 for the proposed method and 0.131 for the MLE, indicating that the proposed method estimates the weights better than the MLE whereas the MLE estimates the propensity scores better than the proposed method. To estimate the quantity of interest through estimating inverse probability weights, the proposed method is better suited than the MLE.

3.1 Setup

Suppose that there is a random sample of n units (i = 1,  2,  …,  n) from a population. Each unit consists of a triplet (Y _i,  R _i,  X _i), where Y _i is an outcome variable, R _i is a binary outcome response indicator variable that takes one when Y _i is observed and zero otherwise, and X _i is a vector of observed covariates whose support is denoted by ${\cal X}$. The propensity score is defined as the conditional response probability given observed covariates and denoted as $\Pr ( R_i = 1 \mid X_i) = \pi ( X_i)$. We impose the overlap assumption that the propensity score is bounded away from 0 and 1: 0 < π(x) < 1 for any ${\bf x} \in {\cal X}$. We also assume that the outcome is missing at random, i,e., the outcome variable does not account for the response conditional on the observed covariates: $Y_i \; \mathop {\perp \!\!\!\!\!\!\perp }\; R_i \mid X_i$ (Little and Rubin, Reference Little and Rubin2019). These assumptions correspond to the conditional ignorability assumption in causal inference.

For simplicity, the parameter of interest in this study is the average outcome $\mu = {\mathbb E}[ Y_i]$ for the target group of R _i = 1 and R _i = 0, and the proposed method is easily extended to other estimands such as the average treatment effect in causal inference (Kang and Schafer, Reference Kang and Schafer2007; Vermeulen and Vansteelandt, Reference Vermeulen and Vansteelandt2015; Zubizarreta, Reference Zubizarreta2015; Wong and Chan, Reference Wong and Chan2017; Seaman and Vansteelandt, Reference Seaman and Vansteelandt2018; Zhao, Reference Zhao2019; Tan, Reference Tan2020; Wang and Zubizarreta, Reference Wang and Zubizarreta2020).

3.2 Problems in the inverse probability weighting with the maximum likelihood estimation

The average outcome can be estimated unbiasedly and consistently by the IPW with true propensity scores: $\sum _{i = 1}^n R_i w_i Y_i = \sum _{i = 1}^n R_i Y_i / \pi ( X_i)$, where w _i = 1/π(X _i) is the true inverse probability weight for unit i. In observational studies, however, neither the true propensity scores nor the true inverse probability weights are known. In practice, researchers use a parametric propensity score model π(X _i,  β), where the most common choice is the logistic regression model $\pi ( X_i,\; \, \beta ) = \exp ( X_i^{\sf T} \beta ) / ( 1 + \exp ( X_i^{\sf T} \beta ) )$. Typically, the nuisance parameters β are estimated using the MLE to obtain the estimated propensity scores for each unit $\hat {\pi }_{\rm MLE}( X_i) = \pi ( X_i,\; \, \hat {\beta }_{\rm MLE})$. The resulting IPW estimator with the estimated propensity scores is consistent if the propensity score model is correct.

However, if the propensity score model is misspecified, the MLE estimator for β and the resulting IPW estimator is not consistent. Generally, the limiting values of $\hat {\beta }_{\rm MLE}$ are the solution to the following problem (Akaike, Reference Akaike1973; Tan, Reference Tan2020; White, Reference White1982):

(1)$$\hat{\beta}_{\rm MLE} = \mathop{\rm arg max}\limits_{\beta} \; \sum_{i = 1}^n \{ R_i \log( \pi( X_i,\; \beta) ) + ( 1 - R_i) \log( 1 - \pi( X_i,\; \beta) ) \} $$

(2)$$ = \mathop{\rm arg min}\limits_{\beta}\; \{ {\rm KL}( R,\; \pi( {\bf X},\; \beta) ) + {\rm KL}( 1 - R,\; 1 - \pi( {\bf X},\; \beta) ) \} .$$

This implies that the MLE estimates $\hat {\beta }_{\rm MLE}$ to minimizes the (generalized) KL divergence between the true responses (R and 1 − R) and their estimates ($\pi ( {\bf X},\; \, \hat {\beta }_{\rm MLE})$ and $1 - \pi ( {\bf X},\; \, \hat {\beta }_{\rm MLE})$) for the target group of R _i = 1 and R _i = 0, where the (generalized) KL divergence is the most common measure of the difference between two points P and Q, which is defined as ${\rm KL}( P,\; \, Q) = \sum _{i = 1}^{n} \{ P_i \log ( P_i / Q_i) - P_i + Q_i\}$. The problem is that this is not the same as the KL divergence between the true weights and estimated inverse probability weights although the latter determines the performance of the IPW estimator.

3.3 Distribution balancing weighting

This study proposes the distribution balancing weighting (DBW), which estimates such propensity scores $\pi ( X_i,\; \, \hat {\beta }_{\rm DBW})$ that minimize the KL divergence between the true and the estimated inverse probability weights. The idea of minimizing the KL divergence with possibly misspecified models originates from the quasi MLE, which minimizes the KL divergence between the true and estimated distribution when the true distribution is unknown (Akaike, Reference Akaike1973; White, Reference White1982). Since the propensity scores and their inverse probability weights have a non-linear relationship, minimizing the KL divergence between the true and estimated response assignment like the MLE does not minimize that for the inverse probability weights in general. The DBW applies the idea to the KL divergence between the true and estimated weights and obtains the following KL divergence:

(3)$${\mathbb E} \left[{\rm KL}\left({R\over \pi( {\bf X}) },\; {R\over \pi( {\bf X},\; \hat{\beta}_{\rm DBW}) } \right)\right] = {\mathbb E} \left[\sum_{i = 1}^n \left\{{R_i\over \pi( X_i,\; \hat{\beta}_{\rm DBW}) } + \log( \pi( X_i,\; \hat{\beta}_{\rm DBW}) ) - 1 - \log( \pi( X_i) ) \right\}\right].$$

The DBW estimates coefficients as the minimizers of the following loss function, which is the sample version of the KL divergence above up to a constant:

(4)$$\hat{\beta}_{\rm DBW} = \mathop{\rm arg min}\limits_{\beta}\; \sum_{i = 1}^n \left\{{R_i\over \pi( X_i,\; \beta) } + \log( \pi( X_i,\; \beta) ) \right\}.$$

As stated formally in proposition 5, the DBW estimates propensity scores and resulting inverse probability weights consistently, like the MLE, when the propensity score model is correctly specified.

3.4 Minimizing an upper bound of bias

The main advantage of the DBW is that it also minimizes the imbalance in the multivariate covariate distribution in the target and weighted groups as stated in proposition 1 below, whose importance has been acknowledged in the existing studies (Hainmueller, Reference Hainmueller2012; Zubizarreta, Reference Zubizarreta2015; Li et al., Reference Li, Morgan and Zaslavsky2018; Zhao, Reference Zhao2019). While the propensity scores estimated by the MLE minimize the distributional imbalance in the limit when the propensity score model is correctly specified, it does not when the model is misspecified. The DBW achieves this even when the model is misspecified. The following proposition summarizes this result.

Proposition 1 (Multivariate distribution balancing property): When the propensity scores and resulting inverse probability weights are estimated using (4), the weights minimize the following KL divergence for the multivariate covariate distribution between target and weighted groups in the limit even when the propensity score model is misspecified.

(5)$$\hat{\beta}_{\rm DBW} \ { {{\,p} \longrightarrow} }\beta_{\rm DBW}^\ast $$

(6)$$\beta_{\rm DBW}^\ast = \mathop{\rm arg min}\limits_{\beta}\; {\rm KL}\left(\,f( {\bf x}) ,\; {\,f_1( {\bf x}) \over \pi( {\bf X},\; \beta) } \right),\; $$

where $\beta _{\rm DBW}^\ast$ is the probability limit of $\hat {\beta }_{\rm DBW}$ and f(x) and f ₁(x) is the multivariate covariate distribution for the target group (${\cal S} = \{ i \mid R_i \in \{ 0,\; \, 1 \} \}$) and the response group (${\cal S}_{1} = \{ i \mid R_i = 1 \}$), respectively.

Proof. For any β, the following equations hold:

(7)$$\sum_{i = 1}^n \left\{{R_i\over \pi( X_i,\; \beta) } + \log( \pi( X_i,\; \beta) ) \right\} \ \ { {p} \longrightarrow} \int_{{\bf x} \in \chi} f( {\bf x}) \left\{{R\over \pi( {\bf x},\; \beta) } + \log( \pi( {\bf x},\; \beta) ) \right\}d{\bf x}$$

(8)$$ \ \ \ \ \ = \int_{{\bf x} \in \chi} \left\{{\,f_1( {\bf x}) \over \pi( {\bf x},\; \beta) } + f( {\bf x}) \log( \pi( {\bf x},\; \beta) ) \right\}d{\bf x}$$

(9)$$ = {\rm KL}\left(\,f( {\bf x}) ,\; {\,f_1( {\bf x}) \over \pi( {\bf X},\; \beta) } \right) + {\rm constant} .$$

By applying the continuous mapping theorem, one can obtain

(10)$$\mathop{\rm arg min}\limits_{\beta}\; \sum_{i = 1}^n \left\{{R_i\over \pi( X_i,\; \beta) } + \log( \pi( X_i,\; \beta) ) \right\} \ { {p} \longrightarrow} \mathop{\rm arg min}\limits_{\beta}\; {\rm KL}\left(\,f( {\bf x}) ,\; {\,f_1( {\bf x}) \over \pi( {\bf X},\; \beta) } \right).$$

□

This property is advantageous because it implies that bias is bounded above as stated in the following proposition.

Proposition 2 (Minimizing a sharp upper bound of bias): Let the absolute value of bias be denoted as B and $C = \max ( \vert {\mathbb E}[ Y_i - \mu \mid R_i = 1,\; \, {\bf x}] \vert )$. A function of imbalance in the multivariate covariate distribution quantified by the KL divergence sharply bounds bias from above as follows:

(11)$$B \leq C \; \; {\rm KL} \left(\,f( {\bf x}) - {\,f_1( {\bf x}) \over \pi( {\bf X},\; \beta_{\rm DBW}^\ast ) } \right)$$

Proof.

(12)$$B = \left\vert \int_{{\bf x} \in \chi} {\mathbb E}[ Y_i - \mu \mid R_i = 1,\; {\bf x}] \left(\,f( {\bf x}) - {\,f_1( {\bf x}) \over \pi( {\bf X},\; \beta_{\rm DBW}^\ast ) } \right)d{\bf x} \right\vert $$

(13)$$\leq C \int_{{\bf x} \in \chi} \left\vert f( {\bf x}) - {\,f_1( {\bf x}) \over \pi( {\bf X},\; \beta_{\rm DBW}^\ast ) } d{\bf x} \right\vert $$

(14)$$\leq C \; \; {\rm KL} \left(\,f( {\bf x}) - {\,f_1( {\bf x}) \over \pi( {\bf X},\; \beta_{\rm DBW}^\ast ) } \right),\; $$

where the last inequality follows from Pinsker's inequality, and the equalities hold when the propensity score model is correct. □

3.5 Minimizing an upper bound of the mean squared error

This section demonstrates that the proposed method sharply bounds an upper bound of the MSE by showing that the proposed method minimizes a sharp upper bound of the relative error of weights (proposition 3) and the relative error of weights bounds the MSE of the IPW estimator (proposition 4).

First, the KL divergence between the weights, which the proposed method minimizes as a loss function, bounds the relative error of the weights $\xi ( \hat {\beta }) = {\mathbb E} [ ( \pi ( X_i) / \pi ( X_i,\; \, \hat {\beta }) - 1 ) ^2 ]$ from above as the following proposition states:

Proposition 3 (Minimizing a sharp upper bound of the relative error of weights): When the propensity scores and resulting inverse probability weights are estimated using (4), a sharp upper bound of the relative error of weights is minimized in the limit even when the propensity score model is misspecified. If $\pi ( X_i,\; \, \hat {\beta }) \geq a \pi ( X_i)$ almost surely for some constant a ∈ (0,  1], the bound is as follows:

(15)$$\xi( \hat{\beta}) \leq 2\; {\mathbb E}\left[{\rm KL}\left({R_i\over \pi( X_i) },\; {R_i\over \pi( X_i,\; \hat{\beta}) } \right)\right]\quad {\rm for\ }a = 1$$

(16)$$\xi( \hat{\beta}) \leq {a^2 - 2a + 1\over a ( 1 - a + a \log( a) ) }\; {\mathbb E}\left[{\rm KL}\left({R_i\over \pi( X_i) },\; {R_i\over \pi( X_i,\; \hat{\beta}) } \right)\right]\quad {\rm for\ }0 < a < 1 ,\; $$

where the equalities hold when $\pi ( X_i,\; \, \hat {\beta }) = a \pi ( X_i)$, for all i. Proof is available in Supplementary Material C.1.

Note that the MLE does not minimize this upper bound as its coefficients converge to different limits than the DBW unless the propensity score model is correctly specified. Note also that this bound is tighter than Tan (Reference Tan2020)'s bound (Supplementary Material C.1), where the leading term is 5/(3a) for a case where $\pi ( X_i,\; \, \hat {\beta }) \geq a \pi ( X_i)$ almost surely for some constant a ∈ (0,  1/2] instead of a ∈ (0,  1] in proposition 3. Tan (Reference Tan2020) neither provide an intuitive interpretation of this relative error nor investigate the estimator minimizing it as the loss function.

Second, the MSE of the IPW estimator is bounded above by the relative error of weights as the following proposition states, which is also shown in Proposition 2 of Tan (2020) (a typo in the original proposition corrected):

Proposition 4 (Minimizing a sharp upper bound of the mean squared error): The MSE of the IPW estimator is sharply bounded above by the relative error of weights as:

(17)$${\mathbb E}\left[\left({1\over n} \sum_{i = 1}^n {R_i Y_i\over \pi( X_i,\; \hat{\beta}) } - \mu \right)^2 \right]\leq {2c\over n \rho} + c\left(1 + {1\over n \rho} \right)\xi( \hat{\beta}) ,\; $$

where ${\mathbb E}[ Y^2 \mid X] \leq c$ and π(X _i) ≥ ρ almost surely for some constants c > 0 and ρ ∈ (0,  1). The equality holds when Y _i = μ and $\pi ( X_i) = \pi ( X_i,\; \, \hat {\beta }) = 0.5$ for all i. Proof is available in Supplementary Material C.2.

Propositions 3 and 4 together imply that the proposed method minimizes the upper bound of the MSE even when the propensity score model is misspecified.

3.6 Comparison with the existing methodologies

This section briefly discusses related literature and compares the proposed methods with existing ones. A more detailed discussion including a comparison with the non-parametric methods is provided in Supplementary Material F.

Table 1 summarizes the loss function and three important properties of proposed and three existing methods, the standard IPW, calibrated weighting (Tan, Reference Tan2020), and entropy balancing (Hainmueller, Reference Hainmueller2012), for the average outcome estimation, where each of the DBW, calibrated weighting, and entropy balancing has two of those three properties. First, the distribution balancing property minimizes the imbalance in the multivariate covariate distribution between the target and weighted groups (Section 3.3) and an upper bound of the mean squared error (MSE) of the parameter of interest (Section 3.5). The DBW and entropy balancing have this property because they use the same loss function (4) with different link functions, where the DBW uses the logistic function and the entropy balancing uses the exponential function (Wang and Zubizarreta, Reference Wang and Zubizarreta2020, See also, Supplementary Material D and F). Second, the calibrated weighting and entropy balancing have the mean balancing property, which is increasingly employed in recently proposed methods (Imai and Ratkovic, Reference Imai and Ratkovic2014; Vermeulen and Vansteelandt, Reference Vermeulen and Vansteelandt2015; Zubizarreta, Reference Zubizarreta2015; Zhao, Reference Zhao2019; Wang and Zubizarreta, Reference Wang and Zubizarreta2020). This property implies that the estimators are sample-bounded (Robins et al., Reference Robins, Sued, Lei-Gomez and Rotnitzky2007; Wang and Zubizarreta, Reference Wang and Zubizarreta2020). Third, the estimated propensity score is bounded between 0 and 1 except for the entropy balancing, whose implicit propensity scores can take larger values than one (See, Wang and Zubizarreta, Reference Wang and Zubizarreta2020). This unboundedness implies that the implicit propensity score model in the entropy balancing is always misspecified, and it cannot be a doubly robust estimator for the average outcome estimation (Zhao, Reference Zhao2019), which contrasts with the double robustness property for the estimation of the average treatment effect for the treated (Zhao and Percival, Reference Zhao and Percival2017).

Table 1. Comparison of the proposed and existing methods for the average outcome estimation

Notes: The standard IPW is the IPW with the MLE.

The comparison shows that the drawback of the DBW is the lack of mean balancing and sample-boundedness properties. However, the DBW can incorporate or substitute for these properties. First, it can substitute the mean balancing with outcome models as explained in Section 4.3, which removes bias under the same assumptions as the mean balancing methods such as the calibrated weighting and entropy balancing. I also demonstrate that incorporating the linear outcome model is equivalent to modify the DBW weights to have the mean balancing property. Second, it can incorporate the normalization constraint to gain the sample-boundedness as explained in Section 4.1.

4.1 Normalization, regularization, and estimation

Like the IPW with the MLE, a drawback of the DBW is that estimated weights may not sum to one. To become normalized and thus sample-bounded, the normalized DBW (nDBW) estimator incorporates the normalization restriction as follows:

(18)$$\hat{\beta}_{\rm nDBW} = \mathop{\rm arg min}\limits_{\beta}\; \sum_{i = 1}^n \log( \pi( X_i,\; \beta) ) \quad{\rm s.t.}\; \; \sum_{i = 1}^n \left({R_i\over \pi( X_i,\; \beta) } - 1 \right) = 0 .$$

This nDBW estimator minimizes the distribution imbalance and an upper bound of the MSE under the normalization constraint. When the propensity score model is correctly specified, the limiting values of the DBW and nDBW estimators are the same because the (unnormalized) DBW estimator asymptotically satisfies the normalization constraint.

Using the Lagrangian function, the solution to this constrained optimization problem is obtained as the solution to the following unconstrained optimization problem:

(19)$$\hat{\beta}_{\rm nDBW} = \mathop{\rm arg min}\limits_{\beta}\; \left\{\sum_{i = 1}^n \log( \pi( X_i,\; \beta) ) + {\sum_{i = 1}^n ( 1 - \pi( X_i,\; \beta) ) \over n_0} \sum_{i = 1}^n \left({R_i\over \pi( X_i,\; \beta) } - 1 \right)\right\}$$

(20)$$n_0 = \sum_{i = 1}^n ( 1 - R_i) .$$

To estimate coefficients for the nDBW in (19), we can use the M-estimation under the regularity conditions that the solution to (19) is unique and an element of the interior of a compact set. The second condition is satisfied when the objective function in (19) is bounded below, which is similar to the nonseparation condition for the MLE (Vermeulen and Vansteelandt, Reference Vermeulen and Vansteelandt2015; Tan, Reference Tan2020). To satisfy this condition, we can utilize the L2 regularization:

(21)$$\hat{\beta}_{\rm nDBW} = \mathop{\rm arg min}\limits_{\beta}\; g_{\rm nDBW}( \beta) $$

(22)$$g_{\rm nDBW}( \beta) = \sum_{i = 1}^n \log( \pi( X_i,\; \beta) ) + {\sum_{i = 1}^n ( 1 - \pi( X_i,\; \beta) ) \over n_0} \sum_{i = 1}^n \left({R_i\over \pi( X_i,\; \beta) } - 1 \right) + {1\over 2} \lambda^{{\sf T}} \Vert \beta_{-1} \Vert ^2 ,\; $$

where λ is a hyper-parameter controlling for the regularization, β ₋₁ is β without the intercept, and || ⋅ || denotes the Euclidean norm.

As the properties of the M-estimator, the following results regarding the consistency and asymptotic normality of the nDBW estimator are obtained (Stefanski and Boos, Reference Stefanski and Boos2002).

Proposition 5 (Large sample properties): Under the conditions that the solution to (21) is unique and that (21) is bounded below, $\hat {\beta }_{\rm nDBW}$ and resulting inverse probability weights $1 / \pi ( X_i,\; \, \hat {\beta }_{\rm nDBW})$ converge in probability to their limiting values $\hat {\beta }_{\rm nDBW}^\ast$ and $1 / \pi ( X_i,\; \, \hat {\beta }_{\rm nDBW}^\ast )$, and since (21) is a smooth function, $\hat {\beta }_{\rm nDBW}$ is asymptotically normally distributed with a certain variance Σ as follows:

(23)$$\sqrt{n} ( \hat{\beta}_{\rm nDBW} - \beta^\ast ) \ \, { {d} \longrightarrow} {\cal N}( 0,\; \Sigma) $$

(24)$$\Sigma = g''_{\rm nDBW}( \beta) ^{-1} {\mathbb E}[ g'_{\rm nDBW}( \beta) g'_{\rm nDBW}( \beta) ^{\sf T}] g''_{\rm nDBW}( \beta) ^{-{\sf T}},\; $$

where β* is the true values of the coefficients, which is the solution for (3), g′_nDBW (β) = ∂g _nDBW(β)/∂β and $g''_{\rm nDBW}( \beta ) = \partial g'_{\rm nDBW}( \beta ) / \partial \beta ^{\sf T}$ are the first and second derivatives of g _nDBW(β) with respect to β. When the propensity score model is correctly specified, $\hat {\beta }_{\rm nDBW}$ and $1 / \pi ( X_i,\; \, \hat {\beta }_{\rm nDBW})$ are consistently estimated.

To demonstrate the benefits of normalization and regularization, I conduct a simulation in Section 2 for the non-normalized DBW and the nDBW with regularization. First, the RMSEs of the estimated weights are 1.23 for the non-normalized DBW against 1.18 for the nDBW, indicating that the normalization improves the weight estimation.

Second, Figure 3 presents the performance of the proposed estimators with different levels of regularization. The optimal value of the regularization hyper-parameter is 0.22 as depicted by the vertical dotted line. Beyond this value, the RMSE of the weights increases as the regularization becomes stronger, approaching the RMSE of the uniform weights shown in the horizontal dotted line. Although tuning hyper-parameters is still a difficult and unsolved open problem (Zubizarreta, Reference Zubizarreta2015; Wong and Chan, Reference Wong and Chan2017; Zhao, Reference Zhao2019), I recommend using regularization mainly to avoid the separation problem by choosing the smallest one that addresses the separation problem as I do in the simulation studies in Section 5. Another way is choosing the one that minimizes the upper bound of bias estimated by the kernel balance method in Hazlett (Reference Hazlett2020), which I adopt in the empirical application in Section 6.

Figure 3. The performance of the proposed estimators with different levels of regularization Notes: The horizontal dotted line indicates the RMSE of the uniform weights.

4.2 Algorithm

In this section, I explain the algorithm for estimating the DBW without the normalization and regularization for simplicity, and the algorithm for the nDBW with regularization is explained in Supplementary Material E. Since the loss function in (4) is non-convex with respect to β, the convex optimization algorithm such as the Newton's method is not applicable. For this non-convex optimization, I develop a majorization-minimization algorithm that iteratively optimizes two convex decompositions of the original function (Wu and Lange, Reference Wu and Lange2010). It exploits the characteristics that the minimization problem of the non-convex loss function of (4) is equivalent to the minimization problem of the difference of the two convex functions as follows:

(25)$$ \hskip-2pc g( \beta) \equiv \sum_{i = 1}^n \left\{{R_i\over \pi( X_i,\; \beta) } + \log( \pi( X_i,\; \beta) ) \right\} = g_1( \beta) - g_2( \beta) $$

(26)$$g_1( \beta) = \sum_{i = 1}^n \left\{{R_i\over \pi( X_i,\; \beta) } + ( 1 - R_i) \log( \pi( X_i,\; \beta) ) - ( 1 - R_i) \log( 1 - \pi( X_i,\; \beta) ) \right\}$$

(27)$$ \hskip-1pc g_2( \beta) = \sum_{i = 1}^n \left\{-R_i \log( \pi( X_i,\; \beta) ) - ( 1 - R_i) \log( 1 - \pi( X_i,\; \beta) ) \right\}.$$

This equivalence implies that the loss function for the DBW (25) is equivalent to the difference of the convex loss functions for the calibrated weighting (26) and the MLE (27). Note that the DBW does not minimize the absolute difference of the two functions but it minimizes the sum of the loss functions of the calibrated weighting and the negative of the loss function of the MLE. Thus, we should not interpret the DBW as the combination of the calibrated weighting and the MLE but rather the calibrated weighting should be regarded as the combination of the DBW and the MLE (Tan, Reference Tan2020). Since the DBW is better at controlling the mean squared error than the MLE as shown in Section 3.5, the DBW is also better than the calibrated weighting in this regard.

The difference of the convex functions algorithm for the DBW repeats the following two steps until convergence. First, as the majorization step, it constructs a surrogate function u(β,  β _t) for iteration t by the Taylor expansion of the second component of (25) as follows:

(28)$$u( \beta,\; \beta_t) = g_1( \beta) - \{ g_2( \beta_t) + g'_2( \beta_t) ^{\sf T} ( \beta - \beta_t) \} ,\; $$

where β _t is the initial values or values obtained in the previous iteration and $g'_2( \beta ) = {\partial g_2( \beta ) \over \partial \beta }$. Next, as the minimization step, it estimates β _t+1 such that minimizes (28) with respect to β while keeping β _t fixed. This minimization is easily conducted, e.g., via the Newton's method, because the objective surrogate function is convex. This algorithm is proved to decrease the loss monotonically in every steps and converge to a stationary point (Wu and Lange, Reference Wu and Lange2010).

4.3 A substitute for the mean balancing condition: doubly robust distribution balancing weighting estimator

Another drawback of the DBW is that it lacks the mean balancing property, which is increasingly employed by recently proposed methods as a safeguard against model misspecification (Hainmueller, Reference Hainmueller2012; Imai and Ratkovic, Reference Imai and Ratkovic2014; Vermeulen and Vansteelandt, Reference Vermeulen and Vansteelandt2015; Zubizarreta, Reference Zubizarreta2015; Zhao, Reference Zhao2019; Tan, Reference Tan2020; Wang and Zubizarreta, Reference Wang and Zubizarreta2020). This removes bias when the true outcome model is a linear regression model with the balanced covariates (Zubizarreta, Reference Zubizarreta2015; Zhao and Percival, Reference Zhao and Percival2017; Zhao, Reference Zhao2019; Fan et al., Reference Fan, Imai, Lee, Liu, Ning and Yang2021).

To substitute the mean balancing property, the (n)DBW can incorporate the outcome regression model. Specifically, this study considers the augmented IPW where the weights are estimated by the (n)DBW instead of the MLE. The augmented IPW combines the outcome model and propensity score model in the following manner (Robins et al., Reference Robins, Rotnitzky and Zhao1994; Seaman and Vansteelandt, Reference Seaman and Vansteelandt2018):

(29)$$\hat{\mu} = {1\over n} \sum_{i = 1}^n \left({R_i Y_i\over \pi( X_i,\; \hat{\beta}) } - {R_i - \pi( X_i,\; \hat{\beta}) \over \pi( X_i,\; \hat{\beta}) } \hat{m}( X_i) \right),\; $$

where $\pi ( X_i,\; \, \hat {\beta })$ is the estimated propensity scores and $\hat {m}( X_i)$ is the estimated conditional expectation function $m( X_i) = {\mathbb E}[ Y_i \mid X_i]$. The (n)DBW DR estimator uses the WLS to estimate m(X _i) because it has better performance than the OLS (Kang and Schafer, Reference Kang and Schafer2007; Imai and Ratkovic, Reference Imai and Ratkovic2014). Under the same condition for the mean balancing methods, the (n)DBW estimator eliminates bias by using a linear regression outcome model. It can also utilize more flexible methods such as machine learning techniques for estimating m(X _i). This estimator is doubly robust, i.e., consistent when either the propensity score or the outcome model is correctly specified (Vermeulen and Vansteelandt, Reference Vermeulen and Vansteelandt2015; Seaman and Vansteelandt, Reference Seaman and Vansteelandt2018). Moreover, the (n)DBW DR estimator achieves semiparametric efficiency bound if both of the models are correctly specified (Vermeulen and Vansteelandt, Reference Vermeulen and Vansteelandt2015; Seaman and Vansteelandt, Reference Seaman and Vansteelandt2018).

Recently, Chattopadhyay and Zubizarreta (Reference Chattopadhyay and Zubizarreta2023) demonstrates that the augmented IPW is expressed as the (non-augmented) IPW with modified weights with the covariate balancing property. Chattopadhyay and Zubizarreta (Reference Chattopadhyay and Zubizarreta2023) also proves that the modified weights minimize the Chi-square-type distance from the (inverse probability) weights used in the augmented IPW with the WLS. This implies that the (n)DBW DR estimator can be expressed as the (non-augmented) IPW estimator with weights $\hat {w}_{\rm nDBWDR}$ satisfying the following conditions:

(30)$$\hat{w}_{\rm nDBWDR} = \mathop{\rm arg min}\limits_{w}\; \sum_{i = 1}^n R_i \left({w_i - 1 / \pi( X_i,\; \hat{\beta}_{\rm nDBW}) \over 1 / \pi( X_i,\; \hat{\beta}_{\rm nDBW}) } \right)^2 \quad{\rm s.t.}\; \; \sum_{i = 1}^n ( R_i w_i - 1) X_i = 0 ,\; $$

which demonstrates that the (n)DBW DR estimator equips the covariate balancing property with the approximate multivariate covariate distribution balancing.

Since the covariate balancing methods satisfy ${1\over n} \sum _{i = 1}^n ( R_i - \pi ( X_i,\; \, \hat {\beta }) ) \, X_i = 0$ by construction, their DR estimators with the linear outcome model of $\hat {m}( X_i)$ cannot improve their non-DR estimators because the augmented term, the second term in (29), is always 0. I will examine the finite-sample performance of the nDBW DR estimators through simulation studies in Section 5.

6.1 Bias due to imbalance in multivariate covariate distribution

While bias induced by multivariate distributional imbalance is difficult to evaluate, an upper bound of bias can be estimated by the kernel technique (Hazlett, Reference Hazlett2020). I use this upper bound to compare the performance of various methods.

Figure 5 presents the upper bounds of bias due to multivariate distributional imbalance for the proposed (shown in red) and other estimators. Those for the calibrated weighting for the distance larger than 30 km are missing because the weights cannot be estimated due to the convergence issue. As the distance from the demarcation line increases, bias due to multivariate distributional imbalance becomes severe, but the proposed estimator mitigates the bias more effectively than the others.

Figure 5. Upper bounds of bias due to distributional imbalance for various methods, Notes: As the distance from the demarcation line increases, bias due to multivariate distributional imbalance becomes severe, but the proposed estimator mitigates it more effectively than the others.