warning: problem initializing stochastic node among other warnings

[Jacob Trowbridge]

Hi all,

I'm trying to make an integrated population model and I'm running into some warnings and errors. Copied below is the code. The model works fine when I don't include and "H" variables (Hf, Hurl, Had). Essentially, H is harvest animals that need to be subtracted from the estimated population size "N" to calculate "N" the next year. The model also works well if I include "H" but make it deterministic rather than stochastic at the Process model spot. Unfortunately, I can't share any data due to a data sharing agreement. I can give some insight to the structures. Nhat, var.Nhat is a vector 25 values long, 1 for each year. "f" is a 1000+ vector that is for 7 years of data. H (harvest) is a matrix 3 rows 25 columns. No harvest occurred for the first 17 years, so there are 0's. Then a known amount was removed for each ageclass (rows) from 18 and on. The errors I got on this code block with the data, initials, and constants below were many but the main ones were: 

1) warning: value of deterministic node ann.growth.rate[7]: value is NA or NaN even after trying to calculate.

2) warning: value of deterministic node Ntot[22]: value is NA or NaN even after trying to calculate.

3) warning: logProb of data node Nhat[21]: logProb is NA or NaN.

4) warning: value of stochastic node N[3, 24]: value is NA or NaN even after trying to simulate.
5) warning: problem initializing stochastic node N[3, 24]: logProb is NA or NaN.

This model gave more errors than I normally got for tinkering on the model. Any help would be greatly appreciated! Thanks!

ipm3 <- nimbleCode({
  # Priors and linear models for population dynamics
  sig.ystar ~ dunif(0, 10000)
  sig2.ystar <- pow(sig.ystar, 2)
  tau.ystar <- pow(sig.ystar, -2)
 
  # Prior for Survival
  mean.sj ~ dunif(0, 1)
  mean.sa ~ dunif(0, 1)

  # Priors for Fecundity
  for(t in 1:nyear){
    mean.f[1,t] ~ dunif(0, 1)
    mean.f[2,t] ~ dunif(0, 4)
    mean.f[3,t] ~ dunif(0, 4)
  }
 
  # Likelihood
  # Model for the initial population size: uniform priors
  for(a in 1:ageclass){
    N[a,1] ~ dunif(100, 2000) # All 3 age classes
  }
 
  # Process Model
  for(t in 2:nyear){
    N[1,t] ~ dpois(round(((N[1,t-1] - Hf[t-1]) * mean.f[1,t-1] / 2 * mean.sj) +
                          (N[2,t-1] * mean.f[2,t-1] / 2 * mean.sa) +
                          (N[3,t-1] * mean.f[3,t-1] / 2 * mean.sa)))
    N[2,t] ~ dbinom(mean.sj, round(N[1,t-1] - Hf[t-1]))
    N[3,t] ~ dbinom(mean.sa, round(N[2,t-1]- Hyrl[t-1] + N[3,t-1]- Had[t-1]))
  }
 
  # Observation process for the abundance estimates with their SEs
  for(t in 1:nyear){
    Nhat[t] ~ dnorm(ystar[t], tau.se[t])
    tau.se[t] <- 1/var.Nhat[t]          # Assumed known and given by variance of Nhat
    ystar[t] ~ dnorm(N[1,t]+N[2,t]+N[3,t], tau.ystar)
  }
 
  # Recruitment
  for(i in 1:nind){
    J[i] ~ dpois(mean.f[age[i], year[i]])
  }
 
  # Derived parameters
  # Annual population growth rate
  for (t in 1:(nyear-1)){
    ann.growth.rate[t] <- (N[1,t+1] + N[2,t+1] + N[3,t+1]) / (N[1,t] + N[2,t] + N[3,t])
  }
  # Total population size
  for (t in 1:nyear){
    Ntot[t] <- N[1,t] + N[2,t] + N[3,t]
  }
})
nimdata = list(J = f$Fetus.Count, Nhat = D$abundace, var.Nhat = D$var,
               Hf = c(rep(0, times = 17), hcopy$n[which(hcopy$ageclass == 1)]),
               Hyrl = c(rep(0, times = 17), hcopy$n[which(hcopy$ageclass == 2)]),
               Had = c(rep(0, times = 17), hcopy$n[which(hcopy$ageclass == 3)]))
params = c("mean.f", "sig2.ystar", "sig.ystar", "ystar", "mean.sj", "mean.sa",
           "N", "ann.growth.rate", "Ntot")
nimconsts = list(nyear = length(D$abundace),
                 nind = length(f$Fetus.Count),
                 year = f$Year-1998,
                 age = f$ageclass,
                 ageclass = length(unique(f$ageclass)))
niminits = list(mean.sj = runif(1, 0, 1), mean.sa = runif(1, 0, 1),
                mean.f = matrix(rep(c(runif(1,0,1),runif(1,0,4)), times = c(25,50)),
                                nrow = 3, ncol = 25, byrow = T))

[Chris Paciorek]

Hi Jacob,

A couple reactions based on a fairly quick look at this:

1) I wonder if a lot of the warnings you list are because you are getting negative values (or zero) for the mean of a Poisson or the sample size of a binomial. Note that you are calculating these based on subtraction, so if you're not careful about initialization, you could get negatives or zero. Furthermore, even if initialized so that everything is valid, in MCMC sampling, one could propose invalid values for the Poisson mean or binomial sample size and, I suspect one would get NA or NaN as log probabilities. 

2) Please see here for some ideas of how to troubleshoot initialization issues. It's framed from the perspective of initializing for running an MCMC but should be relevant for simply understanding what is going on when one first creates a model.

-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/e6b54a8b-ef85-49e7-ad6d-59e12c6cb957n%40googlegroups.com.

[Perry de Valpine]

Hi Jacob,

In case you are still stuck, to follow up a bit more, here is a strategy for tracking down the source of NAs and/or NaNs. You can inspect the model's variables and calculate any subset of the model you want or order to narrow down the problem. For example, you could do

model$ann.growth.rate[7] # look at the first problem node (where `model` is what is returned from a call to nimbleModel, giving the errors you quoted above). 

Since this is NA or NaN, presumably something about its inputs is invalid. You can look at those nodes

model$N[1:3, 8]

model$N[1:3, 7]

etc.; you can keep tracing up the graph (what depends on what) to look at model$Hf, mean.f, etc to track down the problem.

If you want to calculate just some bit of the model, you can do

model$calculate("N[1:3, 8]") # for example

If you want to modify a value and see if that fixes things, you can do say

model$Hf <- <some value you want to try>

 and then calculate nodes that depend on Hf, for example:

model$calculate( model$getDependencies( "Hf" ) ) # or write out the specific nodes by hand as argument to model$calculate

HTH.

Perry

 

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CA%2B7Ts_qnbeRHdUmiL-1rSvV5h3NyUzNZDziuf7oRo2EHgySN%3Dw%40mail.gmail.com.