Introduction

In this demo, we illustrate how to fit multistate models to capture-recapture data. We first consider states as sites and estimate movements. Next we treat states as breeding states for studying life-history trade-offs.

library(tidyverse)
library(nimble)
library(MCMCvis)

Multisite model with 2 sites

We are going to analyze real capture-recapture data on > 20,000 Canada geese (Branta canadensis) marked and resighted in the east coast of the US in the mid-Atlantic (New York, Pennsylvania, New Jersey), Chesapeake (Delaware, Maryland, Virginia), and Carolinas (North and South Carolina). The data were kindly provided by Jay Hestbeck. We analyse a subset of 500 geese from the original dataset, we’ll get back to the whole dataset in the last live demo.

Read in the data.

geese <- read_csv("geese.csv", col_names = TRUE)
y <- as.matrix(geese)
head(y)
     year_1984 year_1985 year_1986 year_1987 year_1988 year_1989
[1,]         0         2         2         0         0         0
[2,]         0         0         0         0         0         2
[3,]         0         0         0         1         0         0
[4,]         0         0         2         0         0         0
[5,]         0         3         0         0         3         2
[6,]         0         0         0         2         0         0

For simplicity, we start by considering only 2 sites, dropping Carolinas.

y2 <- y
y2[y2==3] <- 0 # act as if there was no detection in site 3 Carolinas
mask <- apply(y2, 1, sum)
y2 <- y2[mask!=0,] # remove rows w/ 0s only
head(y2)
     year_1984 year_1985 year_1986 year_1987 year_1988 year_1989
[1,]         0         2         2         0         0         0
[2,]         0         0         0         0         0         2
[3,]         0         0         0         1         0         0
[4,]         0         0         2         0         0         0
[5,]         0         0         0         0         0         2
[6,]         0         0         0         2         0         0

Get occasion of first capture.

get.first <- function(x) min(which(x != 0))
first <- apply(y2, 1, get.first)
first
  [1] 2 6 4 3 6 4 1 5 4 2 2 3 1 5 4 6 3 4 5 2 5 5 4 3 4 4 4 3 4 1 2 4 3 2 2 4 5
 [38] 2 3 2 2 2 1 1 5 5 2 3 3 2 3 3 4 4 4 6 2 3 1 2 2 3 2 2 4 4 3 2 4 3 3 3 4 2
 [75] 2 3 2 4 5 2 3 2 5 4 2 3 3 2 3 1 2 3 4 2 6 3 2 1 4 5 3 1 3 3 5 1 1 6 2 1 4
[112] 2 2 3 2 4 1 5 4 4 5 5 3 1 4 3 5 3 2 5 3 3 3 2 2 4 5 2 3 3 1 5 5 2 3 1 6 3
[149] 3 5 2 3 2 2 3 3 3 1 4 3 2 1 5 3 2 1 3 3 1 4 3 2 4 4 4 2 1 4 1 1 3 3 3 3 3
[186] 3 3 3 1 3 2 3 3 3 4 3 2 1 3 4 4 2 2 2 2 3 3 2 2 5 4 3 1 2 2 3 6 3 3 6 2 1
[223] 4 2 3 6 2 3 5 2 5 1 3 5 6 1 2 3 1 3 5 3 1 3 5 4 2 2 5 3 3 4 1 5 4 2 5 3 6
[260] 4 1 4 2 2 2 4 3 1 3 5 3 2 1 3 1 1 5 3 3 3 3 1 6 4 1 1 5 5 3 2 2 6 4 5 2 4
[297] 2 4 1 3 4 2 4 3 2 4 1 4 1 3 2 1 2 2 2 4 4 5 4 3 3 3 4 4 4 3 5 3 3 2 3 4 3
[334] 4 5 3 6 2 2 2 2 3 2 2 5 3 2 3 3 4 4 2 5 2 3 3 4 5 5 5 4 2 4 2 4 1 3 2 3 6
[371] 2 2 4 3 1 2 2 1 1 3 3 3 3 6 3 1 2 5 1 4 1 2 2 3 4 1 2 3 4 1 2 4 2 5 2 5 2
[408] 5 4 2 6 2 1 1 2 3 4

