I would like to run the following hierarchical models in brms and then conduct kfold on each as follows:

```
f1 <- bf('y1~1+var1+var2+var3+(1+var1+var2+var3|GroupID)')
f2 <- bf('y2~1+var1+var2+var3+(1+var1+var2+var3|GroupID)')
f3 <- bf('y1~1+var1+var2+(1+var1+var2+|GroupID)')
f4 <- bf('y2~1+var1+var2+(1+var1+var2+|GroupID)')
fit1 <- brm(formula = f1,data=df,iter=6000,chains=4,cores=4,control=list(adapt_delta=0.9),family=gaussian(link="identity"),save_all_pars = TRUE))
fit2 <- brm(formula = f2,data=df,iter=6000,chains=4,cores=4,control=list(adapt_delta=0.9),family=gaussian(link="identity"),save_all_pars = TRUE))
fit3 <- brm(formula = f3,data=df,iter=6000,chains=4,cores=4,control=list(adapt_delta=0.9),family=gaussian(link="identity"),save_all_pars = TRUE))
fit4 <- brm(formula = f4,data=df,iter=6000,chains=4,cores=4,control=list(adapt_delta=0.9),family=gaussian(link="identity"),save_all_pars = TRUE))
kfold1 <- brms::kfold(fit1,K=10,folds="stratified",group="GroupID")
kfold2 <- brms::kfold(fit2,K=10,folds="stratified",group="GroupID")
kfold3 <- brms::kfold(fit3,K=10,folds="stratified",group="GroupID")
kfold4 <- brms::kfold(fit4,K=10,folds="stratified",group="GroupID")
```

Then I would like to combine the resulting `elpd_kfold`

estimates from `kfold1`

and `kfold2`

to make a composite value, and compare this composite value with the combined `elpd_kfold`

estimates from `kfold3`

and `kfold4`

. This way I could assess whether, on the whole, the first set of models will generalize better than the second set of models.

I am choosing this method over entering two complex multivariate models into `kfold`

with `compare=TRUE`

for a few reasons: my dataset is large (N=4830) and running `kfold`

on large complex multivariate models takes too much time and often exceeds my application memory; and I often get diagnostic errors from running the models together as multivariate models. It seems like running the models separately and then combining the `elpd_kfold`

estimates may be the best solution as I am not interested in modelling residual correlations.

I am wondering whether 1) this method makes sense, and 2) how to compute a composite standard error for each set of models? Thanks, I appreciate any advice.