library(tidyverse)
library(tidystata)
gss24 <- readRDS("gss24.rds")25 Logit and probit
How is parental schooling associated with the probability of completing a bachelor’s or graduate degree? Chapter 24 described years of schooling with OLS. Here we retain the same respondents and model a binary outcome, then compare the two models on a probability scale.
By the end of this chapter, you will be able to:
- Fit logit and probit models for a binary outcome and interpret their coefficient scales.
- Recreate and describe a common complete-case sample.
- Check predictor support and compare conditional probability predictions.
- State the limits of the sample, model, and interpretation.
Use the table operations from the practical chapters, factors from Chapter 15, and model formulas from Chapter 24. Prepare the saved extract using Appendix A. The setup is repeated here so the chapter runs in a fresh session. stats:: and broom:: identify functions from their packages without attaching them.
25.1 Restore the shared sample
Recreate the sample from Section 24.1. Require recorded ID, schooling, parent schooling, age, sex, and degree, so these models use the same respondents as the preceding OLS comparison. The source variable educ records years of schooling; maeduc records mother’s or younger-parent schooling under the appendix’s parent definition.
The indicator college is 1 for a bachelor’s or graduate degree and 0 for another recorded degree, including associate/junior college. We exclude unknown degrees before using %in%, so an unknown answer cannot become a known zero. The sex factor lists male first to retain the same reference category. The source extract remains unchanged.
model_data <- gss24 |>
select(id, educ, maeduc, age, sex, degree) |>
filter(!is.na(id), !is.na(educ), !is.na(maeduc), !is.na(age),
!is.na(sex), !is.na(degree)) |>
mutate(
sex = factor(sex, levels = c("male", "female")),
college = if_else(degree %in% c("bachelor's", "graduate"), 1, 0)
)
nrow(gss24)[1] 3986nrow(model_data)[1] 3329summ(model_data, educ, maeduc, age)# A tibble: 3 × 8
name type n unique min mean sd max
<chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
1 educ dbl 3329 20 0 14.5 2.82 20
2 maeduc dbl 3329 21 0 12.5 3.44 20
3 age dbl 3329 72 18 49.0 17.5 89
model_description <- model_data |>
summarise(n = n(), mean_educ = mean(educ),
mean_maeduc = mean(maeduc),
college_percent = 100 * mean(college), .by = sex)
model_description# A tibble: 2 × 5
sex n mean_educ mean_maeduc college_percent
<fct> <int> <dbl> <dbl> <dbl>
1 male 1484 14.5 12.7 39.8
2 female 1845 14.4 12.3 39.3
The analysis retains 3329 of the extract’s 3986 respondents. Each row remains one respondent; the summary reports schooling and age in years. The sex-specific rows show the group sizes, mean respondent and parent schooling, and college-attainment percentages within this shared sample. A mean of the 0/1 college indicator is a proportion, multiplied by 100 for a percentage. There is no need for na.rm because the required values are observed.
This restriction makes comparisons transparent, but it also excludes respondents who would be usable in an individual model. Complete-case analysis does not repair missing values or guarantee that included respondents resemble excluded ones. We retain the source table unchanged and make no population correction by restricting it.
All models below are unweighted teaching examples with ordinary model-based uncertainty. They do not account for the GSS sampling design. Its wtssps variable is a post-stratification weight; putting it into lm(weights = ...) would not, by itself, supply survey-design inference. We do not do that here.
25.2 Model a binary outcome and predict probabilities
College attainment takes values zero and one. We now model its probability on the same model_data rows. stats::glm() fits a generalized linear model. family = stats::binomial(link = "logit") specifies a binary outcome with the logit link, which relates a linear combination of the explanatory variables to the log odds of college attainment. For a probability \(p\), odds are \(p/(1-p)\), the ratio of the probability of college attainment to that of non-attainment; log odds are the natural logarithm of that ratio. A logit coefficient is a difference in log odds per unit of a predictor, conditional on the other terms. It is not a percentage-point difference.
logit_model <- stats::glm(
college ~ maeduc + age + sex, data = model_data,
family = stats::binomial(link = "logit")
)
logit_tidy <- broom::tidy(logit_model)
logit_tidy# A tibble: 4 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -3.55 0.235 -15.1 8.92e-52
2 maeduc 0.186 0.0131 14.1 3.19e-45
3 age 0.0151 0.00219 6.86 6.67e-12
4 sexfemale 0.0499 0.0740 0.674 5.00e- 1
The positive parent-schooling coefficient means that the fitted probability increases with parent schooling at the same age and sex. The output’s statistic is a large-sample z statistic for this GLM, rather than the OLS t statistic. We have changed both the outcome and the model, so do not compare the raw schooling-year OLS coefficients with these log-odds coefficients.
A probit model uses the same binary outcome and explanatory terms with link = "probit". Its link is the inverse of the standard normal cumulative distribution, so its coefficients use a different latent scale. Coefficient magnitudes cannot be compared directly with logit magnitudes. Predictions expressed as probabilities are easier to compare.
probit_model <- stats::glm(
college ~ maeduc + age + sex, data = model_data,
family = stats::binomial(link = "probit")
)
probit_tidy <- broom::tidy(probit_model)
probit_tidy# A tibble: 4 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -2.11 0.137 -15.4 9.68e-54
2 maeduc 0.109 0.00762 14.3 1.36e-46
3 age 0.00906 0.00133 6.79 1.11e-11
4 sexfemale 0.0296 0.0454 0.653 5.14e- 1
stats::nobs(logit_model)[1] 3329stats::nobs(probit_model)[1] 3329Both binary models retain the same 3329 respondents as the OLS comparison. They assume the chosen link and additive predictor form are adequate, with independent binary observations for the ordinary model-based inference. Changing the link does not repair missingness, survey dependence, or omitted explanations.
For a concrete comparison, consider female respondents aged 50 whose mother/younger parent completed 12 or 16 years of school. Check that both complete predictor combinations occur in our sample with filter() and count(). Then construct exactly those two rows with tibble(). The two-element columns are aligned by position, and the sex factor uses the model’s level order.
prediction_support <- model_data |>
filter(age == 50, sex == "female", maeduc %in% c(12, 16)) |>
count(maeduc, age, sex)
prediction_support# A tibble: 2 × 4
maeduc age sex n
<dbl> <dbl> <fct> <int>
1 12 50 female 11
2 16 50 female 2
prediction_cases <- tibble(
maeduc = c(12, 16),
age = c(50, 50),
sex = factor(c("female", "female"), levels = c("male", "female"))
)The two combinations occur for 11 and 2 respondents, respectively. The smaller cell is a reminder that observed support can still be sparse. We are not estimating separate means from just those two cells: each fitted model uses the entire common sample.
stats::predict() takes a fitted model and newdata, the table of predictor values. For these GLMs, type = "response" returns probabilities between zero and one rather than values on the link scale. We multiply by 100 to display percentages and store a separate output table, leaving the prediction inputs unchanged. knitr::kable() formats it, with digits = 1 rounding displayed numeric values to one decimal place.
prediction_table <- prediction_cases |>
mutate(
logit_percent = 100 * stats::predict(
logit_model, newdata = prediction_cases, type = "response"
),
probit_percent = 100 * stats::predict(
probit_model, newdata = prediction_cases, type = "response"
)
)
knitr::kable(prediction_table, digits = 1)| maeduc | age | sex | logit_percent | probit_percent |
|---|---|---|---|---|
| 12 | 50 | female | 37.2 | 37.6 |
| 16 | 50 | female | 55.5 | 54.8 |
Table 25.1 gives similar predictions from the two links. Under logit, the fitted percentage rises from 37.2 to 55.5, a difference of 18.2 percentage points. Under probit, the corresponding percentages are 37.6 and 54.8. These are conditional model predictions, not observed cell means or estimated causal effects. The table gives point predictions only; it does not display prediction uncertainty or guarantee individual outcomes.
25.3 Interpret the comparison
The binary models tell a related story for our particular degree definition: at the specified age and sex, predicted college attainment is higher at 16 than at 12 years of parent schooling. Similar logit and probit predictions are reassuring about this particular comparison, not proof that either model is correct. Complete-case selection, ordinary unweighted inference, and unmeasured circumstances limit every result. We can describe conditional associations and fitted probabilities here; we cannot claim a representative national effect or a causal return to parental schooling.
25.4 Review
- Define the outcome and common sample before comparing models; a missing degree is not non-attainment.
glm()fits logit or probit binary models; their coefficient scales differ.tidy()reports model coefficients, andnobs()checks contributing observations.predict()returns conditional probabilities withtype = "response"; state the predictor values and check their support.- Compare probability predictions with their sample and model limitations, without a causal interpretation or survey-design correction.
The optional Programming part, Chapter 27 and Chapter 28, provides additional detail about R expressions and objects when you need it.
25.5 Exercises
- Reconstruct the shared sample, logit model, and probit model. Check support and reproduce the predicted probabilities for female respondents aged 50 at parent schooling 12 and 16. Interpret the probability differences and explain why the raw logit and probit coefficients are not directly comparable.
25.6 Worked answers
NoteCompare probability predictions
Restore the saved data independently:
library(tidyverse)
library(tidystata)
answer_gss24 <- readRDS("gss24.rds")answer_model_data <- answer_gss24 |>
select(id, educ, maeduc, age, sex, degree) |>
filter(!is.na(id), !is.na(educ), !is.na(maeduc), !is.na(age),
!is.na(sex), !is.na(degree)) |>
mutate(
sex = factor(sex, levels = c("male", "female")),
college = if_else(degree %in% c("bachelor's", "graduate"), 1, 0)
)
answer_logit <- stats::glm(
college ~ maeduc + age + sex, data = answer_model_data,
family = stats::binomial(link = "logit")
)
answer_probit <- stats::glm(
college ~ maeduc + age + sex, data = answer_model_data,
family = stats::binomial(link = "probit")
)
answer_model_data |>
filter(age == 50, sex == "female", maeduc %in% c(12, 16)) |>
count(maeduc, age, sex)# A tibble: 2 × 4
maeduc age sex n
<dbl> <dbl> <fct> <int>
1 12 50 female 11
2 16 50 female 2
answer_cases <- tibble(
maeduc = c(12, 16), age = c(50, 50),
sex = factor(c("female", "female"), levels = c("male", "female"))
)
answer_predictions <- answer_cases |>
mutate(
logit_percent = 100 * stats::predict(
answer_logit, newdata = answer_cases, type = "response"
),
probit_percent = 100 * stats::predict(
answer_probit, newdata = answer_cases, type = "response"
)
)
answer_predictions# A tibble: 2 × 5
maeduc age sex logit_percent probit_percent
<dbl> <dbl> <fct> <dbl> <dbl>
1 12 50 female 37.2 37.6
2 16 50 female 55.5 54.8
Logit predicts 37.2% and 55.5% college attainment; probit predicts 37.6% and 54.8%. Each comparison holds modeled age at 50 and sex at female, with four years’ difference in parent schooling. The logit difference is 18.2 percentage points; the probit difference is 17.2. The links use different coefficient scales but give similar probability contrasts here. Both combinations are observed, although sparsely; the predictions pool the model’s entire sample rather than averaging those cells. They remain model-dependent point predictions, with no causal interpretation or survey-design correction.