Извлеките и добавьте во фрейм данных значения сигмы из линейной модели распределения Стэн - PullRequest
Новая
       96

Извлеките и добавьте во фрейм данных значения сигмы из линейной модели распределения Стэн

0 голосов
/ 10 февраля 2019

Учитывая данные выборки 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

Заранее благодарен за любую помощь.

1 Ответ

0 голосов
/ 10 февраля 2019

Как и ожидалось в байесовских моделях, существуют разные способы оценки степени неопределенности.Итак, во-первых, у нас больше нет единственного параметра 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

, которые также являются случайными.

...