Write the model.

multisite <- nimbleCode({
  
  # -------------------------------------------------
  # Parameters:
  # phiA: survival probability site A
  # phiB: survival probability site B
  # psiAB: movement probability from site A to site B
  # psiBA: movement probability from site B to site A
  # pA: recapture probability site A
  # pB: recapture probability site B
  # piA: prop of being in site A at initial capture
  # -------------------------------------------------
  # States (z):
  # 1 alive at A
  # 2 alive at B
  # 3 dead
  # Observations (y):  
  # 1 not seen
  # 2 seen at site A 
  # 3 seen at site B
  # -------------------------------------------------
  
  # priors
  phiA ~ dunif(0, 1)
  phiB ~ dunif(0, 1)
  psiAB ~ dunif(0, 1)
  psiBA ~ dunif(0, 1)
  pA ~ dunif(0, 1)
  pB ~ dunif(0, 1)
  piA ~ dunif(0, 1)
  
  # probabilities of state z(t+1) given z(t)
  gamma[1,1] <- phiA * (1 - psiAB)
  gamma[1,2] <- phiA * psiAB
  gamma[1,3] <- 1 - phiA
  gamma[2,1] <- phiB * psiBA
  gamma[2,2] <- phiB * (1 - psiBA)
  gamma[2,3] <- 1 - phiB
  gamma[3,1] <- 0
  gamma[3,2] <- 0
  gamma[3,3] <- 1
  
  delta[1] <- piA          # Pr(alive in A t = 1)
  delta[2] <- 1 - piA      # Pr(alive in B t = 1)
  delta[3] <- 0            # Pr(dead t = 1) = 0
  
  # probabilities of y(t) given z(t)
  omega[1,1] <- 1 - pA     # Pr(alive A t -> non-detected t)
  omega[1,2] <- pA         # Pr(alive A t -> detected A t)
  omega[1,3] <- 0          # Pr(alive A t -> detected B t)
  omega[2,1] <- 1 - pB     # Pr(alive B t -> non-detected t)
  omega[2,2] <- 0          # Pr(alive B t -> detected A t)
  omega[2,3] <- pB         # Pr(alive B t -> detected B t)
  omega[3,1] <- 1          # Pr(dead t -> non-detected t)
  omega[3,2] <- 0          # Pr(dead t -> detected A t)
  omega[3,3] <- 0          # Pr(dead t -> detected B t)
  
  # likelihood 
  for (i in 1:N){
    # latent state at first capture
    z[i,first[i]] ~ dcat(delta[1:3])
    for (t in (first[i]+1):K){
      # draw z(t) given z(t-1)
      z[i,t] ~ dcat(gamma[z[i,t-1],1:3])
      # draw y(t) given z(t)
      y[i,t] ~ dcat(omega[z[i,t],1:3])
    }
  }
})

Data in a list. Remember to add 1.

my.data <- list(y = y2 + 1)

Constants in a list.

my.constants <- list(first = first, 
                     K = ncol(y2), 
                     N = nrow(y2))

Initial values.

zinits <- y2 # say states are observations, detections in A or B are taken as alive in same sites 
zinits[zinits==0] <- sample(c(1,2), sum(zinits==0), replace = TRUE) # non-detections become alive in site A or B
initial.values <- function(){list(phiA = runif(1, 0, 1), 
                                  phiB = runif(1, 0, 1), 
                                  psiAB = runif(1, 0, 1), 
                                  psiBA = runif(1, 0, 1), 
                                  pA = runif(1, 0, 1), 
                                  pB = runif(1, 0, 1), 
                                  piA = runif(1, 0, 1),
                                  z = zinits)}

Parameters to monitor.

