errors in visualizing results

[cecilia soldatini]

Hi all,

I am running CJS models with and without covariates and with or without age effect, all of them run smoothly and return a WAIC value but when trying to get for example:

MCMCsummary(mcmc.phitrat, round = 2)

MCMCplot(mcmc.phitrat)

MCMCtrace(object = mcmc.phitrat, params = "phi", pdf = FALSE)

the error is invariably:

Error in coda::mcmc.list(lapply(object, function(x) coda::mcmc(x))) : Different start, end or thin values in each chain

for all the models.

Any tip on how to solve this?

thanks,

Cecilia

[Chris Paciorek]

hi Cecilia,

Can you show the code you are using to run the MCMC and produce mcmc.phitrat so we can see more about the number of chains and the start/end/thin/burnin for the chains?

-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/da88d97e-07ad-45c7-8f2f-ea3962f21a3dn%40googlegroups.com.

[cecilia soldatini]

Hi Chris,

thanks for your email, here is the code:

#'## CJS models constant detection survival varying with treatment
##--------------------CJS trat-------------------------------------------------------------------

hmm.phitrat <- nimbleCode({
  delta[1] <- 1          # Pr(alive t = 1) = 1
  delta[2] <- 0          # Pr(dead t = 1) = 0
  for (t in 1:(T-1)){
    logit(phi[t]) <- beta[1] + beta[2] * trat[t]
    gamma[1,1,t] <- phi[t]      # Pr(alive t -> alive t+1)
    gamma[1,2,t] <- 1 - phi[t]  # Pr(alive t -> dead t+1)
    gamma[2,1,t] <- 0           # Pr(dead t -> alive t+1)
    gamma[2,2,t] <- 1           # Pr(dead t -> dead t+1)
  }
  p ~ dunif(0, 1) # prior detection
  omega[1,1] <- 1 - p    # Pr(alive t -> non-detected t)
  omega[1,2] <- p        # Pr(alive t -> detected t)
  omega[2,1] <- 1        # Pr(dead t -> non-detected t)
  omega[2,2] <- 0        # Pr(dead t -> detected t)
  beta[1] ~ dnorm(0, 1.5) # prior intercept
  beta[2] ~ dnorm(0, 1.5) # prior slope
  # likelihood
  for (i in 1:N){
    z[i,first[i]] ~ dcat(delta[1:2])
    for (j in (first[i]+1):T){
      z[i,j] ~ dcat(gamma[z[i,j-1], 1:2, j-1])
      y[i,j] ~ dcat(omega[z[i,j], 1:2])
    }
  }
})


#'
#' Constants in a list.
## ---------------------------------------------------------------
first <- apply(y, 1, function(x) min(which(x !=0)))
trat <- (data1$trat)
trat

my.constants <- list(N = nrow(y),
                     T = ncol(y),
                     first = first,
                     trat= trat)

#'
#' Data in a list.
## ---------------------------------------------------------------
my.data <- list(y = y + 1)

#'
#' Initial values.
## ---------------------------------------------------------------
zinits <- y
zinits[zinits == 0] <- 1
initial.values <- function() list(beta = rnorm(2,0,5),
                                  p = runif(1,0,1),
                                  z = zinits)

#'
#' Parameters to be monitored.
## ---------------------------------------------------------------
parameters.to.save <- c("phi", "p", "z", "beta")
parameters.to.save
#'
#' MCMC details.
## ---------------------------------------------------------------
n.iter <- 5000
n.burnin <- 2500
n.chains <- 2

#'
#' Run nimble.
## ---------------------------------------------------------------
mcmc.phitrat <- nimbleMCMC(code = hmm.phitrat,
                           constants = my.constants,
                           data = my.data,              
                           inits = initial.values,
                           monitors = parameters.to.save,
                           niter = n.iter,
                           nburnin = n.burnin,
                           nchains = n.chains,
                           WAIC = TRUE)

Matt already suggested to subset the MCMC object to refer to the sample, so mcmcphitrat$samples and it worked.

Nonetheless, it gives a Phi for every sample in the graph (see attached). Iraida suggested using set.seed() to say how many phi I want to represent in the graph, but for the moment it didn't work, probably I did something wrong, I'll keep trying.

thanks a lot,

cheers

Cecilia

Remitente notificado con Mailtrack

11/02/22 09:01:38

--

Cecilia Soldatini, Ph.D.

CICESE - Unidad La Paz 

tel: +52 1 6121070538

Skype: cecilia.soldatini

OrcID: 0000-0002-8112-3162 

ResearcherID: F-2950-2011

Research Gate

Scopus

Grupo de Aeroecologia Marina

[Chris Paciorek]

Hi Cecilia,

Yes, the issue here was that if you ask for WAIC then nimbleMCMC returns a list that contains both the samples and the WAIC. In that case when using the MCMCvis functions you need to pass the 'samples' component of the list. If you hadn't asked for WAIC, it would have worked to pass the output of nimbleMCMC.  We're trying to have nimbleMCMC be user-friendly in terms of simply returning the samples if that is all that is asked for, but it can make it a bit confusing what format the output of nimbleMCMC is.

I don't understand what you are saying when you say "Nonetheless, it gives a Phi for every sample in the graph (see attached). Iraida suggested using set.seed() to say how many phi I want to represent in the graph, but for the moment it didn't work"

 - if you monitor 'phi', which has multiple elements, then you'll get back a column of samples for each of the elements of phi. 

 - set.seed has to do with setting the random number seed. It's not related to anything about the model graph. The number of elements of phi is determined by the "[t]" indexing.

-chris

[cecilia soldatini]

Hi Chris,

the problem is that in the graph (and in the results) when considering the effect of a covariate in the model, "trat" in my case, it returns the same value of phi for every sample. So instead of plotting one value of phi,  it gives a series of dots for phi which are all of the same value, like in the graph I was attaching in my previous message. 

If I use a multistate model I have a phi value for each sample and different states have their phi but it is the same value for every sample of the state.

I don't know if I explained clearly, if you want I can send my script and DB.

Thank you, 

best

Cecilia

Remitente notificado con Mailtrack

12/02/22 19:27:50

[Chris Paciorek]

hi Cecilia,

Are you sure that the values in 'trat' are not all the same? Given that beta[1] and beta[2] are scalars, phi should only vary if 'trat' varies.

In general in such a situation one would usually be monitoring beta[1] and beta[2]  since those are parameters, while 'phi' is just the linear predictor that is

a deterministic function of beta[1], beta[2] and 'trat'. 

-chris

[cecilia soldatini]

Hi Chris,

yes you are right, "trat" is the treatment each group of animals is exposed to, so I'd better use beta[1] and beta [2] and "trat" as you suggested. 

Thanks a lot

Cecilia

Remitente notificado con Mailtrack

15/02/22 09:22:44