Robust Estimating Method for Propensity Score Models and its Application to Some Causal Estimands: A review and proposal

In observational study, the propensity score has the central role to estimate causal effects. Since the propensity score is usually unknown, estimating by appropriate procedures is an indispensable step. A point to note that a causal effect estimator might have some bias if a propensity score model was misspecified; valid model construction is important. To overcome the problem, a variety of interesting methods has been proposed. In this paper, we review four methods: using ordinary logistic regression approach; CBPS proposed by Imai and Ratkovic; boosted CART proposed by McCaffrey and colleagues; a semiparametric strategy proposed by Liu and colleagues. Also, we propose the novel robust two step strategy: estimating each candidate model in the first step and integrating them in the second step. We confirm the performance of these methods through simulation examples by estimating the ATE and ATO proposed by Li and colleagues. From the results of the simulation examples, the boosted CART and CBPS with higher-order balancing condition have good properties; both the estimate of the ATE and ATO has the small variance and the absolute value of bias. The boosted CART and CBPS are useful for a variety of estimands and estimating procedures.

To estimate causal effects accurately, adjusting covariates (or confounders; hereafter, “covariates”) is one of the important steps in observational study. When all covariates are observed, the covariates can be adjusted, and an unbiased estimator for causal effects can be obtained; the situation of “no unmeasured confounding” (Hernán and Robins, 2020). In this situation, the propensity score has the central role to estimate causal effects (Rosenbaum and Rubin, 1983). However, the propensity score needs to be estimated by appropriate procedures since it is usually unknown. A point to note that an interested estimator might have some bias if a propensity score model was misspecified (Kang and Schafer, 2007); valid model construction is important.

Recently, some causal estimands, the target population of a causal effects estimation, are considered. One of the well-known estimands is the average treatment effects (ATE; Imbens and Rubin, 2015). The ATE assumes causal effects that the overall population is treated “counterfactual” study treatment and control treatment. In short, the ATE is interested in the overall causal effects. A different estimand called the average treatment effects for the overlap population (ATO) is proposed by Li et al. (2018). The ATO has several attractive features compared with the ATE (Li et al., 2018); for instance, the ATO may become more reasonable causal effects in the sense that subjects who could change their actual treatment to the “counterfactual” treatment have larger weights than subjects who could not change, and is less sensitive to extreme propensity scores (i.e. low propensity scores in the treatment group, or high propensity scores in the control group). Another estimand is defined by the incremental propensity score (Kennedy, 2019). The estimands are defined how change in the intervention probability (i.e., the propensity score) affects the causal effects. There are some attractive estimands, however, there is a crucial problem; some estimands depend on the “true” propensity score. For instance, as mentioned Mao et al. (2019), the ATO has some bias when the propensity score model is misspecified. Therefore, the similar problem mentioned in Kang and Schafer (2007) may be occurred.

There may be a critical problem for not only estimators of causal effects but also estimands itself when the propensity score model is misspecified. To overcome the problem, misspecification of the propensity score model, a variety of interesting methods has been proposed. Covariate balancing methods are fully parametric strategy (Imai and Ratkovic, 2014; Fong et al., 2018). The idea is intuitive and easy to understand; a parametric propensity score model is estimated as the covariate balance between treatment groups is balanced sufficiently. This is somewhat similar procedure as calibration methods (Zubizarreta, 2015; Wang and Zubizarreta, 2020) or multiply robust estimations (Han, 2014; Orihara and Hamada, 2019), however, to my understanding, the philosophy is different; the former is estimating weights nonparametrically so that each group’s weighted expectation of outcome model becomes the marginal expectation of each potential outcome, and the latter is estimating weights so that the weighted score statistics for propensity score models and outcome models become $0$. Nonparametric strategy such as machine learning methods is also considered widely (McCaffrey et al., 2004; Lee et al., 2010; Westreich et al., 2010; Cannas and Arpino, 2019). As mentioned in Cannas and Arpino (2019), by using a sufficient complex model, checking the balancing property (Rosenbaum and Rubin, 1983) may not be that much of a concern. Although the strategy of covariate balancing methods and machine learning methods are different, objective is similar. Recently, a method that is an intermediate between parametric and nonparametric methods has been proposed (Liu et al., 2018). Liu et al. (2018) mentioned that the proposed method is semiparametric in the sense that the method has both sides. In the parametric side, it is assumed that (high dimensional) covariates information is explained by (low dimensional) linear combinations of the covariates. This spirit is the same as a sufficient dimension reduction (Li, 1991). In the nonparametric side, it is assumed that the explicit form of link functions is not assumed.

