The explanandum: Widening gender gap in general life satisfaction

First, we want to reproduce the finding from the press release that states that girls in particular have become unhappier.

load("prep.RData")
library(haven)
library(marginaleffects)
library(ggplot2)

# Simple model that estimates the means of life satisfaction by gender, survey year, and their interaction
explanandum <- lm(satis ~ as.factor(year) + gender + gender:as.factor(year), 
                  data = combined)

# We can take a look at the coefficients which map onto what we are interested in in a rather straightforward manner
summary(explanandum)
## 
## Call:
## lm(formula = satis ~ as.factor(year) + gender + gender:as.factor(year), 
##     data = combined)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.02801 -0.43077  0.08847  0.56923  1.56923 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     3.78805    0.03107 121.929  < 2e-16 ***
## as.factor(year)2015            -0.03921    0.04184  -0.937    0.349    
## as.factor(year)2023            -0.35728    0.04136  -8.638  < 2e-16 ***
## gendermale                      0.21783    0.04431   4.916 9.11e-07 ***
## as.factor(year)2015:gendermale  0.06134    0.06082   1.009    0.313    
## as.factor(year)2023:gendermale  0.26293    0.05889   4.465 8.21e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8237 on 4758 degrees of freedom
## Multiple R-squared:  0.05928,    Adjusted R-squared:  0.05829 
## F-statistic: 59.96 on 5 and 4758 DF,  p-value: < 2.2e-16
# However, do get directly started, with the marginaleffects framework...
# Here we generate predicted values by gender and year
pred <- predictions(explanandum,
                 by = c("gender", "year"))


# Plot the model-implied means
# Figure 1 in the manuscript

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


# Calculating the gender gap by year
# With the help of average counterfactual comparisons
# (the averaging here does not really matter given that our model does not
# allow for any variability within the groups)

comps <- avg_comparisons(explanandum, 
                         variables = "gender", # counterfactual comparisons over gender (what if boy/girl)
                         by = "year") # separately for each year

# Take a look at the resulting gender gaps
print(comps)
## 
##  year Estimate Std. Error     z Pr(>|z|)     S 2.5 % 97.5 %
##  2010    0.218     0.0443  4.92   <0.001  20.1 0.131  0.305
##  2015    0.279     0.0417  6.70   <0.001  35.5 0.198  0.361
##  2023    0.481     0.0388 12.39   <0.001 114.7 0.405  0.557
## 
## Term: gender
## Type: response
## Comparison: male - female
# Test maximum contrast (gender gap in 2010 vs gender gap in 2023)
max_comp <- avg_comparisons(explanandum, variables = "gender", by = "year",
                hypothesis = "b1 = b3") # This tests the first versus the third number in the output
## 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.263     0.0589 -4.46   <0.001 16.9 -0.378 -0.148
## 
## Type: response
# We can also compare the other gender gaps
# Compare gender gap in 2010 with gender gap in 2015
print(avg_comparisons(explanandum, variables = "gender", by = "year",
                hypothesis = "b1 = b2"))
## 
##  Hypothesis Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
##       b1=b2  -0.0613     0.0608 -1.01    0.313 1.7 -0.181 0.0579
## 
## Type: response
# Compare gender gap in 2015 with gender gap in 2023
print(avg_comparisons(explanandum, variables = "gender", by = "year",
                hypothesis = "b2 = b3"))
## 
##  Hypothesis Estimate Std. Error     z Pr(>|z|)    S  2.5 % 97.5 %
##       b2=b3   -0.202     0.0569 -3.54   <0.001 11.3 -0.313  -0.09
## 
## Type: response
# Express the gender gaps in SD units (SD of life satisfaction across all observations)
round(comps$estimate[1]/sd(combined$satis), 2)
## [1] 0.26
round(comps$estimate[2]/sd(combined$satis), 2)
## [1] 0.33
round(comps$estimate[3]/sd(combined$satis), 2)
## [1] 0.57
# Express widening of gender gap in SDs 
# Denominator is SD of life satisfaction across all observations
round(max_comp$estimate/sd(combined$satis), 2)
## [1] -0.31

