Hierarchical Bayesian Distance Sampling based on Sollmann et al 2015

Question

I'm working to adapt code from a paper by Sollman et al (2015) that uses distance sampling to correct for imperfect bird detection using species-specific and community parameter estimates for multiple species (https://doi.org/10.1111/2041-210X.12518). This code is also similar to that of Kery and Royle's HM textbook, though Kery and Royle only estimate single species.

I'm encountering errors with a variable 'fct[s,j,k]', which is a matrix with dimensions equal to the number of species (s), number of sites (j), and number of distance intervals (k), and is used in the code to connect detection probabilities with the observation model. The error I am getting states that node fct[1,3,2] has "Invalid parent values".

The full code (including data prep above the model) is as follows. I'm using a reduced dataset of six sites and five species, modeled by one detection covariate and one process covariate for the purposes of testing:

###########################################################
# Load Packages --------------------------------------------------------------
library(tidyverse)
library(jagsUI)
###########################################################
# Set WD --------------------------------------------------------------
setwd("/Users/colinsweeney/Documents/Documents/PhD/Chapter1/2017_DistanceModel_Test/Models")


Y<- matrix(as.integer(c(0,0,3,0,0,0, 0,0,0,0,3,1, 0,0,0,1,0,0,  1,0,0,0,0,0, 0,1,0,0,0,0)), nrow=6, ncol=5)
colnames(Y)<- c(1,2,3,4,5)

Nin=Y+1 

nsp=5 # number species
species=c(4, 5, 1, 3, 2, 2, 1, 2, 1, 2) # species index
nG=2 # number of distance bins, previously: nD=nD, 
pix=c(0.25, 0.75) # % area covered by each category
midpt=c(62.5, 187.5) # previously: midpt=midpt,
nsites=6
nVAR=1
nOBSVAR=1
VAR = matrix(c( 1.40683409, -0.33717511, -0.91851151,2.07537095,-0.14824078, 0.05522696,-0.91851151,-0.14824078,-0.91851151,-0.14824078), nrow=10, ncol=1)
colnames(VAR)<- c("elevation")
OBSVAR=matrix(c(1.42884613,1.26074658,0.08404977,0.08404977,-1.26074658, 0.75644795, 0.08404977,-1.26074658,0.08404977,-1.26074658),nrow=10, ncol=1)
colnames(OBSVAR)<-c("cloudcover")
dclass=c(2, 2, 1, 1, 1, 2, 1, 1, 1, 1)
y=Y
nind=10
site=c(1, 2, 3, 4, 5, 6, 3, 5, 3, 5)

# Bundle and summarize data set
str(win.df_NEON <- list(spec=nsp,
                        species=species, 
                        nG=nG, 
                        pix=pix, 
                        midpt=midpt,
                        nsites=nsites, 
                        nVAR=nVAR,
                        nOBSVAR=nOBSVAR,
                        VAR = VAR,
                        OBSVAR=OBSVAR,
                        dclass=dclass,
                        y=y,
                        nind=nind,
                        site=site) )