In this paper, we compare the explained three methods: the covariate balancing methods (CBPS), the machine learning methods (boosted CART (bCART) algorithm), and the method proposed in Liu et al. (2018). Also, we propose an another parametric strategy in this paper: a multiply robust estimator for the propensity score (note that “multiply robust” is slightly different from Han, 2014). In brief, we prepare some candidate models, and combine the candidate models all at once (do not have to select just one); we call the method as “integrated method”. To combine the candidate models, we consider two procedures. One is a well-known model averaging approach (MA; c.f. Hoeting et al., 1999; Hjort and Claeskens, 2003; Xie et al., 2019), the other one is novel which has some good theoretical properties. Note that the former is called “MA-based integrated method”, the latter is “proposed integrated method”. For the proposed integrated method, we confirm theoretical properties. If the correct model was included in the candidate models, the parameter estimator would have the consistency. Also, we consider the situation where the true model is not included in the candidate models. Under this situation, we confirm the condition where the estimator becomes a “valid” in the sense of the true Kullback-Leibler (KL) divergence; this property is different from MA. By using the integrated methods, we consider a robust estimating method to estimate a “valid” causal effects without any model/variable selection methods. Therefore, these methods may overcome the problem of a variable selection of a propensity score model (Brookhart et al., 2006; Austin et al., 2007). For instance, a researcher claims that some covariates have to be included in a propensity score model, whereas, another researcher claims that the covariates do not have to be included. Commonly, in this situation, one of the researchers has to compromise, or the researchers use some model/variable selection methods. Meanwhile, integrated procedure provides a simple solution: the both models can be included in a one integrated model. Note that no methods explained here are suffer from the “propensity score tautology” (Imai and Ratkovic, 2014), that is, we need not take care of model adjustment to satisfy the “balancing property” (Rosenbaum and Rubin, 1983). In this sense, the methods have reasonable to apply real-world data since there may be too many covariates to construct an appropriate propensity score model.

The remainder of the paper proceeds as follows. In section 2, we review the previous methods, and explain the proposed methods briefly. Regarding the proposed integrated method, we confirm the theoretical properties as a new robust generalized linear model (GLM) estimator with and without the condition where the true model is included in the candidate models. Since the propensity score model is assumed as not only a logistic regression model but also a regression model with normal error under general treatment regimes (Hirano and Imbens, 2004 and Imai and van Dyk, 2004), considering the general situation is valuable. In section 3, we show simulation data examples when the ATE and ATO is central interest, respectively. The simulation datasets refers to well-established simlation examples (Setoguchi et al., 2008; Setodji et al.,2017; Huang et al., 2022). Some materials are found in Appendix.

Let $n$ be the sample size. $T_{i}\in\{0,\,1\}$, $\boldsymbol{X}_{i}\in\mathbb{R}^{p}$ and $(Y_{1i},\,Y_{0i})\in\mathcal{Y}\subset\mathbb{R}^{2}$ denote the treatment, a vector of covariates measured prior to treatment, and potential outcomes, respectively. $Y_{i}:=T_{i}Y_{1i}+(1-T_{i})Y_{0i}$ denotes an observed outcome and we assume that $i=1,\,2,\,\dots,\,n$ are i.i.d. samples and the stable unit treatment value assumption (c.f. Rosenbaum and Rubin, 1983) holds.

Next, we introduce a general class of estimands called “weighted average treatment effect” (WATE) in Hirano et al. (2003) and Li et al. (2018). The ATE conditional on $\boldsymbol{x}$ is defined as $\tau(\boldsymbol{x}):={\rm E}[Y_{1}-Y_{0}|\boldsymbol{x}]$, and the WATE is defined as

