Convergence issues with default sampler for time stratified CMR model

[guillaum...@gmail.com]

Good afternoon everyone,

 First of, apologies, I deleted my earlier post as I had identified some indexing issues and didn’t want people losing time looking at a model that was incorrect.

My model ‘should’ now be correct.

I provide some context to the model at the end of this message but basically, it’s a variation of a time stratified CMR model. The reason why I’m soliciting feedback here is that I’m very confused by the behavior of the default Nimble sampler.

One set of parameters I'm trying to estimate seems to be generating the issue (at least that’s what I'm observing).

One component of the model looks at the probability that a fish release on day i will be available to recapture on day j, this is done through a discretized distribution function of lognormal distribution.

In this formulation I have parameters lambda[t,i] which are the average log(number of days) until fish cohort released on day i will be available to recapture (t indicates the year).

I can get a version of my model to run correctly i.e., good mixing of chains, low Rhat, good effective size when I set all lambdas to be independent:

lambda(t,i) ~ Normal(mean=log(1), sd = 0.8).

Everything seems fine and values estimated make sense biologically and the chains of all parameters seem fine too. Few examples below:

Fig 1. Examples of total abundances MCMC chains

Fig2. Examples of daily proportions/number of fish starting their migration on a given day (yaxis, first and last day at the top and bottom of the graph respectively) note the larger uncertainties when the trap is not running which are expected.

 

I will note that it seems like this lambda parameter is very sensitive, if I change the prior to something less informative such as Normal(mean=0,tau=0.01) I get some issues with chains not mixing and looking like they are getting stuck on certain values.

 

Similarly, when I modify the model to incorporate a hierarchical structure on lambda such as lambda(t,i) ~ Normal(mu_lambda[t], sigma_lambda[t]) (i.e. in a given year each cohort of fish released on day i, should have similar travel times), I get poor mixing. And chains of a range of parameters, notably the total annual abundances NmTot[t] being stuck on very small range of values.

Fig3 example of total abundances for model with hierarchical structure on lambda, note that two of the chaina (green and red) are moving within a very narrow range of values

 

When digging into it, the output of the models failing to converge seems to find a place in the space of parameters that supports low values of lambda and very high values of catchabilities at the trap i.e., close to 1 when we know from non time stratified models and recaptured/marked fish that the catchability is at best in the low tens (~0.1).

 
Fig 4 examples of chains and poterior distribution of daily catchabilities at the trap

It is a bit frustrating not knowing what is going wrong and It would be great to get any feedback about this. I have attached Rdata with model, constants, Data and inits and script I am using to compile and run.

 

Note: when building model the following warning is printed

Defining model

  [Warning] Multiple definitions for the same node.

            Did you forget indexing with 'k' on the left-hand side of

            `theta[i, j, t] <- phi((log(i - (i - j) + 0.5) - lambda[t, I_lambda_nimble[i, j, t]])/phi_[t]) - nimStep(i - (i - j)) * phi((log(i - (i - j) - 0.5) - lambda[t, I_lambda_nimble[i, j, t]])/phi_[t])`?

  [Warning] Multiple definitions for the same node.

            Did you forget indexing with 'k' on the left-hand side of

            `data_RR[i, j, t] ~ dbin(size = data_RT[t, I_lambda_nimble[i, j, t]], prob = lifted_q_oBt_comma_i_minus_first_day_RT_oBt_cB_plus_1_cB_times_theta_oBi_comma_j_comma_t_cB_L42[t, i, j], lower_ = 0, upper_ = data_RT[t, I_lambda_nimble[i, j, t]])`?

 

This is due to the fact that sometimes there is no likelihood for some combinations of index i,j,t and I need to skip them. From previous discussions I understand that this should be fine (alternative to conditional 'for' loop (google.com)) but maybe this causes some issues…

 

Thanks

 

Guillaume

General context of the model

I am running a time stratified (daily) mark recapture model to estimate abundance of a fish run.

Each a year a trap is operating (almost) daily, fish captured are marked and released a couple of km upstream and we record any recaptures.

Fish migration occur during a short period (a few weeks).

