Class 6 live demo: Transition estimation
last updated: 2023-06-24
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.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 416Second, 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 719Multisite 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.08 0.37 0.51 0.70 1.01 335
pB 0.45 0.05 0.37 0.45 0.55 1.07 349
pC 0.24 0.06 0.14 0.24 0.37 1.18 205
phiA 0.60 0.05 0.50 0.60 0.70 1.01 588
phiB 0.70 0.03 0.63 0.70 0.76 1.04 445
phiC 0.77 0.07 0.64 0.77 0.92 1.25 164
psiA[1] 0.76 0.05 0.64 0.77 0.86 1.00 1028
psiA[2] 0.24 0.05 0.14 0.23 0.36 1.00 1051
psiA[3] 0.00 0.00 0.00 0.00 0.01 1.46 71
psiB[1] 0.07 0.02 0.04 0.07 0.11 1.00 588
psiB[2] 0.85 0.04 0.76 0.85 0.91 1.08 249
psiB[3] 0.09 0.03 0.04 0.08 0.16 1.12 148
psiC[1] 0.01 0.01 0.00 0.01 0.04 1.00 1162
psiC[2] 0.20 0.05 0.11 0.20 0.32 1.15 538
psiC[3] 0.79 0.05 0.67 0.79 0.88 1.16 504Caterpillar 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.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 701Caterpillar 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.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 434Caterpillar 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.8340414Then, 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.