trouble with nimble ecology dDHMMo distribution

[Kirsty Laurenson]

Hi all, 

I have a multistate CMR model that forms part of a larger IPM that I am trying to marginalise using the distributions from NimbleEcology as the run time for my IPM is approaching the time limit for jobs on my University’s HPC. However, I’m running into some issues which I haven’t been able to solve after reading through the introduction to NimbleEcology page, or the relevant conversation threads in this group. 

When I run ‘nimbleModel’ on my model, I get the message telling me my model is not fully initialised, and when I investigate this using model$initializeInfo(), I get: Missing values (NAs) or non-finite values were found in model variables: omega. Omega is my observation matrix. When I look closer at model$omega, the values for omega at the last occasion are all NA. Additionally, after compiling the model and MCMC and trying to run the model, I get the ‘logProb of data node is NA’ warnings, which I suspect is linked to the issues with omega. 

I have attached a script with the model code to run on simulated data. The script includes the non-marginalised model, which runs fine, and a version of the model where I have tried to use the dDHMMo distribution, without success. I get the same warnings when I run the dDHMMo model on the simulated data and my real data. I have a feeling these problems are arising because I haven’t quite coded the model correctly using the dHMMo distribution, but after reading the relevant conversation threads on this forum and the intro to nimble ecology, I still haven’t been able to wrap my head around it! Any help or suggestions would be greatly appreciated!

Best, 

Kirsty

[Thomas Riecke]

Hi Kirsty,

I've struggled a bit with Nimble's marginalized capture history functions in the past, but have some working code I could share, and I'm sure experienced nimble users will have better advice. Before diving down that rabbit hole though, a couple of quick potentially useful thoughts:

If you don't have individual covariates (all covariates are cohort-level or based on t), marginalizing capture histories by individual is generally a very slow way to do things:

It looks like you’re using Schaub-type code. The IPMbook also has examples for multistate m-arrays which typically run much faster (orders of magnitude). You could run this on a good laptop instead of an HPC. Note that there is a minor error in the code as written in the IPM book, it’s corrected in the erratum. Mike Meredith wrote some really nice code here with Michael, Richard Chandler has built some similar code, all based on Williams et al. (2002). You could also marginalize across all individuals with identical capture histories (following Yackulic et al. 2020, EcoApps), but I think that's generally slower than the multinomial approach as well.

You could (and should?) substantially reduce your number of states as well, regardless of whether you have individual covariates or not, by simply making the transition to breeder conditional on time since release.

for (t in f[i]:T){

   psi[pre-breeder, i, t, breeder] <- alpha[t - f[i] + 1] * phi[whatever]

}

This would also greatly reduce runtime. Something like this would still give you different alpha's for each age, rather than having them move forward through a bunch of temporary classes. Then you would have fledglings, non-breeders, and breeders as your states instead of NB2, NB3, .... , NB6+. The absorbing state can be marginalized out. This will substantially reduce the model's complexity. There are many different ways to index these parameters. Dan used a similar trick in one of our papers, there are lots of other good examples:

https://onlinelibrary.wiley.com/doi/pdf/10.1111/gcb.16482

Hope that's helpful and not too off-topic, just a note that rethinking how you write the model may be more helpful here than fixing the code.

Thomas

--
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 visit https://groups.google.com/d/msgid/nimble-users/0ab480ae-a395-49d2-a2de-cc03611e3322n%40googlegroups.com.

[Kirsty Laurenson]

Hi Thomas, 

Thanks for sharing your thoughts! I've had a quick look at your suggestions and those links to papers and code you've shared and I think they will be very helpful!

Best, 

Kirsty

[Perry de Valpine]

Hi Kirsty,

Thanks for posting the question, and Thomas thanks for the useful contribution.

Kirsty, the reproducible example is great, and also the work you show to pin down where the issue is. What I see is that omega is indeed only declared in the model over "t in 1:(K-1)", so that omega[,,30] (where K = 30) was not declared, and that would explain why its values are NAs. I think since the other components being declared within the same for loop only extend over "t in 1:(K-1)", you will want to put omega in a separate for loop over "t in 1:K". (For the technically curious, additional for loops in model code do not add meaningful overhead, but additional declarations do, so it is better to make a separate for loop than to add new, separate, declarations for the final value of t = K). I tried that change, and I also moved p.NB and p.B into the new for loop that covers "t in 1:K", and I also extended their initial values to K instead of K-1. This seems to have resolved the error you were seeing. 

The next issue revealed is that you have some individuals with first[i] == 28 so that the range of data is y[i, 28:29] and the corresponding range of transitions is 29:29 (first[i] + 1):(K-1). You'll need to change those cases to use dHMMo instead of dDHMMo since there is only a single transition so it can't be time-dependent. I see your "first" vector is nicely sorted already so this was easy to do. When I made this change, the model built without error. I stopped there and hope the MCMC will work.

The attached shows my changes.

Perry

To view this discussion visit https://groups.google.com/d/msgid/nimble-users/6f03528c-f47b-4b93-8c9c-35c23d95af6an%40googlegroups.com.

[Kirsty Laurenson]

Hi Perry, 

Thanks very much for your help - I suspected the issue had something to do with the different number of occasions required for the observation matrix and the transition matrix but just couldn't quite figure it out! 

Best, 

Kirsty