Chapter 3: Regression Analysis

Simple and multiple linear regression, OLS, diagnostics, and interpretation using R.
Author

Pawan Jangra

Published

February 5, 2025

Learning Objectives: Fit and interpret simple and multiple linear regression models, check OLS assumptions through diagnostic plots, and communicate results clearly.

1 What is Regression?

Regression analysis quantifies the relationship between variables. The goal is to model how changes in explanatory variables (X) are associated with changes in the outcome variable (Y).

The linear regression model: \(Y_i = \beta_0 + \beta_1 X_{1i} + \cdots + \beta_k X_{ki} + \varepsilon_i\)

Code 
library(tidyverse)
library(knitr)
library(broom)

# Wage dataset with labour economics flavor
set.seed(2025)
n <- 800
wage_data <- data.frame(
  edu_years = sample(0:20, n, replace=TRUE, prob=dpois(0:20, 9)+0.01),
  experience = sample(0:40, n, replace=TRUE),
  female = rbinom(n, 1, 0.4),
  urban = rbinom(n, 1, 0.55),
  union = rbinom(n, 1, 0.18)
) |>
  mutate(
    log_wage = 4.5 + 0.08*edu_years + 0.04*experience - 0.003*experience^2
               - 0.25*female + 0.15*urban + 0.12*union + rnorm(n, 0, 0.35),
    wage = round(exp(log_wage))
  )

2 Simple Linear Regression

Let’s first explore the education-wage relationship:

Code 
# Scatter with regression line
ggplot(wage_data, aes(x=edu_years, y=wage)) +
  geom_point(alpha=0.3, color="#2563eb", size=1.5) +
  geom_smooth(method="lm", color="#e11d48", se=TRUE, linewidth=1.2,
              fill="#fca5a5", alpha=0.3) +
  scale_y_continuous(labels=scales::label_number(prefix="₹")) +
  labs(title="Years of Education vs Daily Wage",
       subtitle="OLS regression line with 95% confidence band",
       x="Years of Education", y="Daily Wage (₹)") +
  theme_minimal(base_size=11) +
  theme(panel.grid.minor=element_blank())

Code 
# Fit simple linear regression
model_simple <- lm(wage ~ edu_years, data=wage_data)
tidy(model_simple) |>
  mutate(across(where(is.numeric), round, 3)) |>
  kable(caption="Simple Linear Regression: Wage ~ Education")
Simple Linear Regression: Wage ~ Education
term estimate std.error statistic p.value
(Intercept) 34.985 10.237 3.417 0.001
edu_years 11.321 1.031 10.980 0.000
Code 
glance(model_simple) |>
  select(r.squared, adj.r.squared, sigma, statistic, p.value, nobs) |>
  mutate(across(where(is.numeric), round, 4)) |>
  kable(caption="Model fit statistics")
Model fit statistics
r.squared adj.r.squared sigma statistic p.value nobs
0.1313 0.1302 110.7371 120.5613 0 800

Interpretation: For each additional year of education, daily wages increase by ₹11 on average, holding all else constant.

3 Multiple Linear Regression

Now let’s add more variables:

Code 
model_full <- lm(wage ~ edu_years + experience + I(experience^2) + 
                   female + urban + union, data=wage_data)

tidy(model_full) |>
  mutate(
    term = recode(term,
      `(Intercept)`="Intercept",
      edu_years="Education (years)",
      experience="Experience (years)",
      `I(experience^2)`="Experience² (quadratic)",
      female="Female (=1)",
      urban="Urban (=1)",
      union="Union member (=1)"
    ),
    across(where(is.numeric), round, 4),
    significance = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      TRUE            ~ ""
    )
  ) |>
  kable(caption="Multiple Regression Results (*** p<0.001, ** p<0.01, * p<0.05)")
