Matching(subclassification/coarsened)

The very beginning idea of estimating ATE (under the Assumption 1, 2 and 3) by first finding pairs of observations with covariates close enough to each other in two arms, and then average them out: $$\tau =E_{\mathcal{X}}\{E[Y|\mathbf{X}\in {\mathcal{X}}_j,T=1]-E[Y|\mathbf{X}\in {\mathcal{X}}_j,T=-1]\}$$ where ${\mathcal{X}}_j$ is a subspace such that $\bigcup _{j=1}^J{\mathcal{X}}_j=\mathcal{X}$ and ${\mathcal{X}}_i\bigcap {\mathcal{X}}_j=\varphi$ for any $i\ne j$. One-to-one paired matching is hard to find, especially when $\mathcal{X}$ is of high dimension. So people propose to use subclassification or coarsened, which also provide more robust estimations.

Inverse Probability Weighting (IPW)

It can be shown that (Proof in Appendix) $$E[Y^{(t)}]=E\left[\frac{1[T=t]Y}{Pr(T=t\mid \mathbf{X})}\right]$$ which leads to $$\begin{array}{}\text{(2)}& \tau =E[Y^{(1)}-Y^{(-1)}]=E\left[\frac{1[T=1]Y}{\pi (\mathbf{X})}\right]-E\left[\frac{1[T=-1]Y}{1-\pi (\mathbf{X})}\right].\end{array}$$ Thus, the estimator of ATE based on observed dataset is $$\hat{\tau }=\frac{1}{n}\sum _{i:t_i=1}\frac{y_i}{\hat{\pi }({\mathbf{x}}_i)}-\frac{1}{n}\sum _{i:t_i=-1}\frac{y_i}{1-\hat{\pi }({\mathbf{x}}_i)}.$$ As this approach can be traced back to Horvitz and Thompson (1952), sometimes it is refered to as Horvitz–Thompson (HT) estimator (Imai and Ratkovic 2014). Since directly weighting on the inverse of propensity score may cause large variability when estimated $\pi ({\mathbf{x}}_i)$ is either close to 1 or 0, Hirano, Imbens, and Ridder (2003) propose to standardize the weight by $$\hat{\tau }=\sum _{i:t_i=1}\frac{y_i}{\hat{\pi }({\mathbf{x}}_i)}/\sum _{i:t_i=1}\frac{1}{\hat{\pi }({\mathbf{x}}_i)}-\sum _{i:t_i=-1}\frac{y_i}{1-\hat{\pi }({\mathbf{x}}_i)}/\sum _{i:t_i=-1}\frac{1}{1-\hat{\pi }({\mathbf{x}}_i)}.$$ In our simulation studies, to distinguish from HT, IPW is refered to this standardized weighting method. Similarly, for ATT, it can be estimated by $$\begin{array}{rl}{\tau }^{\text{ATT}}=& E[Y^{(1)}-Y^{(-1)}\mid T=1]=E[Y^{(1)}\mid T=1]-E[Y^{(-1)}\mid T=1]\\ =& \frac{1}{\pi }(E[1[T=1]Y^{(1)}]+E[1[T=-1]Y^{(-1)}]-E\left[\frac{1[T=-1]Y}{1-\pi (\mathbf{X})}\right])\\ =& \frac{1}{\pi }E\left[(1-\frac{1[T=-1]}{1-\pi (\mathbf{X})})Y\right]=\frac{1}{\pi }E\left[\left(\frac{1[T=1]-\pi (\mathbf{X})}{1-\pi (\mathbf{X})}\right)Y\right]\\ =& \frac{1}{n_1}(\sum _{i:t_i=1}{y}_i-\sum _{i:t_i=-1}\frac{\pi ({\mathbf{x}}_i)}{1-\pi ({\mathbf{x}}_i)}{y}_i)\text{for finite sample estimation}\end{array}$$ where the second equation comes from: $$\begin{array}{rl}& E\left[\frac{1[T=-1]Y}{1-\pi (\mathbf{X})}\right]=E[Y^{(-1)}]=E[Y^{(-1)}\mid T=1]Pr(T=1)+E[Y^{(-1)}\mid T=-1]Pr(T=-1)\\ \Rightarrow & E[Y^{(-1)}\mid T=1]=\frac{1}{\pi }(E\left[\frac{1[T=-1]Y}{1-\pi (\mathbf{X})}\right]-E[1[T=-1]Y^{(-1)}])\end{array}$$

