Logistic Regression
Marc Los Huertos
EA030 · Compiled 2026-04-14
By the end of this module you will be able to
-
Explain why ordinary regression fails for binary outcomes and what
logistic regression does instead.
-
Fit
glm(…, family = binomial) in R and interpret model
output.
-
Convert log-odds to odds ratios and to predicted probabilities.
-
Evaluate model fit with AIC, the null deviance, and a confusion matrix.
-
Visualise a logistic curve over the observed data.
What this module covers. When your outcome is binary
— present/absent, exceeded/not exceeded, survived/died — ordinary linear
regression fails because it can predict probabilities outside [0, 1].
Logistic regression solves this by modeling the log-odds of the
outcome, producing well-calibrated probabilities and interpretable odds
ratios. This tutorial covers the full workflow: fitting with
glm(), interpreting coefficients as odds ratios, assessing
model fit, and evaluating predictive performance with ROC curves.
1. Why Not Just Use Linear Regression?
Suppose you want to predict whether a river site exceeds a pollution
threshold (1 = exceeds, 0 = does not) based on upstream land-use
intensity. If you fit a straight line through 0s and 1s, two problems
arise immediately:
- Out-of-range predictions. A line can predict
probabilities below 0 or above 1, which are nonsensical.
- Non-constant variance. The variance of a binary
outcome is \(p(1-p)\), which changes
with p — violating the homoscedasticity assumption of linear
regression.
Logistic regression maps the linear predictor to a probability using
the logistic (sigmoid) function, guaranteeing
predictions in (0, 1).
2. The Logistic Regression Model
The model links the probability of the outcome p to predictors
through the logit (log-odds) transformation:
\[\log\left(\frac{p}{1-p}\right) = \beta_0
+ \beta_1 x_1 + \cdots + \beta_k x_k\]
Solving for p:
\[p = \frac{1}{1 + e^{-(\beta_0 + \beta_1
x)}} = \text{logistic}(\beta_0 + \beta_1 x)\]
The coefficients \(\beta\) are
estimated by maximum likelihood (not least squares).
The log-likelihood measures how well the model’s predicted probabilities
match the observed 0/1 outcomes.
3. Fitting Logistic Regression in R
We use a dataset of 189 low-birth-weight births from
MASS::birthwt. The outcome is whether birth weight was
below 2.5 kg (low: 1 = low, 0 = normal).
library(MASS)
data(birthwt)
# Recode for clarity
bw <- birthwt |>
mutate(
low = factor(low, labels = c("Normal", "Low")),
smoke = factor(smoke, labels = c("Non-smoker", "Smoker")),
race = factor(race, labels = c("White", "Black", "Other")),
ht = factor(ht, labels = c("No hypertension", "Hypertension")),
low_bin = as.integer(low == "Low") # keep numeric version for glm
)
table(bw$low)
##
## Normal Low
## 130 59
Explore predictors first
bw |>
count(smoke, low) |>
group_by(smoke) |>
mutate(pct = n / sum(n)) |>
filter(low == "Low") |>
ggplot(aes(x = smoke, y = pct, fill = smoke)) +
geom_col(width = 0.5) +
scale_fill_manual(values = c("#4a7c5e", "#c8832a")) +
scale_y_continuous(labels = scales::percent) +
labs(x = NULL, y = "Proportion with low birth weight",
title = "Low birth weight by smoking status") +
theme_minimal(base_size = 13) +
theme(legend.position = "none")
Fit the model
mod_logit <- glm(low_bin ~ smoke + age + lwt + race + ht,
data = bw,
family = binomial(link = "logit"))
The key argument is family = binomial(link = "logit").
This tells glm() that: (1) the outcome is binary (Bernoulli
distribution), and (2) the link function is the logit, connecting the
linear predictor to the probability.
4. Reading summary() Output
summary(mod_logit)
##
## Call:
## glm(formula = low_bin ~ smoke + age + lwt + race + ht, family = binomial(link = "logit"),
## data = bw)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.777274 1.160336 0.670 0.50294
## smokeSmoker 1.064664 0.389159 2.736 0.00622 **
## age -0.021166 0.034960 -0.605 0.54489
## lwt -0.017297 0.006866 -2.519 0.01176 *
## raceBlack 1.238728 0.526253 2.354 0.01858 *
## raceOther 0.919052 0.427195 2.151 0.03145 *
## htHypertension 1.732974 0.687931 2.519 0.01177 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 207.88 on 182 degrees of freedom
## AIC: 221.88
##
## Number of Fisher Scoring iterations: 4
The coefficients table has the same columns as lm() —
Estimate, Std. Error, z value, Pr(>|z|) — but the test statistic is a
z-test (Wald test), not a t-test.
Deviance. The logistic model uses deviance
(−2 × log-likelihood) instead of residual sums of squares. The Null
deviance is for a model with only an intercept; the Residual deviance is
for your model. A substantial drop indicates the predictors explain
something.
AIC. Akaike Information Criterion = deviance + 2k.
Smaller is better; use AIC to compare competing models.
⚠️ Coefficients are on the log-odds scale A
coefficient of 0.72 for smokerSmoker does not mean smoking
increases the probability by 72 percentage points. It means the log-odds
increase by 0.72. To get the odds ratio, exponentiate: exp(0.72) ≈ 2.05.
To get predicted probabilities, use predict(mod, type =
“response”).
5. Odds Ratios
Odds ratios are the most common effect size for logistic regression.
OR > 1 means the predictor increases the odds of the outcome; OR <
1 means it decreases the odds.
# Exponentiate coefficients and confidence intervals
or_table <- exp(cbind(OR = coef(mod_logit), confint(mod_logit)))
round(or_table, 3)
## OR 2.5 % 97.5 %
## (Intercept) 2.176 0.232 22.413
## smokeSmoker 2.900 1.373 6.372
## age 0.979 0.913 1.048
## lwt 0.983 0.969 0.995
## raceBlack 3.451 1.235 9.888
## raceOther 2.507 1.099 5.922
## htHypertension 5.657 1.523 23.882
or_df <- as.data.frame(or_table[-1, ]) # drop intercept
or_df$term <- rownames(or_df)
names(or_df)[2:3] <- c("ci_lo", "ci_hi")
ggplot(or_df, aes(x = term, y = OR)) +
geom_point(size = 3, color = "#4a7c5e") +
geom_errorbar(aes(ymin = ci_lo, ymax = ci_hi),
width = 0.2, color = "#4a7c5e") +
geom_hline(yintercept = 1, linetype = "dashed", color = "#c8832a") +
coord_flip() +
scale_y_log10() +
labs(x = NULL, y = "Odds Ratio (log scale)",
title = "Logistic Regression: Odds Ratios and 95% CI") +
theme_minimal(base_size = 12)
6. Predicted Probabilities
Unlike linear regression, logistic regression predictions are
probabilities, not raw response values. Use
type = "response" with predict().
# Predictions across maternal weight range, by smoking status
new_data <- expand.grid(
lwt = seq(80, 250, by = 2),
smoke = factor(c("Non-smoker", "Smoker"), levels = levels(bw$smoke)),
age = median(bw$age),
race = "White",
ht = "No hypertension"
)
new_data$pred_prob <- predict(mod_logit, newdata = new_data, type = "response")
ggplot(new_data, aes(x = lwt, y = pred_prob, color = smoke)) +
geom_line(linewidth = 1.2) +
scale_color_manual(values = c("#4a7c5e", "#c8832a")) +
labs(x = "Maternal weight (lbs)", y = "P(low birth weight)",
title = "Predicted probability of low birth weight",
color = "Smoking status") +
theme_minimal(base_size = 13)
7. Model Diagnostics
Deviance goodness-of-fit
The likelihood ratio test compares your model to the null
(intercept-only) model.
mod_null <- glm(low_bin ~ 1, data = bw, family = binomial)
anova(mod_null, mod_logit, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: low_bin ~ 1
## Model 2: low_bin ~ smoke + age + lwt + race + ht
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 188 234.67
## 2 182 207.88 6 26.794 0.0001582 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pseudo-R² (McFadden)
pseudo_r2 <- 1 - logLik(mod_logit) / logLik(mod_null)
McFadden pseudo-R² = 0.114.
McFadden values of 0.2–0.4 indicate good fit; above 0.4 is excellent.
This is not directly comparable to OLS R² — it measures relative
improvement over the null model.
Binned residual plot
# Manual binned residual plot (avoids arm package dependency)
pred <- predict(mod_logit, type = "response")
resid <- bw$low_bin - pred
n_bins <- 20
breaks <- quantile(pred, probs = seq(0, 1, length.out = n_bins + 1))
bin_id <- cut(pred, breaks = breaks, include.lowest = TRUE, labels = FALSE)
bin_df <- data.frame(pred, resid, bin_id) |>
group_by(bin_id) |>
summarise(mean_pred = mean(pred), mean_resid = mean(resid),
se = sd(resid) / sqrt(n()), .groups = "drop")
ggplot(bin_df, aes(x = mean_pred, y = mean_resid)) +
geom_ribbon(aes(ymin = -2 * se, ymax = 2 * se), alpha = 0.15, fill = "#4a7c5e") +
geom_point(color = "#1e3a2f", size = 2.5) +
geom_hline(yintercept = 0, linetype = "dashed", color = "#c8832a") +
labs(x = "Predicted probability", y = "Average residual",
title = "Binned residual plot") +
theme_minimal(base_size = 12)
8. ROC Curve and AUC
The receiver operating characteristic (ROC) curve
shows how well the model discriminates between the two outcome classes
across all possible probability thresholds. The area under the
curve (AUC) summarizes discrimination in a single number: 0.5 =
random guessing; 1.0 = perfect discrimination.
# Compute sensitivity and specificity at each threshold
thresholds <- seq(0, 1, length.out = 300)
pred_prob <- predict(mod_logit, type = "response")
obs <- bw$low_bin
roc_df <- data.frame(
threshold = thresholds,
sens = sapply(thresholds, function(t) mean((pred_prob >= t)[obs == 1])),
spec = sapply(thresholds, function(t) mean((pred_prob < t)[obs == 0]))
)
# AUC via trapezoidal rule
auc <- abs(sum(diff(1 - roc_df$spec) * roc_df$sens[-nrow(roc_df)]))
ggplot(roc_df, aes(x = 1 - spec, y = sens)) +
geom_line(color = "#1e3a2f", linewidth = 1.2) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed",
color = "#c8832a", linewidth = 0.8) +
annotate("text", x = 0.7, y = 0.2,
label = paste0("AUC = ", round(auc, 3)),
color = "#1e3a2f", size = 4.5, fontface = "bold") +
labs(x = "1 − Specificity (False Positive Rate)",
y = "Sensitivity (True Positive Rate)",
title = "ROC Curve — Low Birth Weight Model") +
coord_equal() +
theme_minimal(base_size = 12)
Choosing a classification threshold
The default threshold of 0.5 is not always optimal. A confusion
matrix at 0.5 shows the trade-off between sensitivity (catching true
positives) and specificity (avoiding false positives).
pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
table(Predicted = pred_class, Observed = obs)
## Observed
## Predicted 0 1
## 0 122 43
## 1 8 16
💡 When to adjust the threshold In public health
contexts, false negatives (missing a high-risk case) may be more costly
than false positives. Lowering the threshold increases sensitivity at
the cost of specificity. The optimal threshold depends on the relative
costs of the two error types — a decision that requires domain
knowledge, not statistics alone.
🔑 Key takeaways
-
Use
glm(y ~ x, family = binomial) for binary outcomes.
Linear regression will produce nonsensical probability predictions for
extreme x values.
-
Coefficients are on the log-odds scale. Always
exponentiate with
exp(coef(mod)) to get interpretable
odds ratios.
-
OR > 1 increases the odds; OR < 1 decreases them. An OR of 2 means
the odds double for a one-unit increase in x.
-
Use
predict(type = “response”) to get predicted
probabilities.
-
Evaluate discrimination with AUC: above 0.7 is
acceptable; above 0.8 is good. Check calibration with a binned residual
plot.
What’s next?
Exercise 1 — Fit and interpret
Using airquality, create a binary variable
high_ozone (TRUE if Ozone > 60 ppb). Fit a logistic
regression predicting high_ozone from Temp.
(a) Report the coefficient on Temp as a log-odds and as an odds ratio.
(b) What is the predicted probability of a high-ozone day at 90°F?
Show solution
data(airquality)
airquality$high_ozone <- airquality$Ozone > 60
mod <- glm(high_ozone ~ Temp, data = airquality, family = binomial, na.action = na.omit)
summary(mod)
exp(coef(mod)) # odds ratios
predict(mod, newdata = data.frame(Temp = 90), type = "response")
Each 1°F increase in temperature multiplies the odds of a high-ozone day
by OR ≈ 1.15. At 90°F, predicted probability ≈ 0.35 (35%).
