JS model, Bayesian popualtion analysis

[diegoro...@gmail.com]

Hello all, 

I hope you are doing well. 

I am moving some population dynamics code from WinBUGS to NIMBLE. Right now I am working with the book from Kery and Schaub (Bayesian population analysis using WinBUGS). On section 10.3.1 the code from a Jolly Saber model as a Restricted Dynamic Occupacy Model. I manage it to make it work, but I want to model the parameter phi with a logit function and a respective individual covariate, which does not change over time.

What I find interesting is that the model do compile and run with out problem, but the betas values of my logit function are not getting sample. Therefore at the end of the MCMC chain, they keep the same initial value, with an sd=0.0. I already try setting different starting points and different prior distribution. 

Thank you very much for your time in advance. 

Please find attached the code, and the data that I am using. 

For the data, Columns 1:8 are the observed or not observed states

Col 9 is just a bird indicator

Col 10 is the Coeficient that I want to use to model the survival probability.

:)

[Rui Martins]

Hi,

your beta0 and beta1 need to be assigned with a prior distribution: beta1 ~ dist(....)

Cheers

Rui

--
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/782e0c38-4675-4589-ab45-c2795f315bf6n%40googlegroups.com.

[Daniel Turek]

 Rui is correct, for the beta's to be sampled by the MCMC, they need prior distributions.

What was confusing me is that in the attached file JS_NimbleCode.txt, the beta's do have priors:

  beta0~dflat()
  beta1~dflat()

But in the screen shot ("ScreenshotModel.png") we don't see the priors.

If you're still having trouble, then providing a fully reproducible example would be great.

Thanks,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/276114410.3368806.1631191185395%40mail.yahoo.com.

[Daniel Turek]

Thanks for the reproducible example.

It appears this came down to issues of (improper) model initialization.  The warnings you were seeing when you run the MCMC:

> #### Run the model ####
> SJ.samples <- runMCMC(SJ.cMCMC, niter = 1000, nburnin = 500,
+                       nchains = 3, summary=TRUE)
running chain 1...
warning: problem initializing stochastic node z[49, 3]: logProb is -Inf.
warning: problem initializing stochastic node z[53, 4]: logProb is -Inf.
warning: problem initializing stochastic node z[63, 4]: logProb is -Inf.
warning: problem initializing stochastic node z[118, 5]: logProb is -Inf.
warning: problem initializing stochastic node z[120, 5]: logProb is -Inf.
warning: problem initializing stochastic node z[168, 6]: logProb is -Inf.
warning: problem initializing stochastic node z[128, 7]: logProb is -Inf.
warning: problem initializing stochastic node z[149, 7]: logProb is -Inf.
warning: problem initializing stochastic node z[198, 7]: logProb is -Inf.
warning: problem initializing stochastic node z[226, 8]: logProb is -Inf.
warning: problem initializing stochastic node z[256, 8]: logProb is -Inf.

critically indicate that the model is starting at an invalid initial state - in this case, for many of the "z" nodes.  Since the MCMC cannot gain a foothold (does not have a valid starting state) at which to begin MCMC sampling, the sampling never really begins, and is unable to ever transition to a valid state - thus, as you observed, the betas are never updated.

A few changes to your inits list(s):

The "mean.phi" term in your inits lists was never used.

Providing initial values for "gamma" is also a good idea.

Providing an array of all 1's for "z", while not random and perhaps not optimal, is at least a valid starting initial state, from which the MCMC is able to begin sampling, and updating the beta's (and the z's).

You can use an inits list, perhaps as:

inits <- list(
    beta0 = 1,
    beta1 = 1,
    mean.p = 0.5,
    gamma = rep(0.5, constants$Occasions),
    z = array(1, c(constants$M, constants$Occasions))
)

Of course, you could also use runif(...) for the betas, mean.p, and gamma terms.  And also a more clever way of generating random (but *valid*) initial values for "z".

As a quick check, after you build the model you might check that the model is fully initialized.  Here are two ways to do that:

SJ.model$calculate()    ## make sure you get a real number (not NA or -Inf)

SJ.model$initializeInfo()    ## make sure it says that the model is fully initialized

Also, in your code, your inits was a list of 3 initial values lists, as:

inits <- list(  list(...), list(...), list(...))

If you want to do that, that's just fine.  Then use inits[[1]] for building the model, as you did.  However, if you want the 3 MCMC chains to begin using these 3 inits lists, then you also need to supply this to runMCMC:

J.samples <- runMCMC(SJ.cMCMC, niter = 500, nburnin = 200,
                      nchains = 3, summary=TRUE, inits = inits)

That addition will make each of your three chains begin using the 3 lists provided.

I hope this gets it working for you.  Keep in touch if you have other problems.

Cheers,

Daniel

On Thu, Sep 9, 2021 at 2:43 PM diegoro...@gmail.com <diegoro...@gmail.com> wrote:

Hello, 
Thank you for the quick response. 

I have assign dflat() as a prior and also dnorm(0.0,1. 0E-4) (in previous attempts), as show in the file JS_NimbleCode.txt
the .png file idea was to show where I am adding the beta0 and beta1.

Please, find attached all the code, that will be able to make a reproducibe example with the file "Data1.txt" of the first message.

The small amount of iteration is to make sure that the model is working. I have increase it, but the result still the same. 

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/00390c1a-0baa-4f12-b337-b0ae243d4c80n%40googlegroups.com.

[Jose Jimenez Garcia-Herrera]

Hi,

 

You need better starting values to fit the model. Try using *for z*:

 

z.init <- Data1

z.init[z.init==0] <- 1

 

and use:

 

inits<- list(…. z=z.init,…)

 

Good luck!

José

 

De: nimble...@googlegroups.com <nimble...@googlegroups.com> En nombre de diegoro...@gmail.com
Enviado el: jueves, 9 de septiembre de 2021 15:00
Para: nimble-users <nimble...@googlegroups.com>
Asunto: Re: JS model, Bayesian popualtion analysis

 

Hello, 
Thank you for the quick response. 

I have assign dflat() as a prior and also dnorm(0.0,1. 0E-4) (in previous attempts), as show in the file JS_NimbleCode.txt
the .png file idea was to show where I am adding the beta0 and beta1.

Please, find attached all the code, that will be able to make a reproducibe example with the file "Data1.txt" of the first message.

 

The small amount of iteration is to make sure that the model is working. I have increase it, but the result still the same. 

On Thursday, September 9, 2021 at 3:50:38 PM UTC+3 db...@williams.edu wrote:

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/00390c1a-0baa-4f12-b337-b0ae243d4c80n%40googlegroups.com.

[Diego Rondon]

Hello, 
Thank you for your email

It works now :) 

