Statistical Programming: Problem Set 3

Exercise 1: Linear Regression: Prediction

We consider a regression problem where we want to predict our dependent variable \(Y\) in terms of the explanatory variable \(X\). In this exercise we try to create an optimal prediction of \(Y\) given \(X\).

Part 1-2: Load data and fit linear regression

Solution:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

# Load data
data = pd.read_csv('https://raw.githubusercontent.com/JanTeichertKluge/data-statprog/refs/heads/main/data/prediction.csv', sep=',', na_values=".")

# Prepare data
x = pd.DataFrame(data, columns=['x'])
y = pd.DataFrame(data, columns=['y'])

# Fit model
linmod = LinearRegression()
linmod.fit(x, y)

# Predict
ypred_linmod = linmod.predict(x)

# Calculate MSE
MSE_ins = np.mean((ypred_linmod - y)**2)
print(f"In-sample MSE: {MSE_ins:.4f}")

Part 3: Out-of-sample prediction

Solution:

# Load test data
datatest = pd.read_csv('https://raw.githubusercontent.com/JanTeichertKluge/data-statprog/refs/heads/main/data/predictiontest.csv', sep=',', na_values=".")
datatest.rename(columns={'xtest': 'x', 'ytest': 'y'}, inplace=True)

xtest = pd.DataFrame(datatest, columns=['x'])
ytest = pd.DataFrame(datatest, columns=['y'])

# Predict on test data
y_linmod_oos = linmod.predict(xtest)

# Calculate out-of-sample MSE
MSE_oos = np.mean((y_linmod_oos - ytest)**2)
print(f"Out-of-sample MSE: {MSE_oos:.4f}")

Part 4: Polynomial regression

Solution:

from sklearn.preprocessing import PolynomialFeatures

# Create polynomial features
poly_q = 5
poly_transformer = PolynomialFeatures(poly_q, include_bias=False)
x_poly = poly_transformer.fit_transform(x)

# Fit polynomial model
polymod = LinearRegression()
polymod.fit(x_poly, y)

# Predictions
y_poly_ins = polymod.predict(x_poly)
x_poly_test = poly_transformer.transform(xtest)
y_poly_oos = polymod.predict(x_poly_test)

# Calculate MSEs
MSE_poly_ins = np.mean((y_poly_ins - y)**2)
MSE_poly_oos = np.mean((y_poly_oos - ytest)**2)

print(f"In-sample MSE: {MSE_poly_ins:.4f}")
print(f"Out-of-sample MSE: {MSE_poly_oos:.4f}")

Solution to Exercise 1.1.

Solution to Exercise 1.2.

Solution to Exercise 1.3.

Generate out-of-sample predictions

Compute the out-of-sample MSE in prediction.

Solution to Exercise 1.4.

Fit the polynomial model.

Generate in-sample and out-of-sample predictions and compute the in-sample and out-of-sample MSE.

Solution to Exercise 1.5.

Solution to Exercise 1.6.

Exercise 2: Linear Regression: Inference I - Case Study Oregon Health Experiment

In this problem, we analyze a random subset from the Oregon Health Experiment Data. In this large-scale experiment, access to Medicaid (U.S. public health insurance) was provided randomly in a lottery. More information can be found here. The problem set is based on the published study by Finkelstein et al. (2012, Quarterly Journal of Economics).

Load the data and create a pandas data frame.
Our dependent variable is the number of doctor visits in the last 6 months (“docvis”). Inspect the data set and generate a frequency table for the outcome variable”docvis” as well as a barplot.

Hint: For the frequency table use pd.crosstab .

