Check our explanation for different outcomes
Julia Rohrer
2026-03-13
- Check
other outcomes
- Satisfaction with the amount of money with which you have to get by (financial satisfaction)
- Satisfaction with your relationships to friends
- Satisfaction with your relationship to your mother
- Satisfaction with your relationship to your father
- Satisfaction with your leisure time activities
- Satisfaction with your school grades
- Satisfaction with your housing situation
Check other outcomes
Let’s see whether our favored explanation also holds up for the other satisfaction variables. Here, we will also generate all panels of Figure 9.
Satisfaction with the amount of money with which you have to get by (financial satisfaction)
# Fit the model
explanandum_paper_only_money <- lm(satis_money ~ 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[combined$source == "paper",])
# 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)
# Predictions
pred <- predictions(explanandum_paper_only_money,
by = c("gender", "year"),
newdata = all2010,
vcov = ~ unique_classroom)
print(pred)##
## gender year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 2010 3.76 0.0466 80.7 <0.001 Inf 3.67 3.85
## female 2015 3.82 0.0368 103.8 <0.001 Inf 3.75 3.89
## female 2023 3.98 0.0665 59.8 <0.001 Inf 3.85 4.11
## male 2010 3.91 0.0455 85.9 <0.001 Inf 3.82 4.00
## male 2015 3.99 0.0503 79.4 <0.001 Inf 3.89 4.09
## male 2023 4.01 0.1011 39.6 <0.001 Inf 3.81 4.21
##
## Type: response# Gender gap by year
comps <- avg_comparisons(explanandum_paper_only_money, 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.1584 0.0547 2.899 0.00375 8.1 0.0513 0.266
## 2015 0.1782 0.0626 2.847 0.00441 7.8 0.0555 0.301
## 2023 0.0483 0.1261 0.383 0.70175 0.5 -0.1989 0.296
##
## Term: gender
## Type: response
## Comparison: male - female# Test biggest contrast: 2010 versus 2023
max_comp <- avg_comparisons(explanandum_paper_only_money, variables = "gender", by = "year",
hypothesis = "b1 = b3",
newdata = all2010,
vcov = ~ unique_classroom)## 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.11 0.137 0.801 0.423 1.2 -0.159 0.38
##
## Type: responseprint(max_comp$estimate/sd(combined$satis_money, na.rm = TRUE))## [1] 0.1126688# Gender gap by year and Mig back
comps_by_migback <- avg_comparisons(explanandum_paper_only_money, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom)
print(comps_by_migback)##
## year mig_lang Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0 0.1707 0.0581 2.938 0.00331 8.2 0.0568 0.285
## 2010 1 0.0490 0.2133 0.230 0.81832 0.3 -0.3691 0.467
## 2015 0 0.1637 0.0661 2.478 0.01323 6.2 0.0342 0.293
## 2015 1 0.3078 0.1322 2.328 0.01990 5.7 0.0487 0.567
## 2023 0 0.0164 0.1380 0.119 0.90538 0.1 -0.2540 0.287
## 2023 1 0.3338 0.2188 1.526 0.12699 3.0 -0.0949 0.763
##
## Term: gender
## Type: response
## Comparison: male - female# Widening of the gap for no Mig back
# 2023
print(avg_comparisons(explanandum_paper_only_money, 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.154 0.15 1.03 0.303 1.7 -0.139 0.448
##
## Type: response# Widening of the gap for Mig back
# 2023
print(avg_comparisons(explanandum_paper_only_money, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom,
hypothesis = "b2 = b6"))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b6 -0.285 0.306 -0.932 0.351 1.5 -0.884 0.314
##
## Type: responsepred <- avg_predictions(explanandum_paper_only_money,
newdata = all2010,
by = c("gender", "year", "mig_lang"),
vcov = ~ unique_classroom)
# Plot the results
# 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 mig back
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 mig back") +
ylab("Mean financial satisfaction") +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
xlab("Year") +
labs(color = "Gender")
# With Mig back
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("Mig back") +
ylab("Mean financial 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/paper_satis_money.png", width = width_plot, height = height_plot)
# Separate plots for manuscript
ggsave("Plots/paper_satis_money_no.png",
plot = cat_mig_lang_0 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)
ggsave("Plots/paper_satis_money_yes.png",
plot = cat_mig_lang_1 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)Satisfaction with your relationships to friends
# Fit the model
explanandum_paper_only_friends <- lm(satis_friends ~ 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[combined$source == "paper",])
# Predictions
pred <- predictions(explanandum_paper_only_friends,
by = c("gender", "year"),
newdata = all2010,
vcov = ~ unique_classroom)
print(pred)##
## gender year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 2010 4.39 0.0306 143.6 <0.001 Inf 4.33 4.45
## female 2015 4.33 0.0339 127.7 <0.001 Inf 4.26 4.40
## female 2023 4.15 0.0699 59.4 <0.001 Inf 4.01 4.29
## male 2010 4.39 0.0368 119.0 <0.001 Inf 4.31 4.46
## male 2015 4.40 0.0318 138.5 <0.001 Inf 4.34 4.47
## male 2023 4.43 0.0779 56.8 <0.001 Inf 4.27 4.58
##
## Type: response# Gender gap by year
comps <- avg_comparisons(explanandum_paper_only_friends, 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.00367 0.0425 -0.0863 0.9312 0.1 -0.0870 0.0797
## 2015 0.07611 0.0455 1.6714 0.0946 3.4 -0.0131 0.1654
## 2023 0.28095 0.1086 2.5865 0.0097 6.7 0.0681 0.4938
##
## Term: gender
## Type: response
## Comparison: male - female# Test biggest contrast: 2010 versus 2023
max_comp <- avg_comparisons(explanandum_paper_only_friends, 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.285 0.117 -2.44 0.0147 6.1 -0.513 -0.056
##
## Type: responseprint(max_comp$estimate/sd(combined$satis_friends, na.rm = TRUE))## [1] -0.3519398# Gender gap by year and Mig back
comps_by_migback <- avg_comparisons(explanandum_paper_only_friends, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom)
print(comps_by_migback)##
## year mig_lang Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0 -0.0201 0.0448 -0.448 0.6541 0.6 -0.1078 0.0677
## 2010 1 0.1430 0.1568 0.912 0.3616 1.5 -0.1643 0.4504
## 2015 0 0.0774 0.0476 1.626 0.1040 3.3 -0.0159 0.1706
## 2015 1 0.0649 0.1357 0.478 0.6323 0.7 -0.2011 0.3310
## 2023 0 0.2715 0.1214 2.237 0.0253 5.3 0.0336 0.5094
## 2023 1 0.3655 0.2924 1.250 0.2114 2.2 -0.2077 0.9386
##
## Term: gender
## Type: response
## Comparison: male - female# Widening of the gap for no Mig back
print(avg_comparisons(explanandum_paper_only_friends, 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.292 0.129 -2.25 0.0242 5.4 -0.545 -0.038
##
## Type: response# Widening of the gap for Mig back
print(avg_comparisons(explanandum_paper_only_friends, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom,
hypothesis = "b2 = b6"))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b6 -0.222 0.332 -0.67 0.503 1.0 -0.873 0.428
##
## Type: response# Compare widenings
print(avg_comparisons(explanandum_paper_only_friends, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
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.0691 0.372 -0.186 0.853 0.2 -0.798 0.66
##
## Type: responsepred <- avg_predictions(explanandum_paper_only_friends,
newdata = all2010,
by = c("gender", "year", "mig_lang"),
vcov = ~ unique_classroom)
# Plot the results
# 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 mig back
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 mig back") +
ylab("Mean friendship satisfaction") +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
xlab("Year") +
labs(color = "Gender")
# With Mig back
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("Mig back") +
ylab("Mean friendship 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/paper_satis_friends.png", width = width_plot, height = height_plot)
# Separate plots for manuscript
ggsave("Plots/paper_satis_friends_no.png",
plot = cat_mig_lang_0 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)
ggsave("Plots/paper_satis_friends_yes.png",
plot = cat_mig_lang_1 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)Satisfaction with your relationship to your mother
# Fit the model
explanandum_paper_only_mom <- lm(satis_mom ~ 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[combined$source == "paper",])
# Predictions
pred <- predictions(explanandum_paper_only_mom,
by = c("gender", "year"),
newdata = all2010,
vcov = ~ unique_classroom)
print(pred)##
## gender year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 2010 4.28 0.0418 102.3 <0.001 Inf 4.20 4.36
## female 2015 4.27 0.0444 96.3 <0.001 Inf 4.19 4.36
## female 2023 4.33 0.0566 76.5 <0.001 Inf 4.22 4.44
## male 2010 4.42 0.0362 122.0 <0.001 Inf 4.34 4.49
## male 2015 4.40 0.0408 107.9 <0.001 Inf 4.32 4.48
## male 2023 4.52 0.0508 89.0 <0.001 Inf 4.42 4.62
##
## Type: response# Gender gap by year
comps <- avg_comparisons(explanandum_paper_only_mom, 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.141 0.0536 2.64 0.00838 6.9 0.0363 0.246
## 2015 0.135 0.0606 2.23 0.02545 5.3 0.0166 0.254
## 2023 0.196 0.0833 2.35 0.01858 5.8 0.0328 0.359
##
## Term: gender
## Type: response
## Comparison: male - female# Test biggest contrast: 2010 versus 2023
max_comp <- avg_comparisons(explanandum_paper_only_mom, 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.0546 0.0991 -0.552 0.581 0.8 -0.249 0.139
##
## Type: responseprint(max_comp$estimate/sd(combined$satis_mom, na.rm = TRUE))## [1] -0.05792211# Gender gap by year and Mig back
comps_by_migback <- avg_comparisons(explanandum_paper_only_mom, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom)
print(comps_by_migback)##
## year mig_lang Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0 0.15774 0.0547 2.885 0.00391 8.0 0.0506 0.265
## 2010 1 -0.00499 0.1560 -0.032 0.97446 0.0 -0.3107 0.301
## 2015 0 0.15803 0.0633 2.496 0.01256 6.3 0.0339 0.282
## 2015 1 -0.06783 0.1259 -0.539 0.59010 0.8 -0.3146 0.179
## 2023 0 0.14355 0.0915 1.568 0.11680 3.1 -0.0358 0.323
## 2023 1 0.66567 0.2936 2.267 0.02336 5.4 0.0903 1.241
##
## Term: gender
## Type: response
## Comparison: male - female# Widening of the gap for no Mig back
print(avg_comparisons(explanandum_paper_only_mom, 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.0142 0.107 0.133 0.894 0.2 -0.195 0.223
##
## Type: response# Widening of the gap for Mig back
print(avg_comparisons(explanandum_paper_only_mom, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom,
hypothesis = "b2 = b6"))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b6 -0.671 0.332 -2.02 0.0437 4.5 -1.32 -0.0191
##
## Type: response# Compare widenings
print(avg_comparisons(explanandum_paper_only_mom, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
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.685 0.357 1.92 0.0551 4.2 -0.0148 1.38
##
## Type: responsepred <- avg_predictions(explanandum_paper_only_mom,
newdata = all2010,
by = c("gender", "year", "mig_lang"),
vcov = ~ unique_classroom)
# Plot the results
# 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 mig back
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 mig back") +
ylab("Mean satisfaction w mom") +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
xlab("Year") +
labs(color = "Gender")
# With Mig back
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("Mig back") +
ylab("Mean satisfaction w mom") +
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/paper_satis_mother.png", width = width_plot, height = height_plot)
# Separate plots for manuscript
ggsave("Plots/paper_satis_mother_no.png",
plot = cat_mig_lang_0 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)
ggsave("Plots/paper_satis_mother_yes.png",
plot = cat_mig_lang_1 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)Satisfaction with your relationship to your father
# Fit the model
explanandum_paper_only_dad <- lm(satis_dad ~ 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[combined$source == "paper",])
# Predictions
pred <- predictions(explanandum_paper_only_dad,
by = c("gender", "year"),
newdata = all2010,
vcov = ~ unique_classroom)
print(pred)##
## gender year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 2010 3.82 0.0594 64.4 <0.001 Inf 3.70 3.94
## female 2015 3.99 0.0531 75.0 <0.001 Inf 3.88 4.09
## female 2023 3.79 0.0770 49.2 <0.001 Inf 3.64 3.94
## male 2010 4.16 0.0546 76.1 <0.001 Inf 4.05 4.26
## male 2015 4.28 0.0410 104.3 <0.001 Inf 4.20 4.36
## male 2023 4.10 0.0846 48.5 <0.001 Inf 3.94 4.27
##
## Type: response# Gender gap by year
comps <- avg_comparisons(explanandum_paper_only_dad, 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.344 0.0674 5.10 <0.001 21.5 0.2116 0.476
## 2015 0.303 0.0614 4.94 <0.001 20.3 0.1830 0.424
## 2023 0.330 0.1189 2.77 0.0056 7.5 0.0964 0.563
##
## Term: gender
## Type: response
## Comparison: male - female# Test biggest contrast: 2010 versus 2023
max_comp <- avg_comparisons(explanandum_paper_only_dad, 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.0142 0.137 0.104 0.917 0.1 -0.254 0.282
##
## Type: responseprint(max_comp$estimate/sd(combined$satis_dad, na.rm = TRUE))## [1] 0.01193285# Gender gap by year and Mig back
comps_by_migback <- avg_comparisons(explanandum_paper_only_dad, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom)
print(comps_by_migback)##
## year mig_lang Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0 0.347 0.0743 4.68 <0.001 18.4 0.2017 0.493
## 2010 1 0.312 0.1786 1.74 0.0811 3.6 -0.0385 0.662
## 2015 0 0.302 0.0661 4.57 <0.001 17.7 0.1727 0.432
## 2015 1 0.314 0.1741 1.81 0.0709 3.8 -0.0268 0.656
## 2023 0 0.250 0.1381 1.81 0.0706 3.8 -0.0210 0.520
## 2023 1 1.045 0.3027 3.45 <0.001 10.8 0.4513 1.638
##
## Term: gender
## Type: response
## Comparison: male - female# Widening of the gap for no Mig back
print(avg_comparisons(explanandum_paper_only_dad, 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.0977 0.157 0.623 0.533 0.9 -0.21 0.405
##
## Type: response# Widening of the gap for Mig back
print(avg_comparisons(explanandum_paper_only_dad, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom,
hypothesis = "b2 = b6"))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b6 -0.733 0.352 -2.09 0.037 4.8 -1.42 -0.0442
##
## Type: response# Compare widenings
print(avg_comparisons(explanandum_paper_only_dad, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
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.831 0.419 1.99 0.0471 4.4 0.0105 1.65
##
## Type: responsepred <- avg_predictions(explanandum_paper_only_dad,
newdata = all2010,
by = c("gender", "year", "mig_lang"),
vcov = ~ unique_classroom)
# Plot the results
# 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 mig back
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 mig back") +
ylab("Mean satisfaction w dad") +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
xlab("Year") +
labs(color = "Gender")
# With Mig back
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("Mig back") +
ylab("Mean satisfaction w dad") +
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/paper_satis_dad.png", width = width_plot, height = height_plot)
# Separate plots for manuscript
ggsave("Plots/paper_satis_dad_no.png",
plot = cat_mig_lang_0 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)
ggsave("Plots/paper_satis_dad_yes.png",
plot = cat_mig_lang_1 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)Satisfaction with your leisure time activities
# Fit the model
explanandum_paper_only_leisure <- lm(satis_leisure ~ 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[combined$source == "paper",])
# Predictions
pred <- predictions(explanandum_paper_only_leisure,
by = c("gender", "year"),
newdata = all2010,
vcov = ~ unique_classroom)
print(pred)##
## gender year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 2010 3.87 0.0458 84.5 <0.001 Inf 3.78 3.96
## female 2015 3.89 0.0439 88.7 <0.001 Inf 3.81 3.98
## female 2023 3.64 0.0526 69.2 <0.001 Inf 3.54 3.74
## male 2010 3.99 0.0337 118.2 <0.001 Inf 3.92 4.05
## male 2015 4.09 0.0448 91.2 <0.001 Inf 4.00 4.17
## male 2023 3.97 0.0389 102.2 <0.001 Inf 3.90 4.05
##
## Type: response# Gender gap by year
comps <- avg_comparisons(explanandum_paper_only_leisure, 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.135 0.0514 2.62 0.00888 6.8 0.0338 0.235
## 2015 0.197 0.0583 3.38 < 0.001 10.4 0.0825 0.311
## 2023 0.336 0.0541 6.21 < 0.001 30.8 0.2297 0.442
##
## Term: gender
## Type: response
## Comparison: male - female# Test biggest contrast: 2010 versus 2023
max_comp <- avg_comparisons(explanandum_paper_only_leisure, 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.201 0.0746 -2.69 0.00706 7.1 -0.347 -0.0548
##
## Type: responseprint(max_comp$estimate/sd(combined$satis_leisure, na.rm = TRUE))## [1] -0.2148885# Gender gap by year and Mig back
comps_by_migback <- avg_comparisons(explanandum_paper_only_leisure, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom)
print(comps_by_migback)##
## year mig_lang Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0 0.130 0.0592 2.19 0.0285 5.1 0.0136 0.246
## 2010 1 0.179 0.1346 1.33 0.1832 2.4 -0.0847 0.443
## 2015 0 0.176 0.0603 2.91 0.0036 8.1 0.0574 0.294
## 2015 1 0.385 0.1514 2.54 0.0110 6.5 0.0883 0.682
## 2023 0 0.305 0.0627 4.86 <0.001 19.7 0.1818 0.427
## 2023 1 0.613 0.2527 2.43 0.0152 6.0 0.1180 1.109
##
## Term: gender
## Type: response
## Comparison: male - female# Widening of the gap for no Mig back
print(avg_comparisons(explanandum_paper_only_leisure, 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.175 0.0862 -2.03 0.0423 4.6 -0.344 -0.00605
##
## Type: response# Widening of the gap for Mig back
print(avg_comparisons(explanandum_paper_only_leisure, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom,
hypothesis = "b2 = b6"))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b6 -0.434 0.286 -1.52 0.129 2.9 -0.995 0.127
##
## Type: response# Compare widenings
print(avg_comparisons(explanandum_paper_only_leisure, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
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.259 0.322 0.806 0.421 1.2 -0.371 0.89
##
## Type: responsepred <- avg_predictions(explanandum_paper_only_leisure,
newdata = all2010,
by = c("gender", "year", "mig_lang"),
vcov = ~ unique_classroom)
# Plot the results
# 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 mig back
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 mig back") +
ylab("Mean leisure satisfaction") +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
xlab("Year") +
labs(color = "Gender")
# With Mig back
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("Mig back") +
ylab("Mean leisure 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/paper_satis_leisure.png", width = width_plot, height = height_plot)
# Separate plots for manuscript
ggsave("Plots/paper_satis_leisure_no.png",
plot = cat_mig_lang_0 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)
ggsave("Plots/paper_satis_leisure_yes.png",
plot = cat_mig_lang_1 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)Satisfaction with your school grades
# Fit the model
explanandum_paper_only_grades <- lm(satis_grades ~ 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[combined$source == "paper",])
# Predictions
pred <- predictions(explanandum_paper_only_grades,
by = c("gender", "year"),
newdata = all2010,
vcov = ~ unique_classroom)
print(pred)##
## gender year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 2010 3.32 0.0518 64.1 <0.001 Inf 3.22 3.42
## female 2015 3.30 0.0515 64.1 <0.001 Inf 3.20 3.40
## female 2023 3.37 0.0656 51.4 <0.001 Inf 3.24 3.50
## male 2010 3.32 0.0433 76.6 <0.001 Inf 3.23 3.40
## male 2015 3.35 0.0410 81.7 <0.001 Inf 3.27 3.43
## male 2023 3.29 0.1171 28.1 <0.001 574.6 3.06 3.52
##
## Type: response# Gender gap by year
comps <- avg_comparisons(explanandum_paper_only_grades, 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.00416 0.0615 0.0677 0.946 0.1 -0.1164 0.125
## 2015 0.05193 0.0532 0.9758 0.329 1.6 -0.0524 0.156
## 2023 -0.07736 0.1146 -0.6748 0.500 1.0 -0.3020 0.147
##
## Term: gender
## Type: response
## Comparison: male - female# Test biggest contrast: 2010 versus 2023
max_comp <- avg_comparisons(explanandum_paper_only_grades, 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.0815 0.13 0.627 0.531 0.9 -0.173 0.336
##
## Type: responseprint(max_comp$estimate/sd(combined$satis_grades, na.rm = TRUE))## [1] 0.08610917# Gender gap by year and Mig back
comps_by_migback <- avg_comparisons(explanandum_paper_only_grades, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom)
print(comps_by_migback)##
## year mig_lang Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0 0.0426 0.0659 0.646 0.5180 0.9 -0.0865 0.1717
## 2010 1 -0.3397 0.1570 -2.163 0.0305 5.0 -0.6475 -0.0319
## 2015 0 0.0521 0.0550 0.947 0.3437 1.5 -0.0557 0.1599
## 2015 1 0.0505 0.1474 0.343 0.7319 0.5 -0.2385 0.3395
## 2023 0 -0.0975 0.1253 -0.778 0.4364 1.2 -0.3430 0.1480
## 2023 1 0.1030 0.1893 0.544 0.5863 0.8 -0.2679 0.4739
##
## Term: gender
## Type: response
## Comparison: male - female# Widening of the gap for no Mig back
print(avg_comparisons(explanandum_paper_only_grades, 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.14 0.142 0.99 0.322 1.6 -0.137 0.418
##
## Type: response# Widening of the gap for Mig back
print(avg_comparisons(explanandum_paper_only_grades, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom,
hypothesis = "b2 = b6"))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b6 -0.443 0.246 -1.8 0.0718 3.8 -0.925 0.0393
##
## Type: response# Compare widenings
print(avg_comparisons(explanandum_paper_only_grades, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
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.583 0.282 2.07 0.0386 4.7 0.0305 1.14
##
## Type: responsepred <- avg_predictions(explanandum_paper_only_grades,
newdata = all2010,
by = c("gender", "year", "mig_lang"),
vcov = ~ unique_classroom)
# Plot the results
# 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 mig back
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 mig back") +
ylab("Mean grades satisfaction") +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
xlab("Year") +
labs(color = "Gender")
# With Mig back
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("Mig back") +
ylab("Mean grades 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/paper_satis_grades.png", width = width_plot, height = height_plot)
# Separate plots for manuscript
ggsave("Plots/paper_satis_grades_no.png",
plot = cat_mig_lang_0 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)
ggsave("Plots/paper_satis_grades_yes.png",
plot = cat_mig_lang_1 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)Satisfaction with your housing situation
# Fit the model
explanandum_paper_only_dwell <- lm(satis_dwell ~ 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[combined$source == "paper",])
# Predictions
pred <- predictions(explanandum_paper_only_dwell,
by = c("gender", "year"),
newdata = all2010,
vcov = ~ unique_classroom)
print(pred)##
## gender year Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 2010 4.42 0.0325 135.9 <0.001 Inf 4.35 4.48
## female 2015 4.40 0.0342 128.4 <0.001 Inf 4.33 4.46
## female 2023 4.33 0.0664 65.3 <0.001 Inf 4.20 4.46
## male 2010 4.47 0.0289 154.9 <0.001 Inf 4.42 4.53
## male 2015 4.46 0.0316 141.1 <0.001 Inf 4.40 4.52
## male 2023 4.49 0.0520 86.5 <0.001 Inf 4.39 4.60
##
## Type: response# Gender gap by year
comps <- avg_comparisons(explanandum_paper_only_dwell, 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.0596 0.0397 1.50 0.1336 2.9 -0.0183 0.137
## 2015 0.0694 0.0523 1.33 0.1846 2.4 -0.0331 0.172
## 2023 0.1758 0.0981 1.79 0.0732 3.8 -0.0165 0.368
##
## Term: gender
## Type: response
## Comparison: male - female# Test biggest contrast: 2010 versus 2023
max_comp <- avg_comparisons(explanandum_paper_only_dwell, 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.116 0.106 -1.1 0.272 1.9 -0.324 0.0912
##
## Type: responseprint(max_comp$estimate/sd(combined$satis_dwell, na.rm = TRUE))## [1] -0.1408588# Gender gap by year and Mig back
comps_by_migback <- avg_comparisons(explanandum_paper_only_dwell, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom)
print(comps_by_migback)##
## year mig_lang Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 2010 0 0.0503 0.0445 1.132 0.2575 2.0 -0.0368 0.137
## 2010 1 0.1421 0.1373 1.035 0.3007 1.7 -0.1270 0.411
## 2015 0 0.0708 0.0553 1.281 0.2002 2.3 -0.0375 0.179
## 2015 1 0.0574 0.1269 0.452 0.6513 0.6 -0.1914 0.306
## 2023 0 0.1394 0.1041 1.340 0.1803 2.5 -0.0646 0.343
## 2023 1 0.5016 0.1969 2.547 0.0109 6.5 0.1156 0.888
##
## Term: gender
## Type: response
## Comparison: male - female# Widening of the gap for no Mig back
print(avg_comparisons(explanandum_paper_only_dwell, 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.0891 0.113 -0.787 0.431 1.2 -0.311 0.133
##
## Type: response# Widening of the gap for Mig back
print(avg_comparisons(explanandum_paper_only_dwell, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
vcov = ~ unique_classroom,
hypothesis = "b2 = b6"))##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b2=b6 -0.359 0.24 -1.5 0.134 2.9 -0.83 0.111
##
## Type: response# Compare widenings
print(avg_comparisons(explanandum_paper_only_dwell, variables = "gender", by = c("year", "mig_lang"),
newdata = all2010,
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.27 0.26 1.04 0.298 1.7 -0.239 0.78
##
## Type: responsepred <- avg_predictions(explanandum_paper_only_dwell,
newdata = all2010,
by = c("gender", "year", "mig_lang"),
vcov = ~ unique_classroom)
# Plot the results
# 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 mig back
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 mig back") +
ylab("Mean housing satisfaction") +
scale_x_continuous(breaks = c(2010, 2015, 2023)) +
xlab("Year") +
labs(color = "Gender")
# With Mig back
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("Mig back") +
ylab("Mean housing 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/paper_satis_dwell.png", width = width_plot, height = height_plot)
# Separate plots for manuscript
ggsave("Plots/paper_satis_dwell_no.png",
plot = cat_mig_lang_0 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)
ggsave("Plots/paper_satis_dwell_yes.png",
plot = cat_mig_lang_1 + theme(legend.position = "none"),
width = width_plot_single, height = height_plot_single)