Aim of this project

As part of the Multi100 project, I’m reanalyzing the data from Huijts et al. (2013), and investigate one specific hypothesis.

The hypothesis: The disadvantage in psychological well-being of childless people is smaller in countries with tolerant norms towards childlessness (p. 32)

Reference: Huijts, T., Kraaykamp, G., & Subramanian, S. V. (2013). Childlessness and psychological well-being in context: A multilevel study on 24 European countries. European Sociological Review, 29(1), 32–47. https://doi.org/10.1093/esr/jcr037

knitr::opts_chunk$set(
    echo = TRUE,
    message = FALSE,
    warning = FALSE
)

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0      ✔ purrr   1.0.0 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.5.0 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(here)

## here() starts at C:/Users/thoma/OneDrive/Documents/GitHub/collaboration/multi100_2023

library(readxl)
library(lmerTest)

## Loading required package: lme4
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## 
## Attaching package: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step

library(broom)
library(broom.mixed)
library(sjPlot)

theme_set(theme_light())

Reading the data

We made a few changes in the dataset before analysis.

Created a centered version of age (and age^2) to avoid structural multicollinearity.
Set the baseline levels for factor variables.
Found a small rounding error in a small number of cases in childless_disapp_std values for Slovenia that we corrected. In the original dataset, it seemed like as Slovenia having two slightly different values for this variable, although each country should have only one value.

edu_levels <- c("Primary", "Lower secondary", "Upper secondary", "Teritary", "Missing")

wb <- 
  read_csv(here("data/SRD Huijts data final.xls")) %>% 
  mutate(gender = recode(gndr, `1` = "Male", `2` = "Female"),
         # Create centered age (to avoid multicollinearity) and age^2
         age_ctd = scale(agea, scale = FALSE) %>% as.numeric(),
         age_ctd2 = age_ctd^2,
         # Set education baseline to Primary, set explicit NA
         edu = fct_explicit_na(edu, na_level = "Missing") %>% 
               fct_relevel(edu_levels),
         paredu = fct_relevel(paredu, edu_levels),
         # Set parental status baseline to living with children
         par_stat = fct_relevel(par_stat, "Living w children"),
         # Correcting a rounding error for country: SI
         childless_disapp_std = if_else(cntry == "SI", 0.1954625, childless_disapp_std)
         )

Finding the best fitting random structure

We compared the fixed model (multiple linear regression without random components), the random intercept model, and random intercept and slopes model.

This approach is recommended by Zuur et al. (2013).

Zuur, A. F., Ieno, E. N., Walker, N., Saveliev, A. A., & Smith, G. M. (2009). Mixed effects models and extensions in ecology with R. New York, NY: Springer New York. https://doi.org/10.1007/978-0-387-87458-6

# Non-multilevel model
wb_0 <-
  lm(wellbeing ~ age_ctd + age_ctd2 + married + edu + paredu + rel + pdwrk + 
       brncntr + sclmeet + par_stat * childless_disapp_std + gender, data = wb)

# Random intercept model
wb_int <-
  lmer(wellbeing ~ age_ctd + age_ctd2 + married + edu + paredu + rel + pdwrk + 
       brncntr + sclmeet + par_stat * childless_disapp_std + gender +
       (1 | cntry), data = wb)

# Random intercept and slopes model
wb_intslope <-
  lmer(wellbeing ~ age_ctd + age_ctd2 + married + edu + paredu + rel + pdwrk + 
       brncntr + sclmeet + par_stat * childless_disapp_std + gender +
       (1 + par_stat | cntry), data = wb)
     
anova(wb_int, wb_0, wb_intslope)

## Data: wb
## Models:
## wb_0: wellbeing ~ age_ctd + age_ctd2 + married + edu + paredu + rel + pdwrk + brncntr + sclmeet + par_stat * childless_disapp_std + gender
## wb_int: wellbeing ~ age_ctd + age_ctd2 + married + edu + paredu + rel + pdwrk + brncntr + sclmeet + par_stat * childless_disapp_std + gender + (1 | cntry)
## wb_intslope: wellbeing ~ age_ctd + age_ctd2 + married + edu + paredu + rel + pdwrk + brncntr + sclmeet + par_stat * childless_disapp_std + gender + (1 + par_stat | cntry)
##             npar   AIC   BIC logLik deviance   Chisq Df Pr(>Chisq)    
## wb_0          25 36413 36617 -18182    36363                          
## wb_int        26 35686 35897 -17817    35634 729.511  1  < 2.2e-16 ***
## wb_intslope   31 35671 35923 -17805    35609  24.963  5  0.0001416 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We used the BIC for comparison, as we have several predictors, and we wanted to control for complexity. The random intercept model is the best fit.

Fitting the model

We fit two separate LMEMs, one for males and one for females. The names of the predictors in the model output table is unfortunately not in order due to a bug. However, the most important parts are the interaction terms between par_stat and childless_disapp_std. We can see than the standardized betas for these interactions are all negative and significant.

wb_lmer <- 
  wb %>% 
  group_by(gender) %>% 
  nest() %>% 
  ungroup() %>% 
  mutate(model = map(data, 
                     ~lmer(
                       wellbeing ~ age_ctd + age_ctd2 + married + edu + paredu + 
                                   rel + pdwrk + brncntr + sclmeet + 
                                   par_stat * childless_disapp_std + 
                                   (1 | cntry),
                       data = .x)), 
         model_augment = map(model, augment)) %>% 
  arrange(gender)

wb_lmer %>% 
  pull(model) %>% 
  tab_model(
          dv.labels = pull(wb_lmer, gender),
          show.std = TRUE,
          show.est = FALSE,
          show.se = TRUE,
          show.reflvl = TRUE, 
          show.stat = TRUE,
          show.ci = FALSE,
          show.aic = TRUE,
          show.dev = TRUE,
          prefix.labels = "varname")

	Female						Male
Predictors	std. Beta	standardized std. Error	Statistic	std. Statistic	p	std. p	std. Beta	standardized std. Error	Statistic	std. Statistic	p	std. p
(Intercept)	-0.02	0.04	50.59	-0.45	<0.001	0.651	-0.09	0.05	54.83	-1.93	<0.001	0.054
age_ctd	0.00	0.01	0.28	0.28	0.776	0.776	0.07	0.01	5.53	5.53	<0.001	<0.001
age_ctd2	-0.00	0.01	-0.19	-0.19	0.852	0.852	-0.02	0.01	-1.76	-1.76	0.078	0.078
brncntr	0.02	0.01	2.83	2.83	0.005	0.005	0.04	0.01	4.76	4.76	<0.001	<0.001
childless_disapp_std	-0.21	0.04	-5.58	-5.58	<0.001	<0.001	-0.24	0.04	-6.36	-6.36	<0.001	<0.001
married: Married or cohabitating	Reference						Reference
married: Never married	-0.23	0.03	-7.04	-7.04	<0.001	<0.001	-0.25	0.04	-6.71	-6.71	<0.001	<0.001
married: Separated or divorced	-0.29	0.02	-12.32	-12.32	<0.001	<0.001	-0.33	0.03	-10.72	-10.72	<0.001	<0.001
married: Widowed	-0.42	0.02	-19.95	-19.95	<0.001	<0.001	-0.60	0.03	-17.22	-17.22	<0.001	<0.001
edu: Primary	Reference						Reference
edu: Lower secondary	0.12	0.03	4.47	4.47	<0.001	<0.001	0.10	0.03	3.15	3.15	0.002	0.002
edu: Upper secondary	0.25	0.03	9.40	9.40	<0.001	<0.001	0.22	0.03	6.92	6.92	<0.001	<0.001
edu: Teritary	0.33	0.03	11.76	11.76	<0.001	<0.001	0.27	0.03	8.08	8.08	<0.001	<0.001
edu: Missing	0.15	0.15	1.01	1.01	0.315	0.315	0.06	0.16	0.39	0.39	0.699	0.699
paredu: Primary	Reference						Reference
par_statChildless:childless_disapp_std	-0.12	0.03	-4.70	-4.70	<0.001	<0.001	-0.08	0.03	-2.84	-2.84	0.004	0.004
paredu: Lower secondary	0.02	0.02	0.83	0.83	0.404	0.404	0.01	0.03	0.48	0.48	0.630	0.630
par_statEmpty nest:childless_disapp_std	-0.09	0.02	-5.56	-5.56	<0.001	<0.001	-0.06	0.02	-2.99	-2.99	0.003	0.003
paredu: Upper secondary	0.04	0.03	1.51	1.51	0.131	0.131	0.05	0.03	1.89	1.89	0.058	0.058
paredu: Teritary	0.04	0.03	1.24	1.24	0.214	0.214	0.05	0.03	1.51	1.51	0.131	0.131
paredu: Missing	-0.07	0.03	-2.72	-2.72	0.006	0.006	-0.08	0.03	-2.50	-2.50	0.013	0.013
par_stat: Living w children	Reference						Reference
par_stat: Childless	-0.03	0.03	-0.99	-0.88	0.320	0.377	-0.13	0.03	-4.20	-3.91	<0.001	<0.001
par_stat: Empty nest	-0.03	0.02	-1.41	-1.28	0.159	0.200	0.04	0.02	1.21	1.58	0.224	0.114
pdwrk	0.08	0.01	8.39	8.39	<0.001	<0.001	0.16	0.01	13.43	13.43	<0.001	<0.001
rel	0.03	0.01	3.89	3.89	<0.001	<0.001	0.02	0.01	1.69	1.69	0.091	0.091
sclmeet	0.17	0.01	21.69	21.69	<0.001	<0.001	0.12	0.01	13.31	13.31	<0.001	<0.001
Random Effects
σ²	0.26						0.21
τ₀₀	0.01 _cntry						0.01 _cntry
ICC	0.04						0.04
N	24 _cntry						24 _cntry
Observations	14049						10960
Marginal R² / Conditional R²	0.213 / 0.243						0.192 / 0.224
Deviance	21131.589						14284.545
AIC	21358.378						14509.149

