Recursive Function for writing Stan code for Single Release Cormac-Jolly-Seber model of any length.

[Dalton Hance]

I'm working on a project comparing estimates of survival and detection for mark-recapture from Stan to two maximum likelihood methods. Specifically, I'd like to simulate survival and detection typical for mark juvenile Chinook salmon in the Snake River that migrate past eight hydropower projects which serve as detection points. These fish are listed under the endangered species act and mark-recapture models have become an integral part of management and recovery (at the cost of several million dollars a year). My simulation numbers are derived from this report.  For the simulation fish are only marked once at the beginning of the migration. They do not lose marks and we're going to assume independence of each fish, independence of survival and detection at any one project, and independence among projects.

I'm attempting to adapt the code on page 133 of the Stan manual to this situation. But with 128 potential detection histories, that is a lot of code to write! I spent a couple of tedious hours copy and pasting, but I'm still not done yet. I'm also thinking that there's got to be a better way. There is a natural recursive structure to the likelihood for each capture history. So I'm thinking it should be possible to write a function in R to automatically write the code for me for any number of recapture events. I've gotten a start on some code, but I'm getting a little stuck on how to write out the model piece of the code. I've attached my R script that writes out everything through the transformed parameters stage and I've got a place holder for a double for loop that I think can be used to write each increment_log_prob step. I'm hoping for a little help from this community on how to proceed, or maybe for somebody to tell me that I'm thinking about this the wrong way and that there's a better way.

Here is an excerpt from the stan code that show the recursive structure to modelling each detection history. I work the capture history from left to right. We always start with a one. You can see the base 2 structure. I'm sure there's got to be a clever way of stepping through a nested for loop to reproduce this for any number of detection histories, but it's eluded me so far. Any ideas?

