Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

7. Multiple Regression Analysis with Qualitative Information

Learning Objectives

Upon completion of this chapter, readers will demonstrate proficiency in categorical variable analysis using Python by:

7.1 Creating dummy variables from categorical data using pandas get_dummies and patsy C() categorical coding functions.

7.2 Estimating regression models with multiple category indicators and interpreting coefficients as group mean differences using statsmodels.

7.3 Implementing interaction terms between dummy and continuous variables using patsy formula syntax with colon operators.

7.4 Conducting Chow tests for structural breaks through F-tests comparing restricted and unrestricted models with statsmodels.

7.5 Estimating linear probability models (LPM) and computing heteroskedasticity-robust standard errors appropriate for binary outcomes.

7.6 Implementing difference-in-differences estimators for program evaluation using interaction terms between time and treatment dummies.

7.7 Visualizing group differences and interactions through categorical plots using seaborn catplot and interaction_plot functions.

7.8 Conducting post-estimation hypothesis tests for coefficient equality across groups using wald_test methods.

Incorporating categorical variables into regression models enables researchers to analyze the impact of non-numeric factors such as gender, marital status, occupation, or geographic region on economic outcomes. This chapter develops the theory and practice of using dummy variables (indicator variables) to represent qualitative factors numerically, allowing their inclusion alongside quantitative regressors in the OLS framework.

The presentation builds systematically from basic dummy variable techniques to advanced applications. We begin with binary dummy variables and their interpretation as intercept shifts (Section 7.1), extend to categorical variables with multiple categories (Section 7.2), examine interaction terms between qualitative and quantitative variables allowing for slope differences (Section 7.3), address the linear probability model for binary dependent variables (Section 7.4-7.6), develop difference-in-differences methods for program evaluation (Section 7.7), and conclude with the Chow test for structural breaks (Section 7.8). Throughout, we implement these methods using Python’s statsmodels library with real datasets from applied econometric research.

import matplotlib.pyplot as plt
import numpy as np  # noqa: F401
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import wooldridge as wool
from scipy import stats

plt.style.use("theme.mplstyle")

7.1 Linear Regression with Dummy Variables as Regressors

Dummy variables (also called indicator variables or binary variables) take on the value 1 or 0 to indicate the presence or absence of a specific attribute or category. In regression analysis, they are used to include qualitative information in a quantitative framework.

Key Concepts:

  1. Reference Category (Base Category): When including a dummy variable, one category is always omitted and serves as the reference category or base category. All coefficient interpretations are relative to this reference group.

  2. Interpretation: The coefficient on a dummy variable represents the difference in the expected value of the dependent variable between the category indicated by the dummy (when dummy=1) and the reference category (when dummy=0), holding all other variables constant.

  3. Dummy Variable Trap: For a categorical variable with gg categories, we include only g−1g-1 dummy variables in the regression. Including all gg dummies would cause perfect multicollinearity (violating MLR.3) because the dummies would sum to the constant term.

Example: If we have a gender variable with categories “male” and “female”, we create one dummy variable:

Example 7.1: Hourly Wage Equation

In this example, we will investigate how gender affects hourly wages, controlling for education, experience, and tenure. We will use the wage1 dataset from the wooldridge package. The dataset includes information on wages, education, experience, tenure, and gender (female=1 if female, 0 if male).

# Load wage dataset for dummy variable analysis
wage1 = wool.data("wage1")

# Examine the gender variable
# GENDER DISTRIBUTION IN DATASET
# Dataset size
pd.DataFrame(
    {
        "Metric": ["Total observations"],
        "Value": [len(wage1)],
    }
)

# Gender distribution details
gender_dist = pd.DataFrame(
    {
        "Gender": ["Males (female=0)", "Females (female=1)"],
        "Count": [(wage1["female"] == 0).sum(), (wage1["female"] == 1).sum()],
        "Percentage": [
            f"{(wage1['female'] == 0).mean():.1%}",
            f"{(wage1['female'] == 1).mean():.1%}",
        ],
    }
)
gender_dist

# Estimate model with dummy variable for gender
dummy_model = smf.ols(
    formula="wage ~ female + educ + exper + tenure",
    data=wage1,
)
dummy_results = dummy_model.fit()

# Create enhanced results table with interpretations
results_table = pd.DataFrame(
    {
        "Coefficient": dummy_results.params.round(4),
        "Std_Error": dummy_results.bse.round(4),
        "t_statistic": dummy_results.tvalues.round(3),
        "p_value": dummy_results.pvalues.round(4),
        "Interpretation": [
            "Baseline wage ($/hr) for males with zero education/experience",
            "Gender wage gap: females earn $1.81/hr less than males",
            "Return to education: $0.57/hr per year of schooling",
            "Return to experience: $0.03/hr per year",
            "Return to tenure: $0.14/hr per year with current employer",
        ],
    },
)

# REGRESSION RESULTS: WAGE EQUATION WITH GENDER DUMMY
# Dependent Variable: wage (hourly wage in dollars)
# Reference Category: male (female=0)
results_table
# Model statistics
pd.DataFrame(
    {
        "Metric": ["R-squared", "Number of observations"],
        "Value": [f"{dummy_results.rsquared:.4f}", int(dummy_results.nobs)],
    }
)
Loading...

