So far, every model in this book has had a continuous outcome. Now we’re going to look at binary (yes/no) outcomes, and that changes things.
Since Chapter 4, every model has had the shape
\[y_i = (\text{prediction}) + \epsilon_i\]
where \(\epsilon_i\) is each person’s weirdness (how different they are from what the model expects). But that doesn’t really work anymore.
Suppose the model predicts .6 for someone whose answer is 1. Their “weirdness” would be .4. For someone else with the same prediction who answered 0, it would be -.6. Those are the only two possible values, and they’re completely determined by the prediction and the answer. So “weirdness” isn’t really something a person has anymore.
So we change what the model says. Instead of saying the outcome is
a systematic part, plus each person’s own departure from it
we say
a systematic part that determines the probability of each outcome, and then the outcome is a draw from that distribution
We don’t add anything to the prediction. The prediction is a probability, and the randomness comes from the Bernoulli distribution (see Chapter 1), not from an error term.
It turns out the “weirdness” term was never that fundamental. It works for the normal distribution because a normal outcome can be written as its mean plus some random noise. Without normality, that doesn’t work. The more general idea is that the model specifies a distribution for \(y\), and the predictors determine its parameters. The linear models we’ve been using so far were a special case of this all along!
Last time we looked at (in)dependence between two variables. We can also think of this as a simple model:
\[\text{admission} = f(\text{citizenship})\]
If we were only analyzing 2x2 tables, we could stop with tests of independence and odds ratios. But we want to go further. And it turns out the table already has a model in it!
Here’s the Berkeley table from Chapter 15 again, this time as conditional probabilities:
berk <- matrix(c(93, 4, 212, 33),
nrow = 2,
dimnames = list(admitted = c("no", "yes"),
citizenship = c("Other", "U.S.")))
round(prop.table(t(berk), margin = 1), 3) admitted
citizenship no yes
Other 0.959 0.041
U.S. 0.865 0.135Models are functions. Abstractly, \(y\) is some function of \(x\). That is, \(y\) is conditional on \(x\). In high school algebra, you were given the function and asked for the output. Here we have the opposite problem: we don’t have the model, we have the result of an unknown model.
Not \(y = 4x + 12\), but
\[ 16 = ? \times 1 + ? \qquad \text{and} \qquad 12 = ? \times 0 + ? \]
Substituting our conditional probabilities, with \(x = 0\) for other citizenship and \(x = 1\) for U.S.:
\[ .135 = ? \times 1 + ? \qquad \text{and} \qquad .041 = ? \times 0 + ? \]
So the intercept is .041 and the slope is .094. Here’s how we’ll write it from now on:
\[ \pi(y = 1 \mid x) = \alpha + \beta x \]
Choose one of the values of \(X\) to set to 0 (e.g., non-citizen) and one to set to 1 (e.g., citizen). Then \(\pi(y = 1 \mid x = 0) = \alpha\), because \(\beta \times 0\) drops out, and \(\pi(y = 1 \mid x = 1) = \alpha + \beta\). Substituting the conditional probabilities from the table gives \(.041 = \alpha\) and \(.135 = \alpha + \beta\), so \(\beta = .094\).
Remember that \(\alpha\) is a point and \(\beta\) is a distance.
As we discussed in Chapter 15, we don’t use the straight probabilities, but the logit (i.e., log odds). So the model becomes
\[ \log \left( \frac{\pi(y = 1 \mid x)}{1 - \pi(y = 1 \mid x)} \right) = \alpha + \beta x \]
We can do the same thing with the log odds of each row:
alpha <- log(4 / 93) # log odds of admission when x = 0
beta <- log(33 / 212) - log(4 / 93) # difference in log odds
c(alpha = alpha, beta = beta) alpha beta
-3.146305 1.286226 So to make this feel like high school, the data gives us the following model of admissions:
\[ \text{logit}(\text{admission}) = 1.29x - 3.15 \]
where \(x = 0\) if the applicant is not a citizen and \(x = 1\) if the applicant is a citizen.
And \(\exp(\beta)\) is the odds ratio we calculated by hand in Chapter 15:
c(exp_beta = exp(beta),
cross_product = (93 * 33) / (212 * 4)) exp_beta cross_product
3.619104 3.619104 Select one value of \(x\) to be 0 and the other 1. Calculate the log odds of the outcome for each row of the table. Then calculate \(\alpha\) (the log odds that \(y = 1\) when \(x = 0\)) and \(\beta\) (the log odds that \(y = 1\) when \(x = 1\), minus \(\alpha\)).
Country music and education:
country <- matrix(c(387, 209, 819, 157),
nrow = 2,
dimnames = list(BA = c("no", "yes"),
likes_country = c("no", "yes")))
addmargins(country) likes_country
BA no yes Sum
no 387 819 1206
yes 209 157 366
Sum 596 976 1572\(\alpha =\) ______ , \(\beta =\) ______ , \(\exp(\beta) =\) ______
a <- log(819 / 387) # log odds y = 1 when x = 0
b <- log(157 / 209) - a # diff. in log odds (x = 1 vs. x = 0)
c(alpha = a, beta = b, odds_ratio = exp(b)) alpha beta odds_ratio
0.7496594 -1.0357478 0.3549608 Basic interpretation: “The odds that someone with a college degree likes country music are only .35 times as great as for someone without a college degree.”
Let’s show how to do this using a model. Let’s get the data, including a couple of binary variables:
abany: 1 = agrees that a woman should be able to have an abortion for any reason, 0 otherwisenotv: a variable built from tvhours that = 1 if the person watches no TV, 0 otherwised <- readRDS(here::here("data", "gss2024.rds")) |>
haven::zap_labels() |>
select(sex, educ, age, abany, tvhours) |>
mutate(college = factor(if_else(educ >= 16, "Yes", "No")),
abany = if_else(abany == 1, 1, 0),
notv = if_else(tvhours == 0, 1, 0),
female = if_else(sex == 2, 1, 0)) |>
drop_na()
nrow(d)[1] 1014Let’s look at the crosstab of college and abany.
d |> tabyl(college, abany) |> adorn_totals(c("row", "col")) college 0 1 Total
No 283 326 609
Yes 141 264 405
Total 424 590 1014We can convert this into conditional probabilities.
d |> summarize(m_abany = mean(abany), .by = college)# A tibble: 2 x 2
college m_abany
<fct> <dbl>
1 No 0.535
2 Yes 0.652Remember that the GLM has two parts. The family is the probability distribution for the outcome, and for a binary outcome that’s always the binomial (Bernoulli). The link says what function of the outcome’s expected value is linear in the predictors. We get to choose the link, and there are three common ones for binary outcomes.
Identity link. This is the linear probability model, treating the outcome as continuous even when it’s really not.
mod_id <- glm(abany ~ college,
data = d,
family = binomial(link = "identity"))
tidy(mod_id)[, 1:3]# A tibble: 2 x 3
term estimate std.error
<chr> <dbl> <dbl>
1 (Intercept) 0.535 0.0202
2 collegeYes 0.117 0.0311\[ \beta = \pi(y = 1 \mid x = 1) - \pi(y = 1 \mid x = 0) \]
What is the unit of \(\beta\)? Compare these to the conditional probabilities above. \(\alpha\) is the probability for the reference group, and \(\beta\) is the difference in probabilities.
Log link. Now the model is \(\log(\pi_i) = \alpha + \beta x_i\), or equivalently \(\pi_i = e^{\alpha + \beta x_i}\).
mod_log <- glm(abany ~ college,
data = d,
family = binomial(link = "log"))
tidy(mod_log)[, 1:3]# A tibble: 2 x 3
term estimate std.error
<chr> <dbl> <dbl>
1 (Intercept) -0.625 0.0378
2 collegeYes 0.197 0.0524exp(coef(mod_log)["collegeYes"])collegeYes
1.217723 \(\beta\) here is a log risk ratio. Exponentiate it and you get the ratio of the two conditional probabilities, \(.652 / .535\).
Logit link. And finally the one this chapter is named for:
\[ \log \left( \frac{\pi_i}{1 - \pi_i} \right) = \alpha + \beta x_i \qquad \text{or} \qquad \pi_i = \frac{e^{\alpha + \beta x_i}}{1 + e^{\alpha + \beta x_i}} \]
m1 <- glm(abany ~ college,
data = d,
family = binomial(link = "logit"))
tidy(m1)[, 1:3]# A tibble: 2 x 3
term estimate std.error
<chr> <dbl> <dbl>
1 (Intercept) 0.141 0.0812
2 collegeYes 0.486 0.132 We can turn the college coefficient into an odds ratio by exponentiating it.
exp(coef(m1)[2])collegeYes
1.625375 | Term | Identity | Log | Logit |
|---|---|---|---|
| \(\alpha\) | probability | log prob. | log odds |
| \(\beta\) | diff. in prob. | log risk ratio | log odds ratio |
Linear probability models (identity link) and logistic regression models (logit link) are both common in sociology. The main value of the former is simplicity and the main value of the latter is realism. As we’ll see below, sometimes the LPM can create predicted probabilities below 0 or above 1, which is suboptimal!
Just like with the Berkeley table, a logistic regression coefficient is a log odds ratio. With one binary predictor, it’s exactly the log odds ratio from the 2x2 table. With more predictors, it’s the log odds ratio adjusting for the other predictors (like in Chapter 11).
An odds ratio of 1.63 does not mean college graduates are 1.63 times as likely to agree. Odds aren’t probabilities, and the odds ratio is only close to the risk ratio when the outcome is rare. Instead of making readers figure this out in their heads, we can get the quantities we want from marginaleffects.
d |>
summarize(m_abany = mean(abany), .by = college)# A tibble: 2 x 2
college m_abany
<fct> <dbl>
1 No 0.535
2 Yes 0.652# get difference (default)
avg_comparisons(m1, variables = "college")
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
0.117 0.0311 3.74 <0.001 12.4 0.0555 0.178
Term: college
Type: response
Comparison: Yes - No# get risk ratio, not odds ratio!
avg_comparisons(m1, variables = "college", comparison = "ratio")
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
1.22 0.0638 19.1 <0.001 267.4 1.09 1.34
Term: college
Type: response
Comparison: mean(Yes) / mean(No)Notice that the difference is 0.117, which is the \(\beta\) from the identity link model. And the ratio is 1.218, which is the exponentiated \(\beta\) from the log link model. All three models describe the same table, just on different scales. And avg_comparisons() lets you get whichever one you want after fitting the logit.
Odds ratios are what the model estimates, so it’s fine to put them in a table. But they aren’t a great way to communicate results, because hardly anyone thinks in odds, and people often misread them as probabilities. So report the average marginal effect on the probability scale too. It’s just one more line of code!
m2 <- glm(abany ~ college + age, # is this reasonable?
data = d,
family = binomial())
tidy(m2)[, 1:3]# A tibble: 3 x 3
term estimate std.error
<chr> <dbl> <dbl>
1 (Intercept) 0.617 0.197
2 collegeYes 0.492 0.133
3 age -0.00955 0.00360Because the model isn’t linear on the probability scale, the relationship with age isn’t exactly a straight line. We can see this by plotting the predicted probabilities.
curve_data <- predictions(
m2,
newdata = datagrid(age = 18:89, college = c("No", "Yes"))
) |> as.data.frame()
plt(estimate ~ age | college, data = curve_data, type = "l", lwd = 2,
ylim = c(0, 1),
xlab = "Age", ylab = "Predicted probability")
ggpredict(m2, terms = c("age [all]", "college")) |> plot()
And the average marginal effect of college, adjusting for age:
avg_slopes(m2, variables = "college")
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
0.117 0.031 3.78 <0.001 12.7 0.0565 0.178
Term: college
Type: response
Comparison: Yes - NoThis is about the same as the college difference without adjusting for age.
Let’s try a rarer outcome. About 9% of respondents say they don’t watch any television at all.
m3 <- glm(notv ~ college + age, # is this reasonable?
data = d,
family = binomial())
tidy(m3)[, 1:3]# A tibble: 3 x 3
term estimate std.error
<chr> <dbl> <dbl>
1 (Intercept) -1.58 0.324
2 collegeYes 0.299 0.222
3 age -0.0183 0.00646Let’s look at this one on the log odds scale. On this scale, the model is just two parallel straight lines (it’s “linear in the logit”).
link_data <- predictions(
m3,
newdata = datagrid(age = 18:89, college = c("No", "Yes")),
type = "link"
) |> as.data.frame()
plt(estimate ~ age | college, data = link_data, type = "l", lwd = 2,
xlab = "Age", ylab = "Log odds")
ggplot(link_data, aes(x = age, y = estimate, color = college)) +
geom_line(linewidth = 1) +
labs(x = "Age", y = "Log odds", color = "College")
This is the key thing to understand about logistic regression: it’s additive (parallel lines) on the log odds scale, but curved on the probability scale.
This also means that the marginal effects aren’t the same for everyone. slopes() gives one for each respondent:
slopes(m3, variables = "age") |>
as.data.frame() |>
head(4) |>
subset(select = c(age, college, estimate, std.error)) age college estimate std.error
1 19 No -0.002031940 0.0009932325
2 25 No -0.001870056 0.0008687024
3 63 No -0.001052089 0.0002955040
4 31 Yes -0.002155430 0.0009512308avg_slopes() is just the average of these. (Keep in mind that the average slope doesn’t necessarily describe any particular person.)
avg_comparisons(m3, variables = c("college", "age"))
Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
age +1 -0.00147 0.000523 -2.80 0.00505 7.6 -0.00249 -0.000442
college Yes - No 0.02474 0.018725 1.32 0.18639 2.4 -0.01196 0.061444
Type: responseThe college difference is small and its confidence interval includes zero. Now look at the odds ratio:
exp(coef(m3)["collegeYes"])collegeYes
1.349158 An odds ratio of 1.35 sounds pretty big (“35% higher odds”!), but the difference in probability is only about 2.5 percentage points and we can’t tell it apart from zero. When the outcome is rare, a big-sounding odds ratio can go along with a tiny difference in probability.
We don’t have an SSE here, so we can’t use the F test. Instead we use the tools from Chapter 9: deviance and the likelihood ratio test.
m0 <- glm(abany ~ 1, data = d, family = binomial())
anova(m0, m1, m2, test = "LRT")Analysis of Deviance Table
Model 1: abany ~ 1
Model 2: abany ~ college
Model 3: abany ~ college + age
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 1013 1378.4
2 1012 1364.7 1 13.6925 0.0002153 ***
3 1011 1357.7 1 7.0631 0.0078688 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Each line tests whether adding a parameter reduced the deviance more than we’d expect by chance. Both do.
Notice that all three models use the same 1,014 rows. That’s because we used drop_na() on all the variables at the start. If we had dropped missing values separately for each model, the models would be fit to different data and we couldn’t compare them.
Here’s one version of pseudo-\(R^2\) (McFadden’s):
1 - as.numeric(logLik(m2)) / as.numeric(logLik(m0))[1] 0.01505771This is pretty small, as pseudo-\(R^2\) values usually are. It’s not really comparable to an \(R^2\) from a linear model. Appendix C shows how different the various versions can be.
We just spent a whole chapter on logistic regression, so it might be surprising that I think the linear probability model is often the better choice for a lot of sociological problems.
Let’s fit both models with the same predictors and compare the LPM coefficients to the average marginal effects from the logistic model:
lpm_full <- lm(abany ~ college + age, data = d)
logit_full <- glm(abany ~ college + age, data = d, family = binomial())
ame <- as.data.frame(avg_comparisons(logit_full, variables = c("college", "age")))
data.frame(
term = c("college", "age"),
LPM = round(c(coef(lpm_full)["collegeYes"], coef(lpm_full)["age"]), 5),
logit_AME = round(c(ame$estimate[ame$term == "college"],
ame$estimate[ame$term == "age"]), 5),
row.names = NULL
) term LPM logit_AME
1 college 0.11751 0.11731
2 age -0.00230 -0.00228They’re almost identical. This usually happens, because in the middle range of probabilities the logistic curve is pretty close to a straight line.
And the LPM coefficients already are differences in probability. You don’t need to exponentiate anything or calculate marginal effects, and nobody can mistake an odds ratio for a risk ratio.
The usual objection to the LPM is that it can predict values below 0 or above 1. But think about what that tells you.
lpm_bad <- lm(notv ~ educ + age + I(age^2), data = d)
fitted_bad <- predict(lpm_bad)
c(minimum = min(fitted_bad),
maximum = max(fitted_bad),
share_outside = mean(fitted_bad < 0 | fitted_bad > 1)) minimum maximum share_outside
-0.10990297 0.14713108 0.02761341 Here, about 2.8% of the predicted values for people in our data are below zero. That’s a sign that the model is wrong. (We fit a quadratic in age to a rare outcome, and it’s predicting impossible values for real respondents.)
If we fit the same model with a logit, all the predictions will be between 0 and 1, and we’d never know there was a problem. So the thing people like most about logistic regression can actually hide a bad specification. The LPM, on the other hand, lets you know.
The LPM does have one real technical problem. Since the outcome is Bernoulli, its variance is \(\pi(1-\pi)\), which depends on the prediction. That means the errors can’t have constant variance, so the standard errors from lm() are wrong.
The fix is easy: use heteroskedasticity-consistent (“robust”) standard errors.
msummary(list("Classical SEs" = lpm_full, "Robust SEs" = lpm_full),
vcov = c("classical", "HC1"),
gof_map = c("nobs", "r.squared"))| Classical SEs | Robust SEs | |
|---|---|---|
| (Intercept) | 0.649 | 0.649 |
| (0.047) | (0.047) | |
| collegeYes | 0.118 | 0.118 |
| (0.031) | (0.031) | |
| age | -0.002 | -0.002 |
| (0.001) | (0.001) | |
| Num.Obs. | 1014 | 1014 |
| R2 | 0.020 | 0.020 |
The estimates are the same. Only the standard errors change (and here not by much). So if you report the robust standard errors, this problem goes away.
For a lot of the problems in this book (probabilities that aren’t too close to 0 or 1, simple predictors, describing patterns rather than making predictions), I think the LPM is a good default. Its coefficients are easy to interpret, it usually agrees with the logistic model’s marginal effects, and it tells you when something has gone wrong.
Logistic regression makes more sense when probabilities really are near 0 or 1, when the relationship is very nonlinear, when you need predictions for new cases that have to be valid probabilities, or when you want to use the likelihood tools from Chapter 9 to compare models.
But I don’t think “logistic regression is the correct model for a binary outcome” is a good reason. Both models are approximations, and only one of them tells you when it isn’t working.
For a fuller version of this argument, see Hellevik (2009).
avg_comparisons() recovers the difference or the risk ratio from a fitted logit