<jupyter_start><jupyter_text>--- Lecture 2: Potential outcome model---  Introduction Given what we know from the introduction about the potential outcome model, we will initially prepare the Lalonde Dataset to fit the framework and use it as a running example going forward.What are this example's ...* potential outcomes* counterfactual state* treatment<jupyter_code>df = get_lalonde_data()
df.head()<jupyter_output><empty_output><jupyter_text>We are dealing with a binary treatment here: $D = 1$ if the individual did participate in the training program and $D = 0$ if  it did not. However, in practice assigning **treatment** is never that easy. We lump a lot of heterogeneity together (e.g. different sites, content of curriculum) that might affect the success of program participation. Maybe we should stratify the analysis by site? --- Individual-specific effect of treatment--- It would be great if we could get our hands on the individual-specific effect of treatment.\begin{align*}\delta_i = y_i^1 - y_i^0\end{align*}* Why do individuals have potentially different effects of treatment?<jupyter_code>fig, ax = plt.subplots()

x = np.linspace(-5, 5, 5000)
pdf = ss.norm.pdf(x, 0, 1)

ax.plot(x, pdf)

ax.set_xlabel(r"$\delta = Y^1 - Y^0$")
ax.set_ylabel("Density")
ax.set_xticklabels(["", "", "", 0.5, "", "", ""])
ax.set_xlim([-3, 3])
ax.set_ylim([0, 0.5])<jupyter_output><empty_output><jupyter_text>There might be considerable heterogeneity in the benefits of treatment among the population. And summarizing the distribution of benefits with a single number, for example $E[\delta]$, might result in a loss of information.**Examples*** medical treatment* .. Give our definitions of $(Y^1, Y^0, D)$ and their individual realizations $(y^1_i, y^0_i, d_i)$ we can now define the observed outcome $Y$ in terms of them.\begin{align*}Y = \begin{cases}Y^1 & \quad\text{if}\quad D = 1 \\Y^0 & \quad\text{if}\quad D = 0\end{cases}\end{align*}or more compactly in switching-regime notation\begin{align*}Y = D Y^1 + (1 - D) Y^0.\end{align*} This leads Holland (1986) to describe the fundamental problem of causal inference:$\rightarrow$ as only the diagonal of the table is observable we cannot simply compute $\delta_i$ by taking the difference in potential outcomes $(y^1_i, y^0_i)$.<jupyter_code>df.head()<jupyter_output><empty_output><jupyter_text>--- Population-level parameters--- It looks like we need to give up any hope of obtaining the individual-specific effect of treatment. But what can we still hope for?$\rightarrow$ population-level parameters* What are common examples?* What are the policy questions they address?* What is their relationship to each other? \begin{align*}\begin{array}{lll}E[Y^1 - Y^0] & \qquad ATE & \qquad\text{average effect of treatment} \\E[Y^1 - Y^0 \mid D = 1] & \qquad ATT & \qquad\text{average effect on treated} \\E[Y^1 - Y^0 \mid D = 0] & \qquad ATC & \qquad\text{average effect on control}\end{array}\end{align*}<jupyter_code>plot_individual_specific_effects(with_parameters=[0, 0.7, -0.5])
plot_individual_specific_effects(with_parameters=[0, -0.7, 0.5])
plot_individual_specific_effects(with_parameters=[0, 0, 0])<jupyter_output><empty_output><jupyter_text>--- Stable unit treatment value assumption--- The potential outcome model gets its empirical tractability when combined with the **Stable Unit Treatment Value Assumption (SUTVA)** of which there exist many formulations. We will go with the one from Imbens and Rubin (2015):> The potential outcomes for any unit do not vary with the treatments assigned to other units, and, for each unit there are no different forms or versions of each treatment level, which lead to different potential outcomes. The table below shows all possible assignment patterns for a hypothetical treatment where the only constraint is that at least one individual remains in the treatment and control group. As we increase participation from one to two individuals, the potential outcome $Y_1$ declines. * When do we need to expect this is violated?    * **influence patterns** that result from contact across individuals in social or physical space    * **dilution / concentration patterns** that one can assume would result from changes in the prevalence of treatment --- Treatment assignment and observational studies--- * randomized experiment\begin{align*}(Y^0, Y^1)  \perp \!\!\! \perp D\end{align*}* observational study> A *observational study* is an empirical investigation of treatments, policies, or exposures and the effects they cause, but it differs from an experiment in that the investigator cannot control the assignment of treatments to subjects. (Rosenbaum (2002)) --- Naive estimation of average causal effects--- We will now first outline the problem with the naive estimation of average causal effects. Then we take a closer look at the different sources of biases involved and finally discuss the set of assumptions used to ***solve*** these issues.\begin{align*}\hat{\delta}_{NAIVE} \equiv E_N[y_i\mid d_i = 1] - E_N[y_i\mid d_i = 0]\end{align*}We can further decompose the average treatment effect by treatment status as the individual assignment is mutually exclusive. \begin{align*}E[Y^1 - Y^0] & =  E[\delta] & = \{\pi E[Y^1\mid D = 1] + (1 - \pi)E[Y^1\mid D = 0]\} \\             &          & - \{\pi E[Y^0\mid D = 1] + (1 - \pi)E[Y^0\mid D = 0]\}\end{align*}The average treatment effect is a function of five unknowns. Which components can be easily computed from data?\begin{align*}E_N[y_i\mid d_i = 1] \xrightarrow{p} E[Y^1\mid D = 1] \neq E[Y^1]\\E_N[y_i\mid d_i = 0] \xrightarrow{p} E[Y^0\mid D = 0] \neq E[Y^0]\end{align*}  Biases\begin{align*}E[Y^1 \mid D = 1] - E[Y^0 \mid D = 0] & = E[\delta] + \underbrace{\{E[Y^0\mid D= 1] - E[Y^0\mid D= 0]\}}_{\text{Baseline bias}} \\                                      & + (1 - \pi)  \underbrace{\left\{E[\delta \mid D = 1] -E[\delta \mid D = 0]\right\}}_{\text{Differential treatment effect bias}}\end{align*}There are several different representation of the bias when using the naive estimator of mean difference in observed outcomes by treatment status as an estimate for the effect of treatment. We continue with the exposition in Frölich & Sperlich (2019) and Heckman, Urzua, & Vytlacil (2006).\begin{align*}  E[Y\mid D = 1] - E[Y\mid D = 0] & = E[Y^1\mid D = 1] - E[Y^0\mid D = 0] \\                                  & =E[Y^1\mid D = 1]  - E[Y^0\mid D = 1] + E[Y^0\mid D = 1] - E[Y^0\mid D = 0]  \\                                  & = \underbrace{E[Y^1 - Y^0\mid D = 1]}_{TT} + \underbrace{E[Y^0\mid D= 1]- E[Y^0 \mid D = 0]}_{\text{Selection bias}}\end{align*}Now we can simply add and subtract $E[Y_1 - Y_0]$ to get the more economic version.\begin{align*}  E[Y\mid D = 1] - E[Y\mid D = 0] & = \underbrace{E[Y^1 - Y^0]}_{ATE} \\  & + \underbrace{E[Y^1 - Y^0 \mid D = 1] - E[Y^1 - Y^0]}_{\text{Sorting on gains}} \\  & + \underbrace{E[Y^0\mid D = 1] - E[Y^0 \mid D = 0]}_{\text{Sorting on levels}}\end{align*}Sorting on levels is simply a different phrase for selection bias. The exposition in our core textbook is slightly different. Here the term **bias** has two separate components which are **baseline bias** and **differential treatment effect bias**. See the discussion in the book in the subsection on the typical inconsistency and bias of the naive estimator. The term baseline bias refers to the concept of sorting and levels and selection bias.Differential treatment bias is defined as:\begin{align*}E[Y \mid D = 1] - E[Y \mid D = 0] & =  \underbrace{E[\delta]}_{ATE} + \underbrace{\{E[Y^0\mid D= 1] - E[Y^0\mid D= 0]\}}_{\text{Baseline bias}} \\                                      & + \underbrace{(1 - \pi) \{E[\delta \mid D = 1] -E[\delta \mid D = 0]\}}_{\text{Differential treatment effect bias}}\end{align*}The last term is derived derived from the term describing selection on gains by the following decomposition.\begin{align*}E[Y^1 - Y^0] = \pi E [Y^1  - Y^0 \mid D = 1] + (1 - \pi) E [Y^1  - Y^0 \mid D = 0]\end{align*}It is interpreted as the difference in effects between treated and control weighted by  the share of control individuals. It is probably best thought of as an increment to the first term describing the average effect of treatment.  AssumptionsSo, the SUTVA assumption is only necessary but not sufficient to learn about the effect of treatment in light of the biases discussed above. We are still stuck with several unknowns that we need to compute the average effect of treatment.Consider the following two assumptions:\begin{align*}E[Y^1\mid D = 1] = E[Y^1\mid D = 0] \\E[Y^0\mid D = 1] = E[Y^0\mid D = 0] \\\end{align*}and recall our naive estimate \begin{align*}\hat{\delta}_{NAIVE} & = E_N[y_i\mid d_i = 1] - E_N[y_i\mid d_i = 0] \\                     & \xrightarrow{p} E[Y^1\mid D = 1] - E[Y^0\mid D = 0] \end{align*}* What assumptions suffice to estimate the ATE with the naive estimator?    * about potential outcomes for subsets of the population     * about the treatment selection / assignment process  --- Missing data and imputation--- This is an adopted example from Imbens & Rubin (2015).<jupyter_code>df = get_lalonde_data()
df.head()<jupyter_output><empty_output><jupyter_text>We can impute the missing values simply by their average counterpart.<jupyter_code>is_treated = df["D"] == 1

df.loc[~is_treated, "Y_1"] = df["Y_1"].mean()
df.loc[is_treated, "Y_0"] = df["Y_0"].mean()
df.head()
initial_stat = (df["Y_1"] - df["Y_0"]).mean()
print(f"Our estimated treatment effect is {initial_stat:10.2f}")<jupyter_output>Our estimated treatment effect is    1129.49<jupyter_text>However, this does not really account for any uncertainty in our estimate. Can we do better?<jupyter_code>np.random.seed(123)  # set seed to ensure reproducibility
df = get_lalonde_data()  # get the original data

status_counts = df["D"].value_counts().to_dict()

stats = list()
for _ in range(100):
    y_1_sampled = (
        df["Y_1"].dropna().sample(n=status_counts[0], replace=True).values
    )
    y_0_sampled = (
        df["Y_0"].dropna().sample(n=status_counts[1], replace=True).values
    )

    df_boot = df.copy()

    is_treated = df_boot["D"] == 1
    df_boot.loc[is_treated, "Y_0"] = y_0_sampled
    df_boot.loc[~is_treated, "Y_1"] = y_1_sampled

    stat = (df_boot["Y_1"] - df_boot["Y_0"]).mean()
    stats.append(stat)
print(f"Our estimated treatment effect is {np.mean(stats):10.2f}")<jupyter_output>Our estimated treatment effect is    1142.16<jupyter_text>How does the full distribution of estimates look like?<jupyter_code>fig, ax = plt.subplots()
ax.hist(stats)
ax.set_xlabel("Statistic")
ax.set_ylabel("Frequency")
ax.vlines(
    initial_stat, 0, 30, linestyles="--", label="Initial", color="lightgrey"
)
ax.legend()<jupyter_output><empty_output><jupyter_text>Still some limitations remains. For example, we do sample from the empirical distribution of the observed outcomes and not the actual distribution. Phrased differently, we treat the distribution of potential outcomes as known and abstract from any uncertainty in our knowledge about it. --- Extensions of the binary potential outcome model--- * over-time potential outcomes and causal effects    * a single unit over time (time series data)    * many units over time (panel data)    * many-valued treatments  Over-time potential outcomes We explore the case of a single unit over time.**Ingredients** * discrete time periods, $t\in\{1, ..., T\}$* sequence of observed values, $\{y_1, y_2, ..., y_T\}$* treatment initiated in $t^*$ * duration of treatment $k$ Setting up the potential outcome model to explore the basic features of before-and-after designs for a single unit of analsysis.* before the treatment is introduced (for $t < t^*$):\begin{align*}D_t & = 0 \\Y_t & = Y^0_t\\\end{align*}* while the treatment is in place (from $t^*$ through $t^* + k$):\begin{align*}D_t & = 1 \\Y_t & = Y^1_t\\Y^0_t & \text{exists but is counterfactual}\end{align*}* after the treatment ends (for time periods $t > (t^* + k)$):\begin{align*}D_t & = 0 \\Y_t & = Y^0_t\\Y^1_t & \text{exists but is counterfactual}\end{align*} The following is an adapted example from our textbook.  Year of the fire horseWe study the effect of Japanese folk belief that families who give birth to babies will suffer untold miseries. This example does not only illustrative the versatility of the potential outcome framework but also serves as an example that different approaches (informed by domain-expertise) can result in different reasonable imputations for the counterfactual outcome.The treatment indicator is as follows: $D_{1966} = 1$ and $D_{\neq 1966} = 0$ and we are interested in its effect on the birth rate in Japan\begin{align*}\delta_{1966} = y^1_{1966} - y^0_{1966}.\end{align*}A reasonable approach is to estimate it by:\begin{align*}\hat{\delta}_{1966} & = y_{1966} - y^0_{1966}\\\end{align*}<jupyter_code>df = pd.read_csv("material/world_bank.csv", skiprows=4)
df.set_index("Country Code", inplace=True)
df.drop(["Indicator Name", "Indicator Code"], axis=1, inplace=True)

df = df.loc["JPN", "1960":"2017"]
df = df.to_frame()
df.index.name = "Year"
df.columns = ["Birth rate"]

df.sort_index(inplace=True)
df.index = df.index.astype(int)
df.head()<jupyter_output><empty_output><jupyter_text>Let's get to work.<jupyter_code>fig, ax = plt.subplots()
ax.plot(df["Birth rate"].index, df["Birth rate"])
ax.set_ylabel("Birth rate")
ax.set_xlabel("Year")
df.loc[slice(1960, 1970), "Birth rate"]<jupyter_output><empty_output><jupyter_text>We can just take the year before or after treatment?<jupyter_code>estimates = list()
for label, year in [("before", 1965), ("after", 1967)]:
    y_0 = df.loc[year, "Birth rate"]
    y_1 = df.loc[1966, "Birth rate"]
    print(
        f" Using the year {label}, the treatment effect is {y_1 - y_0:10.5f}"
    )
    estimates.append(y_1 - y_0)<jupyter_output>Using the year before, the treatment effect is   -4.90000
 Using the year after, the treatment effect is   -5.60000<jupyter_text>Among demographers, there is the consensus that taking the average of 1963 and 1969 the way to go instead.<jupyter_code>y_0 = df.loc[[1963, 1969], "Birth rate"].mean()
y_1 = df.loc[1966, "Birth rate"]
print(" Another treatment effect is {:10.5f}".format(y_1 - y_0))
estimates.append(y_1 - y_0)<jupyter_output>Another treatment effect is   -4.15000<jupyter_text>Now we have multiple effects of treatment. Which is it?<jupyter_code>labels = ["Before", "After", "Average"]
fig, ax = plt.subplots()
ax.bar(labels, estimates)
ax.set_ylabel("Effect")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12857142857142856, 'fraction_numerical': 0.02506868131868132, 'mean_word_length': 4.281465360899529, 'pattern_counts': {'":': 1, '<': 50, '<?xml version=': 0, '>': 52, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '34626', 'license': 'MIT', 'n_tokens_mistral': 4830, 'n_tokens_neox': 4486, 'n_words': 1676, 'path': 'lectures/02-potential-outcome-model/lecture.ipynb', 'repo_name': 'milakis/microeconometrics', 'size': '14560'}