Developping an age-dependant Hidden Multistate Model in Nimble from scratch
184 views
Skip to first unread message
Etienne Rouby
Aug 17, 2023, 10:46:04 AM8/17/23
to nimble-users
Hello Everyone,
First, I would like to thank you for reading this message. And also the Nimble team to develop such a tool. Which is very powerful for addressing ecological questions.
Mine is to develop a model to estimate the Age at First reproduction in a Long-lived Seabird species. To do so, I have to estimate several parameters, like the survival and reproduction probabilities given multiple Breeding states (PB: Pre-Breeeders, B: Breeders, NB: Non-Breeders), in an age-dependent manner. I choose to use a categorical modeling specification to include age (I may move to a categorial or linear one later, but at the moment, the categorical is good).
I developed some first models that are simple. But now it is necessary to include more complex buildings. And I am a bit lost since I get some warnings but do not know where they come from. And I get some very weird estimations from my parameters. So I decided to simulate some data to investigate how my model behaves to estimate parameters from what I know the value.
At the moment, I did not simulate age-dependent parameters, I would have to, but I do not know how to do. So here is my questions :
1) How to include age-dependent probabilities in my simulated data? For example, say that psiPB should be 0 between 1 and 5 years, then 0.5 between 5 and 10, then 0.1 between 10 and 12+
2) Why do I get these warnings: warning: logProb of data node y[101, 38]: logProb is -Inf.
3) Do you think I approach the modelisation stuff well?
Thank you very much for your help. I would be happy to give you any precision regarding the model.
Here is a reproducible code that turns into 1min. I annotated it as maximum as I could:
######################
######################
# Remove all objects from the environment
rm(list = ls())
# Load the relevant packages
lapply(c("tidyverse", "dplyr", "tidyr", "nimble", "MCMCvis"), library, character.only = TRUE)
## Set parameter values to generate fake data
n_ind <- 200 # number of individuals to generate
n_year <- 40 # number of occasions of samplig
# Data simulation --------------------------------------------------------------------------------------------
## Define the parameters of the transition matrix
phiPB <- 0.9 # Survival probability of birds that are in a Pre-Breeding State (ex: Chicks and Juveniles)
phiB <- 0.95 # Survival probability of birds that are in a Breeding State (Breeding Adults)
phiNB <- 0.92 # Survival probability of birds that are in a Non-Breeding State (Non-Breeding Adults)
psiPB <- 0.1 # Reproduction probability of birds that are in a Pre-Breeding State (ex: Chicks and Juveniles)
psiB <- 0.3 # Reproduction probability of birds that are in a Breeding State (Breeding adults). It is low because birds reproduce then skip the next year of reproduction (there are more probability for them to become Non-Breeder after beeing Breeder)
psiNB <- 0.9 # Reproduction probability of birds that are in a Non-Breeding State (Non-Breeding Adults). Non Breeder birds have an High probability of becoming Breeders the next year.
p <- 0.8 # Detection probability which is the same for each state.
first <- rep(1, n_ind) # First sampling occasion. Here assume a unique cohort. After, I will have to assume several
# Define the transition Marrix
A <- matrix(c(phiPB * (1-psiPB), phiPB * psiPB, 0, 1-phiPB,
0, phiB * psiB, phiB * (1-psiB), 1-phiB,
0, phiNB * psiNB, phiNB * (1-psiNB), 1-phiNB,
0, 0, 0, 1), nrow=4, byrow=TRUE)
namstates <- c("PB","B","NB","D")
dimnames(A) <- list(origin=namstates,destination=namstates)
# Simulate life trajectories for each individuals within the cohort
# Individuals starts as Pre Breeder. They become Breeder with a certain probability. Then, they transit between Breeder and Non-Breeder states each year.
# Ideally, the probabilities should be aged dependant. Because It will be the case in Nature.
z <- matrix(NA, nrow = n_ind, ncol = n_year)
for (i in 1:n_ind){ # loop over individuals
z[i,first[i]] <- 1 # all individuals are alive at first occasion
for (t in (first[i]+1):n_year){ # loop over time
state_draw <- as.data.frame(rmultinom(1, 1, A[z[i,t-1], 1:4])) %>% # draw the state at t given the state at t-1 and transition matrix
tibble::rownames_to_column("state")
if (filter(state_draw, V1 == 1)$state == "PB") { # Assign the good index depending on the state
z[i,t] <- 1
}
if (filter(state_draw, V1 == 1)$state == "B") {
z[i,t] <- 2
}
if (filter(state_draw, V1 == 1)$state == "NB") {
z[i,t] <- 3
}
if (filter(state_draw, V1 == 1)$state == "D") { # D is for Dead
z[i,t] <- 4
}
}
}
# Simulate encounter histories for each individuals within the cohort
y <- matrix(NA, nrow = n_ind, ncol = n_year)
# Observation matrix
B <- matrix(c(p, 0, 0, 1-p,
0, p, 0, 1-p,
0, 0, p, 1-p,
0, 0, 0, 1), nrow=4, byrow=TRUE)
namstates <- c("PB","B","NB","ND") # ND is for "Non-Detected"
dimnames(B) <- list(origin=namstates,destination=namstates)
for (i in 1:n_ind){ # loop over individuals
y[i,first[i]] <- 1 # all individuals are alive at first occasion
for (t in (first[i]+1):n_year){ # loop over time
obs_draw <- as.data.frame(rmultinom(1, 1, B[z[i,t], 1:4])) %>% # The observation state at t depends on the living state (PB, B, NB or Dead) at t.
tibble::rownames_to_column("obs_state")
if (filter(obs_draw, V1 == 1)$obs_state == "PB") {
y[i,t] <- 1
}
if (filter(obs_draw, V1 == 1)$obs_state == "B") {
y[i,t] <- 2
}
if (filter(obs_draw, V1 == 1)$obs_state == "NB") {
y[i,t] <- 3
}
if (filter(obs_draw, V1 == 1)$obs_state == "ND") { # ND is for Non-Detected
y[i,t] <- 4
}
}
}
y[y == 4] <- 0 # Correct the for in order to have 0s, because individuals ar not observed.
# Obtain the age of each monitored individual if they are alive each year.
age <- matrix(NA, nrow = nrow(y), ncol = ncol(y)-1)
for (i in 1:nrow(age)){
for (j in 1:ncol(age)){
if (j == first[i]) age[i,j] <- 1
if (j == first[i]+1) age[i,j] <- 2
if (j == first[i]+2) age[i,j] <- 3
if (j == first[i]+3) age[i,j] <- 4
if (j == first[i]+4) age[i,j] <- 5
if (j == first[i]+5) age[i,j] <- 6
if (j == first[i]+6) age[i,j] <- 7
if (j == first[i]+7) age[i,j] <- 8
if (j == first[i]+8) age[i,j] <- 9
if (j == first[i]+9) age[i,j] <- 10
if (j == first[i]+10) age[i,j] <- 11
if (j >= first[i]+11) age[i,j] <- 12 # I stop at 12 because from 12 to Death, I consider the individuals belong to the same age class ( >= 12)
}
}
# Parameter estimation from simulated data ---------------------------------------------------------------------------------------------------------
# Nimble model
multistate.cate <- nimbleCode({
# -------------------------------------------------
# Parameters:
# phiB: survival probability state B
# phiNB: survival probability state NB
# phiPB: survival probability state PB
# psiPB: reproduction probability PB
# psiB: reproduction probability B
# psiNB: reproduction probability NB
# pPB: recapture probability B
# pB: recapture probability B
# pNB: recapture probability NB
# -------------------------------------------------
# States (z):
# 1 alive PB
# 2 alive B
# 3 alive NB
# 4 dead
# Observations (y):
# 1 not seen
# 2 seen as PB
# 3 seen as B
# 4 seen as NB
# -------------------------------------------------
# Priors
pPB ~ dunif(0,1)
pB ~ dunif(0,1)
pNB ~ dunif(0,1)
# Transition and Observation matrixes
omega[1,1] <- 1-pNB
omega[1,2] <- 0
omega[1,3] <- 0
omega[1,4] <- pNB
omega[2,1] <- 1-pB
omega[2,2] <- 0
omega[2,3] <- pB
omega[2,4] <- 0
omega[3,1] <- 1-pPB
omega[3,2] <- pPB
omega[3,3] <- 0
omega[3,4] <- 0
omega[4,1] <- 1
omega[4,2] <- 0
omega[4,3] <- 0
omega[4,4] <- 0
for(i in 1:n_ind){
for(t in (first[i]):(n_year-1)){
logit(phiPB[i,t]) <- logit(Mu.phiPB[age[i,t]]) # Here I use logit to set up a model that I could use later to include covariates and random effects
logit(phiB[i,t]) <- logit(Mu.phiB[age[i,t]]) # Also, age is consider as a categorical variable at the moment. Maybe I will try to use a function latter
logit(phiNB[i,t]) <- logit(Mu.phiNB[age[i,t]])
logit(psiPB[i,t]) <- logit(Mu.psiPB[age[i,t]])
logit(psiB[i,t]) <- logit(Mu.psiB[age[i,t]])
logit(psiNB[i,t]) <- logit(Mu.psiNB[age[i,t]])
gamma[1,1,i,t] <- phiPB[i,t] * (1-psiPB[i,t])
gamma[1,2,i,t] <- phiPB[i,t] * psiPB[i,t]
gamma[1,3,i,t] <- 0
gamma[1,4,i,t] <- 1-phiPB[i,t]
gamma[2,1,i,t] <- 0
gamma[2,2,i,t] <- phiB[i,t] * psiB[i,t]
gamma[2,3,i,t] <- phiB[i,t] * (1-psiB[i,t])
gamma[2,4,i,t] <- 1-phiB[i,t]
gamma[3,1,i,t] <- 0
gamma[3,2,i,t] <- phiNB[i,t] * psiNB[i,t]
gamma[3,3,i,t] <- phiNB[i,t] * (1-psiNB[i,t])
gamma[3,4,i,t] <- 1-phiNB[i,t]
gamma[4,1,i,t] <- 0
gamma[4,2,i,t] <- 0
gamma[4,3,i,t] <- 0
gamma[4,4,i,t] <- 1
}
}
for(a in 1:Amax){
Mu.psiPB[a] ~ dunif(0, 1)
Mu.psiNB[a] ~ dunif(0, 1)
Mu.psiB[a] ~ dunif(0, 1)
Mu.phiB[a] ~ dunif(0, 1)
Mu.phiNB[a] ~ dunif(0, 1)
Mu.phiPB[a] ~ dunif(0, 1)
}
# Likelihood
for(i in 1:n_ind){
z[i,first[i]] <- y[i,first[i]] - 1
for(t in (first[i]+1):n_year){
z[i,t] ~ dcat(gamma[z[i,t-1],1:4,i,t-1])
y[i,t] ~ dcat(omega[z[i,t],1:4])
}
}
})
my.data <- list(y = y + 1)
my.constants <- list(first = first,
n_year = ncol(y),
n_ind = nrow(y),
Amax = max(age, na.rm = TRUE),
age = age)
zinits <- y
for (i in 1:ncol(zinits)){
if(first[i] > 1) zinits[i, 1:(first[i]-1)] <- NA
}
zinits[is.na(zinits)] <- 1
zinits[zinits == 0] <- sample(c(2,3), sum(zinits == 0), replace = TRUE)
initial.values <- function(){list(phiPB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
phiB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
phiNB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
psiPB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
psiB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
psiNB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
Mu.psiNB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.psiPB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.psiB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiNB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiPB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
z = zinits)}
parameters.to.save <- c("Mu.psiPB", "Mu.psiNB", "Mu.psiB", "Mu.phiPB", "Mu.phiB", "Mu.phiNB")
# parameters.to.save <- c("psiPB", "psiNB", "psiB", "phiPB", "phiB", "phiNB")
# Here there is a problem, when I estimate these parameters they are not estimated. I have NA for the first value and 0 for all the others.
parameters.to.save
n.iter <- 1000
n.burnin <- 500
n.chains <- 3
mcmc.multistate.cate <- nimbleMCMC(code = multistate.cate,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)
summary_2 <- MCMCsummary(mcmc.multistate.cate, round = 2)
### See chains
priors <- runif(3000, 0, 1)
MCMCtrace(object = mcmc.multistate.cate,
ISB = FALSE,
exact = TRUE,
params = c("Mu.psiPB[1]"),
pdf = FALSE,
priors = priors)
First, I would like to thank you for reading this message. And also the Nimble team to develop such a tool. Which is very powerful for addressing ecological questions.
Mine is to develop a model to estimate the Age at First reproduction in a Long-lived Seabird species. To do so, I have to estimate several parameters, like the survival and reproduction probabilities given multiple Breeding states (PB: Pre-Breeeders, B: Breeders, NB: Non-Breeders), in an age-dependent manner. I choose to use a categorical modeling specification to include age (I may move to a categorial or linear one later, but at the moment, the categorical is good).
I developed some first models that are simple. But now it is necessary to include more complex buildings. And I am a bit lost since I get some warnings but do not know where they come from. And I get some very weird estimations from my parameters. So I decided to simulate some data to investigate how my model behaves to estimate parameters from what I know the value.
At the moment, I did not simulate age-dependent parameters, I would have to, but I do not know how to do. So here is my questions :
1) How to include age-dependent probabilities in my simulated data? For example, say that psiPB should be 0 between 1 and 5 years, then 0.5 between 5 and 10, then 0.1 between 10 and 12+
2) Why do I get these warnings: warning: logProb of data node y[101, 38]: logProb is -Inf.
3) Do you think I approach the modelisation stuff well?
Thank you very much for your help. I would be happy to give you any precision regarding the model.
Best,
Etienne
Here is a reproducible code that turns into 1min. I annotated it as maximum as I could:
######################
######################
# Remove all objects from the environment
rm(list = ls())
# Load the relevant packages
lapply(c("tidyverse", "dplyr", "tidyr", "nimble", "MCMCvis"), library, character.only = TRUE)
## Set parameter values to generate fake data
n_ind <- 200 # number of individuals to generate
n_year <- 40 # number of occasions of samplig
# Data simulation --------------------------------------------------------------------------------------------
## Define the parameters of the transition matrix
phiPB <- 0.9 # Survival probability of birds that are in a Pre-Breeding State (ex: Chicks and Juveniles)
phiB <- 0.95 # Survival probability of birds that are in a Breeding State (Breeding Adults)
phiNB <- 0.92 # Survival probability of birds that are in a Non-Breeding State (Non-Breeding Adults)
psiPB <- 0.1 # Reproduction probability of birds that are in a Pre-Breeding State (ex: Chicks and Juveniles)
psiB <- 0.3 # Reproduction probability of birds that are in a Breeding State (Breeding adults). It is low because birds reproduce then skip the next year of reproduction (there are more probability for them to become Non-Breeder after beeing Breeder)
psiNB <- 0.9 # Reproduction probability of birds that are in a Non-Breeding State (Non-Breeding Adults). Non Breeder birds have an High probability of becoming Breeders the next year.
p <- 0.8 # Detection probability which is the same for each state.
first <- rep(1, n_ind) # First sampling occasion. Here assume a unique cohort. After, I will have to assume several
# Define the transition Marrix
A <- matrix(c(phiPB * (1-psiPB), phiPB * psiPB, 0, 1-phiPB,
0, phiB * psiB, phiB * (1-psiB), 1-phiB,
0, phiNB * psiNB, phiNB * (1-psiNB), 1-phiNB,
0, 0, 0, 1), nrow=4, byrow=TRUE)
namstates <- c("PB","B","NB","D")
dimnames(A) <- list(origin=namstates,destination=namstates)
# Simulate life trajectories for each individuals within the cohort
# Individuals starts as Pre Breeder. They become Breeder with a certain probability. Then, they transit between Breeder and Non-Breeder states each year.
# Ideally, the probabilities should be aged dependant. Because It will be the case in Nature.
z <- matrix(NA, nrow = n_ind, ncol = n_year)
for (i in 1:n_ind){ # loop over individuals
z[i,first[i]] <- 1 # all individuals are alive at first occasion
for (t in (first[i]+1):n_year){ # loop over time
state_draw <- as.data.frame(rmultinom(1, 1, A[z[i,t-1], 1:4])) %>% # draw the state at t given the state at t-1 and transition matrix
tibble::rownames_to_column("state")
if (filter(state_draw, V1 == 1)$state == "PB") { # Assign the good index depending on the state
z[i,t] <- 1
}
if (filter(state_draw, V1 == 1)$state == "B") {
z[i,t] <- 2
}
if (filter(state_draw, V1 == 1)$state == "NB") {
z[i,t] <- 3
}
if (filter(state_draw, V1 == 1)$state == "D") { # D is for Dead
z[i,t] <- 4
}
}
}
# Simulate encounter histories for each individuals within the cohort
y <- matrix(NA, nrow = n_ind, ncol = n_year)
# Observation matrix
B <- matrix(c(p, 0, 0, 1-p,
0, p, 0, 1-p,
0, 0, p, 1-p,
0, 0, 0, 1), nrow=4, byrow=TRUE)
namstates <- c("PB","B","NB","ND") # ND is for "Non-Detected"
dimnames(B) <- list(origin=namstates,destination=namstates)
for (i in 1:n_ind){ # loop over individuals
y[i,first[i]] <- 1 # all individuals are alive at first occasion
for (t in (first[i]+1):n_year){ # loop over time
obs_draw <- as.data.frame(rmultinom(1, 1, B[z[i,t], 1:4])) %>% # The observation state at t depends on the living state (PB, B, NB or Dead) at t.
tibble::rownames_to_column("obs_state")
if (filter(obs_draw, V1 == 1)$obs_state == "PB") {
y[i,t] <- 1
}
if (filter(obs_draw, V1 == 1)$obs_state == "B") {
y[i,t] <- 2
}
if (filter(obs_draw, V1 == 1)$obs_state == "NB") {
y[i,t] <- 3
}
if (filter(obs_draw, V1 == 1)$obs_state == "ND") { # ND is for Non-Detected
y[i,t] <- 4
}
}
}
y[y == 4] <- 0 # Correct the for in order to have 0s, because individuals ar not observed.
# Obtain the age of each monitored individual if they are alive each year.
age <- matrix(NA, nrow = nrow(y), ncol = ncol(y)-1)
for (i in 1:nrow(age)){
for (j in 1:ncol(age)){
if (j == first[i]) age[i,j] <- 1
if (j == first[i]+1) age[i,j] <- 2
if (j == first[i]+2) age[i,j] <- 3
if (j == first[i]+3) age[i,j] <- 4
if (j == first[i]+4) age[i,j] <- 5
if (j == first[i]+5) age[i,j] <- 6
if (j == first[i]+6) age[i,j] <- 7
if (j == first[i]+7) age[i,j] <- 8
if (j == first[i]+8) age[i,j] <- 9
if (j == first[i]+9) age[i,j] <- 10
if (j == first[i]+10) age[i,j] <- 11
if (j >= first[i]+11) age[i,j] <- 12 # I stop at 12 because from 12 to Death, I consider the individuals belong to the same age class ( >= 12)
}
}
# Parameter estimation from simulated data ---------------------------------------------------------------------------------------------------------
# Nimble model
multistate.cate <- nimbleCode({
# -------------------------------------------------
# Parameters:
# phiB: survival probability state B
# phiNB: survival probability state NB
# phiPB: survival probability state PB
# psiPB: reproduction probability PB
# psiB: reproduction probability B
# psiNB: reproduction probability NB
# pPB: recapture probability B
# pB: recapture probability B
# pNB: recapture probability NB
# -------------------------------------------------
# States (z):
# 1 alive PB
# 2 alive B
# 3 alive NB
# 4 dead
# Observations (y):
# 1 not seen
# 2 seen as PB
# 3 seen as B
# 4 seen as NB
# -------------------------------------------------
# Priors
pPB ~ dunif(0,1)
pB ~ dunif(0,1)
pNB ~ dunif(0,1)
# Transition and Observation matrixes
omega[1,1] <- 1-pNB
omega[1,2] <- 0
omega[1,3] <- 0
omega[1,4] <- pNB
omega[2,1] <- 1-pB
omega[2,2] <- 0
omega[2,3] <- pB
omega[2,4] <- 0
omega[3,1] <- 1-pPB
omega[3,2] <- pPB
omega[3,3] <- 0
omega[3,4] <- 0
omega[4,1] <- 1
omega[4,2] <- 0
omega[4,3] <- 0
omega[4,4] <- 0
for(i in 1:n_ind){
for(t in (first[i]):(n_year-1)){
logit(phiPB[i,t]) <- logit(Mu.phiPB[age[i,t]]) # Here I use logit to set up a model that I could use later to include covariates and random effects
logit(phiB[i,t]) <- logit(Mu.phiB[age[i,t]]) # Also, age is consider as a categorical variable at the moment. Maybe I will try to use a function latter
logit(phiNB[i,t]) <- logit(Mu.phiNB[age[i,t]])
logit(psiPB[i,t]) <- logit(Mu.psiPB[age[i,t]])
logit(psiB[i,t]) <- logit(Mu.psiB[age[i,t]])
logit(psiNB[i,t]) <- logit(Mu.psiNB[age[i,t]])
gamma[1,1,i,t] <- phiPB[i,t] * (1-psiPB[i,t])
gamma[1,2,i,t] <- phiPB[i,t] * psiPB[i,t]
gamma[1,3,i,t] <- 0
gamma[1,4,i,t] <- 1-phiPB[i,t]
gamma[2,1,i,t] <- 0
gamma[2,2,i,t] <- phiB[i,t] * psiB[i,t]
gamma[2,3,i,t] <- phiB[i,t] * (1-psiB[i,t])
gamma[2,4,i,t] <- 1-phiB[i,t]
gamma[3,1,i,t] <- 0
gamma[3,2,i,t] <- phiNB[i,t] * psiNB[i,t]
gamma[3,3,i,t] <- phiNB[i,t] * (1-psiNB[i,t])
gamma[3,4,i,t] <- 1-phiNB[i,t]
gamma[4,1,i,t] <- 0
gamma[4,2,i,t] <- 0
gamma[4,3,i,t] <- 0
gamma[4,4,i,t] <- 1
}
}
for(a in 1:Amax){
Mu.psiPB[a] ~ dunif(0, 1)
Mu.psiNB[a] ~ dunif(0, 1)
Mu.psiB[a] ~ dunif(0, 1)
Mu.phiB[a] ~ dunif(0, 1)
Mu.phiNB[a] ~ dunif(0, 1)
Mu.phiPB[a] ~ dunif(0, 1)
}
# Likelihood
for(i in 1:n_ind){
z[i,first[i]] <- y[i,first[i]] - 1
for(t in (first[i]+1):n_year){
z[i,t] ~ dcat(gamma[z[i,t-1],1:4,i,t-1])
y[i,t] ~ dcat(omega[z[i,t],1:4])
}
}
})
my.data <- list(y = y + 1)
my.constants <- list(first = first,
n_year = ncol(y),
n_ind = nrow(y),
Amax = max(age, na.rm = TRUE),
age = age)
zinits <- y
for (i in 1:ncol(zinits)){
if(first[i] > 1) zinits[i, 1:(first[i]-1)] <- NA
}
zinits[is.na(zinits)] <- 1
zinits[zinits == 0] <- sample(c(2,3), sum(zinits == 0), replace = TRUE)
initial.values <- function(){list(phiPB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
phiB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
phiNB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
psiPB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
psiB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
psiNB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
Mu.psiNB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.psiPB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.psiB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiNB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiPB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
z = zinits)}
parameters.to.save <- c("Mu.psiPB", "Mu.psiNB", "Mu.psiB", "Mu.phiPB", "Mu.phiB", "Mu.phiNB")
# parameters.to.save <- c("psiPB", "psiNB", "psiB", "phiPB", "phiB", "phiNB")
# Here there is a problem, when I estimate these parameters they are not estimated. I have NA for the first value and 0 for all the others.
parameters.to.save
n.iter <- 1000
n.burnin <- 500
n.chains <- 3
mcmc.multistate.cate <- nimbleMCMC(code = multistate.cate,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)
summary_2 <- MCMCsummary(mcmc.multistate.cate, round = 2)
### See chains
priors <- runif(3000, 0, 1)
MCMCtrace(object = mcmc.multistate.cate,
ISB = FALSE,
exact = TRUE,
params = c("Mu.psiPB[1]"),
pdf = FALSE,
priors = priors)
Daniel Gibson
Aug 17, 2023, 12:06:41 PM8/17/23
to Etienne Rouby, nimble-users
Hi Etienne,
I think there are two distinct problems in your current code, both of which are pretty easy to fix. The first issue is with how the initial values for the latent state are being assigned, and the second issue is how the observation matrix (Omega) is defined.
It looks like the initial values for the latent state are allowing for individuals to 1) switch from the breeding/non-breeding state back to the pre-breeding state and 2) go directly from pre-breeding to non-breeding, both of which the model specifies as impossible. This is happening since the current code is randomly assigning a 2 or 3 to each missed observation, but pre-breeders are being missed in the simulation, which assigns an encounter history of 1-0-1-2 with an initial state assignment of 1-2/3-1-2.
Here is a quick solution to resolve this issue with the model you are working with, but you may have to make some changes here when you are using your actual data.
zinits <- matrix(1, ncol = ncol (y), nrow = nrow(y)) # start everyone as a pre-breeder
get.last_juve <- function(x) max(which(x==1)) # define the last time an individual was seen as a pre-breeder
last <- apply(y, 1, get.last_juve)
for(i in 1:nrow(zinits)){
for(j in (last[i] + 1):ncol(zinits)){
zinits[i,j] <- 2 # all occasions past the last known pre-breeding state are assigned as a breeder (since you can't go from pre-breedinging to non-breeding)
}
}
zinits <- ifelse(y == 3, 3, zinits) # go back and change all known non-breeding states to non-breeding.
get.last_juve <- function(x) max(which(x==1)) # define the last time an individual was seen as a pre-breeder
last <- apply(y, 1, get.last_juve)
for(i in 1:nrow(zinits)){
for(j in (last[i] + 1):ncol(zinits)){
zinits[i,j] <- 2 # all occasions past the last known pre-breeding state are assigned as a breeder (since you can't go from pre-breedinging to non-breeding)
}
}
zinits <- ifelse(y == 3, 3, zinits) # go back and change all known non-breeding states to non-breeding.
The second issue is that I think the pre-breeding and non-breeding assignments are indexed incorrectly (i.e., backwards). If you change the omega's to look like this instead -
omega[1,1] <- 1-pPB
omega[1,2] <- pPB
omega[1,3] <- 0
omega[1,4] <- 0
omega[1,2] <- pPB
omega[1,3] <- 0
omega[1,4] <- 0
omega[2,1] <- 1-pB
omega[2,2] <- 0
omega[2,3] <- pB
omega[2,4] <- 0
omega[3,1] <- 1-pNB
omega[3,2] <- 0
omega[3,3] <- 0
omega[3,4] <- pNB
omega[3,2] <- 0
omega[3,3] <- 0
omega[3,4] <- pNB
omega[4,1] <- 1
omega[4,2] <- 0
omega[4,3] <- 0
omega[4,4] <- 0
- the model should work.
Best,
Dan
--
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/327035f1-be88-4c7b-949d-6f46c92453dbn%40googlegroups.com.
Etienne Rouby
Aug 24, 2023, 11:21:14 AM8/24/23
to nimble-users
Dear Daniel,
Thank you very much for your help. It really helped me to understand better the dynamic and the importance of the initial values. Everything worked with this model and I was able to complexify it. Now, I get the same error with the last model I want to build. I looked at the Omega matrix, which seems correct. Also my initial states Z seems relevant (even if the code is still a bit gross). Here is the model code. Do you have any clues for me to try to debug it ? My feeling is that I miss something with the initial states. Since I added 4 new states with different relations between them.
Error I get : warning: logProb of data node y[46, 40]: logProb is -Inf.
Some precisions :
1 is Pre Breeder (PB): Can only go to 2 or 3
n_year <- 40 # number of occasions of samplig
n_cohorts <- 1 # Number of yearly cohorts
# Data simulation --------------------------------------------------------------------------------------------
## Define the transitions Matrixes age-dependants
#
## --- Survival --- ##
#
# phiPB: Survival probability of birds that are in a Pre-Breeding State (ex: Chicks and Juveniles)
# phiSB: Survival probability of birds that are in a Successful breeder State (Breeding Adults that lay an egg)
# phiFB: Survival probability of birds that are in a Failed breeder State (Breeding Adults that do not lay an egg)
# phiNB: Survival probability of birds that are in a Non-Breeding State (Non-Breeding Adults. Only possible after the sabbatical year)
# phiPSB: Survival probability of birds that are in a Post-successful Breeding State (Sabbatical year at sea after a successful breeding)
# phiPFB: Survival probability of birds that are in a Post-failed Breeding State (Sabbatical year at sea after a fail breeding)
#
## --- Breeding --- ##
#
# psiPB: Breeding probability of birds that are in a Pre-Breeding State
# psiSB: Breeding probability of birds that are in a Successful Breeding State
# psiFB: Breeding probability of birds that are in a Failed Breeding State
# psiNB: Breeding probability of birds that are in a Non-Breeding State
# psiPSB: Breeding probability of birds that are in a Post-successful Breeding State
# psiPFB: Breeding probability of birds that are in a Post-failed Breeding State
#
## --- Breeding success --- ##
#
# rhoPB: Breeding success probability of birds that are in a Pre-Breeding State
# rhoSB: Breeding success probability of birds that are in a Successful Breeding State
# rhoFB: Breeding success probability of birds that are in a Failed Breeding State
# rhoNB: Breeding success probability of birds that are in a Non-Breeding State
# rhoPSB: Breeding success probability of birds that are in a Post-successful Breeding State
# rhoPFB: Breeding success probability of birds that are in a Post-failed Breeding State
#
## Simulate matrixes
# Individuals from 1 to 5 years
psiPB <- 0.00
psiSB <- 0.00
psiFB <- 0.00
psiNB <- 0.00
psiPSB <- 0.00
psiPFB <- 0.00
phiPB <- 0.95
phiSB <- 0.00
phiFB <- 0.00
phiNB <- 0.00
phiPSB <- 0.00
phiPFB <- 0.00
rhoPB <- 0.00
rhoSB <- 0.00
rhoFB <- 0.00
rhoNB <- 0.00
rhoPSB <- 0.00
rhoPFB <- 0.00
A_1_5 <- matrix(c(phiPB * (1-psiPB), phiPB * psiPB * rhoPB, phiPB * psiPB * (1-rhoPB), 0, 0, 0, 1-phiPB,
0, phiSB * psiSB * rhoSB, phiSB * psiSB * (1-rhoSB), 0, phiSB * (1-psiSB), 0, 1-phiSB,
0, phiFB * psiFB * rhoFB, phiFB * psiFB * (1-rhoFB), 0, 0, phiFB * (1-psiFB), 1-phiFB,
0, phiNB * psiNB * rhoNB, phiNB * psiNB * (1-rhoNB), phiNB * (1-psiNB), 0, 0, 1-phiNB,
0, phiPSB * psiPSB * rhoPSB, phiPSB * psiPSB * (1-rhoPSB), phiPSB * (1-psiPSB), 0, 0, 1-phiPSB,
0, phiPFB * psiPFB * rhoPFB, phiPFB * psiPFB * (1-rhoPFB), phiPFB * (1-psiPFB), 0, 0, 1-phiPFB,
0, 0, 0, 0, 0, 0, 1), nrow=7, byrow=TRUE)
namstates <- c("PB","SB","FB","NB","PSB","PFB","D")
dimnames(A_1_5) <- list(origin=namstates,destination=namstates)
# Individuals from 6 to 8 years
psiPB <- 0.80
psiSB <- 0.80
psiFB <- 0.80
psiNB <- 0.80
psiPSB <- 0.80
psiPFB <- 0.80
phiPB <- 0.95
phiSB <- 0.90
phiFB <- 0.90
phiNB <- 0.90
phiPSB <- 0.90
phiPFB <- 0.90
rhoPB <- 0.50
rhoSB <- 0.50
rhoFB <- 0.50
rhoNB <- 0.50
rhoPSB <- 0.50
rhoPFB <- 0.50
A_6_8 <- matrix(c(phiPB * (1-psiPB), phiPB * psiPB * rhoPB, phiPB * psiPB * (1-rhoPB), 0, 0, 0, 1-phiPB,
0, phiSB * psiSB * rhoSB, phiSB * psiSB * (1-rhoSB), 0, phiSB * (1-psiSB), 0, 1-phiSB,
0, phiFB * psiFB * rhoFB, phiFB * psiFB * (1-rhoFB), 0, 0, phiFB * (1-psiFB), 1-phiFB,
0, phiNB * psiNB * rhoNB, phiNB * psiNB * (1-rhoNB), phiNB * (1-psiNB), 0, 0, 1-phiNB,
0, phiPSB * psiPSB * rhoPSB, phiPSB * psiPSB * (1-rhoPSB), phiPSB * (1-psiPSB), 0, 0, 1-phiPSB,
0, phiPFB * psiPFB * rhoPFB, phiPFB * psiPFB * (1-rhoPFB), phiPFB * (1-psiPFB), 0, 0, 1-phiPFB,
0, 0, 0, 0, 0, 0, 1), nrow=7, byrow=TRUE)
namstates <- c("PB","SB","FB","NB","PSB","PFB","D")
dimnames(A_6_8) <- list(origin=namstates,destination=namstates)
# Individuals from 9 to 11 years
psiPB <- 0.80
psiSB <- 0.80
psiFB <- 0.80
psiNB <- 0.80
psiPSB <- 0.80
psiPFB <- 0.80
phiPB <- 0.90
phiSB <- 0.90
phiFB <- 0.95
phiNB <- 0.95
phiPSB <- 0.95
phiPFB <- 0.95
rhoPB <- 0.80
rhoSB <- 0.80
rhoFB <- 0.80
rhoNB <- 0.80
rhoPSB <- 0.80
rhoPFB <- 0.80
A_9_11 <- matrix(c(phiPB * (1-psiPB), phiPB * psiPB * rhoPB, phiPB * psiPB * (1-rhoPB), 0, 0, 0, 1-phiPB,
0, phiSB * psiSB * rhoSB, phiSB * psiSB * (1-rhoSB), 0, phiSB * (1-psiSB), 0, 1-phiSB,
0, phiFB * psiFB * rhoFB, phiFB * psiFB * (1-rhoFB), 0, 0, phiFB * (1-psiFB), 1-phiFB,
0, phiNB * psiNB * rhoNB, phiNB * psiNB * (1-rhoNB), phiNB * (1-psiNB), 0, 0, 1-phiNB,
0, phiPSB * psiPSB * rhoPSB, phiPSB * psiPSB * (1-rhoPSB), phiPSB * (1-psiPSB), 0, 0, 1-phiPSB,
0, phiPFB * psiPFB * rhoPFB, phiPFB * psiPFB * (1-rhoPFB), phiPFB * (1-psiPFB), 0, 0, 1-phiPFB,
0, 0, 0, 0, 0, 0, 1), nrow=7, byrow=TRUE)
namstates <- c("PB","SB","FB","NB","PSB","PFB","D")
dimnames(A_9_11) <- list(origin=namstates,destination=namstates)
# Individuals from 12 to 12+ years
psiPB <- 0.80
psiSB <- 0.80
psiFB <- 0.80
psiNB <- 0.80
psiPSB <- 0.80
psiPFB <- 0.80
phiPB <- 0.85
phiSB <- 0.90
phiFB <- 0.95
phiNB <- 0.95
phiPSB <- 0.95
phiPFB <- 0.95
rhoPB <- 0.80
rhoSB <- 0.80
rhoFB <- 0.80
rhoNB <- 0.80
rhoPSB <- 0.80
rhoPFB <- 0.80
A_12 <- matrix(c(phiPB * (1-psiPB), phiPB * psiPB * rhoPB, phiPB * psiPB * (1-rhoPB), 0, 0, 0, 1-phiPB,
0, phiSB * psiSB * rhoSB, phiSB * psiSB * (1-rhoSB), 0, phiSB * (1-psiSB), 0, 1-phiSB,
0, phiFB * psiFB * rhoFB, phiFB * psiFB * (1-rhoFB), 0, 0, phiFB * (1-psiFB), 1-phiFB,
0, phiNB * psiNB * rhoNB, phiNB * psiNB * (1-rhoNB), phiNB * (1-psiNB), 0, 0, 1-phiNB,
0, phiPSB * psiPSB * rhoPSB, phiPSB * psiPSB * (1-rhoPSB), phiPSB * (1-psiPSB), 0, 0, 1-phiPSB,
0, phiPFB * psiPFB * rhoPFB, phiPFB * psiPFB * (1-rhoPFB), phiPFB * (1-psiPFB), 0, 0, 1-phiPFB,
0, 0, 0, 0, 0, 0, 1), nrow=7, byrow=TRUE)
namstates <- c("PB","SB","FB","NB","PSB","PFB","D")
dimnames(A_12) <- list(origin=namstates,destination=namstates)
# first <- c(rep(1, n_ind), # First sampling occasion. Here assume a unique cohort. After, I will have to assume several
# rep(2, n_ind),
# rep(3, n_ind))
first <- c(rep(1, n_ind))
real_age <- matrix(NA, nrow = n_ind * n_cohorts, ncol = n_year)
for (i in 1:nrow(real_age)){
for (j in 1:ncol(real_age)){
if (j == first[i]) real_age[i,j] <- 1
if (j == first[i]+1) real_age[i,j] <- 2
if (j == first[i]+2) real_age[i,j] <- 3
if (j == first[i]+3) real_age[i,j] <- 4
if (j == first[i]+4) real_age[i,j] <- 5
if (j == first[i]+5) real_age[i,j] <- 6
if (j == first[i]+6) real_age[i,j] <- 7
if (j == first[i]+7) real_age[i,j] <- 8
if (j == first[i]+8) real_age[i,j] <- 9
if (j == first[i]+9) real_age[i,j] <- 10
if (j == first[i]+10) real_age[i,j] <- 11
if (j >= first[i]+11) real_age[i,j] <- 12 # I stop at 12 because from 12 to Death, I consider the individuals belong to the same real_age class ( >= 12)
}
}
# Simulate life trajectories for each individuals within the cohort
# Individuals starts as Pre Breeder. They become Successful or Failed Breeder with a certain probability.
# Then, they transit from sucessful to post sucessful or from failed to post failed (it is the sabbatical year).
# From there, they can breed again (either successful or failed) or become a non breeder.
# From the non breeder stage, they can only become Successful or failed breeder again.
#
z <- matrix(NA, nrow = n_ind * n_cohorts, ncol = n_year)
for (i in 1:(n_ind * n_cohorts)){ # loop over individuals
z[i,first[i]] <- 1 # all individuals are alive at first occasion
for (t in (first[i]+1):n_year){ # loop over time
if(is.na(real_age[i,t])) {
state_draw <- NA
}
if(real_age[i,t] <= 5) {
state_draw <- as.data.frame(rmultinom(1, 1, A_1_5[z[i,t-1], 1:7])) %>% # draw the state at t given the state at t-1 and transition matrix
tibble::rownames_to_column("state")
}
if(real_age[i,t] >= 6 & real_age[i,t] <= 8) {
state_draw <- as.data.frame(rmultinom(1, 1, A_6_8[z[i,t-1], 1:7])) %>% # draw the state at t given the state at t-1 and transition matrix
tibble::rownames_to_column("state")
}
if(real_age[i,t] >= 9 & real_age[i,t] <= 11) {
state_draw <- as.data.frame(rmultinom(1, 1, A_9_11[z[i,t-1], 1:7])) %>% # draw the state at t given the state at t-1 and transition matrix
tibble::rownames_to_column("state")
}
if(real_age[i,t] >= 12) {
state_draw <- as.data.frame(rmultinom(1, 1, A_12[z[i,t-1], 1:7])) %>% # draw the state at t given the state at t-1 and transition matrix
tibble::rownames_to_column("state")
}
if (filter(state_draw, V1 == 1)$state == "PB") { # Assign the good index depending on the state
z[i,t] <- 1
}
if (filter(state_draw, V1 == 1)$state == "SB") {
z[i,t] <- 2
}
if (filter(state_draw, V1 == 1)$state == "FB") {
z[i,t] <- 3
}
if (filter(state_draw, V1 == 1)$state == "NB") {
z[i,t] <- 4
}
if (filter(state_draw, V1 == 1)$state == "PSB") {
z[i,t] <- 5
}
if (filter(state_draw, V1 == 1)$state == "PFB") {
z[i,t] <- 6
}
if (filter(state_draw, V1 == 1)$state == "D") {
z[i,t] <- 7
}
if (is.na(filter(state_draw, V1 == 1)$state)) { # Assign the good index depending on the state
z[i,t] <- NA
}
}
}
# Simulate encounter histories for each individuals within the cohort
y <- matrix(NA, nrow = n_ind * n_cohorts, ncol = n_year)
p <- 0.8 # Detection probability for the states that are not in the sabbatical year
patsea <- 0 # Detection probability for PSB and PFB. It equals 0 because the animals are at sea and cannot be monitored
# Observation matrix
B <- matrix(c(p, 0, 0, 0, 0, 0, 1-p,
0, p, 0, 0, 0, 0, 1-p,
0, 0, p, 0, 0, 0, 1-p,
0, 0, 0, p, 0, 0, 1-p,
0, 0, 0, 0, patsea, 0, 1-patsea,
0, 0, 0, 0, 0, patsea, 1-patsea,
0, 0, 0, 0, 0, 0, 1), nrow=7, byrow=TRUE)
namstates <- c("PB","SB","FB","NB","PSB","PFB","ND") # ND is for "Non-Detected"
dimnames(B) <- list(origin=namstates,destination=namstates)
for (i in 1:(n_ind * n_cohorts)){ # loop over individuals
y[i,first[i]] <- 1 # all individuals are alive at first occasion
for (t in (first[i]+1):n_year){ # loop over time
obs_draw <- as.data.frame(rmultinom(1, 1, B[z[i,t], 1:7])) %>% # The observation state at t depends on the living state (PB, B, NB or Dead) at t.
tibble::rownames_to_column("obs_state")
if (filter(obs_draw, V1 == 1)$obs_state == "PB") {
y[i,t] <- 1
}
if (filter(obs_draw, V1 == 1)$obs_state == "SB") {
y[i,t] <- 2
}
if (filter(obs_draw, V1 == 1)$obs_state == "FB") {
y[i,t] <- 3
}
if (filter(obs_draw, V1 == 1)$obs_state == "NB") { # ND is for Non-Detected
y[i,t] <- 4
}
if (filter(obs_draw, V1 == 1)$obs_state == "PSB") { # ND is for Non-Detected
y[i,t] <- 5
}
if (filter(obs_draw, V1 == 1)$obs_state == "PFB") { # ND is for Non-Detected
y[i,t] <- 6
}
if (filter(obs_draw, V1 == 1)$obs_state == "ND") { # ND is for Non-Detected
y[i,t] <- 7
}
if (is.na(filter(obs_draw, V1 == 1)$obs_state == "ND")) {
y[i,t] <- NA
}
}
}
y[y == 7] <- 0 # Correct the for in order to have 0s, because dead individuals ar not observed.
# Obtain the age of each monitored individual if they are alive each year.
age <- matrix(NA, nrow = nrow(y), ncol = ncol(y)-1)
for (i in 1:nrow(age)){
for (j in 1:ncol(age)){
if (j == first[i]) age[i,j] <- 1
if (j == first[i]+1) age[i,j] <- 2
if (j == first[i]+2) age[i,j] <- 3
if (j == first[i]+3) age[i,j] <- 4
if (j == first[i]+4) age[i,j] <- 5
if (j == first[i]+5) age[i,j] <- 6
if (j == first[i]+6) age[i,j] <- 7
if (j == first[i]+7) age[i,j] <- 8
if (j == first[i]+8) age[i,j] <- 9
if (j == first[i]+9) age[i,j] <- 10
if (j == first[i]+10) age[i,j] <- 11
if (j >= first[i]+11) age[i,j] <- 12 # I stop at 12 because from 12 to Death, I consider the individuals belong to the same age class ( >= 12)
}
}
# Parameter estimation from simulated data ---------------------------------------------------------------------------------------------------------
# Nimble model
multistate.cate <- nimbleCode({
# -------------------------------------------------
# Parameters:
# phiPB: survival probability state Pre-Breeder
# phiSB: survival probability state Successful Breeder
# phiFB: survival probability state Failed-Breeder
# phiNB: survival probability state Non-Breeder
# phiPSB: survival probability state Post-successful Breeder
# phiPFB: survival probability state Post-Fail-Breeder
#
# psiPB: breeding probability state Pre-Breeder
# psiSB: breeding probability state Successful Breeder
# psiFB: breeding probability state Failed-Breeder
# psiNB: breeding probability state Non-Breeder
# psiPSB: breeding probability state Post-successful Breeder
# psiPFB: breeding probability state Post-Fail-Breeder
#
# rhoPB: breeding success probability state Pre-Breeder
# rhoSB: breeding success probability state Successful Breeder
# rhoFB: breeding success probability state Failed-Breeder
# rhoNB: breeding success probability state Non-Breeder
# rhoPSB: breeding success probability state Post-successful Breeder
# rhoPFB: breeding success probability state Post-Fail-Breeder
# -------------------------------------------------
# States (z):
# 1 alive PB
# 2 alive SB
# 3 alive FB
# 4 alive NB
# 5 alive PSB
# 6 alive PFB
# 7 dead
# Observations (y):
# 1 not seen
# 2 seen as PB
# 3 seen as SB
# 4 seen as FB
# 5 seen as NB
# 6 seen as PSB
# 7 seen as PFB
# -------------------------------------------------
# Priors
pPB ~ dunif(0,1)
pSB ~ dunif(0,1)
pFB ~ dunif(0,1)
pNB ~ dunif(0,1)
pPSB ~ dunif(0,1)
pPFB ~ dunif(0,1)
# Transition and Observation matrixes
omega[1,6] <- 0
omega[1,7] <- 0
omega[2,1] <- 1-pSB
omega[2,2] <- 0
omega[2,3] <- pSB
omega[2,4] <- 0
omega[2,5] <- 0
omega[2,6] <- 0
omega[2,7] <- 0
omega[3,1] <- 1-pFB
omega[3,2] <- 0
omega[3,3] <- 0
omega[3,4] <- pFB
omega[3,5] <- 0
omega[3,6] <- 0
omega[3,7] <- 0
omega[4,1] <- 1-pNB
omega[4,2] <- 0
omega[4,3] <- 0
omega[4,4] <- 0
omega[4,5] <- pNB
omega[4,6] <- 0
omega[4,7] <- 0
omega[5,1] <- 1-pPSB
omega[5,2] <- 0
omega[5,3] <- 0
omega[5,4] <- 0
omega[5,5] <- 0
omega[5,6] <- pPSB
omega[5,7] <- 0
omega[6,1] <- 1-pPFB
omega[6,2] <- 0
omega[6,3] <- 0
omega[6,4] <- 0
omega[6,5] <- 0
omega[6,6] <- 0
omega[6,7] <- pPFB
omega[7,1] <- 1
omega[7,2] <- 0
omega[7,3] <- 0
omega[7,4] <- 0
omega[7,5] <- 0
omega[7,6] <- 0
omega[7,7] <- 0
for(i in 1:n_ind){
for(t in (first[i]):(n_year-1)){
logit(phiPB[i,t]) <- logit(Mu.phiPB[age[i,t]]) # Here I use logit to set up a model that I could use later to include covariates and random effects
logit(phiSB[i,t]) <- logit(Mu.phiSB[age[i,t]]) # Also, age is consider as a categorical variable at the moment. Maybe I will try to use a function latter
logit(phiFB[i,t]) <- logit(Mu.phiFB[age[i,t]])
logit(phiNB[i,t]) <- logit(Mu.phiNB[age[i,t]])
logit(phiPSB[i,t]) <- logit(Mu.phiPSB[age[i,t]])
logit(phiPFB[i,t]) <- logit(Mu.phiPFB[age[i,t]])
logit(psiPB[i,t]) <- logit(Mu.psiPB[age[i,t]])
logit(psiSB[i,t]) <- logit(Mu.psiSB[age[i,t]])
logit(psiFB[i,t]) <- logit(Mu.psiFB[age[i,t]])
logit(psiNB[i,t]) <- logit(Mu.psiNB[age[i,t]])
logit(psiPSB[i,t]) <- logit(Mu.psiPSB[age[i,t]])
logit(psiPFB[i,t]) <- logit(Mu.psiPFB[age[i,t]])
logit(rhoPB[i,t]) <- logit(Mu.rhoPB[age[i,t]])
logit(rhoSB[i,t]) <- logit(Mu.rhoSB[age[i,t]])
logit(rhoFB[i,t]) <- logit(Mu.rhoFB[age[i,t]])
logit(rhoNB[i,t]) <- logit(Mu.rhoNB[age[i,t]])
logit(rhoPSB[i,t]) <- logit(Mu.rhoPSB[age[i,t]])
logit(rhoPFB[i,t]) <- logit(Mu.rhoPFB[age[i,t]])
gamma[1,1,i,t] <- phiPB[i,t] * (1-psiPB[i,t])
gamma[1,2,i,t] <- phiPB[i,t] * psiPB[i,t] * rhoPB[i,t]
gamma[1,3,i,t] <- phiPB[i,t] * psiPB[i,t] * (1-rhoPB[i,t])
gamma[1,4,i,t] <- 0
gamma[1,5,i,t] <- 0
gamma[1,6,i,t] <- 0
gamma[1,7,i,t] <- 1-phiPB[i,t]
gamma[2,1,i,t] <- 0
gamma[2,2,i,t] <- phiSB[i,t] * psiSB[i,t] * rhoSB[i,t]
gamma[2,3,i,t] <- phiSB[i,t] * psiSB[i,t] * (1-rhoSB[i,t])
gamma[2,4,i,t] <- 0
gamma[2,5,i,t] <- phiSB[i,t] * (1-psiSB[i,t])
gamma[2,6,i,t] <- 0
gamma[2,7,i,t] <- 1-phiSB[i,t]
gamma[3,1,i,t] <- 0
gamma[3,2,i,t] <- phiFB[i,t] * psiFB[i,t] * rhoFB[i,t]
gamma[3,3,i,t] <- phiFB[i,t] * psiFB[i,t] * (1-rhoFB[i,t])
gamma[3,4,i,t] <- 0
gamma[3,5,i,t] <- 0
gamma[3,6,i,t] <- phiFB[i,t] * (1-psiFB[i,t])
gamma[3,7,i,t] <- 1-phiFB[i,t]
gamma[4,1,i,t] <- 0
gamma[4,2,i,t] <- phiNB[i,t] * psiNB[i,t] * rhoNB[i,t]
gamma[4,3,i,t] <- phiNB[i,t] * psiNB[i,t] * (1-rhoNB[i,t])
gamma[4,4,i,t] <- phiNB[i,t] * (1-psiNB[i,t])
gamma[4,5,i,t] <- 0
gamma[4,6,i,t] <- 0
gamma[4,7,i,t] <- 1-phiNB[i,t]
gamma[5,1,i,t] <- 0
gamma[5,2,i,t] <- phiPSB[i,t] * psiPSB[i,t] * rhoPSB[i,t]
gamma[5,3,i,t] <- phiPSB[i,t] * psiPSB[i,t] * (1-rhoPSB[i,t])
gamma[5,4,i,t] <- phiPSB[i,t] * (1-psiPSB[i,t])
gamma[5,5,i,t] <- 0
gamma[5,6,i,t] <- 0
gamma[5,7,i,t] <- 1-phiPSB[i,t]
gamma[6,1,i,t] <- 0
gamma[6,2,i,t] <- phiPFB[i,t] * psiPFB[i,t] * rhoPFB[i,t]
gamma[6,3,i,t] <- phiPFB[i,t] * psiPFB[i,t] * (1-rhoPFB[i,t])
gamma[6,4,i,t] <- phiPFB[i,t] * (1-psiPFB[i,t])
gamma[6,5,i,t] <- 0
gamma[6,6,i,t] <- 0
gamma[6,7,i,t] <- 1-phiPFB[i,t]
gamma[7,1,i,t] <- 0
gamma[7,2,i,t] <- 0
gamma[7,3,i,t] <- 0
gamma[7,4,i,t] <- 0
gamma[7,5,i,t] <- 0
gamma[7,6,i,t] <- 0
gamma[7,7,i,t] <- 1
}
}
for(a in 1:Amax){
Mu.psiPB[a] ~ dunif(0, 1)
Mu.psiSB[a] ~ dunif(0, 1)
Mu.psiFB[a] ~ dunif(0, 1)
Mu.psiNB[a] ~ dunif(0, 1)
Mu.psiPSB[a] ~ dunif(0, 1)
Mu.psiPFB[a] ~ dunif(0, 1)
Mu.phiPB[a] ~ dunif(0, 1)
Mu.phiSB[a] ~ dunif(0, 1)
Mu.phiFB[a] ~ dunif(0, 1)
Mu.phiNB[a] ~ dunif(0, 1)
Mu.phiPSB[a] ~ dunif(0, 1)
Mu.phiPFB[a] ~ dunif(0, 1)
Mu.rhoPB[a] ~ dunif(0, 1)
Mu.rhoSB[a] ~ dunif(0, 1)
Mu.rhoFB[a] ~ dunif(0, 1)
Mu.rhoNB[a] ~ dunif(0, 1)
Mu.rhoPSB[a] ~ dunif(0, 1)
Mu.rhoPFB[a] ~ dunif(0, 1)
}
# Likelihood
for(i in 1:n_ind){
z[i,first[i]] <- y[i,first[i]] - 1
for(t in (first[i]+1):n_year){
z[i,t] ~ dcat(gamma[z[i,t-1],1:7,i,t-1])
y[i,t] ~ dcat(omega[z[i,t],1:7])
}
}
})
my.data <- list(y = y + 1)
my.constants <- list(first = first,
n_year = ncol(y),
n_ind = nrow(y),
Amax = max(age, na.rm = TRUE),
age = age)
zinits <- matrix(1, ncol = ncol(y), nrow = nrow(y)) # start everyone as a pre-breeder
get.last_juve <- function(x) max(which(x==1)) # define the last time an individual was seen as a pre-breeder
last <- apply(y, 1, get.last_juve)
for(i in 1:nrow(zinits)){
if((last[i] + 1) <= ncol(zinits)) {
for(j in (last[i] + 1):ncol(zinits)){
zinits[i,j] <- sample(2:3, 1) # all occasions past the last known pre-breeding state are assigned as a sucessful-breeder or failed breeder (since you can't go from pre-breeding to non-breeding, post fail or post sucess)
}
}
}
zinits <- ifelse(z == 5, 5, zinits) # go back and change all known post succesful breeding states.
zinits <- ifelse(z == 6, 6, zinits) # go back and change all known post failed breeding states.
for(i in 1:nrow(zinits)) { # then from here, correct the values before the states 5 and 6. Before a 5, there only can be a 2. Before a 6 there only can be a 3.
for(j in 2:ncol(zinits)) {
if(zinits[i,j] == 6 & zinits[i,j-1] == 2) {
zinits[i,j-1] <- 3
}
if(zinits[i,j] == 5 & zinits[i,j-1] == 3) {
zinits[i,j-1] <- 2
}
}
}
zinits <- ifelse(z == 4, 4, zinits) # go back and include all known Non-Breeding states
for(i in 1:nrow(zinits)) { # from here, change the values before the 4. Before a 4 there only can be a 5 or 6. So if no 5 or 6, 4 become a 2 if the next value is a 5, and a 3 if the next value is a 6. If not, 2 or 3 ar randomly drawed.
for(j in 2:ncol(zinits)) {
if((j+1) <= ncol(zinits)) {
if(zinits[i,j] == 4) {
if(zinits[i,j-1] != 5 & zinits[i,j-1] != 6 & zinits[i,j+1] == 5) { #)
zinits[i,j] <- 2
}
if(zinits[i,j-1] != 5 & zinits[i,j-1] != 6 & zinits[i,j+1] == 6) { #)
zinits[i,j] <- 3
}
if(zinits[i,j-1] != 5 & zinits[i,j-1] != 6 & zinits[i,j+1] != 6 & zinits[i,j+1] != 5) { #)
zinits[i,j] <- sample(2:3,1)
}
}
}
}
}
initial.values <- function(){list(phiPB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
phiSB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
phiFB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
phiNB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
phiPSB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
phiPFB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
psiPB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
psiSB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
psiFB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
psiNB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
psiPSB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
psiPFB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
rhoPB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
rhoSB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
rhoFB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
rhoNB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
rhoPSB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
rhoPFB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
Mu.psiFB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.psiNB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.psiPSB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.psiPFB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiPB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiSB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiFB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiNB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiPSB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiPFB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.rhoPB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.rhoSB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.rhoFB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.rhoNB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.rhoPSB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.rhoPFB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
z = zinits)}
parameters.to.save <- c("Mu.psiPB", "Mu.psiSB", "Mu.psiFB", "Mu.psiNB", "Mu.psiPSB", "Mu.psiPFB",
"Mu.phiPB", "Mu.phiSB", "Mu.phiFB", "Mu.phiNB", "Mu.phiPSB", "Mu.phiPFB",
"Mu.rhoPB", "Mu.rhoSB", "Mu.rhoFB", "Mu.rhoNB", "Mu.rhoPSB", "Mu.rhoPFB")
parameters.to.save
n.iter <- 1000
n.burnin <- 500
n.chains <- 3
mcmc.multistate.cate <- nimbleMCMC(code = multistate.cate,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)
Thank you very much for your help. It really helped me to understand better the dynamic and the importance of the initial values. Everything worked with this model and I was able to complexify it. Now, I get the same error with the last model I want to build. I looked at the Omega matrix, which seems correct. Also my initial states Z seems relevant (even if the code is still a bit gross). Here is the model code. Do you have any clues for me to try to debug it ? My feeling is that I miss something with the initial states. Since I added 4 new states with different relations between them.
Error I get : warning: logProb of data node y[46, 40]: logProb is -Inf.
Some precisions :
1 is Pre Breeder (PB): Can only go to 2 or 3
2 is Successful Breeder (SB): Can only go to 2,3 or 5
3 is Failed Breeder (FB): Can only go to 2,3 or 6
4 is Non Breeder (NB): Can only go to 2 or 3
5 is Post successful Breeder (PSB): Can only go to 2,3 or 4
5 is Post Failed breeder (PFB): Can only go to 2,3 or 4
I tried to represent this dynamic in my Zinitis. But I think I miss something. A 0 is taken into account somewhere, that is why I get a logProb_y equals to Inf in my opinion.
I really apologize for this new message since it may not be a pure debugging question. If nobody has time for that, I understand perfectly. But If anybody has even a clue I will take it.
Best regards,
Etienne
I really apologize for this new message since it may not be a pure debugging question. If nobody has time for that, I understand perfectly. But If anybody has even a clue I will take it.
Best regards,
Etienne
######################
######################
# Remove all objects from the environment
rm(list = ls())
# Load the relevant packages
lapply(c("tidyverse", "dplyr", "tidyr", "nimble", "MCMCvis"), library, character.only = TRUE)
## Set parameter values to generate fake data
n_ind <- 100 # number of individuals to generate per cohort######################
# Remove all objects from the environment
rm(list = ls())
# Load the relevant packages
lapply(c("tidyverse", "dplyr", "tidyr", "nimble", "MCMCvis"), library, character.only = TRUE)
## Set parameter values to generate fake data
n_year <- 40 # number of occasions of samplig
# Data simulation --------------------------------------------------------------------------------------------
## Define the transitions Matrixes age-dependants
#
## --- Survival --- ##
#
# phiPB: Survival probability of birds that are in a Pre-Breeding State (ex: Chicks and Juveniles)
# phiSB: Survival probability of birds that are in a Successful breeder State (Breeding Adults that lay an egg)
# phiFB: Survival probability of birds that are in a Failed breeder State (Breeding Adults that do not lay an egg)
# phiNB: Survival probability of birds that are in a Non-Breeding State (Non-Breeding Adults. Only possible after the sabbatical year)
# phiPSB: Survival probability of birds that are in a Post-successful Breeding State (Sabbatical year at sea after a successful breeding)
# phiPFB: Survival probability of birds that are in a Post-failed Breeding State (Sabbatical year at sea after a fail breeding)
#
## --- Breeding --- ##
#
# psiPB: Breeding probability of birds that are in a Pre-Breeding State
# psiSB: Breeding probability of birds that are in a Successful Breeding State
# psiFB: Breeding probability of birds that are in a Failed Breeding State
# psiNB: Breeding probability of birds that are in a Non-Breeding State
# psiPSB: Breeding probability of birds that are in a Post-successful Breeding State
# psiPFB: Breeding probability of birds that are in a Post-failed Breeding State
#
## --- Breeding success --- ##
#
# rhoPB: Breeding success probability of birds that are in a Pre-Breeding State
# rhoSB: Breeding success probability of birds that are in a Successful Breeding State
# rhoFB: Breeding success probability of birds that are in a Failed Breeding State
# rhoNB: Breeding success probability of birds that are in a Non-Breeding State
# rhoPSB: Breeding success probability of birds that are in a Post-successful Breeding State
# rhoPFB: Breeding success probability of birds that are in a Post-failed Breeding State
#
## Simulate matrixes
# Individuals from 1 to 5 years
psiPB <- 0.00
psiSB <- 0.00
psiFB <- 0.00
psiNB <- 0.00
psiPSB <- 0.00
psiPFB <- 0.00
phiPB <- 0.95
phiSB <- 0.00
phiFB <- 0.00
phiNB <- 0.00
phiPSB <- 0.00
phiPFB <- 0.00
rhoPB <- 0.00
rhoSB <- 0.00
rhoFB <- 0.00
rhoNB <- 0.00
rhoPSB <- 0.00
rhoPFB <- 0.00
A_1_5 <- matrix(c(phiPB * (1-psiPB), phiPB * psiPB * rhoPB, phiPB * psiPB * (1-rhoPB), 0, 0, 0, 1-phiPB,
0, phiSB * psiSB * rhoSB, phiSB * psiSB * (1-rhoSB), 0, phiSB * (1-psiSB), 0, 1-phiSB,
0, phiFB * psiFB * rhoFB, phiFB * psiFB * (1-rhoFB), 0, 0, phiFB * (1-psiFB), 1-phiFB,
0, phiNB * psiNB * rhoNB, phiNB * psiNB * (1-rhoNB), phiNB * (1-psiNB), 0, 0, 1-phiNB,
0, phiPSB * psiPSB * rhoPSB, phiPSB * psiPSB * (1-rhoPSB), phiPSB * (1-psiPSB), 0, 0, 1-phiPSB,
0, phiPFB * psiPFB * rhoPFB, phiPFB * psiPFB * (1-rhoPFB), phiPFB * (1-psiPFB), 0, 0, 1-phiPFB,
0, 0, 0, 0, 0, 0, 1), nrow=7, byrow=TRUE)
namstates <- c("PB","SB","FB","NB","PSB","PFB","D")
dimnames(A_1_5) <- list(origin=namstates,destination=namstates)
# Individuals from 6 to 8 years
psiPB <- 0.80
psiSB <- 0.80
psiFB <- 0.80
psiNB <- 0.80
psiPSB <- 0.80
psiPFB <- 0.80
phiPB <- 0.95
phiSB <- 0.90
phiFB <- 0.90
phiNB <- 0.90
phiPSB <- 0.90
phiPFB <- 0.90
rhoPB <- 0.50
rhoSB <- 0.50
rhoFB <- 0.50
rhoNB <- 0.50
rhoPSB <- 0.50
rhoPFB <- 0.50
A_6_8 <- matrix(c(phiPB * (1-psiPB), phiPB * psiPB * rhoPB, phiPB * psiPB * (1-rhoPB), 0, 0, 0, 1-phiPB,
0, phiSB * psiSB * rhoSB, phiSB * psiSB * (1-rhoSB), 0, phiSB * (1-psiSB), 0, 1-phiSB,
0, phiFB * psiFB * rhoFB, phiFB * psiFB * (1-rhoFB), 0, 0, phiFB * (1-psiFB), 1-phiFB,
0, phiNB * psiNB * rhoNB, phiNB * psiNB * (1-rhoNB), phiNB * (1-psiNB), 0, 0, 1-phiNB,
0, phiPSB * psiPSB * rhoPSB, phiPSB * psiPSB * (1-rhoPSB), phiPSB * (1-psiPSB), 0, 0, 1-phiPSB,
0, phiPFB * psiPFB * rhoPFB, phiPFB * psiPFB * (1-rhoPFB), phiPFB * (1-psiPFB), 0, 0, 1-phiPFB,
0, 0, 0, 0, 0, 0, 1), nrow=7, byrow=TRUE)
namstates <- c("PB","SB","FB","NB","PSB","PFB","D")
dimnames(A_6_8) <- list(origin=namstates,destination=namstates)
# Individuals from 9 to 11 years
psiPB <- 0.80
psiSB <- 0.80
psiFB <- 0.80
psiNB <- 0.80
psiPSB <- 0.80
psiPFB <- 0.80
phiPB <- 0.90
phiSB <- 0.90
phiFB <- 0.95
phiNB <- 0.95
phiPSB <- 0.95
phiPFB <- 0.95
rhoPB <- 0.80
rhoSB <- 0.80
rhoFB <- 0.80
rhoNB <- 0.80
rhoPSB <- 0.80
rhoPFB <- 0.80
A_9_11 <- matrix(c(phiPB * (1-psiPB), phiPB * psiPB * rhoPB, phiPB * psiPB * (1-rhoPB), 0, 0, 0, 1-phiPB,
0, phiSB * psiSB * rhoSB, phiSB * psiSB * (1-rhoSB), 0, phiSB * (1-psiSB), 0, 1-phiSB,
0, phiFB * psiFB * rhoFB, phiFB * psiFB * (1-rhoFB), 0, 0, phiFB * (1-psiFB), 1-phiFB,
0, phiNB * psiNB * rhoNB, phiNB * psiNB * (1-rhoNB), phiNB * (1-psiNB), 0, 0, 1-phiNB,
0, phiPSB * psiPSB * rhoPSB, phiPSB * psiPSB * (1-rhoPSB), phiPSB * (1-psiPSB), 0, 0, 1-phiPSB,
0, phiPFB * psiPFB * rhoPFB, phiPFB * psiPFB * (1-rhoPFB), phiPFB * (1-psiPFB), 0, 0, 1-phiPFB,
0, 0, 0, 0, 0, 0, 1), nrow=7, byrow=TRUE)
namstates <- c("PB","SB","FB","NB","PSB","PFB","D")
dimnames(A_9_11) <- list(origin=namstates,destination=namstates)
# Individuals from 12 to 12+ years
psiPB <- 0.80
psiSB <- 0.80
psiFB <- 0.80
psiNB <- 0.80
psiPSB <- 0.80
psiPFB <- 0.80
phiPB <- 0.85
phiSB <- 0.90
phiFB <- 0.95
phiNB <- 0.95
phiPSB <- 0.95
phiPFB <- 0.95
rhoPB <- 0.80
rhoSB <- 0.80
rhoFB <- 0.80
rhoNB <- 0.80
rhoPSB <- 0.80
rhoPFB <- 0.80
A_12 <- matrix(c(phiPB * (1-psiPB), phiPB * psiPB * rhoPB, phiPB * psiPB * (1-rhoPB), 0, 0, 0, 1-phiPB,
0, phiSB * psiSB * rhoSB, phiSB * psiSB * (1-rhoSB), 0, phiSB * (1-psiSB), 0, 1-phiSB,
0, phiFB * psiFB * rhoFB, phiFB * psiFB * (1-rhoFB), 0, 0, phiFB * (1-psiFB), 1-phiFB,
0, phiNB * psiNB * rhoNB, phiNB * psiNB * (1-rhoNB), phiNB * (1-psiNB), 0, 0, 1-phiNB,
0, phiPSB * psiPSB * rhoPSB, phiPSB * psiPSB * (1-rhoPSB), phiPSB * (1-psiPSB), 0, 0, 1-phiPSB,
0, phiPFB * psiPFB * rhoPFB, phiPFB * psiPFB * (1-rhoPFB), phiPFB * (1-psiPFB), 0, 0, 1-phiPFB,
0, 0, 0, 0, 0, 0, 1), nrow=7, byrow=TRUE)
namstates <- c("PB","SB","FB","NB","PSB","PFB","D")
dimnames(A_12) <- list(origin=namstates,destination=namstates)
# first <- c(rep(1, n_ind), # First sampling occasion. Here assume a unique cohort. After, I will have to assume several
# rep(2, n_ind),
# rep(3, n_ind))
first <- c(rep(1, n_ind))
real_age <- matrix(NA, nrow = n_ind * n_cohorts, ncol = n_year)
for (i in 1:nrow(real_age)){
for (j in 1:ncol(real_age)){
if (j == first[i]) real_age[i,j] <- 1
if (j == first[i]+1) real_age[i,j] <- 2
if (j == first[i]+2) real_age[i,j] <- 3
if (j == first[i]+3) real_age[i,j] <- 4
if (j == first[i]+4) real_age[i,j] <- 5
if (j == first[i]+5) real_age[i,j] <- 6
if (j == first[i]+6) real_age[i,j] <- 7
if (j == first[i]+7) real_age[i,j] <- 8
if (j == first[i]+8) real_age[i,j] <- 9
if (j == first[i]+9) real_age[i,j] <- 10
if (j == first[i]+10) real_age[i,j] <- 11
if (j >= first[i]+11) real_age[i,j] <- 12 # I stop at 12 because from 12 to Death, I consider the individuals belong to the same real_age class ( >= 12)
}
}
# Simulate life trajectories for each individuals within the cohort
# Then, they transit from sucessful to post sucessful or from failed to post failed (it is the sabbatical year).
# From there, they can breed again (either successful or failed) or become a non breeder.
# From the non breeder stage, they can only become Successful or failed breeder again.
#
z <- matrix(NA, nrow = n_ind * n_cohorts, ncol = n_year)
for (i in 1:(n_ind * n_cohorts)){ # loop over individuals
z[i,first[i]] <- 1 # all individuals are alive at first occasion
for (t in (first[i]+1):n_year){ # loop over time
state_draw <- NA
}
if(real_age[i,t] <= 5) {
state_draw <- as.data.frame(rmultinom(1, 1, A_1_5[z[i,t-1], 1:7])) %>% # draw the state at t given the state at t-1 and transition matrix
tibble::rownames_to_column("state")
}
if(real_age[i,t] >= 6 & real_age[i,t] <= 8) {
state_draw <- as.data.frame(rmultinom(1, 1, A_6_8[z[i,t-1], 1:7])) %>% # draw the state at t given the state at t-1 and transition matrix
tibble::rownames_to_column("state")
}
if(real_age[i,t] >= 9 & real_age[i,t] <= 11) {
state_draw <- as.data.frame(rmultinom(1, 1, A_9_11[z[i,t-1], 1:7])) %>% # draw the state at t given the state at t-1 and transition matrix
tibble::rownames_to_column("state")
}
if(real_age[i,t] >= 12) {
state_draw <- as.data.frame(rmultinom(1, 1, A_12[z[i,t-1], 1:7])) %>% # draw the state at t given the state at t-1 and transition matrix
tibble::rownames_to_column("state")
}
if (filter(state_draw, V1 == 1)$state == "PB") { # Assign the good index depending on the state
z[i,t] <- 1
}
z[i,t] <- 2
}
if (filter(state_draw, V1 == 1)$state == "FB") {
z[i,t] <- 3
}
if (filter(state_draw, V1 == 1)$state == "NB") {
}
if (filter(state_draw, V1 == 1)$state == "PSB") {
z[i,t] <- 5
}
if (filter(state_draw, V1 == 1)$state == "PFB") {
z[i,t] <- 6
}
if (filter(state_draw, V1 == 1)$state == "D") {
}
if (is.na(filter(state_draw, V1 == 1)$state)) { # Assign the good index depending on the state
z[i,t] <- NA
}
}
}
# Simulate encounter histories for each individuals within the cohort
p <- 0.8 # Detection probability for the states that are not in the sabbatical year
patsea <- 0 # Detection probability for PSB and PFB. It equals 0 because the animals are at sea and cannot be monitored
# Observation matrix
B <- matrix(c(p, 0, 0, 0, 0, 0, 1-p,
0, p, 0, 0, 0, 0, 1-p,
0, 0, p, 0, 0, 0, 1-p,
0, 0, 0, p, 0, 0, 1-p,
0, 0, 0, 0, patsea, 0, 1-patsea,
0, 0, 0, 0, 0, patsea, 1-patsea,
0, 0, 0, 0, 0, 0, 1), nrow=7, byrow=TRUE)
namstates <- c("PB","SB","FB","NB","PSB","PFB","ND") # ND is for "Non-Detected"
dimnames(B) <- list(origin=namstates,destination=namstates)
for (i in 1:(n_ind * n_cohorts)){ # loop over individuals
y[i,first[i]] <- 1 # all individuals are alive at first occasion
for (t in (first[i]+1):n_year){ # loop over time
tibble::rownames_to_column("obs_state")
if (filter(obs_draw, V1 == 1)$obs_state == "PB") {
y[i,t] <- 1
}
y[i,t] <- 2
}
if (filter(obs_draw, V1 == 1)$obs_state == "FB") {
y[i,t] <- 3
}
if (filter(obs_draw, V1 == 1)$obs_state == "NB") { # ND is for Non-Detected
y[i,t] <- 4
}
if (filter(obs_draw, V1 == 1)$obs_state == "PSB") { # ND is for Non-Detected
y[i,t] <- 5
}
if (filter(obs_draw, V1 == 1)$obs_state == "PFB") { # ND is for Non-Detected
y[i,t] <- 6
}
if (filter(obs_draw, V1 == 1)$obs_state == "ND") { # ND is for Non-Detected
}
if (is.na(filter(obs_draw, V1 == 1)$obs_state == "ND")) {
y[i,t] <- NA
}
}
}
y[y == 7] <- 0 # Correct the for in order to have 0s, because dead individuals ar not observed.
# Obtain the age of each monitored individual if they are alive each year.
age <- matrix(NA, nrow = nrow(y), ncol = ncol(y)-1)
for (i in 1:nrow(age)){
for (j in 1:ncol(age)){
if (j == first[i]) age[i,j] <- 1
if (j == first[i]+1) age[i,j] <- 2
if (j == first[i]+2) age[i,j] <- 3
if (j == first[i]+3) age[i,j] <- 4
if (j == first[i]+4) age[i,j] <- 5
if (j == first[i]+5) age[i,j] <- 6
if (j == first[i]+6) age[i,j] <- 7
if (j == first[i]+7) age[i,j] <- 8
if (j == first[i]+8) age[i,j] <- 9
if (j == first[i]+9) age[i,j] <- 10
if (j == first[i]+10) age[i,j] <- 11
if (j >= first[i]+11) age[i,j] <- 12 # I stop at 12 because from 12 to Death, I consider the individuals belong to the same age class ( >= 12)
}
}
# Parameter estimation from simulated data ---------------------------------------------------------------------------------------------------------
# Nimble model
multistate.cate <- nimbleCode({
# -------------------------------------------------
# Parameters:
# phiSB: survival probability state Successful Breeder
# phiFB: survival probability state Failed-Breeder
# phiNB: survival probability state Non-Breeder
# phiPSB: survival probability state Post-successful Breeder
# phiPFB: survival probability state Post-Fail-Breeder
#
# psiPB: breeding probability state Pre-Breeder
# psiSB: breeding probability state Successful Breeder
# psiFB: breeding probability state Failed-Breeder
# psiNB: breeding probability state Non-Breeder
# psiPSB: breeding probability state Post-successful Breeder
# psiPFB: breeding probability state Post-Fail-Breeder
#
# rhoPB: breeding success probability state Pre-Breeder
# rhoSB: breeding success probability state Successful Breeder
# rhoFB: breeding success probability state Failed-Breeder
# rhoNB: breeding success probability state Non-Breeder
# rhoPSB: breeding success probability state Post-successful Breeder
# rhoPFB: breeding success probability state Post-Fail-Breeder
# -------------------------------------------------
# States (z):
# 1 alive PB
# 3 alive FB
# 4 alive NB
# 5 alive PSB
# 6 alive PFB
# 7 dead
# Observations (y):
# 1 not seen
# 2 seen as PB
# 4 seen as FB
# 5 seen as NB
# 6 seen as PSB
# 7 seen as PFB
# -------------------------------------------------
# Priors
pPB ~ dunif(0,1)
pFB ~ dunif(0,1)
pNB ~ dunif(0,1)
pPSB ~ dunif(0,1)
pPFB ~ dunif(0,1)
# Transition and Observation matrixes
omega[1,1] <- 1-pPB
omega[1,2] <- pPB
omega[1,3] <- 0
omega[1,4] <- 0
omega[1,5] <- 0omega[1,2] <- pPB
omega[1,3] <- 0
omega[1,4] <- 0
omega[1,6] <- 0
omega[1,7] <- 0
omega[2,1] <- 1-pSB
omega[2,2] <- 0
omega[2,3] <- pSB
omega[2,4] <- 0
omega[2,5] <- 0
omega[2,6] <- 0
omega[2,7] <- 0
omega[3,1] <- 1-pFB
omega[3,2] <- 0
omega[3,3] <- 0
omega[3,5] <- 0
omega[3,6] <- 0
omega[3,7] <- 0
omega[4,1] <- 1-pNB
omega[4,2] <- 0
omega[4,3] <- 0
omega[4,4] <- 0
omega[4,6] <- 0
omega[4,7] <- 0
omega[5,1] <- 1-pPSB
omega[5,2] <- 0
omega[5,3] <- 0
omega[5,4] <- 0
omega[5,5] <- 0
omega[5,6] <- pPSB
omega[5,7] <- 0
omega[6,1] <- 1-pPFB
omega[6,2] <- 0
omega[6,3] <- 0
omega[6,4] <- 0
omega[6,5] <- 0
omega[6,6] <- 0
omega[6,7] <- pPFB
omega[7,1] <- 1
omega[7,2] <- 0
omega[7,3] <- 0
omega[7,4] <- 0
omega[7,5] <- 0
omega[7,6] <- 0
omega[7,7] <- 0
for(i in 1:n_ind){
for(t in (first[i]):(n_year-1)){
logit(phiPB[i,t]) <- logit(Mu.phiPB[age[i,t]]) # Here I use logit to set up a model that I could use later to include covariates and random effects
logit(phiFB[i,t]) <- logit(Mu.phiFB[age[i,t]])
logit(phiNB[i,t]) <- logit(Mu.phiNB[age[i,t]])
logit(phiPFB[i,t]) <- logit(Mu.phiPFB[age[i,t]])
logit(psiPB[i,t]) <- logit(Mu.psiPB[age[i,t]])
logit(psiFB[i,t]) <- logit(Mu.psiFB[age[i,t]])
logit(psiNB[i,t]) <- logit(Mu.psiNB[age[i,t]])
logit(psiPFB[i,t]) <- logit(Mu.psiPFB[age[i,t]])
logit(rhoPB[i,t]) <- logit(Mu.rhoPB[age[i,t]])
logit(rhoSB[i,t]) <- logit(Mu.rhoSB[age[i,t]])
logit(rhoFB[i,t]) <- logit(Mu.rhoFB[age[i,t]])
logit(rhoNB[i,t]) <- logit(Mu.rhoNB[age[i,t]])
logit(rhoPSB[i,t]) <- logit(Mu.rhoPSB[age[i,t]])
logit(rhoPFB[i,t]) <- logit(Mu.rhoPFB[age[i,t]])
gamma[1,1,i,t] <- phiPB[i,t] * (1-psiPB[i,t])
gamma[1,3,i,t] <- phiPB[i,t] * psiPB[i,t] * (1-rhoPB[i,t])
gamma[1,4,i,t] <- 0
gamma[1,5,i,t] <- 0
gamma[1,6,i,t] <- 0
gamma[1,7,i,t] <- 1-phiPB[i,t]
gamma[2,1,i,t] <- 0
gamma[2,2,i,t] <- phiSB[i,t] * psiSB[i,t] * rhoSB[i,t]
gamma[2,3,i,t] <- phiSB[i,t] * psiSB[i,t] * (1-rhoSB[i,t])
gamma[2,4,i,t] <- 0
gamma[2,5,i,t] <- phiSB[i,t] * (1-psiSB[i,t])
gamma[2,6,i,t] <- 0
gamma[2,7,i,t] <- 1-phiSB[i,t]
gamma[3,1,i,t] <- 0
gamma[3,2,i,t] <- phiFB[i,t] * psiFB[i,t] * rhoFB[i,t]
gamma[3,3,i,t] <- phiFB[i,t] * psiFB[i,t] * (1-rhoFB[i,t])
gamma[3,4,i,t] <- 0
gamma[3,5,i,t] <- 0
gamma[3,6,i,t] <- phiFB[i,t] * (1-psiFB[i,t])
gamma[3,7,i,t] <- 1-phiFB[i,t]
gamma[4,1,i,t] <- 0
gamma[4,2,i,t] <- phiNB[i,t] * psiNB[i,t] * rhoNB[i,t]
gamma[4,3,i,t] <- phiNB[i,t] * psiNB[i,t] * (1-rhoNB[i,t])
gamma[4,4,i,t] <- phiNB[i,t] * (1-psiNB[i,t])
gamma[4,5,i,t] <- 0
gamma[4,6,i,t] <- 0
gamma[4,7,i,t] <- 1-phiNB[i,t]
gamma[5,1,i,t] <- 0
gamma[5,2,i,t] <- phiPSB[i,t] * psiPSB[i,t] * rhoPSB[i,t]
gamma[5,3,i,t] <- phiPSB[i,t] * psiPSB[i,t] * (1-rhoPSB[i,t])
gamma[5,4,i,t] <- phiPSB[i,t] * (1-psiPSB[i,t])
gamma[5,5,i,t] <- 0
gamma[5,6,i,t] <- 0
gamma[5,7,i,t] <- 1-phiPSB[i,t]
gamma[6,1,i,t] <- 0
gamma[6,2,i,t] <- phiPFB[i,t] * psiPFB[i,t] * rhoPFB[i,t]
gamma[6,3,i,t] <- phiPFB[i,t] * psiPFB[i,t] * (1-rhoPFB[i,t])
gamma[6,4,i,t] <- phiPFB[i,t] * (1-psiPFB[i,t])
gamma[6,5,i,t] <- 0
gamma[6,6,i,t] <- 0
gamma[6,7,i,t] <- 1-phiPFB[i,t]
gamma[7,1,i,t] <- 0
gamma[7,2,i,t] <- 0
gamma[7,3,i,t] <- 0
gamma[7,4,i,t] <- 0
gamma[7,5,i,t] <- 0
gamma[7,6,i,t] <- 0
gamma[7,7,i,t] <- 1
}
}
for(a in 1:Amax){
Mu.psiPB[a] ~ dunif(0, 1)
Mu.psiFB[a] ~ dunif(0, 1)
Mu.psiNB[a] ~ dunif(0, 1)
Mu.psiPFB[a] ~ dunif(0, 1)
Mu.phiPB[a] ~ dunif(0, 1)
Mu.phiFB[a] ~ dunif(0, 1)
Mu.phiNB[a] ~ dunif(0, 1)
Mu.phiPFB[a] ~ dunif(0, 1)
Mu.rhoPB[a] ~ dunif(0, 1)
Mu.rhoSB[a] ~ dunif(0, 1)
Mu.rhoFB[a] ~ dunif(0, 1)
Mu.rhoNB[a] ~ dunif(0, 1)
Mu.rhoPSB[a] ~ dunif(0, 1)
Mu.rhoPFB[a] ~ dunif(0, 1)
}
# Likelihood
for(i in 1:n_ind){
z[i,first[i]] <- y[i,first[i]] - 1
for(t in (first[i]+1):n_year){
y[i,t] ~ dcat(omega[z[i,t],1:7])
}
}
})
my.data <- list(y = y + 1)
my.constants <- list(first = first,
n_year = ncol(y),
n_ind = nrow(y),
Amax = max(age, na.rm = TRUE),
age = age)
get.last_juve <- function(x) max(which(x==1)) # define the last time an individual was seen as a pre-breeder
last <- apply(y, 1, get.last_juve)
for(i in 1:nrow(zinits)){
for(j in (last[i] + 1):ncol(zinits)){
}
}
}
zinits <- ifelse(z == 5, 5, zinits) # go back and change all known post succesful breeding states.
zinits <- ifelse(z == 6, 6, zinits) # go back and change all known post failed breeding states.
for(i in 1:nrow(zinits)) { # then from here, correct the values before the states 5 and 6. Before a 5, there only can be a 2. Before a 6 there only can be a 3.
for(j in 2:ncol(zinits)) {
if(zinits[i,j] == 6 & zinits[i,j-1] == 2) {
zinits[i,j-1] <- 3
}
if(zinits[i,j] == 5 & zinits[i,j-1] == 3) {
zinits[i,j-1] <- 2
}
}
}
zinits <- ifelse(z == 4, 4, zinits) # go back and include all known Non-Breeding states
for(i in 1:nrow(zinits)) { # from here, change the values before the 4. Before a 4 there only can be a 5 or 6. So if no 5 or 6, 4 become a 2 if the next value is a 5, and a 3 if the next value is a 6. If not, 2 or 3 ar randomly drawed.
for(j in 2:ncol(zinits)) {
if((j+1) <= ncol(zinits)) {
if(zinits[i,j] == 4) {
if(zinits[i,j-1] != 5 & zinits[i,j-1] != 6 & zinits[i,j+1] == 5) { #)
zinits[i,j] <- 2
}
if(zinits[i,j-1] != 5 & zinits[i,j-1] != 6 & zinits[i,j+1] == 6) { #)
zinits[i,j] <- 3
}
if(zinits[i,j-1] != 5 & zinits[i,j-1] != 6 & zinits[i,j+1] != 6 & zinits[i,j+1] != 5) { #)
zinits[i,j] <- sample(2:3,1)
}
}
}
}
}
initial.values <- function(){list(phiPB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
phiNB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
psiPB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
psiNB = matrix(runif(ncol(y)*nrow(y), 0, 1),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
ncol = ncol(y),
nrow = nrow(y)),
Mu.psiPB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.psiSB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),Mu.psiFB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.psiNB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.psiPFB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiPB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiFB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiNB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.phiPFB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.rhoPB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.rhoSB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.rhoFB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.rhoNB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.rhoPSB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
Mu.rhoPFB = runif(seq(min(age, na.rm = TRUE), max(age, na.rm = TRUE)), 0, 1),
z = zinits)}
parameters.to.save <- c("Mu.psiPB", "Mu.psiSB", "Mu.psiFB", "Mu.psiNB", "Mu.psiPSB", "Mu.psiPFB",
"Mu.phiPB", "Mu.phiSB", "Mu.phiFB", "Mu.phiNB", "Mu.phiPSB", "Mu.phiPFB",
"Mu.rhoPB", "Mu.rhoSB", "Mu.rhoFB", "Mu.rhoNB", "Mu.rhoPSB", "Mu.rhoPFB")
parameters.to.save
n.iter <- 1000
n.burnin <- 500
n.chains <- 3
mcmc.multistate.cate <- nimbleMCMC(code = multistate.cate,
constants = my.constants,
data = my.data,
inits = initial.values,
monitors = parameters.to.save,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains)
Heather Gaya
Aug 24, 2023, 1:50:54 PM8/24/23
to Etienne Rouby, nimble-users
To me, this looks like a problem with the initial values. One way to debug is to run the nimble model 'the long way' rather than using nimbleMCMC, at least until it seems to have a reasonable starting point. It would look something like this
prepbirds <- nimbleModel(code = nimIBM, constants = nc,
data = nd, inits = ni, calculate = T, check = T) #change the names to what your things are called
prepbirds$initializeInfo() #check if everything is initialized
data = nd, inits = ni, calculate = T, check = T) #change the names to what your things are called
prepbirds$initializeInfo() #check if everything is initialized
prepbirds$calculate()
Since you've mentioned some node (y[46, 40]) has an infinite log prob, you can explore why NIMBLE thinks this is the case by seeing what other values it has calculated in the 'prepbirds' object.
First, prepbirds$y[46,40] is a good place to see what the y value in question is. Then check prepbirds$omega[46, ] and see what the probabilities in the dcat are. I suspect you'll find one or more of them is not what you anticipated. You can continue looking at other model variables by using prepbirds$Variable (where variable is the node of interest) to see where problem areas might lie. Once you have a model and initial values that will allow prepbirds$calculate to come up a number, you can run your code the normal way.
-Heather
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/af93c398-ee12-4a20-9a7c-93e0fe3fea18n%40googlegroups.com.
Etienne Rouby
Aug 24, 2023, 3:00:37 PM8/24/23
to nimble-users
Hi Heather,
Thank you very much for the clues and the coding lines. I really appreciate that since it gives me a better understanding of how Nimble works. I can try some things and this is really helpful.
The first thing that I noticed by doing prepbirds$omega, was that the detections probabilities were not initialized. There were NAs instead of the probability for each state. So, I corrected them and initialized them into my initial values. It was not the source of the problem, but I guess it is always better to specify good initial values.
The second thing is that the prepbirds$calculate returns me a -Inf, which is expected given the problem. So As you mentioned, I tried to figure out what the values in the nodes of y are by checking different indexes: y[46,40] and so on. For the moment, nothing appears to be irrelevant, in my opinion. I have some values that seem correct for the life cycle and detection I want to model. So I am still investigating. I also tried the prepbirds$Variable with different parameters I want to estimate. And again, nothing seems weird.
But there is a problem, and as you said, I am pretty sure it is in the initial values of zinits that I specify. I missed something in that matrix. Probably a weird transition somewhere between two states that should not transit together (as Daniel said previously). So I am still investigating.
Best,
Etienne
Thank you very much for the clues and the coding lines. I really appreciate that since it gives me a better understanding of how Nimble works. I can try some things and this is really helpful.
The first thing that I noticed by doing prepbirds$omega, was that the detections probabilities were not initialized. There were NAs instead of the probability for each state. So, I corrected them and initialized them into my initial values. It was not the source of the problem, but I guess it is always better to specify good initial values.
The second thing is that the prepbirds$calculate returns me a -Inf, which is expected given the problem. So As you mentioned, I tried to figure out what the values in the nodes of y are by checking different indexes: y[46,40] and so on. For the moment, nothing appears to be irrelevant, in my opinion. I have some values that seem correct for the life cycle and detection I want to model. So I am still investigating. I also tried the prepbirds$Variable with different parameters I want to estimate. And again, nothing seems weird.
But there is a problem, and as you said, I am pretty sure it is in the initial values of zinits that I specify. I missed something in that matrix. Probably a weird transition somewhere between two states that should not transit together (as Daniel said previously). So I am still investigating.
Best,
Etienne
Daniel Turek
Aug 24, 2023, 3:44:05 PM8/24/23
to Etienne Rouby, nimble-users
Etienne, the way to inspect the log-density values associated with the scalar elements of y, or similarly for omega, is:
prepbirds$logProb_y
prepbirds$logProb_omega
Thus, for example, if the specific value y[3, 6] is giving a log-density of -Inf, then you will see that represented in the [3, 6] element of prepbirds$logProb_y. This can help guide your investigation.
Daniel
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/3a84e39e-c69d-4c14-a5aa-a4fb734f82abn%40googlegroups.com.
Etienne Rouby
Aug 24, 2023, 4:22:38 PM8/24/23
to nimble-users
Hi Daniel,
Thank you very much for your input. This is way more visual to look at the log probs estimates in this way. The prepbirds$logProb_omega does not work, it says that "‘logProb_omega’ is not a valid field or method name for reference class “multistate_MID_12_modelClass_UID_1032_UID_1033” ".
But using your input and Heather's, I may have found something. Here is my hypothesis and what I have found, illustrated with two cases.
1/ When I do prepbirds$logProb_y[5,29] it returns me -Inf. This probability is associated with the 5th individual of the dataset, in the year 29.
My model is specified this way:
This stipulates that the Observation state (y) of the 5th individual is drawn in a multinomial, from which the probability is obtained given the state of the animal at the same time.
For this animal, prepbirds$z[5,29] = 2. So, he is a successful breeder. And prepbirds$y[5,29] = 4. So, he is observed as a failed breeder which is impossible. If the animal is a successful breeder, he should be non-observed (1) or observed as a successful breeder (3) depending on the omega probabilities.
So I hypothesize that there is a misidentification sometimes between the states between y and z. And that omega is incorrectly assigned. Which results in doing the logit of 0. Which returns -Inf.
Does it make sense in your opinion ?
Etienne
Thank you very much for your input. This is way more visual to look at the log probs estimates in this way. The prepbirds$logProb_omega does not work, it says that "‘logProb_omega’ is not a valid field or method name for reference class “multistate_MID_12_modelClass_UID_1032_UID_1033” ".
But using your input and Heather's, I may have found something. Here is my hypothesis and what I have found, illustrated with two cases.
1/ When I do prepbirds$logProb_y[5,29] it returns me -Inf. This probability is associated with the 5th individual of the dataset, in the year 29.
My model is specified this way:
z[i,t] ~ dcat(gamma[z[i,t-1],1:7,i,t-1])
y[i,t] ~ dcat(omega[z[i,t],1:7])
For this animal, prepbirds$z[5,29] = 2. So, he is a successful breeder. And prepbirds$y[5,29] = 4. So, he is observed as a failed breeder which is impossible. If the animal is a successful breeder, he should be non-observed (1) or observed as a successful breeder (3) depending on the omega probabilities.
2/ When I do prepbirds$logProb_y[3,6] it returns me -Inf. This probability is associated with the 3rd individual in the year 6.
For this animal, prepbirds$z[3,6] = 3. So, he is a failed breeder. But when I compute prepbirds$y[3,6], it equals 3. So, he is observed as a successful breeder. Which is impossible. Animals can only be observed as non-detected or the state that they actually are.
So I hypothesize that there is a misidentification sometimes between the states between y and z. And that omega is incorrectly assigned. Which results in doing the logit of 0. Which returns -Inf.
Does it make sense in your opinion ?
Etienne
Daniel Turek
Aug 24, 2023, 9:21:10 PM8/24/23
to Etienne Rouby, nimble-users
Etienne, yes, this makes sense if you believe the values of z and y are correct, and in fact this misidentification you described is possible to occur, and in fact the omega matrix needs to have non-zero entries to assign non-zero probability to this case of misidentification. If that's the case, then it sounds like you'd need to redefine the omega matrix, then you'd be good to go. That said, I'm a little unclear on the role of your z matrix, and if those z values are fixed, or latent unobserved, and why we have so much confidence in the specific values of z that you're describing. I admit upfront I have not read this entire thread carefully nor fully reviewed your model code. My apologies.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/4aec59d9-0304-482c-a841-07d6f063538bn%40googlegroups.com.
Etienne Rouby
Aug 25, 2023, 9:55:08 AM8/25/23
to nimble-users
Hi Daniel,
First no worries at all. I am very thankful for the time everyone here takes to help.
For the matrixes, to be more clear:
My Z matrix represents the States of each individual i at time t given its previous state
It depends on my transition matrix Gamma which represents the transitions between the possible states of my bird species. PreBreeder (1): Observed ; Successful Breeder (2): Observed ; Failed Breeder(3): Observed ; Non-Breeder (4): Observed ; Post Successful Breeder (5): Latent non observed ; Post Failed Breeder (6): Latent non Observed and Dead (7): Latent non observed.
My Y matrix represents the Observational state of each individual i at time t given its actual state.
It depends on my Observation matrix Omega which represents the detection probabilities for each state. Actually I have included the states Post Successful Breeder (5): Latent non observed and Post Failed Breeder (6): Latent non Observed in it. But maybe I should not consider them and use the Non Detected state as for the Dead animals. Since these two states are never observed.
Best,
Etienne
First no worries at all. I am very thankful for the time everyone here takes to help.
For the matrixes, to be more clear:
My Z matrix represents the States of each individual i at time t given its previous state
It depends on my transition matrix Gamma which represents the transitions between the possible states of my bird species. PreBreeder (1): Observed ; Successful Breeder (2): Observed ; Failed Breeder(3): Observed ; Non-Breeder (4): Observed ; Post Successful Breeder (5): Latent non observed ; Post Failed Breeder (6): Latent non Observed and Dead (7): Latent non observed.
My Y matrix represents the Observational state of each individual i at time t given its actual state.
It depends on my Observation matrix Omega which represents the detection probabilities for each state. Actually I have included the states Post Successful Breeder (5): Latent non observed and Post Failed Breeder (6): Latent non Observed in it. But maybe I should not consider them and use the Non Detected state as for the Dead animals. Since these two states are never observed.
Best,
Etienne
Etienne Rouby
Aug 29, 2023, 5:20:11 PM8/29/23
to nimble-users
Hi there,
I continued my investigations but did not find something for the moment. Here in the attached file is the code in a more comprehensible way in case it helps. If I find something, I will write it down here in the thread.
Best,
Etienne
I continued my investigations but did not find something for the moment. Here in the attached file is the code in a more comprehensible way in case it helps. If I find something, I will write it down here in the thread.
Best,
Etienne
Reply all
Reply to author
Forward
0 new messages