I am trying to construct the kind of histogram described here (https://blogs.worldbank.org/impactevaluations/judge-leniency-iv-designs-now-not-just-crime-studies). I don't really feel like my code is accomplishing this because my output seems off. Is what I have yielding the correct histogram?

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
from sklearn.utils import resample

# 1. Load the data
df = pd.read_excel("path")

# 2. Calculate the residuals after accounting for Day * Hospital fixed effects
X = pd.get_dummies(df[['weekday', 'hospital_id']], drop_first=True)
X = sm.add_constant(X)
model = sm.OLS(df['removal_stringency_instrument'], X).fit()
df['residuals'] = model.resid

# 3. Plot the histogram of residuals
plt.figure(figsize=(10, 6))
sns.histplot(df['residuals'], kde=False, bins=30, color='lightblue')

# 4. Calculate and plot non-parametric regression using bootstrap for 95% CI
def bootstrap_ci(data, n_iterations=1000):
    """Function to generate bootstrap CI for non-parametric regression."""
    sample_size = len(data)
    bootstrapped_means = np.zeros(n_iterations)

    for i in range(n_iterations):
        sample = resample(data)
        bootstrapped_means[i] = np.mean(sample['TREATMENT'])

    lower = np.percentile(bootstrapped_means, 2.5)
    upper = np.percentile(bootstrapped_means, 97.5)

    return lower, upper

# Grouping by residuals and getting means of treatment
grouped_means = df.groupby('residuals').agg({'TREATMENT': 'mean'}).reset_index()
sns.lineplot(data=grouped_means, x='residuals', y='TREATMENT', color='black')

# Highlighting the 95% CI using bootstrap
lower, upper = bootstrap_ci(df[['residuals', 'TREATMENT']])
plt.axhline(y=lower, color='red', linestyle='--')
plt.axhline(y=upper, color='red', linestyle='--')

plt.title("Histogram with Non-parametric Regression and 95% CI")




Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.