library(haven) # import Stata data and labels
library(dplyr) # data manipulation
library(ggplot2) # graphs
library(nnet) # multinomial logit: multinom()
library(marginaleffects) # predictions and marginal effects
library(cleanplots) # graph scheme
cleanplots_defaults()
# Center changes in continuous IVs on each observation's value
options(marginaleffects_contrast_direction = "center")
gss <- read_dta("https://tdmize.github.io/data/data/cda_gss.dta")
gss <- gss |> filter(year >= 2016 & year <= 2021)
# Keep model variables; drop any missing data
gss <- gss[c("conserv", "polviews3", "age", "woman", "race3", "college",
"married", "parent", "faminc", "employed")]
gss <- gss[complete.cases(gss[c("conserv", "age", "woman", "race3", "college",
"married", "parent", "faminc")]), ]
# Not conservative = 1; treat binary/nominal variables as factors
gss <- gss |>
mutate(notconserv = 1 - conserv,
across(c(polviews3, woman, race3, college, married, parent, employed), as_factor)) |>
droplevels()
# 4.0.a Conservative = 1
conmod <- glm(conserv ~ age + woman + race3 + college + married + parent +
faminc + employed, data = gss, family = binomial(link = "logit"))
# 4.0.b Not Conservative = 1
notconmod <- glm(notconserv ~ age + woman + race3 + college + married + parent +
faminc + employed, data = gss, family = binomial(link = "logit"))
# Predictions at ages 18 to 88 from each model
pr_age <- bind_rows(
"[4.0.a] Pr(Identifying as Conservative)" = avg_predictions(conmod, by = "age",
newdata = datagrid(age = seq(18, 88, 2), grid_type = "counterfactual")),
"[4.0.b] Pr(NOT Identifying as Conservative)" = avg_predictions(notconmod, by = "age",
newdata = datagrid(age = seq(18, 88, 2), grid_type = "counterfactual")),
.id = "model")
pr_labels <- data.frame(model = rep(unique(pr_age$model), each = 2),
age = c(30, 70, 30, 70), y = c(.2, .3, .8, .7),
label = c("Pr = 0.27", "Pr = 0.38", "Pr = 0.73", "Pr = 0.62"))
ggplot(pr_age, aes(x = age, y = estimate, color = model, linetype = model)) +
geom_line() +
geom_text(data = pr_labels, aes(y = y, label = label), show.legend = FALSE) +
facet_wrap(~ model) +
scale_color_manual(values = c("red", "purple")) +
scale_linetype_manual(values = c("solid", "dashed")) +
scale_x_continuous(breaks = seq(18, 88, 10)) +
scale_y_continuous(breaks = seq(0, 1, .2), limits = c(0, 1)) +
labs(x = "age of respondent", y = "Pr(y = 1)") +
theme(legend.position = "none")
![Two panels of predicted probabilities by age from 18 to 88: [4.0.a] Pr(conservative) rises from about 0.24 to 0.43 (0.27 at age 30, 0.38 at age 70); [4.0.b] Pr(not conservative) falls from about 0.76 to 0.57 (0.73 at 30, 0.62 at 70).](fig/cda-04_01_a-bin_compare-01.png)
![Two panels of predicted probabilities by age from 18 to 88: [4.0.a] Pr(conservative) rises from about 0.24 to 0.43 (0.27 at age 30, 0.38 at age 70); [4.0.b] Pr(not conservative) falls from about 0.76 to 0.57 (0.73 at 30, 0.62 at 70).](4-0-nominal-dvs_files/figure-html/unnamed-chunk-1-1.png)
![[4.0.c] Both lines on one graph by age: Pr(conservative) rises from about 0.24 to 0.43 and Pr(not conservative) falls from about 0.76 to 0.57; the two lines mirror each other around 0.5.](fig/cda-04_01_b-bin_compare-02.png)
![[4.0.c] Both lines on one graph by age: Pr(conservative) rises from about 0.24 to 0.43 and Pr(not conservative) falls from about 0.76 to 0.57; the two lines mirror each other around 0.5.](4-0-nominal-dvs_files/figure-html/unnamed-chunk-2-1.png)
![[4.2.a] Predicted probabilities of each political ideology by age with confidence intervals: Pr(liberal) falls from about 0.35 to 0.27, Pr(moderate) falls from about 0.41 to 0.29, and Pr(conservative) rises from about 0.24 to 0.44; the conservative line crosses the others in the late 40s to 50s.](fig/cda-04_02_a-polviews.png)