Augmented Inverse Probability Weighting (AIPW)

J. M. Robins, Rotnitzky, and Zhao (1995) augmented the IPW by a weighted average of the outcome model $$\begin{array}{rl}{\tau }^{\text{AIPW}}=& \underset{\text{IPW}}{\underbrace{E\left[\frac{1[T=1]Y}{\pi (\mathbf{X})}\right]-E\left[\frac{1[T=-1]Y}{1-\pi (\mathbf{X})}\right]}}\\ \text{(3)}& & -\underset{\text{Augmentation}}{\underbrace{E[\frac{1[T=1]-\pi (\mathbf{X})}{\pi (\mathbf{X})}{\mu }_1(\mathbf{X})+\frac{1[T=1]-\pi (\mathbf{X})}{1-\pi (\mathbf{X})}{\mu }_{-1}(\mathbf{X})]}}\end{array}$$ This AIPW is called “doubly robust” because it is consistent as long as either the treatment assignment mechanism or the outcome model is correctly specified (Kurz 2021). If the propensity score is correctly specified, then $E(1[T=t]-Pr(T=t\mid \mathbf{X}))=E[E(1[T=t]-Pr(T=t\mid \mathbf{X})\mid \mathbf{X})]=0$ which simplifies AIPW to IPW no matter whether response surface functions ${\mu }_T()$ are correct. On the other hand, if response model ${\mu }_T()$ are correctly specified but propensity score is not correctly estimated, AIPW reduces to S- or T-learner (definitions are in the next section. Proof in Appendix)

Covariate Balancing Propensity Score (CBPS)

The very nature of weighting is it helps to make sure the confounding covariates distribution between the observed two arms are balanced. If covariate balancing is guaranteed, then the observed data can mimic the randomized trial data and the treatment effect can be directly obtained from observations. Imai and Ratkovic (2014) propose to directly target on the covariate balancing instead of weighting. So they developed a Covariate Balancing Propensity Score (CBPS) method which is a propensity score estimating method that yields the PSs improve the covariate balancing. The core idea of CBPS is fit a logistic PS model–maximize the likelihood function–subject to covariate balancing constraints. Specifically, the logistic PS model looks like $$\hat{\beta }=\arg \underset{\beta \in \mathrm{\Theta }}{min}-\sum _{i=1}^N{T}_i\log \pi ({\mathbf{X}}_i;\beta )+(1-T_i)\log (1-\pi ({\mathbf{X}}_i;\beta ))$$ which is equivalent to (take derivative w.r.t. $\beta$) $$E[\frac{T_i{\tilde{\mathbf{X}}}_i}{\pi ({\mathbf{X}}_i;\beta )}-\frac{(1-T_i){\tilde{\mathbf{X}}}_i}{1-\pi ({\mathbf{X}}_i;\beta )}]=0$$ where ${\tilde{\mathbf{X}}}_i={\pi }^{\prime }({\mathbf{X}}_i;\beta )$ which is the first derivative of $\pi ({\mathbf{X}}_i;\beta )$ over $\beta$. The balancing conditions are $$E\left[\frac{T_{i}f({\mathbf{X}}_i)}{\pi ({\mathbf{X}}_i;\beta )}\right]=E\left[\frac{(1-T_i)f({\mathbf{X}}_i)}{1-\pi ({\mathbf{X}}_i;\beta )}\right]$$ where $f()$ can be any type of function that intended to be balanced, including first moment $\mathbf{X}$, second moment ${\mathbf{X}}^2$, interaction ${\mathbf{X}}_i{\mathbf{X}}_j,i\ne j$, etc. When $f({\mathbf{X}}_i)={\tilde{\mathbf{X}}}_i$, CBPS is exactly the MLE of logistic regression. If ${\tilde{\mathbf{X}}}_i={\mathbf{X}}_i$ (a linear logistic model), the number of constraints is equal to the number of parameters, which is called exact CBPS. But we can manually add more constraints to make sure higher moments of covariates are balanced as well, then we have more constraints than the covariates, which is call over-identified CBPS. But “overidentifying restrictions generally improve asymptotic efficiency but may result in a poor finite sample performance” (Imai and Ratkovic 2014). The caveat of CBPS is, even PS model is completely misspecificed, the imposed moments (of $X$, $X_i{X}_j$) are still balanced. CBPS largely avoids extremely propensities and thus, is more robust and better balance.