--

Have a great day!

Diego Rondon

[Rémi Patin]

Hi,
I had a quick look to see how the samplers were functionning with this model.
I allowed recording of the RW samplers properties (acceptance rate and scale history) by adding at the beginning of the script

nimbleOptions(buildInterfacesForCompiledNestedNimbleFunctions = TRUE)
nimbleOptions(MCMCsaveHistory = TRUE)

at the beginning of the script.
It appears that the RW samplers for beta0 and beta1 are working, but all proposals are refused (acceptance rate = 0). The sampler is trying to adapt and is decreasing the scale of the RW progressively, but all new values of beta0/beta1 are still refused.

(SJ.conf$getSamplers())[[1]] # to check that this sampler is beta0
(SJ.conf$getSamplers())[[2]] # to check that this sampler is beta1


SJ.samples <- runMCMC(SJ.cMCMC, niter = 100000, nburnin = 500,
                      nchains = 1, summary=TRUE)

# Acceptance history for beta1
> SJ.cMCMC$samplerFunctions[[2]]$acceptanceHistory
 [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[38] 0 0 0 0 0 0 0 0 0 0 0 0 0

# Scale history for beta1
> SJ.cMCMC$samplerFunctions[[2]]$scaleHistory
 [1] 1.000000e+00 2.342284e-01 6.955616e-02 2.435549e-02 9.632360e-03
 [6] 4.184906e-03 1.959725e-03 9.757500e-04 5.113471e-04 2.798741e-04
[11] 1.590176e-04 9.333845e-05 5.637687e-05 3.492686e-05 2.213414e-05
[16] 1.431594e-05 9.431663e-06 6.318920e-06 4.298888e-06 2.966065e-06
[21] 2.073176e-06 1.466555e-06 1.049038e-06 7.581905e-07 5.532984e-07
[26] 4.074402e-07 3.025849e-07 2.265103e-07 1.708381e-07 1.297638e-07
[31] 9.922632e-08 7.635721e-08 5.911277e-08 4.602460e-08 3.602924e-08
[36] 2.835077e-08 2.241895e-08 1.781189e-08 1.421539e-08 1.139404e-08
[41] 9.170403e-09 7.409987e-09 6.010296e-09 4.892817e-09 3.997109e-09
[46] 3.276421e-09 2.694424e-09 2.222767e-09 1.839227e-09 1.526317e-09

Additionaly it appears that there are warning like this when running the model:

warning: problem initializing stochastic node z[49, 3]: logProb is -Inf.
warning: problem initializing stochastic node z[53, 4]: logProb is -Inf.

For different z depending on the seed (I think, but maybe it is the different init values)
And all those z have also contant values, as well as gamma[7] and gamma[8]

Not sure what is happening but the problem is not limited to beta0 and beta1 but also on some z and gamma.
Hope this can help, although that is not a solution to your problem...
Cheers,
Rémi

On 09/09/2021 15:00, diegoro...@gmail.com wrote:

Hello, 
Thank you for the quick response. 

I have assign dflat() as a prior and also dnorm(0.0,1. 0E-4) (in previous attempts), as show in the file JS_NimbleCode.txt
the .png file idea was to show where I am adding the beta0 and beta1.

Please, find attached all the code, that will be able to make a reproducibe example with the file "Data1.txt" of the first message.

The small amount of iteration is to make sure that the model is working. I have increase it, but the result still the same. 

On Thursday, September 9, 2021 at 3:50:38 PM UTC+3 db...@williams.edu wrote:

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/00390c1a-0baa-4f12-b337-b0ae243d4c80n%40googlegroups.com.