Lab: Logit and Probit Models

Overview

In this lab you will model labor force participation as a binary outcome using logit and probit regressions. You will learn why OLS (the linear probability model) has limitations for binary outcomes, how maximum likelihood estimation works conceptually, and — most importantly — how to compute and interpret marginal effects, which are what matter for substantive conclusions.

What you will learn:

• Why the linear probability model can produce impossible predictions
• How to estimate logit and probit models and read the output
• How to compute average marginal effects (AMEs) and marginal effects at the mean (MEMs)
• How to generate predicted probabilities for specific covariate profiles
• How to compare models using pseudo-R-squared, AIC, and classification accuracy

Prerequisites: Familiarity with OLS regression and the concept of probability.

Step 1: Simulate Labor Force Participation Data

We simulate 3,000 individuals with characteristics that determine whether they participate in the labor force.

1library(estimatr)
2library(margins)
3library(modelsummary)
4
5set.seed(42)
6n <- 3000
7
8age <- round(pmin(pmax(rnorm(n, 40, 12), 18), 65))
9educ <- round(pmin(pmax(rnorm(n, 13, 3), 6), 22))
10married <- rbinom(n, 1, 0.55)
11children <- pmin(rpois(n, 1.2), 5)
12spouse_income <- married * rexp(n, rate = 1/30000)
13
14xb <- -3.5 + 0.04 * age - 0.0006 * age^2 + 0.15 * educ -
15    0.25 * children - 0.00001 * spouse_income + 0.3 * married
16
17prob <- 1 / (1 + exp(-xb))
18lfp <- rbinom(n, 1, prob)
19
20df <- data.frame(lfp, age, age_sq = age^2, educ,
21               married, children, spouse_income)
22
23cat("Participation rate:", mean(lfp), "\n")
24summary(df)

Expected output:

Participation rate: 0.612

About 61% of the sample participates in the labor force.

Step 2: The Linear Probability Model (Baseline)

Start with OLS on the binary outcome to see its limitations.

1# Linear Probability Model
2lpm <- lm_robust(lfp ~ age + age_sq + educ + married + children + spouse_income,
3                data = df, se_type = "HC2")
4summary(lpm)
5
6# Check for impossible predictions
7pred_lpm <- predict(lpm, df)
8cat("Predictions < 0:", sum(pred_lpm < 0), "\n")
9cat("Predictions > 1:", sum(pred_lpm > 1), "\n")
10cat("Range:", round(min(pred_lpm), 3), "to", round(max(pred_lpm), 3), "\n")

Expected output:

Predictions < 0: 18
Predictions > 1: 12
Range: [-0.085, 1.062]

The LPM produces 30 predictions outside the [0, 1] range, illustrating its key limitation for binary outcomes.

Step 3: Estimate Logit and Probit Models

1# Logit
2logit <- glm(lfp ~ age + age_sq + educ + married + children + spouse_income,
3           data = df, family = binomial(link = "logit"))
4
5# Probit
6probit <- glm(lfp ~ age + age_sq + educ + married + children + spouse_income,
7            data = df, family = binomial(link = "probit"))
8
9summary(logit)
10summary(probit)
11
12# Coefficient ratio (should be ~1.6)
13cat("\nCoefficient ratio (logit/probit):\n")
14print(round(coef(logit) / coef(probit), 2))

Expected output:

=== Logit Coefficients ===

Pseudo R-squared: 0.0852
Log-likelihood: -1785.2

=== Coefficient Ratio (logit / probit) ===
~1.60 for all coefficients (confirming the logit/probit scaling rule)

The logit/probit coefficient ratio is approximately 1.6, as theory predicts.

Step 4: Compute Marginal Effects

Marginal effects translate logit/probit coefficients into interpretable probability changes.

1# Average Marginal Effects using margins package
2ame_logit <- margins(logit)
3summary(ame_logit)
4
5ame_probit <- margins(probit)
6summary(ame_probit)
7
8# Compare AMEs across models
9cat("\n=== Comparison of AMEs on educ ===\n")
10cat("LPM:", coef(lpm)["educ"], "\n")
11cat("Logit AME:", summary(ame_logit)[summary(ame_logit)$factor == "educ", "AME"], "\n")
12cat("Probit AME:", summary(ame_probit)[summary(ame_probit)$factor == "educ", "AME"], "\n")

Expected output:

