errors in vizualizing results + warnings error

[Céline Artero]

Hi,

I read lots of conversations to see if I could find a solution for the errors I get running nimble. It seems that people had similar problems than me however I wasn't able to solve mine by following the corresponding tips (code and data set at the end).

1. warnings while running the model. 

I have several warnings when I am running the model:

* Warning: R object of different size than NimArrDouble. R obj has size 1 but NimArrDbl has size 4.

Warning: R object of different size than NimArrDouble. R obj has size 1 but NimArrDbl has size 2.
Warning: R object of different size than NimArrDouble. R obj has size 1 but NimArrDbl has size 2.

* warning: value of deterministic node phi[278, 9]: value is NA or NaN even after trying to calculate.

* warning: problem initializing stochastic node z[672, 27]: logProb is -Inf.

* warning: problem initializing stochastic node z[278, 10]: logProb is NA or NaN.

The model results of p and phi are coherent, so it seems that the model is working. Are these warnings a problem or they are just informative ? Do you identify something wrong in my code or these warnings are due to my dataset? 

2. error in visualising results

I can't use the nimble functions to visualise my data.

Neither MCMCsummary(), MCMCtrace(), or plot() are working.

I already tried to look at mcmc.phitpt$samples but it still creates the following error message :

Error in quantile.default(newX[, i], ...) :

  missing values and NaN's not allowed if 'na.rm' is FALSE

It does not surprise me, I have NAs in my data set.

As I was unable to find the solution, I tried to calculate the summary by myself by calculating the mean, sd, quantile of the model output. However, I didn't find any solution to calculate Rhat and n.eff.

I found some function on internet that could have helped but again, I am having errors.

I tried for n.eff:

coda::effectiveSize(mcmc.phitpt$samples)

for Rhat:

coda::gelman.diag(mcmc.phitpt, confidence = 0.95, transform=FALSE, autoburnin=TRUE,
                  multivariate=TRUE)

Error message is the same for both:

Error in mcmc.list(x) : Arguments must be mcmc objects

I converted mcmc.phitpt$samples into a mcmc object using coda::as.mcmc() but I keep having the same error.

Could you help please ?

My code and data :

My data set is a table of 6+28 columns (corresponding to 6 variables + 28 sampling occasions) and 730 individuals.

Each individual was tagged in a specific site (4 sites), at a specific year (2 years). There are 6 or 8 receivers on each site.

I was only able to develop the model by having a specific column for each receiver, therefore I have 28 columns in total with NAs (I created a summarised my data set in joined file - blanks are NAs).

I also tried to have only 8 columns to avoid too much NAs but then I would have had p[site (s), receiver (r), year (t)] instead of p[r, t] that complicated again the extraction of the results. 

Thank you in advance for your help. It would be very helpful as I tried everything I was thinking of.

Céline

### IMPORT DATA ----

Estuary <- read.csv("data set.csv", sep = ",")

site <- Estuary$River
sex <- Estuary$Sex
species <- Estuary$Species
size <- Estuary$ZSize
WL <- Estuary$ZresWL

y <- select(Estuary, -c(1:6))
y <- as.matrix(y)

### MODEL HMM VF ----
#---

phitpt <- nimbleCode({
 
  # coefficients
  alpha ~ dnorm(mean = 0, sd = 1.5)
  beta_site[1] <- 0
  beta_site[2] ~ dnorm(0, 1.5)
  beta_site[3] ~ dnorm(0, 1.5)
  beta_site[4] ~ dnorm(0, 1.5)
  beta_species[1] <- 0
  beta_species[2] ~ dnorm(0, 1.5)
  beta_sex[1] <- 0
  beta_sex[2] ~ dnorm(0, 1.5)
  beta_size ~ dnorm(0, 1.5)
  beta_WL ~ dnorm(0, 1.5)
  beta_dist ~ dnorm(0, 1.5)
  delta[1] <- 1
  delta[2] <- 0
 
   
  # phi: coded to vary in between receivers (r)
  for (i in 1:N) {
    for (r in first[i]:last[i]-1) {
      logit(phi[i, r]) <- alpha + beta_site[site[i]] + beta_species[species[i]] + beta_sex[sex[i]] +
                          (beta_size * size[i]) + (beta_WL * WL[i]) 
    }
  }
 
  # p: vary between receivers (r) and year of study (t)
  for (r in 1:T-1){
    for (t in 1:2){
      p[t, r] ~ dunif(0,1)
    }
  }
 

  # transition matrices
  for (i in 1:N){
    for (r in first[i]:last[i]-1){
     
      gamma[1, 1, i, r] <- phi[i, r]            # alive at t-1 & alive at t
      gamma[1, 2, i, r] <- 1 - phi[i, r]        # alive at t-1 & dead at t
      gamma[2, 1, i, r] <- 0                    # dead at t-1 & alive at t
      gamma[2, 2, i, r] <- 1                    # dead at t-1 & dead at t
     
      omega[1, 1, i, r] <- 1 - p[year[i], r]       # alive & detected at rec a
      omega[1, 2, i, r] <- p[year[i], r]   # alive & not detected at rec a
      omega[2, 1, i, r] <- 1                   # dead & detected at rec a
      omega[2, 2, i, r] <- 0                    # dead & not detected at rec r
     
    }
  }
 
  # state-space model likelihood
  for (i in 1:N) {
    z[i, first[i]] ~ dcat(delta[1:2])
    for (r in (first[i]+1):last[i]) {
      z[i, r] ~ dcat(gamma[z[i, r - 1], 1:2, i, r - 1])
      y[i, r] ~ dcat(omega[z[i, r], 1:2, i, r - 1])
    }
  }
 })

