How to iteratively fit a brms regression model and extract means and sigma to the dataframe

Given the sample data sampleDT below, I would appreciate any help to iteratively fit the brms model below n times, and each time extract the means and sigma from the brmsfit object brm.fit.n and add them to the data frame sampleDT.

If n=10, then there should be 10 columns of means and 10 columns of sigma added to the data frame.

My attempt below does not work as intended. It allows me to run the brms model n times and generate the means and sigma n times, but does not add them to the data frame - one column for each means and one column for each sigma from each run - as intended.

I suspect there might be a natural approach to accomplish this with brms that I am not aware of. Therefore, any help would be much appreciated.

#sample data

sampleDT<-structure(list(id = 1:10, N = c(10L, 10L, 10L, 10L, 10L, 10L, 
    10L, 10L, 10L, 10L), A = c(62L, 96L, 17L, 41L, 212L, 143L, 143L, 
    143L, 73L, 73L), B = c(3L, 1L, 0L, 2L, 170L, 21L, 0L, 33L, 62L, 
    17L), C = c(0.05, 0.01, 0, 0.05, 0.8, 0.15, 0, 0.23, 0.85, 0.23
    ), employer = c(1L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L), F = c(0L, 
    0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L), G = c(1.94, 1.19, 1.16, 
    1.16, 1.13, 1.13, 1.13, 1.13, 1.12, 1.12), H = c(0.14, 0.24, 
    0.28, 0.28, 0.21, 0.12, 0.17, 0.07, 0.14, 0.12), dollar.wage_1 = c(1.94, 
    1.19, 3.16, 3.16, 1.13, 1.13, 2.13, 1.13, 1.12, 1.12), dollar.wage_2 = c(1.93, 
    1.18, 3.15, 3.15, 1.12, 1.12, 2.12, 1.12, 1.11, 1.11), dollar.wage_3 = c(1.95, 
    1.19, 3.16, 3.16, 1.14, 1.13, 2.13, 1.13, 1.13, 1.13), dollar.wage_4 = c(1.94, 
    1.18, 3.16, 3.16, 1.13, 1.13, 2.13, 1.13, 1.12, 1.12), dollar.wage_5 = c(1.94, 
    1.19, 3.16, 3.16, 1.14, 1.13, 2.13, 1.13, 1.12, 1.12), dollar.wage_6 = c(1.94, 
    1.18, 3.16, 3.16, 1.13, 1.13, 2.13, 1.13, 1.12, 1.12), dollar.wage_7 = c(1.94, 
    1.19, 3.16, 3.16, 1.14, 1.13, 2.13, 1.13, 1.12, 1.12), dollar.wage_8 = c(1.94, 
    1.19, 3.16, 3.16, 1.13, 1.13, 2.13, 1.13, 1.12, 1.12), dollar.wage_9 = c(1.94, 
    1.19, 3.16, 3.16, 1.13, 1.13, 2.13, 1.13, 1.12, 1.12), dollar.wage_10 = c(1.94, 
    1.19, 3.16, 3.16, 1.13, 1.13, 2.13, 1.13, 1.12, 1.12)), row.names = c(NA, 
    -10L), class = "data.frame")

#my attempt

map_dfc(1:10, function(i) {
        brm.fit.n <-brm(dollar.wage_1 ~ A + B + C + employer + F + G + H,
                data=sampleDT, iter = 200, family = gaussian())
        sampleDT$mean.n<-fitted(brm.fit.n)[, 1]
        sampleDT$sd.n<-summary(brm.fit.n)$spec_pars[1]
        return(sampleDT)
    })

This question has also been posted here. Thanks in advance for any help.

Did you figure out how to do this? I’m trying to do something similar at the moment…
The reason I want to do it is to get a separate estimate of sigma (for the normal distribution) for each ‘batch’ and I can’t see how to do that in a hierarchical model. Was this also your reason for doing it?

Also, there is this blog post here: https://cran.r-project.org/web/packages/brms/vignettes/brms_distreg.html that shows how to model the group variance separately

Hey all,

I’ve been working through the part in Kruschke’s text where he covers Bayesian power analysis using simulations (section 13.2.4). The form of the model is a little different, but the workflow might be of help.