如何指定嵌套模型

Question

我正在使用 runjags 对一些分层数据建模。我可以对层次结构的一个级别进行建模，但我不知道如何将其扩展到更多级别。我正在尝试使用 Nicky Best 等人使用嵌套循环（而不是嵌套索引）的“使用 WinBUGS 的贝叶斯分层建模”第 24 页的 方法 3 来执行此操作。

我可以建模一个级别

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

但收到错误

Compilation error on line 5.
Attempt to redefine node y[1,1]

我如何扩展它以对第一个嵌套的另一个级别进行建模？谢谢。

下面是一些可重现的数据，显示了数据的结构。我想估计 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)
}"

Answer 1

可以通过使用嵌套索引添加外部组截距，同时仍然使用循环格式。我将使用来自 lme4 的 Pastes 数据集进行比较。

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

如何指定嵌套模型

How to specify nested model

r

jags

runjags