一个 Winbugs 脚本中的两个模型答案

【问题标题】：Two models in one Winbugs script一个 Winbugs 脚本中的两个模型
【发布时间】：2012-12-25 13:04:43
【问题描述】：

我正在使用 R 中的 Winbugs 进行贝叶斯分析。我需要将两个 Winbugs 脚本合并为一个：但是，我收到一条错误消息 (Variable x2 is not defined in model or in data set)。这是winbugs代码：

model{
# Model’s likelihood
 for (i in 1:n) {
    tto[i] ~ dnorm( mu[i], tau ) # stochastic componenent
    b[i] ~ dnorm(0.0, tau2)
    # link and linear predictor
    mu[i] <- 1 - (beta.concern2*concern2[i] + beta.concern3*concern3[i] + b[i])
 }

 for (i in 1:1002) {
    # Linear regression on logit
    logit(p[i]) <- beta.concern2*x2[i,1] + beta.concern2*x2[i,2]

    # Likelihood function for each data point
    y2[i] ~ dbern(p[i])
 }

  s2<-1/tau
  s <-sqrt(s2)

  a2<-1/tau2
  a <-sqrt(a2)

  }

x2 是一个 1002*2 矩阵，y 是一个向量

这是定义数据的 R 代码：

 combined.data <- list(n=n,tto=tto,concern2=concern2,
                       concern3=concern3,y2=y2, x2=x2)

有人知道怎么回事吗？

【问题讨论】：

@Giulio，我认为这个问题是绝对合法的。在大多数 WinBUGS 错误中，您无法进行任何研究。只是尝试和失败。他为什么不缩短刑期并寻求帮助？这就是 stackoverflow 的用途。
@Tomas，非常感谢。你对我很好。我真的尽力了，winbugs对新用户不是很友好。
你能把得到的 combine.data 转储到这里吗？尤其是 n 常量。如果 x2 充满数字（没有 NA）和正确的尺寸，这一点很重要。
@Guilio，当然 - 这对任何错误都有效 - 只需修改您的代码并更正它，仅此而已:-) Stackoverflow 很好，但我认为更有经验的人对其他人的款待可以得到改进。
@James：JAGS 和 WinBUGS 之间并没有那么很大的区别，尽管 JAGS 的错误消息无疑更清晰一些。请阅读 tinyurl.com/reproducible-000 关于构建一个可重现的示例 ...您还没有向我们展示 (1) 您的数据是什么，(2) 您如何将 combined.data 传递给 WinBUGS 模型... （例如，您使用的是 R2WinBUGS 吗？精确的 bugs() 函数调用是什么？）

标签： r winbugs

【解决方案1】：

我将在这里做很多假设......

也许您可以添加一个图表来说明变量之间的关系，这些关系是确定性的还是随机的。在 BUGS 中制作模型时，我发现这很有帮助。此外，了解所有数据的维度、n 的含义以及有关您正在建模的内容和您感兴趣的节点的一些上下文或详细信息会很有帮助。

我猜y 是一个长度为 1002 的二进制 (0,1) 向量，并具有 x2[,1] 和 x2[,2] 的对应值（此处为 x1、x2）和 concern2 、concern3（此处为c2、c3）和tto，即

nrow(x2) == 1002

这里有一些与nrow==10 一起使用的示例数据：

y <- sample(x=c(0,1), size=10, replace=TRUE, prob=c(0.5,0.5))
x2 <- matrix(rnorm(20), nrow=10, ncol=2)
c2 <- rnorm(10)
c3 <- rnorm(10)
tto <- rnorm(10)

您似乎正在尝试确定 logit 中 x2 的两个值的 beta.concern2（此处为 b2）的值。不知道为什么要为两个不同的预测变量使用相同的参数。如果这是一个错字，我将提供 b2 和 b3 作为参数。我希望您能够根据自己的需要进行调整。

b2、b3（随机）和c2、c3（给定）的这些值的乘积用于生成变量mu，该变量也有一个错误项。（我假设b[i]（此处为b1[i]）是一个正态分布的误差项。）那么tto是一个正态分布的变量，它依赖于mu的值，并且它本身有一个误差项。我已将两种情况下的误差项的精度设置为相等。

所以对于这样的模型：

require(rjags)

### The data
dataList <- list(
    x1 = x2[,1],
    x2 = x2[,2],
    y = y,
    c2 = c2,
    c3 = c3,
    tto = tto,
    nRowX = nrow(x2)
    )

### make sure logistic model can be fitted
f1 <- stats::glm(dataList$y ~ dataList$x1 + dataList$x2 -1, family=binomial(logit))
show(f1)

### set some approximate initial values
b1Init <- 0.1 # arbitrary
b2Init <- f1$coef[2]
b3Init <- f1$coef[3]

initsList <-  list(
 b1 = b1Init,
 b2 = b2Init,
 b3 = b3Init)

### Model: varying parameters (b2, b3) per observation; 2x error terms
modelstring <- "
    model {
        for(i in 1:nRowX){
            tto[i] ~ dnorm(mu[i], prec)
            mu[i] <- 1 - (b1 + b2*c2[i] + b3*c3[i])
            y[i] ~ dbern(L[i]) # L for logit
            L[i] <- 1/(1+exp(- ( b2*x1[i] + b3*x2[i]) ))
        }
        b1 ~ dnorm(0, prec) # precision
        prec <- 1/sqrt(SD) # convert to Std Deviation
        SD <- 0.5
        b2 ~ dnorm(0, 1.4) # arbitrary
        b3 ~ dnorm(0, 1.4)
    }
"

writeLines(modelstring,con="model.txt")

parameters <- c("b1","b2","b3") # to monitor
adaptSteps <- 1e4 # "tune in" samplers
burnInSteps <- 2e4 # "burn in" samplers
nChains <- 3
numSavedSteps <-2e3
thinSteps <- 1 # Steps to "thin" (1=keep every step).
nPerChain <- ceiling(( numSavedSteps * thinSteps ) / nChains) # Steps per chain

rm(jagsModel) # in case already present

jagsModel <- rjags::jags.model(
 "model.txt", data=dataList,
 inits=initsList, n.chains=nChains,
 n.adapt=adaptSteps)

stats::update(jagsModel, n.iter=burnInSteps)

### MCMC chain
MCMC1 <- as.matrix(rjags::coda.samples(
 jagsModel, variable.names=parameters,
 n.iter=nPerChain, thin=thinSteps))

### Extract chain values
b2Sample <- as.vector(MCMC1[,grep("b2",colnames(MCMC1))])

【讨论】：