Given the sample data sampleDT
and the brms
models brm.fit
and brm.fit.distr
below, I would like to:
estimate, extract and add to the data frame the values of the standard deviations for each observation from the distributional model
brm.fit.distr
.
I can do this using brm.fit
, but my approach fails when I use brm.fit.distr
.
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 models
library(brms)
brm.fit <-brm(dollar.wage_1 ~ A + B + C + employer + F + G + H,
data=sampleDT, iter = 4000, family = gaussian())
brm.fit.distr <-brm(bf(dollar.wage_1 ~ A + B + C + employer + F + G + H,
sigma ~ A + B + C + employer + F + G + H),
data=sampleDT, iter = 4000, family = gaussian())
My approach for brm.fit
and attempt for brm.fit.distr
sampleDT$sd_brm_fit<-summary(brm.fit)$spec_pars[1] //this works
sampleDT$sd_brm_fit_distr<-summary(brm.fit.distr)$spec_pars[1] //this does not work
Thanks in advance for any help.
As expected in Bayesian models, there are different ways to look at the extent of uncertainty. So, first, we no longer have a single parameter
sigma
; instead there are several standard deviation parameters inand, in particular,
where I use
exp
to make the number positive.Now as an aggregate measure of uncertainty we may look at
Note that those are random (!) In some cases those number are pretty large
but we have that, e.g.,
so that the values for
A
andB
are very large. Similarly you could look atwhich also are random.