Ideally, we would like to run the trap have a couple of days with 0 captures a ramp up to a maximum and a ramp dow followed by again a few 0.

The trick is that the trap is not always operational during the whole migration and as a consequence some years the trap is installed and we get a bunch of fish from day 1.

To compensate for the potential missing part of the run , I am trying to estimate the probability that a fish will start its migration on a given day based on environmental covariates and using it to obtain the proportion of fish migrating on a given day and therefore available to be captured by a trap when operational.

[guillaum...@gmail.com]

Hello again,
just bringing this up a bit to see if anyone has encountered a similar issue and has suggestion to fix this.
Thanks

[Perry de Valpine]

Dear Guillaume,

Thanks for refreshing your question. I am not quite sure if it is a nimble question (how the system works) or a statistical question (you are not getting the estimates you think you should). Our main purpose here is the former, and your model is somewhat complicated for one to dig into the latter unless there is someone very familiar with the model context already. What strikes me from the warning and looking at the relevant code is that either you should remove the k loop or, possibly, you might be hoping your code is doing something that it isn't. I don't see any use of k at all within your k loop, so I am not able to follow what the intention of writing the k loop is. Can you clarify? As written, nimble may be clever enough to ignore the k loop and use only one declaration of each theta[i,j,t] and each data_RR[i,j,t]. However, there is really no possible meaning in the model language of embedding declarations within a loop whose index is never used. Hence I am wondering if this might be a source of confusion that we could help sort out.

Best wishes,

Perry

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/f2cec3e3-74fd-41bd-98ce-9e57eebb156an%40googlegroups.com.

[guillaum...@gmail.com]

Good afternoon Perry,

Thank you for your input, I fully appreciate that this is a complicated modelling object and that it is not an easy thing to dive into. 

re the k-loop, this is a trick to define likelihoods only when they actually exist (in this particular case there are some days in a season when the tagging is interrupted and therefore no fish from that daily cohort can be recaptured).
the loop index I_RR is a Boolean, when ==  0 the loop is ignored  when == 1 the code is read / likelihood defined  .

This trick was available in Bugs and Jags and Daniel Turek confirmed in an earlier conversation that this would work too in nimble

but to your point I am wondering if this loop could potentially impact the parameter space exploration when there is a hierarchical structure on a parameter present in this loop

thanks again

G

[Perry de Valpine]

Thanks for the explanation. That makes sense. I think the idea is that "for(i in 1:0)" should result in nothing. It doesn't look like nimble really supports that. E.g.

m <- nimbleModel(nimbleCode({x ~ dnorm(0,1); for(i in 1:0) y ~ dnorm(0,1)}))

m$getNodeNames()

m$getVarInfo()

#  A scalar node y was actually created

(And on the other hand, this produces an error:

m <- nimbleModel(nimbleCode({x ~ dnorm(0,1); for(i in 1:0) y[i] ~ dnorm(0,1)}))

)

You could check if you have been getting nodes you didn't really want if for example you try:

dailyCMR_dir_model$calculate("theta[1, 2, 3]") 

dailyCMR_dir_model$calculate("data_RR[1, 2, 3]") 

where  you replace "1, 2, 3" with a set of i, j, t values that should not have either of these nodes.

Supporting the trick from BUGS and JAGS is something we should think about. It would provide a path through this need. Meanwhile, here is another way to do it:

If you create in R three vectors, i_values, j_values, and t_values, such that together they contain every triplet of i,j,t values for which you need a model declaration, then you can write (as a single for-loop, not a nested set of for-loops)

for(k in 1:length_of_values) {

 theta[i_values[k], j_values[k], t_values[k] ] <- [etc...]

 data_RR[i_values[k], j_values[k], t_values[k] ] ~ dbin(<etc....>)

}

where length_of_values is length(i_values) (same for the other two vectors). Creating these vectors may seem a little bit silly in terms of redundancy, but if they are passed as constants they shouldn't be a burden in the model bookkeeping.

Will that solution work for you?

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/074288ff-f52c-4eeb-86f5-c024f753397en%40googlegroups.com.

[Guillaume Dauphin]

Hi Perry, 

Thanks for all of that.
I'm going to explore these options and hopefully get a fix!

Thank you!

G