library(haven) # import Stata data and labels
library(dplyr) # data manipulation
library(sjmisc) # frq() frequency tables
library(ggplot2) # graphs
library(nnet) # multinomial logit: multinom()
library(cleanplots) # graph scheme
cleanplots_defaults()
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("polviews3", "polviews", "age", "woman", "race3", "college", "married",
"parent", "faminc", "income", "employed", "conserv", "year")]
gss <- gss[complete.cases(gss[c("conserv", "age", "woman", "race3", "college",
"married", "parent", "faminc")]), ]
# Treat binary/nominal variables as factors
gss <- gss |>
mutate(across(c(polviews3, woman, race3, college, married, parent, employed), as_factor)) |>
droplevels()
# 3 category political views var
pol3mod <- multinom(polviews3 ~ age + woman + race3 + college + married + parent +
faminc + employed, data = gss, trace = FALSE)
# Wald tests for combining each pair of outcomes (like mlogtest, combine):
# all coefficients except the intercepts are the same for the two outcomes
combine_tests <- function(model) {
b <- coef(model) # one row per outcome (other than the base outcome)
V <- vcov(model)
bvec <- setNames(as.vector(t(b)), colnames(V))
outcomes <- model$lev
terms <- colnames(b)[-1]
pairs <- combn(outcomes, 2)
out <- lapply(seq_len(ncol(pairs)), function(i) {
o1 <- pairs[1, i]
o2 <- pairs[2, i]
# differences between the two outcomes' coefficients (base outcome = 0)
L <- matrix(0, length(terms), length(bvec), dimnames = list(terms, names(bvec)))
for (t in terms) {
if (o1 != outcomes[1]) L[t, paste0(o1, ":", t)] <- 1
if (o2 != outcomes[1]) L[t, paste0(o2, ":", t)] <- -1
}
d <- L %*% bvec
chi2 <- as.numeric(t(d) %*% solve(L %*% V %*% t(L)) %*% d)
data.frame(outcomes = paste(o1, "&", o2), chi2 = round(chi2, 3), df = length(terms),
p = round(pchisq(chi2, length(terms), lower.tail = FALSE), 3))
})
do.call(rbind, out)
}
combine_tests(pol3mod)
![[4.4.a] Lowess curve of the proportion saying premarital sex is not wrong at all, by GSS survey year, over jittered points at 0 and 1: it rises from about 0.27 in 1972 to about 0.40 by 1980, stays nearly flat around 0.40 to 0.45 until about 2000, then rises to about 0.66 in 2021.](fig/cda-04_04_a-prems-01.png)
![[4.4.a] Smoothed curve of the proportion saying premarital sex is not wrong at all, by GSS survey year, over jittered points at 0 and 1: it rises from about 0.27 in 1972 to about 0.40 by 1980, stays nearly flat around 0.40 to 0.45 until about 2000, then rises to about 0.66 in 2021.](appendix-4a-combining-outcomes_files/figure-html/unnamed-chunk-4-1.png)