Issue with dynamic indexing in nimbleCode and nimbleMCMC

[Anthony S]

Hello everyone, 

I am trying to set up a model in which I want to save one of my parameters using dynamic indexing. I built the original code from a JAGS script (which works), but I am doing something wrong when translating it to Nimble, and get the following error when running the NimbleMCMC() = “ Dynamic index out of bounds”.

 Below is a (I hope) reproducible example: 

library(plyr)

library(nimble
set.seed(12345)

Let's imagine that I'm sampling individuals in a population, and recording their birth year. I am then doing pairwise comparisons between all individuals, and counting the frequency of the different Ind1 x Ind2 birth year combinations.

Ind_sampled <- floor(runif (200, 1, 15))

df_Ind <- expand.grid(Ind_sampled, Ind_sampled)
names (df_Ind) <- c("Ind1", "Ind2")

counts <- ddply(df_Ind, .(df_Ind$Ind1, df_Ind$Ind2), nrow)
names(counts) <- c("Ind1", "Ind2", "Sampled")

Total <- floor(runif(nrow(counts), 1000, 5000))

df <- cbind(counts, Total)

# Ind1 = Ind1 birth year
# Ind2 = Ind 2 birth year
# Sampled = How many were sampled
# Total = number of combinations for the entire population (sampled and not sampled)

The model uses the Sampled/Total ratio for each line (ie each birth year combination) to calculate N, the total number of individuals in the population.
I show the "simple" model below, where I estimate an average value for N over all birth years of Ind 2.

I then show the dynamic model, in which I want to estimate a unique value for N for each birth year of Ind 2. I attach the rest of the code to this dynamic model.

#### SIMPLE CKMR MODEL (works well, but boring) ####

Model <- nimbleCode({
 
  # priors
 
N ~ dunif(0, 1e5) # Uninformative prior for N
 
  # likelihoods
  for (j in 1:df_length){
    SAMP[j] ~ dbinom((1/N), Total[j])
  }
})

#### DYNAMIC CKMR MODEL ####

Nyears = sort(unique(df$Ind2))
Nindex = match(df$Ind2, Nyears)
Nlen = length(unique(Nindex))

Model <- nimbleCode({
 
  # priors
   for(i in 1:Nlen)
  {
    N[i] ~ dunif(0, 1e5) # Uninformative prior for N
  }
 
   # likelihoods
  for (j in 1:df_length)

{
    SAMP[j] ~ dbinom((1/N[Nindex[j]]), Total[j])
  }
})

###---We put the constants in a list---###

my.constants <- list(df_length = nrow(df))

###---We put the data in a list---###

my.data <- list(SAMP = df$Sampled, # Sampled combinations
                Total = df$Total) # Total combinations (hypothetical)

###---We write a function that generates initial values---###

initial.values <- function() list(N = rnorm(1, mean = 2000, sd = 100))

###---We specify the variables we would like to monitor---###

parameters.to.save <- c("N")

###---We specify a few details about the chains---###

n.iter <- 5000
n.burnin <- 1000
n.chains <- 2
n.thin <- 10

###---And now, we run nimble---###

mcmc.output <- nimbleMCMC(code = Model,
                          constants = my.constants,
                          data = my.data,
                          inits = initial.values,
                          monitors = parameters.to.save,
                          niter = n.iter,
                          thin = n.thin,
                          nburnin = n.burnin,
                          nchains = n.chains)


Ultimately, I wish to have a result that shows a value for C for each Nindex[]
(e.g. N[1] = 6000, N[2] = 3500, ... , N[14] = 3110)

Should I specify dimensions for N in the parameters to save? Initial values? Both? I'm at loss.

Best regards,

Anthony

[Frédéric LeTourneux]

Hello Anthony,

I have not completely gotten the model you're trying to fit, but two 'mechanic' things come to mind when I look at your code. First, if you are using Nindex and Nlen in your nimble code, you must give these vectors to nimble as constants, it is not sufficient that they simply exist in your environment. This should deal with your indexing problem.

The other thing I notice is you're using a dunif(0, 1e5) as a prior for your N. If you are not specifying starting values for N (in your code you are specifying 1 initial value with N = rnorm(1, mean = 2000, sd = 100), but I think this wont specify 1 value for each of your 14 N nodes), then  Nimble generates initial values by simulating from the prior you are giving, (i.e. dunif(0,1e5)).  Technically (although highly unlikely give the large range of possibilities), this means that you could have a 0 as a starting value for N, which you are then using as the denominator of the probability of your binomial distribution. That can't be good... I think. Anyway,  if you want to initialize all N nodes, you should ask for 14 values in your inits function, i.e. N = rnorm(14, mean = 2000, sd = 100).

I'm not yet an expert with Nimble so others might want to weigh in on this, but at least with these two modifications, your code seems to run and generate samples.

Hope that helps!

F

[Anthony S]

Hi Frédéric, 

These work perfectly, thank you for pointing out my errors. I cannot believe I'm still forgetting to give my variables to nimble as constants/data!

Anthony