Problem with unequal time intervals and marginalization
275 views
Skip to first unread message
Erica Mason
Sep 29, 2021, 4:02:10 PM9/29/21
to nimble-users
Hello all,
I have been stuck with marginalization using the dDHMMo function and would really appreciate some help. I have a working model with simulated data (without marginalization) that works with JAGS and Nimble.
Some notes about the model:
--extension of the multistate model in BPA book, section 9.5
--time is a random effect of survival
--I have unequal occasions, in which 1 occasion per year is a long non-survey period when angler resights and recoveries can occur
--real data has ~3,500 fish and 16 occasions per year (x12 years) (the reason I want to try Nimble for marginalization and weighting)
My first problem I run into using Nimble is the model comes up with -Inf in logProby whether I use unmarginalized or marginalized formulations. I think this is a result of the unequal intervals, in which I adjust parameters for the lengthier non-survey occasion by exponentiating (e.g., (R[t]^6 * equals(survey[t],0)), but I'm unclear how to resolve it....I understand the probabilities estimated are likely very small...but why is this handled okay in JAGS and not Nimble?
If I remove the distinction in state and obs probability calculations for survey vs nonsurvey periods, the -Inf problem goes away for the unmarginalized model run, but the new problem I run into with the marginalized model is the following," value in right-hand-side-only variable is NA or NaN, in variable: po" followed by repeated warnings for every single marginalization..."warning: logProb of data node y[1, 1:15]: logProb is NA or NaN.''' I am hoping this is due to some simple error I've made in specifying the likelihood for marginalization? During building, there is an error referring to the rows not summing to 1...missing value where TRUE/FALSE needed...but I think this has more to do with NAs somehow ending up in the prob matrices because I've checked that all rows do sum to 1.
Any ideas/suggestions on how to 1) resolve my unequal occasions problem, and 2) where I've gone wrong with the likelihood/inits for use with dDHMMo?
Code to simulate the data and run the model is attached. I've included a separate file with code using JAGS for comparison. The code file for use with Nimble contains sections with unmarginalized and marginalized models (with and without unequal occasions).
Olivier G. has been trying to help me with this, but we (including my advisor) haven't been able to move forward.
Thank you for your time!
Erica
Ben Goldstein
Sep 30, 2021, 6:24:53 PM9/30/21
to Erica Mason, nimble-users
Hi Erica,
Thanks for providing super clear and reproducible code, which makes giving help a lot easier.
The issue you've run into is a difficult one to identify. It relates to the fact that the observation probabilities and transition probabilities in the dynamic hidden Markov model are of different dimensions. Specifically, in nimbleEcology's dDHMMo, we ask for them in matrices of S x O x N and S x O x N-1, respectively. You've done all the indexing exactly right, which is why the model is compiling and not throwing errors until the MCMC begins.
The problem is that three variables--specifically po, p, and R--have a dimension of length N (n.occasions in the code) but values are only being filled for the first 1 through N-1 elements in a loop that goes from t = 1 to t = N-1. For example the 15th element of the vector p is NA and never defined. This is ok at model build time but when the MCMC starts to run the model notices that these never get filled in by sampling, since nothing is assigned to them--that's the reason for the warning "value in right-hand-side-only variable is NA". I think the reason this came up is that you're using the same loop to define a bunch of variables, some of which have a dimension N-1 and some of which have dimension N.
I solved this by adding the following lines of code just after the end of the for-loop over t. (I recommend you double-check these values as I wasn't thinking about the content of the model):
po[1,1,n.occasions] <- (p[n.occasions])
po[1,2,n.occasions] <- 0
po[1,3,n.occasions] <- ((1-p[n.occasions])*(1-R[n.occasions]))
po[1,4,n.occasions] <- (R[n.occasions]*(1-p[n.occasions]))
po[2,1,n.occasions] <- 0
po[2,2,n.occasions] <- 0
po[2,3,n.occasions] <- 1
po[2,4,n.occasions] <- 0
po[3,1,n.occasions] <- 0
po[3,2,n.occasions] <- 1
po[3,3,n.occasions] <- 0
po[3,4,n.occasions] <- 0
po[4,1,n.occasions] <- 0
po[4,2,n.occasions] <- 0
po[4,3,n.occasions] <- 1
po[4,4,n.occasions] <- 0
p[n.occasions] <- mean.p
R[n.occasions] <- mean.R
For future reference, to figure this out I had to use the longer workflow (using the nimbleModel, buildMCMC, mcmc$run functions). Once the mcmc had run (and, theoretically, sampled values were populated whether or not they were initialized) I was able to query the compiled model using model$p and see that p[15] was NA, which helped me find it in the code.
It's great to see nimbleEcology get some use. Hopefully it will provide some benefit to you. Let me know if you have any questions about the above, if I've gotten this wrong, or if there's anything else I can do to help!
Thanks,
Ben
Thanks for providing super clear and reproducible code, which makes giving help a lot easier.
The issue you've run into is a difficult one to identify. It relates to the fact that the observation probabilities and transition probabilities in the dynamic hidden Markov model are of different dimensions. Specifically, in nimbleEcology's dDHMMo, we ask for them in matrices of S x O x N and S x O x N-1, respectively. You've done all the indexing exactly right, which is why the model is compiling and not throwing errors until the MCMC begins.
The problem is that three variables--specifically po, p, and R--have a dimension of length N (n.occasions in the code) but values are only being filled for the first 1 through N-1 elements in a loop that goes from t = 1 to t = N-1. For example the 15th element of the vector p is NA and never defined. This is ok at model build time but when the MCMC starts to run the model notices that these never get filled in by sampling, since nothing is assigned to them--that's the reason for the warning "value in right-hand-side-only variable is NA". I think the reason this came up is that you're using the same loop to define a bunch of variables, some of which have a dimension N-1 and some of which have dimension N.
I solved this by adding the following lines of code just after the end of the for-loop over t. (I recommend you double-check these values as I wasn't thinking about the content of the model):
po[1,1,n.occasions] <- (p[n.occasions])
po[1,2,n.occasions] <- 0
po[1,3,n.occasions] <- ((1-p[n.occasions])*(1-R[n.occasions]))
po[1,4,n.occasions] <- (R[n.occasions]*(1-p[n.occasions]))
po[2,1,n.occasions] <- 0
po[2,2,n.occasions] <- 0
po[2,3,n.occasions] <- 1
po[2,4,n.occasions] <- 0
po[3,1,n.occasions] <- 0
po[3,2,n.occasions] <- 1
po[3,3,n.occasions] <- 0
po[3,4,n.occasions] <- 0
po[4,1,n.occasions] <- 0
po[4,2,n.occasions] <- 0
po[4,3,n.occasions] <- 1
po[4,4,n.occasions] <- 0
p[n.occasions] <- mean.p
R[n.occasions] <- mean.R
For future reference, to figure this out I had to use the longer workflow (using the nimbleModel, buildMCMC, mcmc$run functions). Once the mcmc had run (and, theoretically, sampled values were populated whether or not they were initialized) I was able to query the compiled model using model$p and see that p[15] was NA, which helped me find it in the code.
It's great to see nimbleEcology get some use. Hopefully it will provide some benefit to you. Let me know if you have any questions about the above, if I've gotten this wrong, or if there's anything else I can do to help!
Thanks,
Ben
--
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/de4fcd4f-08aa-41c4-8dee-b8c8abe8c28bn%40googlegroups.com.
Marc Kery
Oct 1, 2021, 4:19:08 AM10/1/21
to nimble...@googlegroups.com
Dear all,
I have been translating my little blue bugs book from WinBUGS to JAGS and Nimble and last night discovered that the simple linear regression example, which had worked fine the last time I tried it a couple of months ago, no longer does. Now, all Nimble produces
is a bunch of NaNs. Any ideas why this might be so ? Code example attached (and, btw, the problems persist also when the Bayesian p-value stuff is dropped).
Thanks a lot and best regards --- Marc
_____________________________________________________________
Marc Kéry * * * * *
+41 41 462 97 93
marc...@vogelwarte.ch
www.vogelwarte.ch
Swiss Ornithological Institute | Seerose 1 | CH-6204 Sempach | Switzerland
_____________________________________________________________
+41 41 462 97 93
marc...@vogelwarte.ch
www.vogelwarte.ch
Swiss Ornithological Institute | Seerose 1 | CH-6204 Sempach | Switzerland
_____________________________________________________________
Luca
Oct 1, 2021, 4:41:44 AM10/1/21
to marc...@vogelwarte.ch, nimble...@googlegroups.com
Hi Marc,
Not sure I can help, but a couple observations on your code:
– In "model7.5" some of the for-loops go from 1:16 instead of 1:n. I see that n<-16 at the beginning, but be careful in case you change it.
– Personally I'm never sure whether the third argument of dnorm() is tau or sd or var. *If* it is sd, then the original alpha of your dataset (a=40) has basically zero probability according to your prior.
– The priors for alpha and beta have a standard deviation of 0.001, but the init values you generate and from normals with SD=1. I wonder whether this may lead to roundoff errors and effectively zero probability depending on the init draws.
Cheers,
Luca
> To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/ZR0P278MB086948BFE36C79D674F87474EBAB9%40ZR0P278MB0869.CHEP278.PROD.OUTLOOK.COM <https://groups.google.com/d/msgid/nimble-users/ZR0P278MB086948BFE36C79D674F87474EBAB9%40ZR0P278MB0869.CHEP278.PROD.OUTLOOK.COM?utm_medium=email&utm_source=footer>.
Not sure I can help, but a couple observations on your code:
– In "model7.5" some of the for-loops go from 1:16 instead of 1:n. I see that n<-16 at the beginning, but be careful in case you change it.
– Personally I'm never sure whether the third argument of dnorm() is tau or sd or var. *If* it is sd, then the original alpha of your dataset (a=40) has basically zero probability according to your prior.
– The priors for alpha and beta have a standard deviation of 0.001, but the init values you generate and from normals with SD=1. I wonder whether this may lead to roundoff errors and effectively zero probability depending on the init draws.
Cheers,
Luca
> --
> 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 <mailto:nimble-users...@googlegroups.com>.
> You received this message because you are subscribed to the Google Groups "nimble-users" group.
> To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/ZR0P278MB086948BFE36C79D674F87474EBAB9%40ZR0P278MB0869.CHEP278.PROD.OUTLOOK.COM <https://groups.google.com/d/msgid/nimble-users/ZR0P278MB086948BFE36C79D674F87474EBAB9%40ZR0P278MB0869.CHEP278.PROD.OUTLOOK.COM?utm_medium=email&utm_source=footer>.
Marc Kery
Oct 1, 2021, 7:28:46 AM10/1/21
to Luca, nimble...@googlegroups.com
Dear Luca,
thanks for your help. However, I do think the second argument in the Normal distribution function in Nimble is the precision, so this can't be at the root of my problem. Plus: the thing ran fine a short while ago.
Best regards --- Marc
From: nimble...@googlegroups.com <nimble...@googlegroups.com> on behalf of Luca <pgl...@gmail.com>
Sent: Friday, October 1, 2021 10:41
To: Marc Kery <marc...@vogelwarte.ch>
Cc: nimble...@googlegroups.com <nimble...@googlegroups.com>
Subject: Re: problem with simple linear regression
Sent: Friday, October 1, 2021 10:41
To: Marc Kery <marc...@vogelwarte.ch>
Cc: nimble...@googlegroups.com <nimble...@googlegroups.com>
Subject: Re: problem with simple linear regression
Hi Marc,
Not sure I can help, but a couple observations on your code:
– In "model7.5" some of the for-loops go from 1:16 instead of 1:n. I see that n<-16 at the beginning, but be careful in case you change it.
– Personally I'm never sure whether the third argument of dnorm() is tau or sd or var. *If* it is sd, then the original alpha of your dataset (a=40) has basically zero probability according to your prior.
– The priors for alpha and beta have a standard deviation of 0.001, but the init values you generate and from normals with SD=1. I wonder whether this may lead to roundoff errors and effectively zero probability depending on the init draws.
Cheers,
Luca
On 2021-10-01 10:19, Marc Kery wrote:
> Dear all,
>
> I have been translating my little blue bugs book from WinBUGS to JAGS and Nimble and last night discovered that the simple linear regression example, which had worked fine the last time I tried it a couple of months ago, no longer does. Now, all Nimble produces is a bunch of NaNs. Any ideas why this might be so ? Code example attached (and, btw, the problems persist also when the Bayesian p-value stuff is dropped).
>
> Thanks a lot and best regards --- Marc
>
> _____________________________________________________________
>
> Marc Kéry * * * * *
> +41 41 462 97 93
> marc...@vogelwarte.ch
Not sure I can help, but a couple observations on your code:
– In "model7.5" some of the for-loops go from 1:16 instead of 1:n. I see that n<-16 at the beginning, but be careful in case you change it.
– Personally I'm never sure whether the third argument of dnorm() is tau or sd or var. *If* it is sd, then the original alpha of your dataset (a=40) has basically zero probability according to your prior.
– The priors for alpha and beta have a standard deviation of 0.001, but the init values you generate and from normals with SD=1. I wonder whether this may lead to roundoff errors and effectively zero probability depending on the init draws.
Cheers,
Luca
On 2021-10-01 10:19, Marc Kery wrote:
> Dear all,
>
> I have been translating my little blue bugs book from WinBUGS to JAGS and Nimble and last night discovered that the simple linear regression example, which had worked fine the last time I tried it a couple of months ago, no longer does. Now, all Nimble produces is a bunch of NaNs. Any ideas why this might be so ? Code example attached (and, btw, the problems persist also when the Bayesian p-value stuff is dropped).
>
> Thanks a lot and best regards --- Marc
>
> _____________________________________________________________
>
> Marc Kéry * * * * *
> +41 41 462 97 93
> marc...@vogelwarte.ch
> Swiss Ornithological Institute | Seerose 1 | CH-6204 Sempach | Switzerland
> _____________________________________________________________
>
> --
> 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 <mailto:nimble-users...@googlegroups.com>.
> To view this discussion on the web visit
https://che01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fnimble-users%2FZR0P278MB086948BFE36C79D674F87474EBAB9%2540ZR0P278MB0869.CHEP278.PROD.OUTLOOK.COM&data=04%7C01%7Cmarc.kery%40vogelwarte.ch%7Ce5441c6bff7144a9a8f008d984b74c32%7C251ae8d19a904569a35b9d2bc41188d2%7C1%7C0%7C637686745177123965%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=yYDTnZI2QL05DhAFZ%2BphcXEOq2BJ05Ej%2BhQf0L99%2Bak%3D&reserved=0
<https://che01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fnimble-users%2FZR0P278MB086948BFE36C79D674F87474EBAB9%2540ZR0P278MB0869.CHEP278.PROD.OUTLOOK.COM%3Futm_medium%3Demail%26utm_source%3Dfooter&data=04%7C01%7Cmarc.kery%40vogelwarte.ch%7Ce5441c6bff7144a9a8f008d984b74c32%7C251ae8d19a904569a35b9d2bc41188d2%7C1%7C0%7C637686745177123965%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=ohhrmNEV0Of0aOP3tokhpw6%2BKxVzUDqyo25Wir%2Fezzw%3D&reserved=0>.
--
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://che01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fnimble-users%2F2bc50630-aee6-de8f-b681-6a1156d866ed%2540gmail.com&data=04%7C01%7Cmarc.kery%40vogelwarte.ch%7Ce5441c6bff7144a9a8f008d984b74c32%7C251ae8d19a904569a35b9d2bc41188d2%7C1%7C0%7C637686745177123965%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=HGriuOcLiDZhRJrOU2ySWZ5uqDVIJrIB%2BToOSH79q0s%3D&reserved=0.
To view this discussion on the web visit https://che01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fnimble-users%2F2bc50630-aee6-de8f-b681-6a1156d866ed%2540gmail.com&data=04%7C01%7Cmarc.kery%40vogelwarte.ch%7Ce5441c6bff7144a9a8f008d984b74c32%7C251ae8d19a904569a35b9d2bc41188d2%7C1%7C0%7C637686745177123965%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=HGriuOcLiDZhRJrOU2ySWZ5uqDVIJrIB%2BToOSH79q0s%3D&reserved=0.
Ben Goldstein
Oct 1, 2021, 12:27:38 PM10/1/21
to Marc Kery, Luca, nimble...@googlegroups.com
Hi Marc,
I believe the problem can be fixed by providing both "x" and "y" to the model as data rather than as constants. I can't speak to why this was working earlier, but I was able to resolve the problem and get reasonable sampling by making that change.
Luca, dnorm() by default uses precision, but NIMBLE allows named arguments to distributions. So y ~ dnorm(0, 0.01) is equivalent to y ~ dnorm(0, sd = 10). If you're ever stumped on what a distribution's defaults and alternative parameterizations are you can check out page 2 of the NIMBLE cheat sheet!
Ben
I believe the problem can be fixed by providing both "x" and "y" to the model as data rather than as constants. I can't speak to why this was working earlier, but I was able to resolve the problem and get reasonable sampling by making that change.
Luca, dnorm() by default uses precision, but NIMBLE allows named arguments to distributions. So y ~ dnorm(0, 0.01) is equivalent to y ~ dnorm(0, sd = 10). If you're ever stumped on what a distribution's defaults and alternative parameterizations are you can check out page 2 of the NIMBLE cheat sheet!
Ben
Chris Paciorek
Oct 1, 2021, 8:40:41 PM10/1/21
to Ben Goldstein, Marc Kery, Luca, nimble...@googlegroups.com
Marc, thank you for noting this for us.
I've looked into it, and this is a bug that was introduced in version 0.11.0 (released in April 2021) when we improved the efficiency of our conjugate samplers to avoid doing unnecessary calculations in a number of cases. Unfortunately there is glitch in how we determine cases in which things can be simplified, such that if 'x' is specified in constants and if the first x[i] is exactly 1, our processing thinks that it can simplify the conjugate sampler calculation when in fact it can't. (A similar problem would arise if one had `mu[i] = alpha[i] + beta*x[i]` and 'alpha' is specified in constants with alpha[1] is exactly 0.)
As Ben noticed, a work-around is to specify 'x' in 'data' or 'inits', but not 'constants'. (One could also work around it by changing x[1] to be a small epsilon different from one or by reordering the x-y pairs such that the the first x is something other than exactly 1.)
Fortunately in your case it was obvious that something had gone wrong, but there are other situations in which the MCMC would not return obviously wrong results.
We're going to work on a quick fix for this and hopefully roll that out in the coming week. This will be version 0.12.1, which will be layered on top of our 0.12.0 release that is just appearing on CRAN but which we hadn't yet announced.
-Chris
--
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/CAEbSdx39D5UhfzgZFeehYSWX85eAKkZMR92BSFZsLKmLPoiSJQ%40mail.gmail.com.
Marc Kery
Oct 2, 2021, 5:34:02 AM10/2/21
to Chris Paciorek, Ben Goldstein, Luca, nimble...@googlegroups.com
Dear Chris and all others (including Pepe Jimenez, who wrote to me off-list),
thanks for helping so quickly !
Best regards --- Marc
From: Chris Paciorek <christophe...@gmail.com>
Sent: Saturday, October 2, 2021 02:40
To: Ben Goldstein <ben.go...@berkeley.edu>
Cc: Marc Kery <marc...@vogelwarte.ch>; Luca <pgl...@gmail.com>; nimble...@googlegroups.com <nimble...@googlegroups.com>
Sent: Saturday, October 2, 2021 02:40
To: Ben Goldstein <ben.go...@berkeley.edu>
Cc: Marc Kery <marc...@vogelwarte.ch>; Luca <pgl...@gmail.com>; nimble...@googlegroups.com <nimble...@googlegroups.com>
Subject: Re: problem with simple linear regression
Erica Mason
Oct 4, 2021, 7:24:45 PM10/4/21
to nimble-users
Hi Ben,
Thank you so much for your quick response! You are right about p and R requiring an additional dimension. Unfortunately, after some trial and error, I think there is yet another issue. Specifically, when we use the parameterization that assumes equal time intervals, some of the parameters become confounded and thus, the results using the additional code you provided are not adequately representative of the parm values used to simulate the data. Thus, I think the unequal time interval issue needs to be resolved first? Just as an aside, the values you inserted for the po matrix were appropriate when assuming equal time intervals. They would need to be adjusted in the other case.
To explore the unequal time interval issue, I went back to the unmarginalized version using parameterization for unequal time intervals. I looked into what might be going on with the -Inf vals referenced in my original post and using the longer workflow (thank you btw for making that more clear to me). I discovered that the -Inf vals are not coming from extremely small values due to exponentiation as I originally thought, but occurring where a 0 probability is drawn for the likelihood of y[i,t]. This is occurring right after the long nonsurvey period only at specific instances in which a fish is recaptured by a biologist (obs = 1). I think this is because the likelihood for y[i, t] is drawing from po[,,t-1] rather than po[,,t], which yields the impossible event. Is there a different way to specify the transition and observation probabilities during survey and nonsurvey periods than what I am currently using? Or, is this example not capable of translating well into the format necessary for dDHMMo? Also, not sure this even matters, but, in contrast to the dimensions of the transition and obs matrices used to simulate the data, the matrices in the model exclude individual due the dimensionality requirement of dDHMMo...
Thank you again for your time,
Erica
Ben Goldstein
Oct 6, 2021, 2:17:41 PM10/6/21
to Erica Mason, nimble-users
Hi Erica,
When I applied the above solution to the different-survey-lengths version--that is, I commented it back in and added the appropriate lines after the for-loop in the model code, I didn't see any -Inf probabilities in the unmarginalized code. I'll attach the file I'm running in case I've changed something important without realizing it. Let me know if I'm not understanding the problem you're describing; I have a feeling I don't have my head around it fully.
As an aside, I think po may be parameterized slightly differently in the marginalized and non-marginalized versions; I think in dDHMMo the probability of observing state B given true state A is O[A, B, t] not O[A, B, t - 1] as you've written it in the unmarginalized code. This is a matter of parameterization so there isn't a "right" way, but I wanted to point it out in case it causes confusion when developing the two cases simultaneously. You mention this difference in your email, and I think it can be resolved by just changing dcat(po[z[i, t], 1:4, t - 1]) to dcat(po[z[i, t], 1:4, t]), so long as those elements are defined appropriately above. That would at least bring it in line with what dDHMMo is doing internally.
You could conceivably have individual-specific transition and observation probability arrays. dDHMMo needs inputs of a certain dimension, but these could be slices from higher-dimension matrices. For example if ps is indexed by species with dimensions S x S x T x N, you could have
y[i, 1:n.occasions] ~ dDHMMo(init = init[i,1:4], # count data from f[i] + 1
probObs = po[1:4,1:4,1:n.occasions,i],
probTrans = ps[1:4,1:4,1:(n.occasions-1)],
len = n.occasions,
checkRowSums = 0) # if 1, the checker claims that rows don't sum to 1, but they do!
probObs = po[1:4,1:4,1:n.occasions,i],
probTrans = ps[1:4,1:4,1:(n.occasions-1)],
len = n.occasions,
checkRowSums = 0) # if 1, the checker claims that rows don't sum to 1, but they do!
which dDHMMo would handle just fine.
Ben
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/ec470e86-31db-400d-8f5e-01e6327a8c47n%40googlegroups.com.
Erica Mason
Oct 7, 2021, 5:07:14 PM10/7/21
to Ben Goldstein, nimble-users
Hi Ben,
Thanks again for such a quick turnaround. It is so great to have the model running now without yielding errors! Apologies for the confusion over the -Inf vals...I just needed to convince myself that come of the problems I was having weren't due to the exponents. Also, thank you for clarifying that additional dimensions are allowed beyond the stated requirement for dDHMMo.
I ran the model with Nimble using both the marginalized and unmarginalized (discrete version) code you updated for unequal time intervals. Interestingly, the discrete version gets a couple of the posteriors close, with the exception of sdS and mean.r, but the posteriors for the marginalized version are wonky (they don't match the simulated values). Yet, when I run the model using JAGS with the same simulated values (the discrete version as outlined in BPA Ch. 9.5 where y[i,t] ~ dcat(po[z,[i,t], i, t-1,], the posteriors from the JAGS run are close for all params). So it seems the model is still not translating well going from JAGS to Nimble? Here's the output for JAGS and Nimble (red vertical dashed lines on the density plots indicate the simulated values) :
JAGS (discrete, unequal time intervals):
2 chains, each with 7000 iterations (first 500 discarded), n.thin = 3
n.sims = 4332 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
ann.s[1] 0.531 0.058 0.417 0.493 0.532 0.570 0.646 1.001 4300
ann.s[2] 0.499 0.045 0.413 0.468 0.499 0.529 0.588 1.002 1500
mean.R 0.381 0.020 0.343 0.367 0.381 0.394 0.420 1.005 330
mean.RR 0.947 0.009 0.927 0.941 0.947 0.953 0.963 1.001 4300
mean.p 0.275 0.015 0.246 0.264 0.275 0.285 0.307 1.005 350
mean.r 0.611 0.046 0.520 0.579 0.612 0.643 0.698 1.001 4300
mean.s 0.910 0.025 0.852 0.896 0.913 0.928 0.949 1.007 360
s[1] 0.912 0.051 0.788 0.884 0.923 0.950 0.981 1.001 2100
s[2] 0.959 0.027 0.890 0.945 0.965 0.979 0.993 1.001 4300
s[3] 0.893 0.039 0.806 0.868 0.897 0.922 0.956 1.002 960
s[4] 0.952 0.030 0.878 0.937 0.958 0.974 0.991 1.009 300
s[5] 0.678 0.054 0.570 0.642 0.678 0.715 0.781 1.001 2300
s[6] 0.922 0.031 0.851 0.904 0.926 0.945 0.971 1.001 4300
s[7] 0.925 0.027 0.863 0.909 0.929 0.945 0.969 1.006 460
s[8] 0.909 0.029 0.844 0.891 0.913 0.930 0.958 1.001 4300
s[9] 0.943 0.034 0.862 0.923 0.949 0.969 0.989 1.003 860
s[10] 0.646 0.047 0.555 0.614 0.645 0.678 0.736 1.003 810
s[11] 0.940 0.026 0.879 0.924 0.944 0.959 0.980 1.006 360
s[12] 0.905 0.031 0.839 0.886 0.907 0.928 0.956 1.003 780
s[13] 0.905 0.031 0.836 0.885 0.909 0.927 0.954 1.001 4300
s[14] 0.919 0.072 0.725 0.894 0.940 0.968 0.990 1.056 450
sdS 0.542 0.806 0.000 0.054 0.248 0.687 2.739 1.001 4300
deviance 2422.954 21.151 2383.323 2408.372 2421.992 2436.450 2466.201 1.031 56
n.sims = 4332 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
ann.s[1] 0.531 0.058 0.417 0.493 0.532 0.570 0.646 1.001 4300
ann.s[2] 0.499 0.045 0.413 0.468 0.499 0.529 0.588 1.002 1500
mean.R 0.381 0.020 0.343 0.367 0.381 0.394 0.420 1.005 330
mean.RR 0.947 0.009 0.927 0.941 0.947 0.953 0.963 1.001 4300
mean.p 0.275 0.015 0.246 0.264 0.275 0.285 0.307 1.005 350
mean.r 0.611 0.046 0.520 0.579 0.612 0.643 0.698 1.001 4300
mean.s 0.910 0.025 0.852 0.896 0.913 0.928 0.949 1.007 360
s[1] 0.912 0.051 0.788 0.884 0.923 0.950 0.981 1.001 2100
s[2] 0.959 0.027 0.890 0.945 0.965 0.979 0.993 1.001 4300
s[3] 0.893 0.039 0.806 0.868 0.897 0.922 0.956 1.002 960
s[4] 0.952 0.030 0.878 0.937 0.958 0.974 0.991 1.009 300
s[5] 0.678 0.054 0.570 0.642 0.678 0.715 0.781 1.001 2300
s[6] 0.922 0.031 0.851 0.904 0.926 0.945 0.971 1.001 4300
s[7] 0.925 0.027 0.863 0.909 0.929 0.945 0.969 1.006 460
s[8] 0.909 0.029 0.844 0.891 0.913 0.930 0.958 1.001 4300
s[9] 0.943 0.034 0.862 0.923 0.949 0.969 0.989 1.003 860
s[10] 0.646 0.047 0.555 0.614 0.645 0.678 0.736 1.003 810
s[11] 0.940 0.026 0.879 0.924 0.944 0.959 0.980 1.006 360
s[12] 0.905 0.031 0.839 0.886 0.907 0.928 0.956 1.003 780
s[13] 0.905 0.031 0.836 0.885 0.909 0.927 0.954 1.001 4300
s[14] 0.919 0.072 0.725 0.894 0.940 0.968 0.990 1.056 450
sdS 0.542 0.806 0.000 0.054 0.248 0.687 2.739 1.001 4300
deviance 2422.954 21.151 2383.323 2408.372 2421.992 2436.450 2466.201 1.031 56
Gray circles/purple shading = JAGS estimated values
Black circles = values used to simulate data
Nimble (discrete, unequal time intervals):
user system elapsed
852.58 5.12 881.03
852.58 5.12 881.03
mean sd 2.5% 50% 97.5% Rhat n.eff
ann.s[1] 0.551 0.057 0.440 0.552 0.660 1.04 2123
ann.s[2] 0.510 0.043 0.427 0.510 0.599 1.01 2361
mean.R 0.364 0.018 0.328 0.364 0.399 1.04 2129
mean.RR 0.947 0.010 0.926 0.947 0.964 1.00 9975
mean.p 0.265 0.015 0.236 0.265 0.295 1.01 1968
mean.r 0.891 0.015 0.861 0.892 0.918 1.06 1924
mean.s 0.947 0.023 0.889 0.950 0.983 1.03 58
s[1] 0.941 0.046 0.823 0.953 0.995 1.02 1028
s[2] 0.983 0.018 0.932 0.989 1.000 1.01 1314
s[3] 0.921 0.033 0.843 0.925 0.972 1.01 1632
s[4] 0.991 0.009 0.966 0.993 1.000 1.00 1525
s[5] 0.618 0.051 0.515 0.620 0.715 1.02 2344
s[6] 0.940 0.027 0.876 0.944 0.981 1.00 1800
s[7] 0.944 0.023 0.892 0.946 0.981 1.00 1999
s[8] 0.935 0.023 0.882 0.937 0.971 1.00 1916
s[9] 0.993 0.007 0.974 0.995 1.000 1.00 1635
s[10] 0.587 0.043 0.501 0.586 0.669 1.03 1692
s[11] 0.966 0.018 0.921 0.969 0.991 1.00 1651
s[12] 0.937 0.023 0.884 0.939 0.973 1.01 1525
s[13] 0.938 0.021 0.891 0.940 0.972 1.01 1228
s[14] 0.994 0.006 0.978 0.996 1.000 1.00 1571
sdS 1.754 0.480 1.056 1.670 2.930 1.00 608
ann.s[1] 0.551 0.057 0.440 0.552 0.660 1.04 2123
ann.s[2] 0.510 0.043 0.427 0.510 0.599 1.01 2361
mean.R 0.364 0.018 0.328 0.364 0.399 1.04 2129
mean.RR 0.947 0.010 0.926 0.947 0.964 1.00 9975
mean.p 0.265 0.015 0.236 0.265 0.295 1.01 1968
mean.r 0.891 0.015 0.861 0.892 0.918 1.06 1924
mean.s 0.947 0.023 0.889 0.950 0.983 1.03 58
s[1] 0.941 0.046 0.823 0.953 0.995 1.02 1028
s[2] 0.983 0.018 0.932 0.989 1.000 1.01 1314
s[3] 0.921 0.033 0.843 0.925 0.972 1.01 1632
s[4] 0.991 0.009 0.966 0.993 1.000 1.00 1525
s[5] 0.618 0.051 0.515 0.620 0.715 1.02 2344
s[6] 0.940 0.027 0.876 0.944 0.981 1.00 1800
s[7] 0.944 0.023 0.892 0.946 0.981 1.00 1999
s[8] 0.935 0.023 0.882 0.937 0.971 1.00 1916
s[9] 0.993 0.007 0.974 0.995 1.000 1.00 1635
s[10] 0.587 0.043 0.501 0.586 0.669 1.03 1692
s[11] 0.966 0.018 0.921 0.969 0.991 1.00 1651
s[12] 0.937 0.023 0.884 0.939 0.973 1.01 1525
s[13] 0.938 0.021 0.891 0.940 0.972 1.01 1228
s[14] 0.994 0.006 0.978 0.996 1.000 1.00 1571
sdS 1.754 0.480 1.056 1.670 2.930 1.00 608
Gray circles/purple shading = Nimble estimated values
Black circles = values used to simulate data
Nimble (marginalized, unequal time intervals):
user system elapsed
1011.45 2.65 1036.67
1011.45 2.65 1036.67
mean sd 2.5% 50% 97.5% Rhat n.eff
ann.s[1] 0.955 0.020 0.915 0.956 0.994 1.00 2595
ann.s[2] 0.879 0.026 0.826 0.879 0.927 1.00 2860
mean.R 0.096 0.006 0.085 0.096 0.107 1.00 2231
mean.RR 0.947 0.009 0.927 0.947 0.964 1.00 9506
mean.p 0.173 0.007 0.161 0.173 0.187 1.02 2331
mean.r 0.989 0.009 0.966 0.992 1.000 1.01 1157
mean.s 0.984 0.007 0.968 0.985 0.995 1.09 90
s[1] 0.995 0.004 0.985 0.995 0.999 1.00 1925
s[2] 0.997 0.003 0.989 0.998 1.000 1.00 1429
s[3] 0.983 0.007 0.967 0.984 0.994 1.00 2277
s[4] 0.997 0.003 0.989 0.997 1.000 1.01 1472
s[5] 0.931 0.015 0.899 0.932 0.958 1.00 2239
s[6] 0.985 0.007 0.968 0.986 0.995 1.00 2394
s[7] 0.982 0.008 0.964 0.983 0.993 1.00 2553
s[8] 0.974 0.009 0.953 0.975 0.989 1.00 2301
s[9] 0.996 0.003 0.988 0.997 1.000 1.00 1609
s[10] 0.886 0.020 0.844 0.887 0.922 1.00 2259
s[11] 0.986 0.007 0.969 0.987 0.997 1.00 2329
s[12] 0.970 0.011 0.945 0.971 0.987 1.00 2021
s[13] 0.965 0.012 0.938 0.966 0.985 1.00 2342
s[14] 0.996 0.004 0.985 0.997 1.000 1.00 1759
sdS 1.501 0.413 0.878 1.442 2.425 1.02 608
ann.s[1] 0.955 0.020 0.915 0.956 0.994 1.00 2595
ann.s[2] 0.879 0.026 0.826 0.879 0.927 1.00 2860
mean.R 0.096 0.006 0.085 0.096 0.107 1.00 2231
mean.RR 0.947 0.009 0.927 0.947 0.964 1.00 9506
mean.p 0.173 0.007 0.161 0.173 0.187 1.02 2331
mean.r 0.989 0.009 0.966 0.992 1.000 1.01 1157
mean.s 0.984 0.007 0.968 0.985 0.995 1.09 90
s[1] 0.995 0.004 0.985 0.995 0.999 1.00 1925
s[2] 0.997 0.003 0.989 0.998 1.000 1.00 1429
s[3] 0.983 0.007 0.967 0.984 0.994 1.00 2277
s[4] 0.997 0.003 0.989 0.997 1.000 1.01 1472
s[5] 0.931 0.015 0.899 0.932 0.958 1.00 2239
s[6] 0.985 0.007 0.968 0.986 0.995 1.00 2394
s[7] 0.982 0.008 0.964 0.983 0.993 1.00 2553
s[8] 0.974 0.009 0.953 0.975 0.989 1.00 2301
s[9] 0.996 0.003 0.988 0.997 1.000 1.00 1609
s[10] 0.886 0.020 0.844 0.887 0.922 1.00 2259
s[11] 0.986 0.007 0.969 0.987 0.997 1.00 2329
s[12] 0.970 0.011 0.945 0.971 0.987 1.00 2021
s[13] 0.965 0.012 0.938 0.966 0.985 1.00 2342
s[14] 0.996 0.004 0.985 0.997 1.000 1.00 1759
sdS 1.501 0.413 0.878 1.442 2.425 1.02 608
Gray circles/purple shading = JAGS estimated values
Black circles = values used to simulate data
It seems like we are really close....maybe just missing something simple? As you can see in the discrete Nimble version, mean.r is being estimated really high for some reason, as well as sdS (I'm using a gamma prior on it, but the posterior doesn't reflect that distribution). Then, things go really haywire when we try to use dDHMMo with the marginalized version. I should also say that I tried running the marginalized version without a random effect of time on survival and that didn't seem to improve things. I also tried indexing on individual with both the discrete and marginalized versions and that didn't change things either.
I hope this helps to make the problem more clear...ahh, it's difficult to do this over email! :) It would be really great to get the marginalization going as it should theoretically save a lot of computational time when it comes to running the actual data? As the code is currently written, the model run time is longest for marginalization and least for the JAGS run (though to be fair, I had to really bump up the iterations in Nimble to reach convergence).
The version of code I used to generate the results above is attached. The JAGS code I sent earlier remains unchanged.
Erica
Erica Jarvis Mason (she/her/hers)
PhD Candidate, Semmens Lab
Marine Biology Research Division
Scripps Institution of Oceanography, UC San Diego
I respectfully acknowledge that we live and work on traditional, ancestral & unceded land of the Kumeyaay nation. Whose land are you on?
_____________________________________________________
Ben Goldstein
Oct 18, 2021, 6:50:17 PM10/18/21
to Erica Mason, nimble-users
Hi Erica,
I'm finally getting around to taking a look at this. Thanks for your patience. A couple things I'm noticing:
It looks like in the unmarginalized version, you're only looking at the elements in each row of y starting after first capture--specifically, only y[i, f[i]:n] contributes to the likelihood. In the marginalized version, all elements of y[i, 1:n] are contributing to the likelihood. (I think if elements before f[i] are NA, they're ignored, but to be safe I wouldn't count on that.) The substantive difference here is that the marginalized version isn't forcing the state y[i, f[i]] to represent a capture but is treating it the same as any other capture, as the result of a transition from unknown states and also an imperfect detection process.
I think you can achieve equivalence with the unmarginalized version using the following:
y[i, (f[i]+1):n.occasions] ~ dDHMMo(init = po[1, 1:4, f[i]], # use the transition matrix from captured to first unknown state
probObs = po[1:4, 1:4, 1 (f[i]+1):n.occasions],
probTrans = ps[1:4, 1:4, (f[i]+1):(n.occasions-1)],
len = n.occasions - f[i],
probObs = po[1:4, 1:4, 1 (f[i]+1):n.occasions],
probTrans = ps[1:4, 1:4, (f[i]+1):(n.occasions-1)],
len = n.occasions - f[i],
checkRowSums = 0) # if 1, the checker claims that rows don't sum to 1, but they do!
This will cause errors if
f[i]+1 == n.occasions because NIMBLE will insist on the input being a matrix, not a vector. The solution is to handle these cases in a separate for-loop using unmarginalized code (which doesn't require a latent state). But this is a bit of a headache, so if possible could you just prevent this from happening as a condition during simulation (or drop those data)? If not I can provide more info on this.
Unrelated to all this, I see your comment about checkRowSums giving you a problem. There's a phenomenon where R and (compiled) C++ handle rounding differently that can lead to this check passing for uncompiled models but failing for compiled models with the same values. We should have caught this in the latest version of nimbleEcology but if not that's useful to know--are you on the latest version (0.4.0)?
So the above should explain at least in part why the marginalized version looks different from the unmarginalized version.
I also see a potential reason for the discrepancy between JAGS and NIMBLE-unmarginalized. It looks like in JAGS the array is index with t-1, similarly to ps. This is the reason you're getting away with only having 1:(n.occasions-1) time elements in po. In NIMBLE, probObs[i, j, t] is the probability that an individual in state i is observed in state j at time t, not at time t+1 as you've written it in JAGS. If you want to use the same parameterization of po in NIMBLE's unmarginalized version it should work fine with t-1. I think this was caused by some confusing advice I gave in a previous email. If you want to keep your code defining the elements of po[] as written you can stick with the t-1 parameterization in the unmarginalized version. In the marginalized version, you can use this same thing if you substitute
probObs = po[1:4, 1:4, 1
(f[i]):(n.occasions-1)]
for
probObs = po[1:4, 1:4, 1
(f[i]+1):n.occasion]
Hopefully this more comprehensive look is helpful. Let me know if there's further confusion or if you'd like me to elaborate on anything here. There may be more to do when these two points are addressed, and I'm happy to stay looped in.
Ben
Erica Mason
Oct 19, 2021, 12:54:59 PM10/19/21
to Ben Goldstein, nimble-users
Thanks Ben! I appreciate your time looking into this for me. The indexing issues really help to clarify what is going on...I will take some time to digest all of this and happy to keep you looped in!
Reply all
Reply to author
Forward
0 new messages