Explanation:

Interpretation of Results:

This simple example demonstrates how dummy variables can be used to quantify the effect of qualitative factors like gender on a dependent variable like wages in a multiple regression framework.

Example 7.6: Log Hourly Wage Equation

Using the logarithm of the dependent variable, such as wage, is common in economics. When the dependent variable is in log form, the coefficients on the regressors (after multiplying by 100) can be interpreted as approximate percentage changes in the dependent variable for a one-unit change in the regressor.

In this example, we will use the natural logarithm of hourly wage as the dependent variable and explore the interaction effect between marital status and gender, along with education, experience, and tenure (including quadratic terms for experience and tenure to capture potential non-linear relationships).

wage1 = wool.data("wage1")

reg = smf.ols(
    formula="np.log(wage) ~ married*female + educ + exper +"
    "I(exper**2) + tenure + I(tenure**2)",
    data=wage1,
)
results = reg.fit()

# print regression table:
table = pd.DataFrame(
    {
        "b": round(results.params, 4),
        "se": round(results.bse, 4),
        "t": round(results.tvalues, 4),
        "pval": round(results.pvalues, 4),
    },
)
table  # Display regression results
Loading...

Explanation:

Interpretation of Results (with Interaction Terms):

When interpreting interaction terms involving dummy variables, we must be careful to identify the reference category at each step:

Summary of all four groups relative to baseline (single males):

7.2 Boolean variables

In Python (and programming in general), boolean variables are often used to represent true/false or yes/no conditions. In the context of regression with dummy variables, a boolean variable that is True or False can directly correspond to 1 or 0, respectively.

In the wage1 dataset, the female variable is already coded as 1 for female and 0 for male. We can explicitly create a boolean variable, although it is not strictly necessary in this case as statsmodels and regression formulas in general can interpret 1/0 variables as dummy variables.

wage1 = wool.data("wage1")

# regression with boolean variable:
wage1["isfemale"] = wage1["female"] == 1
reg = smf.ols(formula="wage ~ isfemale + educ + exper + tenure", data=wage1)
results = reg.fit()

# print regression table:
table = pd.DataFrame(
    {
        "b": round(results.params, 4),
        "se": round(results.bse, 4),
        "t": round(results.tvalues, 4),
        "pval": round(results.pvalues, 4),
    },
)
table  # Display regression results
Loading...

Explanation:

Interpretation of Results:

The regression table will be virtually identical to the one in Example 7.1. The coefficient, standard error, t-statistic, and p-value for isfemale will be the same as for female in the previous example.

Key takeaway:

Using a boolean variable (isfemale) that evaluates to True or False is functionally equivalent to using a dummy variable coded as 1 or 0 in regression analysis with statsmodels. statsmodels internally handles boolean True as 1 and False as 0 in regression formulas when they are used as regressors. Therefore, creating an explicit boolean variable like isfemale here doesn’t change the regression outcome compared to directly using the 1/0 coded female variable.

7.3 Categorical Variables

Categorical variables represent groups or categories, and they can have more than two categories. Examples include occupation, industry, region, or education levels (e.g., high school, bachelor’s, master’s, PhD). To include categorical variables in a regression model, we typically use dummy variable encoding. If a categorical variable has k categories, we include k-1 dummy variables in the regression to avoid perfect multicollinearity (the dummy variable trap). One category is chosen as the “reference” or “baseline” category, and the coefficients on the dummy variables represent the difference in the dependent variable between each category and the reference category, holding other regressors constant.

In this section, we will use the cps78_85 dataset from the Wooldridge package, which contains data from the Current Population Survey for 1978 and 1985. We will investigate how gender, region (south/non-south), and union membership affect wages, controlling for education and experience.

cps = wool.data("cps78_85")

# Create categorical gender variable for clearer output
cps["gender"] = cps["female"].map({0: "male", 1: "female"})
cps["region"] = cps["south"].map({0: "non-south", 1: "south"})
cps["union_member"] = cps["union"].map({0: "non-union", 1: "union"})

# table of categories and frequencies for categorical variables:
freq_gender = pd.crosstab(cps["gender"], columns="count")
# Gender distribution
freq_gender

freq_region = pd.crosstab(cps["region"], columns="count")
# Region distribution
freq_region

freq_union = pd.crosstab(cps["union_member"], columns="count")
# Union membership distribution
freq_union
Loading...

Explanation:

Output:

The frequency tables show the counts for each category in the dataset, giving us an overview of the distribution of these categorical variables across the sample.

# directly using categorical variables in regression formula:
reg = smf.ols(
    formula="lwage ~ educ + exper + C(gender) + C(region) + C(union_member)",
    data=cps,
)
results = reg.fit()

# print regression table:
table = pd.DataFrame(
    {
        "b": round(results.params, 4),
        "se": round(results.bse, 4),
        "t": round(results.tvalues, 4),
        "pval": round(results.pvalues, 4),
    },
)
table  # Display regression results
Loading...

Explanation:

Interpretation of Results:

# rerun regression with different reference category:
reg_newref = smf.ols(
    formula="lwage ~ educ + exper + "
    'C(gender, Treatment("male")) + '
    'C(region, Treatment("south"))',
    data=cps,
)
results_newref = reg_newref.fit()

# print results:
table_newref = pd.DataFrame(
    {
        "b": round(results_newref.params, 4),
        "se": round(results_newref.bse, 4),
        "t": round(results_newref.tvalues, 4),
        "pval": round(results_newref.pvalues, 4),
    },
)
# Results with new reference category
table_newref
Loading...

Explanation:

Interpretation of Results:

Key takeaway:

By using C() and Treatment(), we can easily incorporate categorical variables into our regression models and control which category serves as the reference group. The interpretation of the coefficients for categorical variables is always in comparison to the chosen reference category.

7.3.1 ANOVA Tables

ANOVA (Analysis of Variance) tables in regression context help to assess the overall significance of groups of regressors. They decompose the total variance in the dependent variable into components attributable to different sets of regressors. In the context of categorical variables, ANOVA tables can be used to test the joint significance of all dummy variables representing a categorical variable.

cps = wool.data("cps78_85")
# Create categorical variables
cps["gender"] = cps["female"].map({0: "male", 1: "female"})
cps["region"] = cps["south"].map({0: "non-south", 1: "south"})

# run regression:
reg = smf.ols(
    formula="lwage ~ educ + exper + gender + region + union",
    data=cps,
)
results = reg.fit()

# print regression table:
table_reg = pd.DataFrame(
    {
        "b": round(results.params, 4),
        "se": round(results.bse, 4),
        "t": round(results.tvalues, 4),
        "pval": round(results.pvalues, 4),
    },
)
table_reg  # Display regression results
Loading...

Explanation:

# ANOVA table:
table_anova = sm.stats.anova_lm(results, typ=2)
# ANOVA table
table_anova
Loading...

Explanation:

Interpretation of ANOVA Table:

The ANOVA table output will typically include columns like:

For example, in the ANOVA table output:

Key takeaway:

ANOVA tables provide F-tests for the joint significance of regressors, especially useful for categorical variables with multiple dummy variables. They help assess the overall contribution of each set of regressors to explaining the variation in the dependent variable.

7.4 Breaking a Numeric Variable Into Categories

Sometimes, it might be useful to convert a continuous numeric variable into a categorical variable. This can be done for several reasons:

Example 7.8: Effects of Law School Rankings on Starting Salaries

In this example, we will examine the effect of law school rankings on starting salaries of graduates. The lawsch85 dataset from wooldridge contains information on law school rankings (rank), starting salaries (salary), LSAT scores (LSAT), GPA, library volumes (libvol), and cost (cost).

lawsch85 = wool.data("lawsch85")

# define cut points for the rank:
cutpts = [0, 10, 25, 40, 60, 100, 175]

# create categorical variable containing ranges for the rank:
lawsch85["rc"] = pd.cut(
    lawsch85["rank"],
    bins=cutpts,
    labels=["(0,10]", "(10,25]", "(25,40]", "(40,60]", "(60,100]", "(100,175]"],
)

# display frequencies:
freq = pd.crosstab(lawsch85["rc"], columns="count")
# Marital status distribution
freq
Loading...

Explanation:

Output:

The freq output shows the number of law schools falling into each defined rank category (e.g., how many schools are ranked between 0 and 10, between 10 and 25, etc.).

# run regression:
reg = smf.ols(
    formula='np.log(salary) ~ C(rc, Treatment("(100,175]")) +'
    "LSAT + GPA + np.log(libvol) + np.log(cost)",
    data=lawsch85,
)
results = reg.fit()

# print regression table:
table_reg = pd.DataFrame(
    {
        "b": round(results.params, 4),
        "se": round(results.bse, 4),
        "t": round(results.tvalues, 4),
        "pval": round(results.pvalues, 4),
    },
)
table_reg  # Display regression results
Loading...

Explanation:

Interpretation of Results:

# ANOVA table:
table_anova = sm.stats.anova_lm(results, typ=2)
# ANOVA table
table_anova
Loading...

Explanation & Interpretation:

As before, sm.stats.anova_lm(results, typ=2) calculates the ANOVA table for this regression model. The row for C(rc) in the ANOVA table will provide an F-test and p-value for the joint significance of all the dummy variables representing the rank categories. This tests the null hypothesis that law school rank (categorized as rc) has no effect on starting salary, once LSAT, GPA, library volumes, and cost are controlled for. A small p-value would indicate that law school rank category is a statistically significant factor in predicting starting salaries.

7.5 Interactions and Differences in Regression Functions Across Groups

We can allow for regression functions to differ across groups by including interaction terms between group indicators (dummy variables for qualitative variables) and quantitative regressors. This allows the slope coefficients of the quantitative regressors to vary across different groups defined by the qualitative variable.

gpa3 = wool.data("gpa3")

# model with full interactions with female dummy (only for spring data):
reg = smf.ols(
    formula="cumgpa ~ female * (sat + hsperc + tothrs)",
    data=gpa3,
    subset=(gpa3["spring"] == 1),
)
results = reg.fit()