##################################################################################################
##################### JAGS model code for community DS model: NEON Adaptated ####################
cat(" # commented out just for the purposes of editing code (delete before running)
model{ 
  #################################### Priors #######################
  mu_s~dnorm(0,0.01) # part of asig: species specific covariate
  sig_s~dunif(0,100) # part of tau_s, which is part of asig: species specific covariate
  tau_s<-1/(sig_s*sig_s) # part of asig: species specific covariate
  
  mu_a~dnorm(0,0.01) # part of alpha: species specific covariate
  tau_a~dgamma(0.1,0.1)# part of sig_a, which is part of alpha: species specific covariate
  sig_a<-1/sqrt(tau_a) # part of alpha: species specific covariate
  
  r.N~dunif(0,100) # community hyperparamter: used in rho.N
  
  ########################################################
  # community hyperpriors
  ########################################################
  
  for (k in 1:nVAR){
    mu_b[k]~dnorm(0,0.01)
    sig_b[k]<-1/sqrt(tau_b[k]) # This doesn't appear to be used anywhere??
    tau_b[k]~dgamma(0.1,0.1)
  }
  
  for (k in 1:nOBSVAR){ # NEW: Trying to rework bsig (Observation covariates): treat it like beta
    mu_bsig[k]~dnorm(0,0.01)
    tau_bsig[k]~dgamma(0.1,0.1)
    sig_bsig[k]<-1/sqrt(tau_bsig[k])
  }
  
  ########################################################
  ### draw species-specific parameters from hyperdistribution
  ########################################################
  for (s in 1:spec){
    for (k in 1:nVAR){
      beta[s,k]~dnorm(mu_b[k],tau_b[k])
    }
    
    ########   ############   ############
    for (k in 1:nOBSVAR){ # NEW Trying to rework bsig (Observation covariates): treat it like beta
         bsig[s,k]~dnorm(mu_bsig[k],tau_bsig[k]) # community hyperparameter for detection model 
       }
    
    asig[s] ~dnorm(mu_s, tau_s) # detection model species specific intercept
    
    alpha[s]~dnorm(mu_a,tau_a)  # ecological model species specific intercept
  }
  
#########################################################################################  
### Main HM - > parallel structure to Kery and Royle 
#########################################################################################
  

#########################################################################################
    for (s in 1:spec){ # New 
    
    for (j in 1:nsites){ # number of sites
    
      ########### linear model for detection (Detection model) #######
      sigma[s,j]<- exp(asig[s] + inprod(bsig[s,],OBSVAR[j,]) )
      
      
      ### construct detection probabilities by distance class w/ half-normal model
      for(k in 1:nG){ # nG: number of distance intervauls (nG is db-1 in length)
        log(p[s,j,k])<- -midpt[k]*midpt[k]/(2*(sigma[s,j]*sigma[s,j])) # half-normal detection
        
        fc[s,j,k]<- p[s,j,k]*pix[k] # pix: proportion of area of each intervaul
        fct[s,j,k]<- fc[s,j,k]/sum(fc[s,j,1:nG])  # swap pcap for it's value AND fsc -> fst
        
      }
      
      pcap[s,j]<-sum(fc[s,j,1:nG])                   # overall detection probability
      
      ############ Part 2 of HM   ################
      y[j,s]~ dbin(pcap[s,j],N[j,s]) 
      
      ############ Part 3 of HM   ################
      N[j,s]~dpois(lambda.star[j,s]) 
      
        ########### linear model for abundance (Ecological model) #######
        log(lambda[j,s])<- alpha[s] + inprod(beta[s,],VAR[j,])

        rho.N[j,s]~dgamma(r.N,r.N) 
        
        lambda.star[j,s]<-lambda[j,s]*rho.N[j,s]
    }
  }
  
  for(i in 1:nind){
    
    ############ Part 1 of HM   ################
    dclass[i] ~ dcat(fct[species[i],site[i],1:nG]) 
    
    # distance class of individual distributed by catagorical distribution of binned probabilites
    # nG: number of distance bins
  }
  ################################################################################################
  ### monitor total abundance
  for (i in 1:spec){
    Nspec[i]<-sum(N[1:nsites,i])
  }
} 
    ",fill=TRUE, file="MultispeciesDist_2017NEON_ver1.txt")

################################## MCMC Settings #########################
# Inits --------------------------------------------------------------
# Initial values to start run at

inits<-function(){list(N=Nin,
                       # Shape paramters
                       mu_a=runif(1,0,1), # alpha
                       tau_a=runif(1,0,1), # alpha
                       
                       mu_b=runif(nVAR,0,1), # beta
                       tau_b=runif(nVAR,0,1), # beta
                       
                       mu_bsig=runif(nOBSVAR,0,1), # NEW: bsig (OBSVAR)
                       tau_bsig=runif(nOBSVAR,0,1), # NEW: bsig (OBSVAR)
                       
                       asig=runif(nsp,0,1), # OBSVAR
                       mu_s = runif(1, 0, 1), # asig (OBSVAR) ~ like alpha
                       sig_s=runif(1) )} # asig (OBSVAR) ~ like alpha

# Params --------------------------------------------------------------
# Params to save: list of parameters to save (variables of interest)

params<-c('asig', # species specific parameter, used in detection model
          'mu_s', # species specific shape parameter for asig (detection model )
          'sig_s', # species specific shape parameter for asig (detection model )
          'mu_a', # species specific shape parameter for alpha (ecological model)
          'alpha', # species specific covariate, ecological model (intercept)
          'sig_a', # species specific shape parameter for alpha (ecological model)
          'mu_b', # species specific shape parameter for beta (ecological model)-VAR
          'sig_b', # species specific: beta parameter(ecological model)-VAR
          'beta', # species specific covariate, ecological model
          'Nspec', # abundance for each species 
          'r.N', # community hyper parameter
          'mu_bsig', # NEW
          'sig_bsig', # NEW
          'bsig' # MODIFIED # new format to match beta's format
)

# MCMC chain settings--------------------------------------------------------------
# ni: number of iterations in MCMC chain
# nb: burn in period
# nt: thinning
# nc: number of chains

ni <- 52000   ;   nb <- 2000   ;   nt <- 2   ;   nc <- 3

# TEST ONLY
ni <- 300   ;   nb <- 200   ;   nt <- 2   ;   nc <- 3

###########################################################
# Make sure wd is correct--------------------------------------------------------------
setwd("/Users/colinsweeney/Documents/Documents/PhD/Chapter1/2017_DistanceModel_Test/Models")
# Make sure the txt file is exported to the correct folder

###########################################################
# Run MCMC (Gibbs Sampler) in JAGS --------------------------------------------------------------
##### Run JAGS and summarize posteriors.
# jagsUI 
MultiSp_FullTest_1 <- jags(win.df_NEON, 
                           inits, 
                           params, 
                           "MultispeciesDist_2017NEON_ver1.txt", 
                           n.thin=nt,
                           n.chains=nc, 
                           n.burnin=nb, 
                           n.iter=ni)

I've tried to modify initial values and priors but that hasn't appeared to make a difference with the "invalid parent values" error. Any thoughts are appreciated.

score 0 · Answer 1 · answered Aug 03 '22 at 07:15

This almost always means you are not specifying priors/joint likelihoods that "respect" the support of your prior distributions/joint likelihood. This can happen, for example, if you (implicitly) try to sample from distributions with negative variances, Poisson distributions with negative lambda, or divide by zero. My understanding is that jags does not allow improper priors as well.

While I don't have an exact solution for you because your example isn't a minimal example, here's what I recommend:

Start with a simpler version of the model and add in complexity one component at a time to ensure you are getting sensible estimates, things converge, etc. This will also help you identify where you are passing unsupported value of parameters or allowing the sampler to move into those spaces.
Move from the top of the hierarchical structure to the unit observations to ensure that everything is supported and appropriately constrained. For example, sig_s ~ dunif(0, 100) imposes a uniform prior on the closed interval [0, 100]. You then divide 1 by this value when defining tau_s as follows: tau_s <- 1 / (sig_s * sig_s). tau_s is then the precision (1/var) of a normal distribution (asig). If sig_s = 0 then you're asking jags to (a) set tau_s to 1/0 and (b) give the normal distribution asig a variance of 1/(1/0). Infinite variance is not a sensible prior and will produce a bunch of issues if jags doesn't just throw an error like this.
Take a look at the jags forums, which will be much better-suited for this question. This error appears in over 1,331 posts on those forums, so I am quite sure you will find a solution there in conjunction with the other parts of this.

Thanks advice! I was able to figure it out in the end. It actually had to do with the midpt values, resulting in the value fct being divided by zero. — colinpsweeney, Aug 08 '22 at 15:22

Hierarchical Bayesian Distance Sampling based on Sollmann et al 2015

1 Answers1