In 2023, the survey was conducted with tablets whenever they were available in the schools. Here, we are going to look at the possibility that this might have affected how students report their life satisfaction.

Does life satisfaction vary by survey mode?

First, we will simply see whether gender and mode jointly predict life satisfaction in 2023.

mode_effects <- lm(satis ~ source*gender, data = combined[combined$year == 2023,])
summary(mode_effects)
## 
## Call:
## lm(formula = satis ~ source * gender, data = combined[combined$year == 
##     2023, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.91186 -0.37681  0.08814  0.62319  1.62319 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             3.37681    0.03298 102.376  < 2e-16 ***
## sourcepaper             0.16990    0.05853   2.903  0.00374 ** 
## gendermale              0.53505    0.04659  11.484  < 2e-16 ***
## sourcepaper:gendermale -0.17098    0.08379  -2.041  0.04143 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.822 on 1799 degrees of freedom
## Multiple R-squared:  0.0829, Adjusted R-squared:  0.08137 
## F-statistic: 54.21 on 3 and 1799 DF,  p-value: < 2.2e-16
pred <- avg_predictions(mode_effects, variables = c("source", "gender"), vcov = ~ unique_classroom)
print(pred)
## 
##  source gender Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
##  online female     3.38     0.0360  93.9   <0.001 Inf  3.31   3.45
##  online male       3.91     0.0274 142.8   <0.001 Inf  3.86   3.97
##  paper  female     3.55     0.0548  64.7   <0.001 Inf  3.44   3.65
##  paper  male       3.91     0.0524  74.7   <0.001 Inf  3.81   4.01
## 
## Type: response
ggplot(pred, aes(x = source, group = gender, color = gender, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_point() +
  geom_line() +
  geom_errorbar(width = .2) +
  coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
  theme_classic() +
  ylab("Mean life satisfaction") +
  xlab("Survey mode") +
  labs(color = "Gender")

ggsave("Plots/mode.png", width = 4, height = 3)


print(avg_predictions(mode_effects, variables = c("source", "gender"),
                hypothesis = "b1 = b3"), vcov = ~ unique_classroom) # Girls: Paper versus online
## Warning: 
## It is essential to check the order of estimates when specifying hypothesis tests using positional indices like b1, b2, etc. The indices of estimates can change depending on the order of rows in the original dataset, user-supplied arguments, model-fitting package, and version of `marginaleffects`.
## 
## It is also good practice to use assertions that ensure the order of estimates is consistent across different runs of the same code. Example:
## 
## ```r
## mod <- lm(mpg ~ am * carb, data = mtcars)
## 
## # assertion for safety
## p <- avg_predictions(mod, by = 'carb')
## stopifnot(p$carb[1] != 1 || p$carb[2] != 2)
## 
## # hypothesis test
## avg_predictions(mod, by = 'carb', hypothesis = 'b1 - b2 = 0')
## ```
## 
## Disable this warning with: `options(marginaleffects_safe = FALSE)`
##  This warning appears once per session.
## 
##  Hypothesis Estimate Std. Error    z Pr(>|z|)   S  2.5 %  97.5 %
##       b1=b3    -0.17     0.0585 -2.9   0.0037 8.1 -0.285 -0.0552
## 
## Type: response
print(avg_predictions(mode_effects, variables = c("source", "gender"),
                hypothesis = "b2 = b4"), vcov = ~ unique_classroom) # Boys: paper versus online
## 
##  Hypothesis Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
##       b2=b4  0.00108       0.06 0.018    0.986 0.0 -0.116  0.119
## 
## Type: response
# Girls are less satisfied on the tablet
# Boys are pretty much the same

# Compare the predictions
print(avg_predictions(mode_effects, variables = c("source", "gender", "year"),
                hypothesis = "b1-b3 = b2-b4", vcov = ~ unique_classroom))
## 
##   Hypothesis Estimate Std. Error     z Pr(>|z|)   S  2.5 %  97.5 %
##  b1-b3=b2-b4   -0.171     0.0772 -2.21   0.0269 5.2 -0.322 -0.0196
## 
## Type: response

Let’s look at whether survey mode covaries with other variables

Maybe survey mode predicts life satisfaction because it covaries with other factors? To check this possibility, let’s see how the samples differ

# Extract 2023 final data for convenience
combined_2023 <- combined[combined$year == 2023,]
prop.table(table(combined_2023$source))
## 
##    online     paper 
## 0.6905158 0.3094842
# Which schools have which assessment mode(s)
has_tablet <- unique(combined_2023$school_id[combined_2023$source == "online" & !is.na(combined_2023$school_id)])
has_paper <- unique(combined_2023$school_id[combined_2023$source == "paper" & !is.na(combined_2023$school_id)])
has_tablet
##  [1] 49  2  3 58 41 60 68 13 20 63 50 72 40 10 66 15 16 71 52 18  6 59 57 55 61
## [26] 53 48 11 14
has_paper
##  [1]  2 60 13 63 40 66 15 16 17 71 46  7 59 57 53 21
# Assessment mode table
assessment_modes <- data.frame(matrix(NA, ncol = 3,nrow = length(unique(combined_2023$school_id[!is.na(combined_2023$school_id)]))))
names(assessment_modes) <- c("school_id", "has_paper", "has_tablet")
assessment_modes$school_id <- unique(combined_2023$school_id[!is.na(combined_2023$school_id)])

for (i in 1:nrow(assessment_modes)) {
  assessment_modes$has_paper[i] <-  assessment_modes$school_id[i] %in% has_paper
  assessment_modes$has_tablet[i] <-  assessment_modes$school_id[i] %in% has_tablet
}

assessment_modes$count <- assessment_modes$has_paper + assessment_modes$has_tablet
prop.table(table(assessment_modes$count))
## 
##         1         2 
## 0.6363636 0.3636364
# 63 percent of schools had only one mode, 36 percent of schools had both

# Is response mode correlated with the type of school?
prop.table(table(combined_2023$schooltype, combined_2023$source), margin = 1)
##    
##        online     paper
##   1 0.8876404 0.1123596
##   2 0.5864407 0.4135593
chisq.test(table(combined_2023$schooltype, combined_2023$source))
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(combined_2023$schooltype, combined_2023$source)
## X-squared = 171.68, df = 1, p-value < 2.2e-16
# Significantly more tablets in Mittelschule

# Is response mode correlated with age?
prop.table(table(combined_2023$age, combined_2023$source), margin = 1)
##     
##         online     paper
##   12 0.6194690 0.3805310
##   13 0.5949008 0.4050992
##   14 0.7326478 0.2673522
##   15 0.7459016 0.2540984
##   16 0.6577181 0.3422819
##   17 0.7028571 0.2971429
##   18 0.8073394 0.1926606
t.test(combined_2023$age[combined_2023$source == "paper"], combined_2023$age[combined_2023$source == "online"])
## 
##  Welch Two Sample t-test
## 
## data:  combined_2023$age[combined_2023$source == "paper"] and combined_2023$age[combined_2023$source == "online"]
## t = -3.4147, df = 1066, p-value = 0.0006626
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4411449 -0.1191695
## sample estimates:
## mean of x mean of y 
##  14.55197  14.83213
mean(combined_2023$age[combined_2023$source == "paper"], na.rm = TRUE)
## [1] 14.55197
mean(combined_2023$age[combined_2023$source == "online"], na.rm = TRUE)
## [1] 14.83213
# Students who filled out on tablet are slightly but significantly older

# Is response mode correlated with migration background?
prop.table(table(combined_2023$mig_lang, combined_2023$source), margin = 1)
##    
##        online     paper
##   0 0.6859873 0.3140127
##   1 0.7210300 0.2789700
chisq.test(table(combined_2023$mig_lang, combined_2023$source))
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(combined_2023$mig_lang, combined_2023$source)
## X-squared = 1.0076, df = 1, p-value = 0.3155

Does life satisfaction vary by response mode? Let’s check with controls

To rule out that the associations between survey mode and life satisfaction are due to confounding with demographics, let’s again adjust for demographics

library(splines)

mode_effects_controlled <- lm(satis ~ source*gender +
                                as.factor(schooltype) + as.factor(schooltype):source + as.factor(schooltype):gender +
                                as.factor(mig_lang) + as.factor(mig_lang):source + as.factor(mig_lang):gender +
                                bs(age, df = 3) + bs(age, df = 3):source + bs(age, df = 3):gender, data = combined[combined$year == 2023,])
summary(mode_effects_controlled)
## 
## Call:
## lm(formula = satis ~ source * gender + as.factor(schooltype) + 
##     as.factor(schooltype):source + as.factor(schooltype):gender + 
##     as.factor(mig_lang) + as.factor(mig_lang):source + as.factor(mig_lang):gender + 
##     bs(age, df = 3) + bs(age, df = 3):source + bs(age, df = 3):gender, 
##     data = combined[combined$year == 2023, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.00649 -0.41965  0.06273  0.56244  1.79345 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         3.62683    0.11261  32.207  < 2e-16 ***
## sourcepaper                         0.19217    0.17882   1.075  0.28268    
## gendermale                          0.47987    0.15684   3.060  0.00225 ** 
## as.factor(schooltype)2             -0.05443    0.06740  -0.807  0.41950    
## as.factor(mig_lang)1               -0.10794    0.09029  -1.195  0.23208    
## bs(age, df = 3)1                   -0.52547    0.28085  -1.871  0.06151 .  
## bs(age, df = 3)2                   -0.07846    0.20143  -0.390  0.69694    
## bs(age, df = 3)3                   -0.18745    0.15120  -1.240  0.21525    
## sourcepaper:gendermale             -0.16204    0.08909  -1.819  0.06908 .  
## sourcepaper:as.factor(schooltype)2  0.06774    0.11775   0.575  0.56517    
## gendermale:as.factor(schooltype)2  -0.01428    0.08887  -0.161  0.87233    
## sourcepaper:as.factor(mig_lang)1    0.15674    0.12806   1.224  0.22112    
## gendermale:as.factor(mig_lang)1     0.20403    0.11550   1.766  0.07749 .  
## sourcepaper:bs(age, df = 3)1        0.05965    0.40285   0.148  0.88230    
## sourcepaper:bs(age, df = 3)2       -0.35985    0.28925  -1.244  0.21364    
## sourcepaper:bs(age, df = 3)3       -0.13804    0.24577  -0.562  0.57441    
## gendermale:bs(age, df = 3)1         0.24909    0.38556   0.646  0.51832    
## gendermale:bs(age, df = 3)2         0.04578    0.26445   0.173  0.86258    
## gendermale:bs(age, df = 3)3        -0.24087    0.22411  -1.075  0.28262    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.818 on 1784 degrees of freedom
## Multiple R-squared:  0.09931,    Adjusted R-squared:  0.09023 
## F-statistic: 10.93 on 18 and 1784 DF,  p-value: < 2.2e-16
# Take counterfactual contrasts: everybody with paper vs everybody online

print(predictions(mode_effects_controlled, by = c("source", "gender"), vcov = ~ unique_classroom))
## 
##  source gender Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
##  online female     3.38     0.0337 100.1   <0.001 Inf  3.31   3.44
##  online male       3.91     0.0247 158.6   <0.001 Inf  3.86   3.96
##  paper  female     3.55     0.0535  66.3   <0.001 Inf  3.44   3.65
##  paper  male       3.91     0.0463  84.5   <0.001 Inf  3.82   4.00
## 
## Type: response
contrasts <- avg_comparisons(mode_effects_controlled, variables = "source", by = "gender", vcov = ~ unique_classroom)
print(contrasts)
## 
##  gender Estimate Std. Error      z Pr(>|z|)   S   2.5 % 97.5 %
##  female  0.15995     0.0671  2.385   0.0171 5.9  0.0285  0.291
##  male   -0.00606     0.0543 -0.112   0.9111 0.1 -0.1125  0.100
## 
## Term: source
## Type: response
## Comparison: paper - online
# Girls are less satisfied on the tablet
# Boys are pretty much the same

print(avg_comparisons(mode_effects_controlled, variables = "source", by = "gender", vcov = ~ unique_classroom,
                hypothesis = "b1 = b2"))
## 
##  Hypothesis Estimate Std. Error    z Pr(>|z|)   S  2.5 % 97.5 %
##       b1=b2    0.166     0.0772 2.15   0.0315 5.0 0.0147  0.317
## 
## Type: response

School fixed effects

Of course, survey mode may also correlate with some other unobserved variable (e.g., SES) that ultimately explains its associations with life satisfaction. As survey mode was randomized, we cannot rule out this possibility. However, let’s assume that within schools, these unobserved variables are approximately constant. In that scenario, we can compare classrooms that used different modes within the same school.

So, if we still find that survey mode is associated with life satisfaction, we can rule out that the association is explained away by unobserved variables on the school level.

# Determine which schools used both assessment modes
schools_with_both <- assessment_modes$school_id[assessment_modes$count == 2]

# Run a fixed effects model
# These model include a school-wise intercept 
# and allow the effect of gender to vary between schools
# (in completely arbitrary ways)

mode_effects_controlled_fe <- lm(satis ~ source*gender +
                                   as.factor(schooltype) + as.factor(schooltype):source + as.factor(schooltype):gender +
                                   as.factor(mig_lang) + as.factor(mig_lang):source + as.factor(mig_lang):gender +
                                   bs(age, df = 3) + bs(age, df = 3):source + bs(age, df = 3):gender +
                                   as.factor(school_id) + as.factor(school_id):gender, 
                                 data = combined[combined$year == 2023 & combined$school_id %in% schools_with_both,])
summary(mode_effects_controlled_fe)
## 
## Call:
## lm(formula = satis ~ source * gender + as.factor(schooltype) + 
##     as.factor(schooltype):source + as.factor(schooltype):gender + 
##     as.factor(mig_lang) + as.factor(mig_lang):source + as.factor(mig_lang):gender + 
##     bs(age, df = 3) + bs(age, df = 3):source + bs(age, df = 3):gender + 
##     as.factor(school_id) + as.factor(school_id):gender, data = combined[combined$year == 
##     2023 & combined$school_id %in% schools_with_both, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1005 -0.4858  0.1053  0.5142  2.1399 
## 
## Coefficients: (2 not defined because of singularities)
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         3.83532    0.31650  12.118   <2e-16 ***
## sourcepaper                         0.29560    0.32807   0.901   0.3679    
## gendermale                          0.53043    0.38650   1.372   0.1703    
## as.factor(schooltype)2             -0.21051    0.31620  -0.666   0.5058    
## as.factor(mig_lang)1               -0.28403    0.16625  -1.708   0.0880 .  
## bs(age, df = 3)1                   -0.74345    0.45824  -1.622   0.1051    
## bs(age, df = 3)2                    0.02417    0.33153   0.073   0.9419    
## bs(age, df = 3)3                   -0.52380    0.25327  -2.068   0.0390 *  
## as.factor(school_id)13              0.21291    0.21502   0.990   0.3224    
## as.factor(school_id)15              0.20731    0.21730   0.954   0.3404    
## as.factor(school_id)16             -0.05594    0.20664  -0.271   0.7867    
## as.factor(school_id)40              0.13049    0.25135   0.519   0.6038    
## as.factor(school_id)53             -0.07619    0.29122  -0.262   0.7937    
## as.factor(school_id)57             -0.55207    0.38593  -1.430   0.1530    
## as.factor(school_id)59              0.17605    0.23201   0.759   0.4482    
## as.factor(school_id)60             -0.14545    0.21087  -0.690   0.4905    
## as.factor(school_id)63              0.14159    0.20314   0.697   0.4860    
## as.factor(school_id)66             -0.09659    0.33586  -0.288   0.7737    
## as.factor(school_id)71                   NA         NA      NA       NA    
## sourcepaper:gendermale             -0.28683    0.15808  -1.814   0.0700 .  
## sourcepaper:as.factor(schooltype)2 -0.16381    0.26015  -0.630   0.5291    
## gendermale:as.factor(schooltype)2   0.03450    0.39647   0.087   0.9307    
## sourcepaper:as.factor(mig_lang)1    0.29052    0.19876   1.462   0.1442    
## gendermale:as.factor(mig_lang)1     0.28909    0.20106   1.438   0.1509    
## sourcepaper:bs(age, df = 3)1        0.26197    0.61441   0.426   0.6700    
## sourcepaper:bs(age, df = 3)2       -0.35685    0.47841  -0.746   0.4560    
## sourcepaper:bs(age, df = 3)3        0.22481    0.44639   0.504   0.6147    
## gendermale:bs(age, df = 3)1        -0.06913    0.59215  -0.117   0.9071    
## gendermale:bs(age, df = 3)2        -0.03315    0.41625  -0.080   0.9365    
## gendermale:bs(age, df = 3)3         0.19402    0.36007   0.539   0.5902    
## gendermale:as.factor(school_id)13  -0.12690    0.29004  -0.438   0.6619    
## gendermale:as.factor(school_id)15  -0.63514    0.30702  -2.069   0.0389 *  
## gendermale:as.factor(school_id)16  -0.02628    0.28229  -0.093   0.9258    
## gendermale:as.factor(school_id)40  -0.03589    0.35950  -0.100   0.9205    
## gendermale:as.factor(school_id)53   0.14303    0.37771   0.379   0.7050    
## gendermale:as.factor(school_id)57   0.23982    0.56431   0.425   0.6710    
## gendermale:as.factor(school_id)59  -0.00223    0.33763  -0.007   0.9947    
## gendermale:as.factor(school_id)60   0.14039    0.30221   0.465   0.6424    
## gendermale:as.factor(school_id)63   0.01959    0.27675   0.071   0.9436    
## gendermale:as.factor(school_id)66  -0.16776    0.43020  -0.390   0.6967    
## gendermale:as.factor(school_id)71        NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8475 on 778 degrees of freedom
## Multiple R-squared:  0.1192, Adjusted R-squared:  0.07615 
## F-statistic:  2.77 on 38 and 778 DF,  p-value: 1.286e-07
pred <- avg_predictions(mode_effects_controlled_fe, variables = c("source", "gender"), vcov = ~ unique_classroom)

ggplot(pred, aes(x = source, group = gender, color = gender, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_point() +
  geom_line() +
  geom_errorbar(width = .2) +
  coord_cartesian(ylim = c(mean(combined$satis) - sd(combined$satis), mean(combined$satis) + sd(combined$satis))) +
  theme_classic() +
  ylab("Mean life satisfaction") +
  xlab("Survey mode") +
  labs(color = "Gender")

ggsave("Plots/mode_fe.png", width = 4, height = 3)

# Let's again look at mode effects by gender
contrasts <- avg_comparisons(mode_effects_controlled_fe, variables = "source", by = "gender", vcov = ~ unique_classroom)
print(contrasts)
## 
##  gender Estimate Std. Error      z Pr(>|z|)   S   2.5 % 97.5 %
##  female   0.2317     0.0813  2.850  0.00437 7.8  0.0723 0.3910
##  male    -0.0545     0.0609 -0.895  0.37089 1.4 -0.1738 0.0648
## 
## Term: source
## Type: response
## Comparison: paper - online
# Once again, girls are less satisfied on tablet
# Boys are pretty much the same

print(avg_comparisons(mode_effects_controlled_fe, variables = "source", by = "gender", vcov = ~ unique_classroom,
                hypothesis = "b1 = b2"))
## 
##  Hypothesis Estimate Std. Error    z Pr(>|z|)   S  2.5 % 97.5 %
##       b1=b2    0.286      0.115 2.49   0.0129 6.3 0.0606  0.512
## 
## Type: response
# This speaks against the idea that survey mode is correlated with some unobserved school-level
# variable that ultimately explains differences in the gender gaps

What happens to the explanandum if we take into account survey mode?

So, let’s assume that the (gender-specific) associations between survey mode and life satisfaction actually reflect mode effects: answering the survey on a tablet makes girls (but not boys) report lower levels of satisfaction. What does the widening of the gender gap look like if we hold survey mode constant? Tablets were only introduced in 2023, and so we cannot see how the gender gap widens for this assessment mode. However, paper was used throughout, and so we can look at the widening of the gender gap for this survey mode

explanandum_paper_only <- 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[combined$source == "paper",])

# We will evaluate this for various data grids
# 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)

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

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

# Demographics like 2023, but only those who actually used paper
temp1 <- combined[combined$year == 2023 & combined$source == "paper",]
temp1$year <- 2010
temp2 <- combined[combined$year == 2023 & combined$source == "paper",]
temp2$year <- 2015
temp3 <- combined[combined$year == 2023 & combined$source == "paper",]
temp3$year <- 2023
all2023_paper <- rbind(temp1, temp2, temp3)


##########################
# Use 2010 distribution
##########################
print(avg_comparisons(explanandum_paper_only, variables = "gender", by = "year",
                 newdata = all2010,
                vcov = ~ unique_classroom))
## 
##  year Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
##  2010    0.223     0.0463 4.82   <0.001 19.4 0.133  0.314
##  2015    0.309     0.0493 6.26   <0.001 31.3 0.212  0.405
##  2023    0.412     0.0946 4.35   <0.001 16.2 0.226  0.597
## 
## Term: gender
## Type: response
## Comparison: male - female
max_comp <- avg_comparisons(explanandum_paper_only, variables = "gender", by = "year",
                hypothesis = "b1 = b3",
                 newdata = all2010,
                vcov = ~ unique_classroom)

# Widening of gender gap
print(max_comp)
## 
##  Hypothesis Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
##       b1=b3   -0.188      0.105 -1.79   0.0742 3.8 -0.395 0.0184
## 
## Type: response
print(max_comp$estimate/sd(combined$satis))
## [1] -0.2216236
##########################
# Use 2015 distribution
##########################
max_comp <- avg_comparisons(explanandum_paper_only, variables = "gender", by = "year",
                hypothesis = "b1 = b3",
                 newdata = all2015,
                vcov = ~ unique_classroom)
# Widening of gender gap
print(max_comp)
## 
##  Hypothesis Estimate Std. Error     z Pr(>|z|)   S  2.5 %  97.5 %
##       b1=b3   -0.246      0.113 -2.17   0.0304 5.0 -0.468 -0.0233
## 
## Type: response
##########################
# Use 2023 distribution
##########################
pred <- predictions(explanandum_paper_only,
                 by = c("gender", "year"),
                 newdata = all2023,
                 vcov = ~ unique_classroom)

comps <- avg_comparisons(explanandum_paper_only, variables = "gender", by = "year",
                 newdata = all2023,
                 vcov = ~ unique_classroom)
print(comps)
## 
##  year Estimate Std. Error    z Pr(>|z|)    S  2.5 % 97.5 %
##  2010    0.169     0.0563 3.00  0.00274  8.5 0.0583  0.279
##  2015    0.279     0.0414 6.74  < 0.001 35.9 0.1981  0.360
##  2023    0.408     0.0829 4.93  < 0.001 20.2 0.2460  0.571
## 
## Term: gender
## Type: response
## Comparison: male - female
max_comp <- avg_comparisons(explanandum_paper_only, variables = "gender", by = "year",
                hypothesis = "b1 = b3",
                 newdata = all2023,
                vcov = ~ unique_classroom)
# Widening of gender gap
print(max_comp)
## 
##  Hypothesis Estimate Std. Error     z Pr(>|z|)   S  2.5 %  97.5 %
##       b1=b3    -0.24        0.1 -2.39   0.0167 5.9 -0.436 -0.0435
## 
## Type: response

Paper only, modification by migback

Before, we had already noted that the widening of the gender gap strongly differed depending on migration background/language at home. So let’s see how migration background modifies the widening of the gender gap if we hold constant survey mode by only analyzing responses collected on paper.

comps_by_migback <- avg_comparisons(explanandum_paper_only, 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.2542     0.0465  5.464  < 0.001 24.4  0.1630  0.345
##  2010        1  -0.0519     0.1771 -0.293  0.76928  0.4 -0.3990  0.295
##  2015        0   0.3095     0.0506  6.115  < 0.001 29.9  0.2103  0.409
##  2015        1   0.3025     0.1399  2.163  0.03056  5.0  0.0284  0.577
##  2023        0   0.3656     0.0999  3.659  < 0.001 11.9  0.1698  0.561
##  2023        1   0.8228     0.2508  3.281  0.00104  9.9  0.3312  1.314
## 
## Term: gender
## Type: response
## Comparison: male - female
# 2010 demographics
# Widening of the gap for no migration background
comp_no <- avg_comparisons(explanandum_paper_only, variables = "gender", by = c("year", "mig_lang"),
                           newdata = all2010,
                           vcov = ~ unique_classroom,
                           hypothesis = "b1 = b5")

print(comp_no)
## 
##  Hypothesis Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
##       b1=b5   -0.111       0.11 -1.01    0.312 1.7 -0.327  0.105
## 
## Type: response
print(comp_no$estimate/sd(combined$satis))
## [1] -0.1312453
# Widening of the gap for migration background
comp_yes <- avg_comparisons(explanandum_paper_only, variables = "gender", by = c("year", "mig_lang"),
                           newdata = all2010,
                           vcov = ~ unique_classroom,
                           hypothesis = "b2 = b6")

print(comp_yes)
## 
##  Hypothesis Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
##       b2=b6   -0.875      0.307 -2.85  0.00438 7.8 -1.48 -0.273
## 
## Type: response
print(comp_yes$estimate/sd(combined$satis))
## [1] -1.030477
# Comparison of the widenings by migration background
print(avg_comparisons(explanandum_paper_only, 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.763      0.321 2.38   0.0173 5.9 0.135   1.39
## 
## Type: response

Plot the implied means of this final model

# Generate predictions for combinations of gender and year and age
pred <- avg_predictions(explanandum_paper_only,
                        newdata = all2023,
                        by = c("gender", "year", "mig_lang"),
                        vcov = ~ unique_classroom)

# 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/paper_satis_by_mig_lang_categorical.png", width = 7, height = 3)