# print regression table:
table = pd.DataFrame(
    {
        "b": round(results.params, 4),
        "se": round(results.bse, 4),
        "t": round(results.tvalues, 4),
        "pval": round(results.pvalues, 4),
    },
)
table  # Display regression results
Loading...

Explanation:

Interpretation of Results:

To get the regression equation for each group:

# F-Test for H0 (the interaction coefficients of 'female' are zero):
hypotheses = ["female:sat = 0", "female:hsperc = 0", "female:tothrs = 0"]
ftest = results.f_test(hypotheses)
fstat = ftest.statistic
fpval = ftest.pvalue

# F-test for interaction terms
pd.DataFrame(
    {
        "Metric": ["F-statistic", "p-value"],
        "Value": [f"{fstat:.4f}", f"{fpval:.4f}"],
    },
)
Loading...

Explanation:

Interpretation:

gpa3 = wool.data("gpa3")

# estimate model for males (& spring data):
reg_m = smf.ols(
    formula="cumgpa ~ sat + hsperc + tothrs",
    data=gpa3,
    subset=(gpa3["spring"] == 1) & (gpa3["female"] == 0),
)
results_m = reg_m.fit()

# print regression table:
table_m = pd.DataFrame(
    {
        "b": round(results_m.params, 4),
        "se": round(results_m.bse, 4),
        "t": round(results_m.tvalues, 4),
        "pval": round(results_m.pvalues, 4),
    },
)
# Regression results for married individuals
table_m
Loading...

Explanation:

Interpretation:

The regression table table_m shows the estimated regression equation specifically for male students. The coefficients are the intercept and slopes of sat, hsperc, and tothrs for males only. These coefficients should be very similar to the coefficients for const, sat, hsperc, and tothrs in the full interaction model (from the first regression in this section), as those were also interpreted as the coefficients for the male group (reference group).

# estimate model for females (& spring data):
reg_f = smf.ols(
    formula="cumgpa ~ sat + hsperc + tothrs",
    data=gpa3,
    subset=(gpa3["spring"] == 1) & (gpa3["female"] == 1),
)
results_f = reg_f.fit()

# print regression table:
table_f = pd.DataFrame(
    {
        "b": round(results_f.params, 4),
        "se": round(results_f.bse, 4),
        "t": round(results_f.tvalues, 4),
        "pval": round(results_f.pvalues, 4),
    },
)
# Regression results for females
table_f
Loading...

Explanation:

Interpretation:

The regression table table_f shows the estimated regression equation specifically for female students. The coefficients are the intercept and slopes of sat, hsperc, and tothrs for females only. These coefficients should correspond to the combined coefficients from the full interaction model. For example:

Key takeaway:

7.6 The Linear Probability Model

When the dependent variable is binary (takes only values 0 or 1), we can still use OLS to estimate the regression model. This is called the linear probability model (LPM) because the fitted values from the regression can be interpreted as estimated probabilities that the dependent variable equals one.

Binary dependent variables in econometrics often represent choices or outcomes:

7.6.1 Model Specification and Interpretation

For a binary outcome yi∈{0,1}y_i \in \{0, 1\}, the linear probability model is:

yi=β0+β1xi1+β2xi2+⋯+βkxik+uiy_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_k x_{ik} + u_i
(1)

Taking expectations conditional on the xx variables:

E(yi∣xi1,…,xik)=β0+β1xi1+⋯+βkxikE(y_i | x_{i1}, \ldots, x_{ik}) = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_k x_{ik}
(2)

Since yiy_i is binary, E(yi∣x)=P(yi=1∣x)E(y_i | \mathbf{x}) = P(y_i = 1 | \mathbf{x}). Therefore, the fitted values represent estimated probabilities:

y^i=β^0+β^1xi1+⋯+β^kxik=P(yi=1∣x)^\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_{i1} + \cdots + \hat{\beta}_k x_{ik} = \widehat{P(y_i = 1 | \mathbf{x})}
(3)

Interpretation of coefficients:

7.6.2 Example: Labor Force Participation

Let’s examine married women’s labor force participation using the mroz dataset:

# Load data
mroz = wool.data("mroz")

# Create binary variable is not already numeric
mroz["inlf_binary"] = (mroz["inlf"] == 1).astype(int)

# Estimate linear probability model
lpm = smf.ols(
    formula="inlf ~ nwifeinc + educ + exper + I(exper**2) + age + kidslt6 + kidsge6",
    data=mroz,
)
results_lpm = lpm.fit()

# Display results
table_lpm = pd.DataFrame(
    {
        "b": round(results_lpm.params, 6),
        "se": round(results_lpm.bse, 6),
        "t": round(results_lpm.tvalues, 4),
        "pval": round(results_lpm.pvalues, 4),
    },
)
print("Linear Probability Model: Labor Force Participation")
table_lpm
Linear Probability Model: Labor Force Participation
Loading...

Interpretation:

Looking at key coefficients:

7.6.3 Predicted Probabilities and Limitations

Source
# Calculate predicted probabilities
mroz["pred_prob"] = results_lpm.fittedvalues

