Как указать вложенную модель - PullRequest
Новая
       50

Как указать вложенную модель

0 голосов
/ 11 июля 2020

Я использую runjags для моделирования некоторых иерархических данных. Я могу смоделировать один уровень иерархии, но не знаю, как расширить его на большее количество уровней. Я пытаюсь сделать это, используя метод 3 со страницы 24 книги «Байесовское иерархическое моделирование с использованием WinBUGS» Ники Бест и др., В котором используется вложенный l oop (в отличие от вложенной индексации).

Для одного уровня я могу смоделировать 

filestring <-
  "model{ 
    for(j in 1:Ninner){
      for(i in 1:N){
        y[j,i] ~ dnorm(beta + alpha[j], py)
      }
      alpha[j] ~ dnorm(0, taua)
    }

  beta ~ dnorm(0, 0.001)
  taua ~ dgamma(0.01, 0.01)
  py ~ dgamma(0.01, 0.1)
}"

INITS <- list(list(.RNG.seed=1, .RNG.name="base::Wichmann-Hill"), 
              list(.RNG.seed=2, .RNG.name="base::Wichmann-Hill"))
results <- run.jags(filestring, monitor=c("py", "beta", "alpha"), data=jags_data, sample=1e3, 
                    n.chains=2, inits=INITS, summarise=FALSE)

Затем я попытался добавить еще один уровень, используя 

for(k in 1:Nouter){
 for(j in 1:Ninner){
  for(i in 1:N){
    y[j,i] ~ dnorm(beta + alpha_in[j] + alpha_out[k], py)
} } }

, но получил сообщение об ошибке

Компиляция ошибка в строке 5. Попытка переопределить узел y [1,1]

Как мне расширить это, чтобы смоделировать другой уровень, первый из которых вложен? Спасибо.

Ниже приведены некоторые воспроизводимые данные, которые показывают структуру данных. Я использую sh для оценки случайных оценок как для outer_grp, так и для inner_grp. 

library(data.table)
library(runjags)

set.seed(12345)
dat <- data.table(outer_grp=rep(1:5, each=10), inner_grp=rep(1:10, each=5), y=rnorm(50), x=rnorm(50), time=1:5)

wdat = dcast(dat, inner_grp + outer_grp ~ time, value.var=c("y", "x"))
jags_data = c(setNames(
  lapply(split.default(wdat, substr(names(wdat), 1, 1)),as.matrix), 
  c("inner_grp", "outer_grp","x", "y")),
  N=5, Nouter=5, Ninner=10)

РЕДАКТИРОВАТЬ

Пожалуй хватит модели вроде ?? 

filestring <-
  "model{ 

     for(j in 1:Ninner){
      for(i in 1:N){
        y[j,i] ~ dnorm(beta + alpha_in[j] + alpha_out[outer_grp[j]], py)
      }
     }

  for(i in 1:Ninner){ alpha_in[i] ~ dnorm(0, taua) }
  for(i in 1:Nouter){ alpha_out[i] ~ dnorm(0, taub) }
  beta ~ dnorm(0, 0.001)
  taua ~ dgamma(0.01, 0.01)
  taub ~ dgamma(0.01, 0.01)
  py ~ dgamma(0.01, 0.1)
}"

1 Ответ

0 голосов
/ 13 июля 2020

Можно добавить перехват внешней группы, используя вложенную индексацию, при этом сохраняя формат l oop. Я буду использовать набор данных Pastes из lme4 для сравнения. 

filestring <-
  "model{ 
     for(j in 1:Ninner){
      for(i in 1:N){
        y[j,i] ~ dnorm(beta + alpha_in[j] + alpha_out[batch[j]], py)
      }
     }

  for(i in 1:Ninner){ alpha_in[i] ~ dnorm(0, taua) }
  for(i in 1:Nouter){ alpha_out[i] ~ dnorm(0, taub) }
  beta ~ dnorm(0, 0.001)
  taua <- 1/(sa*sa)
  sa ~ dunif(0,100)
  taub <- 1/(sb*sb)
  sb ~dunif(0,100)
  py ~ dgamma(0.001, 0.001)
}"
INITS <- list(list(.RNG.seed=1, .RNG.name="base::Wichmann-Hill"), 
              list(.RNG.seed=2, .RNG.name="base::Wichmann-Hill"))
results <- run.jags(filestring, monitor=c("py", "beta", "alpha_in", "alpha_out", "sa", "sb"), 
                    data=jags_data, burnin=1e4, sample=1e4, n.chains=2, 
                    inits=INITS, summarise=0)
summary(results, vars=c("py", "beta", "sa", "sb"))

Сравните с lme4

fm1 <- lmer(strength ~ (1|batch) + (1|sample), Pastes)
print(summary(fm1), corr=FALSE)

Используемые данные 

library(lme4); library(data.table); library(runjags)
data(Pastes); setDT(Pastes)
Pastes[,time := sequence(.N), by=sample]

# Change format to match question
wdat = dcast(Pastes, batch + sample ~ time, value.var="strength")
jags_data = list(y=as.matrix(wdat[,3:4]), batch=wdat$batch, N=2, Ninner=length(unique(wdat$sample)), Nouter=length(unique(wdat$batch)))
...