Add demographic variables
Julia Rohrer
2026-03-16
Add demographic variables
Here, we are going to add in demographic variables in an iterative manner, also looking at their associations with the outcome to gauge whether they could plausibly explain away anything. If you just want to see the results reported in the manuscript, you can immediately skip to the final model at the bottom of the page.
Check for differences between the years
Age
load("prep.RData")
library(marginaleffects)
library(ggplot2)
library(patchwork)
library(splines)
# Are there age differences between the years?
age_diffs <- lm(age ~ as.factor(year)*gender, data = combined)
print(predictions(age_diffs,
by = c("gender", "year")))##
## gender year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 2010 14.1 0.0574 245 <0.001 Inf 13.9 14.2
## female 2015 14.8 0.0518 287 <0.001 Inf 14.7 14.9
## female 2023 14.8 0.0504 293 <0.001 Inf 14.7 14.9
## male 2010 14.2 0.0583 243 <0.001 Inf 14.0 14.3
## male 2015 14.7 0.0569 258 <0.001 Inf 14.6 14.8
## male 2023 14.7 0.0509 289 <0.001 Inf 14.6 14.8
##
## Type: responseprint(avg_comparisons(age_diffs))## Warning: The `year` variable is treated as a categorical (factor) variable, but
## the original data is of class numeric. It is safer and faster to convert such
## variables to factor before fitting the model and calling a `marginaleffects`
## function. This warning appears once per session.##
## Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## gender male - female -0.0334 0.0441 -0.756 0.45 1.2 -0.120 0.0531
## year 2015 - 2010 0.6706 0.0561 11.957 <0.001 107.1 0.561 0.7805
## year 2023 - 2010 0.6446 0.0544 11.850 <0.001 105.2 0.538 0.7512
##
## Type: responseprint(avg_comparisons(age_diffs, variable = "gender", by = "year"))##
## year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0.0931 0.0818 1.138 0.2552 2.0 -0.0673 0.2535
## 2015 -0.1275 0.0769 -1.657 0.0976 3.4 -0.2782 0.0233
## 2023 -0.0481 0.0717 -0.671 0.5023 1.0 -0.1885 0.0924
##
## Term: gender
## Type: response
## Comparison: male - femaleSchooltype
# Recode to turn it into a Gymnasium dummy
combined$gymnasium <- ifelse(combined$schooltype == 2, 1, 0)
# Are there age differences between the years?
# Fit a simple linear model to estimate probabilities
schooltype_diffs <- lm(gymnasium ~ as.factor(year)*gender, data = combined)
# Differences between years
print(predictions(schooltype_diffs,
by = c("year"), vcov = "HC3"))##
## year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0.536 0.0134 39.9 <0.001 Inf 0.509 0.562
## 2015 0.588 0.0124 47.6 <0.001 Inf 0.564 0.612
## 2023 0.654 0.0112 58.4 <0.001 Inf 0.633 0.676
##
## Type: response# Differences between genders
print(predictions(schooltype_diffs,
by = c("gender"), vcov = "HC3"))##
## gender Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 0.623 0.00969 64.2 <0.001 Inf 0.604 0.642
## male 0.571 0.01031 55.4 <0.001 Inf 0.551 0.592
##
## Type: response# Differences between genders and years
print(predictions(schooltype_diffs,
by = c("gender", "year"), vcov = "HC3"))##
## gender year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 2010 0.543 0.0188 28.9 <0.001 606.9 0.507 0.580
## female 2015 0.630 0.0164 38.3 <0.001 Inf 0.597 0.662
## female 2023 0.677 0.0155 43.6 <0.001 Inf 0.647 0.707
## male 2010 0.528 0.0192 27.5 <0.001 552.1 0.490 0.566
## male 2015 0.538 0.0187 28.8 <0.001 602.8 0.501 0.574
## male 2023 0.632 0.0162 39.1 <0.001 Inf 0.600 0.663
##
## Type: response# Gender gap by year
print(avg_comparisons(schooltype_diffs, variable = "gender", by = "year", vcov = "HC3"))##
## year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 -0.0154 0.0269 -0.575 0.565 0.8 -0.0681 0.03720
## 2015 -0.0918 0.0249 -3.688 <0.001 12.1 -0.1406 -0.04303
## 2023 -0.0453 0.0224 -2.024 0.043 4.5 -0.0893 -0.00143
##
## Term: gender
## Type: response
## Comparison: male - female# Are the differences in migration background between the years?
# Simple linear probability model predicting a migration background
mig_lang_diffs <- lm(mig_lang ~ as.factor(year)*gender, data = combined)
# Means by survey year
print(predictions(mig_lang_diffs,
by = c("year"), vcov = "HC3"))##
## year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0.1005 0.00810 12.4 <0.001 115.1 0.0846 0.116
## 2015 0.0938 0.00735 12.8 <0.001 121.5 0.0794 0.108
## 2023 0.1292 0.00791 16.3 <0.001 197.0 0.1137 0.145
##
## Type: response# Means by gender
print(predictions(mig_lang_diffs,
by = c("gender"), vcov = "HC3"))##
## gender Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 0.108 0.00623 17.3 <0.001 220.3 0.0956 0.120
## male 0.111 0.00656 16.9 <0.001 209.7 0.0978 0.123
##
## Type: response# Means by gender and year
print(predictions(mig_lang_diffs,
by = c("gender", "year"), vcov = "HC3"))##
## gender year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 2010 0.1038 0.01152 9.01 <0.001 62.1 0.0813 0.126
## female 2015 0.0914 0.00982 9.31 <0.001 66.1 0.0722 0.111
## female 2023 0.1264 0.01103 11.46 <0.001 98.6 0.1048 0.148
## male 2010 0.0971 0.01137 8.54 <0.001 56.0 0.0748 0.119
## male 2015 0.0966 0.01107 8.73 <0.001 58.4 0.0749 0.118
## male 2023 0.1321 0.01135 11.65 <0.001 101.7 0.1099 0.154
##
## Type: response# Gender gaps per year
print(avg_comparisons(mig_lang_diffs, variable = "gender", by = "year", vcov = "HC3"))##
## year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 -0.00678 0.0162 -0.419 0.675 0.6 -0.0385 0.0249
## 2015 0.00520 0.0148 0.352 0.725 0.5 -0.0238 0.0342
## 2023 0.00577 0.0158 0.364 0.716 0.5 -0.0252 0.0368
##
## Term: gender
## Type: response
## Comparison: male - femaleLooking more closely at age
Age and life satisfaction
# Visualize the association between age and satisfaction, by gender and year
# Fully categorical model
# This allows the effect of every year of age to vary freely
# A lot of flexibility but also huge error bars
add_age <- lm(satis ~ gender*as.factor(year)*as.factor(age),
data = combined)
summary(add_age)##
## Call:
## lm(formula = satis ~ gender * as.factor(year) * as.factor(age),
## data = combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.02198 -0.47059 0.06736 0.52632 1.70238
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 3.82524 0.08070 47.400
## gendermale 0.37476 0.12002 3.122
## as.factor(year)2015 -0.09797 0.19236 -0.509
## as.factor(year)2023 -0.14408 0.12742 -1.131
## as.factor(age)13 -0.01942 0.09884 -0.196
## as.factor(age)14 0.02091 0.10585 0.198
## as.factor(age)15 -0.17818 0.11023 -1.617
## as.factor(age)16 0.05937 0.12293 0.483
## as.factor(age)17 -0.03037 0.15399 -0.197
## as.factor(age)18 -0.49191 0.22635 -2.173
## gendermale:as.factor(year)2015 0.35951 0.31064 1.157
## gendermale:as.factor(year)2023 0.03499 0.19842 0.176
## gendermale:as.factor(age)13 0.01842 0.14491 0.127
## gendermale:as.factor(age)14 -0.33130 0.15314 -2.163
## gendermale:as.factor(age)15 -0.06142 0.16335 -0.376
## gendermale:as.factor(age)16 -0.43848 0.18171 -2.413
## gendermale:as.factor(age)17 -0.38963 0.21218 -1.836
## gendermale:as.factor(age)18 -0.11718 0.29936 -0.391
## as.factor(year)2015:as.factor(age)13 0.18589 0.21084 0.882
## as.factor(year)2023:as.factor(age)13 -0.18806 0.15302 -1.229
## as.factor(year)2015:as.factor(age)14 0.04348 0.21258 0.205
## as.factor(year)2023:as.factor(age)14 -0.25819 0.15604 -1.655
## as.factor(year)2015:as.factor(age)15 0.12768 0.21454 0.595
## as.factor(year)2023:as.factor(age)15 -0.20536 0.16083 -1.277
## as.factor(year)2015:as.factor(age)16 -0.01556 0.22281 -0.070
## as.factor(year)2023:as.factor(age)16 -0.33785 0.17128 -1.972
## as.factor(year)2015:as.factor(age)17 -0.14236 0.24557 -0.580
## as.factor(year)2023:as.factor(age)17 -0.18020 0.20329 -0.886
## as.factor(year)2015:as.factor(age)18 0.57714 0.35164 1.641
## as.factor(year)2023:as.factor(age)18 0.18575 0.26509 0.701
## gendermale:as.factor(year)2015:as.factor(age)13 -0.54720 0.33527 -1.632
## gendermale:as.factor(year)2023:as.factor(age)13 0.12012 0.23146 0.519
## gendermale:as.factor(year)2015:as.factor(age)14 -0.19995 0.33557 -0.596
## gendermale:as.factor(year)2023:as.factor(age)14 0.41032 0.23520 1.745
## gendermale:as.factor(year)2015:as.factor(age)15 -0.33374 0.34017 -0.981
## gendermale:as.factor(year)2023:as.factor(age)15 0.29850 0.24297 1.229
## gendermale:as.factor(year)2015:as.factor(age)16 -0.09230 0.35332 -0.261
## gendermale:as.factor(year)2023:as.factor(age)16 0.40457 0.25882 1.563
## gendermale:as.factor(year)2015:as.factor(age)17 0.14769 0.37907 0.390
## gendermale:as.factor(year)2023:as.factor(age)17 0.35374 0.29212 1.211
## gendermale:as.factor(year)2015:as.factor(age)18 -0.52049 0.52404 -0.993
## gendermale:as.factor(year)2023:as.factor(age)18 -0.10000 0.37687 -0.265
## Pr(>|t|)
## (Intercept) <2e-16 ***
## gendermale 0.0018 **
## as.factor(year)2015 0.6106
## as.factor(year)2023 0.2582
## as.factor(age)13 0.8443
## as.factor(age)14 0.8434
## as.factor(age)15 0.1060
## as.factor(age)16 0.6291
## as.factor(age)17 0.8437
## as.factor(age)18 0.0298 *
## gendermale:as.factor(year)2015 0.2472
## gendermale:as.factor(year)2023 0.8600
## gendermale:as.factor(age)13 0.8988
## gendermale:as.factor(age)14 0.0306 *
## gendermale:as.factor(age)15 0.7069
## gendermale:as.factor(age)16 0.0159 *
## gendermale:as.factor(age)17 0.0664 .
## gendermale:as.factor(age)18 0.6955
## as.factor(year)2015:as.factor(age)13 0.3780
## as.factor(year)2023:as.factor(age)13 0.2191
## as.factor(year)2015:as.factor(age)14 0.8379
## as.factor(year)2023:as.factor(age)14 0.0981 .
## as.factor(year)2015:as.factor(age)15 0.5518
## as.factor(year)2023:as.factor(age)15 0.2017
## as.factor(year)2015:as.factor(age)16 0.9443
## as.factor(year)2023:as.factor(age)16 0.0486 *
## as.factor(year)2015:as.factor(age)17 0.5621
## as.factor(year)2023:as.factor(age)17 0.3754
## as.factor(year)2015:as.factor(age)18 0.1008
## as.factor(year)2023:as.factor(age)18 0.4835
## gendermale:as.factor(year)2015:as.factor(age)13 0.1027
## gendermale:as.factor(year)2023:as.factor(age)13 0.6038
## gendermale:as.factor(year)2015:as.factor(age)14 0.5513
## gendermale:as.factor(year)2023:as.factor(age)14 0.0811 .
## gendermale:as.factor(year)2015:as.factor(age)15 0.3266
## gendermale:as.factor(year)2023:as.factor(age)15 0.2193
## gendermale:as.factor(year)2015:as.factor(age)16 0.7939
## gendermale:as.factor(year)2023:as.factor(age)16 0.1181
## gendermale:as.factor(year)2015:as.factor(age)17 0.6968
## gendermale:as.factor(year)2023:as.factor(age)17 0.2260
## gendermale:as.factor(year)2015:as.factor(age)18 0.3206
## gendermale:as.factor(year)2023:as.factor(age)18 0.7908
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.819 on 4722 degrees of freedom
## Multiple R-squared: 0.07702, Adjusted R-squared: 0.06901
## F-statistic: 9.611 on 41 and 4722 DF, p-value: < 2.2e-16# Generate predictions for combinations of gender and year and age
pred <- predictions(add_age,
by = c("gender", "year", "age"))
# Visualize
# 2010
cat_age_2010 <- ggplot(pred[pred$year == 2010,], aes(x = age, group = gender,
color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point(position = position_dodge(width = .5)) +
geom_line(position = position_dodge(width = .5)) +
geom_errorbar(width = .5, position = position_dodge(width = .5)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("2010") +
ylab("Mean life satisfaction") +
xlab("Age") +
labs(color = "Gender")
# 2015
cat_age_2015 <- ggplot(pred[pred$year == 2015,], aes(x = age, group = gender,
color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point(position = position_dodge(width = .5)) +
geom_line(position = position_dodge(width = .5)) +
geom_errorbar(width = .5, position = position_dodge(width = .5)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("2015") +
ylab("Mean life satisfaction") +
xlab("Age") +
labs(color = "Gender")
# 2023
cat_age_2023 <- ggplot(pred[pred$year == 2023,], aes(x = age, group = gender,
color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point(position = position_dodge(width = .5)) +
geom_line(position = position_dodge(width = .5)) +
geom_errorbar(width = .5, position = position_dodge(width = .5)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("2023") +
ylab("Mean life satisfaction") +
xlab("Age") +
labs(color = "Gender")
# Combine plots
library(patchwork)
combined_cat_age <- (cat_age_2010|cat_age_2015|cat_age_2023) + plot_layout(guides = "collect")
combined_cat_age
ggsave("Plots/satis_by_age_categorical.png", width = 6, height = 3)
# Re-fit the model including only the linear age effect
# This is the least flexible solution
# But also greatly induces the uncertainty in the estimates
add_age_linear <- lm(satis ~ gender*as.factor(year)*age,
data = combined)
summary(add_age_linear)##
## Call:
## lm(formula = satis ~ gender * as.factor(year) * age, data = combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.04914 -0.45410 0.08008 0.51557 1.66722
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.127277 0.286403 14.411 < 2e-16 ***
## gendermale 1.309701 0.402875 3.251 0.00116 **
## as.factor(year)2015 0.374777 0.408654 0.917 0.35914
## as.factor(year)2023 -0.248569 0.373317 -0.666 0.50555
## age -0.024132 0.020256 -1.191 0.23356
## gendermale:as.factor(year)2015 -1.343628 0.604473 -2.223 0.02628 *
## gendermale:as.factor(year)2023 -0.304458 0.538597 -0.565 0.57191
## gendermale:age -0.077005 0.028396 -2.712 0.00671 **
## as.factor(year)2015:age -0.026646 0.028159 -0.946 0.34406
## as.factor(year)2023:age -0.006197 0.025880 -0.239 0.81077
## gendermale:as.factor(year)2015:age 0.097856 0.041602 2.352 0.01870 *
## gendermale:as.factor(year)2023:age 0.041279 0.037247 1.108 0.26782
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8197 on 4752 degrees of freedom
## Multiple R-squared: 0.06961, Adjusted R-squared: 0.06745
## F-statistic: 32.32 on 11 and 4752 DF, p-value: < 2.2e-16# generate predictions for combinations of gender and year and age
pred <- predictions(add_age_linear,
by = c("gender", "year", "age"))
# Plot the results
# 2010
lin_age_2010 <- ggplot(pred[pred$year == 2010,], aes(x = age, group = gender,
color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point(position = position_dodge(width = .5)) +
geom_line(position = position_dodge(width = .5)) +
geom_errorbar(width = .5, position = position_dodge(width = .5)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("2010") +
ylab("Mean life satisfaction") +
xlab("Age") +
labs(color = "Gender")
# 2015
lin_age_2015 <- ggplot(pred[pred$year == 2015,], aes(x = age, group = gender,
color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point(position = position_dodge(width = .5)) +
geom_line(position = position_dodge(width = .5)) +
geom_errorbar(width = .5, position = position_dodge(width = .5)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ylab("Mean life satisfaction") +
xlab("Age") +
labs(color = "Gender") +
ggtitle("2015")
# 2023
lin_age_2023 <- ggplot(pred[pred$year == 2023,], aes(x = age, group = gender,
color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point(position = position_dodge(width = .5)) +
geom_line(position = position_dodge(width = .5)) +
geom_errorbar(width = .5, position = position_dodge(width = .5)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ylab("Mean life satisfaction") +
xlab("Age") +
labs(color = "Gender") +
ggtitle("2023")
# Combine into single plot
combined_lin_age <- (lin_age_2010|lin_age_2015|lin_age_2023) + plot_layout(guides = "collect")
combined_lin_age
ggsave("Plots/satis_by_age_linear.png", width = 6, height = 3)
# Best of both worlds, splines
# Now, we allow for some flexibility
# But we still allow some pooling across years
add_age_splines <- lm(satis ~ gender*as.factor(year)*bs(age, df = 3),
data = combined)
summary(add_age_splines)##
## Call:
## lm(formula = satis ~ gender * as.factor(year) * bs(age, df = 3),
## data = combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.01823 -0.41345 0.09321 0.51493 1.62858
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 3.848565 0.077322 49.773
## gendermale 0.408636 0.114297 3.575
## as.factor(year)2015 -0.019961 0.163407 -0.122
## as.factor(year)2023 -0.160591 0.119535 -1.343
## bs(age, df = 3)1 -0.202410 0.251542 -0.805
## bs(age, df = 3)2 0.189582 0.251148 0.755
## bs(age, df = 3)3 -0.347724 0.207408 -1.677
## gendermale:as.factor(year)2015 0.164888 0.255014 0.647
## gendermale:as.factor(year)2023 0.007901 0.182693 0.043
## gendermale:bs(age, df = 3)1 -0.104318 0.362279 -0.288
## gendermale:bs(age, df = 3)2 -0.615447 0.351343 -1.752
## gendermale:bs(age, df = 3)3 -0.257779 0.278330 -0.926
## as.factor(year)2015:bs(age, df = 3)1 0.338526 0.445707 0.760
## as.factor(year)2023:bs(age, df = 3)1 -0.335873 0.354975 -0.946
## as.factor(year)2015:bs(age, df = 3)2 -0.534900 0.334732 -1.598
## as.factor(year)2023:bs(age, df = 3)2 -0.401170 0.308787 -1.299
## as.factor(year)2015:bs(age, df = 3)3 0.215285 0.320769 0.671
## as.factor(year)2023:bs(age, df = 3)3 0.064902 0.247241 0.263
## gendermale:as.factor(year)2015:bs(age, df = 3)1 -0.828851 0.680909 -1.217
## gendermale:as.factor(year)2023:bs(age, df = 3)1 0.406889 0.527842 0.771
## gendermale:as.factor(year)2015:bs(age, df = 3)2 0.815900 0.486187 1.678
## gendermale:as.factor(year)2023:bs(age, df = 3)2 0.694505 0.439340 1.581
## gendermale:as.factor(year)2015:bs(age, df = 3)3 -0.090340 0.477675 -0.189
## gendermale:as.factor(year)2023:bs(age, df = 3)3 0.056690 0.355331 0.160
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## gendermale 0.000353 ***
## as.factor(year)2015 0.902779
## as.factor(year)2023 0.179186
## bs(age, df = 3)1 0.421046
## bs(age, df = 3)2 0.450370
## bs(age, df = 3)3 0.093701 .
## gendermale:as.factor(year)2015 0.517933
## gendermale:as.factor(year)2023 0.965507
## gendermale:bs(age, df = 3)1 0.773398
## gendermale:bs(age, df = 3)2 0.079890 .
## gendermale:bs(age, df = 3)3 0.354409
## as.factor(year)2015:bs(age, df = 3)1 0.447576
## as.factor(year)2023:bs(age, df = 3)1 0.344101
## as.factor(year)2015:bs(age, df = 3)2 0.110110
## as.factor(year)2023:bs(age, df = 3)2 0.193946
## as.factor(year)2015:bs(age, df = 3)3 0.502156
## as.factor(year)2023:bs(age, df = 3)3 0.792945
## gendermale:as.factor(year)2015:bs(age, df = 3)1 0.223562
## gendermale:as.factor(year)2023:bs(age, df = 3)1 0.440832
## gendermale:as.factor(year)2015:bs(age, df = 3)2 0.093382 .
## gendermale:as.factor(year)2023:bs(age, df = 3)2 0.113993
## gendermale:as.factor(year)2015:bs(age, df = 3)3 0.850003
## gendermale:as.factor(year)2023:bs(age, df = 3)3 0.873248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8197 on 4740 degrees of freedom
## Multiple R-squared: 0.07199, Adjusted R-squared: 0.06748
## F-statistic: 15.99 on 23 and 4740 DF, p-value: < 2.2e-16# generate predictions for combinations of gender and year and age
pred <- predictions(add_age_splines,
by = c("gender", "year", "age"))
# Plot the results
# 2010
spline_age_2010 <- ggplot(pred[pred$year == 2010,], aes(x = age, group = gender,
color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point(position = position_dodge(width = .5)) +
geom_line(position = position_dodge(width = .5)) +
geom_errorbar(width = .5, position = position_dodge(width = .5)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("2010") +
ylab("Mean life satisfaction") +
xlab("Age") +
labs(color = "Gender")
# 2015
spline_age_2015 <- ggplot(pred[pred$year == 2015,], aes(x = age, group = gender,
color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point(position = position_dodge(width = .5)) +
geom_line(position = position_dodge(width = .5)) +
geom_errorbar(width = .5, position = position_dodge(width = .5)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ylab("Mean life satisfaction") +
xlab("Age") +
labs(color = "Gender") +
ggtitle("2015")
# 2023
spline_age_2023 <- ggplot(pred[pred$year == 2023,], aes(x = age, group = gender,
color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point(position = position_dodge(width = .5)) +
geom_line(position = position_dodge(width = .5)) +
geom_errorbar(width = .5, position = position_dodge(width = .5)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ylab("Mean life satisfaction") +
xlab("Age") +
labs(color = "Gender") +
ggtitle("2023")
# Combine into single plot
combined_spline_age <- (spline_age_2010|spline_age_2015|spline_age_2023) + plot_layout(guides = "collect")
combined_spline_age
ggsave("Plots/satis_by_age_splines.png", width = 6, height = 3)Does age in isolation explain away the explanandum?
# Fit model including age
# Categorical version here (but in the final model we will use the spline)
explanandum_age <- lm(satis ~ as.factor(year)*gender*as.factor(age),
data = combined)
# New hypothetical data: demographics like in 2010
# This is our counterfactual world in which the sample
# composition did not change with respects to age
# (and any other variable we will include in the model later on)
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)
# Generate predictions for combinations of gender and year
pred <- predictions(explanandum_age,
by = c("gender", "year"),
newdata = all2010)
# Plot the results
ggplot(pred, aes(x = year, group = gender, color = gender, y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point() +
geom_line() +
geom_errorbar(width = .5) +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
ylab("Mean life satisfaction") +
xlab("Survey year") +
labs(color = "Gender")
ggsave("Plots/explanandum_age_controlled.png", width = 4, height = 3)
# Calculate the gender gaps on the hypothetical data
# In which the demographics remained unchanged across the three survey waves
comps <- avg_comparisons(explanandum_age, variables = "gender", by = "year",
newdata = all2010)
print(comps)##
## year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0.225 0.0442 5.10 <0.001 21.5 0.139 0.312
## 2015 0.313 0.0554 5.65 <0.001 25.9 0.205 0.422
## 2023 0.494 0.0425 11.61 <0.001 101.1 0.410 0.577
##
## Term: gender
## Type: response
## Comparison: male - female# Compare gender gaps 2010 vs 2023
max_comp <- avg_comparisons(explanandum_age, variables = "gender", by = "year",
hypothesis = "b1 = b3",
newdata = all2010)## 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 Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.268 0.0613 -4.38 <0.001 16.4 -0.389 -0.148
##
## Type: response# Compare gender gaps 2010 vs 2015
print(avg_comparisons(explanandum_age, variables = "gender", by = "year",
hypothesis = "b1 = b2"),
newdata = all2010)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.0537 0.0607 -0.885 0.376 1.4 -0.173 0.0652
##
## Type: response# Compare gender gaps 2015 vs 2023
print(avg_comparisons(explanandum_age, variables = "gender", by = "year",
hypothesis = "b2 = b3"),
newdata = all2010)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b3 -0.201 0.0569 -3.53 <0.001 11.3 -0.313 -0.0895
##
## Type: response# Express differences in SDs
round(comps$estimate[1]/sd(combined$satis), 2)## [1] 0.27round(comps$estimate[2]/sd(combined$satis), 2)## [1] 0.37round(comps$estimate[3]/sd(combined$satis), 2)## [1] 0.58# Express widening gender gap in SDs
round(max_comp$estimate/sd(combined$satis), 2)## [1] -0.32# Also calculate widening of gender gap for 2015 and 2023 age distribution
# new data: demographics like in 2015
temp1 <- combined[combined$year == 2015,]
temp1$year <- 2010
temp2 <- combined[combined$year == 2015,]
temp2$year <- 2015
temp3 <- combined[combined$year == 2015,]
temp3$year <- 2023
all2015 <- rbind(temp1, temp2, temp3)
# max comparison for 2015 age distribution
max_comp <- avg_comparisons(explanandum_age, variables = "gender", by = "year",
hypothesis = "b1 = b3")
print(max_comp)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.255 0.0588 -4.33 <0.001 16.1 -0.37 -0.14
##
## Type: response# new data: demographics like in 2023
temp1 <- combined[combined$year == 2023,]
temp1$year <- 2010
temp2 <- combined[combined$year == 2023,]
temp2$year <- 2015
temp3 <- combined[combined$year == 2023,]
temp3$year <- 2023
all2023 <- rbind(temp1, temp2, temp3)
# max comparison for 2023 age distribution
max_comp <- avg_comparisons(explanandum_age, variables = "gender", by = "year",
hypothesis = "b1 = b3",
newdata = all2015)
print(max_comp)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.339 0.064 -5.3 <0.001 23.1 -0.465 -0.214
##
## Type: responseDoes age modify the widening of the gender gap?
# Use the spline model
explanandum_age_spline <- lm(satis ~ as.factor(year)*gender*bs(age, df = 3),
data = combined)
# Calculate the gender gaps for each survey year and each age
all_comps <- avg_comparisons(explanandum_age_spline, variables = "gender", by = c("year", "age"))
print(all_comps)##
## year age Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 12 0.4086 0.1143 3.575 < 0.001 11.5 0.1846 0.633
## 2010 13 0.3285 0.0638 5.151 < 0.001 21.9 0.2035 0.453
## 2010 14 0.2160 0.0633 3.410 < 0.001 10.6 0.0918 0.340
## 2010 15 0.1065 0.0698 1.526 0.12701 3.0 -0.0303 0.243
## 2010 16 0.0355 0.0922 0.386 0.69972 0.5 -0.1451 0.216
## 2010 17 0.0385 0.1101 0.350 0.72654 0.5 -0.1773 0.254
## 2010 18 0.1509 0.2421 0.623 0.53318 0.9 -0.3236 0.625
## 2015 12 0.5735 0.2280 2.516 0.01188 6.4 0.1267 1.020
## 2015 13 0.2618 0.0763 3.432 < 0.001 10.7 0.1123 0.411
## 2015 14 0.1904 0.0627 3.039 0.00237 8.7 0.0676 0.313
## 2015 15 0.2552 0.0566 4.512 < 0.001 17.3 0.1444 0.366
## 2015 16 0.3521 0.0715 4.922 < 0.001 20.2 0.2119 0.492
## 2015 17 0.3769 0.0945 3.990 < 0.001 13.9 0.1917 0.562
## 2015 18 0.2254 0.2594 0.869 0.38481 1.4 -0.2829 0.734
## 2023 12 0.4165 0.1425 2.923 0.00347 8.2 0.1372 0.696
## 2023 13 0.5262 0.0637 8.264 < 0.001 52.7 0.4014 0.651
## 2023 14 0.5611 0.0588 9.539 < 0.001 69.2 0.4458 0.676
## 2023 15 0.5345 0.0536 9.978 < 0.001 75.5 0.4295 0.640
## 2023 16 0.4593 0.0680 6.758 < 0.001 36.1 0.3261 0.593
## 2023 17 0.3486 0.0771 4.521 < 0.001 17.3 0.1975 0.500
## 2023 18 0.2154 0.1522 1.415 0.15695 2.7 -0.0829 0.514
##
## Term: gender
## Type: response
## Comparison: male - female# Calculate the widenings of the gender gap for each age group
# Plot the age trajectory of the widening of the gender gap
widening_per_year <- avg_comparisons(explanandum_age_spline, variables = "gender", by = c("year", "age"),
hypothesis = c("b15 - b1 = 0",
"b16 - b2 = 0",
"b17 - b3 = 0",
"b18 - b4 = 0",
"b19 - b5 = 0",
"b20 - b6 = 0",
"b21 - b7 = 0"))
widening_per_year$age <- 12:18
# Plot widening per year
ggplot(widening_per_year, aes(x = age, y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point() +
geom_line() +
geom_errorbar(width = .5) +
coord_cartesian(ylim = c(-.5, 1)) +
theme_classic() +
ylab("Widening of gender gap") +
xlab("Age")
ggsave("Plots/age_modification.png", width = 4, height = 3)A closer look at Schooltype
Schooltype and life satisfaction
# Schooltype as a categorical predictor
add_schooltype <- lm(satis ~ gender*as.factor(year)*as.factor(schooltype),
data = combined)
summary(add_schooltype)##
## Call:
## lm(formula = satis ~ gender * as.factor(year) * as.factor(schooltype),
## data = combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.07576 -0.43344 0.03343 0.56656 1.57483
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 3.84112 0.04596
## gendermale 0.20872 0.06499
## as.factor(year)2015 -0.13800 0.06504
## as.factor(year)2023 -0.41595 0.06647
## as.factor(schooltype)2 -0.09767 0.06234
## gendermale:as.factor(year)2015 0.16391 0.09163
## gendermale:as.factor(year)2023 0.33267 0.09269
## gendermale:as.factor(schooltype)2 0.01440 0.08881
## as.factor(year)2015:as.factor(schooltype)2 0.17028 0.08516
## as.factor(year)2023:as.factor(schooltype)2 0.10594 0.08540
## gendermale:as.factor(year)2015:as.factor(schooltype)2 -0.17578 0.12277
## gendermale:as.factor(year)2023:as.factor(schooltype)2 -0.10980 0.12065
## t value Pr(>|t|)
## (Intercept) 83.582 < 2e-16 ***
## gendermale 3.212 0.001329 **
## as.factor(year)2015 -2.122 0.033921 *
## as.factor(year)2023 -6.258 4.24e-10 ***
## as.factor(schooltype)2 -1.567 0.117279
## gendermale:as.factor(year)2015 1.789 0.073722 .
## gendermale:as.factor(year)2023 3.589 0.000335 ***
## gendermale:as.factor(schooltype)2 0.162 0.871238
## as.factor(year)2015:as.factor(schooltype)2 2.000 0.045602 *
## as.factor(year)2023:as.factor(schooltype)2 1.240 0.214858
## gendermale:as.factor(year)2015:as.factor(schooltype)2 -1.432 0.152252
## gendermale:as.factor(year)2023:as.factor(schooltype)2 -0.910 0.362830
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8234 on 4752 degrees of freedom
## Multiple R-squared: 0.06129, Adjusted R-squared: 0.05911
## F-statistic: 28.2 on 11 and 4752 DF, p-value: < 2.2e-16# Generate predictions for combinations of gender and year and age
pred <- predictions(add_schooltype,
by = c("gender", "year", "schooltype"))
# Plot the results
# Schooltype 1 (general secondary)
cat_schooltype_1 <- ggplot(pred[pred$schooltype == 1,], aes(x = year, group = gender,
color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point(position = position_dodge(width = .5)) +
geom_line(position = position_dodge(width = .5)) +
geom_errorbar(width = .5, position = position_dodge(width = .5)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("General secondary school") +
ylab("Mean life satisfaction") +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
xlab("Year") +
labs(color = "Gender")
# Schooltype 2 (grammar school)
cat_schooltype_2 <- ggplot(pred[pred$schooltype == 2,], aes(x = year, group = gender,
color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point(position = position_dodge(width = .5)) +
geom_line(position = position_dodge(width = .5)) +
geom_errorbar(width = .5, position = position_dodge(width = .5)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("Grammar school") +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
ylab("Mean life satisfaction") +
xlab("Year") +
labs(color = "Gender")
# Combine into single plot
combined_cat_schooltype <- (cat_schooltype_1|cat_schooltype_2) + plot_layout(guides = "collect")
combined_cat_schooltype
ggsave("Plots/satis_by_schooltype_categorical.png", width = 6, height = 3)Does schooltype explain away the explanandum?
# Fit linear model including both age and schooltype
explanandum_schooltype <- lm(satis ~ as.factor(year) + gender + as.factor(year):gender +
as.factor(age) + as.factor(year):as.factor(age) + gender:as.factor(age) + as.factor(year):gender:as.factor(age) +
as.factor(schooltype) + as.factor(year):as.factor(schooltype) + gender:as.factor(schooltype) + as.factor(year):gender:as.factor(schooltype),
data = combined)
# Generate predictions for combinations of gender and year
pred <- predictions(explanandum_schooltype,
by = c("gender", "year"),
newdata = all2010,
vcov = ~ unique_classroom)
# Plot the results
ggplot(pred, aes(x = year, group = gender, color = gender, y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point() +
geom_line() +
geom_errorbar(width = .5) +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
ylab("Mean life satisfaction") +
xlab("Survey year") +
labs(color = "Gender")
ggsave("Plots/explanandum_age_schooltype_controlled.png", width = 4, height = 3)
# Gender gap by year
# For the 2010 schooltype distribution
comps <- avg_comparisons(explanandum_schooltype, variables = "gender", by = "year",
newdata = all2010,
vcov = ~ unique_classroom)
print(comps)##
## year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0.222 0.0468 4.76 <0.001 19.0 0.131 0.314
## 2015 0.320 0.0490 6.52 <0.001 33.7 0.223 0.416
## 2023 0.501 0.0398 12.59 <0.001 118.4 0.423 0.579
##
## Term: gender
## Type: response
## Comparison: male - female# Test widening from 2010 to 2023
max_comp <- avg_comparisons(explanandum_schooltype, variables = "gender", by = "year",
hypothesis = "b1 = b3", newdata = all2010, vcov = ~ unique_classroom)
print(max_comp)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.279 0.0614 -4.54 <0.001 17.4 -0.399 -0.158
##
## Type: response# Gender gap 2010 vs 2015
print(avg_comparisons(explanandum_schooltype, variables = "gender", by = "year",
hypothesis = "b1 = b2", newdata = all2010, vcov = ~ unique_classroom))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.0972 0.0678 -1.44 0.151 2.7 -0.23 0.0356
##
## Type: response# Gender gap 2015 vs 2023
print(avg_comparisons(explanandum_schooltype, variables = "gender", by = "year",
hypothesis = "b2 = b3", newdata = all2010, vcov = ~ unique_classroom))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b3 -0.182 0.0632 -2.87 0.00405 7.9 -0.305 -0.0578
##
## Type: response# Express differences in SDs
round(comps$estimate[1]/sd(combined$satis), 2)## [1] 0.26round(comps$estimate[2]/sd(combined$satis), 2)## [1] 0.38round(comps$estimate[3]/sd(combined$satis), 2)## [1] 0.59# Express widening gender gap in SDs
round(max_comp$estimate/sd(combined$satis), 2)## [1] -0.33# Widening of the gap from 2010 to 2023 at 2015 schooltype distribution
max_comp <- avg_comparisons(explanandum_schooltype, variables = "gender", by = "year",
hypothesis = "b1 = b3", newdata = all2015, vcov = ~ unique_classroom)
print(max_comp)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.346 0.0703 -4.92 <0.001 20.2 -0.484 -0.208
##
## Type: response# Widening of the gap from 2010 to 2023 at 2023 schooltype distribution
max_comp <- avg_comparisons(explanandum_schooltype, variables = "gender", by = "year",
hypothesis = "b1 = b3", newdata = all2023, vcov = ~ unique_classroom)
print(max_comp)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.305 0.0667 -4.58 <0.001 17.7 -0.436 -0.174
##
## Type: responseDoes schooltype modify the widening of the gender gap?
# Gender gap for each year and each schooltype
all_comps <- avg_comparisons(explanandum_schooltype,
variables = "gender",
by = c("year", "schooltype"),
vcov = ~ unique_classroom)
# Test widening of gender gap for the lower track students
print(avg_comparisons(explanandum_schooltype,
variables = "gender",
by = c("year", "schooltype"),
vcov = ~ unique_classroom,
hypothesis = "b1 = b5")) ##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b5 -0.348 0.0923 -3.76 <0.001 12.5 -0.529 -0.167
##
## Type: response# Test widening of the gender gap for upper track students
print(avg_comparisons(explanandum_schooltype,
variables = "gender",
by = c("year", "schooltype"),
vcov = ~ unique_classroom,
hypothesis = "b2 = b6")) # upper track##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b6 -0.204 0.0772 -2.64 0.00829 6.9 -0.355 -0.0525
##
## Type: response# Compare widening of the gaps between upper and lower track
# Lower track versus upper track
print(avg_comparisons(explanandum_schooltype,
variables = "gender",
by = c("year", "schooltype"),
vcov = ~ unique_classroom,
hypothesis = "b1 - b5 = b2 - b6")) ##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1-b5=b2-b6 -0.144 0.12 -1.19 0.233 2.1 -0.38 0.0924
##
## Type: responseMigration background/German not main language
Migration background and life satisfaction
# Predict life satisfaction from migration background
add_mig_lang <- lm(satis ~ gender*as.factor(year)*as.factor(mig_lang),
data = combined)
summary(add_mig_lang)##
## Call:
## lm(formula = satis ~ gender * as.factor(year) * as.factor(mig_lang),
## data = combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0403 -0.4415 0.1058 0.5585 1.6435
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 3.77619 0.03280 115.133
## gendermale 0.24824 0.04669 5.317
## as.factor(year)2015 -0.01441 0.04403 -0.327
## as.factor(year)2023 -0.33468 0.04391 -7.622
## as.factor(mig_lang)1 0.11422 0.10178 1.122
## gendermale:as.factor(year)2015 0.03029 0.06398 0.473
## gendermale:as.factor(year)2023 0.20444 0.06250 3.271
## gendermale:as.factor(mig_lang)1 -0.30532 0.14742 -2.071
## as.factor(year)2015:as.factor(mig_lang)1 -0.25575 0.14072 -1.817
## as.factor(year)2023:as.factor(mig_lang)1 -0.19921 0.13079 -1.523
## gendermale:as.factor(year)2015:as.factor(mig_lang)1 0.31958 0.20505 1.559
## gendermale:as.factor(year)2023:as.factor(mig_lang)1 0.52153 0.18734 2.784
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## gendermale 1.1e-07 ***
## as.factor(year)2015 0.74355
## as.factor(year)2023 3.0e-14 ***
## as.factor(mig_lang)1 0.26183
## gendermale:as.factor(year)2015 0.63597
## gendermale:as.factor(year)2023 0.00108 **
## gendermale:as.factor(mig_lang)1 0.03840 *
## as.factor(year)2015:as.factor(mig_lang)1 0.06921 .
## as.factor(year)2023:as.factor(mig_lang)1 0.12779
## gendermale:as.factor(year)2015:as.factor(mig_lang)1 0.11917
## gendermale:as.factor(year)2023:as.factor(mig_lang)1 0.00539 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8232 on 4752 degrees of freedom
## Multiple R-squared: 0.0616, Adjusted R-squared: 0.05943
## F-statistic: 28.36 on 11 and 4752 DF, p-value: < 2.2e-16# Generate predictions for combinations of gender and year and age
pred <- predictions(add_mig_lang,
by = c("gender", "year", "mig_lang"))
# Plot the results
# No migration background
cat_mig_lang_0 <- ggplot(pred[pred$mig_lang == 0,], aes(x = year, group = gender,
color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point(position = position_dodge(width = .5)) +
geom_line(position = position_dodge(width = .5)) +
geom_errorbar(width = .5, position = position_dodge(width = .5)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("No migration background") +
ylab("Mean life satisfaction") +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
xlab("Year") +
labs(color = "Gender")
# With migration background
cat_mig_lang_1 <- ggplot(pred[pred$mig_lang == 1,], aes(x = year, group = gender,
color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point(position = position_dodge(width = .5)) +
geom_line(position = position_dodge(width = .5)) +
geom_errorbar(width = .5, position = position_dodge(width = .5)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("Migration background") +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
ylab("Mean life satisfaction") +
xlab("Year") +
labs(color = "Gender")
# Combine into single plot
combined_cat_mig_lang <- (cat_mig_lang_0|cat_mig_lang_1) + plot_layout(guides = "collect")
combined_cat_mig_lang
ggsave("Plots/satis_by_mig_lang_categorical.png", width = 6, height = 3)Does migration background explain away the explanandum?
# Full model including the previous controls
explanandum_mig_lang <- lm(satis ~ as.factor(year) + gender + as.factor(year):gender +
as.factor(age) + as.factor(year):as.factor(age) + gender:as.factor(age) + as.factor(year):gender:as.factor(age) +
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),
data = combined)
# Generate predictions for combinations of gender and year
# Given the 2010 distribution of covariates
pred <- predictions(explanandum_mig_lang,
by = c("gender", "year"),
newdata = all2010,
vcov = ~ unique_classroom)
# Plot the results
ggplot(pred, aes(x = year, group = gender, color = gender, y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point() +
geom_line() +
geom_errorbar(width = .5) +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
ylab("Mean life satisfaction") +
xlab("Survey year") +
labs(color = "Gender")
ggsave("Plots/explanandum_age_schooltype_migrant_controlled.png", width = 4, height = 3)
# Gender gaps by year
comps <- avg_comparisons(explanandum_mig_lang, variables = "gender", by = "year",
newdata = all2010,
vcov = ~ unique_classroom)
print(comps)##
## year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0.222 0.0465 4.78 <0.001 19.1 0.131 0.313
## 2015 0.317 0.0494 6.43 <0.001 32.9 0.221 0.414
## 2023 0.496 0.0391 12.69 <0.001 120.2 0.419 0.572
##
## Term: gender
## Type: response
## Comparison: male - female# Test gender gap 2010 vs 2023
max_comp <- avg_comparisons(explanandum_mig_lang, variables = "gender", by = "year",
hypothesis = "b1 = b3", newdata = all2010, vcov = ~ unique_classroom)
print(max_comp)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.273 0.0607 -4.5 <0.001 17.2 -0.393 -0.154
##
## Type: response# Test gender gap 2010 vs 2015
print(avg_comparisons(explanandum_mig_lang, variables = "gender", by = "year",
hypothesis = "b1 = b2", newdata = all2010, vcov = ~ unique_classroom))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.0952 0.0678 -1.4 0.16 2.6 -0.228 0.0377
##
## Type: response# Test gender gap 2015 vs 2023
print(avg_comparisons(explanandum_mig_lang, variables = "gender", by = "year",
hypothesis = "b2 = b3", newdata = all2010, vcov = ~ unique_classroom))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b3 -0.178 0.0629 -2.83 0.00462 7.8 -0.302 -0.0549
##
## Type: response# Express differences in SDs
round(comps$estimate[1]/sd(combined$satis), 2)## [1] 0.26round(comps$estimate[2]/sd(combined$satis), 2)## [1] 0.37round(comps$estimate[3]/sd(combined$satis), 2)## [1] 0.58# Express widening gender gap in SDs
round(max_comp$estimate/sd(combined$satis), 2)## [1] -0.32# Compare gender gap 2010 vs 2023, this time for the 2015 distribution of covariates
max_comp <- avg_comparisons(explanandum_mig_lang, variables = "gender", by = "year",
hypothesis = "b1 = b3", newdata = all2015, vcov = ~ unique_classroom)
print(max_comp)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.332 0.0693 -4.79 <0.001 19.2 -0.468 -0.196
##
## Type: response# Compare gender gap 2010 vs 2023, this time for the 2023 distribution of covariates
max_comp <- avg_comparisons(explanandum_mig_lang, variables = "gender", by = "year",
hypothesis = "b1 = b3", newdata = all2023, vcov = ~ unique_classroom)
print(max_comp)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.308 0.0656 -4.7 <0.001 18.5 -0.436 -0.179
##
## Type: responseDoes migration background modify the widening of the gender gap?
# Calculate gender gaps for each year and each level of migration background
all_comps <- avg_comparisons(explanandum_mig_lang,
variables = "gender",
by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom)
print(all_comps)##
## year mig_lang Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0 0.2542 0.0463 5.486 <0.001 24.5 0.163 0.345
## 2010 1 -0.0628 0.1813 -0.346 0.7292 0.5 -0.418 0.293
## 2015 0 0.3208 0.0508 6.310 <0.001 31.7 0.221 0.420
## 2015 1 0.2883 0.1394 2.067 0.0387 4.7 0.015 0.562
## 2023 0 0.4778 0.0424 11.274 <0.001 95.5 0.395 0.561
## 2023 1 0.6571 0.1274 5.156 <0.001 21.9 0.407 0.907
##
## Term: gender
## Type: response
## Comparison: male - female# Widening of the gap for no migration background
print(avg_comparisons(explanandum_mig_lang,
variables = "gender",
by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom,
hypothesis = "b1 = b5")) ##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b5 -0.224 0.0628 -3.56 <0.001 11.4 -0.347 -0.101
##
## Type: response# Widening of the gap for migration background
print(avg_comparisons(explanandum_mig_lang,
variables = "gender",
by = c("year", "mig_lang"),
vcov = ~ unique_classroom,
newdata = all2010,
hypothesis = "b2 = b6")) # migback##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b6 -0.72 0.222 -3.25 0.00116 9.8 -1.15 -0.286
##
## Type: response# Compare the widening of the gap between no migration background and migration background
print(avg_comparisons(explanandum_mig_lang,
variables = "gender",
by = c("year", "mig_lang"),
vcov = ~ unique_classroom,
newdata = all2010,
hypothesis = "b1 - b5 = b2 - b6")) ##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1-b5=b2-b6 0.496 0.23 2.16 0.0311 5.0 0.0451 0.948
##
## Type: response# Widening of the gap for no migration background, 2023 demographics
print(avg_comparisons(explanandum_mig_lang,
variables = "gender",
by = c("year", "mig_lang"),
newdata = all2023,
vcov = ~ unique_classroom,
hypothesis = "b1 = b5")) ##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b5 -0.244 0.0697 -3.5 <0.001 11.0 -0.38 -0.107
##
## Type: response# Widening of the gap for migration background, 2023 demographics
print(avg_comparisons(explanandum_mig_lang,
variables = "gender",
by = c("year", "mig_lang"),
vcov = ~ unique_classroom,
newdata = all2023,
hypothesis = "b2 = b6")) # migback##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b6 -0.741 0.219 -3.39 <0.001 10.5 -1.17 -0.312
##
## Type: responseFinal model
In the manuscript, we report a final model to address all questions at once.
explanandum_demo_final <- lm(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),
data = combined)
summary(explanandum_demo_final)##
## Call:
## lm(formula = 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), data = combined)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.03100 -0.42474 0.08081 0.51296 1.68600
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 3.896507 0.085906
## as.factor(year)2015 -0.117415 0.173293
## as.factor(year)2023 -0.209274 0.132052
## gendermale 0.438833 0.127242
## bs(age, df = 3)1 -0.180031 0.251854
## bs(age, df = 3)2 0.102056 0.256868
## bs(age, df = 3)3 -0.321316 0.209255
## as.factor(schooltype)2 -0.101549 0.064027
## as.factor(mig_lang)1 0.145433 0.102242
## as.factor(year)2015:gendermale 0.220952 0.265749
## as.factor(year)2023:gendermale 0.011550 0.200678
## as.factor(year)2015:bs(age, df = 3)1 0.327360 0.446301
## as.factor(year)2023:bs(age, df = 3)1 -0.352580 0.355244
## as.factor(year)2015:bs(age, df = 3)2 -0.436729 0.338988
## as.factor(year)2023:bs(age, df = 3)2 -0.302771 0.313682
## as.factor(year)2015:bs(age, df = 3)3 0.154757 0.322326
## as.factor(year)2023:bs(age, df = 3)3 0.036376 0.249425
## gendermale:bs(age, df = 3)1 -0.133051 0.362450
## gendermale:bs(age, df = 3)2 -0.597044 0.359734
## gendermale:bs(age, df = 3)3 -0.235878 0.280632
## as.factor(year)2015:as.factor(schooltype)2 0.199255 0.087116
## as.factor(year)2023:as.factor(schooltype)2 0.110627 0.087721
## gendermale:as.factor(schooltype)2 -0.005963 0.091168
## as.factor(year)2015:as.factor(mig_lang)1 -0.276221 0.140841
## as.factor(year)2023:as.factor(mig_lang)1 -0.217127 0.131052
## gendermale:as.factor(mig_lang)1 -0.273404 0.148016
## as.factor(year)2015:gendermale:bs(age, df = 3)1 -0.777231 0.681305
## as.factor(year)2023:gendermale:bs(age, df = 3)1 0.395215 0.528223
## as.factor(year)2015:gendermale:bs(age, df = 3)2 0.787005 0.492449
## as.factor(year)2023:gendermale:bs(age, df = 3)2 0.646881 0.446553
## as.factor(year)2015:gendermale:bs(age, df = 3)3 -0.005570 0.481223
## as.factor(year)2023:gendermale:bs(age, df = 3)3 0.037347 0.357913
## as.factor(year)2015:gendermale:as.factor(schooltype)2 -0.179717 0.125545
## as.factor(year)2023:gendermale:as.factor(schooltype)2 -0.061359 0.123751
## as.factor(year)2015:gendermale:as.factor(mig_lang)1 0.279035 0.205460
## as.factor(year)2023:gendermale:as.factor(mig_lang)1 0.484946 0.187673
## t value Pr(>|t|)
## (Intercept) 45.358 < 2e-16 ***
## as.factor(year)2015 -0.678 0.498089
## as.factor(year)2023 -1.585 0.113081
## gendermale 3.449 0.000568 ***
## bs(age, df = 3)1 -0.715 0.474755
## bs(age, df = 3)2 0.397 0.691157
## bs(age, df = 3)3 -1.536 0.124723
## as.factor(schooltype)2 -1.586 0.112796
## as.factor(mig_lang)1 1.422 0.154964
## as.factor(year)2015:gendermale 0.831 0.405773
## as.factor(year)2023:gendermale 0.058 0.954105
## as.factor(year)2015:bs(age, df = 3)1 0.733 0.463293
## as.factor(year)2023:bs(age, df = 3)1 -0.993 0.321003
## as.factor(year)2015:bs(age, df = 3)2 -1.288 0.197694
## as.factor(year)2023:bs(age, df = 3)2 -0.965 0.334485
## as.factor(year)2015:bs(age, df = 3)3 0.480 0.631162
## as.factor(year)2023:bs(age, df = 3)3 0.146 0.884054
## gendermale:bs(age, df = 3)1 -0.367 0.713569
## gendermale:bs(age, df = 3)2 -1.660 0.097045 .
## gendermale:bs(age, df = 3)3 -0.841 0.400657
## as.factor(year)2015:as.factor(schooltype)2 2.287 0.022227 *
## as.factor(year)2023:as.factor(schooltype)2 1.261 0.207327
## gendermale:as.factor(schooltype)2 -0.065 0.947857
## as.factor(year)2015:as.factor(mig_lang)1 -1.961 0.049912 *
## as.factor(year)2023:as.factor(mig_lang)1 -1.657 0.097625 .
## gendermale:as.factor(mig_lang)1 -1.847 0.064793 .
## as.factor(year)2015:gendermale:bs(age, df = 3)1 -1.141 0.254012
## as.factor(year)2023:gendermale:bs(age, df = 3)1 0.748 0.454379
## as.factor(year)2015:gendermale:bs(age, df = 3)2 1.598 0.110078
## as.factor(year)2023:gendermale:bs(age, df = 3)2 1.449 0.147513
## as.factor(year)2015:gendermale:bs(age, df = 3)3 -0.012 0.990765
## as.factor(year)2023:gendermale:bs(age, df = 3)3 0.104 0.916898
## as.factor(year)2015:gendermale:as.factor(schooltype)2 -1.431 0.152355
## as.factor(year)2023:gendermale:as.factor(schooltype)2 -0.496 0.620043
## as.factor(year)2015:gendermale:as.factor(mig_lang)1 1.358 0.174498
## as.factor(year)2023:gendermale:as.factor(mig_lang)1 2.584 0.009796 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8189 on 4728 degrees of freedom
## Multiple R-squared: 0.07611, Adjusted R-squared: 0.06927
## F-statistic: 11.13 on 35 and 4728 DF, p-value: < 2.2e-16Does the model explain away the explanandum?
# Look at gender gaps for 2010 demographics, by year
comps <- avg_comparisons(explanandum_demo_final, variables = "gender", by = "year",
newdata = all2010,
vcov = ~ unique_classroom)
print(comps)##
## year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0.223 0.0463 4.83 <0.001 19.5 0.133 0.314
## 2015 0.309 0.0493 6.27 <0.001 31.3 0.212 0.405
## 2023 0.495 0.0390 12.69 <0.001 120.2 0.419 0.572
##
## Term: gender
## Type: response
## Comparison: male - female# Test widening from 2010 to 2023
max_comp <- avg_comparisons(explanandum_demo_final, variables = "gender", by = "year",
hypothesis = "b1 = b3", newdata = all2010, vcov = ~ unique_classroom)
print(max_comp)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.272 0.0605 -4.49 <0.001 17.1 -0.391 -0.153
##
## Type: response# Express widening gender gap in SDs
round(max_comp$estimate/sd(combined$satis), 2)## [1] -0.32# Widening of the gap from 2010 to 2023 at 2015 distribution
max_comp <- avg_comparisons(explanandum_demo_final, variables = "gender", by = "year",
hypothesis = "b1 = b3", newdata = all2015, vcov = ~ unique_classroom)
print(max_comp)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.329 0.0697 -4.72 <0.001 18.7 -0.466 -0.193
##
## Type: response# Express widening gender gap in SDs
round(max_comp$estimate/sd(combined$satis), 2)## [1] -0.39# Widening of the gap from 2010 to 2023 at 2023 distribution
max_comp <- avg_comparisons(explanandum_demo_final, variables = "gender", by = "year",
hypothesis = "b1 = b3", newdata = all2023, vcov = ~ unique_classroom)
print(max_comp)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.312 0.0661 -4.72 <0.001 18.7 -0.442 -0.182
##
## Type: response# Express widening gender gap in SDs
round(max_comp$estimate/sd(combined$satis), 2)## [1] -0.37Do demographics modify the widening of the gap?
Modification by age
This will give us Figure 5 from the manuscript.
# Calculate margins
all_age_comps <- avg_comparisons(explanandum_demo_final, variables = "gender",
by = c("age", "year"),
vcov = ~ unique_classroom,
newdata = all2010)
print(all_age_comps)##
## age year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 12 2010 0.4135 0.0973 4.252 < 0.001 15.5 0.2229 0.604
## 12 2015 0.5516 0.1740 3.170 0.00152 9.4 0.2106 0.893
## 12 2023 0.4279 0.1368 3.128 0.00176 9.1 0.1598 0.696
## 13 2010 0.3262 0.0649 5.026 < 0.001 20.9 0.1990 0.453
## 13 2015 0.2394 0.0722 3.314 < 0.001 10.1 0.0978 0.381
## 13 2023 0.5172 0.0615 8.412 < 0.001 54.5 0.3967 0.638
## 14 2010 0.2058 0.0616 3.340 < 0.001 10.2 0.0850 0.327
## 14 2015 0.1953 0.0552 3.537 < 0.001 11.3 0.0871 0.304
## 14 2023 0.5594 0.0598 9.356 < 0.001 66.7 0.4422 0.677
## 15 2010 0.1123 0.0891 1.260 0.20756 2.3 -0.0623 0.287
## 15 2015 0.2932 0.0570 5.148 < 0.001 21.9 0.1816 0.405
## 15 2023 0.5347 0.0574 9.316 < 0.001 66.2 0.4222 0.647
## 16 2010 0.0285 0.1130 0.252 0.80100 0.3 -0.1930 0.250
## 16 2015 0.4034 0.0794 5.083 < 0.001 21.4 0.2479 0.559
## 16 2023 0.4809 0.0767 6.270 < 0.001 31.4 0.3306 0.631
## 17 2010 0.0325 0.1270 0.256 0.79770 0.3 -0.2163 0.281
## 17 2015 0.3966 0.0793 5.002 < 0.001 20.8 0.2412 0.552
## 17 2023 0.3629 0.0738 4.919 < 0.001 20.1 0.2183 0.507
## 18 2010 0.1535 0.3635 0.422 0.67292 0.6 -0.5590 0.866
## 18 2015 0.2587 0.1731 1.494 0.13510 2.9 -0.0806 0.598
## 18 2023 0.2279 0.1393 1.636 0.10178 3.3 -0.0451 0.501
##
## Term: gender
## Type: response
## Comparison: male - female# Calculate widening of the gender gap separately per year of age
widening_per_year <- avg_comparisons(explanandum_age_spline, variables = "gender", by = c("age", "year"),
hypothesis = c("b3 - b1 = 0",
"b6 - b4 = 0",
"b9 - b7 = 0",
"b12 - b10 = 0",
"b15 - b13 = 0",
"b18 - b16 = 0",
"b21 - b19 = 0"),
vcov = ~ unique_classroom,
newdata = all2010)
widening_per_year$age <- 12:18
# Plot widening of the gender gap
ggplot(widening_per_year, aes(x = age, y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point() +
geom_line() +
geom_errorbar(width = .5) +
coord_cartesian(ylim = c(-.5, 1)) +
theme_classic() +
ylab("Widening of gender gap") +
xlab("Age")
ggsave("Plots/age_modification_final.png", width = 4, height = 3)
# Test maximum contrast
avg_comparisons(explanandum_age_spline, variables = "gender", by = c("age", "year"),
hypothesis = "b3 - b1 = b12 - b10",
vcov = ~ unique_classroom,
newdata = all2010)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b3-b1=b12-b10 -0.42 0.203 -2.07 0.0385 4.7 -0.818 -0.0223
##
## Type: responseModification by schooltype
# Calculate margins
all_schooltype_comps <- avg_comparisons(explanandum_demo_final, variables = "gender",
vcov = ~ unique_classroom,
by = c("schooltype", "year"),
newdata = all2010)
print(all_schooltype_comps)##
## schooltype year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 1 2010 0.214 0.0671 3.19 0.00141 9.5 0.0827 0.346
## 1 2015 0.403 0.0696 5.78 < 0.001 27.0 0.2663 0.539
## 1 2023 0.535 0.0623 8.59 < 0.001 56.7 0.4130 0.657
## 2 2010 0.231 0.0644 3.59 < 0.001 11.6 0.1052 0.358
## 2 2015 0.227 0.0595 3.82 < 0.001 12.9 0.1108 0.344
## 2 2023 0.461 0.0484 9.52 < 0.001 69.0 0.3662 0.556
##
## Term: gender
## Type: response
## Comparison: male - female# Widening among schooltype 1
avg_comparisons(explanandum_demo_final, variables = "gender",
by = c("schooltype", "year"),
vcov = ~ unique_classroom,
hypothesis = "b1 = b3", newdata = all2010)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.321 0.0915 -3.51 <0.001 11.1 -0.5 -0.141
##
## Type: response# Widening among schooltype 2
avg_comparisons(explanandum_demo_final, variables = "gender",
by = c("schooltype", "year"),
vcov = ~ unique_classroom,
hypothesis = "b4 = b6", newdata = all2010)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b4=b6 -0.23 0.0806 -2.85 0.00438 7.8 -0.388 -0.0717
##
## Type: response# Interaction
avg_comparisons(explanandum_demo_final, variables = "gender",
by = c("schooltype", "year"),
vcov = ~ unique_classroom,
hypothesis = "b1 - b3 = b4 - b6", newdata = all2010)##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1-b3=b4-b6 -0.0912 0.122 -0.747 0.455 1.1 -0.33 0.148
##
## Type: responseModification by migration background
# Calculate margins
all_mig_lang_comps <- avg_comparisons(explanandum_demo_final, variables = "gender",
vcov = ~ unique_classroom,
newdata = all2010,
by = c("mig_lang", "year"))
print(all_mig_lang_comps)##
## mig_lang year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 2010 0.2542 0.0465 5.471 <0.001 24.4 0.1631 0.345
## 0 2015 0.3095 0.0505 6.123 <0.001 30.0 0.2104 0.409
## 0 2023 0.4766 0.0423 11.262 <0.001 95.3 0.3936 0.560
## 1 2010 -0.0519 0.1769 -0.294 0.7690 0.4 -0.3986 0.295
## 1 2015 0.3025 0.1397 2.166 0.0303 5.0 0.0287 0.576
## 1 2023 0.6643 0.1258 5.283 <0.001 22.9 0.4178 0.911
##
## Term: gender
## Type: response
## Comparison: male - female# Widening among mig_lang 0
avg_comparisons(explanandum_demo_final, variables = "gender",
by = c("mig_lang", "year"),
vcov = ~ unique_classroom,
newdata = all2010,
hypothesis = "b1 = b3")##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b3 -0.222 0.0629 -3.54 <0.001 11.3 -0.346 -0.0992
##
## Type: response# Widening among mig_lang 1
avg_comparisons(explanandum_demo_final, variables = "gender",
by = c("mig_lang", "year"),
vcov = ~ unique_classroom,
newdata = all2010,
hypothesis = "b4 = b6")##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b4=b6 -0.716 0.217 -3.3 <0.001 10.0 -1.14 -0.291
##
## Type: response# Interaction
avg_comparisons(explanandum_demo_final, variables = "gender",
by = c("mig_lang", "year"),
vcov = ~ unique_classroom,
newdata = all2010,
hypothesis = "b1 - b3 = b4 - b6")##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1-b3=b4-b6 0.494 0.226 2.18 0.0289 5.1 0.0507 0.937
##
## Type: responseAdjusted widening of gender gap, separated by migration background
Generate Figure 6 which summarizes our preliminary results.
# Generate predictions for combinations of gender and year and age
pred <- predictions(explanandum_demo_final,
vcov = ~ unique_classroom,
newdata = all2010,
by = c("gender", "year", "mig_lang"))
# Duplicate predictions for plotting purposes
pred1 <- pred
pred2 <- pred
# Pred 1: color
pred1$color <- 1
pred1 <- pred1[pred1$year != 2015,]
# Pred 2: grey
pred2$color <- 0
# Manually dodge the years for plotting purposes
pred1$year <- pred1$year + 0.5 * (pred1$gender == "male")
pred2$year <- pred2$year + 0.5 * (pred2$gender == "male")
# Plot the results
# No migration background
cat_mig_lang_0 <- ggplot(pred1[pred1$mig_lang == 0,],
aes(x = year, group = gender, color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
# Parts in grey
geom_point(data = pred2[pred2$mig_lang == 0,],
aes(x = year, group = gender,
y = estimate), color = "lightgrey") +
geom_line(data = pred2[pred2$mig_lang == 0,],
aes(x = year, group = gender,
y = estimate), color = "lightgrey") +
geom_errorbar(data = pred2[pred2$mig_lang == 0,],
aes(x = year, group = gender,
y = estimate, ymin = conf.low, ymax = conf.high),
width = .5, color = "lightgrey") +
# Parts in color
geom_point() +
geom_line() +
geom_errorbar(width = .5) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("No migration background") +
ylab("Mean life satisfaction") +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
xlab("Year") +
labs(color = "Gender")
# With migration background
cat_mig_lang_1 <- ggplot(pred1[pred1$mig_lang == 1,],
aes(x = year, group = gender, color = gender,
y = estimate, ymin = conf.low, ymax = conf.high)) +
# Parts in grey
geom_point(data = pred2[pred2$mig_lang == 1,],
aes(x = year, group = gender,
y = estimate), color = "lightgrey") +
geom_line(data = pred2[pred2$mig_lang == 1,],
aes(x = year, group = gender,
y = estimate), color = "lightgrey") +
geom_errorbar(data = pred2[pred2$mig_lang == 1,],
aes(x = year, group = gender,
y = estimate, ymin = conf.low, ymax = conf.high),
width = .5, color = "lightgrey") +
# Parts in color
geom_point() +
geom_line() +
geom_errorbar(width = .5) +
coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)) +
ggtitle("Migration background") +
ylab("Mean life satisfaction") +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
xlab("Year") +
labs(color = "Gender")
# Combine into single plot
combined_cat_mig_lang <- (cat_mig_lang_0|cat_mig_lang_1) + plot_layout(guides = "collect")
combined_cat_mig_lang
ggsave("Plots/demo_final_by_mig.png", width = 7, height = 3)