Skip to content

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

Merged
strengejacke merged 31 commits intomainfrom
strengejacke/issue1067
Mar 3, 2025
Merged

Add option to pull coef() instead of ranef() for random effects#1068
strengejacke merged 31 commits intomainfrom
strengejacke/issue1067

Conversation

@strengejacke
Copy link
Member

@strengejacke strengejacke commented Mar 2, 2025

Fixes #1067
Fixes #484

@strengejacke
Copy link
Member Author

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

m1 <- suppressWarnings(glmmTMB::glmmTMB(
  count ~ child + camper + (1 | ID),
  data = fish,
  family = poisson()
))

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

m3 <- suppressWarnings(glmmTMB::glmmTMB(
  count ~ child + camper,
  ziformula = ~ child + livebait + (1 | ID),
  data = fish,
  family = glmmTMB::truncated_poisson()
))

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

m5 <- suppressWarnings(lme4::glmer(
  count ~ child + camper + (1 | ID),
  data = fish,
  family = poisson()
))

m6 <- suppressWarnings(lme4::lmer(
  Reaction ~ Days + (1 + Days | Subject),
  data = lme4::sleepstudy
))


parameters::model_parameters(m1, effects = "total")
#> Group | Level |   Parameter | Coefficient |   Component | Effects
#> -----------------------------------------------------------------
#> ID    |     1 | (Intercept) |        0.22 | conditional |   total
#> ID    |     2 | (Intercept) |        0.77 | conditional |   total
#> ID    |     3 | (Intercept) |        0.86 | conditional |   total
#> ID    |     4 | (Intercept) |        1.36 | conditional |   total

parameters::model_parameters(m2, effects = "total")
#> Group   | Level |   Parameter | Coefficient |     Component | Effects
#> ---------------------------------------------------------------------
#> persons |     1 | (Intercept) |       -1.20 |   conditional |   total
#> persons |     1 |          xb |        1.43 |   conditional |   total
#> persons |     2 | (Intercept) |       -0.88 |   conditional |   total
#> persons |     2 |          xb |        1.30 |   conditional |   total
#> persons |     3 | (Intercept) |       -0.88 |   conditional |   total
#> persons |     3 |          xb |        1.30 |   conditional |   total
#> persons |     4 | (Intercept) |        0.23 |   conditional |   total
#> persons |     4 |          xb |        0.87 |   conditional |   total
#> ID      |     1 | (Intercept) |        1.94 |   conditional |   total
#> ID      |     1 |          zg |        0.20 |   conditional |   total
#> ID      |     2 | (Intercept) |        2.06 |   conditional |   total
#> ID      |     2 |          zg |        0.15 |   conditional |   total
#> ID      |     3 | (Intercept) |        2.35 |   conditional |   total
#> ID      |     3 |          zg |        0.04 |   conditional |   total
#> ID      |     4 | (Intercept) |        1.52 |   conditional |   total
#> ID      |     4 |          zg |        0.36 |   conditional |   total
#> ID      |     1 | (Intercept) |        1.92 | zero_inflated |   total
#> ID      |     1 |          zg |       -0.97 | zero_inflated |   total
#> ID      |     1 |     nofish1 |    2.59e-18 | zero_inflated |   total
#> ID      |     2 | (Intercept) |        1.58 | zero_inflated |   total
#> ID      |     2 |          zg |       -0.67 | zero_inflated |   total
#> ID      |     2 |     nofish1 |    1.77e-18 | zero_inflated |   total
#> ID      |     3 | (Intercept) |        1.38 | zero_inflated |   total
#> ID      |     3 |          zg |       -0.50 | zero_inflated |   total
#> ID      |     3 |     nofish1 |    1.33e-18 | zero_inflated |   total
#> ID      |     4 | (Intercept) |        1.63 | zero_inflated |   total
#> ID      |     4 |          zg |       -0.71 | zero_inflated |   total
#> ID      |     4 |     nofish1 |    1.86e-18 | zero_inflated |   total

parameters::model_parameters(m3, effects = "total")
#> Group | Level |   Parameter | Coefficient |     Component | Effects
#> -------------------------------------------------------------------
#> ID    |     1 | (Intercept) |        0.30 | zero_inflated |   total
#> ID    |     2 | (Intercept) |        0.30 | zero_inflated |   total
#> ID    |     3 | (Intercept) |        0.30 | zero_inflated |   total
#> ID    |     4 | (Intercept) |        0.30 | zero_inflated |   total

