Skip to content
Statistical RethinkingStatistical Rethinking

10.1 Binomial regression

Binomial regression models binary and binomial outcomes using a logit link. The logit function maps the real line to the unit interval, ensuring that predicted probabilities are valid.

10.1.1 The logit link

For binary outcomes (yi∈{0,1}y_i \in \{0, 1\}) or binomial outcomes (yi∼Binomial(ni,pi)y_i \sim \text{Binomial}(n_i, p_i)):

logit(pi)=log⁡pi1−pi=α+βxi\text{logit}(p_i) = \log\frac{p_i}{1 - p_i} = \alpha + \beta x_i

pi=logistic(α+βxi)=11+e−(α+βxi)p_i = \text{logistic}(\alpha + \beta x_i) = \frac{1}{1 + e^{-(\alpha + \beta x_i)}}

The coefficients are on the log-odds scale. To interpret them, convert to odds ratios: eβe^\beta is the multiplicative change in odds per unit increase in xx.

# Binomial regression: chimpanzee prosociality example
# Does the prosocial option being on the left increase pulling?

set.seed(42)
n <- 200

# Simulate: prosocial on left increases probability of pulling left

prosoc_left <- sample(c(0, 1), n, replace = TRUE)
logit_p <- -0.5 + 0.8 * prosoc_left
p_true <- 1 / (1 + exp(-logit_p))
pulled_left <- rbinom(n, 1, p_true)

# Fit logistic regression

m <- glm(pulled_left ~ prosoc_left, family = binomial)

cat("Binomial regression (logit link):\n")
cat(" Intercept:", round(coef(m)[1], 3), "(log-odds when prosoc_left = 0)\n")
cat(" prosoc_left:", round(coef(m)[2], 3), "(log-odds change)\n")

# Convert to odds ratio

cat("\n Odds ratio:", round(exp(coef(m)[2]), 3), "\n")

# Convert to probabilities

p0 <- 1 / (1 + exp(-coef(m)[1]))
p1 <- 1 / (1 + exp(-(coef(m)[1] + coef(m)[2])))
cat("\n P(pull left | prosoc_right) =", round(p0, 3), "\n")
cat(" P(pull left | prosoc_left) =", round(p1, 3), "\n")

# Plot

plot(jitter(prosoc_left, 0.1), jitter(pulled_left, 0.1),
pch = 16, col = adjustcolor("steelblue", 0.5),
xlab = "Prosocial on left", ylab = "Pulled left",
main = "Binomial regression: prosociality")
x_seq <- seq(0, 1, length.out = 100)
logit_pred <- coef(m)[1] + coef(m)[2] * x_seq
p_pred <- 1 / (1 + exp(-logit_pred))
lines(x_seq, p_pred, lwd = 2, col = "firebrick")
Click Run to execute this code in your browser with WebR

10.1.2 Aggregated binomial

When each observation has ni>1n_i > 1 trials, we use the aggregated binomial:

yi∼Binomial(ni,pi)y_i \sim \text{Binomial}(n_i, p_i)

# Aggregated binomial regression
set.seed(42)
n_groups <- 20

# Each group has a different number of trials

n_trials <- sample(10:50, n_groups, replace = TRUE)
x <- runif(n_groups, 0, 10)

# True relationship

logit_p <- -1 + 0.3 * x
p_true <- 1 / (1 + exp(-logit_p))
successes <- rbinom(n_groups, n_trials, p_true)

# Fit aggregated binomial regression

m <- glm(cbind(successes, n_trials - successes) ~ x, family = binomial)

cat("Aggregated binomial regression:\n")
cat(" Intercept:", round(coef(m)[1], 3), "\n")
cat(" Slope:", round(coef(m)[2], 3), "\n")

# Predicted probabilities

cat("\nPredicted probabilities:\n")
for (x_val in c(2, 5, 8)) {
p <- 1 / (1 + exp(-(coef(m)[1] + coef(m)[2] * x_val)))
cat(sprintf(" x = %d: p = %.3f\n", x_val, p))
}

# Plot

x_seq <- seq(0, 10, length.out = 100)
logit_pred <- coef(m)[1] + coef(m)[2] * x_seq
p_pred <- 1 / (1 + exp(-logit_pred))

plot(x, successes / n_trials, pch = 16, col = "steelblue",
xlab = "x", ylab = "Proportion of successes",
main = "Aggregated binomial regression")
lines(x_seq, p_pred, lwd = 2, col = "firebrick")
Click Run to execute this code in your browser with WebR
Rethinking

Log-odds are not intuitive. The logit link maps probabilities to the real line, but log-odds are hard to interpret directly. A coefficient of 0.8 does not mean "80% increase" — it means the log-odds increase by 0.8, which corresponds to an odds ratio of e0.8≈2.2e^{0.8} \approx 2.2. Always convert to the probability scale for interpretation. Plot the predicted probabilities at meaningful values of the predictors. The model lives on the logit scale, but the meaning lives on the probability scale.

Summary of Section 10.1

  • Binomial regression uses the logit link: log⁡(p/(1−p))=α+βx\log(p/(1-p)) = \alpha + \beta x.
  • Coefficients are log-odds — convert with eβe^\beta for odds ratios.
  • For aggregated data, use glm(cbind(successes, failures) ~ x, family = binomial).
  • Always interpret on the probability scale — plot predicted probabilities.