WAIC calculation for CJS model in JAGS

[Dilsad Dagtekin]

Hi everyone!

I am trying to fit a CJS model in JAGS using jagsUI package.

I have several models that I would like to compare with WAIC.

For that, I am trying to calculate a log-likelihood inside the model and take it out as a parameter, so I can use the loo package to have a WAIC value. 

However, I encounter with this error for my log likelihood calculation: "Unable to resolve the following parameters:

log_lik0[1,1] (line 114)

log_lik0[1,2] (line 114)

log_lik0[1,3] (line 114)

........ " (rest is the same as well). 

Below is my full JAGS code, I got lost in the many forums to find a solution however still couldn't manage the code to work. Does someone have any idea on what am I doing wrong or what is missing?

I'd appreciate any help!

Thank you so much in advance.

Dilsad

------------------ JAGS CODE------------------ J

cat(file = "p1global.txt"," 

model {

## Priors and constraints

#########################

for (i in 1:n.individuals){

  for (t in f[i]:(n.occasions-1)){

  

    logit(phi[i,t]) <- alphaphi + 

                        betaphi1[sex[i]] + betaphi2[agemat[i,t]] + betaphi3[season[t]] + 

                        betaphi4[sex[i],agemat[i,t]] + betaphi5[sex[i],season[t]] + betaphi6[agemat[i,t],season[t]]

                        

    logit(p[i,t]) <- alphap + 

                        betap1[sex[i]] + betap2[agemat[i,t]] + betap3[season[t]] + 

                        betap4[sex[i],agemat[i,t]] + betap5[sex[i],season[t]] + betap6[agemat[i,t],season[t]] +

                        epsilonp[year[t]]

    

  } # t

} # i

### Intercepts ###

alphaphi <- logit(mean.phi)

mean.phi ~ dunif(0, 1)                        # Prior for mean survival

alphap <- logit(mean.p)

mean.p ~ dunif(0, 1)                          # Prior for mean recapture

### Slopes ###

for (b in 1:2){

  

  # Prior for the slope of each sex class, F=1 , M=2

  betaphi1[b] ~ dunif(0,1)

  betap1[b] ~ dunif(0,1)

  

  # Prior for the slope of each age class, J=1, A=2

  betaphi2[b] ~ dunif(0,1)

  betap2[b] ~ dunif(0,1)

  

  # Prior for the slope of each  season, D=1, W=2

  betaphi3[b] ~ dunif(0,1)

  betap3[b] ~ dunif(0,1)

}

# Prior for the slope of the interaction effects

# order for betas always: 1,1; 1,2; 2,1; 2,2

for (i in 1:2){

  for (t in 1:2) {

  

    # sex:age will have total 4, FJ, FA, MJ, MA

    betaphi4[i,t] ~ dunif(0,1)

    betap4[i,t] ~ dunif(0,1)

    

    # sex:season will have total 4, FD, FW, MD, MW

    betaphi5[i,t] ~ dunif(0,1)

    betap5[i,t] ~ dunif(0,1)

    

    # age:season will have total 4, JD, JW, AD, AW

    betaphi6[i,t] ~ dunif(0,1)

    betap6[i,t] ~ dunif(0,1)

  }

}

## Random effect ##

for (t in 1:(n.years)){

    epsilonp[t] ~ dnorm(0, taup)

}

sigmap ~ dunif(0, 10)                     # Prior for standard deviation

taup <- pow(sigmap, -2)

sigma2p <- pow(sigmap, 2)                  # Temporal variance

## Likelihood 

#########################

for (i in 1:n.individuals){

  # Define latent state at first capture

  z[i,f[i]] <- 1

  

  for (t in (f[i]+1):n.occasions){

    # State process

    z[i,t] ~ dbern(mu1[i,t])

    mu1[i,t] <- phi[i,t-1] * z[i,t-1]

    

    # Observation process

    y[i,t] ~ dbern(mu2[i,t])

    mu2[i,t] <- p[i,t-1] * z[i,t]

  } # t 

} # i

  

## Log-Likelihood for WAIC ##

for (i in 1:n.individuals){

  for (t in (f[i]+1):n.occasions){

  

    # estimate the cell log-likelihood, all distributions in JAGS have a log density function associated with them

    log_lik0[i,t] <- logdensity.bern(y[i,t],p[i,t-1] * z[i,t])

    # get the row log-likelihood

    log_lik[i] <- sum(log_lik0[i,])

    

  } # t 

} # i

}")

Later I want to use log_lik with the loo package

library(loo)

waic(p1global_out$sims.list$log_lik)

[Marc Kery]

Dear Dilsad,

that's a good one ! Would love to be able to use the loo for model selection when fitting models in JAGS. Two perhaps only partly helpful comments:

• I think that in Nimble there is directly an automatic node for the log density so that might help ? BUGS model definition language is of course almost the same as in JAGS.
  
• And for JAGS, I think that the log density is not so easy ... I think that you need to specify the real marginal likelihood of the data. Which in your model would involve integrating out the random effects.  Something for a statistician....

Best regards --- Marc

---

From: hmec...@googlegroups.com <hmec...@googlegroups.com> on behalf of Dilsad Dagtekin <dils...@gmail.com>
Sent: Wednesday, November 24, 2021 15:28
To: hmecology: Hierarchical Modeling in Ecology <hmec...@googlegroups.com>
Subject: WAIC calculation for CJS model in JAGS

 

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2020)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/05a4f3da-4ac2-4cf9-a202-73cc3ec04b6en%40googlegroups.com.