parameters::model_parameters(m4, effects = "total")
#> Group   | Level |   Parameter | Coefficient |   Component | Effects
#> -------------------------------------------------------------------
#> persons |     1 | (Intercept) |       -1.23 | conditional |   total
#> persons |     1 |          xb |        1.35 | conditional |   total
#> persons |     2 | (Intercept) |       -0.87 | conditional |   total
#> persons |     2 |          xb |        1.22 | conditional |   total
#> persons |     3 | (Intercept) |       -1.11 | conditional |   total
#> persons |     3 |          xb |        1.30 | conditional |   total
#> persons |     4 | (Intercept) |        0.10 | conditional |   total
#> persons |     4 |          xb |        0.87 | conditional |   total

parameters::model_parameters(m5, effects = "total")
#> Group | Level |   Parameter | Coefficient | Effects
#> ---------------------------------------------------
#> ID    |     1 | (Intercept) |        0.22 |   total
#> ID    |     2 | (Intercept) |        0.77 |   total
#> ID    |     3 | (Intercept) |        0.86 |   total
#> ID    |     4 | (Intercept) |        1.36 |   total

parameters::model_parameters(m6, effects = "total")
#> Group   | Level |   Parameter | Coefficient | Effects
#> -----------------------------------------------------
#> Subject |   308 | (Intercept) |      253.66 |   total
#> Subject |   308 |        Days |       19.67 |   total
#> Subject |   309 | (Intercept) |      211.01 |   total
#> Subject |   309 |        Days |        1.85 |   total
#> Subject |   310 | (Intercept) |      212.44 |   total
#> Subject |   310 |        Days |        5.02 |   total
#> Subject |   330 | (Intercept) |      275.10 |   total
#> Subject |   330 |        Days |        5.65 |   total
#> Subject |   331 | (Intercept) |      273.67 |   total
#> Subject |   331 |        Days |        7.40 |   total
#> Subject |   332 | (Intercept) |      260.44 |   total
#> Subject |   332 |        Days |       10.20 |   total
#> Subject |   333 | (Intercept) |      268.25 |   total
#> Subject |   333 |        Days |       10.24 |   total
#> Subject |   334 | (Intercept) |      244.17 |   total
#> Subject |   334 |        Days |       11.54 |   total
#> Subject |   335 | (Intercept) |      251.07 |   total
#> Subject |   335 |        Days |       -0.28 |   total
#> Subject |   337 | (Intercept) |      286.30 |   total
#> Subject |   337 |        Days |       19.10 |   total
#> Subject |   349 | (Intercept) |      226.19 |   total
#> Subject |   349 |        Days |       11.64 |   total
#> Subject |   350 | (Intercept) |      238.34 |   total
#> Subject |   350 |        Days |       17.08 |   total
#> Subject |   351 | (Intercept) |      255.98 |   total
#> Subject |   351 |        Days |        7.45 |   total
#> Subject |   352 | (Intercept) |      272.27 |   total
#> Subject |   352 |        Days |       14.00 |   total
#> Subject |   369 | (Intercept) |      254.68 |   total
#> Subject |   369 |        Days |       11.34 |   total
#> Subject |   370 | (Intercept) |      225.79 |   total
#> Subject |   370 |        Days |       15.29 |   total
#> Subject |   371 | (Intercept) |      252.21 |   total
#> Subject |   371 |        Days |        9.48 |   total
#> Subject |   372 | (Intercept) |      263.72 |   total
#> Subject |   372 |        Days |       11.75 |   total

Created on 2025-03-02 with reprex v2.1.1

@DominiqueMakowski
Copy link
Member

I'd call it "random_total" because it still filters out fixed effects (and it's something that applies to random effects)

@strengejacke
Copy link
Member Author

Yeah, the name is not final yet. But what about the general output?

@strengejacke
Copy link
Member Author

And we need methods for rstanarm and brms

@DominiqueMakowski
Copy link
Member

DominiqueMakowski commented Mar 2, 2025

We need to also keep the SE/CI for lme4 and other models, it is a glmmTMB particularity to not yet computing them i think

@strengejacke
Copy link
Member Author

I don't think you can get SE/CI for the coef()?

@DominiqueMakowski
Copy link
Member

