0

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

对于一个级别，我可以建模

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]

如何扩展它以模拟第一个嵌套的另一个级别？谢谢你。

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

1 回答 1

0

可以使用嵌套索引添加外部组拦截，同时仍使用循环格式。我将使用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)))
于 2020-07-12T22:05:24.543 回答