Hierarchical Lincoln-Petersen model: flat discrete prior and start values

[erice...@gmail.com]

Dear all,

I need to estimate the effect of 4 experimental treatments on the size of 12 populations of aquatic snails (n = 3 populations/treatment). Snails were batch-marked and captured only twice (Lincoln-Petersen design). I started from a WinBUGS code provided as teaching material by Simon Bonner ( http://people.stat.sfu.ca/~cschwarz/Consulting/Trinity/Phase2/TrinityWorkshop/Workshop-material-Simon/). However, I am running into multiple problems (I am new to JAGS, population modelling and mark-recapture, so I might ask silly questions...):

1). The model requires sampling recapture probability p from an unknown number U of unmarked individuals (un[i] ~ dbin(p[i],U[i]) in the model below). p and U are the parameters to be estimated (I treat p as a random parameter so that the model is hierarchical). This requires providing a flat discrete prior for U, but I can't figure out how to do that. Any help with that would be greatly appreciated! As a work around, I have used a discrete informative prior produced by dbin, so that I was still able to pursue model fitting. Then, I met problem 2:

2). I can't find proper start values (from the previous posts, I see that this is a recurrent problem with JAGS). I am wondering whether this might come from the data (two populations, #1 and #7, went extinct during the course of the experiment). Would anyone have an opinion on that?

3). Finally, in case the two previous probelms would be solved, I would like to add a covariate (treatment) effect on population size (and ideally also on p), but this is another story (only commented in model code below).

Best wishes!

Eric



## 1) Model
cat(file = "Lincoln_petersen.jags", "
model{
    ## Likelihood function
      for(i in 1:pop){    #loop through populations
        ## Marked individuals
           mk[i] ~ dbin(p[i],n[i])
        ## Unmarked individuals
            un[i] ~ dbin(p[i],U[i])
        ## Population numbers
            N[i] <- n[i] + U[i]
        ##ANOVA for effects of experimental treatments on population size
        #N[i] ~ dpois(lambda)
        #log(lambda) <- alpha[treat[i]]   
      }
    ## Priors for population parameters
      for(i in 1:pop){
                ## Capture probabilities
        logit(p[i]) <- lproba[i]
            lproba[i] ~ dnorm(mu,tau)
            ## Number of unmarked fish
            U[i] ~ dbin(0.1,50)    #CANNOT FIND FLAT DISCRETE DISTRIBUTION!
      }
    ## Hyperpriors for capture probability
              mu ~ dnorm(0,.001)
              tau ~ dgamma(.01,.01)
    ## Priors for ANOVA parameters
      #for(i in 1:4){    #loop through experimental treatments
      #        alpha[i] ~ dnorm(0, 0.001)
      #}
}
")

## 2) Data list
bac <- c(1:12)                            #pop indices
treat <- c(rep(c("C","S","P","SP"), 3))                #experimental treatments
n <-c(0, 20, 19, 7, 1, 20, 0, 135, 30, 190, 14, 17)        #nb captured and marked
recapt <-  c(0, 21, 18, 8, 1, 23, 0, 150, 42, 251, 18, 28)    #nb recaptured
mk <- c(0, 20, 17, 7, 1, 20 , 0, 134, 30, 186, 14, 17)        #nb marked among recaptures
un <- (recapt - mk)                        #nb unmarked among recaptures

data <- list(mk = mk, n = n, un = un, pop = length(bac))

## 3) Initial values
inits <- function(){list(lproba = rep(1, 12), U = rep(1, 12))} #c(0,1,1,1,0,3,0,2,2,2,4,2)

## 4) Parameters monitored
parameters <- c("mu", "tau", "N")
## 5) MCMC settings
ni <- 2000
nt <- 5
nb <- 1000
nc <- 3
#model
library(jagsUI)
mod <- jags(data, inits, parameters, "Lincoln_petersen.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, n.adapt = 2000, parallel = TRUE)

[Daniel A. M. Villela]

Hi Eric,

If you did not capture any snails on any of the two occasions for these two populations (treatments?), there is not much that you can infer for the probability p.

Therefore, the design would not allow you to get an estimate of population size in these cases.  

I believe that you would have to remove these two populations and go ahead with the estimation, and report these two as inconclusive. Of course, other options would be to try multiple capture recaptures to avoid getting into no observations.

Finally, the idea to use the treatment as covariate sounds good.

Best,

Daniel

--
*** 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 (2015)
---
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 post to this group, send email to hmec...@googlegroups.com.
Visit this group at http://groups.google.com/group/hmecology.
To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/42118424-c2ce-405e-946c-8e68f4316ab6%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[John Clare]

Eric,

I believe a flat discrete distribution for capture recapture problems is typically implemented using data augmentation and a uniform (0, 1) prior for individual existence: z[i]~dbern(psi), where psi~Uniform(0, 1), and if i indexes possible individuals.    The idea is that the number of possible individuals is set to some extremely large number M, and the encounter history fed in as data is an M x numsamples matrix (for a Bernoulli encounter model), where n animals have observed encounter histories (010, 110, etc.), and M-n rows in the matrix are all zeroes: the question is how many of all of these unobserved individuals exist.  [Also note that Link (2013) strongly warns against using the discrete uniform, instead suggesting a scale prior approximated by psi~dbeta(1/M, 1); I'm not qualified to comment one way or the other].   Royle & Dorazio (2008) describe the data augmentation process in far more detail than I ever could, so that's likely a useful reference.

It seems to me that you could set up the encounter history for each site in a M rows by 2 column matrix, with the number of individuals initially marked having a 1 in the first column, and a 1 in the second column if recaptured; novel captures on the second occasion would have matrix values of 0, 1.  Technically, this is a multi-strata problem, so you might either choose to augment each site independently, or assign individuals from one large augmentation pool to each site.  Even though you can't really reconcile which individuals initially captured were recaptured, it probably doesn't matter for a straight forward constant detection type model.  I agree with Daniel that the 2 populations with 0 detections will not be particularly useful unless you treat p as constant across sites. 

Anyway, probably an easier way to go about this.

John

[erice...@gmail.com]

Thank you all for your feedbacks! I will try to write a model that works and be back to you asap. Regarding 0 populations, I find it weird to discard these observations because they are still valuable: extinction of a population is a quite strong effect of a treatment! I think that there should be way to account for this, but don't see how yet...

Cheers,

éric