在 R2WinBUGS 中为 n 个数据集运行模拟研究答案

【问题标题】：Run simulation studies in R2WinBUGS for n datasets在 R2WinBUGS 中为 n 个数据集运行模拟研究
【发布时间】：2018-03-31 01:33:09
【问题描述】：

我在弄清楚如何使用 R2WinBUGS 运行一些模拟研究时遇到问题。目的是模拟 n 个数据集（以 1000 个为目标，但从 10 个开始），并将它们全部作为矩阵放入 R2WinBUGS 代码中，以便当它移植到 WinBUGS 时，它将运行生成 n 个数据集的估计值。这是我目前拥有的：

模型：

model{
      alpha0 ~ dnorm(66.6, 0.01)
      alpha1 ~ dnorm(0.3, 0.01)
      alpha2 ~ dnorm(100, 0.01)
      alpha3 ~ dnorm(0.2, 0.01)
      beta0 ~ dnorm(35, 0.01)
      beta1 ~ dnorm(80, 0.01)
      tau ~ dgamma(0.3,1)

    for(k in  1:Ndat) {
      y[k] ~ dnorm(mu[k], tau)
      mu[k] <- ((alpha0/(1 + exp(-alpha1*(28-beta0)))) + (alpha2/(1 + exp(-alpha3*(28-beta1)))))
    }
}

我使用的错误代码是：

grapedat.sim = bugs(data = list('Ndat' = Ndat, 'y' = p.y[,1]),inits,
                model.file="H:/R coding/R2WinBUGS/multsimt1.bug",
             parameters=c("alpha0","alpha1","alpha2","alpha3","beta0","beta1","tau"),
                n.chains=1,n.iter=8000,n.sim = 6000, 
n.burnin=2000,n.thin=1,
                bugs.directory="H:/WinBUGS14",
                codaPkg=FALSE,
                debug = T)

其中 Ndat 是数据集的数量，p.y。是一个 13 x n 矩阵，初始值是：

inits <- function(){
  list(alpha0=70, alpha1=0.4, tau=0.15, alpha2=105, alpha3=0.25,beta0 = 
    40, beta1 = 85)
  }

有什么帮助吗？

【问题讨论】：

标签： r winbugs r2winbugs

【解决方案1】：

理论上，您可以在单个 BUGS 模型中执行此操作，方法是使用第二个（外部）循环，然后将 y 和 mu 索引为矩阵，并将系数作为向量，例如（未经测试）：

model{
    for(d in 1:n){
      alpha0[d] ~ dnorm(66.6, 0.01)
      alpha1[d] ~ dnorm(0.3, 0.01)
      alpha2[d] ~ dnorm(100, 0.01)
      alpha3[d] ~ dnorm(0.2, 0.01)
      beta0[d] ~ dnorm(35, 0.01)
      beta1[d] ~ dnorm(80, 0.01)
      tau[d] ~ dgamma(0.3,1)

      for(k in  1:Ndat) {
        y[k,d] ~ dnorm(mu[k,d], tau[d])
        mu[k,d] <- ((alpha0[d]/(1 + exp(-alpha1[d]*(28-beta0[d])))) +
                    (alpha2[d]/(1 + exp(-alpha3[d]*(28-beta1[d])))))
      }
    }
}

这对于 n=10 是可行的，但对于 n=1000 会非常笨拙，因此通常不是解决问题的好方法。相反，我会在 R 中使用循环独立运行每个数据集/模型。如果您愿意使用 JAGS 而不是 WinBUGS，那么您可以在 runjags 包中查看半自动解决方案 - 即阅读本文的第 2.9 节：https://www.jstatsoft.org/article/view/v071i09

马特

【讨论】：