parameters.to.save <- c("phiA", "phiB","psiAB", "psiBA", "pA", "pB", "piA")
parameters.to.save
[1] "phiA"  "phiB"  "psiAB" "psiBA" "pA"    "pB"    "piA"  

MCMC settings.

n.iter <- 10000
n.burnin <- 5000
n.chains <- 2

Run nimble.

mcmc.multisite <- nimbleMCMC(code = multisite, 
                             constants = my.constants,
                             data = my.data,              
                             inits = initial.values,
                             monitors = parameters.to.save,
                             niter = n.iter,
                             nburnin = n.burnin, 
                             nchains = n.chains)
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|

Let’s inspect the results. First the numerical summaries.

MCMCsummary(mcmc.multisite, round = 2)
      mean   sd 2.5%  50% 97.5% Rhat n.eff
pA    0.56 0.11 0.37 0.55  0.78 1.01   226
pB    0.39 0.04 0.30 0.39  0.48 1.00   207
phiA  0.56 0.08 0.41 0.55  0.72 1.01   208
phiB  0.73 0.05 0.64 0.73  0.82 1.00   180
piA   0.38 0.10 0.20 0.38  0.58 1.00    85
psiAB 0.23 0.10 0.08 0.22  0.47 1.01   191
psiBA 0.06 0.02 0.02 0.06  0.11 1.02   410

Second, a caterpillar plot of the parameter estimates.

MCMCplot(mcmc.multisite)

Actually, the initial state of a goose is known exactly. It is alive at site of initial capture. Therefore, we don’t need to try and estimate the probability of initial states. Let’s rewrite the model.

multisite <- nimbleCode({
  
  # -------------------------------------------------
  # Parameters:
  # phiA: survival probability site A
  # phiB: survival probability site B
  # psiAB: movement probability from site A to site B
  # psiBA: movement probability from site B to site A
  # pA: recapture probability site A
  # pB: recapture probability site B
  # -------------------------------------------------
  # States (z):
  # 1 alive at A
  # 2 alive at B
  # 3 dead
  # Observations (y):  
  # 1 not seen
  # 2 seen at A 
  # 3 seen at B
  # -------------------------------------------------
  
  # priors
  phiA ~ dunif(0, 1)
  phiB ~ dunif(0, 1)
  psiAB ~ dunif(0, 1)
  psiBA ~ dunif(0, 1)
  pA ~ dunif(0, 1)
  pB ~ dunif(0, 1)

  # probabilities of state z(t+1) given z(t)
  gamma[1,1] <- phiA * (1 - psiAB)
  gamma[1,2] <- phiA * psiAB
  gamma[1,3] <- 1 - phiA
  gamma[2,1] <- phiB * psiBA
  gamma[2,2] <- phiB * (1 - psiBA)
  gamma[2,3] <- 1 - phiB
  gamma[3,1] <- 0
  gamma[3,2] <- 0
  gamma[3,3] <- 1
  
  # probabilities of y(t) given z(t)
  omega[1,1] <- 1 - pA     # Pr(alive A t -> non-detected t)
  omega[1,2] <- pA         # Pr(alive A t -> detected A t)
  omega[1,3] <- 0          # Pr(alive A t -> detected B t)
  omega[2,1] <- 1 - pB     # Pr(alive B t -> non-detected t)
  omega[2,2] <- 0          # Pr(alive B t -> detected A t)
  omega[2,3] <- pB         # Pr(alive B t -> detected B t)
  omega[3,1] <- 1          # Pr(dead t -> non-detected t)
  omega[3,2] <- 0          # Pr(dead t -> detected A t)
  omega[3,3] <- 0          # Pr(dead t -> detected B t)
  
  # likelihood 
  for (i in 1:N){
    # latent state at first capture
    z[i,first[i]] <- y[i,first[i]] - 1 # if seen at site A (y = 2), state is alive in A (y - 1 = z = 1) with prob = 1
    for (t in (first[i]+1):K){         # if seen at site B (y = 3), state is alive in A (y - 1 = z = 2) with prob = 1
      # draw (t) given z(t-1)
      z[i,t] ~ dcat(gamma[z[i,t-1],1:3])
      # draw y(t) given z(t)
      y[i,t] ~ dcat(omega[z[i,t],1:3])
    }
  }
})

Initial values without \(p_A\).

zinits <- y2
zinits[zinits==0] <- sample(c(1,2), sum(zinits==0), replace = TRUE)
initial.values <- function(){list(phiA = runif(1, 0, 1), 
                                  phiB = runif(1, 0, 1), 
                                  psiAB = runif(1, 0, 1), 
                                  psiBA = runif(1, 0, 1), 
                                  pA = runif(1, 0, 1), 
                                  pB = runif(1, 0, 1), 
                                  z = zinits)}

Parameters to monitor.

parameters.to.save <- c("phiA", "phiB","psiAB", "psiBA", "pA", "pB")
parameters.to.save
[1] "phiA"  "phiB"  "psiAB" "psiBA" "pA"    "pB"   

Run nimble.

mcmc.multisite <- nimbleMCMC(code = multisite, 
                             constants = my.constants,
                             data = my.data,              
                             inits = initial.values,
                             monitors = parameters.to.save,
                             niter = n.iter,
                             nburnin = n.burnin, 
                             nchains = n.chains)
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|

Have a look to the results. Note that convergence is better.

MCMCsummary(mcmc.multisite, round = 2)
      mean   sd 2.5%  50% 97.5% Rhat n.eff
pA    0.52 0.09 0.35 0.52  0.70 1.01   280
pB    0.40 0.04 0.32 0.40  0.49 1.01   323
phiA  0.61 0.05 0.51 0.61  0.72 1.02   400
phiB  0.69 0.04 0.63 0.69  0.77 1.01   341
psiAB 0.27 0.06 0.16 0.26  0.39 1.00   598
psiBA 0.07 0.02 0.04 0.07  0.12 1.01   488

Multisite model with 3 sites

Now we proceed with the analysis of the 3 sites. Two methods exist to force the movement probabilities from a site to sum to 1 and to be between 0 and 1 at the same time.

Dirichlet prior

Another method is to use Dirichlet priors for the movement parameters. Let’s write the model.