Alternative measures of satisfaction

Let’s repeat the previous analyses for all satisfaction measures in our data that are available for all three years. Here, we also generate the individual panels of Figure 2.

Satisfaction with the amount of money with which you have to get by (financial satisfaction)

# Estimate means
model_money <- lm(satis_money ~ as.factor(year) + gender + gender:as.factor(year), 
                  data = combined)
summary(model_money)
## 
## Call:
## lm(formula = satis_money ~ as.factor(year) + gender + gender:as.factor(year), 
##     data = combined)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9921 -0.7660  0.0885  1.0079  1.2404 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.759602   0.036740 102.331  < 2e-16 ***
## as.factor(year)2015            0.006405   0.049543   0.129  0.89713    
## as.factor(year)2023            0.078504   0.048937   1.604  0.10874    
## gendermale                     0.151903   0.052434   2.897  0.00378 ** 
## as.factor(year)2015:gendermale 0.014294   0.072071   0.198  0.84280    
## as.factor(year)2023:gendermale 0.002109   0.069735   0.030  0.97588    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9741 on 4738 degrees of freedom
##   (20 observations deleted due to missingness)
## Multiple R-squared:  0.007918,   Adjusted R-squared:  0.006871 
## F-statistic: 7.563 on 5 and 4738 DF,  p-value: 4.384e-07
pred <- predictions(model_money,
                 by = c("gender", "year"))

