r - 从 stan 分布线性模型中提取 sigma 的值并将其添加到数据框中
问题描述
鉴于示例数据sampleDT
和brms
模型brm.fit
及brm.fit.distr
以下,我想:
估计、提取分布模型中每个观测值的标准差值并将其添加到数据框中
brm.fit.distr
。
我可以使用 来做到这一点brm.fit
,但是当我使用brm.fit.distr
.
样本数据
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")
我的模型
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())
我的方法brm.fit
和尝试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
提前感谢您的帮助。
解决方案
正如贝叶斯模型中所预期的那样,有不同的方法可以查看不确定性的程度。所以,首先,我们不再有一个参数sigma
;相反,有几个标准偏差参数
summary(brm.fit.distr)$fixed
特别是,
exp(summary(brm.fit.distr)$fixed[, 1])[grep("sigma", rownames(summary(brm.fit.distr)$fixed))]
# sigma_Intercept sigma_A sigma_B sigma_C sigma_employer
# 1.17043390 0.99913160 1.01382623 0.28655150 1.06713923
# sigma_F sigma_G sigma_H
# 0.50428952 0.87669186 0.01203015
我exp
用来使数字为正数的地方。
现在作为不确定性的综合衡量标准,我们可以看看
predict(brm.fit.distr)[, 2]
请注意,这些是随机的(!)在某些情况下,这些数字非常大
predict(brm.fit.distr)[, 2]
# [1] 34.620936 4.456770 2.837869 1.727396 107.116980 2.238100 2.350523 3.037880
# [9] 6.266055 2.517457
但我们有,例如,
sampleDT[5, 1:5]
# id N A B C
# 5 5 10 212 170 0.8
A
因此和的值B
非常大。同样,您可以查看
predict(brm.fit)[, 2]
# [1] 5.203937 4.846928 4.960600 4.827138 4.937323 4.625976 5.122794 4.767257 4.862458 4.219394
这也是随机的。
推荐阅读
- javascript - 如何以角度验证“_”或“.”(任何一个)
- android - 如何使用 adb 命令或类似命令在 3rd 方应用程序上执行命令
- linux - 为什么 kill -2 不杀死进程?
- r - ggplot2 的匹配图例
- html - 如何在手机上以正确的分辨率打开我的 HTML 项目?
- java - 注销后如何让webview记住用户名和密码等用户凭据以在android中自动填充?
- python - lmdb.Error:磁盘空间不足
- sql-server - 在 TSQL 中将 Pivot 函数与 Sum 函数结合使用
- java - 第二种方法的 Java Scanner
- tabulator - 如何在初始设置之外设置 Tabulator groupValues