Иерархическая модель с использованием JAGS для R показывает неправильное количество средних

введите описание изображения здесь

У меня есть 2 столбца данных, y и grp, и я пытаюсь создать модель JAGS, показанную выше. grp - это группа, а у меня 5 групп. Следующий код взят из здесь. Я использую этот код, потому что описание под заголовком Model and Data похоже на эту иерархическую модель.

Но когда я просматриваю сводку, я получаю только один mu. Должно быть 5 mu's, по одному на каждую группу. Может кто поправить код? Вы также можете указать на аналогичный пример, доступный в другом месте, и я мог бы попытаться его изменить. Мне что-то не хватает в коде, и я считаю, что код может быть похож на вопрос, но когда я изменяю его таким образом, мне кажется, что я не получаю нужных средств, даже несмотря на то, что есть 5 средств.

Не уверен, относится ли этот вопрос к math stackexchange.

mod_string = " model {
   for (i in 1:length(y) {
    theta[i] ~ dnorm(mu[grp[i]], invTau2)
    y[i] ~ dnorm(theta[i], 1/sig)
  }
  mu ~ dnorm(0, 1e6)
  invTau2 ~ dgamma(1.0/2.0, 1.3/2.0)
  tau2 <- 1/invTau2

  invgamma2 ~ dgamma(1.0/2.0, 2.1/2.0)
  sig = 1/invgamma2
    } "


summary(mod_sim)

Iterations = 2001:52000
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 50000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

         Mean    SD  Naive SE Time-series SE
mu  5.639e-07 0.001 2.582e-06      2.582e-06
sig 1.570e+00 1.888 4.874e-03      7.068e-03

jags rjags hierarchical-bayesian
person Mohan Radhakrishnan    schedule 26.01.2020    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Похоже, вы используете вложенную индексацию, где grp - вектор длины y и описывает, к какой группе принадлежит каждый элемент y. Если каждой группе нужно собственное среднее значение, вам нужно создать 5 theta априорных значений. Это проще всего сделать в собственном цикле. Таким образом, строка вашей модели будет выглядеть примерно так:

mod_string <- " model {
   # loop for y
   for (i in 1:length(y)) {
     y[i] ~ dnorm(theta[grp[i]], sig)
   }
   # loop for theta
   for (g in 1:ntheta){
     theta[g] ~ dnorm(mu, tau2)
   }
   # other priors. JAGS uses precision so you need to use the reciprocal 
   # of 1e6 for the variance of mu.
   mu ~ dnorm(0, 0.000001)
   invTau2 ~ dgamma(1.0/2.0, 1.3/2.0)
   tau2 <- 1/invTau2

   invgamma2 ~ dgamma(1.0/2.0, 2.1/2.0)
   sig <- 1/invgamma2
} "

Обратите внимание, что вам также необходимо указать в модели скаляр ntheta, который будет числом уникальных групп (в данном случае ntheta = 5). Теперь, если вы отслеживаете theta, вы получите 5 оценок, тогда как mu имеет 1 оценку (это общее оценочное среднее всех групп).

person mfidino    schedule 26.01.2020
comment
y и grp - столбцы данных. Каждый y принадлежит к группе. Я думаю, это то, что вы имели в виду. Когда я выполняю код и отслеживаю тета, я получаю правильные значения. Я думаю. - person Mohan Radhakrishnan; 28.01.2020