# Test maximum contrast
max_comp <- avg_comparisons(model_money, 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.00211     0.0697 -0.0302    0.976 0.0 -0.139  0.135
## 
## Type: response
# Express in SD units
round(max_comp$estimate/sd(combined$satis_money, na.rm = TRUE), 2)
## [1] 0
# Plot
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 financial satisfaction") +
  xlab("Survey year") +
  labs(color = "Gender")

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

Satisfaction with your relationships to friends

# Estimate means
model_friends <- lm(satis_friends ~ as.factor(year) + gender + gender:as.factor(year), 
                  data = combined)
summary(model_friends)
## 
## Call:
## lm(formula = satis_friends ~ as.factor(year) + gender + gender:as.factor(year), 
##     data = combined)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3894 -0.3697 -0.0617  0.6303  0.9383 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     4.389444   0.030239 145.158  < 2e-16 ***
## as.factor(year)2015            -0.110698   0.040729  -2.718 0.006594 ** 
## as.factor(year)2023            -0.327702   0.040263  -8.139 5.05e-16 ***
## gendermale                     -0.004488   0.043126  -0.104 0.917122    
## as.factor(year)2015:gendermale  0.095431   0.059264   1.610 0.107407    
## as.factor(year)2023:gendermale  0.209638   0.057344   3.656 0.000259 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8006 on 4735 degrees of freedom
##   (23 observations deleted due to missingness)
## Multiple R-squared:  0.02094,    Adjusted R-squared:  0.0199 
## F-statistic: 20.25 on 5 and 4735 DF,  p-value: < 2.2e-16
pred <- predictions(model_friends,
                 by = c("gender", "year"))
# Test maximum contrast
max_comp <- avg_comparisons(model_friends, 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.21     0.0573 -3.66   <0.001 11.9 -0.322 -0.0972
## 
## Type: response
# Express in SD units
round(max_comp$estimate/sd(combined$satis_friends, na.rm = TRUE), 2)
## [1] -0.26
# Plot
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 friendship satisfaction") +
  xlab("Survey year") +
  labs(color = "Gender")

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

Satisfaction with your relationship to your mother

# Estimate means
model_mom <- lm(satis_mom ~ as.factor(year) + gender + gender:as.factor(year), 
                  data = combined)
summary(model_mom)
## 
## Call:
## lm(formula = satis_mom ~ as.factor(year) + gender + gender:as.factor(year), 
##     data = combined)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5130 -0.3621  0.4870  0.7171  0.7699 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     4.282857   0.035452 120.809  < 2e-16 ***
## as.factor(year)2015            -0.052717   0.047797  -1.103  0.27011    
## as.factor(year)2023            -0.046393   0.047212  -0.983  0.32582    
## gendermale                      0.133072   0.050541   2.633  0.00849 ** 
## as.factor(year)2015:gendermale -0.001119   0.069473  -0.016  0.98715    
## as.factor(year)2023:gendermale  0.143429   0.067219   2.134  0.03291 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.938 on 4727 degrees of freedom
##   (31 observations deleted due to missingness)
## Multiple R-squared:  0.01261,    Adjusted R-squared:  0.01156 
## F-statistic: 12.07 on 5 and 4727 DF,  p-value: 1.221e-11
pred <- predictions(model_mom,
                 by = c("gender", "year"))
# Test maximum contrast
max_comp <- avg_comparisons(model_mom, 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.143     0.0672 -2.13   0.0329 4.9 -0.275 -0.0117
## 
## Type: response
# Express in SD units
round(max_comp$estimate/sd(combined$satis_mom, na.rm = TRUE), 2)
## [1] -0.15
# Plot
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 satisfaction w mother") +
  xlab("Survey year") +
  labs(color = "Gender")

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

Satisfaction with your relationship to your father

# Estimate means
model_dad <- lm(satis_dad ~ as.factor(year) + gender + gender:as.factor(year), 
                  data = combined)
summary(model_dad)
## 
## Call:
## lm(formula = satis_dad ~ as.factor(year) + gender + gender:as.factor(year), 
##     data = combined)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3107 -0.8080  0.1920  0.8448  1.1920 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     3.82353    0.04504  84.898  < 2e-16 ***
## as.factor(year)2015             0.13232    0.06062   2.183   0.0291 *  
## as.factor(year)2023            -0.01554    0.05966  -0.260   0.7945    
## gendermale                      0.33169    0.06393   5.188 2.21e-07 ***
## as.factor(year)2015:gendermale -0.12328    0.08786  -1.403   0.1607    
## as.factor(year)2023:gendermale  0.17097    0.08474   2.018   0.0437 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.174 on 4659 degrees of freedom
##   (99 observations deleted due to missingness)
## Multiple R-squared:  0.02583,    Adjusted R-squared:  0.02479 
## F-statistic: 24.71 on 5 and 4659 DF,  p-value: < 2.2e-16
pred <- predictions(model_dad,
                 by = c("gender", "year"))

# Test maximum contrast
max_comp <- avg_comparisons(model_dad, 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.171     0.0847 -2.02   0.0436 4.5 -0.337 -0.00488
## 
## Type: response
# Express in SD units
round(max_comp$estimate/sd(combined$satis_dad, na.rm = TRUE), 2)
## [1] -0.14
# Plot
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 satisfaction w father") +
  xlab("Survey year") +
  labs(color = "Gender")

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

Satisfaction with your leisute time activities

# Estimate means
model_leisure <- lm(satis_leisure ~ as.factor(year) + gender + gender:as.factor(year), 
                  data = combined)
summary(model_leisure)
## 
## Call:
## lm(formula = satis_leisure ~ as.factor(year) + gender + gender:as.factor(year), 
##     data = combined)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.07725 -0.59647  0.01182  0.92275  1.40353 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     3.86610    0.03476 111.231  < 2e-16 ***
## as.factor(year)2015            -0.01476    0.04683  -0.315   0.7526    
## as.factor(year)2023            -0.26962    0.04629  -5.824 6.12e-09 ***
## gendermale                      0.12209    0.04961   2.461   0.0139 *  
## as.factor(year)2015:gendermale  0.10383    0.06809   1.525   0.1274    
## as.factor(year)2023:gendermale  0.33336    0.06598   5.053 4.52e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9209 on 4739 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.03236,    Adjusted R-squared:  0.03134 
## F-statistic:  31.7 on 5 and 4739 DF,  p-value: < 2.2e-16
pred <- predictions(model_leisure,
                 by = c("gender", "year"))

# Test maximum contrast
max_comp <- avg_comparisons(model_leisure, 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.333      0.066 -5.05   <0.001 21.1 -0.463 -0.204
## 
## Type: response
# Express in SD units
round(max_comp$estimate/sd(combined$satis_leisure, na.rm = TRUE), 2)
## [1] -0.36
# Plot
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 leisure satisfaction") +
  xlab("Survey year") +
  labs(color = "Gender")

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

Satisfaction with your school grades

# Estimate means
model_grades <- lm(satis_grades ~ as.factor(year) + gender + gender:as.factor(year), 
                  data = combined)
summary(model_grades)
## 
## Call:
## lm(formula = satis_grades ~ as.factor(year) + gender + gender:as.factor(year), 
##     data = combined)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3834 -0.3263 -0.2809  0.6823  1.7191 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     3.317664   0.035729  92.855   <2e-16 ***
## as.factor(year)2015            -0.036778   0.048178  -0.763    0.445    
## as.factor(year)2023             0.008687   0.047588   0.183    0.855    
## gendermale                      0.001391   0.050993   0.027    0.978    
## as.factor(year)2015:gendermale  0.029430   0.070063   0.420    0.674    
## as.factor(year)2023:gendermale  0.055667   0.067772   0.821    0.411    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9467 on 4739 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.001152,   Adjusted R-squared:  9.844e-05 
## F-statistic: 1.093 on 5 and 4739 DF,  p-value: 0.3617
pred <- predictions(model_grades,
                 by = c("gender", "year"))

# Test maximum contrast
max_comp <- avg_comparisons(model_grades, 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.0557     0.0678 -0.821    0.411 1.3 -0.188 0.0772
## 
## Type: response
# Express in SD units
round(max_comp$estimate/sd(combined$satis_grades, na.rm = TRUE), 2)
## [1] -0.06
# Plot
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 school grades satisfaction") +
  xlab("Survey year") +
  labs(color = "Gender")

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

Satisfaction with your housing situation

# Estimate means
model_dwell <- lm(satis_dwell ~ as.factor(year) + gender + gender:as.factor(year), 
                  data = combined)
summary(model_dwell)
## 
## Call:
## lm(formula = satis_dwell ~ as.factor(year) + gender + gender:as.factor(year), 
##     data = combined)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4720 -0.4268  0.5280  0.5843  0.7387 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     4.415714   0.031097 141.999  < 2e-16 ***
## as.factor(year)2015            -0.049009   0.041893  -1.170 0.242119    
## as.factor(year)2023            -0.154413   0.041392  -3.730 0.000193 ***
## gendermale                      0.056262   0.044333   1.269 0.204474    
## as.factor(year)2015:gendermale  0.003793   0.060884   0.062 0.950328    
## as.factor(year)2023:gendermale  0.133615   0.058919   2.268 0.023387 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8227 on 4739 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.007653,   Adjusted R-squared:  0.006606 
## F-statistic: 7.309 on 5 and 4739 DF,  p-value: 7.836e-07
pred <- predictions(model_dwell,
                 by = c("gender", "year"))

# Test maximum contrast
max_comp <- avg_comparisons(model_dwell, 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.134     0.0589 -2.27   0.0233 5.4 -0.249 -0.0181
## 
## Type: response
# Express in SD units
round(max_comp$estimate/sd(combined$satis_dwell, na.rm = TRUE), 2)
## [1] -0.16
# Plot
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 housing satisfaction") +
  xlab("Survey year") +
  labs(color = "Gender")

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