Fit model combining live recapture and dead recoveries in Jags using SSM formulation of multistate models

lifedead_og.R

# Fit model combining live recapture and dead recoveries in Jags using SSM formulation of multistate models
# Warning: if initial values are generated from the ld.init() function, we failed at recovering the values
# used to simulate the data; everything goes well if we use the true latent states (data are simulated) as initial values 

# Define function to simulate multistate capture-recapture data
simul.ms <- function(PSI.STATE, PSI.OBS, marked, unobservable = NA){
   # Unobservable: number of state that is unobservable
   n.occasions <- dim(PSI.STATE)[4] + 1
   CH <- CH.TRUE <- matrix(NA, ncol = n.occasions, nrow = sum(marked))
   # Define a vector with the occasion of marking
   mark.occ <- matrix(0, ncol = dim(PSI.STATE)[1], nrow = sum(marked))
   g <- colSums(marked)
   for (s in 1:dim(PSI.STATE)[1]){
      if (g[s]==0) next  # To avoid error message if nothing to replace
      mark.occ[(cumsum(g[1:s])-g[s]+1)[s]:cumsum(g[1:s])[s],s] <-
      rep(1:n.occasions, marked[1:n.occasions,s])
      } #s
   for (i in 1:sum(marked)){
      for (s in 1:dim(PSI.STATE)[1]){
         if (mark.occ[i,s]==0) next
         first <- mark.occ[i,s]
         CH[i,first] <- s
         CH.TRUE[i,first] <- s
         } #s
      for (t in (first+1):n.occasions){
         # Multinomial trials for state transitions
         if (first==n.occasions) next
         state <- which(rmultinom(1, 1, PSI.STATE[CH.TRUE[i,t-1],,i,t-1])==1)
         CH.TRUE[i,t] <- state
         # Multinomial trials for observation process
         event <- which(rmultinom(1, 1, PSI.OBS[CH.TRUE[i,t],,i,t-1])==1)
         CH[i,t] <- event
         } #t
      } #i
   # Replace the NA and the highest state number (dead) in the file by 0
   CH[is.na(CH)] <- 0
   CH[CH==dim(PSI.STATE)[1]] <- 0
   CH[CH==unobservable] <- 0
   id <- numeric(0)
   for (i in 1:dim(CH)[1]){
      z <- min(which(CH[i,]!=0))
      ifelse(z==dim(CH)[2], id <- c(id,i), id <- c(id))
      }
   return(list(CH=CH[-id,], CH.TRUE=CH.TRUE[-id,]))
   # CH: capture histories to be used
   # CH.TRUE: capture histories with perfect observation
   }




# Define mean survival, transitions, recapture, as well as number of occasions, states, observations and released individuals 
ss <- 0.8
ff <- 0.6
rr <- 0.1
pp <- 0.5
n.occasions <- 10  
n.states <- 4
n.obs <- 3
marked <- matrix(0, ncol = n.states, nrow = n.occasions)
marked[,1] <- rep(100, n.occasions)	# Releases in study area

# Define matrices with survival, transition and recapture probabilities
# These are 4-dimensional matrices, with 
   # Dimension 1: state of departure
   # Dimension 2: state of arrival
   # Dimension 3: individual
   # Dimension 4: time
# 1. State process matrix
totrel <- sum(marked)*(n.occasions-1)
PSI.STATE <- array(NA, dim=c(n.states, n.states, totrel, n.occasions-1))
for (i in 1:totrel){
   for (t in 1:(n.occasions-1)){
      PSI.STATE[,,i,t] <- matrix(c(
      ss*ff, ss*(1-ff), 1-ss, 0,
      0,   ss,       1-ss, 0,
      0,   0,       0,   1,
      0,   0,       0,   1), nrow = n.states, byrow = TRUE)
      } #t
   } #i

# 2.Observation process matrix
PSI.OBS <- array(NA, dim=c(n.states, n.obs, totrel, n.occasions-1))
for (i in 1:totrel){
   for (t in 1:(n.occasions-1)){
      PSI.OBS[,,i,t] <- matrix(c(
      pp, 0, 1-pp,
      0, 0, 1,
      0, rr, 1-rr,
      0, 0, 1), nrow = n.states, byrow = TRUE)
      } #t
   } #i

# Execute simulation function
sim <- simul.ms(PSI.STATE, PSI.OBS, marked)
CH <- sim$CH
zsimul <- sim$CH.TRUE

# Compute date of first capture
get.first <- function(x) min(which(x!=0))
first <- apply(CH, 1, get.first)

# Recode CH matrix: note, a 0 is not allowed!
# 1 = alive and in study are, 2 = recovered dead, 3 = not seen or recovered
rCH <- CH  # Recoded CH
rCH[rCH==0] <- 3