Visualization of the association between country well being, parental status, and childlessness disaproval

The interactions can be best seen in the following plot.

# Run separate linear regressions by country by gender
wb_lm <- 
  wb %>% 
  group_by(country = cntry, gender, childless_disapp_std) %>% 
  nest() %>% 
  mutate(model = map(data, 
                      ~lm(wellbeing ~ par_stat + age_ctd + age_ctd2 + married + 
                                      edu + paredu + rel + pdwrk + brncntr +
                                      sclmeet,
                                      data = .x)),
         model_tidy = map(model, tidy, conf.int = TRUE),
         n = map_int(data, nrow))

wb_lm %>% 
    unnest(model_tidy) %>% 
    ungroup() %>% 
    filter(str_detect(term, "par_stat")) %>% 
    mutate(term = str_remove(term, "par_stat")) %>% 
    select(childless_disapp_std:country, term, estimate, n) %>% 
    ggplot() +
    aes(x = childless_disapp_std, y = estimate, color = term, group = term) +
    geom_point(aes(size = n), alpha = .5) +
    # Fit wls regression
    geom_smooth(method = "lm", mapping = aes(weight = n), se = FALSE) + 
    geom_text(aes(label = country), check_overlap = TRUE, color = "black") +
    facet_wrap(~gender) +
    labs(title = "Association between effect sizes (beta) of parental status and well-being and childlessness disapproval by gender",
         subtitle = "Effects are in comparison with living with children. Lines show weighted regressions.",
         x = "Childlessness disapproval",
         y = "Effect of parental status on well-being",
         color = "Parental status",
         size = "Sample size")

Alternative visualization

In this plot all three groups are shown (but the interactions are less visible). The two plots become basically the same if we rotate the previous plot by the slope of “Living w children” group.

wb_lmer %>% 
  unnest(model_augment) %>% 
  group_by(gender, cntry, par_stat) %>% 
  summarise(estimate = mean(.fitted), 
            .groups = "drop") %>% 
  left_join(count(wb, 
                  cntry, childless_disapp_std),
            by = "cntry") %>% 
  ggplot() +
  aes(x = childless_disapp_std,
      y = estimate,
      color = par_stat,
      group = par_stat) +
    geom_point(aes(size = n), alpha = .5) +
    # Fit wls regression
    geom_smooth(method = "lm", mapping = aes(weight = n), se = FALSE) + 
    geom_text(aes(label = cntry), check_overlap = TRUE, color = "black") +
    facet_wrap(~gender) +
    labs(title = "Predicted wellbeing by gender, country, and the association with childlessness disapproval",
         subtitle = "Lines show weighted regressions.",
         x = "Childlessness disapproval",
         y = "Average prediction",
         color = "Parental status")

Single model including both genders, with a 3-way interaction.

The 3-way interactions show that there is no significant gender difference in the effect of childlessness disapproval on well-being for people with childless and empty nest parental status (as compared to those living with children). Note that a bug makes the order of predictors slightly confusing.

