NAs produced during simulation from nimbleCode

[Marie Ozanne]

Hello! I am attempting to simulate from a SIRS model in NIMBLE with a seasonal force of infection (beta[t] in the second set of code below), and all of the random distributions generate NAs (e.g.,  Warning messages: 1: In rbinom(1, size = model$I[getNodeFunctionIndexedInfo(INDEXEDNODEINFO_, ... : NAs produced ) . 

The first set of code below does not include a seasonal component and works. The second set is the simplest example I could come up with that throws NA warnings. Do you have any suggestions on how to approach this problem? Thank you in advance for any help you can provide.

Code that works:

SIRS_code <-  nimbleCode({
 
  S[1] <- N - I0 - R0
  I[1] <- I0
  R[1] <- R0
 
  probIR <- 1 - exp(-gamma)
  probRS <- 1 - exp(-phi)
 
  ### loop over time
  for(t in 1:tau) {
   
    probSI[t] <- 1 - exp(- beta1 * I[t] / N)
   
    Istar[t] ~ dbin(probSI[t], S[t])
    Rstar[t] ~ dbin(probIR, I[t])
    Sstar[t] ~ dbin(probRS, R[t])
   
    # update S, I, R
    S[t + 1] <- S[t] - Istar[t] + Sstar[t]
    I[t + 1] <- I[t] + Istar[t] - Rstar[t]
    R[t + 1] <- R[t] + Rstar[t] - Sstar[t]
   
    # data not perfectly observed
    nu[t] <- logit(Istar[t]/N+0.0001)
    logit(Y[t]) ~ dnorm(nu[t], 0.1)
   
  }
 
  # priors
  beta1 ~ dgamma(0.1, 0.1)
  gamma ~ dgamma(0.1, 0.1)
  phi ~ dgamma(0.1, 0.1)
})

n_weeks <- 52
n_seasons <- 2
constantsList <- list(N = 10000,
                      I0 = 5,
                      R0 = 0,
                      tau = n_weeks*n_seasons)


sirsModel <- nimbleModel(SIRS_code,
                         constants = constantsList)

# exclude data from parent nodes
dataNodes <- c('Istar', 'Rstar', 'Sstar', 'logit_Y')
dataNodes <- sirsModel$expandNodeNames(dataNodes, returnScalarComponents = TRUE)
parentNodes <- sirsModel$getParents(dataNodes, stochOnly = TRUE)
parentNodes <- parentNodes[-which(parentNodes %in% dataNodes)]
parentNodes <- sirsModel$expandNodeNames(parentNodes, returnScalarComponents = TRUE)
nodesToSim <- sirsModel$getDependencies(parentNodes, self = FALSE, downstream = T)

initsList <- list(beta1 = 0.6,
                  gamma = 0.2,
                  phi = 0.1)
sirsModel$setInits(initsList)

set.seed(1)
sirsModel$simulate(nodesToSim, includeData = TRUE)

Minimal example that throws NA warnings:

SIRS_code <-  nimbleCode({
 
  S[1] <- N - I0 - R0
  I[1] <- I0
  R[1] <- R0
 
  probIR <- 1 - exp(-gamma)
  probRS <- 1 - exp(-phi)
 
  ### loop over time
  for(t in 1:tau) {
   
    probSI[t] <- 1 - exp(- beta[t] * I[t] / N)
    ## need to add seasonal beta to this next
    beta[t] <- 0.5*(1+0.9*sin((2*pi/52)*(w[t]-52/4)))
   
    Istar[t] ~ dbin(probSI[t], S[t])
    Rstar[t] ~ dbin(probIR, I[t])
    Sstar[t] ~ dbin(probRS, R[t])
   
    # update S, I, R
    S[t + 1] <- S[t] - Istar[t] + Sstar[t]
    I[t + 1] <- I[t] + Istar[t] - Rstar[t]
    R[t + 1] <- R[t] + Rstar[t] - Sstar[t]
   
    # data not perfectly observed
    nu[t] <- logit(Istar[t]/N+0.0001)
    logit(Y[t]) ~ dnorm(nu[t], 0.1)
   
  }
 
  # priors
  gamma ~ dgamma(0.1, 0.1)
  phi ~ dgamma(0.1, 0.1)
})

n_weeks <- 52
n_seasons <- 2
constantsList <- list(N = 10000,
                      I0 = 5,
                      R0 = 0,
                      tau = n_weeks*n_seasons,
                      w = rep(1:n_weeks, n_seasons))


sirsModel <- nimbleModel(SIRS_code,
                         constants = constantsList)

# exclude data from parent nodes
dataNodes <- c('Istar', 'Rstar', 'Sstar', 'logit_Y')
dataNodes <- sirsModel$expandNodeNames(dataNodes, returnScalarComponents = TRUE)
parentNodes <- sirsModel$getParents(dataNodes, stochOnly = TRUE)
parentNodes <- parentNodes[-which(parentNodes %in% dataNodes)]
parentNodes <- sirsModel$expandNodeNames(parentNodes, returnScalarComponents = TRUE)
nodesToSim <- sirsModel$getDependencies(parentNodes, self = FALSE, downstream = T)

initsList <- list(gamma = 0.2,
                  phi = 0.1)
sirsModel$setInits(initsList)

set.seed(1)
sirsModel$simulate(nodesToSim, includeData = TRUE)

[Perry de Valpine]

Hi Marie,

Thanks for the question and especially the simplified and reproducible example.  That does make it much easier to take a look.

I see two things going on.  The first is that there is some kind of glitch with using the constant 'pi' in the model code.  This shouldn't be the case, so we need to take a look into that.  As a workaround, please replace it with something like "my_pi" and then in your constantsList include "my_pi = pi", so the actual value of pi is inserted into the model as a constant without relying on the special name "pi".

Then I see that even with that workaround, your code throws errors.  I am not completely following what you want to end up with from the model structure queries, but what I can see is that "Istar[1]" is not included in nodesToSim, and having that not simulated seems to set off a cascade of NAs through the time steps of the model.  I was just inspecting with code like:

sirsModel$Y[1]
sirsModel$nu[1]
sirsModel$Istar[1]
sirsModel$S[1]
sirsModel$probSI[1]
sirsModel$w[1]
sirsModel$beta[1]
sirsModel$I[1]
constantsList$w[1]

This issue was undoubtedly somewhat masked from your debugging efforts by the problem with pi, so I am hoping you'll be able to work out what you want from here.  I wonder if some of the options in model$getNodeNames would be helpful, e.g. "topOnly = TRUE" will give all the nodes with no stochastic parents.  Is something like setdiff(sirsModel$getNodeNames(), c("gamma", "phi")) useful?

HTH.

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/29dd0a2e-e684-4d39-a378-a53f1490bc57n%40googlegroups.com.

[Marie Ozanne]

Thank you for your response, Perry! I explicitly included "Istar[1]" in nodesToSim and then the simulation worked with no issues. I am still puzzled as to why just "Istar[1]" was excluded in nodesToSim in the given setup. Is there a straightforward answer to this?

Thank you for your help!

Marie 

[Perry de Valpine]

Hi Marie,

If I follow, after this step

parentNodes <- parentNodes[-which(parentNodes %in% dataNodes)]

parentNodes contains "gamma" and "phi".  And then Istar[1] does not depend on gamma or phi, so it is not included in 

nodesToSim <- sirsModel$getDependencies(parentNodes, self = FALSE, downstream = T)

Rstar[1] and Sstar[1] depend on gamma and phi, and then via those nodes Istar[2] and later ones are all in the downstream results.  It's just that Istar[1] starts the sequence and does not depend on gamma or phi.

I need to amend my reply about pi. I was forgetting that the choice to not support pi as a constant in model code, the way it is in R or nimbleFunction code, was deliberate because pi is not uncommonly used as a variable name in models. So instead what I think we should look into is providing a helpful error or warning message.  If one follows the output of nimbleModel and checks sirsModel$initializeInfo(), one can see that pi is indeed part of what is unknown, and this can be a flag that it is not supported as a constant.  But it is somewhat buried.

Best wishes,

Perry