estimate_contrasts() for estimate_relation() etc

[strengejacke]

When we want contrasts for random effects, marginaleffects returns too small standard errors (at least for glmmTMB), while estimate_relation() relies on the SEs from predict() (which are the best we can get, according to the package authors).

Thus, we implement a very basic contrast method, which, however works for random effects, and thereby allows MAIHDA analysis with the modelbased package.

[DominiqueMakowski]

MAIDHA stands for?

[strengejacke]

https://www.sciencedirect.com/science/article/pii/S235282732400065X

Latest shit in epidemiology

[strengejacke]

See also https://strengejacke.github.io/ggeffects/articles/practical_intersectionality.html

[strengejacke]

See also https://strengejacke.github.io/ggeffects/articles/practical_intersectionality.html

You can now fully replace ggeffects with modelbased in that vignette. 🥂

[strengejacke]

Are you ok with this "feature"? (another hidden one actually, but... else we don't get uncertainties for contrasts of random effects, neither for glmmTMB nor lme4:

library(modelbased)
data(efc, package = "modelbased")

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

# Quality of Life score ranges from 0 to 25
m_null <- glmmTMB::glmmTMB(qol ~ 1 + (1 | gender:employed:age), data = d)
out1 <- estimate_relation(m_null, by = c("gender", "employed", "age"))
out2 <- estimate_means(m_null, by = c("gender", "employed", "age"))

out1
#> Model-based Predictions
#> 
#> gender | employed | age   | Predicted |   SE |         95% CI | Residuals
#> -------------------------------------------------------------------------
#> Male   | no       | -40   |     15.24 | 0.88 | [13.52, 16.96] |     -0.87
#> Female | no       | -40   |     15.17 | 0.69 | [13.83, 16.52] |     -0.80
#> Male   | yes      | -40   |     16.12 | 0.79 | [14.58, 17.66] |     -1.75
#> Female | yes      | -40   |     16.04 | 0.61 | [14.84, 17.23] |     -1.67
#> Male   | no       | 41-64 |     14.94 | 0.66 | [13.65, 16.23] |     -0.57
#> Female | no       | 41-64 |     13.75 | 0.33 | [13.11, 14.40] |      0.61
#> Male   | yes      | 41-64 |     15.47 | 0.56 | [14.38, 16.56] |     -1.10
#> Female | yes      | 41-64 |     14.10 | 0.35 | [13.42, 14.77] |      0.27
#> Male   | no       | 65+   |     14.41 | 0.63 | [13.18, 15.65] |     -0.05
#> Female | no       | 65+   |     13.53 | 0.44 | [12.68, 14.39] |      0.83
#> Male   | yes      | 65+   |     15.34 | 1.09 | [13.22, 17.47] |     -0.97
#> Female | yes      | 65+   |     14.85 | 1.04 | [12.82, 16.89] |     -0.49
#> 
#> Variable predicted: qol
#> Predictors modulated: gender, employed, age
out2
#> Estimated Marginal Means
#> 
#> gender | employed | age   |  Mean |   SE |         95% CI |     z
#> -----------------------------------------------------------------
#> Male   | no       | -40   | 15.24 | 0.40 | [14.46, 16.02] | 38.23
#> Female | no       | -40   | 15.17 | 0.40 | [14.39, 15.95] | 38.06
#> Male   | yes      | -40   | 16.12 | 0.40 | [15.34, 16.90] | 40.43
#> Female | yes      | -40   | 16.04 | 0.40 | [15.25, 16.82] | 40.22
#> Male   | no       | 41-64 | 14.94 | 0.40 | [14.16, 15.72] | 37.47
#> Female | no       | 41-64 | 13.75 | 0.40 | [12.97, 14.54] | 34.50
#> Male   | yes      | 41-64 | 15.47 | 0.40 | [14.69, 16.25] | 38.80
#> Female | yes      | 41-64 | 14.10 | 0.40 | [13.31, 14.88] | 35.36
#> Male   | no       | 65+   | 14.41 | 0.40 | [13.63, 15.19] | 36.15
#> Female | no       | 65+   | 13.53 | 0.40 | [12.75, 14.31] | 33.95
#> Male   | yes      | 65+   | 15.34 | 0.40 | [14.56, 16.12] | 38.48
#> Female | yes      | 65+   | 14.85 | 0.40 | [14.07, 15.63] | 37.26
#> 
#> Variable predicted: qol
#> Predictors modulated: gender, employed, age

estimate_contrasts(out1, "gender", by = c("employed", "age"))
#> Model-based Contrasts Analysis
#> 
#> Level1 | Level2 | employed | age   | Difference |   SE |        95% CI
#> ----------------------------------------------------------------------
#> Male   | Female | no       | -40   |       0.07 | 1.11 | [-2.11, 2.25]
#> Male   | Female | yes      | -40   |       0.08 | 0.99 | [-1.86, 2.03]
#> Male   | Female | no       | 41-64 |       1.18 | 0.74 | [-0.26, 2.63]
#> Male   | Female | yes      | 41-64 |       1.37 | 0.66 | [ 0.09, 2.66]
#> Male   | Female | no       | 65+   |       0.88 | 0.77 | [-0.62, 2.38]
#> Male   | Female | yes      | 65+   |       0.49 | 1.50 | [-2.45, 3.43]
#> 
#> Level1 | Statistic |     p
#> --------------------------
#> Male   |      0.06 | 0.951
#> Male   |      0.08 | 0.932
#> Male   |      1.60 | 0.109
#> Male   |      2.09 | 0.036
#> Male   |      1.15 | 0.250
#> Male   |      0.33 | 0.745
#> 
#> Variable predicted: qol
#> Predictors contrasted: gender
estimate_contrasts(m_null, "gender", by = c("employed", "age"))
#> Level1 | Level2 | employed | age   | Difference
#> -----------------------------------------------
#> Female | Male   | no       | -40   |      -0.07
#> Female | Male   | yes      | -40   |      -0.08
#> Female | Male   | no       | 41-64 |      -1.18
#> Female | Male   | yes      | 41-64 |      -1.37
#> Female | Male   | no       | 65+   |      -0.88
#> Female | Male   | yes      | 65+   |      -0.49

^{Created on 2025-01-30 with reprex v2.1.1}

library(modelbased)
data(efc, package = "modelbased")

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

# Quality of Life score ranges from 0 to 25
m_null <- lme4::lmer(qol ~ 1 + (1 | gender:employed:age), data = d)
out1 <- estimate_relation(m_null, by = c("gender", "employed", "age"))
out2 <- estimate_means(m_null, by = c("gender", "employed", "age"))

out1
#> Model-based Predictions
#> 
#> gender | employed | age   | Predicted |   SE |         95% CI | Residuals
#> -------------------------------------------------------------------------
#> Male   | no       | -40   |     15.28 | 0.41 | [14.48, 16.08] |     -0.91
#> Female | no       | -40   |     15.19 | 0.41 | [14.39, 15.99] |     -0.83
#> Male   | yes      | -40   |     16.20 | 0.41 | [15.40, 17.00] |     -1.83
#> Female | yes      | -40   |     16.08 | 0.41 | [15.28, 16.88] |     -1.72
#> Male   | no       | 41-64 |     14.94 | 0.41 | [14.15, 15.74] |     -0.58
#> Female | no       | 41-64 |     13.74 | 0.41 | [12.94, 14.54] |      0.63
#> Male   | yes      | 41-64 |     15.49 | 0.41 | [14.69, 16.29] |     -1.13
#> Female | yes      | 41-64 |     14.09 | 0.41 | [13.29, 14.89] |      0.28
#> Male   | no       | 65+   |     14.40 | 0.41 | [13.60, 15.20] |     -0.03
#> Female | no       | 65+   |     13.51 | 0.41 | [12.71, 14.31] |      0.86
#> Male   | yes      | 65+   |     15.42 | 0.41 | [14.62, 16.21] |     -1.05
#> Female | yes      | 65+   |     14.86 | 0.41 | [14.06, 15.66] |     -0.49
#> 
#> Variable predicted: qol
#> Predictors modulated: gender, employed, age
out2
#> Estimated Marginal Means
#> 
#> gender | employed | age   |  Mean |   SE |         95% CI | t(892)
#> ------------------------------------------------------------------
#> Male   | no       | -40   | 15.28 | 0.41 | [14.48, 16.08] |  37.49
#> Female | no       | -40   | 15.19 | 0.41 | [14.39, 15.99] |  37.28
#> Male   | yes      | -40   | 16.20 | 0.41 | [15.40, 17.00] |  39.74
#> Female | yes      | -40   | 16.08 | 0.41 | [15.28, 16.88] |  39.46
#> Male   | no       | 41-64 | 14.94 | 0.41 | [14.15, 15.74] |  36.67
#> Female | no       | 41-64 | 13.74 | 0.41 | [12.94, 14.54] |  33.71
#> Male   | yes      | 41-64 | 15.49 | 0.41 | [14.69, 16.29] |  38.01
#> Female | yes      | 41-64 | 14.09 | 0.41 | [13.29, 14.89] |  34.56
#> Male   | no       | 65+   | 14.40 | 0.41 | [13.60, 15.20] |  35.32
#> Female | no       | 65+   | 13.51 | 0.41 | [12.71, 14.31] |  33.14
#> Male   | yes      | 65+   | 15.42 | 0.41 | [14.62, 16.21] |  37.82
#> Female | yes      | 65+   | 14.86 | 0.41 | [14.06, 15.66] |  36.46
#> 
#> Variable predicted: qol
#> Predictors modulated: gender, employed, age

estimate_contrasts(out1, "gender", by = c("employed", "age"))
#> Model-based Contrasts Analysis
#> 
#> Level1 | Level2 | employed | age   | Difference |   SE |        95% CI
#> ----------------------------------------------------------------------
#> Male   | Female | no       | -40   |       0.09 | 0.58 | [-1.04, 1.22]
#> Male   | Female | yes      | -40   |       0.11 | 0.58 | [-1.02, 1.25]
#> Male   | Female | no       | 41-64 |       1.20 | 0.58 | [ 0.07, 2.34]
#> Male   | Female | yes      | 41-64 |       1.41 | 0.58 | [ 0.28, 2.54]
#> Male   | Female | no       | 65+   |       0.89 | 0.58 | [-0.24, 2.02]
#> Male   | Female | yes      | 65+   |       0.55 | 0.58 | [-0.58, 1.68]
#> 
#> Level1 | Statistic |     p
#> --------------------------
#> Male   |      0.15 | 0.879
#> Male   |      0.20 | 0.843
#> Male   |      2.09 | 0.037
#> Male   |      2.44 | 0.015
#> Male   |      1.54 | 0.123
#> Male   |      0.96 | 0.338
#> 
#> Variable predicted: qol
#> Predictors contrasted: gender
estimate_contrasts(m_null, "gender", by = c("employed", "age")) |> as.data.frame()
#>   Level1 Level2 employed   age  Difference SE CI_low CI_high  t  df  p
#> 1 Female   Male       no   -40 -0.08757017 NA     NA      NA NA 892 NA
#> 2 Female   Male       no 41-64 -1.20390797 NA     NA      NA NA 892 NA
#> 3 Female   Male       no   65+ -0.88945903 NA     NA      NA NA 892 NA
#> 4 Female   Male      yes   -40 -0.11429974 NA     NA      NA NA 892 NA
#> 5 Female   Male      yes 41-64 -1.40747801 NA     NA      NA NA 892 NA
#> 6 Female   Male      yes   65+ -0.55253657 NA     NA      NA NA 892 NA

^{Created on 2025-01-30 with reprex v2.1.1}