"convergence" for a derived quantity in JAGS/R2Jags

Question

UPDATE: Now with Traceplot example

UPDATE: Now with new traceplot

I am trying to adapt Outhwaite et. als 2018 code for occupancy modelling and have a couple of questions that I just can't seem to find an answer for...

Code used to create model

cat(
  "model{ 
  
  ### Model ###

        # State model
        for (i in 1:nsite){ 
          for (t in 1:nyear){   
            z[i,t] ~ dbern(psi[i,t]) 
            logit(psi[i,t])<- b[t] + u[i] 
          }}   
        
        # Observation model 
        for(j in 1:nvisit) {
          y[j] ~  dbern(Py[j]+0.0001)
          Py[j]<- z[Site[j],Year[j]]*p[j]      
          logit(p[j]) <- a[Year[j]] + c*logL[j]
        }
        
   ### Priors ###

      # State model priors
          for(t in 1:nyear){
            b[t] ~ dunif(-10,10)        # fixed year effect
          }                 
          
          for (i in 1:nsite) {
            u[i] ~ dnorm(0, tau.u)  # random site effect      
          } 
          
          tau.u <- 1/(sd.u * sd.u)  
          sd.u ~ dunif(0, 5)        # half-uniform hyperpriors
          
        
      # Observation model priors 
          for (t in 1:nyear) {
            a[t] ~ dnorm(mu.a, tau.a)   # random year effect           
          }
          
          mu.a ~ dnorm(0, 0.01)                         
          tau.a <- 1 / (sd.a * sd.a)    
          sd.a ~ dunif(0, 5)        # half-uniform hyperpriors             
          
          c ~ dunif(-10, 10)        # sampling effort effect
  
   
  ### Derived parameters ###
  # Finite sample occupancy - proportion of occupied sites
  for (t in 1:nyear) {  
    psi.fs[t] <- sum(z[1:nsite,t])/nsite
  } 
  
  #data# nyear, nsite, nvisit, y, logL, Site, Year
  
}", file="bmmodel.txt"
)

Note that dbern(Py[j]+0.0001) includes a correction factor since dbern(0) is not supported in JAGS.

I am running the model on some plant data just basically trying it out to see if it runs and converges and behaves as I would expect it to.

Question number 1(ANSWERED): I am interested in the quantity psi.fs[t]. But since the model calculates this quantity after the actual modelling process, can convergence be assessed for psi.fs[t]?

R code for running model with R2JAGS

jagsrespsi<-jags(data.list, inits=test.inits,
     n.chains=2, n.iter=15000, n.thin=3,
     DIC=T, 
     model.file=paste0(modeltype,"model.txt"), parameters.to.save=c("psi.fs"))

Question number 2: When I use traceplot(jagsrespsi) to plot the traceplot seems all over the place but the Rhat for jagsrespsi$BUGSoutput is 1 for all my years? gelman.diag(as.mcmc(jagsrespsi)) also indicates convergence. Same goes for monitoring psi!

I am very astonished by this model behaviour and am suspecting there is something wrong... but no idea where to look

Answer 1

Yes, you can check psi.ft[] for convergence in exactly the same way as you check the convergence of the model's parameters. That's exactly what happens, for example, in a logistic regression, where the fitted probabilities of response are calculated as exp(z)/(1 + exp(z)) for some linear predictor z.

When you say the traceplot is "all over the place", what do you mean? This could be either good or bad. Can you show an example? A "good" traceplot looks like a "fat, hairy caterpillar": consecutive samples taken from all regions of the sample space, a horizontal hair ball. Although written for SAS, this page gives a reasonable high level description of what a good trace plot looks like, and what problems might be indicated by less-than-ideal examples.

In response to your edit to include the trace plot...

That doesn't look like a particularly good traceplot to me: there seems to be some negative autocorrelation between successive samples. Have you calculated the effective sample size [ESS]?

But the plot may look a little odd because your chain is very short, IMHO. You can use the ESS to provide a very rough approximation for the accuracy of an estimated probability: the worst case half width CI of a binomial proportion is +/-2 * sqrt(0.5*0.5/N), where N is the sample size (or ESS in this case). So even if the efficiency of your MCMC process is 1 - so that the ESS is equal to the chain length - then the accuracy of your estimates is only +/-0.02. To estimate a probability to 2 decimal places (so that the half width of the CI is no more than 0.005), you need an ESS of 40,000.

There's nothing wrong with using short chain lengths during testing, but for "production" runs, then I would always use a chan length much greater than 2,500. (And I'd also use multiple chains so that I can use Gelman-Rubin statistics to test for convergence.)