We are interested in estimating an intent-to-treatment effect, i.e. how does access to Medicaid increase the number of doctor visits. Thus, our regression model is \[ DOCVIS_i = \beta_0 + \beta_1 * TREATMENT_i + X_i'\beta + \varepsilon_i,\] where \(DOCVIS_i\) is the number of doctor visits, \(TREATMENT_i\) is access to health insurance. The controls \(X_i\) contain the variables:

dddraw_sur_2 ddddraw_sur_3 ddddraw_sur_4 ddddraw_sur_5 ddddraw_sur_6 ddddraw_sur_7 dddnumhh_li_2 dddnumhh_li_3 ddddraXnum_2_2 ddddraXnum_2_3 ddddraXnum_3_2 ddddraXnum_3_3 ddddraXnum_4_2 ddddraXnum_5_2 ddddraXnum_6_2 ddddraXnum_7_2

Hint: At the very end of this exercise, you’ll find the list with column names that you need to build your data frame.
Is there a significant effect of \(TREATMENT\), i.e., access to Medicaid?
Now, we would like to extend the model and include variables on gender, education and race/ethnicity in addition to the technical regressors

female_list, hhinc_pctfpl_12m, age2008, edu_12m_2, edu_12m_3, edu_12m_4, english_list, zip_msa, race_white_12m, race_black_12m, race_hisp_12m

How do you interpret the coefficients on female, age and income? Does your conclusion from part 2. change?
(optional) Does the treatment dummy vary with the observable characteristics? Run a regression that includes the interaction of the treatment variable with the non-technical regressors. Interpret your findings.

Hint: Use the formula interface of the statsmodel.formula.api to generate interactions.
How would you summarize your results? Are the OLS conditions met in the Oregon Health Example?

Parts 1-2: Data Loading and Exploration

Solution:

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

# Load data
Oregdata = pd.read_csv('https://raw.githubusercontent.com/JanTeichertKluge/data-statprog/refs/heads/main/data/Oregon.csv', sep=',', na_values=".")

# Frequency table
freq_table = pd.crosstab(columns=Oregdata["docvis"], index="count")
print(freq_table)

# Histogram
docvis = Oregdata['docvis'].values
plt.hist(docvis, bins=30, density=True)
plt.xlabel('Doctor Visits')
plt.ylabel('Density')
plt.title('Distribution of Doctor Visits')
plt.show()

Part 3: Basic Regression with Technical Controls

Solution:

import pandas as pd
import statsmodels.api as sm

columnnames = ['treatment', 'ddddraw_sur_2', 'ddddraw_sur_3', 'ddddraw_sur_4',
               'ddddraw_sur_5', 'ddddraw_sur_6', 'ddddraw_sur_7',
               'dddnumhh_li_2', 'dddnumhh_li_3', 'ddddraXnum_2_2', 'ddddraXnum_2_3',
               'ddddraXnum_3_2', 'ddddraXnum_3_3', 'ddddraXnum_4_2',
               'ddddraXnum_5_2', 'ddddraXnum_6_2', 'ddddraXnum_7_2']

X1c = pd.DataFrame(Oregdata, columns=columnnames)
X1c = sm.add_constant(X1c)

ols1 = sm.OLS(docvis, X1c)
ols1 = ols1.fit(cov_type="HC3")
print(ols1.summary())

Part 4: Extended Model with Demographic Controls

Solution:

import pandas as pd
import statsmodels.api as sm

# All columns
X2c = pd.DataFrame(Oregdata, columns=[
    'treatment', 'female_list', 'hhinc_pctfpl_12m', 'age2008',
    'edu_12m_2', 'edu_12m_3', 'edu_12m_4', 'english_list', 'zip_msa',
    'race_white_12m', 'race_black_12m', 'race_hisp_12m',
    'ddddraw_sur_2', 'ddddraw_sur_3', 'ddddraw_sur_4',
    'ddddraw_sur_5', 'ddddraw_sur_6', 'ddddraw_sur_7',
    'dddnumhh_li_2', 'dddnumhh_li_3', 'ddddraXnum_2_2', 'ddddraXnum_2_3',
    'ddddraXnum_3_2', 'ddddraXnum_3_3', 'ddddraXnum_4_2',
    'ddddraXnum_5_2', 'ddddraXnum_6_2', 'ddddraXnum_7_2'
])

X2c = sm.add_constant(X2c)
ols2 = sm.OLS(docvis, X2c).fit(cov_type="HC3")
print(ols2.summary())

# Interpretation:
# - Female: women have more doctor visits than men
# - Age: older people have more doctor visits
# - Income: higher income associated with more visits
# - Treatment: still significant, similar magnitude to basic model

Alternative Solutions to Exercise 2.1.

Alternative Solution to Exercise 2.2.

(+) Boxplot.

We are interested in estimating an intent-to-treatment effect, i.e. how does access to Medicaid increase the number of doctor visits. Thus, our regression model is

\[ DOCVIS_i = \beta_0 + \beta_1 * TREATMENT_i + X_i'\beta + \varepsilon_i,\]

where \(DOCVIS_i\) is the number of doctor visits, \(TREATMENT_i\) is access to health insurance. The controls \(X_i\) contain the variables:

dddraw_sur_2 ddddraw_sur_3 ddddraw_sur_4 ddddraw_sur_5 ddddraw_sur_6 ddddraw_sur_7 dddnumhh_li_2 dddnumhh_li_3 ddddraXnum_2_2 ddddraXnum_2_3 ddddraXnum_3_2 ddddraXnum_3_3 ddddraXnum_4_2 ddddraXnum_5_2 ddddraXnum_6_2 ddddraXnum_7_2

Solution to Exercise 2.3.

Solution to Exercise 2.4.

Solution to Exercise 2.5. (optional)

Solution to Exercise 2.6.

DOCVIS is a count variable \(\rightarrow\) normality approximation might be questionnable. Asymptotically, confidence bands (etc.) are valid, but not exactly.

Exercise 3: Linear Regression: Inference II - Simulation for the Linear Regression Model (OLS)

For the beginning, let us consider a univariate regression model

\[Y = X\beta + \varepsilon,\] where \(Y\) is the outcome variable and \(X\) is a regressor. \(\varepsilon\) is drawn from a Normal distribution with variance \(\mu = 0\) and \(\sigma = 1\) and independent from \(X\). We are interested in inference on \(\beta\).

In this exercise, we simulate data according to the regression model above in order to provide evidence in favor of the theoretical results we learned in the lecture.

Set up a data generating process (DGP) according to the regression model above. Write a function with inputs \(\beta\) and \(n\) and outputs \(Y\) and \(X\). Generate \(X\) as \(X \sim_{i.i.d} N(0,1)\). Start with \(n=20\) observations first and use \(\beta = 1\) in your example.
Set up a simulation study to estimate the Bias and Standard Error of \(\hat{\beta}_1\). How do the results change if \(n\) increases? Illustrate your results with an appropriate graphic. Do your result support the claim that the OLS estimator is an unbiased estimator? What can you say about estimation uncertainty?
Repeat the simulation from before \(R=100\) times and report the average results (i.e. the average Bias and Standard Error over the \(R\) repetitions).

Hint: Try to write a function that executes the estimation from before automatically and that repeats the calculations \(R\) times.
In the lecture, we have learned that the t-statistic \(t\) for the regression coefficient \(\beta\) is asymptotically normal if \(n\) is large.

Use the simulation study to provide evidence that in the setting above it holds that: \[\sqrt{n}(\hat{\beta} - \beta) \xrightarrow[]{d} N(0,1)\] Hint: Repeat the simulation from above. Now, we need to save the \(\hat{\beta}\) so that we can compute \(\sqrt{n}(\hat{\beta}-\beta)\). We can illustrate the results by generating a bar plot similar to that in Lecture 1.
(optional) Illustrate that using a t-test allows to control the type-1-error at level of α or below.

Part 1: Data Generating Process (DGP)

Solution:

import numpy as np
import statsmodels.api as sm

def DGP(beta, n):
    X = np.random.normal(0, 1, n)
    epsilon = np.random.normal(0, 1, n)
    Y = np.dot(X, beta) + epsilon
    return Y, X

# Test
n = 20
beta = 1
outcome, regressor = DGP(beta, n)

# Run regression
ols_test = sm.OLS(outcome, regressor).fit()
print(ols_test.summary())

Parts 2-3: Bias and Standard Error Simulation

Solution:

import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

beta = 1
nobs = [10, 20, 30, 40, 50, 100, 200, 400]

Bias = np.zeros(len(nobs))
SE = np.zeros(len(nobs))
Estim = np.zeros(len(nobs))

np.random.seed(1234)

for i in range(len(nobs)):
    n = nobs[i]
    outcome, regressor = DGP(beta, n)
    ols = sm.OLS(outcome, regressor).fit()

    Estim[i] = ols.params[0]
    Bias[i] = beta - ols.params[0]
    SE[i] = ols.bse[0]

print("Results:")
for i in range(len(nobs)):
    print(f"n={nobs[i]:3d}: β̂={Estim[i]:.4f}, Bias={Bias[i]:.4f}, SE={SE[i]:.4f}")

# Plot
plt.figure(figsize=(10, 6))
plt.plot(nobs, Estim, 'o-', label='Estimate')
plt.axhline(y=beta, color='r', linestyle='--', label='True β')
plt.plot(nobs, SE, 's-', label='SE')
plt.legend()
plt.show()

Alternative Solutions to Exercise 3.1.

Solution to Exercise 3.2.

Solution to Exercise 3.3.

Solution to Exercise 3.4.

We can already show that \(\hat{\beta}\) is \(\hat{\beta} \sim N(\beta,\sigma^2_{\hat{\beta}})\)

We want to illustrate that \(\sqrt{n}(\hat{\beta} - \beta) \xrightarrow[]{d} N(0,1)\).

Solution to Exercise 3.5. (optional)

The type-1-error control says that, if the \(H_0\) is true, we reject the \(H_0\) with less than probability \(\alpha\).

Hence, we need to simulate a situation where the null hypothesis is true and count how many times we reject the \(H_0\). We should not reject more often than in \(\alpha\) percent of the cases.

In our example, we test \(H_0: \beta = 0\) vs. \(H_A: \beta \neq 0\).

The test is rejected if the p-value is smaller than \(\alpha\).