CBPS is just one of the Covariate-Balancing-type method to obtain the weights. There are many other methods following this idea (Hainmueller 2012; Chan, Yam, and Zhang 2016; Pirracchio and Carone 2018; Huling and Mak 2020) but will not be introduced here.

TMLE

Another estimator that is receiving more attention for ATE estimation is the targeted maximum likelihood estimator (TMLE) (Van Der Laan and Rubin 2006), which also enjoys the double robustness as AIPW along with other desired features like efficient and aysmptotically linear (allows for construction of valid Wald-type confidence intervals). Either ${\hat{\mu }}_i(\cdot )$ and $\hat{\pi }(\cdot )$ are correctly specified leads to consistent estimator while if both are consistently estimated, TMLEs achieve the semi-parametric efficiency bound (Van Der Laan and Rubin 2006).

TMLE is a two-stage procedure with the first stage to obtain an initial estimate of response surface/conditional mean response ${\hat{\mu }}_1(\mathbf{X})$ and ${\hat{\mu }}_{-1}(\mathbf{X})$. If the initial estimator response surfaces is consistent, the TMLE remains consistent; but if the initial estimator is not consistent, the subsequent targeting step provides an opportunity for TMLE to reduce any residual bias in the estimate of the parameter of interest. Specifically, TMLE solves the efficient influence curve estimating equation for the target parameter (ATE here) in the second stage, which fluctuate/correct the initial estimation on the correct direction as a one-step estimator. The whole idea is supported by the semiparametric theory and here are some references: Semiparametric theory; Nonparametric efficiency theory in causal inference; influence functions.

$${\hat{\tau }}^{\text{TMLE}}(\mathbf{X})={\hat{\mu }}_1(\mathbf{X})-{\hat{\mu }}_{-1}(\mathbf{X})+\hat{\epsilon }(\frac{1[T=1]}{\hat{\pi }(\mathbf{X})}-\frac{1[T=-1]}{1-\hat{\pi }(\mathbf{X})})$$ where $\epsilon$ is a flunctuation parameter to update the initial estimation of response surfaces ${\hat{\mu }}_1(\mathbf{X})-{\hat{\mu }}_{-1}(\mathbf{X})$ using the information from propensity score, just like other doubly-robust estimator or X-learner. $\epsilon$ can be estimated by fitting ‘clever covariate’ $H(T,\mathbf{X})=\frac{1[T=1]}{\hat{\pi }(\mathbf{X})}-\frac{1[T=-1]}{1-\hat{\pi }(\mathbf{X})}$ on the residual $Y-\mathbb{E}[Y\mid T,\mathbf{X}]\sim H(T,\mathbf{X})$.

Interaction-oriented Methods

Without loss of generality, outcome $Y$ can be written by $$E(Y)=f(\mathbf{X})+1[T=1]g(\mathbf{X}),$$ then $E(Y\mid \mathbf{X}=\mathbf{x},T=1)=f(\mathbf{x})+g(\mathbf{x})$, $E(Y\mid \mathbf{X}=\mathbf{x},T=-1)=f(\mathbf{x})$, and therefore, $\tau (\mathbf{x})=g(\mathbf{x})$. So the problem of estimating treatment effect, or subgroup identification turns to be figuring out $g(\mathbf{x})$. Notice that here $1(T=1)$ is independent of covariates $\mathbf{X}$, reflecting the nature of randomized trial.

