Class 6 live demo: Transition estimation
The team
last updated: 2021-05-14
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 0For 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 0Get 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 4Write 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 <- 2Run 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 410Second, 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 488Multisite 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.
Multinomial logit link
Get the date of first capture.
get.first <- function(x) min(which(x != 0))
first <- apply(y, 1, get.first)Write the model.
multisite <- nimbleCode({
# -------------------------------------------------
# Parameters:
# phiA: survival probability site A
# phiB: survival probability site B
# phiC: survival probability site B
# psiAA: movement probability from site A to site A (reference)
# psiAB = psiA[1]: movement probability from site A to site B
# psiAC = psiA[2]: movement probability from site A to site C
# psiBA = psiB[1]: movement probability from site B to site A
# psiBB: movement probability from site B to site B (reference)
# psiBC = psiB[2]: 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: 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: multinomial logit
# normal priors on logit of all but one transition probs
for (i in 1:2){
lpsiA[i] ~ dnorm(0, sd = 1000)
lpsiB[i] ~ dnorm(0, sd = 1000)
lpsiC[i] ~ dnorm(0, sd = 1000)
}
# constrain the transitions such that their sum is < 1
for (i in 1:2){
psiA[i] <- exp(lpsiA[i]) / (1 + exp(lpsiA[1]) + exp(lpsiA[2]))
psiB[i] <- exp(lpsiB[i]) / (1 + exp(lpsiB[1]) + exp(lpsiB[2]))
psiC[i] <- exp(lpsiC[i]) / (1 + exp(lpsiC[1]) + exp(lpsiC[2]))
}
# last transition probability
psiA[3] <- 1 - psiA[1] - psiA[2]
psiB[3] <- 1 - psiB[1] - psiB[2]
psiC[3] <- 1 - psiC[1] - psiC[2]
# 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])
}
}
})List of data.
my.data <- list(y = y + 1)List of constants.
my.constants <- list(first = first,
K = ncol(y),
N = nrow(y))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),
lpsiA = rnorm(2, 0, 10),
lpsiB = rnorm(2, 0, 10),
lpsiC = rnorm(2, 0, 10),
pA = runif(1, 0, 1),
pB = runif(1, 0, 1),
pC = runif(1, 0, 1),
z = zinits)}Parameters to monitor.
parameters.to.save <- c("phiA", "phiB", "phiC", "psiA", "psiB", "psiC","pA", "pB", "pC")
parameters.to.save[1] "phiA" "phiB" "phiC" "psiA" "psiB" "psiC" "pA" "pB" "pC" 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 the results.
MCMCsummary(mcmc.multisite, round = 2) mean sd 2.5% 50% 97.5% Rhat n.eff
pA 0.52 0.09 0.36 0.52 0.70 1.00 292
pB 0.45 0.05 0.35 0.45 0.54 1.04 300
pC 0.26 0.06 0.16 0.26 0.41 1.01 215
phiA 0.60 0.05 0.49 0.60 0.70 1.01 514
phiB 0.70 0.04 0.63 0.70 0.77 1.03 293
phiC 0.76 0.06 0.64 0.76 0.88 1.01 274
psiA[1] 0.76 0.06 0.64 0.76 0.85 1.00 560
psiA[2] 0.24 0.06 0.14 0.24 0.36 1.00 567
psiA[3] 0.00 0.00 0.00 0.00 0.01 1.32 28
psiB[1] 0.07 0.02 0.03 0.07 0.11 1.00 644
psiB[2] 0.85 0.04 0.77 0.86 0.91 1.03 255
psiB[3] 0.08 0.03 0.03 0.07 0.16 1.07 139
psiC[1] 0.01 0.01 0.00 0.01 0.04 1.01 1132
psiC[2] 0.21 0.06 0.12 0.21 0.34 1.02 468
psiC[3] 0.77 0.06 0.65 0.78 0.87 1.02 445Caterpillar plot.
MCMCplot(mcmc.multisite)
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 553Caterpillar 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 1The 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 <- 2Run 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 447Caterpillar 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.7731614Then, 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.