# Examine range of predicted probabilities
print(f"Range of predicted probabilities:")
print(f"Minimum: {mroz['pred_prob'].min():.4f}")
print(f"Maximum: {mroz['pred_prob'].max():.4f}")
print(f"\nProportion of predictions outside [0, 1]:")
outside_range = ((mroz["pred_prob"] < 0) | (mroz["pred_prob"] > 1)).sum()
print(f"{outside_range} out of {len(mroz)} ({100*outside_range/len(mroz):.2f}%)")

# Visualize predicted probabilities
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Histogram of predicted probabilities
ax1.hist(mroz["pred_prob"], bins=30, edgecolor="black", alpha=0.7)
ax1.axvline(0, color="red", linestyle="--", linewidth=2, label="Valid range bounds")
ax1.axvline(1, color="red", linestyle="--", linewidth=2)
ax1.set_xlabel("Predicted Probability")
ax1.set_ylabel("Frequency")
ax1.set_title("Distribution of LPM Predicted Probabilities")
ax1.legend()

# Predicted vs actual
ax2.scatter(
    mroz["educ"],
    mroz["inlf"],
    alpha=0.3,
    label="Actual (0 or 1)",
    s=20,
)
# Sort by educ for smooth line
mroz_sorted = mroz.sort_values("educ")
ax2.plot(
    mroz_sorted["educ"],
    mroz_sorted["pred_prob"],
    color="red",
    linewidth=2,
    label="Predicted probability",
)
ax2.axhline(0, color="gray", linestyle="--", alpha=0.5)
ax2.axhline(1, color="gray", linestyle="--", alpha=0.5)
ax2.set_xlabel("Years of Education")
ax2.set_ylabel("Probability of Labor Force Participation")
ax2.set_title("LPM: Education and Labor Force Participation")
ax2.legend()
ax2.set_ylim(-0.1, 1.1)

plt.tight_layout()
plt.show()
Range of predicted probabilities:
Minimum: -0.3451
Maximum: 1.1272

Proportion of predictions outside [0, 1]:
33 out of 753 (4.38%)
<Figure size 1200x500 with 2 Axes>

Key Limitations of the LPM:

  1. Nonsensical predictions: Fitted values can fall outside the [0,1][0, 1] interval, which is problematic for probabilities

  2. Heteroskedasticity: The error variance is inherently heteroskedastic: Var(ui∣x)=pi(1−pi)\text{Var}(u_i | \mathbf{x}) = p_i(1 - p_i) where pip_i is the true probability

  3. Constant marginal effects: The LPM assumes the effect of xjx_j on the probability is constant across all values of xx, which may be unrealistic

  4. Linearity assumption: True relationship between covariates and probabilities may be nonlinear

When to use LPM:

Addressing heteroskedasticity:

Since LPM has heteroskedastic errors by construction, we should use heteroskedasticity-robust standard errors (covered in Chapter 8):

# Re-estimate with robust standard errors
results_lpm_robust = lpm.fit(cov_type="HC3")

# Compare standard errors
comparison = pd.DataFrame(
    {
        "Coefficient": round(results_lpm.params, 6),
        "SE (usual)": round(results_lpm.bse, 6),
        "SE (robust)": round(results_lpm_robust.bse, 6),
    },
)
print("\nComparison of Standard Errors:")
comparison

Comparison of Standard Errors:
Loading...

Interpretation:

The robust standard errors account for heteroskedasticity and are generally larger than the usual OLS standard errors. When using LPM, always report robust standard errors for valid inference.

Connection to later chapters:

7.7 Policy Analysis and Program Evaluation

Dummy variables are essential tools for evaluating the effects of policy interventions and programs. When units (individuals, firms, regions) receive a treatment or participate in a program, we can use a treatment dummy to estimate the average treatment effect while controlling for observable characteristics.

7.7.1 Regression Adjustment for Treatment Effects

The basic program evaluation model is:

yi=β0+δ0Di+β1xi1+⋯+βkxik+uiy_i = \beta_0 + \delta_0 D_i + \beta_1 x_{i1} + \cdots + \beta_k x_{ik} + u_i
(4)

where:

Interpretation:

7.7.2 Example: Job Training Program Evaluation

Consider evaluating a job training program’s effect on earnings. We use hypothetical data similar to Wooldridge Example 7.9:

# Simulate job training data (similar structure to real evaluations)
np.random.seed(42)
n = 500

# Generate covariates (pre-treatment characteristics)
educ = np.random.randint(10, 17, n)  # Years of education (10-16)
age = np.random.randint(22, 55, n)  # Age (22-54)
prior_earn = np.random.gamma(3, 5, n)  # Prior earnings (in thousands)

# Treatment assignment (not fully random - selection bias)
# People with lower prior earnings more likely to participate
prob_treat = 0.3 + 0.2 * (prior_earn < prior_earn.mean())
treatment = np.random.binomial(1, prob_treat)

# Generate outcome (post-program earnings)
# True treatment effect is $3,000
earnings = (
    2
    + 1.5 * educ
    + 0.3 * age
    + 0.8 * prior_earn
    + 3 * treatment  # True ATE = 3 (thousands)
    + np.random.normal(0, 4, n)
)