multisite <- nimbleCode({
  
  # -------------------------------------------------
  # Parameters:
  # phiA: survival probability site A
  # phiB: survival probability site B
  # phiC: survival probability site B
  # psiAA = psiA[1]: movement probability from site A to site A (reference)
  # psiAB = psiA[2]: movement probability from site A to site B
  # psiAC = psiA[3]: movement probability from site A to site C 
  # psiBA = psiB[1]: movement probability from site B to site A
  # psiBB = psiB[2]: movement probability from site B to site B (reference)
  # psiBC = psiB[3]: movement probability from site B to site C
  # psiCA = psiC[1]: movement probability from site C to site A
  # psiCB = psiC[2]: movement probability from site C to site B
  # psiCC = psiC[3]: movement probability from site C to site C (reference)
  # pA: recapture probability site A
  # pB: recapture probability site B
  # pC: recapture probability site C
  # -------------------------------------------------
  # States (z):
  # 1 alive at A
  # 2 alive at B
  # 2 alive at C
  # 3 dead
  # Observations (y):  
  # 1 not seen
  # 2 seen at A 
  # 3 seen at B
  # 3 seen at C
  # -------------------------------------------------
  
  # Priors
  phiA ~ dunif(0, 1)
  phiB ~ dunif(0, 1)
  phiC ~ dunif(0, 1)
  pA ~ dunif(0, 1)
  pB ~ dunif(0, 1)
  pC ~ dunif(0, 1)
  # transitions: Dirichlet priors
  psiA[1:3] ~ ddirch(alpha[1:3])
  psiB[1:3] ~ ddirch(alpha[1:3])
  psiC[1:3] ~ ddirch(alpha[1:3])

  # probabilities of state z(t+1) given z(t)
  gamma[1,1] <- phiA * psiA[1]
  gamma[1,2] <- phiA * psiA[2]
  gamma[1,3] <- phiA * psiA[3]
  gamma[1,4] <- 1 - phiA
  gamma[2,1] <- phiB * psiB[1]
  gamma[2,2] <- phiB * psiB[2]
  gamma[2,3] <- phiB * psiB[3]
  gamma[2,4] <- 1 - phiB
  gamma[3,1] <- phiC * psiC[1]
  gamma[3,2] <- phiC * psiC[2]
  gamma[3,3] <- phiC * psiC[3]
  gamma[3,4] <- 1 - phiC
  gamma[4,1] <- 0
  gamma[4,2] <- 0
  gamma[4,3] <- 0
  gamma[4,4] <- 1
  
  # probabilities of y(t) given z(t)
  omega[1,1] <- 1 - pA     # Pr(alive A t -> non-detected t)
  omega[1,2] <- pA         # Pr(alive A t -> detected A t)
  omega[1,3] <- 0          # Pr(alive A t -> detected B t)
  omega[1,4] <- 0          # Pr(alive A t -> detected C t)
  omega[2,1] <- 1 - pB     # Pr(alive B t -> non-detected t)
  omega[2,2] <- 0          # Pr(alive B t -> detected A t)
  omega[2,3] <- pB         # Pr(alive B t -> detected B t)
  omega[2,4] <- 0          # Pr(alive B t -> detected C t)
  omega[3,1] <- 1 - pC     # Pr(alive C t -> non-detected t)
  omega[3,2] <- 0          # Pr(alive C t -> detected A t)
  omega[3,3] <- 0          # Pr(alive C t -> detected B t)
  omega[3,4] <- pC         # Pr(alive C t -> detected C t)
  omega[4,1] <- 1          # Pr(dead t -> non-detected t)
  omega[4,2] <- 0          # Pr(dead t -> detected A t)
  omega[4,3] <- 0          # Pr(dead t -> detected B t)
  omega[4,4] <- 0          # Pr(dead t -> detected C t)
  
  # likelihood 
  for (i in 1:N){
    # latent state at first capture
    z[i,first[i]] <- y[i,first[i]] - 1
    for (t in (first[i]+1):K){
      # z(t) given z(t-1)
      z[i,t] ~ dcat(gamma[z[i,t-1],1:4])
      # y(t) given z(t)
      y[i,t] ~ dcat(omega[z[i,t],1:4])
    }
  }
})

Data in a list.

my.data <- list(y = y + 1)

Constants in a list.

my.constants <- list(first = first, 
                     K = ncol(y), 
                     N = nrow(y),
                     alpha = c(1, 1, 1))

Initial values.

zinits <- y
zinits[zinits==0] <- sample(c(1,2,3), sum(zinits==0), replace = TRUE)
initial.values <- function(){list(phiA = runif(1, 0, 1), 
                                  phiB = runif(1, 0, 1), 
                                  phiC = runif(1, 0, 1), 
                                  psiA = rep(1/3, 3),
                                  psiB = rep(1/3, 3),
                                  psiC = rep(1/3, 3),
                                  pA = runif(1, 0, 1), 
                                  pB = runif(1, 0, 1), 
                                  pC = runif(1, 0, 1),
                                  z = zinits)}

Run nimble.

mcmc.multisite <- nimbleMCMC(code = multisite, 
                             constants = my.constants,
                             data = my.data,              
                             inits = initial.values,
                             monitors = parameters.to.save,
                             niter = n.iter,
                             nburnin = n.burnin, 
                             nchains = n.chains)
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|

Inspect results. Compare to the estimates we obtained with the multinomial logit link method.

MCMCsummary(mcmc.multisite, round = 2)
        mean   sd 2.5%  50% 97.5% Rhat n.eff
