2

我正在使用 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)

有谁知道出了什么问题？

4

1 回答 1

1

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

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

我猜这是一个长度为 1002的二进制( 0,1 )向量，并且有对应的y值x2[,1]x2[,2]x1x2concern2concern3c2c3tto

nrow(x2) == 1002

以下是一些与 of 一起使用的示例数据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 中的两个值确定beta.concern2(herein b2) 的值。x2不知道为什么要为两个不同的预测变量使用相同的参数。如果这是我给出的错字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))])
于 2013-01-25T22:42:22.990 回答