Bayes Factor in R with Jaggs

Question

I am trying to decide between linear regression models using the Bayes Factor in Jaggs through R. For simplicity's sake, I am working on the mtcars dataset and the models are:

(1) mpg on only intercept (2) mpg on wt (3) mpg on disp (4) mpg on wt and disp

I am utilizing flat priors and try to implement the model comparison with the following code, but run into some problems and errors:

reg.bf.model <- "model{
    M ~ dcat(model.p[])
    model.p[1] <- 0.25
    model.p[2] <- 0.25
    model.p[3] <- 0.25
    model.p[4] <- 0.25

  # Likelihood
  for (i in 1:n){
      y[i] ~ dnorm(mu[i,M],tau[M])
        mu[i,1] <- alpha[1]
        mu[i,2] <- alpha[2] + beta1[1] * wt[i]
        mu[i,3] <- alpha[3] + beta1[2] * disp[i]
        mu[i,4] <- alpha[4] + beta1[3] * disp[i]+ beta2* wt[i]
  }

  for (m in 1:4){
      alpha[m] ~ dnorm(0,0.0001)
      sigma2[m] <- (1/tau[m])
  }

  for (j in 1:3) {
      beta1[j] ~ dnorm(0,0.0001)
  }

    beta2 ~ dnorm(0,0.0001)

    tau[1] ~ dgamma(prior.shape1[M],prior.rate1[M])
    tau[2] ~ dgamma(prior.shape2[M],prior.rate2[M])
    tau[3] ~ dgamma(prior.shape3[M],prior.rate3[M])
    tau[4] ~ dgamma(prior.shape4[M],prior.rate4[M])
}"

jags.data <- list(y=mtcars$mpg,n=nrow(mtcars), 
                  wt = mtcars$wt, disp = mtcars$disp,
                  prior.shape1=c(0.0001,1),prior.rate1=c(0.0001,1),
                  prior.shape2=c(1,0.0001),prior.rate2=c(1,0.0001),
                  prior.shape3=c(1,0.0001),prior.rate3=c(1,0.0001),
                  prior.shape4=c(1,0.0001),prior.rate4=c(1,0.0001))

jags.bf <- jags.model(file=textConnection(reg.bf.model),
                      #                      inits=jags.inits,
                      data=jags.data, n.chains=3)

update(jags.bf, 100)
jags.bf.out <- coda.samples(jags.bf,
                            variable.names=c("alpha","beta1","beta2"
                                             ,"sigma2","M"),
                            n.iter=50000, thin=400)

But I get the following error messages:

"Error in update.jags(jags.bf, 100) : LOGIC ERROR: SimpleRange::leftOffset. Index outside of allowed range Please send a bug report to martyn_plummer@users.sourceenter code hereforge.net"

and

"Error in jags.samples(model, variable.names, n.iter, thin, type = "trace", : No valid monitors set In addition: Warning messages: 1: In FUN(X[[i]], ...) : Can't set monitor. No model!

2: In FUN(X[[i]], ...) : Can't set monitor. No model!

3: In FUN(X[[i]], ...) : Can't set monitor. No model!

4: In FUN(X[[i]], ...) : Can't set monitor. No model!

5: In FUN(X[[i]], ...) : Can't set monitor. No model!"

Any ideas or nudges in the right direction are greatly appreciated!

I think the problem arises from the model parameter `M` and it's usage. I have to admit, while I do have some experience with `rjags` I'm having trouble understanding what the line `M ~ dcat(model.p[])` and the subsequent usage of `M` as an index does. Anyhow, the model param `prior.shape1` has only 2 elements ( in `jags.data`) and the error tells me that it tried to access element [3]... — dario, Feb 18 '20 at 10:45

m1gnoc · Answer 1 · 2020-02-18T14:57:08.690

As an update, apparently the Bayes factor needs to be used as an alternative to the frequentist's hypothesis test between two models. Therefore the following code adjustment got things to work.

reg.bf.model <- "model{
  M ~ dcat(model.p[])
  model.p[1] <- prior1
  model.p[2] <- 1-prior1


  # Likelihood
  for (i in 1:n){
    y[i] ~ dnorm(mu[i,M],tau[M])
    mu[i,1] <- alpha[1]
    mu[i,2] <- alpha[2] + beta1[1] * wt[i]
  }

  for (m in 1:2){
      alpha[m] ~ dnorm(0,0.0001)
      sigma2[m] <- (1/tau[m])
  }

  for (j in 1:1) {
    beta1[j] ~ dnorm(0,0.0001)
  }

  beta2 ~ dnorm(0,0.0001)

  tau[1] ~ dgamma(prior.shape1[M],prior.rate1[M])
  tau[2] ~ dgamma(prior.shape2[M],prior.rate2[M])

}"
prior1 = 0.5
jags.data <- list(y=mtcars$mpg,n=nrow(mtcars), prior1 = prior1,
                  wt = mtcars$wt, disp = mtcars$disp,
                  prior.shape1=c(0.0001,1),prior.rate1=c(0.0001,1),
                  prior.shape2=c(1,0.0001),prior.rate2=c(1,0.0001))

jags.bf <- jags.model(file=textConnection(reg.bf.model),
                      #                      inits=jags.inits,
                      data=jags.data, n.chains=3)

update(jags.bf, 100)
jags.bf.out <- coda.samples(jags.bf,
                            variable.names=c("alpha","beta1","beta2"
                                             ,"sigma2","M"),
                            n.iter=50000, thin=400)

# calculating bayes factor
posterior.M <- unlist(jags.bf.out[,"M"])
posterior.odds <- mean(posterior.M==1)/mean(posterior.M==2)
prior.odds <- prior1/(1-prior1)
bayes.factor <- posterior.odds/prior.odds
bayes.factor

Bayes Factor in R with Jaggs

1 Answers1