We’re going to look at models with multiple categorical predictors, including their interactions. The approach here is dummy coding and model comparison rather than contrast coding.
We are going to use the full GSS data for this example. That is the cumulative file you built in the setup chapter, because we want variation across half a century and no single survey can provide it.
d <- readRDS(here::here("data", "gss-1972-2024.rds")) |>
haven::zap_labels() |>
select(tvhours, degree, year, sex) |>
drop_na() |>
mutate(female = if_else(sex == 2, 1, 0), # female
degree_fac = factor(degree), # degree as factor
year_fac = factor(year)) # svy year as factor
c(respondents = nrow(d), survey_years = nlevels(d$year_fac)) respondents survey_years
46000 29 We’re going to consider a few categorical predictors of tvhours:
female: whether the respondent is female (0, 1)degree_fac: highest degree earned from none (0) to graduate degree (4); stored as a factoryear_fac: a factor encoding the survey year1 from 1975 to 2024 (with gaps)1 year is clearly a “numeric” variable in the obvious sense. We have year_fac stored as a factor so that we make no assumptions about its functional (e.g., straight line) relationship to the outcome when we use it in a model. That is exactly the trade-off from the ladder in Chapter 13, decided here in favor of flexibility because we have tens of thousands of observations to spend.
Let’s consider a simple model that uses all of these additively.
m1 <- lm(tvhours ~ degree_fac + female + year_fac,
data = d)The output here is super unwieldy but I’ll put it here if you want to see the summary().
summary(m1)
Call:
lm(formula = tvhours ~ degree_fac + female + year_fac, data = d)
Residuals:
Min 1Q Median 3Q Max
-4.7364 -1.3895 -0.4297 0.9042 21.9266
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.54935 0.06894 51.484 < 2e-16 ***
degree_fac1 -0.72826 0.03205 -22.725 < 2e-16 ***
degree_fac2 -1.19014 0.05425 -21.937 < 2e-16 ***
degree_fac3 -1.61124 0.04060 -39.681 < 2e-16 ***
degree_fac4 -1.91856 0.04910 -39.074 < 2e-16 ***
female 0.18315 0.02345 7.810 5.83e-15 ***
year_fac1977 -0.10716 0.09099 -1.178 0.238899
year_fac1978 -0.21179 0.09091 -2.330 0.019833 *
year_fac1980 -0.04476 0.09208 -0.486 0.626872
year_fac1982 0.18080 0.08694 2.079 0.037577 *
year_fac1983 0.03034 0.09000 0.337 0.736058
year_fac1985 0.05173 0.09099 0.569 0.569654
year_fac1986 0.15821 0.09187 1.722 0.085063 .
year_fac1988 0.29295 0.10274 2.851 0.004357 **
year_fac1989 0.14147 0.10232 1.383 0.166757
year_fac1990 0.04205 0.10475 0.401 0.688087
year_fac1991 0.24879 0.10182 2.444 0.014547 *
year_fac1993 0.09034 0.09017 1.002 0.316402
year_fac1994 0.04430 0.08605 0.515 0.606636
year_fac1996 0.19042 0.08619 2.209 0.027160 *
year_fac1998 0.09076 0.08311 1.092 0.274845
year_fac2000 0.18342 0.08746 2.097 0.035989 *
year_fac2002 0.20883 0.10536 1.982 0.047477 *
year_fac2004 0.15770 0.10567 1.492 0.135608
year_fac2006 0.21870 0.08590 2.546 0.010896 *
year_fac2008 0.25951 0.09455 2.745 0.006061 **
year_fac2010 0.30289 0.09275 3.266 0.001092 **
year_fac2012 0.37884 0.09506 3.985 6.75e-05 ***
year_fac2014 0.30112 0.08934 3.371 0.000751 ***
year_fac2016 0.35250 0.08698 4.053 5.07e-05 ***
year_fac2018 0.25950 0.09086 2.856 0.004294 **
year_fac2021 1.00387 0.08194 12.251 < 2e-16 ***
year_fac2022 0.74754 0.08329 8.975 < 2e-16 ***
year_fac2024 0.72786 0.08482 8.582 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.492 on 45966 degrees of freedom
Multiple R-squared: 0.05467, Adjusted R-squared: 0.05399
F-statistic: 80.55 on 33 and 45966 DF, p-value: < 2.2e-16There are three sets of coefficients:
degree_fac ones that show how each level of degree_fac are different from the “none” reference categoryfemale one that shows how female respondents are different from malesyear_fac ones that show how different each survey is from the 1975 referenceEvery one is a difference from a reference category, exactly as in Chapter 13. There are simply a lot of them. Coefficient tables stop being a useful way to read a model at about this size.
This is easier to see in picture form. I will use plot_predictions() from the marginaleffects package. Using the newdata = "balanced" argument means that the predictions are averaged over equal values of the other predictors (i.e, in the first plot: half male, equal representation from all survey years).
by_degree <- as.data.frame(
avg_predictions(m1, by = "degree_fac", newdata = "balanced")
)
plt(estimate ~ degree_fac, data = by_degree,
ymin = conf.low, ymax = conf.high,
type = "pointrange", lwd = 2,
xlab = "degree_fac", ylab = "tvhours")
plot_predictions(m1,
by = "degree_fac",
newdata = "balanced")
And here it is for the other ones. (I messed with the angle of the x-axis labels to avoid a mess on year.)
by_year <- as.data.frame(
avg_predictions(m1, by = "year_fac", newdata = "balanced")
)
by_year$yr <- as.numeric(as.character(by_year$year_fac))
plt(estimate ~ yr, data = by_year,
ymin = conf.low, ymax = conf.high,
type = "pointrange",
xlab = "year_fac", ylab = "tvhours")
plot_predictions(m1,
by = "year_fac",
newdata = "balanced") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
by_sex <- as.data.frame(
avg_predictions(m1, by = "female", newdata = "balanced")
)
plt(estimate ~ factor(female), data = by_sex,
ymin = conf.low, ymax = conf.high,
type = "pointrange", lwd = 2,
xlab = "female", ylab = "tvhours")
plot_predictions(m1,
by = "female",
newdata = "balanced")
Averaging over an equal mix of survey years treats 1975 and 2024 as equally representative, though the GSS did not sample equally across them. Averaging over half men and half women is similarly a construction. Balanced predictions describe a hypothetical population you have specified, which is often what you want for comparison (it holds the composition fixed while one variable moves), but it is a choice you are making, and worth stating when you report it.
The model above says that there are educational differences and year differences and sex differences. But the model also assumes that those differences are constant. That is, m1 assumes that, for example, the educational differences don’t differ by sex or year.
We can relax that assumption in several ways and compare the results. We can allow any pair of those differences to moderate each other (3 options) or all three to affect each other.
m_deg_sex <- update(m1, tvhours ~ degree_fac * female + year_fac)
m_deg_yr <- update(m1, tvhours ~ degree_fac * year_fac + female)
m_sex_yr <- update(m1, tvhours ~ female * year_fac + degree_fac)
m_all <- update(m1, tvhours ~ degree_fac * female * year_fac)Now we can compare their performance.2
2 The AIC_wt and BIC_wt columns convert the differences into something like relative probabilities: exponentiate half the negative difference from the best model, then divide by the total. Chapter 10 introduced them.
mods <- list(additive = m1,
"degree x sex" = m_deg_sex,
"degree x year" = m_deg_yr,
"sex x year" = m_sex_yr,
"all three" = m_all)
cmp <- data.frame(
model = names(mods),
parameters = map_int(mods, \(m) length(coef(m))),
PRE = round(map_dbl(mods, \(m) summary(m)$r.squared), 4),
AIC = round(map_dbl(mods, AIC), 1),
BIC = round(map_dbl(mods, BIC), 1),
row.names = NULL
)
wt <- \(x) round(exp(-0.5 * (x - min(x))) / sum(exp(-0.5 * (x - min(x)))), 3)
cmp$AIC_wt <- wt(cmp$AIC)
cmp$BIC_wt <- wt(cmp$BIC)
cmp[, c("model", "parameters", "PRE", "AIC", "AIC_wt", "BIC", "BIC_wt")] model parameters PRE AIC AIC_wt BIC BIC_wt
1 additive 34 0.0547 214593.7 0 214899.5 0.997
2 degree x sex 38 0.0553 214570.3 1 214911.0 0.003
3 degree x year 146 0.0580 214657.5 0 215941.7 0.000
4 sex x year 62 0.0556 214602.9 0 215153.3 0.000
5 all three 290 0.0616 214766.6 0 217308.9 0.000PRE rises with every added interaction, as always, and the criteria do not follow it: the three-way model, with 290 parameters, is the worst by both.
With 46,000 cases, it’s not too surprising that the BIC prefers the simple additive model. The AIC, however, prefers the model where degree and sex moderate each other. Let’s take a look at that.
This disagreement is not a malfunction. Chapter 10 explained that the two answer different questions, and that BIC’s penalty grows with the sample size. Here \(\log(n)\) is about 10.7, so BIC charges more than three times AIC’s rate per parameter. A modest interaction can easily be worth keeping for prediction and not worth asserting as real.
preds <- avg_predictions(m_deg_sex,
variables = c("degree_fac", "female"),
newdata = "balanced")plt(estimate ~ degree_fac | factor(female),
data = preds,
type = "pointrange",
ymax = estimate + 1.96 * std.error,
ymin = estimate - 1.96 * std.error,
lw = 2)
ggplot(as.data.frame(preds),
aes(x = degree_fac, y = estimate, color = factor(female))) +
geom_pointrange(aes(ymin = estimate - 1.96 * std.error,
ymax = estimate + 1.96 * std.error),
position = position_dodge(width = 0.3), linewidth = 1) +
labs(x = "degree_fac", y = "tvhours", color = "female")
cells <- as.data.frame(preds)
gaps <- cells |>
select(degree_fac, female, estimate) |>
tidyr::pivot_wider(names_from = female, values_from = estimate,
names_prefix = "f") |>
mutate(gap = round(f1 - f0, 3))
gaps# A tibble: 5 x 4
degree_fac f0 f1 gap
<fct> <dbl> <dbl> <dbl>
1 0 3.63 4.05 0.424
2 1 3.04 3.22 0.185
3 2 2.64 2.71 0.076
4 3 2.24 2.26 0.02
5 4 1.91 1.96 0.049Among respondents with no degree, women are predicted to watch about 0.424 hours more television per day than men. Among those with graduate degrees, the gap is about 0.049 hours (essentially nothing). The sex difference narrows as education rises, and an additive model cannot express that; it would have imposed a single sex gap at every education level, and reported the average of gaps that are in fact quite different.
It is tempting to write that “education reduces the sex gap in television watching.” It does not follow. We have compared people at different education levels in a survey, not intervened on anyone’s schooling, and everything from birth cohort to employment to household composition differs across those groups as well.
The defensible statement: among respondents with more education, the predicted difference between women and men is smaller.