DominiqueMakowski commented Mar 3, 2025

It seems you're right, somehow I (mis)remembered that it was possible for lme4

(Unless with bootstrapping... that could be a longer-term feature)

@DominiqueMakowski
Copy link
Member

Note: In https://hypatia.math.ethz.ch/pipermail/r-help/2007-November/145442.html, Doug Bates calles the "default" random effects (from ranef()) "conditional"

So we could have:

  • effects="random" (which is technically "random_conditional")
  • effects="random_total" OR what about effects="random_unbiased" (to link it to the BLUP terminology)

@bwiernik
Copy link
Contributor

bwiernik commented Mar 3, 2025

I'm not a fan of "random_total". I think that's kind of confusing. If something other than "total", my preference is "group_total".

@strengejacke
Copy link
Member Author

strengejacke commented Mar 3, 2025

brms works (at least this complex example, I bet there are dozens of not yet considered exceptions)

library(brms)
library(parameters)

x <- insight::download_model("brms_zi_4")

### Model

#>  Family: zero_inflated_poisson 
#>   Links: mu = log; zi = logit 
#> Formula: count ~ child + camper + (1 + xb | persons) + (1 + zg | ID) 
#>          zi ~ child + livebait + (1 + zg + nofish | ID)
#>    Data: fish (Number of observations: 250) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000

### Here are the over all coefficients

coef(x)
#> $ID
#> , , Intercept
#> 
#>     Estimate Est.Error      Q2.5    Q97.5
#> 1 -0.7872119  1.042311 -2.700848 1.478536
#> 2 -0.6960301  1.036887 -2.593541 1.606411
#> 3 -0.3863236  1.028798 -2.281176 1.859450
#> 4 -1.1993232  1.035853 -3.146711 1.089575
#> 
#> , , zg
#> 
#>     Estimate  Est.Error        Q2.5     Q97.5
#> 1 0.15976724 0.08921813 -0.01071167 0.3196564
#> 2 0.24963216 0.08264031  0.08816210 0.4097783
#> 3 0.08353011 0.07200988 -0.04812716 0.2280746
#> 4 0.36529945 0.04519791  0.27709860 0.4533866
#> 
#> , , zi_Intercept
#> 
#>    Estimate Est.Error      Q2.5      Q97.5
#> 1 -7.296146  6.817151 -22.84229  0.5702874
#> 2 -9.930848  8.403730 -29.46093 -0.7190626
#> 3 -9.532446  8.373872 -28.81624  0.2686403
#> 4 -9.502189  8.415209 -30.05622  0.1605108
#> 
#> , , zi_zg
#> 
#>    Estimate Est.Error      Q2.5      Q97.5
#> 1 -8.514177  6.120933 -25.27990 -1.9568175
#> 2 -6.563639  4.558390 -16.05121 -1.8399470
#> 3 -8.958289  5.989657 -21.74333 -2.6511637
#> 4 -4.571465  3.546923 -13.44316 -0.8261003
#> 
#> , , zi_nofish1
#> 
#>      Estimate Est.Error      Q2.5    Q97.5
#> 1 -0.69059771  2.558553 -6.878615 3.391735
#> 2 -0.66113253  1.886437 -5.312619 2.324461
#> 3 -0.55772179  2.194537 -6.212487 3.026720
#> 4 -0.08506966  1.992193 -4.468332 4.176201
#> 
#> , , child
#> 
#>     Estimate Est.Error      Q2.5      Q97.5
#> 1 -0.4544075 0.1220026 -0.698228 -0.2138667
#> 2 -0.4544075 0.1220026 -0.698228 -0.2138667
#> 3 -0.4544075 0.1220026 -0.698228 -0.2138667
#> 4 -0.4544075 0.1220026 -0.698228 -0.2138667
#> 
#> , , camper1
#> 
#>    Estimate Est.Error        Q2.5     Q97.5
#> 1 0.1078043 0.1110276 -0.09501827 0.3259446
#> 2 0.1078043 0.1110276 -0.09501827 0.3259446
#> 3 0.1078043 0.1110276 -0.09501827 0.3259446
#> 4 0.1078043 0.1110276 -0.09501827 0.3259446
#> 
#> , , zi_child
#> 
#>     Estimate Est.Error      Q2.5    Q97.5
#> 1 -0.7358728  3.085961 -7.267424 4.137575
#> 2 -0.7358728  3.085961 -7.267424 4.137575
#> 3 -0.7358728  3.085961 -7.267424 4.137575
#> 4 -0.7358728  3.085961 -7.267424 4.137575
#> 
#> , , zi_livebait1
#> 
#>   Estimate Est.Error    Q2.5    Q97.5
#> 1 2.377975  5.663644 -5.9263 16.07785
#> 2 2.377975  5.663644 -5.9263 16.07785
#> 3 2.377975  5.663644 -5.9263 16.07785
#> 4 2.377975  5.663644 -5.9263 16.07785
#> 
#> 
#> $persons
#> , , Intercept
#> 
#>     Estimate Est.Error      Q2.5      Q97.5
#> 1 -1.4726950 0.5559270 -2.558071 -0.4086148
#> 2 -1.1021684 0.5289535 -2.150847 -0.1149768
#> 3 -1.2686065 0.5395648 -2.328640 -0.2619923
#> 4 -0.1746891 0.5134628 -1.198719  0.7418030
#> 
#> , , xb
#> 
#>    Estimate  Est.Error      Q2.5     Q97.5
#> 1 1.5626166 0.20201624 1.1732615 1.9904038
#> 2 1.2953383 0.09746668 1.0999387 1.4829048
#> 3 1.3196510 0.09650849 1.1370000 1.5147897
#> 4 0.8891383 0.04823144 0.7929538 0.9809271
#> 
#> , , zi_Intercept
#> 
#>    Estimate Est.Error      Q2.5    Q97.5
#> 1 -3.854074  4.813789 -14.58513 4.052241
#> 2 -3.854074  4.813789 -14.58513 4.052241
#> 3 -3.854074  4.813789 -14.58513 4.052241
#> 4 -3.854074  4.813789 -14.58513 4.052241
#> 
#> , , child
#> 
#>     Estimate Est.Error      Q2.5      Q97.5
#> 1 -0.4544075 0.1220026 -0.698228 -0.2138667
#> 2 -0.4544075 0.1220026 -0.698228 -0.2138667
#> 3 -0.4544075 0.1220026 -0.698228 -0.2138667
#> 4 -0.4544075 0.1220026 -0.698228 -0.2138667
#> 
#> , , camper1
#> 
#>    Estimate Est.Error        Q2.5     Q97.5
#> 1 0.1078043 0.1110276 -0.09501827 0.3259446
#> 2 0.1078043 0.1110276 -0.09501827 0.3259446
#> 3 0.1078043 0.1110276 -0.09501827 0.3259446
#> 4 0.1078043 0.1110276 -0.09501827 0.3259446
#> 
#> , , zi_child
#> 
#>     Estimate Est.Error      Q2.5    Q97.5
#> 1 -0.7358728  3.085961 -7.267424 4.137575
#> 2 -0.7358728  3.085961 -7.267424 4.137575
#> 3 -0.7358728  3.085961 -7.267424 4.137575
#> 4 -0.7358728  3.085961 -7.267424 4.137575
#> 
#> , , zi_livebait1
#> 
#>   Estimate Est.Error    Q2.5    Q97.5
#> 1 2.377975  5.663644 -5.9263 16.07785
#> 2 2.377975  5.663644 -5.9263 16.07785
#> 3 2.377975  5.663644 -5.9263 16.07785
#> 4 2.377975  5.663644 -5.9263 16.07785

### formatted output from parameters

model_parameters(x, effects = "total")
#> Group   | Level | Parameter | Median |          95% CI |     pd |     Component | Effects
#> -----------------------------------------------------------------------------------------
#> persons |     1 | Intercept |  -1.45 | [ -2.56, -0.41] | 99.62% |   conditional |   total
#> persons |     2 | Intercept |  -1.07 | [ -2.15, -0.11] | 98.35% |   conditional |   total
#> persons |     3 | Intercept |  -1.27 | [ -2.33, -0.26] | 98.72% |   conditional |   total
#> persons |     4 | Intercept |  -0.13 | [ -1.20,  0.74] | 62.60% |   conditional |   total
#> persons |     1 |        xb |   1.57 | [  1.17,  1.99] |   100% |   conditional |   total
#> persons |     2 |        xb |   1.30 | [  1.10,  1.48] |   100% |   conditional |   total
#> persons |     3 |        xb |   1.32 | [  1.14,  1.51] |   100% |   conditional |   total
#> persons |     4 |        xb |   0.89 | [  0.79,  0.98] |   100% |   conditional |   total
#> persons |     1 | Intercept |  -3.56 | [-14.59,  4.05] | 83.65% | zero_inflated |   total
#> persons |     2 | Intercept |  -3.56 | [-14.59,  4.05] | 83.65% | zero_inflated |   total
#> persons |     3 | Intercept |  -3.56 | [-14.59,  4.05] | 83.65% | zero_inflated |   total
#> persons |     4 | Intercept |  -3.56 | [-14.59,  4.05] | 83.65% | zero_inflated |   total
#> ID      |     1 | Intercept |  -0.89 | [ -2.70,  1.48] | 79.97% |   conditional |   total
#> ID      |     2 | Intercept |  -0.81 | [ -2.59,  1.61] | 78.05% |   conditional |   total
#> ID      |     3 | Intercept |  -0.53 | [ -2.28,  1.86] | 70.40% |   conditional |   total
#> ID      |     4 | Intercept |  -1.32 | [ -3.15,  1.09] | 87.22% |   conditional |   total
#> ID      |     1 |        zg |   0.16 | [ -0.01,  0.32] | 96.67% |   conditional |   total
#> ID      |     2 |        zg |   0.25 | [  0.09,  0.41] | 99.85% |   conditional |   total
#> ID      |     3 |        zg |   0.08 | [ -0.05,  0.23] | 86.28% |   conditional |   total
#> ID      |     4 |        zg |   0.37 | [  0.28,  0.45] |   100% |   conditional |   total
#> ID      |     1 | Intercept |  -5.98 | [-22.84,  0.57] | 95.97% | zero_inflated |   total
#> ID      |     2 | Intercept |  -7.43 | [-29.46, -0.72] | 98.72% | zero_inflated |   total
#> ID      |     3 | Intercept |  -7.23 | [-28.82,  0.27] | 96.85% | zero_inflated |   total
#> ID      |     4 | Intercept |  -7.13 | [-30.06,  0.16] | 97.20% | zero_inflated |   total
#> ID      |     1 |        zg |  -7.23 | [-25.28, -1.96] |   100% | zero_inflated |   total
#> ID      |     3 |        zg |  -7.13 | [-21.74, -2.65] |   100% | zero_inflated |   total
#> ID      |     2 |   nofish1 |  -0.23 | [ -5.31,  2.32] | 61.22% | zero_inflated |   total
#> ID      |     4 |   nofish1 |  -0.08 | [ -4.47,  4.18] | 54.62% | zero_inflated |   total

@DominiqueMakowski
Copy link
Member

In this example the zi term is on the conditional component, but what if you have it on say another auxiliarity component, in a model like:

bf(Resp ~ Formula1,
   sigma ~ Formula_with_zi)

@strengejacke
Copy link
Member Author

From ?lme4:::coef.merMod:

image

I still would prefer "total" for now.

@strengejacke
Copy link
Member Author

In this example the zi term is on the conditional component, but what if you have it on say another auxiliarity component, in a model like:

bf(Resp ~ Formula1,
   sigma ~ Formula_with_zi)

How can we make a "nested" formula? (i.e. sigma ~ zi_formula)?

@bwiernik
Copy link
Contributor

bwiernik commented Mar 3, 2025

In this example the zi term is on the conditional component, but what if you have it on say another auxiliarity component, in a model like:

bf(Resp ~ Formula1,

sigma ~ Formula_with_zi)

How can we make a "nested" formula? (i.e. sigma ~ zi_formula)?

You need to use brms' nonlinear formula syntax. See the last example here https://rdrr.io/cran/brms/man/brmsformula.html

@DominiqueMakowski
Copy link
Member

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

m4 <- brms::brm(brms::bf(Sepal.Width ~ Petal.Width + (Petal.Width | Group),
     sigma ~ Petal.Width + (Petal.Width | Group)),
    data = data, refresh=0)

@strengejacke
Copy link
Member Author

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

m4 <- brms::brm(brms::bf(Sepal.Width ~ Petal.Width + (Petal.Width | Group),
     sigma ~ Petal.Width + (Petal.Width | Group)),
    data = data, refresh=0)

That's no zero-inflation formula.

@DominiqueMakowski
Copy link
Member

I wonder if the zi parts should not be a different Parameter rather than Component

@strengejacke
Copy link
Member Author

Assuming that auxiliary components have no zero-inflation components, that stuff should work now, too. I bet there are more cases that won't work and which are hard to anticipate, so we can "fix on demand".

@strengejacke
Copy link
Member Author

I wonder if the zi parts should not be a different Parameter rather than Component

That would be a different behaviour than what we usually have across the easystats verse.

@strengejacke
Copy link
Member Author

strengejacke commented Mar 3, 2025

x <- insight::download_model("brms_zi_4")
parameters::model_parameters(x)
#> # Fixed Effects
#> 
#> Parameter   | Median |         95% CI |     pd |  Rhat |     ESS
#> ----------------------------------------------------------------
#> (Intercept) |  -0.67 | [-2.54,  1.78] | 72.47% | 1.020 |  214.00
#> child       |  -0.45 | [-0.70, -0.21] |   100% | 1.001 | 2571.00
#> camper1     |   0.11 | [-0.10,  0.33] | 81.30% | 1.025 |  174.00
#> 
#> # Zero-Inflation
#> 
#> Parameter   | Median |          95% CI |     pd |  Rhat |     ESS
#> -----------------------------------------------------------------
#> (Intercept) |  -3.56 | [-14.59,  4.05] | 83.65% | 1.003 | 1061.00
#> child       |  -0.63 | [ -7.27,  4.14] | 59.80% | 1.035 |  103.00
#> livebait1   |   1.52 | [ -5.93, 16.08] | 67.30% | 1.006 |  648.00
#> 
#> Uncertainty intervals (equal-tailed) computed using a MCMC distribution
#>   approximation.
#> 
#> The model has a log- or logit-link. Consider using `exponentiate =
#>   TRUE` to interpret coefficients as ratios.

Created on 2025-03-03 with reprex v2.1.1

@strengejacke strengejacke merged commit 90b6d4c into main Mar 3, 2025
16 of 23 checks passed
@strengejacke strengejacke deleted the strengejacke/issue1067 branch March 3, 2025 17:59
## Changes

* The `effects` argument in `model_parameters()` for classes `merMod`, `glmmTMB`,
`brmsfit` and `stanreg` gets an additional `"total"` option, to return the
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think now that "random_unbiased" would be the best name...

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, it's not finalized yet, we can change the name, if we find a name we're all happy with.

@ your component comment: We always had it composed like this, where the parameter names are "clean" and the related component is indicated in the Component column.

Should we now revise estimate_grouplevel()? It seems it could just be a wrapper around model_parameters(..., effects = "random_<whatever>"), or is estimate_grouplevel() doing more? The implementation here supports lme4/glmmTMB/rstanarm/brms.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But random_unbiased suggests the other effects are biased. But they're just deviations from the overall mean. Just curious, why are you so hesitant with "total" or "overall", since it's really just summing up fixed+random effects...

And: we should have the same options for the type argument in estimate_grouplevel(), or not?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we now revise estimate_grouplevel()?

Yes!

But random_unbiased

The main thing is to avoid effects="total" because it's too vague. Since it concerns random effects, I'd at least name it "random_total" so that it's clear what's it about

Now about unbiased, well I'm not an expert of the underlying mathematics, but:

  • "It simply computes the "total" of random + fixed": in general yes, but I'm not 100% that is always exactly the case. As @bwiernik mentioned, sometimes there is no fixed effect and I'm not sure what happens then. Also, In the coef.brmsfit code, there seem to be some edgecases with some models etc. etc. So I'm not sure it is always just as simple as the sum.
  • "unbiased": initially I wasn't a fan either. But 1) it has the benefit of directly linking to the now well-established literature about BLUPs (best linear unbiased). And one could understand "biased" as "relative to the fixed effect", whereas "unbiased" means "absolute". Technically we could even consider "random_relative" and "random_absolute" 🤷

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do BLUPs refer to the sum of fixed and random effects? I thought it was just a description of the estimation method of random effects in general?

When there are not fixef, it's still the sum, with fixef's being zero:

m1 <- lme4::lmer(Reaction ~ Days + (Days | Subject), lme4::sleepstudy)
m2 <- lme4::lmer(Reaction ~ 1 + (Days | Subject), lme4::sleepstudy)

