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.36 0.56  0.77 1.04   190
pB    0.39 0.05 0.31 0.39  0.49 1.12   195
phiA  0.56 0.09 0.41 0.56  0.76 1.11   155
phiB  0.72 0.05 0.63 0.72  0.82 1.01   182
piA   0.37 0.11 0.17 0.37  0.59 1.27    73
psiAB 0.24 0.09 0.09 0.23  0.42 1.01   248
psiBA 0.06 0.02 0.02 0.06  0.11 1.02   416

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.53 0.09 0.37 0.53  0.71 1.04   284
pB    0.40 0.04 0.32 0.40  0.48 1.00   309
phiA  0.61 0.05 0.51 0.60  0.71 1.06   425
phiB  0.70 0.03 0.63 0.70  0.77 1.00   374
psiAB 0.27 0.06 0.17 0.26  0.40 1.01   549
psiBA 0.07 0.02 0.04 0.07  0.11 1.03   719

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.51 0.08 0.36 0.50  0.68 1.03   324
pB      0.45 0.05 0.37 0.45  0.55 1.05   309
pC      0.25 0.06 0.14 0.24  0.38 1.01   270
phiA    0.61 0.05 0.50 0.61  0.71 1.01   556
phiB    0.70 0.04 0.63 0.70  0.77 1.04   372
phiC    0.76 0.06 0.64 0.76  0.89 1.02   273
psiA[1] 0.74 0.06 0.62 0.74  0.84 1.03   858
psiA[2] 0.24 0.05 0.15 0.23  0.35 1.03  1004
psiA[3] 0.02 0.02 0.00 0.01  0.08 1.00   419
psiB[1] 0.07 0.02 0.04 0.07  0.12 1.03   772
psiB[2] 0.84 0.04 0.74 0.84  0.90 1.01   334
psiB[3] 0.09 0.03 0.04 0.09  0.17 1.00   317
psiC[1] 0.02 0.02 0.00 0.02  0.06 1.00  1272
psiC[2] 0.21 0.05 0.12 0.21  0.33 1.04   709
psiC[3] 0.77 0.06 0.65 0.77  0.86 1.03   701

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.59 0.03 0.54 0.59  0.65 1.00   264
pNB    0.57 0.03 0.51 0.57  0.62 1.00   241
phiB   0.80 0.02 0.77 0.80  0.83 1.00   537
phiNB  0.85 0.02 0.82 0.85  0.88 1.00   353
psiBNB 0.25 0.02 0.21 0.25  0.29 1.01   363
psiNBB 0.24 0.02 0.20 0.24  0.29 1.01   434

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.8160603
2  phiB 0.8160603
3  phiB 0.8340414
4  phiB 0.8340414
5  phiB 0.8340414
6  phiB 0.8340414

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.