Warning message: Failed to set trace monitor for gamma1 Node not found

Help
2014-03-26
2014-03-28
  • Nicolas Fasel

    Nicolas Fasel - 2014-03-26

    Hi every one,

    I have the exact question as alon.benari:

    http://r.789695.n4.nabble.com/Error-in-JAGS-cannot-monitor-z-td4333685.html

    I am also trying to adapt my model to the example given by Zuur et al. 2012.

    My code:

     Noff.ij <- VecMat(factor(total$ID),total$OFFSPRING)
     Tobs.ij <- VecMat(factor(total$ID),total$Tobs)
     Statut.ij<- VecMat(factor(total$ID), total$STATUT)
     iStatut.ij <-Statut.ij-1
     NMales <- length(levels(factor( total$ID)))
     Nobs.ij <- as.numeric(tapply(total$OFFSPRING,FUN=length,INDEX= total$ID))
     win.data <- list(Noffspring = Noff.ij, Statut= iStatut.ij, NMales= NMales, Tobs=Tobs.ij,Nobs.ij = Nobs.ij)
     sink("modelglmm.txt")
     cat("
    model{
    #Priors
    beta1~dnorm(0,0.001)
    for(i in 1:NMales){a[i]~dnorm(0,tau.Male)}
    alpha~dnorm(0,0.001)
    gamalpha~dnorm(0,0.001)
    gambeta~dnorm(0,0.001)
    tau.Male <- 1/pow(sigma.Male,2)
    sigma.Male~dunif(0,10)
    
    #Likelihood
    for(i in 1:NMales){
        for(j in 1:Nobs.ij[i]){
    
            #logit part
            W[j,i]~dbern(psi.min1[j,i])
            psi.min1[j,i]<- min(0.99999,max(0.00001,(1-psi[j,i])))
            gamma1[j,i] <-gamalpha+gambeta*Tobs[j,i] 
            logit(psi[j,i]) <- max(-20,min(20,gamma1[j,i]))
    
            #Poisson part
            Noffspring[j,i] ~ dpois(eff.mu[j,i])
            eff.mu[j,i] <- W[j,i]*mu[j,i] 
            log(mu[j,i]) <- max(-20,min(20,eta[j,i]))
            eta[j,i] <- alpha+ beta1 * Statut[j,i]+a[i]
    
            #Extract Pearson residuals
            EZip[j,i] <- mu[j,i]*(1-psi[j,i])
            VarZip[j,i] <- (1-psi[j,i])*(mu[j,i]+psi[j,i]*pow(mu[j,i],2))
            PRes[j,i] <- (Noffspring[j,i]-EZip[j,i])/sqrt(VarZip[j,i])
                }
    }
    }
    ")
    sink()
    
    MyW <- Noff.ij
    MyW[Noff.ij>0] <- 1
    # Params monitored
    params <- c("alpha","beta1","gambeta","gamma1","PRes") 
    
    # MCMC settings
    ni <- 3000
    nt <- 1 
    nb <- 500
    nc <- 3 
    
    #inits
     inits <- function(){
        list(alpha=rnorm(1),gambeta=rnorm(1),gamalpha=rnorm(1),beta1=rnorm(1),a=rnorm(length(levels(total$ID)),0,1),sigma.Male=rlnorm(1),W = MyW)}
    
    out <- jags(win.data, inits, params, "modelglmm.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
    

    I did not do long chains in order to save time, utile models run correctly...

    Jags run well, but cannot monitor deterministic parameter such as gamma1 or PRes.

    I do not think that I have missed the determination of any of them...

    Any idea would be very useful!

    Thank you

     
    Last edit: Martyn Plummer 2014-03-27
  • Jonas Lindeløv

    Jonas Lindeløv - 2014-03-26

    A quick guess is that gamma1 is not "filled out" on the j-dimension, i.e. the loop does not go to max(Nobs.ij) every time. JAGS won't return the object if this is the case. I spent quite some time discovering this myself :-)

    There's two solutions to this. The most secure is to do a second "j-loop" after the first, filling in the remaining places.

    I found that just looping to max(Nobs.ij) and having NA's for non-existing points/subjects doesn't change the posterior. This gives simpler code but I cannot recommend it because I do not know if this behavior is guaranteed.

     
  • Nicolas Fasel

    Nicolas Fasel - 2014-03-27

    Sorry to ask that, but if it would not have go until the required Nobs.ij. How could it have succeed to calculate at all? I would have receive a error message.

    Then I did not understand how to do the second j-loop without specifying the upper limit. Or how do you think to fix this upper limit? With the max(Nobs.ij) from your second option?

    thanks for your help

     
  • Martyn Plummer

    Martyn Plummer - 2014-03-27

    Like this:

    for(i in 1:NMales){
        for(j in (Nobs.ij[i]+1):max(Nobs.ij)){
            gamma1[j,i] <- 0
        }
    }
    

    When you monitor gamma1, you will see a lot of zeros corresponding to the elements that you did not define before in your model. These zeros fill out the array and allow you to monitor it.

    In the next version, JAGS should fill out incomplete arrays like this automatically with missing values when you want to monitor them.

     
  • Nicolas Fasel

    Nicolas Fasel - 2014-03-28

    Great! Thank you

     

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks