# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Potential outcome model

# +
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

pd.options.display.float_format = "{:,.2f}".format
# -

# The National Health Interview Survey (NHIS) data is collected on U.S. households since 1957. It covers a broad range of health-related topics from medical conditions, health insurance, and the number of doctor visits to measures of physical activity. Here we focus on indicators relevant for the Potential outcome model (POM) framework. In particular, we will compare the health status of hospitalized and non-hospitalized individuals in 2018. For this purpose, we use answers to the survey question During the past 12 months, has the respondent been hospitalized overnight? with potential answers Yes and No which we code as one and zero. Further, we consider answers to the questions Would you say your health, in general, is excellent, very good, good, fair, poor? where responses are coded as one for poor health up to five for excellent health. The survey also collects data on relevant characteristics as sex, age, level of education, hours worked last week, and total earnings.

# ## Task A.1
#
# Open a Jupyter Notebook and import the data set nhis-initial.xslx (available at https://bit.ly/nhis-initial). Try to think of ways to answer the following questions: Are there more females or males? Are there more individuals who hold a degree or not?. Now try to relate individual characteristics to the hospitalization status. Are high or low earners/old or young people more often hospitalized?

df = pd.read_excel("data/nhis-initial.xls")
df.index.set_names("Individual", inplace=True)
df.head()

# We want to study average age and working hours.

stat = df["age"].mean()
print(f"Average age in the sample is {stat:.2f}")

stat = df["hours"].mean()
print(f"Average of working hours per week in the sample is {stat:.0f}")

for column in ["sex", "education", "earnings", "health"]:

    fig, ax = plt.subplots()

    info = df[column].value_counts(normalize=True).to_frame()
    x, y = info.index, info.to_numpy().flatten()

    ax.bar(x, y)

    ax.set_xlabel(column.capitalize())

    ax.set_ylim(None, 1)
    ax.set_ylabel("Share")

    fig.show()

# Now try to relate individual characteristics to the hospitalization status. 

df.groupby("hospitalized")[["age", "hours"]].mean()

for column in ["sex", "education", "earnings"]:

    fig, ax = plt.subplots()

    rslt = df.groupby(column)["hospitalized"].mean()
    x, y = rslt.index, rslt.to_numpy()

    ax.bar(x, y)

    ax.set_ylabel("Share")

    fig.show()

# ## Task A.2
#
# Compute the average health status of hospitalized and non-hospitalized individuals. Who is healthier on average? What could be a reason for this difference?

df.groupby("hospitalized")["health"].mean().to_frame()

# ## Task A.3
#
# Adjust the data set for the POM framework, with health status as the outcome and hospitalization as the treatment status.

# +
df = pd.read_excel("data/nhis-initial.xls")

df.rename(columns={"health": "Y", "hospitalized": "D"}, inplace=True)

df["Y_1"] = np.where(df["D"] == 1, df["Y"], np.nan)
df["Y_0"] = np.where(df["D"] == 0, df["Y"], np.nan)

df.head()
# -

# ## Task A.4
#
# Compute the naive estimate for the average treatment effect (ATE)

stat = df["Y_1"].mean() - df["Y_0"].mean()
print(f"Our naive estimate is {stat:.1f}")

# ## Task B.1
#
# As we’ve seen in the lecture, in reality, we can only ever observe one counterfactual; however, when simulating data, we can bypass this problem. The (simulated) data set nhis-simulated.xslx (available at https://bit.ly/nhis-simulated)
# contains counterfactual outcomes, i.e., outcomes under control for individuals assigned to the treatment group and vice versa. Derive and compute the average outcomes in the two observable and two unobservables states. Design them similar
# to Table 2.3 in Morgan & Winship (2014).

df = pd.read_excel("data/nhis-simulated.xls")
df.index.set_names("Individual", inplace=True)
df.head()

# +
rslt = df.groupby("D")[["Y_1", "Y_0"]].mean()

rslt.columns = ["E[Y_1|D]", "E[Y_0|D]"]
rslt.index = ["Untreated", "Treated"]

rslt
# -

# ## Task B.2 
#
# From here on we assume that 5% of the population take the treatment. Derive and explain Equation 2.12 from Morgan & Winship (2014) for the naive estimator as a decomposition of true ATE, baseline bias, and differential treatment effect bias.

# This derivation is straightforward.

# ## Task B.3
#
# Compute the naive estimate and true value of the ATE for the simulated data. Is the naive estimator upwardly or downwardly biased? Calculate the baseline bias and differential treatment effect bias. How could we interpret these biases in our framework of health status of hospitalized and non-hospitalized respondents?

# +
pi = 0.05

# naive estimate
naive = rslt.loc["Treated", "E[Y_1|D]"] - rslt.loc["Untreated", "E[Y_0|D]"]

# baseline bias
base = rslt.loc["Treated", "E[Y_0|D]"] - rslt.loc["Untreated", "E[Y_0|D]"]

# differential effect
diff = 0
diff += rslt.loc["Treated", "E[Y_1|D]"] - rslt.loc["Treated", "E[Y_0|D]"]
diff -= rslt.loc["Untreated", "E[Y_1|D]"] - rslt.loc["Untreated", "E[Y_0|D]"]
diff *= 1 - pi

# true average treatment effect
true = 0
true += pi * (rslt.loc["Treated", "E[Y_1|D]"] - rslt.loc["Treated", "E[Y_0|D]"])
true += (1 - pi) * (rslt.loc["Untreated", "E[Y_1|D]"] - rslt.loc["Untreated", "E[Y_0|D]"])
print(f"naive: {naive:.2f}, base: {base:.2f}, diff: {diff:.2f}, true: {true:.2f}")

# We can also test the relationships just to be sure.
np.testing.assert_almost_equal(true, naive - (base + diff), decimal=10)
# -

# ## Task B.4 
#
# Under which assumptions does the naive estimator provide the ATE?

# We need the *stable unit treatment value assumption* and independence between potential outcomes and the treatment.

# ## References
#
# * <NAME>., and <NAME>. (2014). [*Counterfactuals and causal inference: Methods and principles for social research*](https://www.cambridge.org/de/academic/subjects/sociology/sociology-general-interest/counterfactuals-and-causal-inference-methods-and-principles-social-research-2nd-edition?format=PB). Cambridge, England: Cambridge University Press.
#
# * <NAME>., and <NAME>. (2009). [*Mostly harmless econometrics: An empiricists companion*](https://press.princeton.edu/titles/8769.html). Princeton, NJ: Princeton University Press.
#
# * [National Health Interview Survey.](https://www.cdc.gov/nchs/nhis/index.htm) (2018). National Center for Health Statistics.

StarcoderdataJupyterScriptsDedupFiltered