Problem with initializing first occurences in CJS multistate (age-class) model

[Oldřich Tomášek]

Hi all,

I am trying to fit a relatively simple CJS multistate model for survival - specifically modelling the transition to the next age class or death.

My states (z matrix) and observations (y matrix) are as follows:

# states (z):
# 1-6 - age 1-6 years
# 7 - age 7 years or older
# 8 - dead
# observations (y):
# 1-7 - observed as 1-7+ years old
# 8 - unobserved

The complete code is at the end of my post.

When running the model, I receive the following warnings:

problem initializing stochastic node z[i, j]: logProb is NA or NaN

categorical sampler for 'z[m, n]' encountered an invalid model density, and sampling results are likely invalid

When checking the positions indexed by i and j in the transition matrix z in the first warning type, these refer to the first captures of the individuals.

The indexes m and n in the second warning refer to the first captures too, but only to those when the individual was first captured as 2 years old instead of 1 year old.

Although the model is quite simple, I cannot figure out what could cause the problem. Could my initial values for z matrix be the cause?

I am providing the code below and the data set is attached.

I would greatly appreciate any help.

Many thanks,

Olda

library(nimble)
library(MCMCvis)
library(basicMCMCplots)


y <- as.matrix(read.csv(file = "data/y_test.csv"))

get.first <- function(x) min(which(x != 8))
first <- apply(y, 1, get.first)


#### the model code ####
# states (Z):
# 1-6 - age 1-6 years
# 7 - age 7 years and older
# 8 - dead
# observations (y):
# 1-7 - observed as 1-7+ years old
# 8 - unobserved

hmm_age_as_known <- nimbleCode({
  # priors
  pi1 ~ dunif(0, 1)
  pi2 ~ dunif(0, 1)
  pi3 ~ dunif(0, 1)
  pi4 ~ dunif(0, 1)
  pi5 ~ dunif(0, 1)
  pi6 ~ dunif(0, 1)

  phi1 ~ dunif(0, 1)
  phi2 ~ dunif(0, 1)
  phi3 ~ dunif(0, 1)
  phi4 ~ dunif(0, 1)
  phi5 ~ dunif(0, 1)
  phi6 ~ dunif(0, 1)
  phi7 ~ dunif(0, 1)

  p1 ~ dunif(0, 1)
  p2 ~ dunif(0, 1)
  p3 ~ dunif(0, 1)
  p4 ~ dunif(0, 1)
  p5 ~ dunif(0, 1)
  p6 ~ dunif(0, 1)
  p7 ~ dunif(0, 1)

  # initial state probs
  delta[1:8] <- c(pi1, pi2, pi3, pi4, pi5, pi6, 1-pi1-pi2-pi3-pi4-pi5-pi6, 0)
 
  # transition matrix - probabilities of state z(t+1) given z(t)
  gamma[1, 1:8] <- c(0, phi1, 0, 0, 0, 0, 0, 1-phi1) # transition from age 1
  gamma[2, 1:8] <- c(0, 0, phi2, 0, 0, 0, 0, 1-phi2) # transition from age 2
  gamma[3, 1:8] <- c(0, 0, 0, phi3, 0, 0, 0, 1-phi3) # transition from age 3
  gamma[4, 1:8] <- c(0, 0, 0, 0, phi4, 0, 0, 1-phi4) # transition from age 4
  gamma[5, 1:8] <- c(0, 0, 0, 0, 0, phi5, 0, 1-phi5) # transition from age 5
  gamma[6, 1:8] <- c(0, 0, 0, 0, 0, 0, phi6, 1-phi6) # transition from age 6
  gamma[7, 1:8] <- c(0, 0, 0, 0, 0, 0, phi7, 1-phi7) # transition from age 7
  gamma[8, 1:8] <- c(0, 0, 0, 0, 0, 0, 0, 1) # transition from dead
 
  # observation matrix - probabilities of y(t) given z(t)
  # p1-p7 - probability of being detected given the specific age group
  omega[1, 1:8] <- c(p1, 0, 0, 0, 0, 0, 0, 1-p1)
  omega[2, 1:8] <- c(0, p2, 0, 0, 0, 0, 0, 1-p2)
  omega[3, 1:8] <- c(0, 0, p3, 0, 0, 0, 0, 1-p3)
  omega[4, 1:8] <- c(0, 0, 0, p4, 0, 0, 0, 1-p4)
  omega[5, 1:8] <- c(0, 0, 0, 0, p5, 0, 0, 1-p5)
  omega[6, 1:8] <- c(0, 0, 0, 0, 0, p6, 0, 1-p6)
  omega[7, 1:8] <- c(0, 0, 0, 0, 0, 0, p7, 1-p7)
  omega[8, 1:8] <- c(0, 0, 0, 0, 0, 0, 0, 1)

  # likelihood
  for (i in 1:N){
    # latent and observed states at the first capture
    z[i, first[i]] ~ dcat(delta[1:8])
    for (t in (first[i]+1):K){
      # z(t) given z(t-1)
      z[i, t] ~ dcat(gamma[z[i, t-1], 1:8])
      # y(t) given z(t)
      y[i, t] ~ dcat(omega[z[i, t], 1:8])
    }
  }
})