# Create DataFrame
train_data = pd.DataFrame(
    {
        "earnings": earnings,
        "treatment": treatment,
        "educ": educ,
        "age": age,
        "prior_earn": prior_earn,
    },
)

print(f"Sample size: {n}")
print(f"Treatment group: {treatment.sum()} ({100*treatment.mean():.1f}%)")
print(f"Control group: {(1-treatment).sum()} ({100*(1-treatment).mean():.1f}%)")
Sample size: 500
Treatment group: 198 (39.6%)
Control group: 302 (60.4%)

Naive comparison (without regression adjustment):

# Simple mean difference (biased if groups differ in observables)
mean_treat = train_data[train_data["treatment"] == 1]["earnings"].mean()
mean_control = train_data[train_data["treatment"] == 0]["earnings"].mean()
naive_effect = mean_treat - mean_control

print(f"\nNaive Comparison (Simple Difference in Means):")
print(f"Average earnings (treatment group): ${mean_treat:.2f}k")
print(f"Average earnings (control group): ${mean_control:.2f}k")
print(f"Naive treatment effect: ${naive_effect:.2f}k")

Naive Comparison (Simple Difference in Means):
Average earnings (treatment group): $45.98k
Average earnings (control group): $45.82k
Naive treatment effect: $0.16k

This naive comparison may be biased if treatment and control groups differ in characteristics like education, age, or prior earnings.

Regression-adjusted estimate:

# Estimate with regression adjustment
eval_model = smf.ols(
    formula="earnings ~ treatment + educ + age + prior_earn",
    data=train_data,
)
results_eval = eval_model.fit()

# Display results
table_eval = pd.DataFrame(
    {
        "b": round(results_eval.params, 4),
        "se": round(results_eval.bse, 4),
        "t": round(results_eval.tvalues, 4),
        "pval": round(results_eval.pvalues, 4),
    },
)
print("\nRegression-Adjusted Treatment Effect:")
table_eval

Regression-Adjusted Treatment Effect:
Loading...

Interpretation:

7.7.3 Difference-in-Differences with Dummies

When we have data before and after program implementation, we can use difference-in-differences (DD) to control for time-invariant unobservables:

yit=β0+δ0(Di×Aftert)+β1Di+β2Aftert+xit′β+uity_{it} = \beta_0 + \delta_0 (D_i \times \text{After}_t) + \beta_1 D_i + \beta_2 \text{After}_t + \mathbf{x}_{it}'\boldsymbol{\beta} + u_{it}
(5)

where:

Intuition:

# Simulate difference-in-differences data
np.random.seed(123)
n_units = 250  # Units per group
n_total = n_units * 2 * 2  # 2 groups * 2 time periods

# Create panel structure
unit_id = np.repeat(np.arange(n_units * 2), 2)
time = np.tile([0, 1], n_units * 2)  # 0 = before, 1 = after
treated = np.repeat([0, 1], n_units * 2)  # First 250 control, next 250 treated

# Interaction term
treat_after = treated * time

# Generate outcome with true DD effect of 5
# treated units have higher baseline (3)
# time trend is 2 for everyone
# treatment effect is 5
y_dd = (
    10
    + 3 * treated  # Pre-existing difference
    + 2 * time  # Time trend
    + 5 * treat_after  # Treatment effect
    + np.random.normal(0, 3, n_total)
)

dd_data = pd.DataFrame(
    {
        "y": y_dd,
        "treated": treated,
        "after": time,
        "treat_after": treat_after,
        "unit_id": unit_id,
    },
)

# Estimate DD model
dd_model = smf.ols(formula="y ~ treated + after + treat_after", data=dd_data)
results_dd = dd_model.fit()

# Display results
table_dd = pd.DataFrame(
    {
        "b": round(results_dd.params, 4),
        "se": round(results_dd.bse, 4),
        "t": round(results_dd.tvalues, 4),
        "pval": round(results_dd.pvalues, 4),
    },
)
print("\nDifference-in-Differences Estimation:")
table_dd

Difference-in-Differences Estimation:
Loading...

Interpretation:

Key assumptions:

7.8 Testing for Structural Breaks: The Chow Test

When we suspect that a regression relationship differs across groups or time periods, we can formally test for structural breaks using the Chow test. This F-test determines whether coefficients differ significantly between subsamples.

7.8.1 The Chow Test for Different Regression Functions

Suppose we have data from two groups (or time periods) and want to test whether the regression function is the same for both.

Pooled model (assumes same coefficients):

yi=β0+β1xi1+⋯+βkxik+uifor all iy_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_k x_{ik} + u_i \quad \text{for all } i
(6)

Separate models (allows different coefficients):

yi=β0,1+β1,1xi1+⋯+βk,1xik+uifor group 1yi=β0,2+β1,2xi1+⋯+βk,2xik+uifor group 2\begin{aligned} y_i &= \beta_{0,1} + \beta_{1,1} x_{i1} + \cdots + \beta_{k,1} x_{ik} + u_i \quad \text{for group 1} \\ y_i &= \beta_{0,2} + \beta_{1,2} x_{i1} + \cdots + \beta_{k,2} x_{ik} + u_i \quad \text{for group 2} \end{aligned}
(7)

