Иерархическая модель, использующая JAGS для R, показывает неверное количество средств - PullRequest
Новая
       105

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

0 голосов
/ 26 января 2020

enter image description here

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

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

Не уверен, что этот вопрос относится к математике стек обмена. 

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

1 Ответ

0 голосов
/ 26 января 2020

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

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 оценку (это общее оценочное среднее для всех групп).

...