[Ben Goldstein]

Hi Dilsad and all,

I think NIMBLE is a good option for you here. NIMBLE has an option for automatic WAIC estimation that was just improved in efficiency in the most recent major update. Also, we have an auxiliary R package, nimbleEcology, which provides a CJS distribution (dCJS) used for marginalizing over the latent state z, saving some memory and in some cases improving mixing time.

Happy to help if this is something you want to learn more about or if you have follow-up questions.

Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/ZR0P278MB0869A200EC39F76B1C24E91EEB619%40ZR0P278MB0869.CHEP278.PROD.OUTLOOK.COM.

[Dilsad Dagtekin]

Thanks all for the suggestions and comments. 

I think I'll try to convert to nimble since it looks like the most logical thing to do for the future as well I suppose. :)

Best,

Dilsad

--

A. Dilsad Dagtekin

PhD Student at Population Ecology Research Group

University of Zurich, Department of Evolutionary Biology and Environmental Studies 

dilsad....@ieu.uzh.ch | www.popecol.org

Ecology and Evol. Biology Society of Turkey: www.ekoevo.org
ITU Alumni '16 & '18 | AFS Alumni '12

[Luke Wilde]

Hi Marc and Dilsad, 

I came upon this thread trying to answer the same question, and seeing the content of the thread was discouraged from trying to calculate WAIC (opting instead for Link and Barker's indicator variable approach with Bernoulli weights). However, I came upon Olivier Gimenez's code (https://gist.github.com/oliviergimenez/68ad17910a62635ff6a062f8ec34292f) and a discussion thread looking at WAIC for models with Jags where Martyn Plummer introduces JAGS functionality to estimate WAIC automatically (https://sourceforge.net/p/mcmc-jags/discussion/610036/thread/8211df61/#ea5c). I tried it with my CJS model (without adding in log_lik terms to the likelihood statement, mind you), and it produced WAIC estimates. 

Would you be inclined to distrust these, Marc?

Thanks so much,

Luke

> samples.mod <- jags.samples(model_jags$model, c("WAIC", "deviance"), type = "mean", + n.iter = 10000, + n.burnin = 3000, + n.thin = 10) > samples.mod$p_waic <- samples.mod$WAIC > samples.mod$waic <- samples.mod$p_waic + samples.mod$deviance > tmp <- sapply(samples.mod, sum) > waic.mod <- round(c(waic = tmp[["waic"]], p_waic = tmp[["p_waic"]]),1) > waic.mod waic p_waic 2144.5 26.6