library(tidyverse)
library(modelsummary)
library(marginaleffects)
library(broom)Session 3: Modelling Binary and Categorical Outcomes
Session 3 · Quantitative Data Analysis
1 Introduction
1.1 Business question
Can we predict which firms are likely to exit the market in the next year?
This question is highly relevant in business practice: investors, banks, and regulators all have an interest in understanding which firms are at risk of failure. But answering it with standard regression runs into an immediate problem — the outcome we care about, firm exit, is a binary variable: a firm either exits (1) or it doesn’t (0).
In this session, we work through why ordinary least squares (OLS) breaks down for such outcomes, develop the intuition for logistic regression, and learn how to estimate, interpret, and communicate logistic regression results in R.
1.2 Packages and data
We use the bisnode-firms dataset from the open-access textbook by Békés & Kézdi (2021). It contains 19,036 Hungarian manufacturing and services firms, observed over time. Our outcome variable is exit, a binary indicator equal to 1 if the firm closed in the following year.
firms <- read_csv("data/bisnode-firms.csv")
glimpse(firms)Rows: 17,506
Columns: 11
$ comp_id <dbl> 1002029, 1011889, 1014183, 1018301, 1022796, 1035705, 10…
$ exit <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0,…
$ sales_log <dbl> 13.943477, 12.980031, 11.773208, 8.813174, 10.169549, 10…
$ employees_log <dbl> -0.78015860, 0.48317409, -0.33508429, -1.87877086, -2.48…
$ profit_margin <dbl> 0.008554417, 0.223101892, -0.018128871, -0.164738285, -0…
$ age <dbl> 6, 20, 11, 8, 11, 1, 14, 5, 2, 4, 2, 1, 2, 4, 16, 1, 23,…
$ liquidity <dbl> 1.69277357, 8.46944789, 24.16025018, 3.05624997, 0.90785…
$ foreign <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,…
$ female_ceo <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1,…
$ industry <dbl> 27, 55, 55, 56, 56, 56, 29, 56, 56, 27, 41, 55, 56, 28, …
$ urban <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0,…Let us first look at the outcome variable:
firms |>
count(exit) |>
mutate(share = n / sum(n))# A tibble: 2 × 3
exit n share
<dbl> <int> <dbl>
1 0 16294 0.931
2 1 1212 0.0692
Note
Notice that exit events are relatively rare — this is typical of survival data and has implications for model evaluation that we will revisit in Session 9 (prediction vs. inference).
A quick look at the main predictors:
firms |>
select(exit, sales_log, employees_log, profit_margin, age) |>
summary() exit sales_log employees_log profit_margin
Min. :0.00000 Min. : 1.309 Min. :-2.485 Min. :-1.0000000
1st Qu.:0.00000 1st Qu.: 9.828 1st Qu.:-2.331 1st Qu.:-0.2542683
Median :0.00000 Median :10.895 Median :-1.473 Median : 0.0001222
Mean :0.06923 Mean :10.967 Mean :-1.256 Mean :-0.1428446
3rd Qu.:0.00000 3rd Qu.:12.021 3rd Qu.:-0.626 3rd Qu.: 0.0322967
Max. :1.00000 Max. :18.478 Max. : 5.019 Max. : 1.0000000
age
Min. : 0.000
1st Qu.: 3.000
Median : 8.000
Mean : 8.991
3rd Qu.:15.000
Max. :34.000 2 Why OLS Fails for Binary Outcomes
2.1 The Linear Probability Model
The simplest approach is to regress the binary outcome directly on predictors using OLS. This is called the Linear Probability Model (LPM):
\[\text{exit}_i = \beta_0 + \beta_1 \cdot \text{sales\_log}_i + \beta_2 \cdot \text{profit\_margin}_i + \varepsilon_i\]
lpm <- lm(exit ~ sales_log + profit_margin + employees_log,
data = firms)
summary(lpm)
Call:
lm(formula = exit ~ sales_log + profit_margin + employees_log,
data = firms)
Residuals:
Min 1Q Median 3Q Max
-0.27563 -0.07935 -0.05414 -0.03527 1.04989
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.216657 0.022001 9.848 < 2e-16 ***
sales_log -0.013902 0.001729 -8.041 9.45e-16 ***
profit_margin -0.087885 0.005908 -14.875 < 2e-16 ***
employees_log 0.005976 0.002556 2.338 0.0194 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2498 on 17502 degrees of freedom
Multiple R-squared: 0.03178, Adjusted R-squared: 0.03162
F-statistic: 191.5 on 3 and 17502 DF, p-value: < 2.2e-16The coefficient interpretation is straightforward: a one-unit increase in profit_margin is associated with a -0.088 change in the probability of exit.
2.1.1 The problem: out-of-bounds predictions
lpm_fitted <- augment(lpm, data = firms)
lpm_fitted |>
summarise(
min_pred = min(.fitted),
max_pred = max(.fitted),
n_below_zero = sum(.fitted < 0),
n_above_one = sum(.fitted > 1)
)# A tibble: 1 × 4
min_pred max_pred n_below_zero n_above_one
<dbl> <dbl> <int> <int>
1 -0.0651 0.277 338 0We have 338 predicted values below zero and 0 above one. These are mathematically impossible as probabilities.
ggplot(lpm_fitted, aes(x = sales_log, y = .fitted)) +
geom_point(alpha = 0.2, size = 0.8, color = "#2c7bb6") +
geom_hline(yintercept = c(0, 1), linetype = "dashed", color = "firebrick") +
annotate("label", x = max(firms$sales_log, na.rm = TRUE) - 1,
y = -0.05, label = "Impossible (P < 0)",
color = "firebrick", size = 3) +
labs(
title = "LPM: Fitted probabilities of firm exit",
x = "Log sales", y = "Predicted P(exit = 1)"
) +
theme_minimal()
2.1.2 The problem: heteroskedastic errors
When Y is binary, the variance of the error term depends on X by construction:
\[\text{Var}(\varepsilon_i | x_i) = P_i(1 - P_i)\]
This means OLS standard errors are not valid, so hypothesis tests and confidence intervals are unreliable. We can verify this visually — the residuals should show a clear fan pattern.
ggplot(lpm_fitted, aes(x = .fitted, y = .resid)) +
geom_point(alpha = 0.2, size = 0.8) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "LPM residuals vs. fitted values",
x = "Fitted value", y = "Residual") +
theme_minimal()
Warning
Despite these problems, the LPM is widely used in economics, especially with panel data. It can serve as a useful benchmark. But you should always be aware of its limitations.
3 Logistic Regression: Intuition
3.1 The logistic function
We need a function that maps any real number to a probability in [0, 1]. The logistic (sigmoid) function does exactly this:
\[\sigma(z) = \frac{e^z}{1 + e^z}\]
sigma <- function(z) exp(z) / (1 + exp(z))
tibble(z = seq(-6, 6, by = 0.05)) |>
mutate(prob = sigma(z)) |>
ggplot(aes(z, prob)) +
geom_line(color = "#1a9641", linewidth = 1.3) +
geom_hline(yintercept = c(0, 0.5, 1), linetype = "dashed", color = "grey60") +
labs(title = "The logistic function",
x = "z (linear combination of predictors)",
y = "Predicted probability") +
theme_minimal()
The idea: let \(z = \beta_0 + \beta_1 x_1 + \ldots\) be the linear combination, and then:
\[\hat{P}(\text{exit} = 1 \mid \mathbf{x}) = \sigma(z) = \frac{e^{\mathbf{X}\boldsymbol{\beta}}}{1 + e^{\mathbf{X}\boldsymbol{\beta}}}\]
3.2 From the sigmoid to the logit
Rearranging the expression above gives the logit link:
\[\ln\left(\frac{P}{1-P}\right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots\]
The left-hand side is the log-odds. The relationship between probability, odds, and log-odds:
| P(exit) | Odds | Log-odds |
|---|---|---|
| 0.10 | 0.11 | −2.20 |
| 0.25 | 0.33 | −1.10 |
| 0.50 | 1.00 | 0.00 |
| 0.75 | 3.00 | 1.10 |
| 0.90 | 9.00 | 2.20 |
3.3 Maximum likelihood estimation
Logistic regression is not estimated by OLS but by maximum likelihood estimation (MLE). The idea: find the coefficient values β that make the observed data most probable. In practice, R’s glm() handles this for you.
4 Logistic Regression in R
4.1 Fitting the model
The syntax is almost identical to lm(). The key difference is family = binomial:
logit_model <- glm(exit ~ sales_log + profit_margin + employees_log,
family = binomial,
data = firms)
summary(logit_model)
Call:
glm(formula = exit ~ sales_log + profit_margin + employees_log,
family = binomial, data = firms)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.048479 0.331230 -3.165 0.00155 **
sales_log -0.168696 0.026067 -6.472 9.7e-11 ***
profit_margin -1.047630 0.085995 -12.183 < 2e-16 ***
employees_log 0.006127 0.043090 0.142 0.88694
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 8810.8 on 17505 degrees of freedom
Residual deviance: 8333.5 on 17502 degrees of freedom
AIC: 8341.5
Number of Fisher Scoring iterations: 6
Important
If you forget family = binomial, R will silently fit a Gaussian (OLS) model on your binary outcome. Always check your model specification.
4.2 Comparing LPM and logit side by side
modelsummary(
list("LPM" = lpm, "Logit" = logit_model),
stars = TRUE,
gof_map = c("nobs", "r.squared", "logLik", "AIC"),
title = "Firm exit: LPM vs. logistic regression"
)| LPM | Logit | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | 0.217*** | -1.048** |
| (0.022) | (0.331) | |
| sales_log | -0.014*** | -0.169*** |
| (0.002) | (0.026) | |
| profit_margin | -0.088*** | -1.048*** |
| (0.006) | (0.086) | |
| employees_log | 0.006* | 0.006 |
| (0.003) | (0.043) | |
| Num.Obs. | 17506 | 17506 |
| R2 | 0.032 | |
| Log.Lik. | -556.349 | -4166.748 |
Note: the logit model reports log-likelihood and AIC rather than R². These are the relevant goodness-of-fit statistics for maximum likelihood models.
5 Interpretation
5.1 Level 1: Log-odds (raw coefficients)
coef(logit_model) (Intercept) sales_log profit_margin employees_log
-1.048478951 -0.168695919 -1.047629865 0.006126538 The coefficient on profit_margin is -1.048.
Interpretation: a one-unit increase in profit_margin is associated with a change of -1.048 in the log-odds of exit. This tells us the direction clearly, but the magnitude is not intuitive.
5.2 Level 2: Odds ratios
Exponentiate the coefficients to get odds ratios:
exp(coef(logit_model)) (Intercept) sales_log profit_margin employees_log
0.3504704 0.8447657 0.3507681 1.0061453 The odds ratio for profit_margin is 0.351.
Interpretation: a one-unit increase in profit margin multiplies the odds of exit by 0.351, i.e., it reduces the odds by 64.9%.
We can add confidence intervals:
exp(confint(logit_model)) 2.5 % 97.5 %
(Intercept) 0.1821600 0.6674216
sales_log 0.8029502 0.8893463
profit_margin 0.2961005 0.4147961
employees_log 0.9243338 1.09444165.3 Level 3: Predicted probabilities (most useful)
5.3.1 For observed data
firms <- firms |>
mutate(pred_prob = predict(logit_model, type = "response"))
firms |>
group_by(exit) |>
summarise(mean_pred_prob = mean(pred_prob))# A tibble: 2 × 2
exit mean_pred_prob
<dbl> <dbl>
1 0 0.0670
2 1 0.0994Firms that actually exited had, on average, a predicted exit probability of 0.099, compared to 0.067 for firms that survived. The model is doing something sensible.
5.3.2 For specific scenarios
This is often the most useful form for communicating results:
scenarios <- tibble(
label = c("Profitable, large", "Loss-making, small"),
sales_log = c(11, 9),
employees_log = c(4, 2),
profit_margin = c(0.10, -0.15)
)
scenarios <- scenarios |>
mutate(
pred_prob = predict(logit_model, newdata = scenarios, type = "response")
)
scenarios# A tibble: 2 × 5
label sales_log employees_log profit_margin pred_prob
<chr> <dbl> <dbl> <dbl> <dbl>
1 Profitable, large 11 4 0.1 0.0481
2 Loss-making, small 9 2 -0.15 0.0834A profitable, large firm has a predicted exit probability of 0.048, compared to 0.083 for a small, loss-making firm.
5.3.3 Visualising predicted probabilities
# How does exit probability vary with profit margin?
grid <- tibble(
profit_margin = seq(-0.5, 0.5, by = 0.01),
sales_log = mean(firms$sales_log, na.rm = TRUE),
employees_log = mean(firms$employees_log, na.rm = TRUE)
)
grid <- grid |>
mutate(pred_prob = predict(logit_model, newdata = grid, type = "response"))
ggplot(grid, aes(x = profit_margin, y = pred_prob)) +
geom_line(color = "#1a9641", linewidth = 1.2) +
geom_hline(yintercept = c(0, 1), linetype = "dashed", color = "grey60") +
labs(
title = "Predicted probability of exit vs. profit margin",
subtitle = "Other variables held at their mean",
x = "Profit margin", y = "P(exit = 1)"
) +
theme_minimal()
5.4 Average Marginal Effects
A cleaner summary of “the effect of X on P(exit)” is the Average Marginal Effect (AME). Because logit is non-linear, marginal effects differ across observations. The AME averages them:
avg_slopes(logit_model)
Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
employees_log 0.000381 0.00268 0.142 0.887 0.2 -0.00487 0.00563
profit_margin -0.065177 0.00546 -11.927 <0.001 106.5 -0.07589 -0.05447
sales_log -0.010495 0.00163 -6.445 <0.001 33.0 -0.01369 -0.00730
Type: response
Comparison: dY/dXThe AME for profit_margin tells us: on average across all firms, a one-unit increase in profit margin is associated with a -0.065 change in the probability of exit. This is directly comparable to the LPM coefficient.
6 Summary
- OLS on binary outcomes produces out-of-bounds predictions and invalid standard errors
- Logistic regression uses the sigmoid function to ensure \(\hat{P} \in [0, 1]\)
- Estimation is via maximum likelihood using
glm(..., family = binomial) - Raw coefficients are in log-odds space — hard to interpret directly
- Odds ratios \(= e^\beta\) — easier, but still require explanation
- Predicted probabilities — most useful for business communication; compute with
predict(..., type = "response") - Average marginal effects via
marginaleffects::avg_slopes()— a cleaner, directly comparable effect size - Inline R code ensures your text and numbers stay in sync
7 Further Reading
- Békés, G. & Kézdi, G. (2021). Data Analysis for Business, Economics, and Policy. Ch. 11. Available at gabors-data-analysis.com
- James, G., Witten, D., Hastie, T., & Tibshirani, R. (2021). An Introduction to Statistical Learning. Ch. 4 (Classification). Available at statlearning.com