Models for categorical data
Packages
Loading packages:
library(lme4)
library(MASS)
library(stargazer)
library(dplyr)
library(ggplot2)
library(ggeffects)
library(nnet)Data
Loading the data:
load('/Users/davidhamad/Documents/Cand.Scient.Pol/2. Semester/Statistical models beyond linear regression - applied statistics for political scientists/9-10) Categorical outcomes/kap6.rda')Rename
Renaming for convenience:
df <- kap6 %>%
mutate(Education = Uddannelse,
Income = Indtaegt,
Female = Kvinde,
Burden = Belastning,
Immigrant = Innvandrer)Exercise 1
Q1: Begin by choosing a reference category. Which one did you choose? Why? You can use the functions as.factor() and relevel().
A1: I chose Socialdemokraterne as the reference category because it is a large category with a familiar baseline.
table(df$Parti)
Andet Dansk Folkeparti Det Konservative Folkeparti
17 143 65
Det Radikale Venstre Enhedslisten Kristendemokraterne
134 77 8
Liberal Alliance Socialdemokraterne Socialistisk Folkeparti
48 268 108
Venstre
311 df <- df %>%
mutate(Parti = as.factor(Parti),
Parti = relevel(Parti, ref = "Socialdemokraterne"))table(df$Parti)
Socialdemokraterne Andet Dansk Folkeparti
268 17 143
Det Konservative Folkeparti Det Radikale Venstre Enhedslisten
65 134 77
Kristendemokraterne Liberal Alliance Socialistisk Folkeparti
8 48 108
Venstre
311 Q2: Estimate a multinomial model where choice of party is a function of education and income.
A2: The coefficients on income support the redistribution hypothesis. The bourgeois parties (V, KF, LA) have positive coefficients, meaning that higher income increases the probability of voting for these parties rather than S. Education is also interesting. It seems that education is especially important when voting for RV (and also SF).
mod1.nom <- multinom(Parti ~ Education + Income, data = df)# weights: 40 (27 variable)
initial value 2516.725507
iter 10 value 2297.906803
iter 20 value 2182.108328
iter 30 value 2111.031390
final value 2110.013593
convergedsummary(mod1.nom)Call:
multinom(formula = Parti ~ Education + Income, data = df)
Coefficients:
(Intercept) Education Income
Andet -3.47649926 0.02586823 0.05825070
Dansk Folkeparti -0.02936041 -0.03946825 -0.02262418
Det Konservative Folkeparti -2.49096556 0.02823281 0.13025930
Det Radikale Venstre -3.14790162 0.11733117 0.12734286
Enhedslisten -1.51281720 0.02578023 -0.01401907
Kristendemokraterne -3.00310057 -0.01311761 -0.05207996
Liberal Alliance -3.09991119 0.03387547 0.14070260
Socialistisk Folkeparti -1.82459370 0.07611650 -0.01328518
Venstre -1.06960446 0.02104517 0.15310841
Std. Errors:
(Intercept) Education Income
Andet 0.8716862 0.05979221 0.09979420
Dansk Folkeparti 0.3071891 0.02289756 0.04053397
Det Konservative Folkeparti 0.4756282 0.03150253 0.05332830
Det Radikale Venstre 0.4264347 0.02693849 0.04309100
Enhedslisten 0.4198162 0.02998338 0.05049790
Kristendemokraterne 1.0260423 0.07701172 0.13630737
Liberal Alliance 0.5904935 0.03844719 0.06496080
Socialistisk Folkeparti 0.3920117 0.02703611 0.04384018
Venstre 0.2869829 0.01951165 0.03339975
Residual Deviance: 4220.027
AIC: 4274.027 Q3: The model effectively estimates the effect of income at equal levels of education. Which party has the voters with the highest/lowest revenue? In your opinion, does our knowledge of respondents’ income help us understand how voters discriminate between parties?
A3: Looking at the income coefficients (all relative to S), four parties have voters with significantly higher incomes than S voters: V (β = 0.15), LA (β = 0.14) KF (β = 0.13), and RV (β = 0.13). All four are statistically significant. The parties with the lowest-income voters are K, DF, Ø, and SF, all with small negative coefficients, But none of these are statistically significant, meaning we cannot distinguish their voters’ income from S voters in this data. Income therefore helps us discriminate the bourgeois bloc (V, KF, LA) from S, and it interestingly also separates RV from S.
coefs <- summary(mod1.nom)$coefficients
ses <- summary(mod1.nom)$standard.errors
z <- coefs / ses
p <- (1 - pnorm(abs(z))) * 2
data.frame(
coef = round(coefs[, "Income"], 3),
se = round(ses[, "Income"], 3),
z = round(z[, "Income"], 2),
p = round(p[, "Income"], 4),
sig = ifelse(p[, "Income"] < 0.001, "***",
ifelse(p[, "Income"] < 0.01, "**",
ifelse(p[, "Income"] < 0.05, "*",
ifelse(p[, "Income"] < 0.1, ".", ""))))
) |>
(\(x) x[order(x$coef), ])()Q4: Interpret the marginal effects (by exponentiating the slope coefficients). Which parties are similar to your reference category with respect to the effect of income? For the party-pairs that are significantly different, what is the effect of income on the respondents’ choice of party?
A4: Exponentiating the income coefficients converts the log-odds into odds ratios, which show how a one-unit increase in income changes the odds of voting for a given party rather than for Socialdemokraterne. Five parties are not significantly different from S with respect to the effect of income: Kristendemokraterne, Dansk Folkeparti, Enhedslisten, Socialistisk Folkeparti, and Andet. Substantively, this means that the voters of these parties have income profiles similar to S voters, once education is controlled for. Three of these (V, LA, KF) are classic bourgeois parties, and the income effect fits the redistribution hypothesis. Higher-income voters prefer parties that favor less redistribution. The fourth, RV, is more surprising, but it attracts high-income voters. This is consistent with its appeal running through the educational-cultural dimension rather than the economic one.
or <- exp(coefs[, "Income"])data.frame(
coef = round(coefs[, "Income"], 3),
odds_ratio = round(or, 3),
pct_change = round((or - 1) * 100, 1),
p = ifelse(p[, "Income"] < 0.0001, "<0.0001",
format(round(p[, "Income"], 4), nsmall = 4)),
sig = ifelse(p[, "Income"] < 0.001, "***",
ifelse(p[, "Income"] < 0.01, "**",
ifelse(p[, "Income"] < 0.05, "*",
ifelse(p[, "Income"] < 0.1, ".", ""))))
) |>
(\(x) x[order(x$coef), ])()Q5: Make predictions for a scenario where a respondent is rich (Income = 9) by filling in the equation for Dansk Folkeparti (DF), Venstre (V) and Socialistisk Folkeparti (SF). What do you find?
A5: For a rich respondent (Income = 9) with average education, the predicted probabilities for the three parties of interest are: Venstre 33.9%, Dansk Folkeparti 8.9%, and Socialistisk Folkeparti 7.3%. Venstre is clearly the most likely choice. V voters have the highest-income profile of any party. DF and SF are both well below Venstre, and interestingly they are fairly close to each other, despite being on opposite ends of the cultural spectrum.
scenario_rich <- data.frame(
Education = mean(df$Education, na.rm = TRUE),
Income = 9
)preds_rich <- predict(mod1.nom, newdata = scenario_rich, type = "probs")
round(preds_rich, 3) Socialdemokraterne Andet Dansk Folkeparti
0.189 0.014 0.089
Det Konservative Folkeparti Det Radikale Venstre Enhedslisten
0.073 0.119 0.051
Kristendemokraterne Liberal Alliance Socialistisk Folkeparti
0.005 0.047 0.073
Venstre
0.339 preds_rich[c("Dansk Folkeparti", "Venstre", "Socialistisk Folkeparti")] Dansk Folkeparti Venstre Socialistisk Folkeparti
0.08938308 0.33904848 0.07327798 Q6 + Q7: Reestimate your model with a new reference category? What happened? Based on your new model fit, make predictions for the same scenarios. What did you find? Why?
A6 + A7: I reestimated the model with Venstre as the reference category. The coefficients change. However, when I compute predictions for the same rich-voter scenario under both models and subtract them, all differences are zero. The reference category only affects how coefficients are presented, not what the model estimates.
df_v <- df %>%
mutate(Parti = relevel(Parti, ref = "Venstre"))mod2.nom <- multinom(Parti ~ Education + Income, data = df_v, trace = FALSE)summary(mod2.nom)$coefficients (Intercept) Education Income
Socialdemokraterne 1.0696836 -0.021043714 -0.15311519
Andet -2.4070773 0.004844270 -0.09487471
Dansk Folkeparti 1.0402804 -0.060511852 -0.17573784
Det Konservative Folkeparti -1.4214332 0.007194775 -0.02285379
Det Radikale Venstre -2.0781289 0.096283116 -0.02578075
Enhedslisten -0.4430975 0.004728908 -0.16713024
Kristendemokraterne -1.9334603 -0.034279840 -0.20497515
Liberal Alliance -2.0300525 0.012823998 -0.01243127
Socialistisk Folkeparti -0.7549107 0.055072451 -0.16640355preds_rich_s <- predict(mod1.nom, newdata = scenario_rich, type = "probs")preds_rich_v <- predict(mod2.nom, newdata = scenario_rich, type = "probs")round(preds_rich_s - preds_rich_v[names(preds_rich_s)], 4) Socialdemokraterne Andet Dansk Folkeparti
0 0 0
Det Konservative Folkeparti Det Radikale Venstre Enhedslisten
0 0 0
Kristendemokraterne Liberal Alliance Socialistisk Folkeparti
0 0 0
Venstre
0 Exercise 2
Q1: Calculate the regression coefficients in an intercept-only model by hand.
A1: The log-odds show the log-ratio of each party’s vote share relative to S. Venstre is the only party with more votes than S, so it has a positive value (+0.15). All other parties have fewer votes than S and therefore have negative values. S itself has a value of 0 because it is the reference category.
tab <- df %>%
filter(!is.na(Parti)) %>%
group_by(Parti) %>%
reframe(n = n()) %>%
mutate(
N = sum(n),
p = n / N,
p_ref = p[Parti == "Socialdemokraterne"],
odds_vs_S = p / p_ref,
logodds = log(odds_vs_S)
) %>%
arrange(logodds)tabQ2: Calculate the parameters of a multinomial model of party choice with a binary predictor by hand. The data contains a predictor that flags the respondent’s gender (Kvinde / Female).
A2: Among men, V is the only party with positive odds against S, so V is the most popular party for men relative to S. Women are more likely to vote SF and Ø, relative to S.
df0 <- df %>% filter(!is.na(Parti), Female == 0)
df1 <- df %>% filter(!is.na(Parti), Female == 1) calc_logodds <- function(d, ref = "Socialdemokraterne") {
d %>%
group_by(Parti) %>%
reframe(n = n()) %>%
mutate(
p = n / sum(n),
odds_vs_ref = p / p[Parti == ref],
logodds = log(odds_vs_ref)
)
}men_tab <- calc_logodds(df0)
women_tab <- calc_logodds(df1)combined <- men_tab %>%
dplyr::select(Parti, logodds_men = logodds, odds_men = odds_vs_ref) %>%
left_join(
women_tab %>% dplyr::select(Parti, logodds_women = logodds, odds_women = odds_vs_ref),
by = "Parti"
) %>%
mutate(
intercept_a = logodds_men,
slope_b = log(odds_women / odds_men)
) %>%
filter(Parti != "Socialdemokraterne")combined %>% dplyr::select(Parti, intercept_a, slope_b) %>% arrange(intercept_a)Q3: How would you explain the slope coefficient to your neighbour/classmate?
A3: The slope tells you how the odds of voting for a party (rather than for Socialdemokraterne) change when we move from a man to a woman. A negative slope means women are less likely than men to vote for that party rather than S. A positive slope means women are more likely. To get a more intuitive number, we exponentiate the slope: exp(b) gives the odds-ratio. For example, Liberal Alliance has a slope of −1.61, so exp(−1.61) = 0.20, meaning women have only a fifth of the odds that men have of voting LA rather than S.
Exercise 3
Q1: Recode the variable in R into an ordinal variable.
A1: See the recoding below.
df <- df %>%
mutate(Burden_rev = 10 - Burden,
Burden_ord = ordered(Burden_rev, levels = 0:10))class(df$Burden_ord)[1] "ordered" "factor" levels(df$Burden_ord) [1] "0" "1" "2" "3" "4" "5" "6" "7" "8" "9" "10"table(df$Burden_ord)
0 1 2 3 4 5 6 7 8 9 10
23 36 161 193 148 389 124 161 116 52 72 Q2: Fit an ordinal model where agreement to the statement is a function of the respondent’s income, education and immigrant background.
A2: The model shows three negative and statistically significant effects: higher income, higher education and immigrant background all decrease the probability of being in a higher Burden category; that is, they reduce the likelihood of viewing immigration as a burden. The 10 cutpoints separate the 11 ordered categories on the latent scale.
mod.ord <- polr(Burden_ord ~ Income + Education + Immigrant, Hess = T, df)summary(mod.ord)Call:
polr(formula = Burden_ord ~ Income + Education + Immigrant, data = df,
Hess = T)
Coefficients:
Value Std. Error t value
Income -0.06457 0.01773 -3.643
Education -0.08159 0.01060 -7.695
Immigrant -0.53010 0.15567 -3.405
Intercepts:
Value Std. Error t value
0|1 -5.7690 0.2768 -20.8454
1|2 -4.7739 0.2167 -22.0302
2|3 -3.3085 0.1803 -18.3485
3|4 -2.4726 0.1699 -14.5523
4|5 -1.9814 0.1650 -12.0057
5|6 -0.8658 0.1570 -5.5146
6|7 -0.4622 0.1563 -2.9571
7|8 0.2229 0.1591 1.4011
8|9 0.9805 0.1706 5.7486
9|10 1.5980 0.1895 8.4351
Residual Deviance: 5526.293
AIC: 5552.293
(203 observations deleted due to missingness)Q3: Present the results in stargazer.
A3: See the results below.
stargazer(mod.ord,
ord.intercepts = TRUE,
type = "text")
========================================
Dependent variable:
---------------------------
Burden_ord
----------------------------------------
Income -0.065***
(0.018)
Education -0.082***
(0.011)
Immigrant -0.530***
(0.156)
0| 1 -5.769***
(0.277)
1| 2 -4.774***
(0.217)
2| 3 -3.308***
(0.180)
3| 4 -2.473***
(0.170)
4| 5 -1.981***
(0.165)
5| 6 -0.866***
(0.157)
6| 7 -0.462***
(0.156)
7| 8 0.223
(0.159)
8| 9 0.980***
(0.171)
9| 10 1.598***
(0.189)
----------------------------------------
Observations 1,299
========================================
Note: *p<0.1; **p<0.05; ***p<0.01Q4: What is the marginal effect of income? Can you make a sentence that states the finding?
A4: A one-decile increase in income lowers the odds of being in a higher Burden category by approximately 6.3% (OR = 0.937), holding education and immigrant background constant. Substantively, wealthier respondents are less likely to see immigration as a burden.
exp(coef(mod.ord)["Income"]) Income
0.9374682 Q5: Calculate the probability that an observation is below the first – then the second – cutpoint in the model.
A5: For an average respondent (mean income, mean education, non-immigrant), the predicted probability of being in the lowest category (Burden_rev = 0, i.e. not viewing immigration as a burden at all) is approximately X%. The probability of being in either of the two lowest categories is approximately Y%.
b_inc <- coef(mod.ord)["Income"]
b_edu <- coef(mod.ord)["Education"]
b_imm <- coef(mod.ord)["Immigrant"]tau1 <- mod.ord$zeta["0|1"]
tau2 <- mod.ord$zeta["1|2"]xb <- b_inc * mean(df$Income, na.rm = TRUE) +
b_edu * mean(df$Education, na.rm = TRUE) +
b_imm * 0xb Income
-1.441242 logodds_tau1 <- tau1 - xb
logodds_tau2 <- tau2 - xbp_below_tau1 <- plogis(logodds_tau1)
p_below_tau2 <- plogis(logodds_tau2)p_below_tau1 0|1
0.0130248 p_below_tau2 1|2
0.03446845 Q6: Using the same scenario, calculate the probability that an observation falls between cutpoints “0|1” and “1|2” by subtracting one cut point from the other.
A6: Subtracting the probability of being below cutpoint 0|1 from the probability of being below cutpoint 1|2 isolates the probability of being in category 1 specifically. For an average respondent, this gives 3.6%. So a typical Dane has approximately a 3.6% probability of placing themselves in category 1 on the Burden scale.
p_kat1 <- p_below_tau2 - p_below_tau1p_kat1 1|2
0.02144365 Q7: Check your results by having R calculate the predicted probabilities. Remember to keep your scenario constant.
A7: See the predicted probabilities below.
scenario <- data.frame(
Income = mean(df$Income, na.rm = TRUE),
Education = mean(df$Education, na.rm = TRUE),
Immigrant = 0
)predict(mod.ord, newdata = scenario, type = "probs") 0 1 2 3 4 5 6 7 8
0.01302480 0.02144365 0.09939390 0.12895222 0.10534712 0.27186722 0.08688645 0.11387419 0.07767918
9 10
0.03584824 0.04568303 Q8: Calculate the cumulative sum of the probabilities.
A8: Calling cumsum() on the predicted probabilities gives the cumulative probability. So the first cumsum value (0.046) matches my hand-calculated. The difference is just perspective. Predict() shows the probability of each specific category, while cumsum() shows the probability of being in that category or any lower one.
predict(mod.ord, newdata = scenario, type = "probs") %>% cumsum() 0 1 2 3 4 5 6 7 8
0.01302480 0.03446845 0.13386235 0.26281457 0.36816169 0.64002891 0.72691536 0.84078955 0.91846873
9 10
0.95431697 1.00000000 Exercise 4
Q1: Create 10 binary variables (one for each of the cut points in the “Burden” variable), then run 10 binomial regressions.
A1: These 10 binary logits test the parallel regressions assumption. Income and Education are stable across cutpoints, so the assumption holds for them. Immigrant varies more and gets stronger at higher cutpoints, so the assumption is questionable for that variable.
binary_models <- list()for (k in 0:9) {
df$y_temp <- as.numeric(df$Burden_rev < (k + 1))
binary_models[[paste0("<", k+1)]] <- glm(
y_temp ~ Income + Education + Immigrant,
data = df,
family = "binomial"
)
}coef_table <- sapply(binary_models, function(m) coef(m)[c("Income", "Education", "Immigrant")])round(coef_table, 3) <1 <2 <3 <4 <5 <6 <7 <8 <9 <10
Income 0.056 0.026 0.052 0.082 0.070 0.062 0.068 0.062 0.048 0.087
Education 0.038 0.068 0.100 0.088 0.098 0.082 0.077 0.063 0.064 0.065
Immigrant 0.895 0.651 0.735 0.642 0.623 0.418 0.439 0.278 0.177 0.697Q2: Display the 10 binary models and the ordinal model from the previous exercise in a stargazer() results table. Are the regression coefficients very different?
A2: Income and Education coefficients in the binary models are close to those in the ordinal model, supporting parallel regressions. Immigrant varies much more (from 0.18 to 0.90) and the ordinal model’s value of 0.53 is essentially an average of these which is masking that the effect is much weaker at lower cutpoints and much stronger at upper cutpoints. The assumption is therefore questionable for immigrant.
stargazer(binary_models[1:5],
type = "text",
column.labels = paste0("<", 1:5),
model.numbers = FALSE)
===================================================================
Dependent variable:
-------------------------------------------------
y_temp
<1 <2 <3 <4 <5
-------------------------------------------------------------------
Income 0.056 0.026 0.052* 0.082*** 0.070***
(0.082) (0.052) (0.029) (0.023) (0.022)
Education 0.038 0.068** 0.100*** 0.088*** 0.098***
(0.049) (0.032) (0.019) (0.015) (0.014)
Immigrant 0.895* 0.651* 0.735*** 0.642*** 0.623***
(0.526) (0.353) (0.211) (0.185) (0.182)
Constant -5.135*** -4.345*** -3.542*** -2.699*** -2.261***
(0.787) (0.507) (0.301) (0.232) (0.211)
-------------------------------------------------------------------
Observations 1,299 1,299 1,299 1,299 1,299
Log Likelihood -105.448 -222.875 -534.369 -741.870 -820.797
Akaike Inf. Crit. 218.896 453.751 1,076.738 1,491.741 1,649.593
===================================================================
Note: *p<0.1; **p<0.05; ***p<0.01stargazer(binary_models[6:10],
type = "text",
column.labels = paste0("<", 6:10),
model.numbers = FALSE)
=================================================================
Dependent variable:
-----------------------------------------------
y_temp
<6 <7 <8 <9 <10
-----------------------------------------------------------------
Income 0.062*** 0.068*** 0.062** 0.048 0.087*
(0.022) (0.023) (0.028) (0.037) (0.048)
Education 0.082*** 0.077*** 0.063*** 0.064*** 0.065**
(0.013) (0.014) (0.016) (0.021) (0.027)
Immigrant 0.418** 0.439** 0.278 0.177 0.697
(0.197) (0.217) (0.256) (0.334) (0.528)
Constant -0.858*** -0.424** 0.473** 1.302*** 1.671***
(0.189) (0.195) (0.223) (0.286) (0.361)
-----------------------------------------------------------------
Observations 1,299 1,299 1,299 1,299 1,299
Log Likelihood -806.452 -727.047 -556.568 -364.433 -238.816
Akaike Inf. Crit. 1,620.903 1,462.093 1,121.137 736.867 485.632
=================================================================
Note: *p<0.1; **p<0.05; ***p<0.01stargazer(mod.ord,
type = "text",
ord.intercepts = TRUE)
========================================
Dependent variable:
---------------------------
Burden_ord
----------------------------------------
Income -0.065***
(0.018)
Education -0.082***
(0.011)
Immigrant -0.530***
(0.156)
0| 1 -5.769***
(0.277)
1| 2 -4.774***
(0.217)
2| 3 -3.308***
(0.180)
3| 4 -2.473***
(0.170)
4| 5 -1.981***
(0.165)
5| 6 -0.866***
(0.157)
6| 7 -0.462***
(0.156)
7| 8 0.223
(0.159)
8| 9 0.980***
(0.171)
9| 10 1.598***
(0.189)
----------------------------------------
Observations 1,299
========================================
Note: *p<0.1; **p<0.05; ***p<0.01Q3: Create a coefplot for each of the 10 coefficients on income. To do so, you will want to make a small data frame with one observation per model and variables reporting the slope coefficient, its standard error and the cutvalue associated with each model. What do you see?
A3: The coefficient on income is stable across the 10 cutpoints. Confidence intervals overlap substantially, especially in the middle of the scale where most data are concentrated. The wider intervals at cutpoints 1, 9 and 10 reflect the smaller number of observations in the tails of the Burden distribution, not a real change in the underlying effect. The plot supports the parallel regressions assumption for income: the effect is essentially constant across the Burden scale, justifying its use in the ordered logit model.
income_df <- data.frame(
cutpoint = 1:10,
estimate = sapply(binary_models, function(m) coef(m)["Income"]),
se = sapply(binary_models, function(m) summary(m)$coefficients["Income", "Std. Error"])
)income_df$lower <- income_df$estimate - 1.96 * income_df$se
income_df$upper <- income_df$estimate + 1.96 * income_df$seggplot(income_df, aes(x = cutpoint, y = estimate)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
geom_point(size = 3) +
scale_x_continuous(breaks = 1:10) +
labs(x = "Cutpoint",
y = "Income coefficient",
title = "Effect of income across cutpoints") +
theme_minimal()
Q4: Create a similar coefplot for the intercepts/cutpoints of the model. There is no need to backtransform them.
A4: The intercepts increase monotonically. There are no reversals, confirming that the categories of Burden are genuinely ordered. Together, this supports the second assumption of the ordinal model: the categories are both ordered and separable.
intercept_df <- data.frame(
cutpoint = 1:10,
estimate = sapply(binary_models, function(m) coef(m)["(Intercept)"]),
se = sapply(binary_models, function(m) summary(m)$coefficients["(Intercept)", "Std. Error"])
)intercept_df$lower <- intercept_df$estimate - 1.96 * intercept_df$se
intercept_df$upper <- intercept_df$estimate + 1.96 * intercept_df$seggplot(intercept_df, aes(x = cutpoint, y = estimate)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
geom_point(size = 3) +
scale_x_continuous(breaks = 1:10) +
labs(x = "Cutpoint",
y = "Intercept (log-odds)",
title = "Intercepts across cutpoints") +
theme_minimal()