wb_lmer2 <- 
  lmer(wellbeing ~ age_ctd + age_ctd2 + married + edu + paredu + 
                   rel + pdwrk + brncntr + sclmeet + 
                   gender * par_stat * childless_disapp_std + 
                   (1 | cntry),
       data = wb)

wb_lmer2 %>% 
  tab_model(
          show.std = TRUE,
          show.est = FALSE,
          show.se = TRUE,
          show.stat = TRUE,
          show.reflvl = TRUE, 
          show.ci = FALSE,
          show.aic = TRUE,
          show.dev = TRUE,
          prefix.labels = "varname")

	wellbeing
Predictors	std. Beta	standardized std. Error	Statistic	std. Statistic	p	std. p
(Intercept)	-0.10	0.04	62.58	-2.52	<0.001	0.012
age_ctd	0.03	0.01	3.86	3.86	<0.001	<0.001
age_ctd2	-0.01	0.01	-1.59	-1.59	0.112	0.112
brncntr	0.03	0.01	5.08	5.08	<0.001	<0.001
childless_disapp_std	-0.21	0.04	-5.98	-5.98	<0.001	<0.001
edu: Primary	Reference
edu: Lower secondary	0.12	0.02	5.75	5.75	<0.001	<0.001
edu: Upper secondary	0.24	0.02	12.00	12.00	<0.001	<0.001
edu: Teritary	0.30	0.02	14.46	14.46	<0.001	<0.001
edu: Missing	0.12	0.11	1.13	1.13	0.257	0.257
genderMale	0.12	0.02	6.62	6.63	<0.001	<0.001
genderMale:childless_disapp_std	-0.01	0.02	-0.60	-0.60	0.552	0.552
genderMale:par_statChildless	-0.09	0.04	-2.52	-2.67	0.012	0.008
genderMale:par_statChildless:childless_disapp_std	0.04	0.04	1.14	1.14	0.253	0.253
genderMale:par_statEmpty nest	0.07	0.02	3.05	2.93	0.002	0.003
genderMale:par_statEmpty nest:childless_disapp_std	0.04	0.02	1.52	1.52	0.129	0.129
marriedNever married	-0.24	0.02	-9.94	-9.94	<0.001	<0.001
marriedSeparated or divorced	-0.31	0.02	-16.71	-16.71	<0.001	<0.001
marriedWidowed	-0.47	0.02	-26.83	-26.83	<0.001	<0.001
paredu: Primary	Reference
par_statChildless:childless_disapp_std	-0.13	0.03	-5.05	-5.05	<0.001	<0.001
paredu: Lower secondary	0.01	0.02	0.76	0.76	0.447	0.447
par_statEmpty nest:childless_disapp_std	-0.09	0.02	-5.72	-5.72	<0.001	<0.001
paredu: Upper secondary	0.04	0.02	2.19	2.19	0.029	0.029
paredu: Teritary	0.04	0.02	1.69	1.69	0.091	0.091
paredu: Missing	-0.08	0.02	-3.79	-3.79	<0.001	<0.001
par_stat: Living w children	Reference
par_stat: Childless	-0.03	0.03	-1.45	-1.08	0.146	0.281
par_stat: Empty nest	-0.03	0.02	-2.22	-1.80	0.027	0.071
pdwrk	0.11	0.01	15.11	15.11	<0.001	<0.001
rel	0.02	0.01	4.09	4.09	<0.001	<0.001
sclmeet	0.15	0.01	25.41	25.41	<0.001	<0.001
Random Effects
σ²	0.24
τ₀₀ _cntry	0.01
ICC	0.04
N _cntry	24
Observations	25009
Marginal R² / Conditional R²	0.220 / 0.248
Deviance	35602.697
AIC	35897.039

Summary of the results

Childlessness disapproval in a country is negatively associated with the well-being.
For childless people, the effect of childlessness disapproval is more negative than for those who have children.
This is also true for those people whose children moved out (i.e., empty nest) but to a lesser - but still statistically significant - extent.
These associations are of similar magnitude for males and females, with betas between -0.12 and -0.06 (SEs = 0.02-0.03), while controlling for several variables.
Thus it is true that in countries which are more tolerant towards childlessness, the negative psychological effect of childlessness is smaller.

Conclusion

We found support for the investigated hypothesis ( The disadvantage in psychological well-being of childless people is smaller in countries with tolerant norms towards childlessness (p. 32) ).

Multi100 reanalysis of ‘Childlessness and Psychological Well-Being in Context: A Multilevel Study on 24 European Countries’

2023-02-28