$$\tau_{h}:=\frac{{\rm E}[\tau(\boldsymbol{X})h(\boldsymbol{X})]}{{\rm E}[h(\boldsymbol{X})]},$$

where $h(\cdot)$ is a weight function for $\boldsymbol{x}$. When $h(\boldsymbol{x})\equiv 1$, the WATE becomes the ATE: $\tau_{ATE}:={\rm E}[Y_{1}-Y_{0}]$ obviously. Under the strongly ignorable treatment assignment assumption $T\mathop{\perp\!\!\!\perp}(Y_{1},\,Y_{0})|\boldsymbol{X}$, we can estimate the ATE using the propensity score $e(\boldsymbol{X}_{i}):={\rm P}(T_{i}=1|\boldsymbol{X}_{i})$. When $h(\boldsymbol{x})=e(\boldsymbol{x})(1-e(\boldsymbol{x}))$, the WATE becomes the ATO (Li et al., 2018):

$$\tau_{ATO}:=\frac{{\rm E}[\tau(\boldsymbol{X})e(\boldsymbol{X})(1-e(\boldsymbol{X}))]}{{\rm E}[e(\boldsymbol{X})(1-e(\boldsymbol{X}))]}.$$

Since the function $e(1-e),\,e\in(0,1)$ is the convex and symmetric function at $1/2$, subjects who could change their actual treatment to the “counterfactual” treatment (i.e. $e\approx 1/2$) have larger weights than subjects who could not change (i.e. $e\approx 0$ or $1$).

As mentioned in Introduction, the propensity score is commonly unknown and need to be estimated. If a propensity score model is misspecified, an estimated ATE may have serious bias. Also, as obviously, the ATO depends on the true propensity score; if the estimated propensity score model is misspecified, the ATO has some bias (see also Mao et al., 2019). In the following subsections, we introduce and propose robust propensity score estimating methods.

As mentioned in the previous section, there are many robust methods for propensity score model misspecification. In this section, we confirm their performance through simulation data examples. About the example data set, we confirm the performance when we are interested in the ATE and ATO, respectively. Note that an iteration time of all simulations is 500. Some simulation results are appeared in Appendix C.

At first, we explain the data generating mechanism in this simulation. The setting is referering to Wyss et al. (2014); Setodji et al. (2017); Mao et al. (2019); Huang et al. (2022). There are 10 covariates generated from the following distributions:

	$\displaystyle X_{i1},\ X_{i3}$	$\displaystyle\stackrel{{\scriptstyle i.i.d.}}{{\sim}}Ber(0.5),$
	$\displaystyle\left(X_{i2},X_{i4},X_{i7},X_{i10}\right)^{\top}$	$\displaystyle\stackrel{{\scriptstyle i.i.d.}}{{\sim}}N(\boldsymbol{0}_{4},{\rm I}_{4}),$

$$X_{i5}\sim Ber(p_{i5}),\ X_{i6}\sim Ber(p_{i6}),\ X_{i8}\sim Ber(p_{i8}),\ X_{i9}\sim Ber(p_{i9}),$$

where

$$p_{i5}=expit\left\{0.4(X_{i1}-1)\right\},\ p_{i6}=expit\left\{5X_{i2}\right\},\ p_{i8}=expit\left\{0.4(X_{i3}-1)\right\},\ p_{i9}=expit\left\{5X_{i4}\right\}.$$

From the setting, the correlations between the covariates become approximately as follows:

$$Corr(X_{1},X_{5})=Corr(X_{3},X_{8})=0.2,\ Corr(X_{2},X_{6})=Corr(X_{4},X_{9})=0.8.$$

Note that the correlations of the other combinations are approximately $0$. Next, we introduce the true propensity score model. The model is the same as Option E of Setodji et al. (2017):

$$e(\boldsymbol{X}_{i})={\rm P}\left(T=1|\boldsymbol{X}_{i}\right)=expit\left\{h_{T}^{\top}\left(\boldsymbol{X}_{i}\right)\boldsymbol{\beta}\right\},$$