model {

  

  increment_log_prob(history[1] * ( #10000000

                     log(chi[1])));

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

  increment_log_prob(history[2] * ( #11000000

                     log(phi[1]) + log(p[2]) + 

                     log(chi[2])));

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

  increment_log_prob(history[3] * ( #10100000

                     log(phi[1]) + log1m(p[2]) +

                     log(phi[2]) + log(p[3]) + 

                     log(chi[3])));                     

  increment_log_prob(history[4] * ( #11100000

                     log(phi[1]) + log(p[2]) +

                     log(phi[2]) + log(p[3]) + 

                     log(chi[3])));

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

  increment_log_prob(history[5] * ( #10010000

                     log(phi[1]) + log1m(p[2]) +

                     log(phi[2]) + log1m(p[3]) +

                     log(phi[3]) + log(p[4]) +

                     log(chi[4])));

  increment_log_prob(history[6] * ( #11010000

                     log(phi[1]) + log(p[2]) +

                     log(phi[2]) + log1m(p[3]) +

                     log(phi[3]) + log(p[4]) +

                     log(chi[4])));

                     

  increment_log_prob(history[7] * ( #10110000

                     log(phi[1]) + log1m(p[2]) +

                     log(phi[2]) + log(p[3]) +

                     log(phi[3]) + log(p[4]) +

                     log(chi[4])));

  increment_log_prob(history[8] * ( #11110000

                     log(phi[1]) + log(p[2]) +

                     log(phi[2]) + log(p[3]) +

                     log(phi[3]) + log(p[4]) +

                     log(chi[4])));

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

 

[Christopher Peterson]

I've been working on some code to handle imputing missing categorical covariates, which also involves marginalizing multiple binary variables in combination.  The project is still very much a work in progress, but the marginalization part of it appears to be working.  Not sure how relevant it will be to your problem, but the link is here:

https://github.com/Christopher-Peterson/stan-imputeDiscreteMissingData/blob/dev/discreteMissingDataFunctions.stan

Essentially, it takes a vector of (0,1) responsibility parameters and computes every combination of them; then calculates the likelihood of each combination, and uses them as a mixture. 

[Dalton Hance]

I will attempt to answer my own question. But I think it's still up in the air as to whether my approach is the correct one, or if there is a more elegant way to implement what I'm trying to do without resorting to goofy paste and cat calls.

I wrote a recursive function that translates a capture history vector to Stan code chunks in character vector. This should work for any length of capture history for the single release study design. It shouldn't be too hard to modify to match the classic CJS situation.

The recursive function and function which writes out a Stan model for the correct length of capture history is attached. 

[Bob Carpenter]

The CJS model with arbitrary numbers of capture times
is here:

https://github.com/stan-dev/example-models/blob/master/misc/ecology/mark-recapture/cjs-K.stan

It relies on precomputing first and last capture time and
then looping through the middle.

An even better way to do this would be to write up a bunch
of these pieces using functions. I had thought I or someone
else had done that at some point, but I can't find it.

Was that the kind of model you were looking for? If not, there
are also a bunch of CJS models in Hiroki Itô's translation of
some of the examples from Kéry and Schaub's book:

https://github.com/stan-dev/example-models/blob/master/example-models/BPA

If not, if you can specify the functions you need in math
notation without trying to generate Stan code, we can probably
figure out how to do it either with a loop or recursively.

For instance, if you're not working at the individual level, you can
code up the histories as input data and key off that in the code.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <cjs_write_stan_functions.R>

[Dalton Hance]

Thanks Bob.

I was attempting to use the aggregated capture histories rather than working at the individual level. My motivation for doing so is purely for computational efficiency though. When I was playing around with the model from the Stan manual with only three capture events after marking (8 possible detection histories), it ran blazingly fast; about a minute or so for 50000 iterations with 17500 individuals. I was hoping to see not much more of a penalty for bumping up to seven capture events after marking (128 capture histories), but I'm looking at about 30 minutes or so for 50000 iterations with 17500 individuals. Maybe I have unrealistic hopes, but I'd like to cut that down. 

I'll play around with your code and Hiroki Itô's version to see how they compare. Maybe there is a better way to work out the marginal likelihood in Stan as well.

I'm not sure I know what you mean by code up the histories as input data and key off that in the code. I think that is what I'm doing. The capture histories translate nicely as binary representation of integers. I relied on that structure to write recursive function in R (moving from right to left) to come up with the increment_log_probability chunks. But maybe there is a way to do a similar thing in Stan without relying on my goofy R function. I'm still learning the Stan language though, so I'm sure I'm making some amateur errors.

As a side note: if the group is interested I did code up a multi-level version of the model on page 133 (multiple releases with 4 capture events each) that I'll clean up and post to my github for sharing.

[Ben Goodrich]

On Thursday, February 18, 2016 at 12:13:13 PM UTC-5, Dalton Hance wrote:

I was hoping to see not much more of a penalty for bumping up to seven capture events after marking (128 capture histories), but I'm looking at about 30 minutes or so for 50000 iterations with 17500 individuals. Maybe I have unrealistic hopes, but I'd like to cut that down. 

The real issue is not that your hopes are unrealistic but that your number of iterations is too large. For the vast majority of well-written Stan models, the effective sample sample size is going to be at least 10% of the nominal post-warmup sample size. And you do not need an effective sample size of 2500 to estimate a mean sufficiently well. Just using the default values of 2000 iterations (of which half are discarded as warmup) with 4 chains should get you an effective sample size of a few hundred. And if you had a few hundred independent draws from the posterior, you would be fine in terms of estimating the posterior means. There are some exceptions to all this, but you shouldn't start by assuming your research project is one of the exceptions.

Ben

[Dalton Hance]

Thanks Ben. That is very good point that I had not considered. But my goal here is specifically to compare the frequentist properties of the Bayesian version of CJS to available MLE methods. I'm interested in the coverage and width of credible intervals in addition to the posterior mean. Actually I was thinking of using the posterior mode and HPD intervals, but if there are pitfalls there I'd be happy to be educated. I'd like to directly compare the credible intervals to confidence intervals coming out of Jeff Laake's marked package for example.

I'll play around with number of iterations and see how that affects the stability of credible intervals. 

[Ben Goodrich]

On Thursday, February 18, 2016 at 12:45:38 PM UTC-5, Dalton Hance wrote:

Thanks Ben. That is very good point that I had not considered. But my goal here is specifically to compare the frequentist properties of the Bayesian version of CJS to available MLE methods. I'm interested in the coverage and width of credible intervals in addition to the posterior mean. Actually I was thinking of using the posterior mode and HPD intervals, but if there are pitfalls there I'd be happy to be educated. I'd like to directly compare the credible intervals to confidence intervals coming out of Jeff Laake's marked package for example.

I'll play around with number of iterations and see how that affects the stability of credible intervals. 

Here you have the additional consideration that the endpoints of the 95% credible intervals are not estimated very precisely for the default settings. It would be better to compare the 68% intervals or something.

Ben

[Bob Carpenter]

You can do the aggregated histories the same way. You just
need to store the aggregation counts and multiply by that
inside increment_log_prob.

And yes, please post example models if you have them.

- Bob