Before diving into logistic regression, it helps to think carefully about what we’re trying to model when we have two binary variables. The fundamental question is: are they associated? And if so, how do we measure the strength of that association?
Chapter 3 built marginal, joint, and conditional probabilities out of a contingency table and left the story there. This chapter finishes it.
We’ll work through an example: PhD sociology admissions at Berkeley in 2009, broken down by citizenship status.
berk <- matrix(c(93, 4, 212, 33),
nrow = 2,
dimnames = list(admitted = c("no", "yes"),
citizenship = c("Other", "U.S.")))
addmargins(berk) citizenship
admitted Other U.S. Sum
no 93 212 305
yes 4 33 37
Sum 97 245 342The marginal probability of admission is the overall rate, ignoring citizenship:
p_admit <- rowSums(berk)["yes"] / sum(berk)
p_admit yes
0.1081871 Similarly, \(P(\text{U.S.}) = 245/342 = .716\). These are “marginal” because they come from the margins of the table—they don’t condition on anything.
More interesting are conditional probabilities—the probability of admission given citizenship status:
p_admit_us <- berk["yes", "U.S."] / colSums(berk)["U.S."]
p_admit_other <- berk["yes", "Other"] / colSums(berk)["Other"]
c(US = p_admit_us, Other = p_admit_other) US.U.S. Other.Other
0.13469388 0.04123711 U.S. applicants are admitted at 13.5%; non-U.S. applicants at 4.1%. You’ll often see this written as \(\pi_{\text{admitted}|\text{U.S.}} = .135\) or \(\Pr(Y = 1 \mid X = 1) = .135\).
This can also be flipped to ask a different question. What is \(\Pr(X = 1 \mid Y = 1)\), the probability that someone who is admitted is a U.S. citizen? And the probability that someone who is not admitted is a U.S. citizen?
c(admitted = berk["yes", "U.S."] / rowSums(berk)["yes"],
rejected = berk["no", "U.S."] / rowSums(berk)["no"])admitted.yes rejected.no
0.8918919 0.6950820 Neither is the “right” conditional probability; which one you want depends on what you are asking.
The key question underneath all of this: is the shape of the outcome distribution conditional on the values of some other variable?
Two variables are independent if the conditional probabilities equal the marginal:
\[ \pi_{\text{admitted}|\text{U.S.}} = \pi_{\text{admitted}|\text{Other}} = \pi_{\text{admitted}} \]
We can think about this in terms of expected cell counts: if the two variables were unrelated, how many people would we expect in each cell?
\[ E_{ij} = \frac{(\text{row total}_i)(\text{col total}_j)}{n} \]
outer(rowSums(berk), colSums(berk)) / sum(berk) Other U.S.
no 86.50585 218.49415
yes 10.49415 26.50585How many of those four numbers did you actually have to calculate? Only one. Once you know a single cell, the margins fix the other three. That is the sense in which a 2x2 table has one degree of freedom, and it will come back when we test for independence below.
The observed counts depart from these expected values, which tells us there’s some association. But how much? Four measures below.
\[ \pi_{\text{admitted}|\text{U.S.}} - \pi_{\text{admitted}|\text{Other}} = .135 - .041 = .094 \]
p_admit_us - p_admit_other U.S.
0.09345676 Simple and interpretable—about a 9 percentage-point gap. Greater difference means greater departure from independence.
The problem is the same difference means different things at different baselines. Going from .001 to .020 is a much bigger deal than going from .500 to .519, even though both are “.019 differences.”
Instead of differencing, we can divide the probabilities. This is called the relative risk ratio.
\[ \text{RRR} = \frac{\pi_{\text{admitted}|\text{U.S.}}}{\pi_{\text{admitted}|\text{Other}}} = \frac{.135}{.041} = 3.3 \]
p_admit_us / p_admit_other U.S.
3.266327 This is more informative than a raw difference when probabilities are far from .5. But the RRR for rejection—keeping non-U.S. applicants in the numerator for consistency—is only:
(1 - p_admit_other) / (1 - p_admit_us) Other
1.108004 The same association (\(3.3\times\) vs \(1.1\times\)) looks very different depending on whether you frame it as success or failure. This asymmetry is a real problem for interpretation.
First, what are odds? One in six is the probability of rolling a 6. One to five are the odds of rolling a 6. The odds of an event is \(\pi / (1 - \pi)\): probability of success divided by probability of failure. If you’ve ever heard “the odds are 5 to 1 against,” that’s what it means.
odds_us <- p_admit_us / (1 - p_admit_us)
odds_other <- p_admit_other / (1 - p_admit_other)
OR <- odds_us / odds_other
OR U.S.
3.619104 A U.S. applicant’s odds of admission are .156, or about 6.4 to 1 against. A non-U.S. applicant’s are .043, or about 23.4 to 1 against. The ratio of those two odds is 3.6.
There’s a convenient shortcut for 2×2 tables using cell counts directly—the product of the two diagonals divided by the product of the other two:
(berk["no", "Other"] * berk["yes", "U.S."]) /
(berk["no", "U.S."] * berk["yes", "Other"])[1] 3.619104The odds ratio solves the success-failure asymmetry: the OR for rejection is also 3.6.1 But the OR is still asymmetric in one way: it ranges from 0 to \(\infty\), and the same relationship looks like 3.6 or 0.28 depending on which group you put in the numerator.
1 Odds of rejection for non-U.S.: \((.959/.041) = 23.4\). For U.S.: \((.865/.135) = 6.4\). Ratio: \(23.4/6.4 = 3.6\).
Taking the natural log solves the remaining problems:
log(OR) U.S.
1.286226 The log odds ratio:
c(log_OR = log(OR), reversed = log(1 / OR)) log_OR.U.S. reversed.U.S.
1.286226 -1.286226 This is why logistic regression is logistic regression—it models the log odds of the outcome as a linear function of predictors. Every coefficient you estimate is, at heart, a log odds ratio. The key transformation is the logit:
\[ \text{logit}(\pi) = \log\left(\frac{\pi}{1 - \pi}\right) \]
p_grid <- data.frame(p = seq(0.01, 0.99, by = 0.001))
p_grid$log_odds <- log(p_grid$p / (1 - p_grid$p))plt(log_odds ~ p, data = p_grid, type = "l", lwd = 2,
xlab = "Probability", ylab = "Log odds",
main = "The logit transformation")
abline(h = 0, lty = 2, col = "gray60")
abline(v = 0.5, lty = 2, col = "gray60")
ggplot(p_grid, aes(p, log_odds)) +
geom_line() +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray60") +
geom_vline(xintercept = 0.5, linetype = "dashed", color = "gray60") +
labs(x = "Probability", y = "Log odds",
title = "The logit transformation")
The curve is symmetric around \(p = .5\), where \(\text{logit}(.5) = 0\), and it maps the bounded interval \((0, 1)\) to the full real line. That stretching is exactly what a linear model needs: a probability cannot be modeled directly as \(b_0 + b_1 x\), because a straight line will eventually predict values below 0 or above 1. Log odds have no such boundaries. We take that up in Chapter 16.
| Measure | Formula | Value | Independence |
|---|---|---|---|
| Difference | \(\pi_1 - \pi_2\) | \(.094\) | \(0\) |
| RRR | \(\pi_1 / \pi_2\) | \(3.3\) | \(1\) |
| Odds ratio | \((\pi_1/\bar\pi_1) / (\pi_2/\bar\pi_2)\) | \(3.6\) | \(1\) |
| Log odds ratio | \(\log(\text{OR})\) | \(1.28\) | \(0\) |
Chapter 17 will build models of the dependence between variables. Stand-alone tests are common, however, and you should be able to read them.
If you’ve taken any stats before, you’ve probably seen the good old chi-square test. You generally use it to analyze categorical data (where none of the variables involved is continuous). It asks: are these two variables independent? It calculates this as follows:
\[ \chi^2 = \sum \frac{(O_i - E_i)^2}{E_i} \]
where \(i\) indexes table cells, \(O\) means “observed,” and \(E\) means “expected values if independent.” Or in English: square the observed-expected difference and divide by the expected frequency for each cell, then add up all these values.
We can do this easily in R.
chisq.test(berk, correct = FALSE)
Pearson's Chi-squared test
data: berk
X-squared = 6.2905, df = 1, p-value = 0.01214Fisher’s exact test answers the same question a different way. For all possible tables with the same marginals, we ask: how many of these give at least as much evidence of association as the observed table? The p-value is the proportion of hypothetical tables where the association is at least as strong.
fisher.test(berk)
Fisher's Exact Test for Count Data
data: berk
p-value = 0.01146
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
1.232012 14.416963
sample estimates:
odds ratio
3.608676 Fisher’s test is the safer choice when cells are small, as they are here. Only four non-U.S. applicants were admitted.
Tables need not be 2x2. Let’s look at the contingency table of degree and sex in the 2024 GSS. We will use the janitor::tabyl() function as a more capable version of the base R table() function.
d <- readRDS(here::here("data", "gss2024.rds")) |>
haven::zap_labels() |>
select(tvhours, degree, sex) |>
drop_na() |>
mutate(sex = factor(sex,
labels = c("M", "F")), # sex as a factor
degree_fac = factor(degree,
labels = c("none", "HS", "AA", "BA", "grad")))
xtab <- d |>
tabyl(sex, degree_fac)
xtab |> adorn_totals(c("row", "col")) sex none HS AA BA grad Total
M 66 471 66 212 130 945
F 117 516 120 260 186 1199
Total 183 987 186 472 316 2144xtab |>
chisq.test(tabyl_results = TRUE)
Pearson's Chi-squared test
data: xtab
X-squared = 16.893, df = 4, p-value = 0.002027This doesn’t seem like a “model”—because it’s not. There is no coefficient, no effect size, no prediction, only a verdict on independence. But it turns out we can make it into a model. We’ll need to first consider a new probability distribution, which is Chapter 17.
Here’s a table from the GSS on whether spanking is acceptable, by sex:
spank <- matrix(c(461, 226, 489, 132),
nrow = 2,
dimnames = list(spank_ok = c("yes", "no"),
sex = c("Male", "Female")))
addmargins(spank) sex
spank_ok Male Female Sum
yes 461 489 950
no 226 132 358
Sum 687 621 1308Calculate all four measures of dependence. Which sex is more likely to say spanking is OK, and by how much? Specifically, work out: