Have coef return all nested random effects, not just one level

[bakaburg1]

Dear Paul,

I have a nested model of this form y ~ (1 | A/B/C).

I'd like to get the baseline estimates for each level of the model, i.e.:

1. Intercept + A
2. Intercept + A + A:B
3. Intercept + A + A:B + A:B:C

I thought of using coef() since it already gives you the join estimate for the fixed and random effect, but then I notice that it won't manage nested effects: coef(model)$A = Intercept + A (OK!) coef(model)$A:B = Intercept + A:B (Not OK!) coef(model)$A:B:C = Intercept + A:B:C (Not OK!)

(I wouldn't even know how to interpret these estimates)

Is there a way to have the correct estimates without having to build them manually? (once you add slopes it gets quite messy)

The only easy alternative I could find is to use posterior_epred but of course, it has much more overhead (I loop over several models).

If at the moment there aren't any alternatives, this would be a very much appreciated new feature.

Thank you for your work as always.

[paul-buerkner]

Yes, this would be a sensible feature indeed. I will think of a way to make that possible in the future.

[bakaburg1]

I reimplemented coef() to add up nested RE. At the moment is tested and works for linear models like y ~ Pred + (Pred | A/B/../N) with Pred or intercept only, but maybe it could help you when adding the functionality. working with multidimensional arrays brought me quite the headaches...

nested_coefs <- function(model, pred = NULL, summary = TRUE, robust = FALSE,
                         probs = c(0.025, 0.975)) {
  fe <- fixef(model, summary = F)
  
  re <- ranef(model, summary = F)
  
  rnd_lvls <- names(re)
  
  draws <- lapply(1:length(rnd_lvls), \(i) {
    vars <- str_split(rnd_lvls[i], ":")[[1]]
    cases <- select(model$data, any_of(vars)) %>% distinct()
    
    re_nested <- lapply(1:i, \(j) {
      re[[j]][,do.call(paste, c(as.data.frame(cases[,1:j]), sep = "_")),] # sum along nested RE up to the current one
    }) %>% Reduce(f = '+')
    
    fe <- replicate(nrow(cases), fe, simplify = F) %>% # replicate FE a number of time equal to the n of random effects ...
      abind::abind(along = 1.5) %>%  # ... along a new dimension (1.5 is for that) and ...
      array(dim(re_nested)) # ... make it conformable with RE

    draws <- fe + re_nested
    
    dimnames(draws) <- dimnames(re[[i]])[1:length(dim(draws))]
    
    if (summary) {
      draws <- posterior_summary(draws, probs, robust)
    }
    
    draws
  }) %>% setNames(rnd_lvls)
}

[bakaburg1]

UPDATE: the way I propagate dimnames is wrong. Trying to fixing it now.

[bakaburg1]

Updated version. fixes label propagation.

nested_coefs <- function(model, pred = NULL, summary = TRUE, robust = FALSE,
                         probs = c(0.025, 0.975)) {
  fe <- fixef(model, summary = F)
  
  re <- ranef(model, summary = F)
  
  rnd_lvls <- names(re)
  
  draws <- lapply(1:length(rnd_lvls), \(i) {
    vars <- str_split(rnd_lvls[i], ":")[[1]]
    cases <- select(model$data, any_of(vars)) %>% distinct()
    
    re_nested <- lapply(i:1, \(j) { # From the innest to the outest level to keep the names during the reduction
      re[[j]][,do.call(paste, c(as.data.frame(cases[,1:j]), sep = "_")),] # sum along nested RE levels
    }) %>% Reduce(f = '+')
    
    fe_rep <- replicate(nrow(cases), fe, simplify = F) %>% # replicate FE a number of time equal to the n of random effects ...
      abind::abind(along = 1.5) %>%  # ... along a new dimension (1.5 is for that) and ...
      array(dim(re_nested)) # ... make it conformable with RE
    
    draws <- fe_rep + re_nested
    
    if (summary) {
      draws <- posterior_summary(draws, probs,robust)
    }
    
    draws
  }) %>% setNames(rnd_lvls)
}