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

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

оценить, извлечь и добавить к кадру данных значения функции плотности для условного нормального распределения, оцененного на наблюдаемом уровне переменной dollar.wage_1.

Я могу сделать это, используя частотную линейную регрессию lm.fit и dnorm, но моя попытка сделать то же самое, используя байесовскую brm.fit модель, терпит неудачу. Поэтому любая помощь будет принята с благодарностью.

## пример данных

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")

## частотная модель: это работает

lm.fit <-lm(dollar.wage_1 ~ A + B + C + employer + F + G + H,
            data=sampleDT)

sampleDT$dens1 <-dnorm(sampleDT$dollar.wage_1,mean=lm.fit$fitted,
sd=summary(lm.fit)$sigma)

## байесовская модель: это моя попытка - не работает

//this works
brm.fit <-brm(dollar.wage_1 ~ A + B + C + employer + F + G + H,
            data=sampleDT, iter = 4000, family = gaussian())

//this does not work
 sampleDT$dens1_bayes <-dnorm(sampleDT$dollar.wage_1, mean = fitted(brm.fit), sd=summary(brm.fit)$sigma)

Ошибка в dnorm (sampleDT $ dollar.wage_1, mean = brm.fit $ fit, sd = summary (brm.fit) $ sigma): нечисловой аргумент математической функции

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


r rstan bayesian probability-density stan
person Krantz    schedule 08.02.2019    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Теперь у нас есть fitted(brm.fit) - это матрица, поэтому мы хотим использовать только ее первый столбец - столбец оценок. Кроме того, поскольку нет причин для одинаковой структуры объекта, summary(brm.fit)$sigma ничего не дает. Вместо этого вы хотите summary(brm.fit)$spec_pars[1]. Следовательно, вы можете использовать

sampleDT$dens1_bayes <- dnorm(sampleDT$dollar.wage_1,
                              mean = fitted(brm.fit)[, 1],
                              sd = summary(brm.fit)$spec_pars[1])
person Julius Vainora    schedule 08.02.2019
comment
Отлично, @JuliusVainora. спасибо за ответ. Очень полезно. Но я немного обеспокоен. Почему эти sampleDT $ dens1_bayes сильно отличаются от sampleDT $ dens1? Я заметил, что > sampleDT$dens1_freq [1] 0.5313967 0.4377899 0.5309715 0.4308041 0.5297744 0.5247409 0.5275020 0.4069652 0.5295822 0.3930264, тогда как > sampleDT$dens1_bayes [1] 0.1644518 0.1613566 0.1644267 0.1621689 0.1644273 0.1641519 0.1642465 0.1591944 0.1642170 0.1601089. Разве эти двое не должны быть как-то примерно равны, по крайней мере, не слишком сильно различаться? - person Krantz; 08.02.2019
comment
@Krantz, разница исходит от sd, который в байесовской модели вдвое больше. Я не совсем уверен, как читать (части 0 и 10) априорное значение sigma в brm.fit$prior, но 3 степени свободы для t-распределения в сочетании с очень маленькой выборкой могут означать, что апостериорная неопределенность останется высокой. Единственное другое объяснение было бы то, что этот sigma - какой-то совершенно другой параметр, но я сомневаюсь в этом. - person Julius Vainora; 08.02.2019
comment
да. Как вы говорите, одна возможность состоит в том, что это из-за very small sample, а другая - из-за the prior of sigma in brm.fit$prior (...) 3 degrees of freedom for a t distribution. sd очень разные: DT$sd_freq [1] 0.7506886 тогда как sampleDT$sd_bayes [1] 2.425812. Данные совпадают, поэтому результаты байесовского и частотного анализа не должны сильно отличаться, как при использовании default configurations пакетов. Есть предположения? - person Krantz; 08.02.2019
comment
@Krantz, я бы предложил сгенерировать некоторые данные и попытаться оценить модель, используя небольшую подвыборку и полную (относительно большую) выборку. И сравнить полученное среднее sd с истинным. - person Julius Vainora; 08.02.2019
comment
Спасибо, @JuliusVainora. Я сделаю это. - person Krantz; 08.02.2019
comment
Привет, @JuliusVainora. Связанный вопрос размещен на странице stackoverflow.com/questions/54615821/. Заранее благодарю за любую помощь. - person Krantz; 10.02.2019