appendix_paper_IHCMR_partB[rmd2r].R

#' ---	
#' title: 'Appendix B: Individual heterogeneity in capture-recapture models - Bayesian	
#'   approach using Jags'	
#' output:	
#'   word_document: default	
#'   html_notebook: default	
#'   pdf_document: default	
#'   html_document: default	
#' ---	
#' 	
#' ## Introduction	
#' 	
#' In this appendix, we introduce three methods to cope with individual heterogeneity in capture-recapture models, which we implement in a Bayesian framework using MCMC methods. First, we present multistate models in which heterogeneity is measured on individuals using states. Then, we illustrate models with individual random effects and finite mixtures that can help in dealing with hidden heterogeneity. We refer to the paper for a formal presentation of these models and a list of references using them. Throughout this appendix, we use R to simulate data and program Jags is called from R using package Rjags to fit models. We do our best to ensure reproducibility. Note that the frequentist approach can be used instead, and implemented either with program E-SURGE (appendix C) or program Mark (appendix A).	
#' 	
#' ## Multistate models	
#' 	
#' In this section, we aim at illustrating how not accounting for individual heterogeneity may obscure the detection of life-history tradeoffs. 	
#' In details, we consider two states for the individuals of our fake population, non-breeding (NB) and breeding (B). To mimic 	
#' individual heterogeneity, we simulate a bunch of good individuals with survival $\phi_{NB}=0.7$ and $\phi_{B}=0.8$ and a bunch of bad individuals with survival $\phi_{NB}=0.7$ and $\phi_{B}=0.6$. Overall, the cost of breeding on survival should be detected only in bad individuals after accounting for individual heterogeneity through quality. For each group of bad vs. good individuals, we consider the same detection probability $p=0.9$, the same transition probabilities between breeding states $\psi_{NB,B}=0.8$ and $\psi_{B,NB}=0.3$, and 100 newly marked individuals for each group in each year of the 6-year study.	
#' 	
#' ### Data simulation	
#' 	
#' Using R code from [Kéry and Schaub (2012)](http://www.vogelwarte.ch/de/projekte/publikationen/bpa/) book (chapter 9), we first define a function to simulate multistate capture-recapture data:	
#' 	
# 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	
}	
#' 	
#' 	
#' Second, we use this function to simulate the two datasets of good and bad individuals:	
#' 	
set.seed(1) # for reproducibility	
p = 0.9	
R = 100	
#------------------------------	
#---- good quality individuals	
#------------------------------	
# Define mean survival, transitions, recapture, as well as number of occasions, states, observations and released individuals	
phiA <- 0.7	
phiB <- 0.8	
psiAB <- 0.8	
psiBA <- 0.3	
pA <- p	
pB <- p	
n.occasions <- 6	
n.states <- 3	
n.obs <- 3	
marked <- matrix(NA, ncol = n.states, nrow = n.occasions)	
marked[,1] <- rep(R, n.occasions)	
marked[,2] <- rep(R, n.occasions)	
marked[,3] <- rep(0, n.occasions)	
# Define matrices with survival, transition and recapture probabilities	
# 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(	
      phiA*(1-psiAB), phiA*psiAB, 1-phiA,	
      phiB*psiBA, phiB*(1-psiBA), 1-phiB,	
      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(	
      pA, 0, 1-pA,	
      0, pB, 1-pB,	
      0, 0, 1 ), nrow = n.states, byrow = TRUE)	
  } #t	
} #i	
	
# Execute function	
sim <- simul.ms(PSI.STATE, PSI.OBS, marked)	
CH <- sim$CH	
his1 = CH[!apply(CH,1,sum)==0,] # remove lines of 0s	
	
#------------------------------	
#---- bad quality individuals	
#------------------------------	
# Define mean survival, transitions, recapture, as well as number of occasions, states, observations and released individuals	
phiA <- 0.7	
phiB <- 0.6	
psiAB <- 0.8	
psiBA <- 0.3	
pA <- p	
pB <- p	
n.occasions <- 6	
n.states <- 3	
n.obs <- 3	
marked <- matrix(NA, ncol = n.states, nrow = n.occasions)	
marked[,1] <- rep(R, n.occasions)	
marked[,2] <- rep(R, n.occasions)	
marked[,3] <- rep(0, n.occasions)	
# Define matrices with survival, transition and recapture probabilities	
# 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(	
    phiA*(1-psiAB), phiA*psiAB, 1-phiA,	
    phiB*psiBA, phiB*(1-psiBA), 1-phiB,	
    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(	
    pA, 0, 1-pA,	
    0, pB, 1-pB,	
    0, 0, 1 ), nrow = n.states, byrow = TRUE)	
  } #t	
} #i	
	
