CODEWITHSUNDEEP
rjagsrjags
Me28
Me28

Reputation: 107

jags priors in for loop

I'm trying to minimize my jags code by work with for loop. my original code is

data<- list(r1=c(16, 62, 14, 23, 570, 63, 63, 116),
            r2=c(10, 66, 20, 27, 522, 31, 31, 95),
            n1=c(53, 106, 56, 82, 1012, 201, 201, 2049),
            n2=c(53, 107, 94, 82, 1081, 149, 149, 3097))
NS <- length(data$r1) 
library(rjags) 
library(coda)

ma <- function(o){
  out <- "
  model{
  for(i in 1:NS){
  delta[i,1] <- 0
  mu[i] ~ dnorm(0, 0.0001) 
  for(k in 1:2){
  r[i,k] ~ dbin(p[i,k], n[i,k])  
  logit(p[i,k]) <- mu[i] + delta[i,k]  
  }
  delta[i,2] ~ dnorm(lor, prec)  
  }
  lor ~ dnorm(0, 0.0001)  
  ############## 1  inverse Gamma Dist.
   tau2 <- 1/prec
   taus <- sqrt(1/prec)

   prec ~ dgamma(0.1,0.1)
  # prec ~ dgamma(0.01,0.01)
  # prec ~ dgamma(0.001,0.001)

  ############## 2 Uniform Dist.
  # tau2 <- tau*tau
  # prec <- 1/tau2
  # taus <- tau

  # tau ~ dunif(0, 2)
  # tau ~ dunif(0, 5)
  # tau ~ dunif(0, 10)

  ### Transform ln(OR) into OR
  OR <- exp ( lor )
  }"
  return(out)
}

as you can see I run each prior manually instead of that i would like to run them in a for loop my try;

# hyper-parameters
IG <- c(a1=0.1, a2=0.01, a3=0.001,b1=0.1, b2=0.01, b3=0.001)
U <- c(a1=0, a2=0, a3=0, b1=2, b2=5, b3=10)
# priors
prior <- data.frame(IG,U)

if (prior == "IG"){
  for (i in 1:3){
    tau2[i] <- 1/prec[i]
    taus[i] <- sqrt(1/prec[i])
    prec[i] ~ dgamma(a[i],b[i])}
} 

if (prior == "U"){
  for (i in 1:3){
    tau2[i] <- tau[i]*tau[i]
    prec[i] <- 1/tau2[i]
    tau[i] ~ dunif(a[i],b[i])}
}

is there any way to do it in rjags, because im not sure if there is an if statement in JAGS software.

Upvotes: 0

Views: 579

Answers (1)

mfidino
mfidino

Reputation: 3055

There really are not any if statements that works that way in JAGS, and you will have an error if you are trying to create two tau, taus, and prec objects. There would be two ways to do this. The way I prefer would be to tack on the parameter name to the object. By this I mean:

for(i in 1:3){
    # IG priors
    tau2_IG[i] <- 1/prec_IG[i]
    taus_IG[i] <- sqrt(1/prec_IG[i])
    prec_IG[i] ~ dgamma(a_IG[i],b_IG[i])
    # U priors
    tau2_U[i] <- 1/prec_U[i]
    taus_U[i] <- sqrt(1/prec_U[i])
    prec_U[i] ~ dgamma(a_U[i],b_U[i])
}

The other way to do it would be to just determine the total number of priors you need for both of these (which in this example would be 6). However, that requires you to make sure you input them in the correct place throughout the model (and then you have to remember what parameters tau[1:3] are tied to compared to tau[1:6].

I cannot really see how you plan to relate your second code block (with all of your hyperparameters) back to the JAGS model (which only has 1 prec prior), but hopefully this example helps.

Upvotes: 0

Related Questions