pA      0.52 0.08 0.37 0.52  0.69 1.01   327
pB      0.45 0.05 0.37 0.45  0.55 1.00   250
pC      0.25 0.06 0.14 0.24  0.37 1.06   225
phiA    0.60 0.05 0.51 0.60  0.70 1.00   566
phiB    0.70 0.04 0.62 0.70  0.77 1.01   326
phiC    0.76 0.07 0.63 0.77  0.89 1.13   203
psiA[1] 0.74 0.06 0.62 0.74  0.84 1.00   661
psiA[2] 0.24 0.05 0.14 0.23  0.36 1.00   765
psiA[3] 0.02 0.02 0.00 0.02  0.08 1.01   473
psiB[1] 0.07 0.02 0.04 0.07  0.12 1.01   715
psiB[2] 0.84 0.04 0.75 0.84  0.90 1.03   462
psiB[3] 0.09 0.03 0.04 0.09  0.16 1.02   386
psiC[1] 0.02 0.01 0.00 0.02  0.06 1.00  1072
psiC[2] 0.21 0.05 0.12 0.21  0.33 1.04   736
psiC[3] 0.77 0.05 0.64 0.77  0.86 1.04   553

Caterpillar plots.

MCMCplot(mcmc.multisite)

Multistate model

Now we illustrate how one can replace site by breeding state and address another question in evolution ecology: life-history trade-offs. We use data on Soory shearwaters that were kindly provided by David Fletcher. Birds are seen as breeders or non-breeders.

titis <- read_csv2("titis.csv", col_names = FALSE)
y <-  as.matrix(titis)
head(y)
     X1 X2 X3 X4 X5 X6 X7
[1,]  0  0  0  0  0  0  1
[2,]  0  0  0  0  0  0  1
[3,]  0  0  0  0  0  0  1
[4,]  0  0  0  0  0  0  1
[5,]  0  0  0  0  0  0  1
[6,]  0  0  0  0  0  0  1

The code is very similar to that for the model above with two sites, i.e. the structure is the same, but the transition probabilities have different interpretations.

multistate <- nimbleCode({
  
  # -------------------------------------------------
  # Parameters:
  # phiB: survival probability state B
  # phiNB: survival probability state NB
  # psiBNB: transition probability from B to NB
  # psiNBB: transition probability from NB to B
  # pB: recapture probability B
  # pNB: recapture probability NB
  # -------------------------------------------------
  # States (z):
  # 1 alive B
  # 2 alive NB
  # 3 dead
  # Observations (y):  
  # 1 not seen
  # 2 seen as B 
  # 3 seen as NB
  # -------------------------------------------------
  
  # priors
  phiB ~ dunif(0, 1)
  phiNB ~ dunif(0, 1)
  psiBNB ~ dunif(0, 1)
  psiNBB ~ dunif(0, 1)
  pB ~ dunif(0, 1)
  pNB ~ dunif(0, 1)
  
  # probabilities of state z(t+1) given z(t)
  gamma[1,1] <- phiB * (1 - psiBNB)
  gamma[1,2] <- phiB * psiBNB
  gamma[1,3] <- 1 - phiB
  gamma[2,1] <- phiNB * psiNBB
  gamma[2,2] <- phiNB * (1 - psiNBB)
  gamma[2,3] <- 1 - phiNB
  gamma[3,1] <- 0
  gamma[3,2] <- 0
  gamma[3,3] <- 1
  
  # probabilities of y(t) given z(t)
  omega[1,1] <- 1 - pB    # Pr(alive B t -> non-detected t)
  omega[1,2] <- pB        # Pr(alive B t -> detected B t)
  omega[1,3] <- 0         # Pr(alive B t -> detected NB t)
  omega[2,1] <- 1 - pNB   # Pr(alive NB t -> non-detected t)
  omega[2,2] <- 0         # Pr(alive NB t -> detected B t)
  omega[2,3] <- pNB       # Pr(alive NB t -> detected NB t)
  omega[3,1] <- 1         # Pr(dead t -> non-detected t)
  omega[3,2] <- 0         # Pr(dead t -> detected N t)
  omega[3,3] <- 0         # Pr(dead t -> detected NB t)
  
  # likelihood 
  for (i in 1:N){
    # latent state at first capture
    z[i,first[i]] <- y[i,first[i]] - 1
    for (t in (first[i]+1):K){
      # z(t) given z(t-1)
      z[i,t] ~ dcat(gamma[z[i,t-1],1:3])
      # y(t) given z(t)
      y[i,t] ~ dcat(omega[z[i,t],1:3])
    }
  }
})