Chow test statistic:

F=(SSRpooled−SSR1−SSR2)/(k+1)(SSR1+SSR2)/(n−2(k+1))F = \frac{(\text{SSR}_{\text{pooled}} - \text{SSR}_1 - \text{SSR}_2) / (k+1)}{(\text{SSR}_1 + \text{SSR}_2) / (n - 2(k+1))}
(8)

where:

Under H0:β0,1=β0,2,β1,1=β1,2,…,βk,1=βk,2H_0: \beta_{0,1} = \beta_{0,2}, \beta_{1,1} = \beta_{1,2}, \ldots, \beta_{k,1} = \beta_{k,2}, the statistic follows Fk+1,n−2(k+1)F_{k+1, n-2(k+1)}.

7.8.2 Implementation using Interaction Terms

An equivalent way to conduct the Chow test is using an interaction model:

yi=β0+δ0Di+β1xi1+δ1(Di×xi1)+⋯+βkxik+δk(Di×xik)+uiy_i = \beta_0 + \delta_0 D_i + \beta_1 x_{i1} + \delta_1 (D_i \times x_{i1}) + \cdots + \beta_k x_{ik} + \delta_k (D_i \times x_{ik}) + u_i
(9)

where Di=1D_i = 1 for group 2, 0 for group 1.

Chow test as F-test:

H0:δ0=δ1=⋯=δk=0(all interactions zero)H_0: \delta_0 = \delta_1 = \cdots = \delta_k = 0 \quad \text{(all interactions zero)}
(10)

If we reject H0H_0, the regression functions differ significantly between groups.

7.8.3 Example: Testing for Gender Differences in Returns to Education

Let’s test whether the wage-education relationship differs by gender:

# Use wage1 data
wage1 = wool.data("wage1")

# Pooled model (assumes same coefficients for males and females)
pooled_model = smf.ols(
    formula="np.log(wage) ~ educ + exper + I(exper**2) + tenure + I(tenure**2)",
    data=wage1,
)
results_pooled = pooled_model.fit()

# Separate regressions
# Male model
male_model = smf.ols(
    formula="np.log(wage) ~ educ + exper + I(exper**2) + tenure + I(tenure**2)",
    data=wage1,
    subset=(wage1["female"] == 0),
)
results_male = male_model.fit()

# Female model
female_model = smf.ols(
    formula="np.log(wage) ~ educ + exper + I(exper**2) + tenure + I(tenure**2)",
    data=wage1,
    subset=(wage1["female"] == 1),
)
results_female = female_model.fit()

# Compute Chow test statistic
SSR_pooled = results_pooled.ssr
SSR_male = results_male.ssr
SSR_female = results_female.ssr
SSR_unrestricted = SSR_male + SSR_female

n = len(wage1)
k = 5  # Number of slope coefficients (excluding intercept)
df1 = k + 1  # Restrictions (all coefficients)
df2 = n - 2 * (k + 1)  # Unrestricted model df

chow_stat = ((SSR_pooled - SSR_unrestricted) / df1) / (SSR_unrestricted / df2)
chow_pval = 1 - stats.f.cdf(chow_stat, df1, df2)

print("\nChow Test for Structural Break (Gender Differences):")
print(f"SSR pooled: {SSR_pooled:.4f}")
print(f"SSR unrestricted (male + female): {SSR_unrestricted:.4f}")
print(f"Chow F-statistic: {chow_stat:.4f}")
print(f"Degrees of freedom: ({df1}, {df2})")
print(f"P-value: {chow_pval:.6f}")
print(
    f"Critical value (5% level): {stats.f.ppf(0.95, df1, df2):.4f}",
)

if chow_pval < 0.05:
    print("\nConclusion: REJECT null hypothesis - significant structural break")
    print("The wage-education relationship differs significantly by gender")
else:
    print("\nConclusion: FAIL TO REJECT null hypothesis")
    print("No evidence of structural break by gender")

Chow Test for Structural Break (Gender Differences):
SSR pooled: 93.9113
SSR unrestricted (male + female): 79.9204
Chow F-statistic: 14.9968
Degrees of freedom: (6, 514)
P-value: 0.000000
Critical value (5% level): 2.1162

Conclusion: REJECT null hypothesis - significant structural break
The wage-education relationship differs significantly by gender

Alternative: F-test on interaction model

# Full interaction model
interact_model = smf.ols(
    formula="np.log(wage) ~ female * (educ + exper + I(exper**2) + tenure + I(tenure**2))",
    data=wage1,
)
results_interact = interact_model.fit()

# Test all interaction terms jointly
# This tests H0: all interaction coefficients = 0
# Equivalent to Chow test
interaction_terms = [
    "female:educ",
    "female:exper",
    "female:I(exper ** 2)",
    "female:tenure",
    "female:I(tenure ** 2)",
    "female",  # Different intercept
]

# F-test for joint significance
f_test_result = results_interact.f_test(
    " + ".join([f"{term} = 0" for term in interaction_terms]),
)

