Add option to pull coef() instead of ranef() for random effects

[DominiqueMakowski]

Initially I thought this should fit in modelbased, but I realize that this is a relatively low-level thing which makes more sense to be supported here.

Basically, for all mixed models, we can extract the random effects as ranef() (which is what parameters() does) or coef() (which corresponds to random + fixed).

I would add an argument to model_parameters random="deviation" (random effects are expressed as deviation from the fixed effect) vs. random="total" (it corresponds to the absolute effect).

Fro brms, we could look at the source code of coef.brmsfit to see how it does the transformation

  data <- iris
  data$Group <- as.factor(rep(c("G1", "G2", "G3"), each = 50))

  m1 <- glmmTMB::glmmTMB(Sepal.Width ~ Petal.Width + (Petal.Width | Group), data = data)
  lme4::ranef(m1)
#> $Group
#>    (Intercept)  Petal.Width
#> G1   1.0758012 -0.005194106
#> G2  -0.4001235  0.014754929
#> G3  -0.6756778 -0.009560821
  coef(m1)
#> $Group
#>    (Intercept) Petal.Width
#> G1    3.243472   0.7369045
#> G2    1.767547   0.7568536
#> G3    1.491993   0.7325378
  
  parameters::model_parameters(m1, effects = "random", group_level = TRUE)
#> # Random Effects
#> 
#> Parameter        | Coefficient |   SE |        95% CI
#> -----------------------------------------------------
#> (Intercept) [G1] |        1.08 | 0.47 | [ 0.15, 2.00]
#> (Intercept) [G2] |       -0.40 | 0.57 | [-1.51, 0.71]
#> (Intercept) [G3] |       -0.68 | 0.55 | [-1.75, 0.40]
#> Petal.Width [G1] |   -5.19e-03 | 0.29 | [-0.58, 0.57]
#> Petal.Width [G2] |        0.01 | 0.29 | [-0.56, 0.59]
#> Petal.Width [G3] |   -9.56e-03 | 0.21 | [-0.42, 0.40]
#> 
  
  
  m3 <- brms::brm(Sepal.Width ~ Petal.Width + (Petal.Width | Group), data = data, refresh=0)
  lme4::ranef(m3)
#> $Group
#> , , Intercept
#> 
#>      Estimate Est.Error       Q2.5    Q97.5
#> G1  1.1400671 0.8776806 -0.5852008 3.093138
#> G2 -0.5236839 0.8781311 -2.3384922 1.340037
#> G3 -0.4704580 0.8858766 -2.2499332 1.474791
#> 
#> , , Petal.Width
#> 
#>        Estimate Est.Error       Q2.5     Q97.5
#> G1  0.009368953 0.3797845 -0.8376778 0.8034064
#> G2  0.104031171 0.3383296 -0.6117355 0.8604620
#> G3 -0.136885669 0.3261730 -0.9531138 0.4863215
  coef(m3)
#> $Group
#> , , Intercept
#> 
#>    Estimate  Est.Error      Q2.5    Q97.5
#> G1 3.225730 0.08978179 3.0362333 3.398977
#> G2 1.561979 0.27815857 0.9840993 2.059511
#> G3 1.615205 0.30416862 1.0410086 2.220016
#> 
#> , , Petal.Width
#> 
#>     Estimate Est.Error      Q2.5     Q97.5
#> G1 0.8172649 0.3210640 0.1884875 1.4853629
#> G2 0.9119271 0.2069303 0.5402409 1.3365941
#> G3 0.6710102 0.1486471 0.3783872 0.9519822
  
  parameters::model_parameters(m3, effects = "random", group_level = TRUE)
#> # Fixed Effects
#> 
#> Parameter               |   Median |        95% CI |     pd |  Rhat |     ESS
#> -----------------------------------------------------------------------------
#> Group.G1                |     1.11 | [-0.59, 3.09] | 92.47% | 1.002 |  794.00
#> Group.G2                |    -0.52 | [-2.34, 1.34] | 76.92% | 1.001 |  919.00
#> Group.G3                |    -0.49 | [-2.25, 1.47] | 74.60% | 1.001 |  869.00
#> Group.G1                | 5.29e-03 | [-0.84, 0.80] | 51.73% | 1.002 |  752.00
#> Group.G2                |     0.06 | [-0.61, 0.86] | 63.62% | 1.001 |  634.00
#> Group.G3                |    -0.09 | [-0.95, 0.49] | 68.88% | 1.003 |  457.00
#> (Intercept)             |     1.40 | [ 0.56, 4.06] |   100% | 1.002 | 1072.00
#> Petal.Width             |     0.34 | [ 0.02, 1.62] |   100% | 1.004 |  807.00
#> Intercept ~ Petal.Width |    -0.06 | [-0.95, 0.93] | 53.02% | 1.003 | 2342.00
#> 
#> Uncertainty intervals (equal-tailed) computed using a MCMC distribution
#>   approximation.

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