my_constants <- list(first = first,
                     K = ncol(y),
                     N = nrow(y))

my_data <- list(y = y)

zinit <- y
for (i in 1:nrow(y)) {
  if (first[i] > 1 & zinit[i, first[i]] > 1 & first[i] <= zinit[i, first[i]]) {zinit[i, 1] <- zinit[i, first[i]]-first[i]+1}
    for (j in 2:ncol(y)) {
    if (j < first[i] & zinit[i, first[i]] > 1 & j > first[i] - zinit[i, first[i]]) {zinit[i, j] <- zinit[i, first[i]]-first[i]+j}
    if (j > first[i] & zinit[i, j-1] < 7 & zinit[i, j] == 8) {zinit[i, j] <- zinit[i, j-1]+1}
    if (j > first[i] & zinit[i, j-1] == 7) {zinit[i, j] <- 7}
  }
}
zinit <- as.matrix(zinit)

initial_values <- function() list(
  pi1 = runif(1, 0, 1),
  pi2 = runif(1, 0, 1),
  pi3 = runif(1, 0, 1),
  pi4 = runif(1, 0, 1),
  pi5 = runif(1, 0, 1),
  pi6 = runif(1, 0, 1),
  phi1 = runif(1, 0, 1),
  phi2 = runif(1, 0, 1),
  phi3 = runif(1, 0, 1),
  phi4 = runif(1, 0, 1),
  phi5 = runif(1, 0, 1),
  phi6 = runif(1, 0, 1),
  phi7 = runif(1, 0, 1),
  p1 = runif(1, 0, 1),
  p2 = runif(1, 0, 1),
  p3 = runif(1, 0, 1),
  p4 = runif(1, 0, 1),
  p5 = runif(1, 0, 1),
  p6 = runif(1, 0, 1),
  p7 = runif(1, 0, 1),
  z = zinit)

parameters_to_save <- c(
  "pi1",
  "pi2",
  "pi3",
  "pi4",
  "pi5",
  "pi6",
  "phi1",
  "phi2",
  "phi3",
  "phi4",
  "phi5",
  "phi6",
  "phi7",
  "p1",
  "p2",
  "p3",
  "p4",
  "p5",
  "p6",
  "p7"
)

# MCMC settings
n_iter <- 5000
n_burnin <- 2500
n_chains <- 2

mcmc_age_as_known <- nimbleMCMC(
  code = hmm_age_as_known,
  constants = my_constants,
  data = my_data,
  inits = initial_values,
  monitors = parameters_to_save,
  niter = n_iter,
  nburnin = n_burnin,
  nchains = n_chains)

[Chris Paciorek]

Hi Oldřich ,

It could indeed be an initial value problem. If you haven't already, please take a look at the section in the NIMBLE manual on initializing MCMC. That 

has some suggestions for how to narrow down to the node(s) that may have NaN or -Inf values as their log probability. If you can identify those

then you may be able to figure out if some initial value(s) are invalid.

Let us know if that doesn't help you get unstuck.

-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/2c020773-33c3-4523-a6f7-9e69b84e1b9an%40googlegroups.com.

[Radovan Václav]

Hi Olda,

I think your problem is related to defining the latent state at the first capture.

# latent and observed states at the first capture
    z[i, first[i]] ~ dcat(delta[1:8])

You may want to try to define the latent state at the first capture this way: 

z[i, first[i]] <- y(i, first[i])

Rado

[Perry de Valpine]

Also just pointing you to nimbleEcology and the vignette there in case they are helpful.

-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/e1a29299-4a44-40ee-a294-b2b92e433093n%40googlegroups.com.