where $\beta_{k}\sim U(-0.8,0.8)$ ($k=1,2,\dots,12$), and

$\displaystyle h_{T}^{\top}\left(\boldsymbol{X}_{i}\right)$

$\displaystyle=\left(X_{i1},\dots,X_{i7},X_{i1}\times X_{i3},X_{i2}\times X_{i4},X_{i4}\times X_{i5},X_{i5}\times X_{i6},X_{i2}^{2}\right).$

Therefore, there is obviously nonlinear structure in the propensity score; the standard linear logistic regression may have some bias. Finally, we introduce the true outcome model:

$$Y_{i}=\left(0.5+X_{i2}+X_{i1}\times X_{i2}\right)T_{i}+h_{Y}^{\top}\left(\boldsymbol{X}_{i}\right)\boldsymbol{\xi}+\varepsilon_{i},\ \ \ \varepsilon_{i}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}N(0,1),$$

where $\xi_{k}\sim U(-0.73,0.73)$ ($k=1,2,\dots,8$), and

$$h_{Y}^{\top}\left(\boldsymbol{X}_{i}\right)=\left(1,\exp\left\{X_{i1}\right\},\exp\left\{X_{i2}\right\},\exp\left\{X_{i3}\right\},\exp\left\{1.3X_{i4}\right\},X_{i8},X_{i9},X_{i10}\right).$$

Therefore, there are the heterogeneous treatment effects, and nonlinear structure in the outcome model.

From the settings, $X_{i1}$, $X_{i2}$, $X_{i3}$, and $X_{i4}$ are confounders; it is at least necessary to adjust the covariates for estimating unbiased causal effects. The true causal effect for the ATE is approximately 0.500, and for the ATO is approximately 0.430; the methods which estimate the true causal effect acculately are more appropriate.

In this paper, we review the four methods, propose the novel robust two step strategy called “integreted method”. The proposed strategies have two steps: estimating each candidate model in the first step and integrating them in the second step. Regarding the proposed integrated method, it is proved that the estimator has the $\sqrt{n}$–consistency when the true model is included in the candidate. In this sense, the proposed method has multiple robustness. Also, the KL divergence of the integrated model is smaller than each candidates when the true model is not included. In the simulation examples, the performance of these methods are compared in the situation where the estimand is the ATE and ATO, respectively. From the results of the simulation examples, the boosted CART and CBPS with higher-order balancing condition have good properties; both the estimate of the ATE and ATO has the small variance and the absolute value of bias. The boosted CART and CBPS are useful for a variety of estimands and estimating procedures.

Properties of boosted CART and CBPS are discussed in the various previous researches (e.g., Lee et al., 2010; Westreich et al., 2010; Wyss et al., 2014; Huang et al., 2022), however, to the best of our knowledge, the estimator for the ATO is not considered. Mao et al. (2019) mentioned that the impact of the proponsity score model misspecification is slight compared with the ATE. However, it is appeared that there are some impact to estimates, and the impact is different between estimating methods through our simulation result. Therefore, in spite of the estimand or estimating methods, it is appeared again that the propensity score model consideration is one of the important process in causal inference.

A different integrated method is considered in Xie et al. (2019), however, the motivation of their proposed method is slight different from our proposal. In Xie’s methods, only one parametric model and one nonparametric model are integrated; in other words, taking good point of the two models. Whereas, our proposed methods are consist of parametric models only; this is one of the strong points. As mentioned in Section 2, only related to GLM procedure such as “glm” function in R, and “GENMOD” procedure in SAS is used. In other words, any special training is not necessary; only “integration” of existing knowledge for many statistics users is necessary.

In this paper, we consider robust estimation for the propensity score misspecification. On the other hand, propensity score subclassification (e.g., Wang et al., 2016; Orihara and Hamada, 2021) is another solution to ease the propensity score misspecicication problem since the propensity scores in each subclass are regarded as “smoothed” (Imbens and Rubin, 2015). In other words, the impact of misspecification is also smoothed. Subclassification estimator is used only for the ATE or ATT (e.g., Antonelli et al., 2018). By extending subclassification estimator, another robust ATO estimator may be derived.