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

Учитывая примеры данных sampleDT и модели brmsbrm.fit и brm.fit.distr ниже, я хотел бы:

estimate, extract and add to the data frame the values of the standard deviations for each observation from the distributional model 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

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

Я считаю, что в brm.fit.distr больше нет единственного sigma; сравните его вывод с brm.fit. Это правильно? ?bf также говорит По умолчанию стандартные ошибки заменяют параметр сигма.

Julius Vainora 10.02.2019 12:47

Вот и все коэффициенты с sigma от summary(brm.fit.distr)$fixed; так что ни одного числа за наблюдение. Возможно, вы хотите predict(brm.fit.distr)[, 2]?

Julius Vainora 10.02.2019 12:50
Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
1
2
122
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

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

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

@Krantz, я намеренно был несколько расплывчатым в своем ответе, так как не уверен, какие вещи вы хотите и даже что правильно, но я надеюсь, что это поможет.

Julius Vainora 10.02.2019 13:21

Но я не уверен в масштабе, потому что эти числа из predict(brm.fit.distr)[, 2], как вы говорите, очень большие, и когда я делаю log(predict(brm.fit.distr)[, 2]), я получаю более приемлемые значения по сравнению с числами из summary(brm.fit)$spec_pars[1]. Есть предположения?

Krantz 10.02.2019 13:36

@Krantz, я считаю, что единицы действительно должны быть правильными, и я объясняю в своем ответе, почему иногда число велико: A и B равны 212 и 170 в 5-м наблюдении, например. Также эти числа являются случайными и модель не сошлась; Rhat для sigma_A и sigma_B особенно велики.

Julius Vainora 10.02.2019 13:38

Большое спасибо, @JuliusVainora, за пояснение.

Krantz 10.02.2019 13:40

Другие вопросы по теме