We’re going to look at models with multiple categorical predictors, including their interactions. I’m going to use dummy coding and model comparison rather than contrast coding.
We are going to use the full GSS data for this example. That’s the cumulative file you saved in the setup chapter. We need it here because we want to look at about fifty years of surveys.
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. This is the same trade-off we saw in Chapter 13. Here we can afford the extra parameters because we have tens of thousands of cases.
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 referenceEach one is a difference from a reference category, just like in Chapter 13. There are just a lot of them! With this many coefficients, the table isn’t a very useful way to understand the model.
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 every year as equally important, even though the GSS didn’t interview the same number of people every year. And averaging over half men and half women is also something we made up. Balanced predictions describe a hypothetical population that you’ve specified. That’s often what you want for comparisons (because it holds the composition constant), but you should say that’s what you did when you report the results.
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 (AIC and BIC “weights”) turn the differences into something like relative probabilities. To get them, you take the difference between each model and the best one, multiply by \(-\frac{1}{2}\), exponentiate, and then divide by the total. The BIC weights are the same idea as the posterior probabilities in Chapter 10.
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 goes up with every interaction we add (as always), but the three-way model with 290 parameters is the worst according to both AIC and BIC.
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 isn’t a problem. As we saw in Chapter 10, AIC and BIC are trying to answer different questions, and BIC’s penalty gets bigger with the sample size. Here \(\log(n)\) is about 10.7, so BIC’s penalty per parameter is more than three times as big as AIC’s. An interaction can help with prediction without being something we’d want to claim is “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 (basically nothing). So the gap between men and women gets smaller as education goes up. The additive model can’t show that, since it assumes the sex gap is the same at every level of education.
It is tempting to write that “education reduces the sex gap in television watching.” But that’s not what we found. We compared people with different amounts of education in a survey. We didn’t change anyone’s education, and lots of other things (birth cohort, employment, household composition, etc.) are different across these groups too.
A better way to say it: among respondents with more education, the predicted difference between women and men is smaller.