Also be aware that the covariates of $f(\cdot )$ and $g(\cdot )$ could be different. Denote $\mathcal{Z}$ and $\mathcal{V}$ are two subspaces of origincal covariate space $\mathcal{X}$, and $$E(Y)=f(\mathbf{V})+1[T=1]g(\mathbf{Z}),$$ where $\mathbf{Z}\in \mathcal{Z},\mathbf{V}\in \mathcal{V}$, and $\mathcal{V}$ and $\mathcal{Z}$ can have overlap. Generally, covariates in $\mathbf{V}$ are called diagnostic variables and $\mathbf{Z}$ are predictive variables (Imai, Ratkovic, et al. 2013). Some variables are predictive only, indicating that they only react to the treatment; some variables are diagnostic only, meaning that they do nothing with treatment but impact the outcome; some variables are both predictive and diagnostic, while some variables play no effect at all. See Figure (1) for illustration.

A bunch of methods are developed following this idea: Interaction trees (Su et al. 2008, 2009), GUIDE (Generalized unbiased interaction detection and estimation) (Loh 2002; Loh, He, and Man 2015), MOB (Model-Based Recursive Partitioning) (Zeileis, Hothorn, and Hornik 2008; Seibold, Zeileis, and Hothorn 2016).

Figure 1: Illustration of predictive and diagnostic variable

Q-learning

Q-learning gets its name because its objective function plays a role similar to that of the Q or reward function in reinforcement learnin (Li, Wang, and Tu 2021) and is first used in estimating optimal dynamic treatment regime (Murphy 2003; J. M. Robins 2004; Schulte et al. 2014). The basic idea is focusing on the estimation of response surfaces $E[Y\mid \mathbf{X},T]$, which is also known as g-computation (J. Robins 1986). So this approach is also called the parametric g-formula (Hernán and Robins 2010) in the literature.

A-learning

The basic idea of advantage-learning or A-learning, is to estimate treatment effect $\tau (\mathbf{x})$ directly, rather than estimating the reponse surfaces like Q-learning which can be deemed as an indirect way. Here are some selected approaches in this camp.

Individualized Treatment Rule

Most of the previous methods target on CATE estimation, then the optimal treatment can be found automatically by looking at the estimates. Whereas, obtaining the consistent and unbiased estimates for CATE is difficult in practice. However, if we only need to know which treatment is the optimal, we actually do not need to estimate the CATE but only the treatment assignment rule, or so-called individual-treatment-rule (ITR), $d(\mathbf{x})=sgn\{\tau (\mathbf{x})\}$, which is a mapping $d:\mathcal{X}\to T$. In practice, we do not need to estimate $\tau (\cdot )$ but $d(\cdot )$ directly. See Figure (2) for a demonstration of the difference between Q-learning, A-learning, and ITR.

Figure 2: Illustration of causal inference approaches. (A) Q-learning aims at estimating the response surfaces of each treatment group; (B) A-learning aims at estimation of treatment effect directly; (C) ITR aims at estimation of treatment region rather than the treatment effect.

Recall the ITR is indeed a mapping from the covariate space into the treatment space $d:\mathcal{X}\to T$. According to Qian and Murphy (2011) and Zhao et al. (2012), the value function under the ITR $d(\mathbf{X})$ is $$V(d)=\mathbb{E}[Y\mid d(\mathbf{X})=t]=\mathbb{E}\left[\frac{1[T=d(\mathbf{X})]Y}{\pi (\mathbf{X}\mid T=t)}\right]$$ and the optimal ITR is the one corresponds to the largest value function $$d^{\ast }()=\arg \underset{d()}{max}V(d).$$ Notice that solving this objective function requires discontinuous loss functions, which may incur certain numerical difficulty. So far, various methods have been proposed in the literature to obtain the ITR following this idea (Qian and Murphy 2011; Zhao et al. 2012; Xu et al. 2015; Zhu et al. 2017; X. Zhou et al. 2017; Zhang and Zhang 2018; Qi et al. 2020). Due to its structure similar to propensity score weighting but only on outcome $Y$, it is also known as outcome weighted learning, or O-learning (Zhao et al. 2012).

Other tree-based methods

There are many tree-based methods such as virtual trees (Foster 2013), causal tree/forest (Athey and Imbens 2016; Wager and Athey 2018), causal boosting/PTO forest/causal MARS (Powers et al. 2018). Actually, most of them can be categorized as S- or T-learners on a high-level.

(For more details on CATE/HTE estimation, please refer to another post Meta-Learners)