Example for building intersectional strata

intersectional.R

easystats::install_latest(force= TRUE)
library(datawizard)  # data wrangling and preparation
library(parameters)  # model summaries
library(glmmTMB)     # multilevel modelling

# sample data set
data(efc, package = "ggeffects")

efc <- efc |>
  # numeric to factors, set labels as levels
  to_factor(select = c("c161sex", "c172code", "c175empl")) |>
  # recode age into three groups
  recode_values(
    select = "c160age",
    recode = list(`1` = "min:40", `2` = 41:64, `3` = "65:max")
  ) |>
  # rename variables
  data_rename(
    pattern = c("c161sex", "c160age", "quol_5", "c175empl"),
    replacement = c("gender", "age", "qol", "employed")
  ) |>
  # age into factor, set levels, and change labels for education
  data_modify(age = factor(age, labels = c("-40", "41-64", "65+")))

efc$strata <- ifelse(
  is.na(efc$employed) | is.na(efc$gender) | is.na(efc$age),
  NA_character_,
  paste0(efc$gender, ", ", efc$employed, ", ", efc$age)
  )
efc$strata <- factor(efc$strata)

data_tabulate(efc$strata)
#> efc$strata <categorical>
#> # total N=908 valid N=900
#> 
#> Value              |   N | Raw % | Valid % | Cumulative %
#> -------------------+-----+-------+---------+-------------
#> Female, no, -40    |  37 |  4.07 |    4.11 |         4.11
#> Female, no, 41-64  | 238 | 26.21 |   26.44 |        30.56
#> Female, no, 65+    | 135 | 14.87 |   15.00 |        45.56
#> Female, yes, -40   |  63 |  6.94 |    7.00 |        52.56
#> Female, yes, 41-64 | 210 | 23.13 |   23.33 |        75.89
#> Female, yes, 65+   |   3 |  0.33 |    0.33 |        76.22
#> Male, no, -40      |  15 |  1.65 |    1.67 |        77.89
#> Male, no, 41-64    |  42 |  4.63 |    4.67 |        82.56
#> Male, no, 65+      |  50 |  5.51 |    5.56 |        88.11
#> Male, yes, -40     |  34 |  3.74 |    3.78 |        91.89
#> Male, yes, 41-64   |  70 |  7.71 |    7.78 |        99.67
#> Male, yes, 65+     |   3 |  0.33 |    0.33 |       100.00
#> <NA>               |   8 |  0.88 |    <NA> |         <NA>

m1 <- glmmTMB(qol ~ 1 + (1 | strata), data = efc)
m2 <- glmmTMB(qol ~ 1 + (1 | gender:employed:age), data = efc)

model_parameters(m1)
#> # Fixed Effects
#> 
#> Parameter   | Coefficient |   SE |         95% CI |     z |      p
#> ------------------------------------------------------------------
#> (Intercept) |       14.91 | 0.40 | [14.13, 15.70] | 37.41 | < .001
#> 
#> # Random Effects
#> 
#> Parameter              | Coefficient |       95% CI
#> ---------------------------------------------------
#> SD (Intercept: strata) |        1.03 | [0.56, 1.89]
#> SD (Residual)          |        5.23 | [4.99, 5.48]
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald z-distribution approximation.
model_parameters(m2)
#> # Fixed Effects
#> 
#> Parameter   | Coefficient |   SE |         95% CI |     z |      p
#> ------------------------------------------------------------------
#> (Intercept) |       14.91 | 0.40 | [14.13, 15.70] | 37.41 | < .001
#> 
#> # Random Effects
#> 
#> Parameter                           | Coefficient |       95% CI
#> ----------------------------------------------------------------
#> SD (Intercept: gender:employed:age) |        1.03 | [0.56, 1.89]
#> SD (Residual)                       |        5.23 | [4.99, 5.48]
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald z-distribution approximation.