Hi!

Honestly, I am not sure. I have played with simulation to explore the possibility that a misspecified hierarchical effect distribution biases the population-level effect estimates, but it is not really conclusive… It seemed that it depends greatly on the particular dataset.

I noted that student hierarchical effects were harder to sample from, with sometimes poorly recovered degrees of freedom.

What I would suggest doing:

- A simulation-based calibration of population-level effects, fitting models to data generated by the other hierarchical effect distribution. I’d like to test that in a case study, but it would take some time because I shall concentrate on something else.
- Compare the two models using loo, and see if one provides better out-of-sample predictions.

I wish you good success!

Lucas

EDIT : I join the code for the simulations

```
# Empty the environment
rm(list = ls())
# Load libraries
library(brms)
library(useful)
library(tidyverse)
# Set seed
set.seed(77)
# Number of categories
K = 40
# Population level parameter
alpha = 7
beta = 2
# Individual level parameters
alpha_K = rstudent_t(K,3, 0, 1)
beta_K = rstudent_t(K,3, 0, 1)
# Create predictor data frame
D = data.frame(time = rep(1:5, K),
cat = factor(rep(1:K, each = 5)))
# Create model matrices
u_int = build.x(formula = ~ cat, data = D, contrasts = F)[,-1]
u_slope <- build.x(formula = ~ cat:time, data = D, contrasts = F)[,-1]
# Compute linear predictor
mu <- alpha + beta * D$time + u_int %*% alpha_K + u_slope %*% beta_K
# Sample observed values
D$y <- rnorm(length(mu), mu, 1)
# Observe simulated data
D %>% ggplot(aes(x = time, y = y)) +
geom_point(aes(color = cat)) +
geom_smooth(se = F, method = "lm") +
geom_smooth(aes(color = cat), se = F, method = "lm")
# Fit a normal model to t data
form_norm <- bf(y ~ time + (1 + time||cat))
# Observe priors
get_prior(form_norm, data = D)
# Set prior
prior_norm <- c(prior(normal(0,5), class = b),
prior(normal(0,3), class = sd))
# Sample from the posterior
fit_norm <- brm(form_norm,
prior = prior_norm,
data = D,
chains = 3, cores = 3)
fit_norm
# Fit a t model to t data
form_t <- bf(y ~ time + (1 + time||gr(cat, dist = "student")))
# Observe priors
get_prior(form_t, data = D)
# Set prior
prior_t <- c(prior(normal(0,5), class = b),
prior(normal(0,3), class = sd),
prior(normal(0,10), class = df, group= cat))
# Sample from the posterior
fit_t <- brm(form_t,
prior = prior_t,
data = D,
chains = 3, cores = 3)
fit_t
```