=== Average Marginal Effects (Logit) ===

=== Comparison of AMEs ===

The AMEs from logit, probit, and LPM are nearly identical (within 0.005 of each other). An additional year of education raises participation probability by about 3.3 percentage points. An additional child reduces it by about 5.4 percentage points.

Step 5: Predicted Probabilities

Generating predicted probabilities for specific profiles is essential for communicating your results.

1# Predicted probability for specific profiles
2profiles <- data.frame(
3age = c(25, 25, 45, 45),
4age_sq = c(625, 625, 2025, 2025),
5educ = c(12, 16, 12, 16),
6married = c(0, 0, 1, 1),
7children = c(0, 0, 2, 2),
8spouse_income = c(0, 0, 40000, 40000)
9)
10
11profiles$pred_logit <- predict(logit, profiles, type = "response")
12profiles$pred_probit <- predict(probit, profiles, type = "response")
13
14print(profiles[, c("age", "educ", "married", "pred_logit", "pred_probit")])
15
16# Plot predicted probability vs education
17educ_range <- 6:21
18pred_df <- data.frame(age = 40, age_sq = 1600, educ = educ_range,
19                     married = 0, children = 1, spouse_income = 0)
20pred_df$prob <- predict(logit, pred_df, type = "response")
21
22plot(educ_range, pred_df$prob, type = "l", lwd = 2, col = "blue",
23   xlab = "Years of Education", ylab = "Predicted Probability",
24   main = "Predicted LFP Probability by Education")

Expected output:

Predicted Probabilities:

Step 6: Goodness of Fit

1# Pseudo R-squared (McFadden)
2null_loglik <- logLik(glm(lfp ~ 1, data = df, family = binomial))
3cat("McFadden Pseudo R-squared:", 1 - logLik(logit)/null_loglik, "\n")
4
5# AIC / BIC
6cat("Logit AIC:", AIC(logit), " BIC:", BIC(logit), "\n")
7cat("Probit AIC:", AIC(probit), " BIC:", BIC(probit), "\n")
8
9# Classification accuracy
10pred_class <- as.integer(predict(logit, type = "response") > 0.5)
11accuracy <- mean(pred_class == df$lfp)
12cat("\nClassification accuracy:", accuracy, "\n")
13
14# Confusion matrix
15table(Actual = df$lfp, Predicted = pred_class)
16
17# ROC AUC (requires pROC)
18# library(pROC)
19# roc_obj <- roc(df$lfp, predict(logit, type = "response"))
20# auc(roc_obj)

Expected output:

Logit pseudo R-squared:  0.0852
Probit pseudo R-squared: 0.0848

Logit  AIC: 3584.4  BIC: 3627.5
Probit AIC: 3585.1  BIC: 3628.2

Classification accuracy (logit): 0.682

Confusion Matrix:
            Pred 0    Pred 1
Actual 0       652       514
Actual 1       442      1392

ROC AUC: 0.7215

Logit and probit yield virtually identical fit. The model correctly classifies about 68% of observations, with an AUC of ~0.72 indicating good discrimination.

Step 7: Exercises

Try these on your own:

1. Odds ratios. Exponentiate the logit coefficients to obtain odds ratios. Interpret the odds ratio on children in plain language.

2. Interaction effects. Add an interaction between married and children. Compute the marginal effect of an additional child for married vs. unmarried individuals. (Warning: interaction effects in nonlinear models require careful treatment — read Ai and Norton 2003.)

3. Multinomial logit. Extend the model to three outcomes: not in labor force, employed part-time, employed full-time. Use statsmodels.discrete.discrete_model.MNLogit (Python), nnet::multinom (R), or mlogit (Stata).

4. Out-of-sample prediction. Split the data 80/20 and evaluate the model's out-of-sample AUC. Does overfitting appear to be a problem?

Summary

In this lab you learned:

• Logit and probit model binary outcomes using S-shaped link functions that keep predicted probabilities in [0,1]
• The raw coefficients are in log-odds (logit) or z-score (probit) units — in most cases, marginal effects are reported instead
• Average marginal effects (AMEs) are generally preferred over marginal effects at the mean (MEMs)
• Logit and probit typically produce very similar marginal effects in practice
• Pseudo R-squared values are not comparable to OLS R-squared; focus on meaningful marginal effects and adequate discrimination (AUC)
• The linear probability model remains a useful quick check but has known limitations at the tails