Ordinal analysis: Could this pattern be a scaling artifact?
Julia Rohrer
2026-03-13
Ordinal analysis
Plotting the histogram
First, let’s look at the actual distribution of the responses (Figure 3 in the manuscript).
load("prep.RData")
library(marginaleffects)
library(ggplot2)
library(MASS)
# Plot histogram
ggplot(data = combined, aes(x = satis)) +
geom_histogram(bins = 5, fill = "white", color = "black") +
theme_classic() +
xlab("Life satisfaction score") +
ylab("Frequency")
ggsave("Plots/histogram.png", width = 4, height = 3)Plotting the hypothetical example
Let’s generate the plot of an hypotetical example in which the widening gender gap was just a scaling artifact (Left panel of Figure 4).
# Plot fictious thresholds and equal latent gender
# gaps that imply different manifest gender differences
# Range for which we are going to generate values
lower_limit <- -3
upper_limit <- 3
# Thresholds for illustrative purposes
thresholds <- c(-1.5, -0.25, 0.25, 1.5)
# helper data frames with the areas to-be-shaded
# each represents the area under the curve
# for one of the response options
shade_1 <- data.frame(x = seq(lower_limit, thresholds[1], length.out = 100))
shade_1$y <- dnorm(shade_1$x)
shade_2 <- data.frame(x = seq(thresholds[1], thresholds[2], length.out = 100))
shade_2$y <- dnorm(shade_2$x)
shade_3 <- data.frame(x = seq(thresholds[2], thresholds[3], length.out = 100))
shade_3$y <- dnorm(shade_3$x)
shade_4 <- data.frame(x = seq(thresholds[3], thresholds[4], length.out = 100))
shade_4$y <- dnorm(shade_4$x)
shade_5 <- data.frame(x = seq(thresholds[4], upper_limit, length.out = 100))
shade_5$y <- dnorm(shade_5$x)
# Hypothetical latent values
# for boys and girls
girls_2010 <- 0.3 # girls in 2010
boys_2010 <- girls_2010 + 0.7 # boys in 2010
girls_2023 <- -0.35 # girls in 2023
boys_2023 <- girls_2023 + 0.7 # boys in 2023
# Plot the implied scenario
ggplot(data.frame(x = c(lower_limit, upper_limit)), aes(x = x)) +
stat_function(fun = dnorm, color = "grey") +
theme_classic() +
# Thresholds of the ordinal model
geom_area(data = shade_1, aes(x = x, y = y), fill = "skyblue", alpha = 0.0) +
geom_area(data = shade_2, aes(x = x, y = y), fill = "skyblue", alpha = 0.25) +
geom_area(data = shade_3, aes(x = x, y = y), fill = "skyblue", alpha = 0.5) +
geom_area(data = shade_4, aes(x = x, y = y), fill = "skyblue", alpha = 0.75) +
geom_area(data = shade_5, aes(x = x, y = y), fill = "skyblue", alpha = 1) +
# Gender gap 2010
geom_segment(aes(x = girls_2010, xend = boys_2010, y = 0.25, yend = 0.25)) +
# Gender gap 2023
geom_segment(aes(x = girls_2023, xend = boys_2023, y = 0.20, yend = 0.20)) +
xlab("Latent life satisfaction") +
ylab("Density")## Warning in geom_segment(aes(x = girls_2010, xend = boys_2010, y = 0.25, : All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.## Warning in geom_segment(aes(x = girls_2023, xend = boys_2023, y = 0.2, yend = 0.2)): All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
ggsave("Plots/ordinal_hypothetical.png", width = 4, height = 3)## Warning in geom_segment(aes(x = girls_2010, xend = boys_2010, y = 0.25, : All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
## All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.Fitting the ordinal model on our data
Now let’s analyze our actual data to check whether such scaling artifacts should keep us up at night.
# The model we fit wants a factor as outcome
combined$satis_factor <- as.factor(combined$satis)
# Simple ordered probit model that estimates the means of life satisfaction by gender and survey year
ordinal <- polr(satis_factor ~ as.factor(year) + gender + gender:as.factor(year),
method = "probit",
Hess = TRUE,
data = combined)
# Loot at coefficients (which happen to be sufficient for our purposes here)
summary(ordinal)## Call:
## polr(formula = satis_factor ~ as.factor(year) + gender + gender:as.factor(year),
## data = combined, Hess = TRUE, method = "probit")
##
## Coefficients:
## Value Std. Error t value
## as.factor(year)2015 -0.06711 0.05482 -1.224
## as.factor(year)2023 -0.46720 0.05412 -8.632
## gendermale 0.28455 0.05861 4.855
## as.factor(year)2015:gendermale 0.08843 0.08032 1.101
## as.factor(year)2023:gendermale 0.32845 0.07751 4.237
##
## Intercepts:
## Value Std. Error t value
## 1|2 -2.3149 0.0632 -36.6400
## 2|3 -1.5299 0.0470 -32.5615
## 3|4 -0.5206 0.0424 -12.2789
## 4|5 0.9324 0.0433 21.5264
##
## Residual Deviance: 11285.03
## AIC: 11303.03# polr in combination with marginaleffects does not natively support predictions and comparisons on the latent scale
# we don't need those here yet so we're fine (later we will switch to brms because of that issue)
# However, in principle, we could still look at effects on the "natural" scale
# Which are probabilities of belonging to particular "groups" (response categories)
# Look at gender gaps in probabilities for the individual response categories
comps <- avg_comparisons(ordinal, variables = "gender", by = "year")
comps##
## Group year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 1 2010 -0.00564 0.00141 -4.01 < 0.001 14.0 -0.008398 -0.00288
## 2 2010 -0.02258 0.00486 -4.65 < 0.001 18.2 -0.032095 -0.01306
## 3 2010 -0.06274 0.01288 -4.87 < 0.001 19.8 -0.087984 -0.03750
## 4 2010 0.00798 0.00399 2.00 0.04564 4.5 0.000155 0.01581
## 5 2010 0.08297 0.01719 4.83 < 0.001 19.5 0.049288 0.11666
## 1 2015 -0.00791 0.00157 -5.04 < 0.001 21.1 -0.010981 -0.00483
## 2 2015 -0.03066 0.00480 -6.39 < 0.001 32.5 -0.040070 -0.02125
## 3 2015 -0.08226 0.01208 -6.81 < 0.001 36.6 -0.105942 -0.05859
## 4 2015 0.01412 0.00475 2.97 0.00294 8.4 0.004817 0.02343
## 5 2015 0.10671 0.01612 6.62 < 0.001 34.7 0.075118 0.13830
## 1 2023 -0.02539 0.00358 -7.10 < 0.001 39.6 -0.032396 -0.01838
## 2 2023 -0.07168 0.00690 -10.38 < 0.001 81.5 -0.085205 -0.05815
## 3 2023 -0.12907 0.01095 -11.79 < 0.001 104.1 -0.150527 -0.10761
## 4 2023 0.09119 0.00929 9.82 < 0.001 73.2 0.072989 0.10940
## 5 2023 0.13494 0.01195 11.30 < 0.001 95.9 0.111522 0.15835
##
## Term: gender
## Type: probs
## Comparison: male - female# The interaction of interest here should, on this scale, result in differences in the gender differences in the probabilities of belongig into certain response categories
# (not exactly intuitive)
# Group 1
# Effect of gender on probability of reporting lowest outcome, 2010 vs 2023
print(avg_comparisons(ordinal, variables = "gender", by = "year",
hypothesis = "b1 = b11"))## Warning:
## It is essential to check the order of estimates when specifying hypothesis tests using positional indices like b1, b2, etc. The indices of estimates can change depending on the order of rows in the original dataset, user-supplied arguments, model-fitting package, and version of `marginaleffects`.
##
## It is also good practice to use assertions that ensure the order of estimates is consistent across different runs of the same code. Example:
##
## ```r
## mod <- lm(mpg ~ am * carb, data = mtcars)
##
## # assertion for safety
## p <- avg_predictions(mod, by = 'carb')
## stopifnot(p$carb[1] != 1 || p$carb[2] != 2)
##
## # hypothesis test
## avg_predictions(mod, by = 'carb', hypothesis = 'b1 - b2 = 0')
## ```
##
## Disable this warning with: `options(marginaleffects_safe = FALSE)`
## This warning appears once per session.##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b11 0.0197 0.00342 5.77 <0.001 26.9 0.013 0.0265
##
## Type: probs# Group 2
# Effect of gender on probability of reporting second lowest outcome, 2010 vs 2023
print(avg_comparisons(ordinal, variables = "gender", by = "year",
hypothesis = "b2 = b12"))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b12 0.0491 0.00806 6.09 <0.001 29.7 0.0333 0.0649
##
## Type: probs# Group 3
# Effect of gender on probability of reporting middle outcome, 2010 vs 2023
print(avg_comparisons(ordinal, variables = "gender", by = "year",
hypothesis = "b3 = b13"))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b3=b13 0.0663 0.0165 4.03 <0.001 14.1 0.0341 0.0986
##
## Type: probs# Group 4
# Effect of gender on probability of reporting second highest outcome, 2010 vs 2023
print(avg_comparisons(ordinal, variables = "gender", by = "year",
hypothesis = "b4 = b14"))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b4=b14 -0.0832 0.00997 -8.35 <0.001 53.7 -0.103 -0.0637
##
## Type: probs# Group 5
# Effect of gender on probability of reporting highest outcome, 2010 vs 2023
print(avg_comparisons(ordinal, variables = "gender", by = "year",
hypothesis = "b5 = b15"))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b5=b15 -0.052 0.0209 -2.48 0.013 6.3 -0.093 -0.011
##
## Type: probsIllustrating the model results
Let’s plot these results, see Panel B of Figure 4.
# Plot the model-implied thresholds and the model-implied means
# Here are the estimated thresholds
ordinal$zeta## 1|2 2|3 3|4 4|5
## -2.3149303 -1.5298668 -0.5205751 0.9324412# Range of the normal distribution that we will plot
lower_limit <- -3
upper_limit <- 3
# helper data frames with the areas to-be-shaded
# between the thresholds
shade_1 <- data.frame(x = seq(lower_limit, ordinal$zeta[1], length.out = 100))
shade_1$y <- dnorm(shade_1$x)
shade_2 <- data.frame(x = seq(ordinal$zeta[1], ordinal$zeta[2], length.out = 100))
shade_2$y <- dnorm(shade_2$x)
shade_3 <- data.frame(x = seq(ordinal$zeta[2], ordinal$zeta[3], length.out = 100))
shade_3$y <- dnorm(shade_3$x)
shade_4 <- data.frame(x = seq(ordinal$zeta[3], ordinal$zeta[4], length.out = 100))
shade_4$y <- dnorm(shade_4$x)
shade_5 <- data.frame(x = seq(ordinal$zeta[4], upper_limit, length.out = 100))
shade_5$y <- dnorm(shade_5$x)
# model-implied latent values
girls_2010 <- 0 # girls in 2010, reference group
boys_2010 <- as.numeric(coefficients(ordinal)["gendermale"]) # boys in 2010
girls_2023 <- as.numeric(coefficients(ordinal)["as.factor(year)2023"]) # girls in 2023
boys_2023 <- as.numeric(coefficients(ordinal)["gendermale"] + coefficients(ordinal)["as.factor(year)2023"] + coefficients(ordinal)["as.factor(year)2023:gendermale"]) # boys in 2023
# Plot the distribution
ggplot(data.frame(x = c(lower_limit, upper_limit)), aes(x = x)) +
stat_function(fun = dnorm, color = "grey") +
theme_classic() +
# Thresholds of the ordinal model
geom_area(data = shade_1, aes(x = x, y = y), fill = "skyblue", alpha = 0.0) +
geom_area(data = shade_2, aes(x = x, y = y), fill = "skyblue", alpha = 0.25) +
geom_area(data = shade_3, aes(x = x, y = y), fill = "skyblue", alpha = 0.5) +
geom_area(data = shade_4, aes(x = x, y = y), fill = "skyblue", alpha = 0.75) +
geom_area(data = shade_5, aes(x = x, y = y), fill = "skyblue", alpha = 1) +
#geom_vline(xintercept = ordinal$zeta, color = "lightgrey") + # optional, add line for the actual thresholds
# Gender gap 2010
geom_segment(aes(x = girls_2010, xend = boys_2010, y = 0.25, yend = 0.25)) +
# Gender gap 2023
geom_segment(aes(x = girls_2023, xend = boys_2023, y = 0.20, yend = 0.20)) +
xlab("Latent life satisfaction") +
ylab("Density")## Warning in geom_segment(aes(x = girls_2010, xend = boys_2010, y = 0.25, : All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.## Warning in geom_segment(aes(x = girls_2023, xend = boys_2023, y = 0.2, yend = 0.2)): All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
ggsave("Plots/ordinal_empirical.png", width = 4, height = 3)## Warning in geom_segment(aes(x = girls_2010, xend = boys_2010, y = 0.25, : All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.
## All aesthetics have length 1, but the data has 2 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
## a single row.Additionally freeing thresholds
In footnote 4, we mention that it is possible to free the thresholds across gender and still arrive at the same results. Here are the corresponding analyses. We can also free the thresholds by year (not reported in the manuscript) and arrive at similar conclusions, however, it’s a bit harder to defend why one would free thresholds across year but not across gender.
Note that in either case, this introduces huge uncertainty about the “main effect” of the variable for which we are freeing the thresholds. That is because any affect on the mean of the latent variable could be arbitrarily re-expressed as a shifting of all thresholds; essentially the model is underidentified (there is an infinite number of solutions that would produce precisely the same observed data). We can still get point estimates because of Bayesian magic and priors, but the uncertainty explodes. However, this does not really affect the uncertainty regarding the interaction of interest. Here’s an intuition for that: No matter which combination of thresholds and main effect we pick, once we commit ourselves, it’s the same across years (if freeing by gender) or the same across genders (if freeing by year). Thus, the contrasts across that second variable are still identified.
But notice that we could not possibly free the threesholds by both gender and both year; then the uncertainty explosion would affect the interaction of interest. One way to think about this: If we say “oh, in every year, girls and boys may use the scale in arbitrarily different ways that change over time”, that means we essentially cannot say anything anymore about the gender differences and how they change over time. They could be literally anything; we can just freely pick any pattern of gender gaps and then just adjust the gender and year specific thresholds to produce the observed data.
# We switch to a Bayesian framework to make the magic happen
library(brms)## Loading required package: Rcpp## Loading 'brms' package (version 2.22.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following object is masked from 'package:stats':
##
## ar# Refit simple ordinal model in brms just to check its essentially the same
# ordinal_b <- brm(satis ~ as.factor(year) + gender + gender:as.factor(year),
# family = cumulative(probit),
# data = combined,
# cores = 4,
# seed = 1)
# saveRDS(ordinal_b, file = "Models/ordinal_brms")
ordinal_b <- readRDS(file = "Models/ordinal_brms")
# Loot at coefficients
summary(ordinal_b)## Family: cumulative
## Links: mu = probit; disc = identity
## Formula: satis ~ as.factor(year) + gender + gender:as.factor(year)
## Data: combined (Number of observations: 4764)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1] -2.32 0.06 -2.44 -2.20 1.00 1826
## Intercept[2] -1.53 0.05 -1.62 -1.44 1.00 2153
## Intercept[3] -0.52 0.04 -0.60 -0.43 1.00 2030
## Intercept[4] 0.93 0.04 0.85 1.02 1.00 2001
## as.factoryear2015 -0.07 0.06 -0.17 0.04 1.00 1750
## as.factoryear2023 -0.47 0.05 -0.57 -0.36 1.00 2042
## gendermale 0.29 0.06 0.18 0.40 1.00 1612
## as.factoryear2015:gendermale 0.09 0.08 -0.07 0.24 1.00 1609
## as.factoryear2023:gendermale 0.32 0.08 0.17 0.47 1.00 1767
## Tail_ESS
## Intercept[1] 2207
## Intercept[2] 2619
## Intercept[3] 2643
## Intercept[4] 2758
## as.factoryear2015 2497
## as.factoryear2023 2650
## gendermale 1750
## as.factoryear2015:gendermale 2253
## as.factoryear2023:gendermale 1857
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Free thresholds by gender
# ordinal_b_free <- brm(satis| thres(gr = gender) ~ as.factor(year) + gender + gender:as.factor(year),
# family = cumulative(probit),
# data = combined,
# cores = 4,
# seed = 1)
# saveRDS(ordinal_b_free, file = "Models/ordinal_brms_free_gender")
ordinal_b_free <- readRDS(file = "Models/ordinal_brms_free_gender")
summary(ordinal_b_free)## Family: cumulative
## Links: mu = probit; disc = identity
## Formula: satis | thres(gr = gender) ~ as.factor(year) + gender + gender:as.factor(year)
## Data: combined (Number of observations: 4764)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[female,1] -2.35 0.07 -2.50 -2.21 1.00 1848
## Intercept[female,2] -1.53 0.05 -1.64 -1.43 1.00 2031
## Intercept[female,3] -0.45 0.04 -0.54 -0.36 1.00 2053
## Intercept[female,4] 0.87 0.05 0.78 0.96 1.00 2423
## Intercept[male,1] -1.41 1.29 -4.02 1.02 1.00 1721
## Intercept[male,2] -0.68 1.29 -3.27 1.73 1.00 1723
## Intercept[male,3] 0.20 1.29 -2.39 2.62 1.00 1722
## Intercept[male,4] 1.80 1.29 -0.80 4.21 1.00 1723
## as.factoryear2015 -0.06 0.06 -0.16 0.05 1.00 1987
## as.factoryear2023 -0.45 0.06 -0.56 -0.34 1.00 1932
## gendermale 1.09 1.29 -1.51 3.50 1.00 1745
## as.factoryear2015:gendermale 0.08 0.08 -0.08 0.24 1.00 2037
## as.factoryear2023:gendermale 0.30 0.08 0.15 0.46 1.00 1987
## Tail_ESS
## Intercept[female,1] 2191
## Intercept[female,2] 2507
## Intercept[female,3] 2811
## Intercept[female,4] 2679
## Intercept[male,1] 2326
## Intercept[male,2] 2265
## Intercept[male,3] 2266
## Intercept[male,4] 2306
## as.factoryear2015 2399
## as.factoryear2023 2457
## gendermale 2281
## as.factoryear2015:gendermale 2299
## as.factoryear2023:gendermale 2561
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# now there is huge uncertainties about the thresholds in men and the
# effect of male gender
# Free the thresholds by year
# ordinal_b_free_year <- brm(satis| thres(gr = as.factor(year)) ~ as.factor(year) + gender + gender:as.factor(year),
# family = cumulative(probit),
# data = combined,
# cores = 4,
# seed = 1)
# saveRDS(ordinal_b_free_year, file = "Models/ordinal_brms_free_year")
ordinal_b_free_year <- readRDS(file = "Models/ordinal_brms_free_year")
summary(ordinal_b_free_year)## Family: cumulative
## Links: mu = probit; disc = identity
## Formula: satis | thres(gr = as.factor(year)) ~ as.factor(year) + gender + gender:as.factor(year)
## Data: combined (Number of observations: 4764)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[2010,1] -2.17 0.10 -2.38 -1.99 1.00 2167
## Intercept[2010,2] -1.38 0.06 -1.50 -1.27 1.00 3809
## Intercept[2010,3] -0.51 0.05 -0.61 -0.43 1.00 2816
## Intercept[2010,4] 0.86 0.05 0.77 0.96 1.00 2644
## Intercept[2015,1] -1.41 1.33 -4.09 1.12 1.00 2435
## Intercept[2015,2] -0.69 1.33 -3.34 1.81 1.00 2438
## Intercept[2015,3] 0.35 1.32 -2.29 2.87 1.00 2434
## Intercept[2015,4] 1.84 1.32 -0.80 4.37 1.00 2428
## Intercept[2023,1] -1.52 1.30 -4.05 1.10 1.00 2336
## Intercept[2023,2] -0.69 1.30 -3.24 1.89 1.00 2327
## Intercept[2023,3] 0.39 1.30 -2.15 2.97 1.00 2328
## Intercept[2023,4] 1.89 1.30 -0.64 4.47 1.00 2335
## as.factoryear2015 0.82 1.33 -1.83 3.35 1.00 2436
## as.factoryear2023 0.43 1.30 -2.10 3.02 1.00 2319
## gendermale 0.27 0.06 0.15 0.38 1.00 1762
## as.factoryear2015:gendermale 0.11 0.08 -0.05 0.27 1.00 2147
## as.factoryear2023:gendermale 0.37 0.08 0.22 0.52 1.00 2146
## Tail_ESS
## Intercept[2010,1] 2305
## Intercept[2010,2] 3372
## Intercept[2010,3] 3034
## Intercept[2010,4] 2546
## Intercept[2015,1] 2336
## Intercept[2015,2] 2233
## Intercept[2015,3] 2173
## Intercept[2015,4] 2198
## Intercept[2023,1] 2224
## Intercept[2023,2] 2226
## Intercept[2023,3] 2225
## Intercept[2023,4] 2239
## as.factoryear2015 2198
## as.factoryear2023 2233
## gendermale 2272
## as.factoryear2015:gendermale 2338
## as.factoryear2023:gendermale 2547
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).