```{r set-options, echo=FALSE, cache=FALSE} options(width = 1000) knitr::opts_chunk$set(fig.width=7, fig.height=5, out.width = "70%", fig.align = 'center') knitr::opts_chunk$set(class.source = "fold-hide") knitr::opts_chunk$set(collapse=TRUE) knitr::opts_chunk$set(message = FALSE) ``` ## Overview As we discussed previously, when the treatment assignment is not independent of the potential outcomes, the difference-in-means estimator can be biased. In this section, we will discuss the propensity score and the inverse probability weighting (IPW) method to adjust for this bias. ## Example: Simpson’s Paradox Suppose we split the population into two subgroups, say those with age $\leq 50$ and those with age $> 50$. We are interested in the effect of a new treatment on recovery from a disease. Within each group, we randomly assign treatment, but with different probabilities. The treatment is binary, and the outcome indicates whether the patient recovered from the disease. Hence, the treatment effect is the difference in recovery probabilities. We are interested in the average treatment effect (ATE) for the entire population. Here are the data we observed: | Group / Treatment | Treatment 1 | Treatment 0 | |---------------|---------------------|---------------------| | Age $\leq 50$ | __18 out of 20 (90\%)__ | 150 out of 180 (83.3\%) | | Age $> 50$ | __50 out of 80 (62.5\%)__ | 10 out of 20 (50\%) | | Pooled | 68 out of 100 (68\%) | __160 out of 200 (80\%)__ | We can see that the treatment effect is positive (higher probability of recovery) in both groups, but negative (worse than control) overall. This is known as [Simpson's paradox](https://en.wikipedia.org/wiki/Simpson%27s_paradox). __Even though we used random treatment assignment within both groups__, the pooled difference-in-means estimator is biased for the overall ATE. Which assumption in our previous lecture is violated? * The unconditional random treatment assignment assumption is violated in this case. Treatment is independent of the potential outcomes within each age group, but not in the pooled population. In the younger group, the potential outcomes are higher overall, but the treatment assignment is more likely to be 0. In the older group, the potential outcomes are lower, but the treatment assignment is more likely to be 1. Although this dependency may not be intentional, it is not uncommon in practice. If we know this ahead of time, we could estimate the average treatment effect separately in each group, where the chance of treatment assignment remains constant. We can then average the effects across groups based on their respective sample sizes. In this case, the estimated ATE would be: \[ \begin{aligned} \tau &= \frac{200}{300} \times \left(\frac{18}{20} - \frac{150}{180}\right) + \frac{100}{300} \times \left(\frac{50}{80} - \frac{10}{20}\right) \\ & \approx 8.63\% \end{aligned} \] ## A Different View Let's consider a different view of the same formula given above. Let's group the two parts associated with treatment 1 as \[ \begin{aligned} & \frac{200}{300} \times \frac{18}{20} + \frac{100}{300} \times \frac{50}{80} \\ =& \frac{1}{300} \times \left( \frac{18}{20 / 200} + \frac{50}{80 / 100} \right) \end{aligned} \] Let $G_i \in \{1,2\}$ denote the group membership. The expression above becomes: \[ \frac{1}{n} \left( \sum_{i: G_i=1} \frac{A_i Y_i}{ \widehat \Pr(A_i = 1| G_i=1)} + \sum_{i: G_i=2} \frac{A_i Y_i}{ \widehat \Pr(A_i = 1| G_i=2)} \right) \] If we further encode the group label as a covariate $X_i$, where $X_i = 1$ if the patient is in group 1 and $X_i = 0$ if the patient is in group 2, then the above formula can be written as \[ \frac{1}{n} \sum_{i = 1}^{n} \frac{A_i Y_i}{ \Pr(A_i = 1| X_i )} \] And similarly, we can write the part associated with treatment 0, and combine the two as an estimator of the ATE: \[ \widehat\tau_{\text{ipw}} = \frac{1}{n} \sum_{i = 1}^{n} \left( \frac{A_i Y_i}{ \widehat e(X_i) } - \frac{(1 - A_i) Y_i}{ 1 - \widehat e(X_i) } \right) \] where $e(X_i) = \Pr(A_i = 1| X_i )$ is the probability of receiving treatment given the covariate $X_i$. $e(X_i)$ is also known as the __propensity score__. Compared with the difference-in-means estimator, the above formula is a weighted average estimator, with each observation weighted by the inverse probability of receiving its observed treatment label so that the sum scales back to the total of $n$ observations. This technique is known as __inverse propensity weighting__ (IPW). The estimator was originally proposed by @horvitz1952generalization in survey sampling and then used in causal inference. ## Assumptions Recall that our ultimate goal is to estimate the __average treatment effect__ (ATE), following the potential outcome notation: \[ \tau = \E[Y_i(1) - Y_i(0)] \] Our goal here is to establish some properties of the new estimator $\widehat\tau_{\text{ipw}}$. To do so, we need to establish some assumptions. Previously, we had two assumptions, Independence and SUTVA. It is obvious that the independence assumption is no longer valid. But some new assumptions could help us to establish the validity of the IPW estimator, besides the SUTVA assumption. > Unconfoundedness: $A_i \perp \{Y_i(0), Y_i(1)\} | X_i$ Essentially this is the same as the independence assumption, but now we condition on the covariate $X_i$. Why do we call it unconfoundedness? Because conditioning on $X_i$ is sufficient to eliminate the treatment-selection bias in our estimation. The assumption has many different names, such as __conditional ignorability__ or __conditional independence__. Together with positivity, it is also called __strong ignorability__. Readers are referred to [@rubin1978bayesian] and [@ding2023first] for more details. The difference is rather technical and is not the focus of this lecture. > Positivity: The treatment assignment is not deterministic for any value of $X$, i.e., $0 < P(A = 1 | X = x) \doteq e(x) < 1$. This is to ensure that we can, in principle, compare both treatment labels at each value of $X$. For stable estimation, we often use the stronger overlap condition $\eta \leq e(x) \leq 1-\eta$ for some $\eta>0$. This may seem to be pretty natural, but when the cardinality of $X$ is large, it could be violated. ## Inverse-propensity Weighting Because estimating $\widehat e(X_i)$ comes with additional randomness, and could potentially harm the performance of the estimator, let's first consider a version with $e(X_i)$ known (e.g., by design). Then the estimator becomes \[ \widehat\tau_\text{ipw}^\ast = \frac{1}{n} \sum_{i = 1}^{n} \left( \frac{A_i Y_i}{ e(X_i) } - \frac{(1 - A_i) Y_i}{ 1 - e(X_i) } \right) \] First, we can see that this estimator is an unbiased estimator of the ATE: \[ \begin{aligned} \E[ \widehat\tau_\text{ipw}^\ast ] =& \E \left[ \frac{A_i Y_i}{ e(X_i) } - \frac{(1 - A_i) Y_i}{ 1 - e(X_i) } \right] \quad \text{by IID} \\ =& \E \left[ \frac{A_i Y_i(1)}{ e(X_i) } - \frac{(1 - A_i) Y_i(0)}{ 1 - e(X_i)} \right] \quad \text{by SUTVA} \\ =& \E \left[ \E \left[ \frac{A_i Y_i(1)}{ e(X_i) } \Biggm| X_i \right] - \E \left[ \frac{(1 - A_i) Y_i(0)}{ 1 - e(X_i)} \Biggm| X_i \right] \right] \\ =& \E \left[ \frac{e(X_i)}{ e(X_i) } \E[ Y_i(1) | X_i ] - \frac{1 - e(X_i)}{1- e(X_i) } \E [Y_i(0) | X_i ] \right] \quad \text{by Unconfoundedness} \\ =& \E \, [ Y_i(1) - Y_i(0) ] \\ \end{aligned} \] Now, to bridge the gap between the known $e(X_i)$ and the estimated $\widehat e(X_i)$, we can further assume that * \( |Y_i| \leq M \) * \( \sup_{x} \left| e(x) - \widehat e(x) \right| \rightarrow 0 \) Under the overlap condition, it is easy to see that \( \left| \widehat\tau_\text{ipw}^\ast - \widehat\tau_\text{ipw} \right| \) also converges to 0 as \( \widehat e(x) \) converges uniformly to \( e(x) \). Furthermore, \( \widehat\tau_\text{ipw}^\ast \) is an average of i.i.d. terms and is therefore consistent. Hence \( \widehat\tau_\text{ipw} \) is also consistent.[^estimated] Obtaining a $\sqrt{n}$ rate for the estimator with estimated propensity scores requires additional rate and regularity conditions. ## The Density Ratio View The idea of propensity score weighting is essentially manipulating the distribution of the covariate $X$ in the treated and control groups so that they resemble the overall distribution of $X$. To see this, let's consider having two samples * Treated group: $\{X_i\}_{i: A_i = 1}$, with density $f_1(x)$ * Control group: $\{X_i\}_{i: A_i = 0}$, with density $f_0(x)$ Then the overall distribution of $X$ can be written as a mixture of the two groups: \[ f(x) = \pi_1 f_1(x) + \pi_0 f_0(x) \] where $\pi_1$ and $\pi_0$ are the proportions of treated and control subjects in the overall population. Note that using the Bayes identity, we have \[ \begin{aligned} f_1(x) &= f(x | A = 1) = \frac{e(x) f(x)}{\Pr(A = 1)} \\ \rightarrow \quad f(x) &= \frac{\Pr(A = 1)}{e(x)} f_1(x) \end{aligned} \] Now, if we look at the expected potential outcome under treatment for the entire population, we have \[ \begin{aligned} E_{f}[Y(1)] &= \int E[Y(1) | X = x] f(x) dx \\ &= \int E[Y(1) | X = x] \frac{\Pr(A = 1)}{e(x)} f(x | A = 1) dx \\ &= \E_{X \sim f_1} \left[ E[Y(1) | X] \frac{\Pr(A = 1)}{e(X)} \right] \\ &= \E_{X \sim f_1} \left[ E[Y(1) | X, A = 1] \frac{\Pr(A = 1)}{e(X)} \right] \\ &= \E \left[ \frac{A Y}{e(X)} \right] \quad \text{by unconfoundedness and SUTVA} \end{aligned} \] Hence a natural estimator of $E_f[Y(1)]$ is \[ \widehat E_f[Y(1)] = \frac{1}{n} \sum_{i = 1}^{n} \frac{A_i Y_i}{e(X_i)} \] This means that we are reweighting the treated group by the density ratio $\Pr(A = 1) / e(X)$ to recover the overall average potential outcome under treatment. Similarly, we can reweight the control group by $\Pr(A = 0) / (1 - e(X))$ to recover the overall average potential outcome under control. Combining the two, we arrive at the IPW estimator again.[^radon] ## Examples ### Example 1: Direct Implementation The `lalonde` dataset in the `Matching` package is a popular dataset for causal effect estimation. It contains the following variables: * `treat`: a binary variable indicating whether the subject received the treatment * `age`: the age of the subject * `educ`: the education level of the subject * `black`: a binary variable indicating whether the subject is black * `hisp`: a binary variable indicating whether the subject is Hispanic * `married`: a binary variable indicating whether the subject is married * `nodegr`: a binary variable indicating whether the subject has a degree * `re74`: the real earnings in 1974 * `re75`: the real earnings in 1975 * `re78`: the real earnings in 1978 LaLonde (1986) was interested in the causal effect of a job training program on earnings. These are experimental data, meaning that the participants were randomly assigned to the treatment and control groups. Propensity adjustment is therefore not required for unbiasedness here, but we can still use this dataset to illustrate the IPW calculation. If we consider a naive difference-in-means estimator, we would have: ```{r} library(Matching) data(lalonde) head(lalonde) # Calculate the average earnings for treated and control groups avg_treated <- mean(lalonde$re78[lalonde$treat == 1]) avg_control <- mean(lalonde$re78[lalonde$treat == 0]) # Calculate the difference in means ate_diff_means <- avg_treated - avg_control ate_diff_means ``` We can now apply the IPW estimator to illustrate how the weighting calculation works. ```{r} ps_model <- glm(treat ~ age + educ + black + hisp + married + nodegr + re74 + re75, family = binomial(link = "logit"), data = lalonde) # Extract the propensity scores lalonde$pscore <- predict(ps_model, type = "response") lalonde$psweight <- ifelse(lalonde$treat == 1, 1 / lalonde$pscore, 1 / (1 - lalonde$pscore)) mean(lalonde$treat*lalonde$re78*lalonde$psweight - (1 - lalonde$treat)*lalonde$re78*lalonde$psweight) ``` ### Example 2: Propensity Score Stratification An interesting fact is that we do not need to know the full information of $X$ to satisfy the conditional independence assumption. We only need to know the propensity score $e(X)$. Formally, this means that \[ \text{If} \quad A_i \perp \{Y_i(0), Y_i(1)\} | X_i \quad \text{then} \quad A_i \perp \{Y_i(0), Y_i(1)\} | e(X_i) \] The proof is relatively straightforward, by showing that the conditional distribution of $A_i$ given $e(X_i)$, $Y_i(1)$ and $Y_i(0)$ is the same as the conditional distribution of $A_i$ given $e(X_i)$. This leads to a new idea: we can first fit a model to estimate the propensity score and then stratify the data based on the estimated score. We can estimate the ATE within each stratum and average the stratum-specific ATEs to obtain the overall ATE. This is known as the __propensity score stratification__ method. Compared with IPW, stratification is often less sensitive to the exact magnitudes of the propensity scores because they are used to form groups. However, accurate rankings alone are not sufficient: treated and control covariates should still be adequately balanced within each stratum. There is a natural trade-off here caused by the number of strata. If we have too many strata, then the variance of the estimator will be large. If we have too few strata, then the bias of the estimator will be large. The following code shows how to implement the method with the `lalonde` dataset. We use an ad hoc choice of 7 strata, which ensures a good number of observations in each stratum. ```{r} lalonde$pscore_strata <- cut(lalonde$pscore, breaks = 7) lalonde$pscore_strata <- as.numeric(lalonde$pscore_strata) # Estimate the ATE within each stratum ate_strata <- sapply(split(lalonde, lalonde$pscore_strata), function(dat) mean(dat$re78[dat$treat == 1]) - mean(dat$re78[dat$treat == 0])) means1_within_strata <- aggregate(re78 ~ pscore_strata, data = lalonde[lalonde$treat == 1,], FUN = mean) means0_within_strata <- aggregate(re78 ~ pscore_strata, data = lalonde[lalonde$treat == 0,], FUN = mean) cbind('strata' = means1_within_strata$pscore_strata, 'treated' = means1_within_strata$re78, 'control' = means0_within_strata$re78, 'strata size' = table(lalonde$pscore_strata)) sum((means1_within_strata$re78 - means0_within_strata$re78)* table(lalonde$pscore_strata)) / nrow(lalonde) ``` ### Example 3: Matching and the Average Treatment Effect of the Treated Up to now, we have not explicitly discussed the population over which the treatment effect is being defined. In the previous examples, we implicitly estimated the average treatment effect (ATE) for the entire population that generated our observed sample of covariates $X_i$'s. In that setting, the sample is assumed to be randomly drawn from a single population that includes both treated and untreated individuals. However, in some studies this assumption does not hold, and the researcher’s interest may lie in a more specific population. For example, treated subjects may come from a particular hospital, while control subjects are obtained from an external public database. In such cases, the goal is often to estimate the treatment effect for the treated-source population, i.e., the population represented by patients in that hospital, rather than for the combined population. This is known as the average treatment effect of the treated (ATT). Formally, our target of interest is (for whatever implicit population that $A_i = 1$ refers to) \[ \tau_{\text{ATT}} = \E[Y_i(1) - Y_i(0) | A_i = 1] \] Note that we can further condition on the propensity score $e(X_i)$ to get independence between the treatment assignment $A$ and the potential outcomes: \[ \begin{aligned} \tau_{\text{ATT}} &= \E\left[ \E[Y_i(1) - Y_i(0) \mid e(X_i), A_i = 1] \mid A_i = 1\right] \\ &= \E\left[ \E[Y_i(1) \mid e(X_i), A_i = 1] - \E[Y_i(0) \mid e(X_i), A_i = 1] \mid A_i = 1\right] \\ &= \E\left[ \E[Y_i(1) \mid e(X_i)] - \E[Y_i(0) \mid e(X_i)] \mid A_i = 1\right] \quad \text{(by unconfoundedness)} \\ &= \E\left[ \E[Y_i(1) \mid e(X_i), A_i = 1] - \E[Y_i(0) \mid e(X_i), A_i = 0] \mid A_i = 1\right] \\ &= \E\left[ \E[Y_i \mid e(X_i), A_i = 1] - \E[Y_i \mid e(X_i), A_i = 0] \mid A_i = 1\right] \quad \text{(by SUTVA)} \\ \end{aligned} \] Note that the first term can be estimated by directly averaging all treated subjects, while the second term, for subject $i$, can be estimated by (locally) averaging all control subjects with the same propensity score $e(X_i)$. Thus, an estimator of ATT is \[ \widehat\tau_{\text{ATT}} = \frac{1}{n_T} \sum_{i: A_i = 1} \left( Y_i - \widehat E[Y \mid e(X_i), A = 0] \right) \] We can use the same stratification idea by binning the propensity scores and treating the subjects within each stratum as having the same propensity score[^covmath]. Alternatively, we can match each treated subject with the closest control subject, or a set of the closest $M$ control subjects, based on the propensity score $e(X_i)$. The difference for each matched set can then be viewed as an estimate of the treatment effect given $e(X_i)$. The following code is an implementation of this idea with the `lalonde` dataset. ```{r} matching_control = data.frame(pscore = lalonde$pscore[lalonde$treat == 0], outcome = lalonde$re78[lalonde$treat == 0]) matching_treat = data.frame(pscore = lalonde$pscore[lalonde$treat == 1], outcome = lalonde$re78[lalonde$treat == 1]) ntreat = nrow(matching_treat) matched_outcome = rep(NA, ntreat) # number of matches M = 5 for (i in 1:ntreat) { abs_diffs = abs(matching_treat$pscore[i] - matching_control$pscore) idmatch = order(abs_diffs)[1:M] matched_outcome[i] = mean(matching_control$outcome[idmatch]) } mean(matching_treat$outcome - matched_outcome) ``` We can see that the estimated ATT is different from the ATE estimate. The ATT focuses on the treated population rather than the entire population represented by the data, so these two targets need not be equal. *** [^estimated]: In some regular parametric settings, using an estimated propensity score can reduce the asymptotic variance relative to directly using the true propensity score because propensity-score estimation induces useful calibration. This improvement is not universal and depends on the propensity model and estimator. [^radon]: This view is closely related to Radon-Nikodym theorem, which states that \[ \mathbb{E}_P[g(X)] = \mathbb{E}_Q\left[g(X) \frac{dP}{dQ}(X)\right] \] where $P$ and $Q$ are two probability measures, and $\frac{dP}{dQ}(X)$ is the Radon-Nikodym derivative (or density ratio) of $P$ with respect to $Q$. Here, we can view the overall distribution of $X$ as $P$, and the distribution of $X$ in the treated group as $Q$. Then the density ratio is exactly $\Pr(A = 1) / e(X)$. This idea will also be used later in the outcome weighted learning method for personalized medicine. [^covmath]: Another alternative we can consider is the covariate matching, which matches the treated and control subjects based on their covariates $X_i$ directly. The idea is similar, and we leave it as an exercise for the readers.