# Specify model in BUGS language
sink("lifedead.jags")
cat("
model {

# -------------------------------------------------
# Parameters:
# ss: true survival probability
# ff: fidelity probability
# rr: recovery probability
# pp: recapture/resighting probability
# -------------------------------------------------
# States (S):
# 1 alive in study area
# 2 alive outside study area
# 3 recently dead and recovered
# 4 recently dead, but not recovered, or dead (absorbing)
# Observations (O):
# 1 seen alive
# 2 recovered dead
# 3 neither seen nor recovered
# -------------------------------------------------

# Priors
ss ~ dunif(0, 1)     # Prior for mean survival
ff ~ dunif(0, 1)     # Prior for mean fidelity
rr ~ dunif(0, 1)     # Prior for mean recovery
pp ~ dunif(0, 1)     # Prior for mean recapture

# Define state-transition and observation matrices 	
   # Define probabilities of state S(t+1) given S(t)
      ps[1,1] <- ss * ff
      ps[1,2] <- ss * (1 - ff)
      ps[1,3] <- 1 - ss
      ps[1,4] <- 0
      ps[2,1] <- 0
      ps[2,2] <- ss
      ps[2,3] <- 1 - ss
      ps[2,4] <- 0
      ps[3,1] <- 0
      ps[3,2] <- 0
      ps[3,3] <- 0
      ps[3,4] <- 1
      ps[4,1] <- 0
      ps[4,2] <- 0
      ps[4,3] <- 0
      ps[4,4] <- 1

      # Define probabilities of O(t) given S(t)
      po[1,1] <- pp
      po[1,2] <- 0
      po[1,3] <- 1 - pp
      po[2,1] <- 0
      po[2,2] <- 0
      po[2,3] <- 1
      po[3,1] <- 0
      po[3,2] <- rr
      po[3,3] <- 1 - rr
      po[4,1] <- 0
      po[4,2] <- 0
      po[4,3] <- 1

# Likelihood 
for (i in 1:nind){
   # Define latent state at first capture
   z[i,first[i]] <- y[i,first[i]]
   for (t in (first[i]+1):n.occasions){
      # State process: draw S(t) given S(t-1)
      z[i,t] ~ dcat(ps[z[i,t-1], 1:4])
      # Observation process: draw O(t) given S(t)
      y[i,t] ~ dcat(po[z[i,t], 1:3])
      } #t
   } #i
}
",fill = TRUE)
sink()


# Bundle data
jags.data <- list(y = rCH, 
                  first = first, 
                  n.occasions = dim(rCH)[2], 
                  nind = dim(rCH)[1])

# Initial values
ld.init <- function(ch, f){
   ch[ch==3] <- NA
   v2 <- which(ch==2, arr.ind = T)
   ch[v2] <- 3
   for (i in 1:nrow(v2)){
      ifelse(v2[i,2]!=ncol(ch), ch[v2[i,1], (v2[i,2]+1):ncol(ch)] <- 4, next)}
   for (i in 1:nrow(ch)){
      m <- max(which(ch[i,]==1))
      ch[i,f[i]:m] <- 1
      }
   for (i in 1:nrow(v2)){
      u1 <- min(which(ch[v2[i,1],]==1))
      ch[v2[i],u1:(v2[i,2]-1)] <- 1
      }
   for (i in 1:nrow(ch)){
      for (j in f[i]:ncol(ch)){
         if(is.na(ch[i,j])==1) ch[i,j] <- 1
         }
      ch[i,f[i]] <- NA
      }
   return(ch)
   }

inits <- function(){list(ss = runif(1, 0, 1), 
                         ff = runif(1, 0, 1), 
                         pp = runif(1, 0, 1),
                         rr = runif(1, 0, 1), 
                         z = ld.init(rCH, first))}  

# instead of using the ld.init() function to generate inits for the latent states
# we use the true latent states

zinit <- zsimul
for (i in 1:nrow(zinit)){
   zinit[i, first[i]] <- NA
}

inits <- function(){list(ss = runif(1, 0, 1), 
                         ff = runif(1, 0, 1), 
                         pp = runif(1, 0, 1),
                         rr = runif(1, 0, 1), 
                         z = zinit)}  

# Parameters monitored
parameters <- c("ss", "ff", "rr", "pp")

# MCMC settings
ni <- 2000
nt <- 1
nb <- 1000
nc <- 3

# Call JAGS from R
library(jagsUI)
lifedead <- jags(data = jags.data, 
                 inits = inits, 
                 parameters.to.save = parameters, 
                 model.file = "lifedead.jags", 
                 n.chains = nc, 
                 n.thin = nt, 
                 n.iter = ni, 
                 n.burnin = nb, 
                 parallel = TRUE)

print(lifedead, digit = 3)

traceplot(lifedead)


# JAGS output for model 'lifedead.jags', generated by jagsUI.
# Estimates based on 3 chains of 2000 iterations,
# adaptation = 100 iterations (sufficient),
# burn-in = 1000 iterations and thin rate = 1,
# yielding 3000 total samples from the joint posterior. 
# MCMC ran in parallel for 1.723 minutes at time 2020-11-26 09:57:41.
# 
# mean     sd     2.5%      50%    97.5% overlap0 f  Rhat n.eff
# ss          0.788  0.008    0.772    0.788    0.805    FALSE 1 1.001  1461
# ff          0.581  0.024    0.535    0.580    0.631    FALSE 1 1.001  1308
# rr          0.102  0.013    0.078    0.102    0.127    FALSE 1 1.000  3000
# pp          0.485  0.031    0.423    0.485    0.546    FALSE 1 1.002  1020
# deviance 1336.051 48.203 1246.511 1335.452 1432.940    FALSE 1 1.001  1505
# 
# Successful convergence based on Rhat values (all < 1.1). 
# Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
# For each parameter, n.eff is a crude measure of effective sample size. 
# 
# overlap0 checks if 0 falls in the parameter's 95% credible interval.
# f is the proportion of the posterior with the same sign as the mean;
# i.e., our confidence that the parameter is positive or negative.
# 
# DIC info: (pD = var(deviance)/2) 
# pD = 1161 and DIC = 2497.048 
# DIC is an estimate of expected predictive error (lower is better).
#