## data

my.data <- list(y = y + 1)

## constants
first <- apply(y, 1, function(x) min(which(x !=0)))
last <- apply(y, 1, function(x) max(which(! is.na(x))))

my.constants <- list(N = nrow(y),
                     T = ncol(y),
                     site = site,
                     year = year,
                     first = first,
                     size = size,
                     species = species,
                     sex = sex,
                     WL = WL,
                     last = last)

## Set initial values
zinits <- y + 1 # non-detection -> alive
zinits[zinits == 2] <- 1 # dead -> alive
initial.values <- function() list(alpha = rnorm(1,0,1),
                                  beta_site = rnorm(1,0,1),
                                  beta_species = rnorm(1,0,1),
                                  beta_sex = rnorm(1,0,1),
                                  beta_size = rnorm(1,0,1),
                                  beta_WL = rnorm(1,0,1),
                                  beta_dist = rnorm(1,0,1),
                                  z = zinits)

# Define parameters to save
#---
parameters.to.save <- c("phi","p", "z", "beta_site", "beta_species", "beta_sex", "beta_size",
                        "beta_WL", "beta_dist")

## mcmc details
n.iter <- 5000
n.burnin <- 1000
n.chains <- 2
n.thin <- 1

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

[Daniel Turek]

Céline, thanks for the question and the detailed information.  There are some issues with initial values in your model (both the sizes in the initial.values function, and also the values for z.init), which is causing the problems with looking at the samples.

Because of the initial values issues, some entries in the 1st row of your MCMC samples (first iteration of the MCMC) are NA.  Therefore, trying to find quantiles (without using na.rm = TRUE) will give the error you're seeing.  You could include na.rm = TRUE, or fix the initialization issues discussed below.

As for the error from coda::gelman.diag function, this function requires an coda::mcmc.list object, where each element is a coda::mcmc object.  You can re-cast the list of samples arrays yourself, or easier, supply this argument to nimbleMCMC: samplesAsCodaMCMC = TRUE.  That said, if the samples arrays still contain NA's then the gelman.diag function will likely still fail.

So, it might be in your best interests to clear up the initialization issues:  That said, I haven't gone ahead and fixed these issues yet, so I'm not sure whether fixing these will take care of everything.  But for starters:

(1) In your initial.values function, note that for beta_site, beta_species, and beta_sex, you're giving initial values of the wrong dimensions.  For example, in your model code:

  beta_sex[1] <- 0
  beta_sex[2] ~ dnorm(0, 1.5)

even though beta_sex[1] is hard-coded as 0, the beta_sex variable is still length-2.  Thus, an appropriate initial value would be:

beta_sex = c(0, rnorm(1, 0, 1))

Similar for the length-4 beta_site variable, where an appropriate initial value might be:

beta_site = c(0, rnorm(3, 0, 1))

And, also for the length-2 beta_species.

(2) The second issue is causing the "dynamic index out of bounds" errors that you're seeing.  In particular, the problem is for individuals with NA's in the capture history *between* non-NA observations, for example:

> unname(y[23,])
 [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA  1  1 NA  1  1  1  1  1 NA NA NA
[26] NA NA NA

and

> unname(y[40,])
 [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA  1  1 NA
[26]  1 NA NA

 Among others, individuals # (24, 25, 26, 27, 28, 29, 30, 46).

In the setup of the initial values for z, zinit:

zinits <- y + 1 # non-detection -> alive
zinits[zinits == 2] <- 1 # dead -> alive

These NA's in y translate to NA's in zinit, which is then later used for indexing the gamma[] and omega[] arrays.  These missing values of zinit are taking a while to recover from.  I would suggest you change how zinit is constructed, so it has valid (non-NA) initial values for all values of z which will be used for indexing either gamma[] or omega[], so that is for all values of z[i, first[i]:last[i]].    This should be a pretty easy change, with a few lines of code.

Assuming these changes work, then you should *not* see all the warning messages you're seeing, and I think that will clear up other problems.

Give it a go, and let me know if it works, or you have other questions.

Cheers,

Danel

--
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/797afd4c-4033-4c2f-8698-18c17fcf3242n%40googlegroups.com.