I think your approach should work for hierarchical models as well. One thing i noticed is that you dont use your priors in `simulModel`

which is probably just a typo.

##
Long answer + example

**Data generation**

```
library(sjPlot)
library(ggplot2)
library(MASS)
library(brms)
library(plyr)
library(gridExtra)
# data generation function
gen_design <- function(n_participant, n_stimulus, n_levels){
participants <- rep(1:n_participant, each = n_stimulus)
#stims <- rep(1:n_stimulus, times = n_participant)
x <- c()
stims <- as.vector(sapply(seq(1, n_participant), function(x) sample(seq(1,n_stimulus), size = n_stimulus)))
for(i in 1:n_participant){
x <- append(x, seq(0,(n_levels-1)))
#x <- append(x, sapply(seq(1,n_stimulus), function(x) sample(seq(1,n_levels), size = n_levels)))
}
expdesign <- data.frame(participant = participants, stimulus = stims, x = x)
return(expdesign)
}
sp <- list(n_participant = 30, n_stimulus = 20, n_levels =2) # data with 30 participants each seeing 20 stimuli (which will not be modelled in this example) with 2 levels (10 obs per participant per cell)
exp1 <- gen_design(n_participant = sp$n_participant, n_stimulus = sp$n_stimulus, n_levels = sp$n_levels)
```

**priors**

```
all_priors <- c(set_prior("normal(10,1)", class = "Intercept"),
set_prior("normal(1,1)", class = "b"),
set_prior("normal(5,1)", class = "sd", coef = "Intercept", group = "participant"),
set_prior("normal(2,1)", class = "sd", coef = "x", group = "participant"),
set_prior("normal(10,1)", class = "sigma"))
```

**fit intial model**

```
y <- rnorm(1:nrow(exp1), 1, 0) # just some arbitrary values for this column
m1 <- brm(y ~ x + (1 + x | participant), data = exp1, sample_prior = "only", prior = all_priors)
```

**make and fit predictions**

```
new_y_sum <- predict(m1)[,1]
new_y<- predict(m1, newdata = exp1, summary = F, nsamples = 1)[1,]
exp1$new_y_sum <- new_y_sum
exp1$new_y <- new_y
m1_new_sum <- brm(new_y_sum ~ x + (1 + x | participant), data = exp1, chains = 4, cores = 4)
m1_new <- brm(new_y ~ x + (1 + x | participant), data = exp1, chains = 4, cores = 4)
```

**model output for predict(m1) - too small SEs**

```
fixef(m1_new_sum)
```

```
Estimate Est.Error Q2.5 Q97.5
Intercept 9.488095 0.01631371 9.4562210 9.520331
x 1.009386 0.01550281 0.9784959 1.039939
```

**model output for predict(m1, summary = F) intended SEs**

```
fixef(m1_new)
```

```
Estimate Est.Error Q2.5 Q97.5
Intercept 10.0320100 1.1200567 7.892897 12.145149
x 0.3468883 0.8816665 -1.388313 2.085174
```

##
More formal simualtion of SEs for predict(m1, summary = F)

```
set.seed(1)
fixef_vec_int <- list()
fixef_vec_x <- list()
fixef_vec_x_lb <- c()
for(i in 1:1000){
new_y <- predict(m1, newdata = exp1, summary = F, nsamples = 1)[1,]
exp1$new_y <- new_y
m1_tmp <- update(m1_new, newdata = exp1, chains = 1)
fixef_vec_int[[i]] <- c(fixef(m1_tmp)[1], fixef(m1_tmp)[3])
fixef_vec_x[[i]] <- c(fixef(m1_tmp)[2], fixef(m1_tmp)[4])
fixef_vec_x_lb[i] <- fixef(m1_tmp)[6]
}
fixef_x_means <- unlist(lapply(1:100, function(x) fixef_vec_x[[x]][1]))
fixef_x_sds <- unlist(lapply(1:100, function(x) fixef_vec_x[[x]][2]))
mean(fixef_x_means) # the mean of the coef-prior
mean(fixef_x_sds) # the mean of the sampled SE
```