Home | Cleaning and preparation of the Youth Surveys | Explanandum and alternative satisfaction measures | Ordinal analysis | Add demographic variables | Survey mode effects | Re-check with ordinal model | Returning to alternative satisfaction measures | Miscellaneous

Check whether results are robust with ordinal model

Do we get the same pattern if we once again use an ordinal model?

# Fit model in ordinal
# We are using brms here simply because polr() does not support
# the postestimation routine of marginaleffects on the link scale
library(splines)
# explanandum_paper_only_ordinal <- brm(satis ~ as.factor(year) + gender + as.factor(year):gender +
#                                     bs(age, df = 3) + as.factor(year):bs(age, df = 3) + gender:bs(age, df = 3) + as.factor(year):gender:bs(age, df = 3) +
#                                     as.factor(schooltype) + as.factor(year):as.factor(schooltype) + gender:as.factor(schooltype) + as.factor(year):gender:as.factor(schooltype) +
#                                     as.factor(mig_lang) + as.factor(year):as.factor(mig_lang) + gender:as.factor(mig_lang) + as.factor(year):gender:as.factor(mig_lang),
#                                   family = cumulative("probit"),
#                                   cores = 4,
#                    data = combined[combined$source == "paper",])
#  
# 
# saveRDS(explanandum_paper_only_ordinal, file = "explanandum_paper_only_ordinal")
explanandum_paper_only_ordinal <- readRDS("explanandum_paper_only_ordinal")

# Demographics like in 2010
temp1 <- combined[combined$year == 2010,]
temp1$year <- 2010
temp2 <- combined[combined$year == 2010,]
temp2$year <- 2015
temp3 <- combined[combined$year == 2010,]
temp3$year <- 2023
all2010 <- rbind(temp1, temp2, temp3)

# Predicted values for all combinations of gender and year
pred <- avg_predictions(explanandum_paper_only_ordinal,
                 by = c("gender", "year"),
                 newdata = all2010,
                 vcov = ~ unique_classroom,
                 type = "link")

## Warning: The `vcov` argument is not supported for models of this class.

print(pred)

## 
##  gender year Estimate   2.5 %  97.5 %
##  female 2010  -0.1448 -0.3508  0.0632
##  female 2015  -0.1820 -0.4153  0.0659
##  female 2023  -0.4459 -0.7272 -0.1758
##  male   2010   0.1413 -0.0892  0.3807
##  male   2015   0.2346 -0.0107  0.4874
##  male   2023   0.0903 -0.2155  0.3932
## 
## Type: link

# Counterfactual comparisons: latent gender gaps per year
comps <- avg_comparisons(explanandum_paper_only_ordinal, variables = "gender", by = "year",
                 newdata = all2010,
                 vcov = ~ unique_classroom,
                 type = "link")

## Warning: The `vcov` argument is not supported for models of this class.

print(comps)

## 
##  year Estimate 2.5 % 97.5 %
##  2010    0.293 0.179  0.411
##  2015    0.424 0.280  0.566
##  2023    0.543 0.279  0.824
## 
## Term: gender
## Type: link
## Comparison: male - female

# Let's compare the 2010 and 2023 gender gaps
max_comp <- avg_comparisons(explanandum_paper_only_ordinal, variables = "gender", by = "year",
                hypothesis = "b1 = b3",
                 newdata = all2010,
                vcov = ~ unique_classroom,
                 type = "link")

## Warning: The `vcov` argument is not supported for models of this class.

## 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.

print(max_comp)

## 
##  Hypothesis Estimate  2.5 % 97.5 %
##       b1=b3    -0.25 -0.545 0.0357
## 
## Type: link

# And lets look at how these vary by migration background

# Widening of gender gap for students without migration background
print(avg_comparisons(explanandum_paper_only_ordinal, variables = "gender", by = c("year", "mig_lang"),
                 newdata = all2010,
                 vcov = ~ unique_classroom,
                 hypothesis = "b1 = b5",
                 type = "link"))

## Warning: The `vcov` argument is not supported for models of this class.

## 
##  Hypothesis Estimate  2.5 % 97.5 %
##       b1=b5   -0.153 -0.456  0.154
## 
## Type: link

# Widening of gender gap for students with migration background
print(avg_comparisons(explanandum_paper_only_ordinal, variables = "gender", by = c("year", "mig_lang"),
                 newdata = all2010,
                 vcov = ~ unique_classroom,
                 hypothesis = "b2 = b6",
                 type = "link"))

## Warning: The `vcov` argument is not supported for models of this class.

## 
##  Hypothesis Estimate 2.5 % 97.5 %
##       b2=b6    -1.16 -1.79 -0.493
## 
## Type: link

# Comparing the widening (triple-interaction gender, year, migration background)
print(avg_comparisons(explanandum_paper_only_ordinal, variables = "gender", by = c("year", "mig_lang"),
                 newdata = all2010,
                 vcov = ~ unique_classroom,
                 hypothesis = "b1 - b5 = b2 - b6",
                 type = "link"))

## Warning: The `vcov` argument is not supported for models of this class.

## 
##   Hypothesis Estimate 2.5 % 97.5 %
##  b1-b5=b2-b6        1 0.293   1.69
## 
## Type: link