```{r set-options, echo=FALSE, cache=FALSE}
  options(width = 1000)
  knitr::opts_chunk$set(fig.width=7, fig.height=5, out.width = "60%", fig.align = 'center')
  knitr::opts_chunk$set(class.source = "fold-hide")
  knitr::opts_chunk$set(collapse=TRUE)
  knitr::opts_chunk$set(message = FALSE)
```

## Overview

For the second half of this course, the central theme is __decision making__. This is a very broad topic and has many applications. Two particular settings we will focus on are __personalized medicine__ [@kosorok2019precision] and __reinforcement learning__ [@sutton2018reinforcement]. The two settings have many overlaps, but we will mostly discuss single-time-point decision making for personalized medicine and sequential decision making under the Markov decision process (MDP) framework for reinforcement learning. However, in a broader view, many practical problems in personalized medicine are also sequential decision-making problems and therefore require reinforcement learning tools to solve.

For both settings, we will put them into a __causal inference__ framework. Causal inference is a branch of statistics that focuses on understanding the cause-and-effect relationships between variables. Naturally, we view the treatment decision as an action that needs to be understood or optimized, and the health outcome is something we want to improve. Hence the essential question is to understand what would happen if we take a certain action (treatment) on a certain individual. This is often referred to as the __counterfactual outcome__. Besides, decision-making problems are frequently encountered in many fields, ranging from economics and marketing to cognitive science, environmental science and sports analytics. Essentially, anytime a decision needs to be made based on individual characteristics, such models can be useful.


## Causal Treatment Effects

Before finding the best treatment for an individual, one fundamental concept is understanding the treatment effect. Let's define some notation. Suppose we have a binary treatment label $A \in \{0, 1\}$, for which 1 indicates a type of treatment and 0 indicates another treatment or simply the placebo control. We also have a response variable $Y$, which is the outcome of interest. For an individual $i$, there are potentially two outcomes, $Y_i(0)$ and $Y_i(1)$, and the difference between the two is the __causal treatment effect__ ^[In this notation, we usually consider each $\Delta_i$ as a random variable. They may or may not have the same mean value across different individuals (e.g., if they depend on additional covariates), but that will not change the main result of this section. In our later lectures, we could consider using additional covariates $X$ to improve our estimation if they are available.]:

\[
\Delta_i = Y_i(1) - Y_i(0)
\]

The fundamental difficulty in estimating the treatment effect is that we can only observe one of the two outcomes. For example, suppose you receive offers from two universities, UIUC and Harvard, and can only choose one of them. If you choose UIUC and observe your salary after graduation, then you will never observe the outcome of going to Harvard. The two outcomes essentially belong to two parallel universes. This dilemma between the __actual__ (or realized) and __counterfactual__ outcomes is known as the __fundamental problem of causal inference__ by @holland1986statistics. However, statistical frameworks allow us to estimate related treatment-effect quantities. The most commonly accepted view is the __potential outcome__ framework, which originated in @neyman1923application, who considered a randomized experiment, and @rubin1974estimating, who considered observational (non-randomized) studies. In this section, let us explore this framework using the randomized control trial and discuss some of its crucial assumptions.

## Randomized Control Trial

Let's consider the __average treatment effect__ (ATE), which is the treatment effect averaged over the entire population:

\[
\tau = \E[\Delta_i] = \E[Y_i(1)] - \E[Y_i(0)]
\]

To entertain the idea of potential outcomes, let's suppose that we could actually visit both parallel universes. Then we would observe both $Y_i(1)$ and $Y_i(0)$. For the observed sample, a natural quantity would be