Multiple Regression Results (*** p<0.001, ** p<0.01, * p<0.05)
term estimate std.error statistic p.value significance
Intercept 136.9419 9.6012 14.2629 0.0000 ***
Education (years) 11.0924 0.6676 16.6161 0.0000 ***
Experience (years) -2.4950 0.8015 -3.1131 0.0019 **
Experience² (quadratic) -0.1063 0.0196 -5.4143 0.0000 ***
Female (=1) -30.0160 5.2114 -5.7597 0.0000 ***
Urban (=1) 15.0396 5.1328 2.9301 0.0035 **
Union member (=1) 38.7287 7.0472 5.4956 0.0000 ***
Code 
# Compare models
bind_rows(
  glance(model_simple) |> mutate(Model="Simple (Education only)"),
  glance(model_full)   |> mutate(Model="Full Model")
) |>
  select(Model, r.squared, adj.r.squared, sigma, AIC) |>
  mutate(across(where(is.numeric), round, 3)) |>
  kable(caption="Model comparison")
Model comparison
Model r.squared adj.r.squared sigma AIC
Simple (Education only) 0.131 0.130 110.737 9805.753
Full Model 0.639 0.637 71.572 9112.400

4 Regression Diagnostics

Always check the four core OLS assumptions:

Code 
par(mfrow=c(2,2), mar=c(4,4,3,1), family="sans")
plot(model_full, which=1:4, 
     col=alpha("#2563eb", 0.4), pch=20, cex=0.7,
     col.smooth="#e11d48", lwd.smooth=2)

Regression diagnostic plots — checking OLS assumptions
Code 
par(mfrow=c(1,1))

What to look for: - Residuals vs Fitted: Should show no pattern (random scatter around 0) - Q-Q Plot: Points should follow the diagonal line (normality of residuals) - Scale-Location: Horizontal band (homoscedasticity — equal variance) - Cook’s Distance: Points above 0.5 are influential observations

5 Coefficient Plot

Code 
tidy(model_full, conf.int=TRUE) |>
  filter(term != "(Intercept)") |>
  mutate(term = recode(term,
    edu_years="Education", experience="Experience",
    `I(experience^2)`="Experience²",
    female="Female", urban="Urban", union="Union")) |>
  ggplot(aes(x=estimate, y=reorder(term, estimate))) +
  geom_vline(xintercept=0, linetype="dashed", color="#94a3b8") +
  geom_errorbarh(aes(xmin=conf.low, xmax=conf.high), height=0.25,
                 color="#7c3aed", linewidth=0.8) +
  geom_point(size=4, color="#2563eb") +
  labs(title="Regression Coefficients with 95% CI",
       x="Coefficient Estimate (₹/unit change)", y=NULL) +
  theme_minimal(base_size=11) +
  theme(panel.grid.minor=element_blank(),
        panel.grid.major.y=element_blank())

6 Log-Linear Wage Equation (Mincer Model)

A standard approach in labour economics — using log wages gives percentage interpretations:

Code 
model_log <- lm(log_wage ~ edu_years + experience + I(experience^2) + 
                  female + urban + union, data=wage_data)

tidy(model_log) |>
  filter(term != "(Intercept)") |>
  mutate(
    term = recode(term,
      edu_years="Education (years)",
      experience="Experience",
      `I(experience^2)`="Experience²",
      female="Female (=1)",
      urban="Urban (=1)",
      union="Union (=1)"),
    pct_effect = round((exp(estimate)-1)*100, 2),
    across(c(estimate, std.error, p.value), round, 4)
  ) |>
  select(term, estimate, std.error, p.value, pct_effect) |>
  rename(`% Effect` = pct_effect) |>
  kable(caption="Mincer Wage Equation (Log-Linear): % change in wage per unit")
Mincer Wage Equation (Log-Linear): % change in wage per unit
term estimate std.error p.value % Effect
Education (years) 0.0769 0.0031 0 7.99
Experience 0.0372 0.0038 0 3.79
Experience² -0.0029 0.0001 0 -0.29
Female (=1) -0.2593 0.0245 0 -22.84
Urban (=1) 0.1452 0.0241 0 15.62
Union (=1) 0.2042 0.0331 0 22.65

Key finding: The gender wage gap is approximately 22.8% — women earn about this much less than men after controlling for education, experience, and other factors.

7 Practice Exercises

TipTry It Yourself
  1. Add state fixed effects to the model using lm(wage ~ ... + factor(state))
  2. Test for heteroscedasticity using lmtest::bptest()
  3. Compute robust standard errors using sandwich::vcovHC()
  4. Interpret the experience² coefficient — at what experience level do wages peak?

Next: Chapter 4 — Time Series Analysis →