Reading in model in JAGS

2012-11-23
2012-11-28
  • Hi,
    I am trying to create a model in JAGS where my process model is based on the output of a large process model. I have a table of output from this model that I want to combine in a certain way before comparing the results with data. I have tried to do this by reading in the table as data and then do the manipulations in JAGS. However, I do not manage to reproduce known parameters when I test with synthetic data. Is it possible to read in the model as data in JAGS without errors? I see in the manual (par 2.1.2) that "it is an error to supply a data value for a deterministic node". Not sure if this means that it is not possible to do what I outline above. If it is not possible are there workarounds?

    Cheers,
    Øystein

     
  • Martyn Plummer
    Martyn Plummer
    2012-11-28

    It sounds like you are talking about data transformation. You want to read in some data, transform it in a deterministic way and then model the transformed data, e.g. you read in skewed data x then log-transform it to a variable y and model y as Gaussian.

    In JAGS, this can be done by putting the parameter transformation in a separate data block

    data {
       for (i in 1:length(y)) {
          y[i] <- log(x[i])
       }
    }
    model {
       for (i in 1:length(y)) {
          y[i] ~ dnorm(mu, tau)
       }
       mu ~ dnorm(0, 1.0E-3)
       tau ~ dgamma(0.001, 0.001)
    }
    

    See the User manual for details.

     
  • Hi,
    Unfortunately I don't think it is that simple. Because I need access to the parameter values (priors) when manipulating the data (which are my model somehow). Here is my code trying to do so, but it does not return right values for synthetic dataset.
    hits contain the model output from a big process model run, and it is this I want to manipulate to model the log-transformed data (ld). dist[1-4] are the parameters I want to estimate.

    model{
    sigma.o~dunif(0,5)
    tau.o <- 1/sigma.o^2
    mu<-30
    std<-14.0
    
    dist1~dunif(0,1000)
    dist2~dunif(0,1000)
    dist3~dunif(0,1000)
    dist4~dunif(0,1000)
    
    scale<-1.0
    
    dist[1] <- dist1
    dist[2] <- dist1
    dist[3] <- dist1
    dist[4] <- dist1
    dist[5] <- dist1
    dist[6] <- dist1
    dist[7] <- dist1
    dist[8] <- dist1
    dist[9] <- dist2
    dist[10] <- dist2
    dist[11] <- dist2
    dist[12] <- dist2
    dist[13] <- dist3
    dist[14] <- dist3
    dist[15] <- dist3
    dist[16] <- dist4
    dist[17] <- dist4
    
    mort.egg <- 0.17
    
    for(i in 1:nh) {  
      gg[i] <- dnorm(hits[i,6],mu,std)
      surv[i] <- exp(-mort.egg*hits[i,3])
      fract[i] <- dist[hits[i,4]+1]  #*ssb[hits[i,2]-1958,2]
    
      sim[i,1] <- gg[i]*surv[i]*fract[i]*scale*equals(hits[i,7],1)
      sim[i,2] <- gg[i]*surv[i]*fract[i]*scale*equals(hits[i,7],2)
      sim[i,3] <- gg[i]*surv[i]*fract[i]*scale*equals(hits[i,7],3)
      sim[i,4] <- gg[i]*surv[i]*fract[i]*scale*equals(hits[i,7],4)
    }
    
    for(k in 1:nn) {
        for(j kn 1:4){
          td[k,j] <- sum(sim[(nrh[k,1]+1):nrh[k+1,1],j])
          ld[k,j]~dnorm(log(td[k,j]+1),tau.o)
        }
      }
    }
    

    This code is runned with r2jags using the following commands in r:
    data <- read.table('data.txt')
    "dd"<-as.matrix(data)[,6:9]
    "ld" <- log(dd+1)
    "nn" <- length(dd[,1])
    "hits"<-as.matrix(read.table('hits.txt'))
    "nh" <- length(hits[,1])
    "nrh" <- as.matrix(read.table('nr_hits.txt'))
    calinit <- list(list(mu=30,std=14.0,dist1=0.1,dist2=100,dist3=50,dist4=0.1,scale=1),
    list(mu=30,std=14.0,dist1=0.1,dist2=100,dist3=50,dist4=0.1,scale=1),
    list(mu=30,std=14.0,dist1=0.1,dist2=100,dist3=50,dist4=0.1,scale=1)
    )
    caldat <- list("ld"=ld,"nn"=nn,"hits"=hits,"nh"=nh,"nrh"=nrh)
    calmod <- "calibrate_model.r"
    calpara <- c("dist1","dist2","dist3","dist4","sigma.o")
    set.seed(123456)

    calfit <- jags(data=caldat, inits=calinit, parameters.to.save=calpara,         jags.seed=1234,model=calmod, n.chains=3, n.iter=3000,n.burnin=1000,n.thin=10,refresh = 10)
    

    I am very happy for any feedback.

    Cheers,
    Øystein