[bwiernik]

I would make this another option for effects, "total" or similar. I don't like "deviation" because models can include random effects without corresponding fixed components.

Is this what estimate_grouplevel() does?

[DominiqueMakowski]

because models can include random effects without corresponding fixed components

Good point! we could have "random" instead of "deviation"

Is this what estimate_grouplevel() does?

estimate_grouplevel pulls and filters the parameters from model_parameters(), and then "simply" does a bit of cleaning and formatting (in particular, it adds a couple of columns to separate the group, the levels, the components, the effects etc. (+ bespoke plotting

[strengejacke]

Possible duplicate to #484

[DominiqueMakowski]

The magic probably happens here: https://github.com/paul-buerkner/brms/blob/2800c0823fe466fe86566b958ba3f45404c6c505/R/brmsfit-methods.R#L241-L259

[DominiqueMakowski]

Might be tedious to duplicate the implementation (requires to find all ranefs and add their corresponding fixef), might be easier to directly call coef(, ...summary=FALSE) and get the correct draws

[strengejacke]

we could have "random" instead of "deviation"

We already have that option for the effects argument. If we want to use the effects argument, I'd say "total" would probably the best option?

[DominiqueMakowski]

I'm not super sure about having that as an option of effect, given that this argument is mostly used for filtering effects, and the option discussed here changes the values of random effects...
Maybe another argument like random=c("effects", "total")? Althoug there we end up with a ridiculous parameters(..., effects="random", random="effects") 😅

[strengejacke]

option discussed here changes the values of random effects...

Unless you see this as not changing values of random effects, but rather returning other effects than isolated fixed or isolated random.

[DominiqueMakowski]

fair, then maybe effects="random_total" ?

[strengejacke]

Here's an example that should work, which probably covers one of the more complex cases. I have an idea in mind how to implement this...

data("fish", package = "insight")

m <- suppressWarnings(glmmTMB::glmmTMB(
  count ~ child + camper + (1 + xb | persons) + (1 | ID),
  ziformula = ~ child + livebait + (1 + zg + nofish | ID),
  dispformula = ~xb,
  data = fish,
  family = glmmTMB::truncated_poisson()
))

coef(m)
#> $cond
#> $persons
#>         xb (Intercept)     child   camper1
#> 1 1.393807 -1.32143624 -1.134178 0.2968573
#> 2 1.185343 -0.78802935 -1.134178 0.2968573
#> 3 1.400908 -1.33961436 -1.134178 0.2968573
#> 4 0.892148 -0.03771508 -1.134178 0.2968573
#> 
#> $ID
#>   xb (Intercept)     child   camper1
#> 1  0    2.039620 -1.134178 0.2968573
#> 2  0    2.383518 -1.134178 0.2968573
#> 3  0    2.338501 -1.134178 0.2968573
#> 4  0    2.233461 -1.134178 0.2968573
#> 
#> attr(,"class")
#> [1] "coef.mer"
#> 
#> $zi
#> $ID
#>           zg     nofish1 (Intercept)      child  livebait1
#> 1 -1.0098205  0.42640593    1.761620 -0.2135642 -0.8754563
#> 2 -0.6800334 -0.08774915    1.427959 -0.2135642 -0.8754563
#> 3 -0.4952986  0.13642795    1.238911 -0.2135642 -0.8754563
#> 4 -0.7351137  0.19240968    1.482873 -0.2135642 -0.8754563
#> 
#> attr(,"class")
#> [1] "coef.mer"
#> 
#> $disp
#> list()

^{Created on 2025-03-02 with reprex v2.1.1}

[strengejacke]

fair, then maybe effects="random_total" ?

For coef.glmmTMB: a similar list, but containing the overall coefficient value for each level, i.e., the sum of the fixed effect estimate and the random effect value for that level.

Or "overall"?

[bwiernik]

I like "total", "group_total", or "group_specific" as names for this option of "effects"

[DominiqueMakowski]

effects="random_raw" (= effects="random", default) vs. effects="random_total"?