Data in a list.

my.data <- list(y = y + 1)

Get data of first capture.

get.first <- function(x) min(which(x != 0))
first <- apply(y, 1, get.first)

Constants in a list.

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

Initial values. Same as before, we initialize the latent states with the detections and non-detections. Then non-detections are replaced by a 1 or a 2 for breeder on non-breeder.

zinits <- y
zinits[zinits == 0] <- sample(c(1,2), sum(zinits == 0), replace = TRUE)
initial.values <- function(){list(phiNB = runif(1, 0, 1), 
                                  phiB = runif(1, 0, 1), 
                                  psiNBB = runif(1, 0, 1), 
                                  psiBNB = runif(1, 0, 1), 
                                  pNB = runif(1, 0, 1), 
                                  pB = runif(1, 0, 1), 
                                  z = zinits)}

Parameters to be monitored.

parameters.to.save <- c("phiNB", "phiB","psiNBB", "psiBNB", "pNB", "pB")
parameters.to.save
[1] "phiNB"  "phiB"   "psiNBB" "psiBNB" "pNB"    "pB"    

MCMC details.

n.iter <- 5000
n.burnin <- 2500
n.chains <- 2

Run nimble.

mcmc.multistate <- nimbleMCMC(code = multistate, 
                             constants = my.constants,
                             data = my.data,              
                             inits = initial.values,
                             monitors = parameters.to.save,
                             niter = n.iter,
                             nburnin = n.burnin, 
                             nchains = n.chains)
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|

Let’s inspect the results.

MCMCsummary(mcmc.multistate, round = 2)
       mean   sd 2.5%  50% 97.5% Rhat n.eff
pB     0.60 0.03 0.55 0.60  0.65 1.00   251
pNB    0.56 0.03 0.51 0.56  0.62 1.01   257
phiB   0.80 0.02 0.77 0.80  0.83 1.01   462
phiNB  0.85 0.02 0.82 0.85  0.88 1.02   468
psiBNB 0.25 0.02 0.21 0.25  0.30 1.00   309
psiNBB 0.24 0.02 0.20 0.24  0.28 1.01   447

Caterpillar plot of the parameter estimates.

MCMCplot(mcmc.multistate)

Non-breeder individuals seem to have a survival higher than breeder individuals, suggesting a trade-off between reproduction and survival. Let’s compare graphically the survival of breeder and non-breeder individuals. First we gather the values generated for \(\phi_B\) and \(\phi_{NB}\) for the two chains.

phiB <- c(mcmc.multistate$chain1[,"phiB"], mcmc.multistate$chain2[,"phiB"])
phiNB <- c(mcmc.multistate$chain1[,"phiNB"], mcmc.multistate$chain2[,"phiNB"])
df <- data.frame(param = c(rep("phiB", length(phiB)), rep("phiNB", length(phiB))), 
                 value = c(phiB, phiNB))
head(df)
  param     value
1  phiB 0.8347852
2  phiB 0.7826723
3  phiB 0.7762608
4  phiB 0.7731614
5  phiB 0.7731614
6  phiB 0.7731614

Then, we plot a histogram of the two posterior distributions.

df %>%
  ggplot(aes(x = value, fill = param)) +
  geom_density(color = "white", alpha = 0.6, position = 'identity') +
  scale_fill_manual(values = c("#69b3a2", "#404080")) +
  labs(fill = "", x = "survival")

To formally test this difference, one could fit a model with no state effect on survival and compare the two models with WAIC.