\[
\tau_n = \frac{1}{n} \sum_{i=1}^n \Delta_i = \frac{1}{n} \sum_{i=1}^n \big( Y_i(1) - Y_i(0) \big) (\#eq:sate)
\]

This quantity is called the __sample average treatment effect__ (SATE). It is the finite-sample treatment-effect target. Under i.i.d. sampling from the population, $\E(\tau_n) = \tau$, so the SATE is also unbiased for the population ATE. However, we cannot observe $\tau_n$ because one potential outcome is missing for every individual.

In reality, we could randomly assign the treatment to a group of individuals and compare the outcomes between the two groups. This is known as the __randomized control trial__ (RCT). In this case, we only observe one of the two potential outcomes for each individual:

\[
\begin{aligned}
Y_i &= A_i Y_i(1) + (1 - A_i) Y_i(0)\\
&= \begin{cases}
Y_i(1) & \text{if } A_i = 1 \\
Y_i(0) & \text{if } A_i = 0
\end{cases}
\end{aligned}
\]

where $A_i$ is the treatment assignment for individual $i$.

## Difference-in-Means Estimator

A naive (but pretty good) idea is to estimate the treatment effect using the mean difference between the two groups, which is called the __difference-in-means estimator__:

\[
\begin{aligned}
\widehat\tau &= \frac{1}{n_1} \sum_{A_i = 1} Y_i - \frac{1}{n_0} \sum_{A_i = 0} Y_i \\
&= \frac{1}{n_1} \sum_{i = 1}^n A_i Y_i - \frac{1}{n_0} \sum_{i = 1}^n (1 - A_i) Y_i
\end{aligned}
\]

where $n_1$ and $n_0$ are the sample sizes of the two groups. A natural question is whether this estimator is a "good" estimator. The answer is, only in some restricted settings. Here is a numerical example, same as the one in section \@ref(discussion1) later. Suppose the potential outcomes are generated from a nonlinear model:

```{r fig.width=8, fig.height=6}
  # setting parameters
  n <- 200
  nsim = 100
  
  set.seed(1)
  
  tauhat <- rep(NA, nsim)
  sate <- rep(NA, nsim)
  
  for (i in 1:nsim) {
    
    # generate covariates
    X <- rnorm(n)
    A <- rbinom(n, 1, exp(X) / (1 + exp(X)))
    Y <- rnorm(n, mean = 0.5*X + A*X^2)
    
    tauhat[i] <- mean(Y[A == 1]) - mean(Y[A == 0])
    
    sate[i] <- mean(X^2)
  }
  
  # barplot of the two estimators
  df <- data.frame(
    Method = rep(c("Difference-in-means", "SATE"), each = nsim),
    Estimate = c(tauhat, sate)
  )
  
  boxplot(Estimate ~ Method, data = df, 
          ylab = "Estimated ATE", xlab = "Method")
  
```

Hence, this estimator can be biased without random treatment assignment. However, under complete randomization, it is unbiased for the SATE and, under i.i.d. sampling, for the ATE. This is the main result of this section. The key advantage of RCT is that this difference-in-means estimator is __unbiased even when there are other covariates__ that affect the outcome. We will discuss the utilization of additional covariates later.

## Assumptions

The unbiasedness of this estimator can be established by connecting it with the SATE in \@ref(eq:sate). Complete randomization allows us to target the SATE for a fixed sample, while i.i.d. sampling connects the SATE to the population ATE. We use two important assumptions:

> Independence: $A_i \perp \{Y_i(0), Y_i(1)\}$

This is a pretty natural assumption, which suggests that the assignment of treatment has nothing to do with all the possible potential outcomes. What situation would violate this assumption? For example, a patient may choose what he/she believes to be the better treatment. In this case, the treatment assignment is not independent of the potential outcomes.

> SUTVA (Stable Unit Treatment Values Assumption): $Y_i = Y_i(A_i)$

This assumption is a bit subtle. It has two components.

  * __No interference__: the outcome of individual $i$ is not affected by other individuals. 
  * __Consistency__: There are no hidden forms of treatment, i.e., the outcome of individual $i$ is the same as the outcome of individual $i$ under the same treatment assignment.

The first part is relatively easy to understand. For example, if you are a patient in a clinical trial, your outcome should not be affected by the treatment of other patients. However, if all patients share and are competing on the same medical resources, then the outcome of one patient may be affected by the treatment of others. The second part is somewhat philosophical and it can be closely related to how we define the treatment. One could simply say that this is automatically satisfied. But if the patient does not take the treatment as prescribed, then the outcome may be different. In that case, there is an entire area of research on __noncompliance__ and __intention-to-treat__ [@gupta2011intention] analysis. 

## Unbiasedness

As we discussed previously, we could consider the difference-in-means estimator. In fact, under complete randomization with fixed $n_1$, this estimator is unbiased for the SATE. Under i.i.d. sampling from the population, it is also unbiased for the ATE. To see this, let's consider the conditional expectation of the estimator given the potential outcomes and the number assigned to treatment^[We condition on the potential outcomes because this avoids making assumptions about how they are generated. Random assignment then determines the expectation over treatment labels.]:

\[
\begin{aligned}
& \, \E \left[ \frac{1}{n_1} \sum_{i = 1}^n A_i Y_i \Biggm| \{Y_i(0), Y_i(1)\}_{i = 1}^n, n_1 \right] \\
=& \, \E \left[ \frac{1}{n_1} \sum_{i = 1}^n A_i Y_i(1) \Biggm| \{Y_i(0), Y_i(1)\}_{i = 1}^n, n_1 \right] \quad \text{by SUTVA} \\
=& \, \frac{1}{n_1} \sum_{i = 1}^n Y_i(1) \, \E \Big[ A_i \Bigm| \{Y_i(0), Y_i(1)\}_{i = 1}^n, n_1 \Big] \quad Y_i(1)'s \text{ are constant} \\
=& \, \frac{1}{n_1} \sum_{i = 1}^n Y_i(1) \frac{n_1}{n} \quad \text{by complete randomization} \\
=& \, \frac{1}{n} \sum_{i = 1}^n Y_i(1) \\
\end{aligned}
\]

Here, complete randomization gives $\E[A_i \mid \{Y_i(0),Y_i(1)\}_{i=1}^n,n_1] = n_1/n$ for every individual. The same argument can be applied to the second term. Therefore, the difference-in-means estimator is conditionally unbiased for the SATE. Averaging over i.i.d. samples then gives unbiasedness for the ATE.

## Asymptotic Inference

We can then analyze the variance of this estimator for the population ATE under i.i.d. sampling, conditional on the sample size of both arms.^[For finite-sample randomization inference about the SATE, the exact variance also contains a term involving the finite-sample variance of $Y_i(1)-Y_i(0)$. Here we use the superpopulation perspective for the ATE.]

\[
\Var(\widehat \tau | n_0, n_1) = \frac{1}{n_1} \Var\big[Y_i(1)\big] + \frac{1}{n_0} \Var\big[Y_i(0)\big] (\#eq:var)
\]

By the central limit theorem, $\widehat \tau$ is asymptotically normal:

\[
\frac{ \widehat \tau - \tau }{ \sqrt{\sigma_1^2/n_1 + \sigma_0^2/n_0} } \xrightarrow{d} \cN \left(0, 1\right)
\]

where $\sigma_1^2 = \Var\big[Y_i(1)\big]$ and $\sigma_0^2 = \Var\big[Y_i(0)\big]$. Sample version estimates of these two quantities are:

\[
\begin{aligned}
\widehat \sigma_1^2 &= \frac{1}{n_1 - 1} \sum_{A_i = 1} \left( Y_i - \bar Y_1\right)^2 \\
\widehat \sigma_0^2 &= \frac{1}{n_0 - 1} \sum_{A_i = 0} \left( Y_i - \bar Y_0\right)^2
\end{aligned}
\]

where $\bar Y_1 = \frac{1}{n_1} \sum_{A_i = 1} Y_i$ and $\bar Y_0 = \frac{1}{n_0} \sum_{A_i = 0} Y_i$ are the sample means of the two groups, respectively. Let $\widehat V = \widehat\sigma_1^2/n_1 + \widehat\sigma_0^2/n_0$. We may construct the confidence interval for the ATE as:

\[
\widehat \tau \pm z_{1 - \alpha / 2} \sqrt{\widehat V}
\]

where $z_{1 - \alpha / 2}$ is the $(1 - \alpha / 2)$-th quantile of the standard normal distribution. You can also use the $t$ distribution if the sample size is relatively small, but you need to calculate the appropriate standard error and degrees of freedom.

## Numerical Example

Let's use a simulation study to illustrate the properties of the difference-in-means estimator. Suppose we have $n = 200$ patients, and we randomly assign $n_1 = 100$ patients to the treatment group and $n_0 = 100$ patients to the control group. Suppose our outcomes are generated from a nonlinear model

\[
\E(Y | X = x, A = a) = 0.5 \times x + a \times x^2. (\#eq:model)
\]

In this case, the potential outcomes for a subject with covariate value $x$ are $0.5\times x + x^2$ for treatment 1 and $0.5\times x$ for treatment 0. If $X$ follows a standard normal distribution, the expected treatment effect is $\tau = E(X^2) = 1$. In a regression problem, we typically want to observe $X$ and model the relationship between $X$ and $Y$ and then infer the treatment effect. However, in an RCT, this is not necessary. Let's generate the observed outcomes and estimate the ATE using the difference-in-means estimator without using $X$. We repeat this process 500 times and examine the point estimates and their confidence intervals.

```{r fig.width=12, fig.height=5}
  # setting parameters
  n <- 200
  n1 <- 100
  n0 <- 100
  tau <- 1 # true ATE
  nsim = 500
  
  set.seed(1)
  
  tauhat <- rep(NA, nsim)
  tauhatsd <- rep(NA, nsim)
  
  for (i in 1:nsim) {
    # generate potential outcomes
    X <- rnorm(n)
    Y1 <- rnorm(n, mean = 0.5*X+X^2)
    Y0 <- rnorm(n, mean = 0.5*X)
    
    # treatment label
    A = sample(c(rep(1, n1), rep(0, n0)))
    
    # observed outcomes
    Y <- A * Y1 + (1 - A) * Y0
    
    tauhat[i] <- mean(Y[A == 1]) - mean(Y[A == 0])
    tauhatsd[i] <- sqrt(var(Y[A == 1]) / n1 + var(Y[A == 0]) / n0)
  }

  # make two plots on the same row
  par(mfrow = c(1, 2))
  
  # set margin of figure
  par(mar = c(4, 4, 1, 1))
  plot(tauhat[1:50], pch = 19, ylim = c(0, 2),
       xlab = "Simulation Runs", ylab = "Estimated ATE")
  abline(h = tau, col = "red")
  
  # adding confidence intervals
  for (i in 1:50) {
    ci_lower <- tauhat[i] - 1.96 * tauhatsd[i]
    ci_upper <- tauhat[i] + 1.96 * tauhatsd[i]
    arrows(x0 = i, y0 = ci_lower, x1 = i, y1 = ci_upper, angle = 90, code = 3, length = 0.05)
  }
  
  coverage = sum((tauhat - 1.96 * tauhatsd < tau) & (tauhat + 1.96 * tauhatsd > tau)) / nsim
  legend("topright", paste("Coverage probability of", nsim, "runs:", coverage))
  
  boxplot(tauhat, xlab = "Boxplot of Estimated ATE", ylab = "Estimated ATE")
  abline(h = tau, col = "red")
```

The result roughly demonstrates that the point estimation is unbiased and the confidence interval has the correct coverage probability. We should note that in this estimation, we __did not use the covariate information__. The unbiasedness comes purely from the randomization of the treatment assignment.

## Discussion

### Discussion 1: Observational Study {#discussion1}

In fact, the difference-in-means estimator is numerically the same statistic used in a two-sample comparison. Suppose we observe outcomes $Y_{11}, \ldots, Y_{1n_1}$ in one group and $Y_{01}, \ldots, Y_{0n_0}$ in the other, and we want to test whether the two groups have the same mean. The test statistic is the difference in means, while the asymptotic variance is the same as in \@ref(eq:var). Here we use the asymptotic variance rather than a $t$ distribution.

However, the way the data are generated determines whether this difference has a causal interpretation. In a randomized trial, the treatment label is independent of the potential outcomes. Therefore, if treatment has no effect, any systematic difference between the two groups cannot be attributed to pre-treatment differences. Without randomization, the two groups may have different outcome distributions even when the treatment has no causal effect. Hence, the difference-in-means estimator may not be valid for estimating the ATE. This is essentially an observational-study situation, which can be illustrated by the following example.

In a randomized trial design, our analysis shows that, even when we do not utilize the covariate information, we can still obtain an unbiased estimate of the average treatment effect. This is a very important property of RCT, and it is the reason why RCT is considered the gold standard for estimating the treatment effect. However, in many cases, the treatment cannot be assigned randomly for practical and ethical reasons. Such a scenario is an __observational study__. In that case, the difference-in-means estimator can be biased.^[A famous example concerns vitamin E supplements and the risk of coronary heart disease. Two papers in the New England Journal of Medicine in 1993 [@rimm1993vitamin; @stampfer1993vitamin] reported that vitamin E supplements were associated with a reduced risk of coronary heart disease in observational studies. However, later randomized controlled trials [@eidelman2004randomized; @lee2005vitamin] reported no significant effect.] Here is a simulation study. In this case, the treatment label is generated by a logistic model with covariate $X$, which makes higher $X$ values more likely to take treatment 1. Some correction will be needed to estimate the treatment effect in observational studies. We will discuss this in the next lecture.

```{r fig.width=12, fig.height=5}
  # setting parameters
  n <- 200
  # unlike the randomized trial above, the realized
  # group sizes vary across simulation runs
  tau <- 1 # true ATE
  nsim = 500
  
  set.seed(1)
  
  tauhat <- rep(NA, nsim)
  tauhatsd <- rep(NA, nsim)
  
  for (i in 1:nsim) {
    # generate potential outcomes
    X <- rnorm(n)
    Y1 <- rnorm(n, mean = 0.5*X+X^2)
    Y0 <- rnorm(n, mean = 0.5*X)
    
    # treatment label
    A = rbinom(n, 1, exp(X) / (1 + exp(X)))
    
    # observed outcomes
    Y <- A * Y1 + (1 - A) * Y0
    
    tauhat[i] <- mean(Y[A == 1]) - mean(Y[A == 0])
    tauhatsd[i] <- sqrt(var(Y[A == 1]) / sum(A == 1) + var(Y[A == 0]) / sum(A == 0))
  }

  # make two plots on the same row
  par(mfrow = c(1, 2))
  
  # set margin of figure
  par(mar = c(4, 4, 1, 1))
  plot(tauhat[1:50], pch = 19, ylim = c(0, 2),
       xlab = "Simulation Runs", ylab = "Estimated ATE")
  abline(h = tau, col = "red")
  
  # adding confidence intervals
  for (i in 1:50) {
    ci_lower <- tauhat[i] - 1.96 * tauhatsd[i]
    ci_upper <- tauhat[i] + 1.96 * tauhatsd[i]
    arrows(x0 = i, y0 = ci_lower, x1 = i, y1 = ci_upper, angle = 90, code = 3, length = 0.05)
  }
  
  coverage = sum((tauhat - 1.96 * tauhatsd < tau) & (tauhat + 1.96 * tauhatsd > tau)) / nsim
  legend("topright", paste("Coverage probability of", nsim, "runs:", coverage))
  
  boxplot(tauhat, xlab = "Boxplot of Estimated ATE", ylab = "Estimated ATE")
  abline(h = tau, col = "red")
```

In this lecture, we did not discuss the estimation of the treatment effect using covariates. Additionally, in many cases, we are not only interested in the average treatment effect over the entire population, but also in treatment effects that may differ across individuals. This is often referred to as the __heterogeneous treatment effect__. Estimating heterogeneous treatment effects is a more challenging problem and requires covariate information. Additional assumptions are also needed.

### Discussion 2: Regression Based Estimator

The following analysis is adapted from @wager2020stats. Let's consider a case where both potential outcomes are generated from a linear regression:

\[
Y_i(a) = c(a) + X_i \beta_{(a)} + \epsilon_i(a), \quad a = 0, 1
\]

where $c(a)$ is the intercept, $\beta_{(a)}$ is the coefficient, and $\epsilon_i(a)$ is the error term satisfying

\[
\E[ \epsilon_i(a) | X_i] = 0, \quad \Var[ \epsilon_i(a) | X_i ] = \sigma^2
\]

We can also assume, without loss of generality, that the randomized trial is balanced with $\Pr(A_i = 1) = \Pr(A_i = 0) = 0.5$ and that the covariate $X$ has mean 0 and variance 1. Under this model, the population ATE is $\tau = c(1)-c(0)$.

\[
\E[ X ] = 0 \quad \text{and} \quad \Var[ X ] = 1.
\]

If we analyze the difference-in-means estimator, we can see that

\[
\begin{aligned}
n \Var(\widehat \tau | n_0, n_1) &= 2 \Var[Y_i(1)] + 2 \Var[Y_i(0)] \\
&= 4 \sigma^2 + 2 \Var[X_i \beta_{(1)}] + 2\Var[X_i \beta_{(0)}] \\
&= 4 \sigma^2 + 2 \beta_{(1)}^2 + 2 \beta_{(0)}^2 \\
&= 4 \sigma^2 + ( \beta_{(1)} + \beta_{(0)} )^2 + ( \beta_{(1)} - \beta_{(0)} )^2  \\
\end{aligned}
\]

Now if we fit linear regression models to estimate these parameters, we would obtain $\widehat c(1)$ and $\widehat \beta_{(1)}$ from the samples that received treatment 1, and $\widehat c(0)$ and $\widehat \beta_{(0)}$ from the samples that received treatment 0. The regression-adjusted estimator can be written as

\[
\widehat \tau_\text{OLS} = \widehat c(1) - \widehat c(0) + \bar X \left( \widehat \beta_{(1)} - \widehat \beta_{(0)} \right)
\]

By analyzing the difference between $\widehat \tau_\text{OLS}$ and $\tau$, the average treatment effect, we have

\[
\begin{aligned}
\widehat \tau_\text{OLS} - \tau &= \left( \widehat c(1) - c(1) \right) - \left( \widehat c(0) - c(0) \right) + \bar X \left( \beta_{(1)} - \beta_{(0)} \right) \\
& \quad + \bar X \left( \widehat \beta_{(1)} - \beta_{(1)} \right) - \bar X \left( \widehat \beta_{(0)} - \beta_{(0)} \right)
\end{aligned} (\#eq:olserror)
\]

Keep in mind that for a regression model, the asymptotic variance of parameter estimates is \( \sigma^2 (Z^\T Z)^{-1} \) where $Z = \text{cbind}(1, X)$ is the design matrix. In our case, $X$ is independent of the intercept covariate, so \( (Z^\T Z)^{-1} \) is a diagonal matrix. In our case, the two diagonal terms are in the rate of $1/n_a$ since $X$ has variance 1.  

Back to \@ref(eq:olserror), the difference between the first two terms is asymptotically normal with variance $4\sigma^2/n$. The third term relies on the rate of $\bar X$ converging to 0, making it asymptotically normal with variance \( (\beta_{(1)} - \beta_{(0)})^2/n \). The last two terms converge to zero at a much faster rate because each is a product of \(\bar X\) and a small estimation error in $\beta$. Hence, overall we have

\[
\sqrt{n} (\widehat \tau_\text{OLS} - \tau ) \xrightarrow{d} \cN \left(0, 4\sigma^2 + (\beta_{(1)} - \beta_{(0)})^2 \right)
\]

Compared with the asymptotic variance of the difference-in-means estimator, the OLS estimator is more efficient because the \( ( \beta_{(1)} + \beta_{(0)} )^2 \) term disappears. Let's use a simulation to compare the two estimators.

```{r fig.width=12, fig.height=5}
  # setting parameters
  n <- 500
  n1 <- n/2
  n0 <- n/2
  tau <- 0.3 # true ATE  
  nsim = 500
  
  set.seed(1)
  
  tauhat_dim <- rep(NA, nsim)
  tauhat_ols <- rep(NA, nsim)
    
  for (i in 1:nsim) {
    # generate potential outcomes
    X <- rnorm(n)
    Y1 <- rnorm(n, mean = 0.5 + tau + X)
    Y0 <- rnorm(n, mean = 0.5 + 2*X) # you can change this coefficient to see the difference
    
    # treatment label
    A = sample(c(rep(1, n1), rep(0, n0)))
    
    # observed outcomes
    Y <- A * Y1 + (1 - A) * Y0
    
    # Difference-in-means estimator
    tauhat_dim[i] <- mean(Y[A == 1]) - mean(Y[A == 0])
    
    # OLS estimator
    X1 <- X[A == 1]
    Y_1 <- Y[A == 1]
    X0 <- X[A == 0]
    Y_0 <- Y[A == 0]
    
    model1 = lm(Y_1 ~ X1)
    model0 = lm(Y_0 ~ X0)
    
    tauhat_ols[i] <- coef(model1)[1] - coef(model0)[1] + mean(X) * (coef(model1)[2] - coef(model0)[2])
  }

  # compare the performance of the two estimators
  mean(tauhat_dim)
  mean(tauhat_ols)
  sd(tauhat_dim)
  sd(tauhat_ols)
```


***