# Execute function	
sim <- simul.ms(PSI.STATE, PSI.OBS, marked)	
CH <- sim$CH	
his2 = CH[!apply(CH,1,sum)==0,] # remove lines of 0s	
#' 	
#' 	
#' Last, we pool these two datasets together:	
#' 	
his = rbind(his1,his2) 	
head(his) # display first lines	
tail(his) # display last lines	
#' 	
#' 	
#' ### Model fitting	
#' 	
#' We adapted R code from [Kéry and Schaub (2012)](http://www.vogelwarte.ch/de/projekte/publikationen/bpa/) book (chapter 9) to fit multistate capture-recapture data.	
#' 	
#' Let us get the occasion of first capture for each individual:	
#' 	
get.first <- function(x) min(which(x!=0))	
f <- apply(his, 1, get.first)	
f	
#' 	
#' 	
#' Recode the data such that 1 = seen alive in A, 2 = seen alive in B, 3 = not seen:	
#' 	
his_recoded <- his	
his_recoded[his_recoded==0] <- 3	
#' 	
#' 	
#' Now we fit a multistate model: we assume that survival depends on the breeding states, 	
#' transition probabilities are constant over time, as well as the detection probability:	
#' 	
sink("state_on_survival.jags")	
cat("	
model {	
	
# -------------------------------------------------	
# Parameters:	
# phiA: survival probability at site A	
# phiB: survival probability at site B	
# psiAB: movement probability from site A to site B	
# psiBA: movement probability from site B to site A	
# p: recapture probability at site A or B	
# -------------------------------------------------	
# States (S):	
# 1 alive at A	
# 2 alive at B	
# 3 dead	
# Observations (O):	
# 1 seen at A	
# 2 seen at B	
# 3 not seen	
# -------------------------------------------------	
	
# Priors	
   phiA ~ dunif(0, 1)	
   phiB ~ dunif(0, 1)	
   psiAB ~ dunif(0, 1)	
   psiBA ~ dunif(0, 1)	
   p ~ dunif(0, 1)	
	
# Define state-transition and observation matrices	
for (i in 1:nind){	
   # Define probabilities of state S(t+1) given S(t)	
   for (t in f[i]:(n.occasions-1)){	
      ps[1,i,t,1] <- phiA * (1-psiAB)	
      ps[1,i,t,2] <- phiA * psiAB	
      ps[1,i,t,3] <- 1-phiA	
      ps[2,i,t,1] <- phiB * psiBA	
      ps[2,i,t,2] <- phiB * (1-psiBA)	
      ps[2,i,t,3] <- 1-phiB	
      ps[3,i,t,1] <- 0	
      ps[3,i,t,2] <- 0	
      ps[3,i,t,3] <- 1	
	
      # Define probabilities of O(t) given S(t)	
      po[1,i,t,1] <- p	
      po[1,i,t,2] <- 0	
      po[1,i,t,3] <- 1-p	
      po[2,i,t,1] <- 0	
      po[2,i,t,2] <- p	
      po[2,i,t,3] <- 1-p	
      po[3,i,t,1] <- 0	
      po[3,i,t,2] <- 0	
      po[3,i,t,3] <- 1	
      } #t	
   } #i	
	
# Likelihood	
for (i in 1:nind){	
   # Define latent state at first capture	
   z[i,f[i]] <- y[i,f[i]]	
   for (t in (f[i]+1):n.occasions){	
      # State process: draw S(t) given S(t-1)	
      z[i,t] ~ dcat(ps[z[i,t-1], i, t-1,])	
      # Observation process: draw O(t) given S(t)	
      y[i,t] ~ dcat(po[z[i,t], i, t-1,])	
      } #t	
   } #i	
}	
",fill = TRUE)	
sink()	
#' 	
#' 	
#' 	
# Function to create known latent states z	
known.state.ms <- function(ms, notseen){	
   # notseen: label for ënot seení	
   state <- ms	
   state[state==notseen] <- NA	
   for (i in 1:dim(ms)[1]){	
      m <- min(which(!is.na(state[i,])))	
      state[i,m] <- NA	
      }	
   return(state)	
   }	
	
# Function to create initial values for unknown z	
ms.init.z <- function(ch, f){	
   for (i in 1:dim(ch)[1]){ch[i,1:f[i]] <- NA}	
   states <- max(ch, na.rm = TRUE)	
   known.states <- 1:(states-1)	
   v <- which(ch==states)	
   ch[-v] <- NA	
   ch[v] <- sample(known.states, length(v), replace = TRUE)	
   return(ch)	
   }	
	
# Bundle data	
jags.data <- list(y = his_recoded, f = f, n.occasions = dim(his_recoded)[2], nind = dim(his_recoded)[1], z = known.state.ms(his_recoded, 3))	
	
# Initial values	
inits <- function(){list(phiA = runif(1, 0, 1), psiAB = runif(1, 0, 1), p = runif(1, 0, 1), z = ms.init.z(his_recoded, f))}  	
	
# Parameters monitored	
parameters <- c("phiA", "phiB", "psiAB", "psiBA", "p")	
	
# MCMC settings	
ni <- 5000	
nb <- 2000	
nc <- 2	
	
# Call JAGS from R	
library(R2jags)	
ms <- jags(jags.data, inits, parameters, "state_on_survival.jags", n.chains = nc, n.iter = ni, n.burnin = nb)	
print(ms, digits = 3)	
#' 	
#' 	
#' 	
#' Run same model without state effect on survival:	
#' 	
sink("no_state_on_survival.jags")	
cat("	
model {	
	
# -------------------------------------------------	
# Parameters:	
# phi: survival probability at site A or B	
# psiAB: movement probability from site A to site B	
# psiBA: movement probability from site B to site A	
# p: recapture probability at site A or B	
# -------------------------------------------------	
# States (S):	
# 1 alive at A	
# 2 alive at B	
# 3 dead	
# Observations (O):	
# 1 seen at A	
# 2 seen at B	
# 3 not seen	
# -------------------------------------------------	
	
# Priors	
   phi ~ dunif(0, 1)	
   psiAB ~ dunif(0, 1)	
   psiBA ~ dunif(0, 1)	
   p ~ dunif(0, 1)	
	
# Define state-transition and observation matrices	
for (i in 1:nind){	
   # Define probabilities of state S(t+1) given S(t)	
   for (t in f[i]:(n.occasions-1)){	
      ps[1,i,t,1] <- phi * (1-psiAB)	
      ps[1,i,t,2] <- phi * psiAB	
      ps[1,i,t,3] <- 1-phi	
      ps[2,i,t,1] <- phi * psiBA	
      ps[2,i,t,2] <- phi * (1-psiBA)	
      ps[2,i,t,3] <- 1-phi	
      ps[3,i,t,1] <- 0	
      ps[3,i,t,2] <- 0	
      ps[3,i,t,3] <- 1	
	
      # Define probabilities of O(t) given S(t)	
      po[1,i,t,1] <- p	
      po[1,i,t,2] <- 0	
      po[1,i,t,3] <- 1-p	
      po[2,i,t,1] <- 0	
      po[2,i,t,2] <- p	
      po[2,i,t,3] <- 1-p	
      po[3,i,t,1] <- 0	
      po[3,i,t,2] <- 0	
      po[3,i,t,3] <- 1	
      } #t	
   } #i	
	
# Likelihood	
for (i in 1:nind){	
   # Define latent state at first capture	
   z[i,f[i]] <- y[i,f[i]]	
   for (t in (f[i]+1):n.occasions){	
      # State process: draw S(t) given S(t-1)	
      z[i,t] ~ dcat(ps[z[i,t-1], i, t-1,])	
      # Observation process: draw O(t) given S(t)	
      y[i,t] ~ dcat(po[z[i,t], i, t-1,])	
      } #t	
   } #i	
}	
",fill = TRUE)	
sink()	
#' 	
#' 	
#' 	
# Function to create known latent states z	
known.state.ms <- function(ms, notseen){	
   # notseen: label for ënot seení	
   state <- ms	
   state[state==notseen] <- NA	
   for (i in 1:dim(ms)[1]){	
      m <- min(which(!is.na(state[i,])))	
      state[i,m] <- NA	
      }	
   return(state)	
   }	
	
# Function to create initial values for unknown z	
ms.init.z <- function(ch, f){	
   for (i in 1:dim(ch)[1]){ch[i,1:f[i]] <- NA}	
   states <- max(ch, na.rm = TRUE)	
   known.states <- 1:(states-1)	
   v <- which(ch==states)	
   ch[-v] <- NA	
   ch[v] <- sample(known.states, length(v), replace = TRUE)	
   return(ch)	
   }	
	
# Bundle data	
jags.data <- list(y = his_recoded, f = f, n.occasions = dim(his_recoded)[2], nind = dim(his_recoded)[1], z = known.state.ms(his_recoded, 3))	
	
# Initial values	
inits <- function(){list(phi = runif(1, 0, 1), psiAB = runif(1, 0, 1), p = runif(1, 0, 1), z = ms.init.z(his_recoded, f))}  	
	
# Parameters monitored	
parameters <- c("phi", "psiAB", "psiBA", "p")	
	
# MCMC settings	
ni <- 5000	
nb <- 2000	
nc <- 2	
	
# Call JAGS from R	
library(R2jags)	
ms_without <- jags(jags.data, inits, parameters, "no_state_on_survival.jags", n.chains = nc, n.iter = ni, n.burnin = nb)	
print(ms_without, digits = 3)	
#' 	
#' 	
#' If you compare DIC values, it sounds like the difference in survival of breeding vs. non-breeding individuals is not detected.	
#' 	
#' Let's add individual heterogeneity through an individual covariate for bad vs. good individuals:	
#' 	
quality=c(rep(0,nrow(his1)),rep(1,nrow(his2))) # 0 for good, 1 for bad	
quality	
#' 	
#' 	
#' Now we fit again the two models from above, including the effect of individual heterogeneity.	
#' 	
#' 	
sink("state_on_survival_withquality.jags")	
cat("	
model {	
	
# -------------------------------------------------	
# Parameters:	
# alpha, beta: regression parameters for survival probability at site A/B	
# psiAB: movement probability from site A to site B	
# psiBA: movement probability from site B to site A	
# p: recapture probability	
# -------------------------------------------------	
# States (S):	
# 1 alive at A	
# 2 alive at B	
# 3 dead	
# Observations (O):	
# 1 seen at A	
# 2 seen at B	
# 3 not seen	
# -------------------------------------------------	
	
# Priors	
   alphaA ~ dnorm(0, 0.01)	
   betaA ~ dnorm(0, 0.01)	
   alphaB ~ dnorm(0, 0.01)	
   betaB ~ dnorm(0, 0.01)	
   for (i in 1:nind){	
   logit(phiA[i]) <- alphaA + betaA * quality[i]	
   logit(phiB[i]) <- alphaB + betaB * quality[i]}	
   psiAB ~ dunif(0, 1)	
   psiBA ~ dunif(0, 1)	
   p ~ dunif(0, 1)	
	
# Define state-transition and observation matrices	
for (i in 1:nind){	
   # Define probabilities of state S(t+1) given S(t)	
   for (t in f[i]:(n.occasions-1)){	
      ps[1,i,t,1] <- phiA[i] * (1-psiAB)	
      ps[1,i,t,2] <- phiA[i] * psiAB	
      ps[1,i,t,3] <- 1-phiA[i]	
      ps[2,i,t,1] <- phiB[i] * psiBA	
      ps[2,i,t,2] <- phiB[i] * (1-psiBA)	
      ps[2,i,t,3] <- 1-phiB[i]	
      ps[3,i,t,1] <- 0	
      ps[3,i,t,2] <- 0	
      ps[3,i,t,3] <- 1	
	
      # Define probabilities of O(t) given S(t)	
      po[1,i,t,1] <- p	
      po[1,i,t,2] <- 0	
      po[1,i,t,3] <- 1-p	
      po[2,i,t,1] <- 0	
      po[2,i,t,2] <- p	
      po[2,i,t,3] <- 1-p	
      po[3,i,t,1] <- 0	
      po[3,i,t,2] <- 0	
      po[3,i,t,3] <- 1	
      } #t	
   } #i	
	
# Likelihood	
for (i in 1:nind){	
   # Define latent state at first capture	
   z[i,f[i]] <- y[i,f[i]]	
   for (t in (f[i]+1):n.occasions){	
      # State process: draw S(t) given S(t-1)	
      z[i,t] ~ dcat(ps[z[i,t-1], i, t-1,])	
      # Observation process: draw O(t) given S(t)	
      y[i,t] ~ dcat(po[z[i,t], i, t-1,])	
      } #t	
   } #i	
}	
",fill = TRUE)	
sink()	
	
# Function to create known latent states z	
known.state.ms <- function(ms, notseen){	
   # notseen: label for ënot seení	
   state <- ms	
   state[state==notseen] <- NA	
   for (i in 1:dim(ms)[1]){	
      m <- min(which(!is.na(state[i,])))	
      state[i,m] <- NA	
      }	
   return(state)	
   }	
	
# Function to create initial values for unknown z	
ms.init.z <- function(ch, f){	
   for (i in 1:dim(ch)[1]){ch[i,1:f[i]] <- NA}	
   states <- max(ch, na.rm = TRUE)	
   known.states <- 1:(states-1)	
   v <- which(ch==states)	
   ch[-v] <- NA	
   ch[v] <- sample(known.states, length(v), replace = TRUE)	
   return(ch)	
   }	
	
# Bundle data	
jags.data <- list(y = his_recoded, f = f, quality = quality, n.occasions = dim(his_recoded)[2], nind = dim(his_recoded)[1], z = known.state.ms(his_recoded, 3))	
	
# Initial values	
inits <- function(){list(alphaA = rnorm(1, 0, 1), betaA = rnorm(1, 0, 1), psiAB = runif(1, 0, 1), p = runif(1, 0, 1), z = ms.init.z(his_recoded, f))}  	
	
# Parameters monitored	
parameters <- c("alphaA", "alphaB","betaA", "betaB", "psiAB", "psiBA", "p")	
	
# MCMC settings	
ni <- 5000	
nb <- 2000	
nc <- 2	
	
# Call JAGS from R	
library(R2jags)	
ms_quality <- jags(jags.data, inits, parameters, "state_on_survival_withquality.jags", n.chains = nc, n.iter = ni, n.burnin = nb)	
print(ms_quality, digits = 3)	
#' 	
#' 	
#' Same model without state effect on survival:	
#' 	
sink("no_state_on_survival_withquality.jags")	
cat("	
model {	
	
# -------------------------------------------------	
# Parameters:	
# alpha, beta: regression parameters for survival probability at site A/B	
# psiAB: movement probability from site A to site B	
# psiBA: movement probability from site B to site A	
# pA: recapture probability	
# -------------------------------------------------	
# States (S):	
# 1 alive at A	
# 2 alive at B	
# 3 dead	
# Observations (O):	
# 1 seen at A	
# 2 seen at B	
# 3 not seen	
# -------------------------------------------------	
	
# Priors	
   alpha ~ dnorm(0, 0.01)	
   beta ~ dnorm(0, 0.01)	
   for (i in 1:nind){logit(phi[i]) <- alpha + beta * quality[i]}	
   psiAB ~ dunif(0, 1)	
   psiBA ~ dunif(0, 1)	
   p ~ dunif(0, 1)	
	
# Define state-transition and observation matrices	
for (i in 1:nind){	
   # Define probabilities of state S(t+1) given S(t)	
   for (t in f[i]:(n.occasions-1)){	
      ps[1,i,t,1] <- phi[i] * (1-psiAB)	
      ps[1,i,t,2] <- phi[i] * psiAB	
      ps[1,i,t,3] <- 1-phi[i]	
      ps[2,i,t,1] <- phi[i] * psiBA	
      ps[2,i,t,2] <- phi[i] * (1-psiBA)	
      ps[2,i,t,3] <- 1-phi[i]	
      ps[3,i,t,1] <- 0	
      ps[3,i,t,2] <- 0	
      ps[3,i,t,3] <- 1	
	
      # Define probabilities of O(t) given S(t)	
      po[1,i,t,1] <- p	
      po[1,i,t,2] <- 0	
      po[1,i,t,3] <- 1-p	
      po[2,i,t,1] <- 0	
      po[2,i,t,2] <- p	
      po[2,i,t,3] <- 1-p	
      po[3,i,t,1] <- 0	
      po[3,i,t,2] <- 0	
      po[3,i,t,3] <- 1	
      } #t	
   } #i	
	
# Likelihood	
for (i in 1:nind){	
   # Define latent state at first capture	
   z[i,f[i]] <- y[i,f[i]]	
   for (t in (f[i]+1):n.occasions){	
      # State process: draw S(t) given S(t-1)	
      z[i,t] ~ dcat(ps[z[i,t-1], i, t-1,])	
      # Observation process: draw O(t) given S(t)	
      y[i,t] ~ dcat(po[z[i,t], i, t-1,])	
      } #t	
   } #i	
}	
",fill = TRUE)	
sink()	
	
# Function to create known latent states z	
known.state.ms <- function(ms, notseen){	
   # notseen: label for ënot seení	
   state <- ms	
   state[state==notseen] <- NA	
   for (i in 1:dim(ms)[1]){	
      m <- min(which(!is.na(state[i,])))	
      state[i,m] <- NA	
      }	
   return(state)	
   }	
	
# Function to create initial values for unknown z	
ms.init.z <- function(ch, f){	
   for (i in 1:dim(ch)[1]){ch[i,1:f[i]] <- NA}	
   states <- max(ch, na.rm = TRUE)	
   known.states <- 1:(states-1)	
   v <- which(ch==states)	
   ch[-v] <- NA	
   ch[v] <- sample(known.states, length(v), replace = TRUE)	
   return(ch)	
   }	
	
# Bundle data	
jags.data <- list(y = his_recoded, f = f, quality = quality, n.occasions = dim(his_recoded)[2], nind = dim(his_recoded)[1], z = known.state.ms(his_recoded, 3))	
	
# Initial values	
inits <- function(){list(alpha = rnorm(1, 0, 1), beta = rnorm(1, 0, 1), psiAB = runif(1, 0, 1), p = runif(1, 0, 1), z = ms.init.z(his_recoded, f))}  	
	
# Parameters monitored	
parameters <- c("alpha","beta", "psiAB", "psiBA", "p")	
	
# MCMC settings	
ni <- 5000	
nb <- 2000	
nc <- 2	
	
# Call JAGS from R	
library(R2jags)	
ms_quality_without <- jags(jags.data, inits, parameters, "no_state_on_survival_withquality.jags", n.chains = nc, n.iter = ni, n.burnin = nb)	
print(ms_quality_without, digits = 3)	
#' 	
#' 	
#' Clearly, the inclusion of quality improves the DIC. Also, the model with a difference in survival between breeders and non-breeders is better supported by the data when individual heterogeneity is accounted for.	
#' 	
#' # Models with individual random effects	
#' 	
#' Here, we aim at illustrating how not accounting for individual heterogeneity may obscure the detection of senescence in survival. 	
#' More specifically, we consider a single cohort of 500 individuals with survival decreasing as they age over a 20-year study. 	
#' We also add a frailty for each individual under the form of a normal distribution. Specifically, we specify $logit(\phi_i(a))=\beta_0+\beta_1 a + \varepsilon_i$ where $\varepsilon_i \sim N(0,\sigma^2)$. We use $\beta_0 = 1$, $\beta_1 = -0.05$ and $\sigma = 1$. If we condition upon the random effect, survival is decreasing as age increases. Note that we consider the same detection probability $p=0.5$ for all individuals.	
#' 	
#' ### Data simulation	
#' 	
#' First, we simulate survival for each individual then plot the individual trajectories (in grey) as well as survival conditional on the random effect (in red):	
#' 	
rm(list=ls())	
r = set.seed(3) # for reproducibility	
p = 0.5 # detection	
intercept_phi = 1 	
slope_phi = -0.05	
sigmaphi = 1	
nind = 500 # nb of individuals	
nyear = 20 # duration of the study	
expit<-function(x){exp(x)/(1+exp(x))} # reciprocal logit function	
z<-data<-x<-matrix(NA,nrow=nind,ncol=nyear)	
first<-rep(1,nind)	
age = matrix(NA,nind,nyear)	
phi = matrix(NA,nind,nyear)	
# simulate age-varying survival for each individual	
for (i in 1:nind){	
  mask <- first[i]:nyear	
  age[i,mask] <- mask - first[i] + 1	
  phi[i,mask] <- expit(intercept_phi + slope_phi * age[i,mask] + rnorm(1,0,sigmaphi))	
}	
plot(age[1,],phi[1,],type='l',col='grey',ylim=c(0,1),xlab='age',ylab='estimated survival')	
for (i in 2:nind){	
  lines(age[i,],phi[i,],type='l',col='grey')	
}	
lines(1:nyear,expit(intercept_phi + slope_phi * 1:nyear),col='red',lwd=2)	
#' 	
#' 	
#' Now simulate the encounter histories:	
#' 	
for(i in 1:nind){	
  z[i,first[i]] <- x[i,first[i]] <- 1	
  for(j in (first[i]+1):nyear){	
    z[i,j]<-rbinom(1,1,phi[i,j-1]*z[i,j-1])	
    x[i,j]<-rbinom(1,1,z[i,j]*p)	
  }	
}	
his = x	
his[is.na(his)]=0 # remove lines with 0's	
#' 	
#' 	
#' ### Model fitting	
#' 	
#' We fit the model with an age effect but no individual heterogeneity to the simulated dataset:	
#' 	
#' 	
sink("cjs_age.jags")	
cat("	
model {	
	
# Priors and constraints	
for (i in 1:nind){	
   for (t in f[i]:(n.occasions-1)){	
      logit(phi[i,t]) <- alpha + beta * t # a single cohort is used here, so that age = time elapsed since first capture = time	
      p[i,t] <- mean.p	
      } #t	
   } #i	
alpha ~ dnorm(0, 0.01)                   # Prior for intercept of age effect	
beta ~ dnorm(0, 0.01)                   # Prior for slope of age effect	
mean.p ~ dunif(0, 1)                     # Prior for mean recapture	
	
# Likelihood	
for (i in 1:nind){	
   # Define latent state at first capture	
   z[i,f[i]] <- 1	
   for (t in (f[i]+1):n.occasions){	
      # State process	
      z[i,t] ~ dbern(mu1[i,t])	
      mu1[i,t] <- phi[i,t-1] * z[i,t-1]	
      # Observation process	
      y[i,t] ~ dbern(mu2[i,t])	
      mu2[i,t] <- p[i,t-1] * z[i,t]	
      } #t	
   } #i	
}	
",fill = TRUE)	
sink()	
	
# Bundle data	
	
known.state.cjs <- function(ch){	
   state <- ch	
   for (i in 1:dim(ch)[1]){	
      n1 <- min(which(ch[i,]==1))	
      n2 <- max(which(ch[i,]==1))	
      state[i,n1:n2] <- 1	
      state[i,n1] <- NA	
      }	
   state[state==0] <- NA	
   return(state)	
   }	
	
jags.data <- list(y = his, f = first, nind = dim(his)[1], n.occasions = dim(his)[2], z = known.state.cjs(his))	
	
# Initial values 	
	
cjs.init.z <- function(ch,f){	
   for (i in 1:dim(ch)[1]){	
      if (sum(ch[i,])==1) next	
      n2 <- max(which(ch[i,]==1))	
      ch[i,f[i]:n2] <- NA	
      }	
   for (i in 1:dim(ch)[1]){	
   ch[i,1:f[i]] <- NA	
   }	
   return(ch)	
   }	
	
inits <- function(){list(z = cjs.init.z(his, first), alpha = rnorm(1, 0, 1), mean.p = runif(1, 0, 1), beta = rnorm(1, 0, 1))}  	
	
# Parameters monitored	
parameters <- c("alpha", "beta", "mean.p")	
	
# MCMC settings	
ni <- 5000	
nb <- 2000	
nc <- 2	
	
# Call JAGS from R	
library(R2jags)	
cjs_age <- jags(jags.data, inits, parameters, "cjs_age.jags", n.chains = nc, n.iter = ni, n.burnin = nb)	
	
# Summarize posteriors	
print(cjs_age, digits = 3)	
#' 	
#' 	
#' Having a look to the parameter estimates, it sounds like the slope of the age effect on survival is estimated positive... 	
#' 	
#' Which means that at the population level, whenever individual heterogeneity is ignored, then senescence (in red) is completely masked. Even worse, survival is increasing with increasing age (in blue).	
#' 	
plot(age[1,],phi[1,],type='l',col='grey',ylim=c(0,1),xlab='age',ylab='estimated survival')	
for (i in 2:nind){	
  lines(age[i,],phi[i,],type='l',col='grey')	
}	
lines(1:nyear,expit(intercept_phi + slope_phi * 1:nyear),col='red',lwd=2)	
lines(1:nyear,expit(mean(cjs_age$BUGSoutput$sims.matrix[,'alpha']) + mean(cjs_age$BUGSoutput$sims.matrix[,'beta']) * 1:nyear),col='blue',lwd=2)	
#' 	
#' 	
#' Now we fit the model with a random effect in the survival process. 	
#' 	
#' 	
sink("cjs_age_re.jags")	
cat("	
model {	
	
# Priors and constraints	
for (i in 1:nind){	
   for (t in f[i]:(n.occasions-1)){	
      logit(phi[i,t]) <- alpha + beta * t + epsilon[i] # a single cohort is used here, so that age = time elapsed since first capture = time	
      p[i,t] <- mean.p	
      } #t	
   } #i	
for (i in 1:nind){	
   epsilon[i] ~ dnorm(0, tau)	
   }	
alpha ~ dnorm(0, 0.01)                   # Prior for intercept of age effect	
beta ~ dnorm(0, 0.01)                   # Prior for slope of age effect	
sigma ~ dunif(0, 5)                      # Prior for standard deviation	
tau <- pow(sigma, -2)	
sigma2 <- pow(sigma, 2)	
mean.p ~ dunif(0, 1)                     # Prior for mean recapture	
	
# Likelihood	
for (i in 1:nind){	
   # Define latent state at first capture	
   z[i,f[i]] <- 1	
   for (t in (f[i]+1):n.occasions){	
      # State process	
      z[i,t] ~ dbern(mu1[i,t])	
      mu1[i,t] <- phi[i,t-1] * z[i,t-1]	
      # Observation process	
      y[i,t] ~ dbern(mu2[i,t])	
      mu2[i,t] <- p[i,t-1] * z[i,t]	
      } #t	
   } #i	
}	
",fill = TRUE)	
sink()	
	
# Bundle data	
jags.data <- list(y = his, f = first, nind = dim(his)[1], n.occasions = dim(his)[2], z = known.state.cjs(his))	
	
# Initial values 	
inits <- function(){list(z = cjs.init.z(his, first), alpha = rnorm(1, 0, 1),  beta = rnorm(1, 0, 1),mean.p = runif(1, 0, 1), sigma = runif(1, 0, 2))}  	
	
# Parameters monitored	
parameters <- c("alpha", "beta", "mean.p", "sigma2")	
	
# MCMC settings	
ni <- 5000	
nb <- 2000	
nc <- 2	
	
# Call JAGS from R	
cjs_age_re <- jags(jags.data, inits, parameters, "cjs_age_re.jags", n.chains = nc, n.iter = ni, n.burnin = nb)	
	
# Summarize posteriors	
print(cjs_age_re, digits = 3)	
#' 	
#' 	
#' The intercept and slope of the age-survival relationship are quite close to the values we used to simulate the data.	
#' 	
#' Now we add to our previous plot the survival as estimated when individual heterogeneity is explicitely accounted for using individual random effects (in green):	
#' 	
plot(age[1,],phi[1,],type='l',col='grey',ylim=c(0,1),xlab='age',ylab='estimated survival')	
for (i in 2:nind){	
  lines(age[i,],phi[i,],type='l',col='grey')	
}	
lines(1:nyear,expit(intercept_phi + slope_phi * 1:nyear),col='red',lwd=2)	
lines(1:nyear,expit(mean(cjs_age$BUGSoutput$sims.matrix[,'alpha']) + mean(cjs_age$BUGSoutput$sims.matrix[,'beta']) * 1:nyear),col='blue',lwd=2)	
lines(1:nyear,expit(mean(cjs_age_re$BUGSoutput$sims.matrix[,'alpha']) + mean(cjs_age_re$BUGSoutput$sims.matrix[,'beta']) * 1:nyear),col='green',lwd=2)	
#' 	
#' 	
#' ## Models with finite mixtures	
#' 	
#' Here, we again aim at illustrating how not accounting for individual heterogeneity may obscure the detection of senescence in survival. In contrast with the previous section, we now use finite mixtures to deal with heterogeneity. More specifically, we consider a cohort of 1000 individuals that are split into a group of robust individuals in proportion $\pi$ with constant high survival $\phi_R$ and a group of frail individuals with survival $\phi_F$ that senesce over the 20 years of the study according to the relationship $logit(\phi_F(a))=\beta_0+\beta_1 a$. We use $\pi = 0.3$, $\phi_R = 0.85$, $\beta_0 = 0$ and $\beta_1 = -0.07$. Note that we consider the same detection probability $p=0.5$ for all individuals.	
#' 	
#' ### Data simulation	
#' 	
#' First simulate data	
#' 	
rm(list=ls())	
r = set.seed(3) # for reproducibility	
p = 0.5 # detection	
prop_class1 = 0.3 # pi	
phi_class1 = 0.85 # survival or robust ind	
intercept_phi_class2 = 0 #beta_0	
slope_phi_class2 = -0.05 # beta_1	
nind = 1000 # nb of ind	
nyear = 20 # duration of the study	
expit<-function(x){exp(x)/(1+exp(x))} # reciprocal of the logit function	
z<-data<-x<-matrix(NA,nrow=nind,ncol=nyear)	
first<-rep(1,nind)	
age = matrix(NA,nind,nyear)	
phi = matrix(NA,nind,nyear)	
which_mixture = rep(NA,nind)	
# simulate age-varying survival for each individual, 	
# by first assigning them to the robust or frail class, then using the corresponding 	
# survival 	
for (i in 1:nind){	
  mask <- first[i]:nyear	
  age[i,mask] <- mask - first[i] + 1	
  which_mixture[i] <- rbinom(1,1,prop_class1) # assign ind i to a class with prob pi	
  if (which_mixture[i] == 1){	
    phi[i,mask] <- phi_class1 # robust	
  } else { 	
  phi[i,mask] <- expit(intercept_phi_class2 + slope_phi_class2 * age[i,mask])} # frail	
}	
#' 	
#' 	
#' Represent graphically survival over time in the two classes:	
#' 	
plot(age[1,],phi[1,],type='l',col='grey',ylim=c(0,1),xlab='age',ylab='estimated survival')	
for (i in 2:nind){	
  lines(age[i,],phi[i,],type='l',col='grey')	
}	
#' 	
#' 	
#' Now simulate the encounter histories:	
#' 	
for(i in 1:nind){	
  z[i,first[i]] <- x[i,first[i]] <- 1	
  for(j in (first[i]+1):nyear){	
    z[i,j]<-rbinom(1,1,phi[i,j-1]*z[i,j-1])	
    x[i,j]<-rbinom(1,1,z[i,j]*p)	
  }	
}	
his = x	
his[is.na(his)]=0	
#' 	
#' 	
#' ### Model fitting	
#' 	
#' Let's fit two models assuming homogeneity, first one with constant survival probability, second one with an age effect: 	
#' 	
#' 	
sink("Mct.jags")	
cat("	
model	
{	
	
# notation used	
# STATES	
# V for alive	
# D for dead	
	
# OBSERVATIONS	
# 0 = non-observed (coded 1)	
# 1 = observed (coded 2)	
	
# prior on survival for V1 and V2	
phi ~ dunif(0,1)	
	
# prior on detection	
p ~ dunif(0,1)	
	
# probabilities for each initial state	
px0[1] <- 1 # prob. of being in initial state V	
px0[2] <- 0 # prob. of being in initial state D	
	
# define probabilities of observations at t given states at t	
po[1,1] <- 1-p	
po[1,2] <- p	
po[2,1] <- 1	
po[2,2] <- 0	
	
po.init[1,1] <- 0	
po.init[1,2] <- 1	
po.init[2,1] <- 1	
po.init[2,2] <- 0	
	
# define probabilities of states at t given states at t-1	
px[1,1] <- phi	
px[1,2] <- 1-phi	
px[2,1] <- 0	
px[2,2] <- 1	
	
for (i in 1:nind){  # for each indiv	
  # initial states	
  z[i,first[i]] ~ dcat(px0[1:2])	
  x[i,first[i]] ~ dcat(po.init[z[i,first[i]],1:2])	
  for (j in (first[i]+1):nyear){  # loop over time	
	
    # state equations	
    z[i,j] ~ dcat(px[z[i,j-1],1:2])	
	
	  # observation equations	
	  x[i,j] ~ dcat(po[z[i,j],1:2])	
    }	
}	
	
}	
	
",fill=TRUE)	
sink()	
	
mydatax <- list(x=his+1,first=first,nind=nind,nyear=nyear)	
	
# initial values	
x.init <- his	
for (i in 1:nind){x.init[i,1:first[i]] <- NA}	
x.init[x.init==0] <- 1	
z = x.init	
init1 <- list(phi=0.3,p=.2,z=z)	
init2 <- list(phi=0.7,p=.8,z=z)	
inits <- list(init1,init2)	
	
	
# parameters to be monitored	
parameters <- c("phi","p")	
	
# MCMC settings	
ni <- 5000	
nb <- 2000	
nc <- 2	
	
# Call JAGS from R	
library(R2jags)	
start<-as.POSIXlt(Sys.time())	
no_het_no_age <- jags(mydatax,inits,parameters,"Mct.jags",n.chains = nc, n.iter = ni,n.burnin = nb)	
end <-as.POSIXlt(Sys.time())	
duration = end-start	
	
# Summarize posteriors	
print(no_het_no_age, digits = 3)	
#' 	
#' 	
#' 	
sink("Mage.jags")	
cat("	
model	
{	
	
# notation used	
# STATES	
# V for alive	
# D for dead	
	
# OBSERVATIONS	
# 0 = non-observed (coded 1)	
# 1 = observed (coded 2)	
	
# prior on survival	
alpha ~ dnorm(0,0.01)	
beta ~ dnorm(0,0.01)	
for (i in 1:nind){  # for each indiv	
  for (j in 1:(nyear-1)){  # loop over time	
    logit(phi[i,j]) <- alpha + beta * j # single cohort, hence age = time	
  }	
}	
# prior on detection	
p ~ dunif(0,1)	
	
# probabilities for each initial state	
px0[1] <- 1 # prob. of being in initial state V	
px0[2] <- 0 # prob. of being in initial state D	
	
# define probabilities of observations at t given states at t	
po[1,1] <- 1-p	
po[1,2] <- p	
po[2,1] <- 1	
po[2,2] <- 0	
	
po.init[1,1] <- 0	
po.init[1,2] <- 1	
po.init[2,1] <- 1	
po.init[2,2] <- 0	
	
# define probabilities of states at t given states at t-1	
for (i in 1:nind){  # for each indiv	
  for (j in 1:(nyear-1)){  # loop over time	
    px[1,i,j,1] <- phi[i,j]	
    px[1,i,j,2] <- 1-phi[i,j]	
    px[2,i,j,1] <- 0	
    px[2,i,j,2] <- 1	
  }	
}	
for (i in 1:nind){  # for each indiv	
  # initial states	
  z[i,first[i]] ~ dcat(px0[1:2])	
  x[i,first[i]] ~ dcat(po.init[z[i,first[i]],1:2])	
  for (j in (first[i]+1):nyear){  # loop over time	
	
    # state equations	
    z[i,j] ~ dcat(px[z[i,j-1],i,j-1,1:2])	
	
	  # observation equations	
	  x[i,j] ~ dcat(po[z[i,j],1:2])	
    }	
}	
	
}	
	
",fill=TRUE)	
sink()	
	
mydatax <- list(x=his+1,first=first,nind=nind,nyear=nyear)	
	
# initial values	
x.init <- his	
for (i in 1:nind){x.init[i,1:first[i]] <- NA}	
x.init[x.init==0] <- 1	
z = x.init	
init1 <- list(alpha=0.3,p=.2,z=z)	
init2 <- list(alpha=0.7,p=.8,z=z)	
inits <- list(init1,init2)	
	
	
# parameters to be monitored	
parameters <- c("alpha","beta","p")	
	
# MCMC settings	
ni <- 5000	
nb <- 2000	
nc <- 2	
	
# Call JAGS from R	
library(R2jags)	
start<-as.POSIXlt(Sys.time())	
no_het_with_age <- jags(mydatax,inits,parameters,"Mage.jags",n.chains = nc, n.iter = ni,n.burnin = nb)	
end <-as.POSIXlt(Sys.time())	
duration = end-start	
	
# Summarize posteriors	
print(no_het_with_age, digits = 3)	
#' 	
#' 	
#' Graphically, we have the estimate from the model with constant survival (in red) vs. age-varying survival (in blue):	
#' 	
plot(age[1,],phi[1,],type='l',col='grey',ylim=c(0,1),xlab='age',ylab='estimated survival')	
for (i in 2:nind){	
  lines(age[i,],phi[i,],type='l',col='grey')	
}	
lines(1:nyear,rep(mean(no_het_no_age$BUGSoutput$sims.matrix[,'phi']),nyear),lwd=2,col='red') # add survival from constant model	
lines(1:nyear,expit(mean(no_het_with_age$BUGSoutput$sims.matrix[,'alpha'])+mean(no_het_with_age$BUGSoutput$sims.matrix[,'beta'])*(1:nyear)),lwd=2,col='blue') # add survival from age model	
#' 	
#' 	
#' Again, as in the previous section, it's striking to see that survival is increasing when age increases if individual heterogeneity is ignored. In other words, senescence is masked. 	
#' 	
#' Now let's fit a model with heterogeneity in the survival probability, with constant parameters over time. Constant survival, two-finite mixture on the survival probability and a constant proportion of individual in each class:	
#' 	
#' 	
sink("M2class.jags")	
cat("	
model	
{	
	
# notation used	
# STATES	
# V1 for alive in class 1	
# V2 for alive in class 2	
# D for dead	
	
# OBSERVATIONS	
# 0 = non-observed (coded 1)	
# 1 = observed (coded 2)	
	
# prior on survival for V1 and V2	
phi1 ~ dunif(0,1)	
phi2 ~ dunif(0,1)	
	
# prior on detection	
p ~ dunif(0,1)	
	
# prior on initial state prob	
pi ~ dunif(0,1)	
	
# probabilities for each initial state	
px0[1] <- pi # prob. of being in initial state V1	
px0[2] <- 1-pi # prob. of being in initial state V2	
px0[3] <- 0 # prob. of being in initial state dead	
	
# define probabilities of observations at t given states at t	
po[1,1] <- 1-p	
po[1,2] <- p	
po[2,1] <- 1-p	
po[2,2] <- p	
po[3,1] <- 1	
po[3,2] <- 0	
	
po.init[1,1] <- 0	
po.init[1,2] <- 1	
po.init[2,1] <- 0	
po.init[2,2] <- 1	
po.init[3,1] <- 1	
po.init[3,2] <- 0	
	
# define probabilities of states at t given states at t-1	
px[1,1] <- phi1	
px[1,2] <- 0	
px[1,3] <- 1-phi1	
px[2,1] <- 0	
px[2,2] <- phi2	
px[2,3] <- 1-phi2	
px[3,1] <- 0	
px[3,2] <- 0	
px[3,3] <- 1	
	
for (i in 1:nind){  # for each indiv	
  # initial states	
  z[i,first[i]] ~ dcat(px0[1:3])	
  x[i,first[i]] ~ dcat(po.init[z[i,first[i]],1:2])	
  for (j in (first[i]+1):nyear){  # loop over time	
	
    # state equations	
    z[i,j] ~ dcat(px[z[i,j-1],1:3])	
	
	  # observation equations	
	  x[i,j] ~ dcat(po[z[i,j],1:2])	
    }	
}	
	
}	
	
",fill=TRUE)	
sink()	
	
mydatax <- list(x=his+1,first=first,nind=nind,nyear=nyear)	
	
# initial values	
x.init <- his	
for (i in 1:nind){x.init[i,1:first[i]] <- NA}	
x.init[x.init==0] <- 1	
z = x.init	
init1 <- list(phi1=0.3,phi2=.2,p=.2,z=z)	
init2 <- list(phi1=0.7,phi2=.8,p=.8,z=z)	
inits <- list(init1,init2)	
	
# load package R2jags	
library(R2jags)	
	
# parameters to be monitored	
parameters <- c("phi1","phi2","p")	
	
# MCMC settings	
ni <- 5000	
nb <- 2000	
nc <- 2	
	
# Call JAGS from R	
start<-as.POSIXlt(Sys.time())	
with_het_no_age <- jags(mydatax,inits,parameters,"M2class.jags",n.chains = nc, n.iter = ni,n.burnin = nb)	
end <-as.POSIXlt(Sys.time())	
duration = end-start	
	
# Summarize posteriors	
print(with_het_no_age, digits = 3)	
#' 	
#' 	
#' Let's have a look graphically:	
#' 	
plot(age[1,],phi[1,],type='l',col='grey',ylim=c(0,1),xlab='age',ylab='estimated survival')	
for (i in 2:nind){	
  lines(age[i,],phi[i,],type='l',col='grey')	
}	
lines(1:nyear,rep(mean(with_het_no_age$BUGSoutput$sims.matrix[,'phi1']),nyear),lwd=2,col='green') # add survival from first class	
lines(1:nyear,rep(mean(with_het_no_age$BUGSoutput$sims.matrix[,'phi2']),nyear),lwd=2,col='green') # add survival from second class	
#' 	
#' 	
#' Not too bad. Obviously, for frail individuals, we miss the age effect to be able to detect senescence. Now let's add age to this model:	
#' 	
#' 	
sink("M2classage.jags")	
cat("	
model	
{	
	
# notation used	
# STATES	
# V1 for alive in class 1	
# V2 for alive in class 2	
# D for dead	
	
# OBSERVATIONS	
# 0 = non-observed (coded 1)	
# 1 = observed (coded 2)	
	
# prior on survival for V1 and V2	
alpha1 ~ dnorm(0,0.01)	
beta1 ~ dnorm(0,0.01)	
alpha2 ~ dnorm(0,0.01)	
beta2 ~ dnorm(0,0.01)	
for (i in 1:nind){  # for each indiv	
  for (j in 1:(nyear-1)){  # loop over time	
    logit(phi1[i,j]) <- alpha1 + beta1 * j # single cohort, hence age = time	
    logit(phi2[i,j]) <- alpha2 + beta2 * j # single cohort, hence age = time	
  }	
}	
	
# prior on detection	
p ~ dunif(0,1)	
	
# prior on initial state prob	
pi ~ dunif(0,1)	
	
# probabilities for each initial state	
px0[1] <- pi # prob. of being in initial state V1	
px0[2] <- 1-pi # prob. of being in initial state V2	
px0[3] <- 0 # prob. of being in initial state dead	
	
# define probabilities of observations at t given states at t	
po[1,1] <- 1-p	
po[1,2] <- p	
po[2,1] <- 1-p	
po[2,2] <- p	
po[3,1] <- 1	
po[3,2] <- 0	
	
po.init[1,1] <- 0	
po.init[1,2] <- 1	
po.init[2,1] <- 0	
po.init[2,2] <- 1	
po.init[3,1] <- 1	
po.init[3,2] <- 0	
	
# define probabilities of states at t given states at t-1	
for (i in 1:nind){  # for each indiv	
  for (j in 1:(nyear-1)){  # loop over time	
    px[1,i,j,1] <- phi1[i,j]	
    px[1,i,j,2] <- 0	
    px[1,i,j,3] <- 1-phi1[i,j]	
    px[2,i,j,1] <- 0	
    px[2,i,j,2] <- phi2[i,j]	
    px[2,i,j,3] <- 1-phi2[i,j]	
    px[3,i,j,1] <- 0	
    px[3,i,j,2] <- 0	
    px[3,i,j,3] <- 1	
  }	
}	
	
for (i in 1:nind){  # for each indiv	
  # initial states	
  z[i,first[i]] ~ dcat(px0[1:3])	
  x[i,first[i]] ~ dcat(po.init[z[i,first[i]],1:2])	
  for (j in (first[i]+1):nyear){  # loop over time	
	
    # state equations	
    z[i,j] ~ dcat(px[z[i,j-1],i,j-1,1:3])	
	
	  # observation equations	
	  x[i,j] ~ dcat(po[z[i,j],1:2])	
    }	
}	
	
}	
	
",fill=TRUE)	
sink()	
	
mydatax <- list(x=his+1,first=first,nind=nind,nyear=nyear)	
	
# initial values	
x.init <- his	
for (i in 1:nind){x.init[i,1:first[i]] <- NA}	
x.init[x.init==0] <- 1	
z = x.init	
init1 <- list(alpha1=1,alpha2=-1,p=.2,z=z)	
init2 <- list(alpha1=-1,alpha2=1,p=.8,z=z)	
inits <- list(init1,init2)	
	
# load package R2jags	
library(R2jags)	
	
# parameters to be monitored	
parameters <- c("alpha1","alpha2","beta1","beta2","p")	
	
# MCMC settings	
ni <- 5000	
nb <- 2000	
nc <- 2	
	
# Call JAGS from R	
start<-as.POSIXlt(Sys.time())	
with_het_with_age <- jags(mydatax,inits,parameters,"M2classage.jags",n.chains = nc, n.iter = ni,n.burnin = nb)	
end <-as.POSIXlt(Sys.time())	
duration = end-start	
	
# Summarize posteriors	
print(with_het_with_age, digits = 3)	
#' 	
#' 	
#' Let's have a look graphically:	
#' 	
plot(age[1,],phi[1,],type='l',col='grey',ylim=c(0,1),xlab='age',ylab='estimated survival')	
for (i in 2:nind){	
  lines(age[i,],phi[i,],type='l',col='grey')	
}	
phi1 = 1/(1+exp(-mean(with_het_with_age$BUGSoutput$sims.matrix[,'alpha1']) + mean(with_het_with_age$BUGSoutput$sims.matrix[,'beta1']) * 1:nyear))	
phi2 = 1/(1+exp(-mean(with_het_with_age$BUGSoutput$sims.matrix[,'alpha2']) + mean(with_het_with_age$BUGSoutput$sims.matrix[,'beta2']) * 1:nyear))	
lines(1:nyear,phi1,lwd=2,col='green') # add survival from first class	
lines(1:nyear,phi2,lwd=2,col='green') # add survival from second class	
lines(1:nyear,expit(mean(no_het_with_age$BUGSoutput$sims.matrix[,'alpha'])+mean(no_het_with_age$BUGSoutput$sims.matrix[,'beta'])*(1:nyear)),lwd=2,col='blue') # add survival from age model	
#' 	
#' 	
#' Seems like we've managed to capture the main patterns in the simulated data.