print("\nF-test on Interaction Terms (Equivalent to Chow Test):")
print(f"F-statistic: {f_test_result.fvalue:.4f}")
print(f"P-value: {f_test_result.pvalue:.6f}")
print("\nNote: This F-statistic should match the Chow test statistic")

F-test on Interaction Terms (Equivalent to Chow Test):
F-statistic: 14.9968
P-value: 0.000000

Note: This F-statistic should match the Chow test statistic

7.8.4 Chow Test with Different Intercepts Only

Sometimes we want to test only if intercepts differ (not slopes). This is a restricted version:

H0:δ0=0(only intercept difference, slopes same)H_0: \delta_0 = 0 \quad \text{(only intercept difference, slopes same)}
(11)

This is simply a t-test on a single dummy variable:

# Model with only intercept difference (no slope interactions)
intercept_only = smf.ols(
    formula="np.log(wage) ~ female + educ + exper + I(exper**2) + tenure + I(tenure**2)",
    data=wage1,
)
results_intercept = intercept_only.fit()

# Extract female coefficient (intercept shift)
female_coef = results_intercept.params["female"]
female_se = results_intercept.bse["female"]
female_t = results_intercept.tvalues["female"]
female_p = results_intercept.pvalues["female"]

print("\nTest for Different Intercepts Only (Same Slopes):")
print(f"Female coefficient (intercept shift): {female_coef:.6f}")
print(f"Standard error: {female_se:.6f}")
print(f"t-statistic: {female_t:.4f}")
print(f"p-value: {female_p:.6f}")

if female_p < 0.05:
    print("\nConclusion: Intercepts differ significantly by gender")
    print("(but this doesn't test whether slopes also differ)")

Test for Different Intercepts Only (Same Slopes):
Female coefficient (intercept shift): -0.296511
Standard error: 0.035805
t-statistic: -8.2812
p-value: 0.000000

Conclusion: Intercepts differ significantly by gender
(but this doesn't test whether slopes also differ)

Interpretation:

Key insight:

The Chow test is essentially the same as testing the joint significance of all interaction terms in a fully interacted model. The interaction approach is more flexible because you can:

  1. Test for different intercepts only (one dummy)

  2. Test for different slopes only (slope interactions)

  3. Test for any differences (all interactions including intercept)

This provides a comprehensive framework for testing structural stability across groups or time periods.

Chapter Summary

Chapter 7 introduced qualitative (categorical) regressors into the multiple regression framework, expanding our ability to model economic relationships and evaluate policies.

Key Concepts

Dummy Variables:

Multiple Categories:

Interactions:

Linear Probability Model (LPM):

Policy Evaluation:

Chow Test:

Python Implementation Patterns

Creating dummy variables:

Interactions in formulas:

Testing with F-tests:

Robust standard errors:

Common Pitfalls

  1. Dummy variable trap: Including all categories plus intercept causes perfect multicollinearity. Always exclude one reference group.

  2. Misinterpreting log models: In log⁡(y)=β0+δD+⋯\log(y) = \beta_0 + \delta D + \cdots, the coefficient δ\delta is NOT the percentage difference. Use 100×[exp⁡(δ)−1]100 \times [\exp(\delta) - 1] for exact percentage.

  3. Ignoring heteroskedasticity in LPM: LPM has heteroskedastic errors by construction. Always use robust standard errors.

  4. Extrapolating LPM: Predictions outside [0,1][0,1] are nonsensical. LPM works best when most fitted values are in (0.2,0.8)(0.2, 0.8).

  5. Confusing treatment assignment with treatment effect: Selection bias occurs when treatment group differs from control group in unobservables. Regression adjustment only controls for observables.

  6. Forgetting interaction interpretation: When x1x_1 and x2x_2 interact, the marginal effect of x1x_1 depends on the value of x2x_2 (and vice versa).

Connections

To previous chapters:

To later chapters:

Practical Guidance

When to use dummy variables:

Choosing reference group:

Interaction modeling strategy:

  1. Start simple: Test dummy main effect only

  2. If theoretically motivated, add interactions

  3. Test joint significance of interactions (F-test)

  4. If interactions significant, interpret using predicted values or marginal effects at different values

Program evaluation best practices:

Reporting results:

Learning Objectives Covered

  1. 7.1: Define dummy variables and provide examples

  2. 7.2: Interpret dummy variable coefficients in various functional forms

  3. 7.3: Understand interactions among dummy variables

  4. 7.4: Interpret interactions between dummy and continuous variables

  5. 7.5: Incorporate multiple categorical variables using multiple dummies

  6. 7.6: Understand purpose of Chow test for parameter stability

  7. 7.7: Implement Chow test with different intercepts (restricted version)

  8. 7.8: Define and interpret linear probability model for binary outcomes

  9. 7.9: Use regression adjustment for program evaluation

  10. 7.10: Interpret regression results with discrete dependent variables

This chapter provided the tools to handle qualitative information in regression analysis, enabling richer modeling of economic relationships and rigorous evaluation of policies and programs.

merino
6. Multiple Regression Analysis: Further Issues
merino
8. Heteroskedasticity