Two models in one Winbugs script

Question

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

  }

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?

Answer 1

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