[精彩WinBUGS答问]Two models in one Winbugs script

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

I am conducting a Bayesian analysis using Winbugs from R. I need to combine two Winbugs scripts into one: however, I am receiving an error message (Variable x2 is not defined in model or in data set). Here is the winbugs code:

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

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

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

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

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

  }
where x2 is a 1002*2 matrix and y is a vector

This is the R code definining the data:

combined.data <- list(n=n,tto=tto,concern2=concern2,
                       concern3=concern3,y2=y2, x2=x2)
Anyone know what is wrong?

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：winbugs WINBUG script models model receiving message however error

I'm going to be making quite a few assumptions here...

Perhaps you could add a diagram illustrating the relationships between the variables, and which are deterministic vs stochastic. I find this helpful when making models in BUGS. Also, it would be helpful to have the dimensions of all your data, the meaning of n and perhaps some context or detail on what you're modelling and the nodes in which you're interested.

I'm guessing that y is a binary (0,1) vector of length 1002, and has corresponding values for x2[,1] and x2[,2] (herein x1, x2) and concern2, concern3 (herein c2, c3) and tto i.e.

nrow(x2) == 1002
Here's some sample data with of nrow==10 to work with:

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)
It appears that you're trying to determine the values of beta.concern2 (herein b2) for both values of x2 in the logit. Not sure why you'd want to fit it with the same parameter for two different predictors. In case this is a typo I'm giving b2 and b3 as parameters instead. I hope you'll be able to adapt this to your needs.

The product of these values of b2, b3 (stochastic) and c2, c3 (given) are used to generate a variable mu, which also has an error term. (I'm presuming b[i] (herein b1[i]) is a normally distributed error term.) Then tto is a normally distributed variable which depends on the value of mu, and itself has an error term. I have set the precision of the error terms as being equal in both cases.

So for such a model:

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