CODEWITHSUNDEEP
rfunctiondata.tabletidyrrstan
Krantz
Krantz

Reputation: 1503

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.

#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.

Upvotes: 1

Views: 452

Answers (1)

Marius
Marius

Reputation: 60170

The two things you need to do to adapt your existing code into a function are:

  • Repeat the fit n times
  • Save the output in a nice structure

There are lots of ways to do both parts, one option is purrr::map_dfr which can do both, applying the fit multiple times and creating a dataframe.

Instead of a brm model, which takes some time to fit, I've fit a simple linear model to random data instead, you would just have to replace those parts with your fitting code and save the mean and sd instead of the intercept and slope:

library(purrr)

# 1:10 - will repeat 10 times
map_dfr(1:10, function(i) {
    random_data = data.frame(
        x = rnorm(30),
        y = rnorm(30)
    )
    fit = lm(y ~ x, data = random_data)
    intercept = coef(fit)[1]
    slope = coef(fit)[2]
    return(data.frame(intercept, slope))
}, .id = "sim_num")

Which gives a dataframe as output:

   sim_num    intercept       slope
1        1  0.274903632 -0.03529736
2        2 -0.005134599 -0.22063748
3        3 -0.134999713  0.12090366
4        4 -0.216886033  0.21265679
5        5  0.261365432  0.02434036
6        6  0.067069791  0.23180334
7        7 -0.235138217  0.28360061
8        8 -0.117489553  0.10781101
9        9 -0.150288480  0.03086797
10      10 -0.031814194 -0.04075479

Upvotes: 1

Related Questions