coef(m1)
#> $Subject
#>     (Intercept)       Days
#> 308    253.6637 19.6662617
#> 309    211.0064  1.8476053
#> 310    212.4447  5.0184295
#> 330    275.0957  5.6529356
#> 331    273.6654  7.3973743
#> 332    260.4447 10.1951090
#> 333    268.2456 10.2436499
#> 334    244.1725 11.5418676
#> 335    251.0714 -0.2848792
#> 337    286.2956 19.0955511
#> 349    226.1949 11.6407181
#> 350    238.3351 17.0815038
#> 351    255.9830  7.4520239
#> 352    272.2688 14.0032871
#> 369    254.6806 11.3395008
#> 370    225.7921 15.2897709
#> 371    252.2122  9.4791297
#> 372    263.7197 11.7513080
#> 
#> attr(,"class")
#> [1] "coef.mer"
lme4::ranef(m1)
#> $Subject
#>     (Intercept)        Days
#> 308   2.2585509   9.1989758
#> 309 -40.3987381  -8.6196806
#> 310 -38.9604090  -5.4488565
#> 330  23.6906196  -4.8143503
#> 331  22.2603126  -3.0699116
#> 332   9.0395679  -0.2721770
#> 333  16.8405086  -0.2236361
#> 334  -7.2326151   1.0745816
#> 335  -0.3336684 -10.7521652
#> 337  34.8904868   8.6282652
#> 349 -25.2102286   1.1734322
#> 350 -13.0700342   6.6142178
#> 351   4.5778642  -3.0152621
#> 352  20.8636782   3.5360011
#> 369   3.2754656   0.8722149
#> 370 -25.6129993   4.8224850
#> 371   0.8070461  -0.9881562
#> 372  12.3145921   1.2840221
#> 
#> with conditional variances for "Subject"
lme4::fixef(m1)
#> (Intercept)        Days 
#>   251.40510    10.46729

coef(m2)
#> $Subject
#>           Days (Intercept)
#> 308 20.6010269    249.4572
#> 309  0.2135211    218.3598
#> 310  3.8918664    217.5142
#> 330  4.1171150    282.0069
#> 331  6.1589652    279.2383
#> 332  9.5231905    263.4683
#> 333  9.5157838    271.5210
#> 334 11.2238002    245.6038
#> 335 -2.5963842    261.4733
#> 337 19.6695601    283.7125
#> 349 11.4860944    226.8907
#> 350 17.7185058    235.4685
#> 351  6.3674017    260.8638
#> 352 13.8582027    272.9217
#> 369 10.9022191    256.6483
#> 370 15.7360708    223.7837
#> 371  8.7573841    255.4600
#> 372 11.3074308    265.7171
#> 
#> attr(,"class")
#> [1] "coef.mer"
lme4::ranef(m2)
#> $Subject
#>     (Intercept)       Days
#> 308   -8.304993 20.6010269
#> 309  -39.402398  0.2135211
#> 310  -40.247926  3.8918664
#> 330   24.244775  4.1171150
#> 331   21.476108  6.1589652
#> 332    5.706141  9.5231905
#> 333   13.758848  9.5157838
#> 334  -12.158381 11.2238002
#> 335    3.711084 -2.5963842
#> 337   25.950356 19.6695601
#> 349  -30.871497 11.4860944
#> 350  -22.293637 17.7185058
#> 351    3.101614  6.3674017
#> 352   15.159485 13.8582027
#> 369   -1.113833 10.9022191
#> 370  -33.978438 15.7360708
#> 371   -2.302158  8.7573841
#> 372    7.954974 11.3074308
#> 
#> with conditional variances for "Subject"
lme4::fixef(m2)
#> (Intercept) 
#>    257.7622

Created on 2025-03-03 with reprex v2.1.1

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So why is it called BLUPs then? Why such a complicated name for "simply" taking the random parameters

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

They are the best linear unbiased predictions of the unknown parameters (random or total slope coefficients), the same way that the OLS coefficients are the best linear unbiased estimates (BLUEs) of those unknown coefficients.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That all makes sense. Incredible that I never found any clear clarification in a lot of reading. AI companies should trains their LLMs on the easystats' github issues & PRs it's a goldmine of knowledge lol

Then:

  • effects="random"
  • effects="random_total"

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it's a goldmine of knowledge lol

😂

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

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

3 participants