appendix_paper_IHCMR_partB.Rmd

---
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:
```{r}
# 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:
```{r}
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:
```{r}
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:
```{r}
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:
```{r}
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:
```{r}
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()
```

```{r}
# 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:
```{r}
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()
```

```{r}
# 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:
```{r}
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.

```{r}
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:
```{r}
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):
```{r}
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:
```{r}
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:

```{r}
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).
```{r}
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. 

```{r}
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):
```{r}
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
```{r}
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:
```{r}
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:
```{r}
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: 

```{r}
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)
```

```{r}
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):
```{r}
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:

```{r}
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:
```{r}
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:

```{r}
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:
```{r}
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.