Workflow for Ridgway’s Hawk Integrated Population Model
Brian W. Rolek
2025-02-15
- 1. Integrated Population Model
(IPM)
- 1.1 Model Code
- 1.2 Results- Figures and Tables
- 1.2.1 Population Trends and Structure
- 1.2.2 Apparent Survival, Recruitment, and Detection
- 1.2.3 Productivity
- 1.2.4 Correlations Among Parameters
- 1.2.5 Population Growth Rates Over Time
- 1.2.6 Correlations Between Population Growth Rates and Demographics
- 1.2.7 Transient Life Table Response Experiments (tLTREs)
- 2. Population Viability Analysis (PVA)
- 3. Additional IPM plots
Supplemental materials for: Brian W. Rolek, Leah Dunn, Marta Curti, Russell Thorstrom, Evan Buechley, Carolina Granthon, Jeff Johnson, Carlos Suárez, Gabriela Diaz, and Christopher J. W. McClure. 2025. Management influences population persistence of the Critically Endangered Ridgway’s Hawk (Buteo ridgwayi) and creates tradeoffs between sites
Contact information: rolek.brian@peregrinefund.org
A permanent archive and DOI is available at: https://zenodo.org/doi/XXXXX
Metadata, data, and scripts used in analyses can be found at https://github.com/The-Peregrine-Fund/Ridgways-Hawk-IPM.
The full workflow below is visible as a html website at: https://the-peregrine-fund.github.io/Ridgways-Hawk-IPM/.
1. Integrated Population Model (IPM)
1.1 Model Code
Code for the Integrated Population Model.
library('nimble')
library('parallel')
library ('coda')
#load("/bsuscratch/brianrolek/riha_ipm/data.rdata")
load("data/data.rdata")
cpus <- 5
#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
# Observations (po) = y
# 1 seen first-year (age=0, just before 1st b-day)
# 2 seen nonbreeder
# 3 seen breeder
# 4 not seen
# States (ps)
# 1 alive first-year
# 2 alive nonbreeder
# 3 alive breeder
# 4 dead
# Groups
# 1 wild-born
# 2 translocated and hacked
###################################################
# PARAMETERS
# phiFY: survival probability first year
# phiA: survival probability nonbreeders
# phiB: survival probability breeders
# psiFYB: recruitment probability from first-year to breeder
# psiAB: recruitment probability from nonbreeders to breeder
# psiBA: recruitment probability from breeder to nonbreeders
# pFY: resight probability first-years
# pA: resight probability nonbreeders
# pB: resight probability breeders
# lmu.prod: mean productivity (males and females) per territory on the log scale
# sds: standard deviations for multivariate normal random effects for time and site
# R: correlation coefficients for multivariate normal random effects for time and site
# lambda: population growth rate (derived)
# extinct: binary indicator of extirpation at a site
# gamma: coefficient of effect from nest treatments
# betas: coefficient of effect from translocations
# deltas: coefficient of effect from survey effort
# mus: overall means for survival, recruitment, and detections
# r and rr: "r" parameter for negative binomial distribution
# also described as omega in manuscript
#**********************
#* Model code
#**********************
mycode <- nimbleCode(
{
###################################################
# Priors and constraints
###################################################
# survival, recruitment, and detection can be correlated
for (k in 1:8){
betas[k] ~ dnorm(0, sd=20) # prior for translocations coefficients
deltas[k] ~ dnorm(0, sd=10) # prior for survey effort coefficients
} # k
for (j in 1:8){
for (s in 1:nsite){
lmus[j,s] <- logit(mus[j,s])
mus[j,s] ~ dbeta(1,1) # prior for overall means
}} #
# Temporal random effects and correlations between sites
# Non-centered parameterization of the multivariate normal distribution to improve convergence
for (jj in 1:p){ sds[jj] ~ dexp(1) }# prior for temporal variation
R[1:p,1:p] <- t(Ustar[1:p,1:p]) %*% Ustar[1:p,1:p] # calculate rhos, correlation coefficients
Ustar[1:p,1:p] ~ dlkj_corr_cholesky(eta=1.1, p=p) # Ustar is the Cholesky decomposition of the correlation matrix
for (t in 1:nyr){ # survival params only have nyr-1, no problem to simulate from however
for (s in 1:nsite){
eta[1:p,s,t] <- diag(sds[1:p]) %*% t(Ustar[1:p,1:p]) %*% z.score[1:p,s,t]
for(j in 1:p){
z.score[j,s,t] ~ dnorm(0, sd=1) # z-scores
} # j
} } # s t
#######################
# Derived params
#######################
for (s in 1:nsite){
for (t in 1:(nyr-1)){
lambda[t, s] <- Ntot[t+1, s]/(Ntot[t, s])
loglambda[t, s] <- log(lambda[t, s])
}} #t
###############################
# Likelihood for productivity
###############################
# Priors
for (s in 1:nsite){
lmu.prod[s] ~ dnorm(0, sd=5)
} # s
gamma ~ dnorm(0, sd=10)
rr ~ dexp(0.05)
# Productivity likelihood
for (k in 1:npairsobs){
prod[k] ~ dnegbin(ppp[k], rr)
ppp[k] <- rr/(rr+mu.prod[k])
log(mu.prod[k]) <- lmu.prod[site.pair[k]] +
gamma*treat.pair[k] +
#eps[9, year.pair[k] ] +
eta[9, site.pair[k], year.pair[k] ]
} # k
# Derive yearly productivity for population model
# need to reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.pair is a matrix of indices for each site
for (t in 1:nyr){
for (s in 1:nsite){
for (xxx in 1:pair.end[t,s]){
prodmat[t,s,xxx] <- mu.prod[ yrind.pair[xxx,t,s] ]
} # xxx
mn.prod[t,s] <- mean( prodmat[t,s,1:pair.end[t,s]] )
}} # s t
# GOF for productivity
for (k in 1:npairsobs){
f.obs[k] <- prod[k] # observed counts
f.exp[k] <- mu.prod[k] # expected counts adult breeder
f.rep[k] ~ dnegbin(ppp[k], rr) # expected counts
f.dssm.obs[k] <- abs( ( f.obs[k] - f.exp[k] ) / (f.obs[k]+0.001) )
f.dssm.rep[k] <- abs( ( f.rep[k] - f.exp[k] ) / (f.rep[k]+0.001) )
} # k
f.dmape.obs <- sum(f.dssm.obs[1:npairsobs])
f.dmape.rep <- sum(f.dssm.rep[1:npairsobs])
f.tvm.obs <- sd(brood[1:npairsobs])^2/mean(brood[1:npairsobs])
f.tvm.rep <- sd(f.rep[1:npairsobs])^2/mean(f.rep[1:npairsobs])
################################
# Likelihood for counts
################################
# Abundance for year=1
for (v in 1:7){
for (s in 1:nsite){
# subtract one to allow dcat to include zero
N[v, 1, s] <- N2[v, 1, s] - 1
N2[v, 1, s] ~ dcat(pPrior[v, 1:s.end[v,s], s]) # Priors differ for FYs and adults
}} # s t
# Abundance for years > 1
for (t in 1:(nyr-1)){
for (s in 1:nsite){
# Number of wild born juvs
N[1, t+1, s] ~ dpois( (NFY[t, s]*mn.phiFY[t, s]*mn.psiFYB[t, s] + # first year breeders
NF[t, s]*mn.phiA[t, s]*mn.psiAB[t, s] + # nonbreeders to breeders
NB[t, s]*mn.phiB[t, s]*(1-mn.psiBA[t, s])) # breeders remaining
*mn.prod[t+1, s]/2 ) # end Poisson
# Abundance of nonbreeders
N[2, t+1, s] ~ dbin(mn.phiFY[t, s]*(1-mn.psiFYB[t, s]), NFY[t, s]) # Nestlings to nonbreeders
N[3, t+1, s] ~ dbin(mn.phiA[t, s]*(1-mn.psiAB[t, s]), NF[t, s]) # Nonbreeders to nonbreeders
N[4, t+1, s] ~ dbin(mn.phiB[t, s]*mn.psiBA[t, s], NB[t, s]) # Breeders to nonbreeders
# Abundance of breeders
N[5, t+1, s] ~ dbin(mn.phiFY[t, s]*mn.psiFYB[t, s], NFY[t, s]) # Nestlings to second year breeders
N[6, t+1, s] ~ dbin(mn.phiA[t, s]*mn.psiAB[t, s], NF[t, s]) # Nonbreeder to breeder
N[7, t+1, s] ~ dbin(mn.phiB[t, s]*(1-mn.psiBA[t, s]), NB[t, s]) # Breeder to breeder
}} # s t
# Count likelihoods, state-space model, and observation process
for (t in 1:nyr){
for (s in 1:nsite){
countsAdults[t, s] ~ dpois(lamAD[t, s]) # adult females
constraint_data[t, s] ~ dconstraint( (N[1, t, s] + hacked.counts[t, s]) >= 0 ) # constrain N1+hacked.counts to be >=0
NFY[t, s] <- N[1, t, s] + hacked.counts[t, s] # Transfers translocated first-year females
NF[t, s] <- sum(N[2:4, t, s]) # number of adult nonbreeder females
NB[t, s] <- sum(N[5:7, t, s]) # number of adult breeder females
NAD[t, s] <- NF[t, s] + NB[t, s] # number of adults
Ntot[t, s] <- sum(N[1:7, t, s]) # total number of females
log(lamAD[t, s]) <- log(NAD[t, s]) + deltas[1]*effort2[t, s] + deltas[2]*effort2[t, s]^2
log(lamFY[t, s]) <- log(N[1, t, s]) + deltas[3]*effort2[t, s] + deltas[4]*effort2[t, s]^2
}# s
# First-years at different sites have different distributions
# for better model fit
countsFY[t, 1] ~ dpois(lamFY[t, 1]) # doesn't have any zeroes so poisson fits
countsFY[t, 2] ~ dnegbin(pp[t], r) # first year females, includes translocated/hacked
pp[t] <- r/(r+(lamFY[t, 2] ))
} # t
r ~ dexp(0.05)
###################
# Assess GOF of the state-space models for counts
# Step 1: Compute statistic for observed data
# Step 2: Use discrepancy measure: mean absolute error
# Step 3: Use test statistic: number of turns
###################
for (t in 1:nyr){
c.repFY[t, 1] ~ dpois( lamFY[t, 1] )
c.repFY[t, 2] ~ dnegbin( pp[t], r )
for (s in 1:nsite){
c.expAD[t, s] <- lamAD[t, s] # expected counts adult breeder
c.expFY[t, s] <- lamFY[t, s]
c.obsAD[t, s] <- countsAdults[t, s]
c.obsFY[t, s] <- countsFY[t, s] # first year
c.repAD[t, s] ~ dpois( lamAD[t, s] ) # simulated counts
dssm.obsAD[t, s] <- abs( ( (c.obsAD[t, s]) - (c.expAD[t, s]) ) / (c.obsAD[t, s]+0.001) )
dssm.obsFY[t, s] <- abs( ( (c.obsFY[t, s]) - (c.expFY[t, s]) ) / (c.obsFY[t, s]+0.001) )
dssm.repAD[t, s] <- abs( ( (c.repAD[t, s]) - (c.expAD[t, s]) ) / (c.repAD[t, s]+0.001) )
dssm.repFY[t, s] <- abs( ( (c.repFY[t, s]) - (c.expFY[t, s]) ) / (c.repFY[t, s]+0.001) )
}} # t
dmape.obs[1] <- sum(dssm.obsAD[1:nyr, 1:nsite])
dmape.obs[2] <- sum(dssm.obsFY[1:nyr, 1:nsite])
dmape.rep[1] <- sum(dssm.repAD[1:nyr, 1:nsite])
dmape.rep[2] <- sum(dssm.repFY[1:nyr, 1:nsite])
tvm.obs[1] <- sd(c.obsB[1:nyr, 1:nsite])^2/mean(c.obsB[1:nyr, 1:nsite])
tvm.obs[2] <- sd(c.obsFY[1:nyr, 1:nsite])^2/mean(c.obsFY[1:nyr, 1:nsite])
tvm.rep[1] <- sd(c.repB[1:nyr, 1:nsite])^2/mean(c.repB[1:nyr, 1:nsite])
tvm.rep[2] <- sd(c.repFY[1:nyr, 1:nsite])^2/mean(c.repFY[1:nyr, 1:nsite])
# # Test statistic for number of turns
for (s in 1:nsite){
for (t in 1:(nyr-2)){
tt1.obsAD[t, s] <- step(c.obsAD[t+2, s] - c.obsAD[t+1, s])
tt2.obsAD[t, s] <- step(c.obsAD[t+1, s] - c.obsAD[t, s])
tt3.obsAD[t, s] <- equals(tt1.obsAD[t, s] + tt2.obsAD[t, s], 1)
tt1.obsFY[t, s] <- step(c.obsFY[t+2, s] - c.obsFY[t+1, s])
tt2.obsFY[t, s] <- step(c.obsFY[t+1, s] - c.obsFY[t, s])
tt3.obsFY[t, s] <- equals(tt1.obsFY[t, s] + tt2.obsFY[t, s], 1)
}} # t
tturn.obs[1] <- sum(tt3.obsAD[1:(nyr-2), 1:nsite])
tturn.obs[2] <- sum(tt3.obsFY[1:(nyr-2), 1:nsite])
for (s in 1:nsite){
for (t in 1:(nyr-2)){
tt1.repAD[t, s] <- step(c.repAD[t+2, s] - c.repAD[t+1, s])
tt2.repAD[t, s] <- step(c.repAD[t+1, s] - c.repAD[t, s])
tt3.repAD[t, s] <- equals(tt1.repAD[t, s] + tt2.repAD[t, s], 1)
tt1.repFY[t, s] <- step(c.repFY[t+2, s] - c.repFY[t+1, s])
tt2.repFY[t, s] <- step(c.repFY[t+1, s] - c.repFY[t, s])
tt3.repFY[t, s] <- equals(tt1.repFY[t, s] + tt2.repFY[t, s], 1)
}} # t
tturn.rep[1] <- sum(tt3.repAD[1:(nyr-2), 1:nsite])
tturn.rep[2] <- sum(tt3.repFY[1:(nyr-2), 1:nsite])
################################
# Likelihood for survival
################################
# Calculate averages for sites each year for integration
for (t in 1:nyr){
for (s in 1:nsite){
for (xxxx in 1:surv.end[t,s]){
# Reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.surv is a matrix of indices for each site
phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
} # xxxx
mn.phiFY[ t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] )
mn.phiA[ t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
mn.phiB[ t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
mn.psiFYB[ t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
mn.psiAB[ t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
mn.psiBA[ t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
mn.pA[ t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
mn.pB[ t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
}}
for (i in 1:nind){
for (t in 1:nyr){
#Survival
logit(phiFY[i,t]) <- eta[1, site[i,t],t] +
lmus[1, site[i,t]] + betas[1]*hacked[i] # first year
logit(phiA[i,t]) <- eta[2, site[i,t],t] +
lmus[2, site[i,t]] #+ betas[2]*hacked[i] # nonbreeder
logit(phiB[i,t]) <- eta[3, site[i,t],t] +
lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
#Recruitment
logit(psiFYB[i,t]) <- eta[4, site[i,t],t] +
lmus[4, site[i,t]] + betas[4]*hacked[i] # first year to breeder
logit(psiAB[i,t]) <- eta[5, site[i,t],t] +
lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
logit(psiBA[i,t]) <-
lmus[6, site[i,t]] # breeder to nonbreeder
#Re-encounter
logit(pA[i,t]) <- eta[7, site[i,t],t] +
lmus[7, site[i,t]] + betas[7]*hacked[i] +
deltas[5]*effort2[t, site[i,t]] + deltas[6]*effort2[t, site[i,t]]^2# resight of nonbreeders
logit(pB[i,t]) <- eta[8, site[i,t],t] +
lmus[8, site[i,t]] + betas[8]*hacked[i] +
deltas[7]*effort2[t, site[i,t]] + deltas[8]*effort2[t, site[i,t]]^2# resight of breeders
}#t
}#i
# Define state-transition and observation matrices
for (i in 1:nind){
# Define probabilities of state S(t+1) given S(t)
for (t in first[i]:(nyr-1)){
ps[1,i,t,1] <- 0
ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t])
ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
ps[1,i,t,4] <- (1-phiFY[i,t])
ps[2,i,t,1] <- 0
ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
ps[2,i,t,4] <- (1-phiA[i,t])
ps[3,i,t,1] <- 0
ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
ps[3,i,t,4] <- (1-phiB[i,t])
ps[4,i,t,1] <- 0
ps[4,i,t,2] <- 0
ps[4,i,t,3] <- 0
ps[4,i,t,4] <- 1
# Define probabilities of O(t) given S(t)
po[1,i,t,1] <- 1
po[1,i,t,2] <- 0
po[1,i,t,3] <- 0
po[1,i,t,4] <- 0
po[2,i,t,1] <- 0
po[2,i,t,2] <- pA[i,t]
po[2,i,t,3] <- 0
po[2,i,t,4] <- (1-pA[i,t])
po[3,i,t,1] <- 0
po[3,i,t,2] <- 0
po[3,i,t,3] <- pB[i,t]
po[3,i,t,4] <- (1-pB[i,t])
po[4,i,t,1] <- 0
po[4,i,t,2] <- 0
po[4,i,t,3] <- 0
po[4,i,t,4] <- 1
} #t
} #i
# Likelihood
for (i in 1:nind){
# Define latent state at first capture
z[i,first[i]] <- y[i,first[i]]
for (t in (first[i]+1):nyr){
# State process: draw S(t) given S(t-1)
z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
# Observation process: draw O(t) given S(t)
y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
} #t
} #i
} )
#**********************
#* Function to run model in NIMBLE
#**********************
run_ipm <- function(info, datl, constl, code){
library('nimble')
library('coda')
params <- c(
# pop growth
"lambda",
# fecundity
"lmu.prod", "gamma", "rr", "mn.prod",
# survival
"mus", "lmus", "betas", "deltas",
# abundance
"NB", "NF", "NFY", "N", "NAD",
"r",
"N", "Ntot",
# error terms
"eta", "sds", "Ustar", "R", "z.score",
# yearly summaries
'mn.phiFY', 'mn.phiA', 'mn.phiB',
'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
'mn.pA', 'mn.pB',
# goodness of fit
"f.dmape.obs", "f.dmape.rep",
"f.tvm.obs", "f.tvm.rep",
"dmape.obs", "dmape.rep",
"tvm.obs", "tvm.rep",
"tturn.obs", "tturn.rep"
)
n.chains=1; n.thin=100; n.iter=100000; n.burnin=50000
mod <- nimbleModel(code,
constants = constl,
data = datl,
inits = info$inits,
buildDerivs = FALSE, # doesn't work when TRUE, no hope for HMC
calculate=T
)
cmod <- compileNimble(mod, showCompilerOutput = TRUE)
confhmc <- configureMCMC(mod)
confhmc$setMonitors(params)
hmc <- buildMCMC(confhmc)
chmc <- compileNimble(hmc, project = mod,
resetFunctions = TRUE,
showCompilerOutput = TRUE)
post <- runMCMC(chmc,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
thin = n.thin,
samplesAsCodaMCMC = T,
setSeed = info$seed)
return(post)
} # run_ipm function end
#*****************
#* Run chains in parallel
#*****************
this_cluster <- makeCluster(cpus)
post <- parLapply(cl = this_cluster,
X = par_info[1:cpus],
fun = run_ipm,
datl = datl,
constl = constl,
code = mycode)
stopCluster(this_cluster)
# save(post, mycode,
# file="/bsuscratch/brianrolek/riha_ipm/outputs/ipm.rdata")1.2 Results- Figures and Tables
Simple summary statistics
library ('MCMCvis')
library ('coda')
library ('ggplot2')
library('reshape2')
library ('tidybayes')
library ('bayestestR')
library ('ggpubr')
library('viridis')
library ('HDInterval')
library ('abind')
load("data/data.rdata")
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_longrun.rdata")
out <- lapply(post[1:4], as.mcmc) # omit chain 5, causing lack of convergence in 1 param
outp <- MCMCpstr(out, type="chains")
#********************
#* Calculate summary stats
#********************
#* Change in abundance at both sites
# 2011 to 2023 pop growth,
# growth between start and end years of monitoring
Nall <- apply(outp$Ntot, c(1,3) , sum)
dimnames(Nall) <- list( Year=2011:2023, Iter=1:2000)
l.all.post2 <- Nall[13,] / Nall[1,]
median((l.all.post2-1)*100)## [1] 135.6286HDInterval::hdi((l.all.post2-1)*100)## lower upper
## 100.0000 173.2026
## attr(,"credMass")
## [1] 0.95# Annual average
l.all.post <- array(NA, dim=c(constl$nyr-1, ncol(Nall)) )
for (t in 2:13){
l.all.post[t-1, ] <- Nall[t,] / Nall[t-1,]
}
# Annual average pop growth
# geometric mean is required here because
# population growth rate multiplies across years
geom.mean <- function(x){exp(mean(log(x)))}
l.mn <- (apply(l.all.post, 2, geom.mean)-1)*100
# Average annual pop growth
median(l.mn)## [1] 7.40364HDInterval::hdi(l.mn)## lower upper
## 5.946309 8.736096
## attr(,"credMass")
## [1] 0.95# Population growth rates, site-specific
# percent change from start and end years on monitoring
l.all.post3 <- outp$Ntot[13,,] / outp$Ntot[1,,]
((apply(l.all.post3, 1, median)-1)*100)## [1] 102.9586 1155.0000((apply(l.all.post3, 1, HDInterval::hdi)-1)*100)## [,1] [,2]
## lower 69.61326 436.3636
## upper 138.56209 2366.6667mn.lambda <- apply(outp$lambda, c(2,3), geom.mean)
(md.lambda <- ((apply(mn.lambda, 1, median)-1)*100))## [1] 6.076037 23.467706(hdi95.lambda <- ((apply(mn.lambda, 1, HDInterval::hdi)-1)*100))## [,1] [,2]
## lower 4.577884 15.77612
## upper 7.573550 30.91046# Overall averages
# from survival model
mn.mus <- apply(outp$mus, c(1,3), mean)
(md.mus <- apply(mn.mus, 1, median)) |> round(2)## [1] 0.35 0.94 0.91 0.11 0.44 0.10 0.18 0.79(hdi.mus <- apply(mn.mus, 1, HDInterval::hdi)) |> round(2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## lower 0.20 0.83 0.84 0.02 0.30 0.05 0.02 0.56
## upper 0.55 1.00 0.97 0.27 0.55 0.14 0.38 0.95# Does detection differ by life stage?
pd(mn.mus[8,]-mn.mus[7,])## Probability of Direction: 1.00# Compare between sites
df.sites <- data.frame(mus=1:8, pd=rep(NA,8), greater=NA)
for (i in 1:8){
first <- median(outp$mus[i, 1, ]) > median(outp$mus[i, 2, ])
df.sites$greater[i] <- ifelse(first, "LH", "PC")
if (first){
diff <- outp$mus[i, 1, ]-outp$mus[i, 2, ]
df.sites$pd[i] <- mean(diff>0) |> round(2)
} else{
diff <- outp$mus[i, 2, ]-outp$mus[i, 1, ]
df.sites$pd[i] <- mean(diff>0) |> round(2)
}
}
df.sites(md.mus <- apply(outp$mus, c(1,2), median)) |> round(2)## [,1] [,2]
## [1,] 0.42 0.26
## [2,] 0.99 0.88
## [3,] 0.95 0.88
## [4,] 0.13 0.07
## [5,] 0.15 0.73
## [6,] 0.15 0.04
## [7,] 0.04 0.29
## [8,] 0.72 0.88(hdi.mus <- apply(outp$mus, c(1,2), HDInterval::hdi)) |> round(2)## , , 1
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## lower 0.20 0.97 0.89 0.01 0.07 0.08 0.0 0.37
## upper 0.69 1.00 0.99 0.34 0.22 0.23 0.2 0.96
##
## , , 2
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## lower 0.06 0.67 0.74 0.00 0.49 0.00 0.04 0.57
## upper 0.54 1.00 0.98 0.35 0.96 0.08 0.69 1.00# productivity
# see section below
# Compare survival, recruitment, and detection rates
# among life stages within each site
df.comp <- data.frame(param1 = c('phiA', 'phiA', 'phiB', 'psiA', 'pB'),
param2 = c('phiB', 'phiFY', 'phiFY', 'psiFY', 'pA'),
pd_LH = NA, pd_PC = NA)
## Los Haitises, survival
df.comp$pd_LH[1] <- (outp$lmus[2,1,]-outp$lmus[3,1,]) |> pd() |> round(2)
df.comp$pd_LH[2] <- (outp$lmus[2,1,]-outp$lmus[1,1,]) |> pd() |> round(2)
df.comp$pd_LH[3] <- (outp$lmus[3,1,]-outp$lmus[1,1,]) |> pd() |> round(2)
## Los Haitises, recruitment rates
df.comp$pd_LH[4] <- (outp$lmus[5,1,]-outp$lmus[4,1,]) |> pd() |> round(2)
## Los Haitises, resight probability
df.comp$pd_LH[5] <- (outp$lmus[8,1,]-outp$lmus[7,1,]) |> pd() |> round(2)
## Punta Canas, survival
df.comp$pd_PC[1] <- (outp$lmus[2,2,]-outp$lmus[3,2,]) |> pd() |> round(2)
df.comp$pd_PC[2] <- (outp$lmus[2,2,]-outp$lmus[1,2,]) |> pd() |> round(2)
df.comp$pd_PC[3] <- (outp$lmus[3,2,]-outp$lmus[1,2,]) |> pd() |> round(2)
## Punta Cana, recruitment rates
df.comp$pd_PC[4] <- (outp$lmus[5,2,]-outp$lmus[4,2,]) |> pd() |> round(2)
## Punta Cana, resight probability
df.comp$pd_PC[5] <- (outp$lmus[8,2,]-outp$lmus[7,2,]) |> pd() |> round(2)
df.comp #*******************
#* Plot IPM results
#*******************1.2.1 Population Trends and Structure
#Plot population size
lNall <- melt(Nall, value.name="Abundance")
lNall$Site <- "Both sites"
lNall <- lNall[,c(1,4,2,3)]
Ntot <- outp$Ntot[1:13,,]
dimnames(Ntot) <- list( Year=2011:2023, Site=c("Los Haitises", "Punta Cana"),
Iter=1:2000)
lN <- melt(Ntot, value.name="Abundance")
lN <- rbind(lN, lNall)
med_hdis <- function(x, cm){
df <- data.frame(y=median(x),
ymin=HDInterval::hdi(x, credMass=cm)[1],
ymax=HDInterval::hdi(x, credMass=cm)[2] )
}
counts.tot <- datl$countsAdults+datl$countsFY+constl$hacked.counts[,-3]
df.counts <- melt(counts.tot)
colnames(df.counts) <- c("Year", "Site_ab", "Count")
df.counts$Site <- ifelse(df.counts$Site_ab=="LHNP", "Los Haitises", "Punta Cana")
df.counts.all <- data.frame(Year=2011:2023, Site_ab="Both sites", Count=counts.tot[,1] + counts.tot[,2], Site="Both sites" )
df.counts <- rbind(df.counts, df.counts.all)
p1 <- ggplot() + theme_minimal(base_size=14) +
geom_line(data=lN, aes(x=Year, y=Abundance, group=Iter),
color="gray30", linewidth=0.1, alpha=0.01 ) +
stat_summary(data=lN, aes(x=Year, y=Abundance),
geom="errorbar", width=0,
fun.data =med_hdis, fun.args=list(cm=0.95) ) +
stat_summary(data=lN, aes(x=Year, y=Abundance),
geom="errorbar", width=0, linewidth=1,
fun.data =med_hdis, fun.args=list(cm=0.80) ) +
stat_summary(data=lN, aes(x=Year, y=Abundance),
geom="line", linewidth=1,
fun.data =function(x){data.frame(y=median(x))} ) +
geom_line(data=df.counts, aes(x=Year, y=Count),
col="black", linewidth=1, linetype=2 ) +
facet_wrap("Site", scales="free_y") +
xlab("Year") + ylab("Number")
p1
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Abundance and counts.tiff",
# plot=p1,
# device="tiff",
# width=9, height=4, dpi=300)# Plot life stage structure
mdFY <- apply(outp$NFY[1:13,,], c(1,2), median)
mdB <- apply(outp$NB[1:13,,], c(1,2), median)
mdF <- apply(outp$NF[1:13,,], c(1,2), median)
lFY <- melt(mdFY)
lB <- melt(mdB)
lF <- melt(mdF)
lFY$Stage <- "First-year"
lB$Stage <- "Breeder"
lF$Stage <- "Nonbreeder"
ldat <- rbind(lFY, lB, lF)
colnames(ldat)[1:3] <- c("Year", "Sitenum", "Number")
ldat$Site <- ifelse(ldat$Sitenum==1, "Los Haitises", "Punta Cana")
ldat$year <- ldat$Year+2010
# Use median number of females in each stage
# to plot an approximate population structure
p3 <- ggplot(ldat, aes(fill=Stage, y=as.numeric(Number), x=year)) +
geom_bar(position="fill", stat="identity") +
ylab("Proportion of population") + xlab("Year") +
facet_wrap("Site") +
theme_minimal(base_size=14)+ theme(legend.position="top") +
scale_fill_viridis_d(option = "mako")
p4 <- ggplot(ldat, aes(fill=Stage, y=as.numeric(Number), x=year)) +
geom_bar(position="stack", stat="identity") +
ylab("Median number of females") + xlab("Year") +
facet_wrap("Site", scales = "free_y") +
theme_minimal(base_size=14 ) + theme(legend.position="top") +
scale_fill_viridis_d(option = "mako")
# plot with uncertainty
lFY <- melt(outp$NFY)
lB <- melt(outp$NB)
lF <- melt(outp$NF)
colnames(lFY) <- colnames(lB) <- colnames(lF) <- c("Time", "Site", "Iteration", "Abundance")
lFY$Year <- lFY$Time + 2010
lB$Year <- lB$Time + 2010
lF$Year <- lF$Time + 2010
lFY$Stage <- "First-year"
lB$Stage <- "Breeder"
lF$Stage <- "Nonbreeder"
lALL <- rbind(lFY, lB, lF)
lALL$Site <- ifelse(lALL$Site==1, "Los Haitises", "Punta Cana")
mdFY <- apply(outp$NFY[1:13,,], c(1,2), median)
mdB <- apply(outp$NB[1:13,,], c(1,2), median)
mdF <- apply(outp$NF[1:13,,], c(1,2), median)
hdiFY <- apply(outp$NFY[1:13,,], c(1,2), HDInterval::hdi)
hdiB <- apply(outp$NB[1:13,,], c(1,2), HDInterval::hdi)
hdiF <- apply(outp$NF[1:13,,], c(1,2), HDInterval::hdi)
medFY <- melt(mdFY)
medB <- melt(mdB)
medF <- melt(mdF)
hdisFY <- melt(hdiFY)
hdisB <- melt(hdiB)
hdisF <- melt(hdiF)
dfstages <- rbind(medFY, medB, medF)
dfhdis <- rbind(hdisFY, hdisB, hdisF)
dfstages$LHDI <- dfhdis[dfhdis$Var1=="lower",]
dfstages$UHDI <- dfhdis[dfhdis$Var1=="upper",]
cols <- mako(3)
med_hdis <- function(x, cm){
df <- data.frame(y=median(x),
ymin=HDInterval::hdi(x, credMass=cm)[1],
ymax=HDInterval::hdi(x, credMass=cm)[2] )
}
p5 <- ggplot() + theme_minimal(base_size=14) +
theme(legend.position="none") +
geom_line(data=lALL, aes(x=Year, y=Abundance,
group=Iteration, color=Stage,
alpha=Stage),
linewidth=0.1) +
stat_summary(data=lALL, aes(x=Year, y=Abundance),
geom="errorbar", width=0,
fun.data =med_hdis, fun.args=list(cm=0.95) ) +
stat_summary(data=lALL, aes(x=Year, y=Abundance),
geom="errorbar", width=0, linewidth=1,
fun.data =med_hdis, fun.args=list(cm=0.80) ) +
stat_summary(data=lALL, aes(x=Year, y=Abundance),
geom="line", linewidth=1,
fun.data =function(x){data.frame(y=median(x))} ) +
scale_color_manual( values = c("First-year"=cols[2], "Breeder"=cols[1], "Nonbreeder"=cols[3]) ) +
scale_alpha_manual( values = c("First-year"=0.015, "Breeder"=0.0075, "Nonbreeder"=0.15) ) +
facet_wrap(Stage~Site, scales="free_y",
shrink=TRUE, ncol=2) +
xlab("Year") + ylab("Number of females") +
theme(
strip.background = element_blank(),
strip.text.x = element_blank()
)
p35 <- ggarrange(p3, p5, ncol=1, nrow=2)
p35
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Stage structure_abundance.tiff",
# plot=p35,
# device="tiff",
# width=6, height=8, dpi=300)1.2.2 Apparent Survival, Recruitment, and Detection
# plot survival, recruitment, and detection
mus.mat <- outp$mus
dimnames(mus.mat)[[1]] <- c("FY",
"A\nSurvival", "B",
"FY:B",
"A:B\nRecruitment",
"B:A",
"A", "B\nDetection")
dimnames(mus.mat)[[2]] <- c("Los Haitises", "Punta Cana")
lmus <- melt(mus.mat)
colnames(lmus) <- c("Parameter", "Site", "Iter", "value")
p6 <- ggplot(data= lmus, aes(x=value, y=Parameter )) +
stat_pointinterval(point_interval=median_hdci, .width = c(.95, .80)) +
coord_flip() +
facet_wrap("Site") +
xlab("Probability") + ylab("Parameter") +
theme_bw(base_size=14) +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())
p6
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Survival Recruitment Detection.tiff",
# plot=p6,
# device="tiff",
# width=8, height=4, dpi=300)
# print estimates
mds <- tapply(lmus$value, list(lmus$Parameter, lmus$Site), median) |> round(2)
hdi_LH <- tapply(lmus$value, list(lmus$Parameter, lmus$Site), HDInterval::hdi)[,1]
hdi_LH2 <- do.call(rbind, hdi_LH) |> round(2)
hdi_PC <- tapply(lmus$value, list(lmus$Parameter, lmus$Site), HDInterval::hdi)[,2]
hdi_PC2 <- do.call(rbind, hdi_PC) |> round(2)
df <- data.frame(param= c("phiFY", "phiA", "phiB", "psiFY:B", "psiA:B", "psiB:A", "pA", "pB"),
md_LH=mds[,1], lhdi_LH=hdi_LH2[,1], uhdi_LH=hdi_LH2[,2],
md_PC=mds[,2], lhdi_PC=hdi_PC2[,1], uhdi_PC=hdi_PC2[,2])
df# Plot predicted survival, recruitment, and detection
# with effects from translocations and hacking
# Birds were only sourced from LHNP and
# translocated to PC here
# so we only predict values for PC
betas.pd <- apply(outp$betas, 1, pd)
ind <- which (betas.pd>0.975) # which are significant
pred.mus <- array(NA, dim=c(dim(outp$mus)))
for (m in ind){
for (tr in 1:2){
pred.mus[m,tr,] <- plogis( outp$lmus[m,2,] + outp$betas[m,]*c(0,1)[tr] )
}}
# treatment, mu, site, iter
mus.md <- apply(pred.mus, c(1,2), median)
mus.HDI80 <- apply(pred.mus, c(1,2), HDInterval::hdi, credMass=0.8)
mus.HDI95 <- apply(pred.mus, c(1,2), HDInterval::hdi)
dimnames(pred.mus) <- list(Params=c("First-year\nSurvival",
"Nonbreeder\nSurvival",
"Breeder\nSurvival",
"First-year\nRecruitment",
"Nonbreeder\nRecruitment",
"Breeder to Nonbreeder\nTransition",
"Nonbreeder\nDetection",
"Breeder\nDetection"),
Translocated=c("Not Translocated",
"Translocated"),
Iter=1:dim(pred.mus)[[3]])
lpms <- melt(pred.mus)
colnames(lpms)[[4]] <- "Probability"
lpms <- lpms[!is.na(lpms$Probability),]
p.betas <- ggplot(lpms, aes(x=Translocated, y=Probability) ) +
stat_pointinterval(point_interval=median_hdci, .width = c(.95, .80)) +
facet_wrap("Params", nrow=1) +
theme_bw(base_size=14) +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()) +
xlab("") +
scale_y_continuous(breaks = seq(0, 1, by = 0.2),
minor_breaks = seq(0, 1, by = 0.1)) +
scale_x_discrete(labels=c("Not Translocated"="Not\nTranslocated",
"Translocated"="Trans-\nlocated"))
p.betas
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Translocation effects_ggplot.tiff",
# plot=p.betas,
# device="tiff",
# width=8, height=4, dpi=300)1.2.3 Productivity
# Plot fecundity in response to nest treatments
f <- exp(outp$lmu.prod)
f.md <- apply(f, 1, median)
f.HDI80 <- apply(f, 1, HDInterval::hdi, credMass=0.8)
f.HDI95 <- apply(f, 1, HDInterval::hdi)
rownames(f) <- c("Los Haitises", "Punta Cana")
lf <- melt(f)
colnames(lf) <- c("Site", "Iter", "Productivity")
lf$Treatment <- "Untreated"
# Calculate treatment effects on fecundity
f.pred <- array(NA, dim=dim(outp$lmu.prod))
for (s in 1:2){
f.pred[s,] <- exp(outp$lmu.prod[s,] + outp$gamma[,1])
} # s
f2.md <- apply(f.pred, 1, median)
f2.HDI80 <- apply(f.pred, 1, HDInterval::hdi, credMass=0.8)
f2.HDI95 <- apply(f.pred, 1, HDInterval::hdi)
rownames(f.pred) <- c("Los Haitises", "Punta Cana")
lf.pred <- melt(f.pred)
colnames(lf.pred) <- c("Site", "Iter", "Productivity")
lf.pred$Treatment <- "Treated"
lprod <- rbind(lf,lf.pred)
lprod$Treatment <- factor(lprod$Treatment, levels=c('Untreated', 'Treated'))
p.prod <- ggplot(lprod, aes(x=Treatment, y=Productivity) ) +
stat_pointinterval(point_interval=median_hdci, .width = c(.95, .80), expand=T) +
facet_wrap("Site") +
theme_bw(base_size=14) +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()) +
xlab("")+
scale_y_continuous(breaks = seq(0, 1.5, by = 0.4),
minor_breaks = seq(0, 1.5, by = 0.1))
p.prod
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Fecundity_treatment_ggplot.tiff",
# plot=p.prod,
# device="tiff",
# width=4, height=3, dpi=300)
# Calculate probability of direction
f.diff <- f[1,]-f[2,]
pd(f.diff) |> round(2)## Probability of Direction: 1.00f.md |> round(2)## lmu.prod[1] lmu.prod[2]
## 0.37 0.23f.HDI95 |> round(2)## lmu.prod[1] lmu.prod[2]
## lower 0.30 0.16
## upper 0.45 0.29f.diff2 <- f.pred[1,]-f.pred[2,]
pd(f.diff2) |> round(2)## Probability of Direction: 1.00f2.md |> round(2)## [1] 1.16 0.71f2.HDI95 |> round(2)## [,1] [,2]
## lower 0.93 0.50
## upper 1.41 0.91median(f.pred[1,]/f[1,])## [1] 3.116722median(f.pred[2,]/f[2,])## [1] 3.1167221.2.4 Correlations Among Parameters
R is the correlation coefficent between rates. Indices for R are
demographics and detection parameters in both dimensions. For example,
R[1,2] is the correlation coefficient between first-year survival and
nonbreeder survival.
1. First-year Survival
2. Nonbreeder Survival
3. Breeder Survival
4. First-year recruitment to breeder
5. Nonbreeder recruitment to breeder
6. Breeder transition to nonbreeder
7. Nonbreeder detection
8. Breeder detection
9. Productivity
#*******************
# Correlations among demographic rates and detection
#*******************
pdR <- apply(outp$R, c(1,2), pd) |> round(2)
max(pdR[pdR<1])## [1] 0.82apply(outp$R, c(1,2), median) |> round(2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 1.00 0.01 0.23 0.10 0.05 0.01 0.08 0.15 -0.08
## [2,] 0.01 1.00 0.04 -0.06 -0.01 0.00 0.00 0.07 -0.10
## [3,] 0.23 0.04 1.00 -0.04 -0.01 0.00 0.25 0.03 -0.13
## [4,] 0.10 -0.06 -0.04 1.00 -0.02 -0.01 0.20 -0.18 -0.15
## [5,] 0.05 -0.01 -0.01 -0.02 1.00 -0.03 -0.01 0.09 -0.12
## [6,] 0.01 0.00 0.00 -0.01 -0.03 1.00 -0.02 0.05 0.04
## [7,] 0.08 0.00 0.25 0.20 -0.01 -0.02 1.00 -0.22 -0.20
## [8,] 0.15 0.07 0.03 -0.18 0.09 0.05 -0.22 1.00 0.09
## [9,] -0.08 -0.10 -0.13 -0.15 -0.12 0.04 -0.20 0.09 1.00apply(outp$R, c(1,2), HDInterval::hdi)[1,,] |> round(2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 1.00 -0.61 -0.31 -0.52 -0.56 -0.61 -0.52 -0.35 -0.55
## [2,] -0.61 1.00 -0.53 -0.60 -0.59 -0.59 -0.58 -0.47 -0.69
## [3,] -0.31 -0.53 1.00 -0.53 -0.59 -0.60 -0.30 -0.40 -0.59
## [4,] -0.52 -0.60 -0.53 1.00 -0.56 -0.55 -0.30 -0.68 -0.63
## [5,] -0.56 -0.59 -0.59 -0.56 1.00 -0.67 -0.52 -0.46 -0.69
## [6,] -0.61 -0.59 -0.60 -0.55 -0.67 1.00 -0.54 -0.61 -0.53
## [7,] -0.52 -0.58 -0.30 -0.30 -0.52 -0.54 1.00 -0.60 -0.65
## [8,] -0.35 -0.47 -0.40 -0.68 -0.46 -0.61 -0.60 1.00 -0.42
## [9,] -0.55 -0.69 -0.59 -0.63 -0.69 -0.53 -0.65 -0.42 1.00apply(outp$R, c(1,2), HDInterval::hdi)[2,,] |> round(2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 1.00 0.58 0.70 0.63 0.58 0.57 0.62 0.63 0.40
## [2,] 0.58 1.00 0.60 0.59 0.59 0.60 0.61 0.69 0.51
## [3,] 0.70 0.60 1.00 0.53 0.55 0.58 0.72 0.48 0.33
## [4,] 0.63 0.59 0.53 1.00 0.53 0.60 0.71 0.34 0.44
## [5,] 0.58 0.59 0.55 0.53 1.00 0.55 0.60 0.58 0.40
## [6,] 0.57 0.60 0.58 0.60 0.55 1.00 0.62 0.62 0.64
## [7,] 0.62 0.61 0.72 0.71 0.60 0.62 1.00 0.23 0.32
## [8,] 0.63 0.69 0.48 0.34 0.58 0.62 0.23 1.00 0.54
## [9,] 0.40 0.51 0.33 0.44 0.40 0.64 0.32 0.54 1.00plt <- function(object, params,...) {
MCMCplot(object=out,
params=params,
guide_axis=TRUE,
HPD=TRUE, ci=c(80, 95), horiz=FALSE,
#ylim=c(-10,10),
...)
}
ind <- 1
Rs <- c()
for (i in 1:(nrow(outp$R)-1)){
for (j in (i+1):nrow(outp$R)){
Rs[ind] <- paste0("R[",i,", ", j, "]")
ind <- ind+1
}}
par(mfrow=c(2,1))
plt(object=out, params=Rs[1:18], exact=TRUE, ISB=FALSE,
main="Correlations btw demographic rates\n over time and sites",
xlab = "Rhos", guide_lines=TRUE)
plt(object=out, params=Rs[19:36], exact=TRUE, ISB=FALSE,
main="Correlations btw demographic rates\n over time and sites, continued ...",
xlab = "Rhos", guide_lines=TRUE)
1.2.5 Population Growth Rates Over Time
#*******************
# Population growth rates over time
#******************
lam.m <- apply(outp$lambda[1:12,,], c(1,2), mean)
lam.hdi <- apply(outp$lambda[1:12,,], c(1,2), HDInterval::hdi)
par(mfrow=c(1,2))
plot(2012:2023, lam.m[,1], type="b", pch=1,
ylab="Population growth rate", xlab="Year",
ylim=c(min(lam.hdi[,,1]), max(lam.hdi[,,1])),
main="Los Haitises")
abline(h=1, lty=2)
segments(x0=2012:2023, x1=2012:2023,
y0 = lam.hdi[1,,1], y1= lam.hdi[2,,1])
plot(2012:2023, lam.m[,2], type="b", pch=1, lty=1,
ylab="Population growth rate", xlab="Year",
ylim=c(min(lam.hdi[,,2]), max(lam.hdi[,,2])),
main="Punta Cana")
abline(h=1, lty=2)
segments(x0=2012:2023, x1=2012:2023,
y0 = lam.hdi[1,,2], y1= lam.hdi[2,,2])
1.2.6 Correlations Between Population Growth Rates and Demographics
#*******************
# Correlations between demographics and population growth rates
#*******************
# create a function to plot correlations
# between demographics and population growth rates
plot.cor <- function (lam.post, x.post, x.lab, ind.x=1:12, ind.y=1:12, site=1, yaxt="b"){
# calculate the correlation coefficient
# over each iteration to propagate error
lambda <- lam.post[ind.y, site,]
x <- x.post[ind.x, site, ]
df <- data.frame(ats1=c(0.8, 0.9, 1.0, 1.1, 1.2, 1.3),
ats2=c(1.0, 1.5, 2.0, 2.5, 3.0, 3.5),
labs1=c("0.8", NA, "1.0", NA, "1.2", NA),
labs2=c("1.0", NA, "2.0", NA, "3.0", NA)
)
if(site==1){ df <- df[, c(1,3)]} else{
df <- df[, c(2,4)]
}
lwd <- 1
cor.post <- array(NA, dim=dim(lambda)[-1])
for (i in 1:dim(lambda)[[2]]){
cor.post[i] <- cor(lambda[,i], x[,i])
}
cor.df <- data.frame(median= median(cor.post) |> round(digits=2),
ldi=HDInterval::hdi(cor.post)[1] |> round(digits=2),
hdi=HDInterval::hdi(cor.post)[2] |> round(digits=2),
pd = pd(cor.post) |> round(digits=2))
lam.m <- apply(lambda, 1, geom.mean)
lam.hdi <- apply(lambda, 1, HDInterval::hdi)
x.m <- apply(x, 1, median)
x.hdi <- apply(x, 1, HDInterval::hdi)
x.lims <- c(min(x.hdi[,]), max(x.hdi[,]))
y.lims <- c(min(lam.hdi[,]), max(lam.hdi[,]))
if(yaxt=="n"){
plot(x.m, lam.m[ind.y],
xlim= x.lims,
ylim= y.lims,
type="n",
ylab="", xlab=x.lab, cex.lab=1.5,
main=c("", "")[site],
yaxt=yaxt, xaxt="n")
axis(2, at=df[,1], labels=c(rep(NA, 5)))
} else {
plot(x.m, lam.m[ind.y],
xlim= x.lims,
ylim= y.lims,
type="n",
ylab="", xlab=x.lab, cex.lab=1.5,
main=c("", "")[site],
yaxt="n", xaxt="n")
axis(2, at=df[,1], labels=df[,2], las=1, cex.axis=1.5)
}
axis(1, cex.axis=1.5)
points(x.m, lam.m, pch=1)
segments(x0=x.hdi[1,], x1=x.hdi[2,],
y0 = lam.m, y1= lam.m, lwd=lwd)
segments(x0=x.m, x1=x.m,
y0 = lam.hdi[1,], y1= lam.hdi[2,], lwd=lwd)
xs <- c(x.lims[1], x.lims[1]+(x.lims[2]-x.lims[1]*1.5))
ys <- c( (y.lims[2]-y.lims[1])+y.lims[1], (y.lims[2]-y.lims[1])+y.lims[1],
(y.lims[2]-y.lims[1])*0.6+y.lims[1], (y.lims[2]-y.lims[1])*0.6+y.lims[1])
polygon(x=c(xs, rev(xs)), y=ys, col=alpha("white", 0.7), border=NA )
text(x = x.lims[1], y = (y.lims[2]-y.lims[1])*0.9+y.lims[1],
paste("r = ", cor.df$median, " (",cor.df$ldi,", ", cor.df$hdi, ")", sep=""),
pos = 4, font = 3, cex = 1.5)
text(x = x.lims[1], y = (y.lims[2]-y.lims[1])*0.7+y.lims[1], paste("pd = ", cor.df$pd, sep=""),
pos = 4, font = 3, cex = 1.5)
}
# Plot correlations between population growth rates
# and demographics.
# "r" is a correlation coefficient and represents
# the magnitude of the correlation
# P(r>0) is the probability of direction (similar to p-values)
# that is, the probability that an effect exists
par(mfrow=c(4,4), mar=c(4,1,3,1), oma=c(1,5,1,1))
plot.cor(outp$lambda, outp$mn.prod, x.lab="Productivity", ind.x=2:13)
plot.cor(outp$lambda, outp$mn.phiFY, x.lab="First-year Survival", yaxt="n")
mtext("Los Haitises", side=3, line=1, cex=2)
plot.cor(outp$lambda, outp$mn.phiA, x.lab="Nonbreeder Survival", yaxt="n")
plot.cor(outp$lambda, outp$mn.phiB, x.lab="Breeder Survival", yaxt="n")
plot.cor(outp$lambda, outp$mn.psiFYB, x.lab="First-year to Breeder")
plot.cor(outp$lambda, outp$mn.psiAB, x.lab="Nonbreeder to Breeder", yaxt="n")
plot.cor(outp$lambda, array( abs(constl$hacked.counts[,1]), c(13, 2, 2000)), x.lab="Number translocated", yaxt="n")
plot.new()
# tiff(height=4, width=8, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//PopGrowthand Demographics2.tiff")
# par(mfrow=c(2,3), mar=c(4,1,1,1), oma=c(0,5,3,0))
plot.cor(outp$lambda, outp$mn.prod, x.lab="Productivity", ind.x=2:13, site=2)
plot.cor(outp$lambda, outp$mn.phiFY, x.lab="First-year Survival", yaxt="n", site=2)
mtext("Punta Cana", side=3, line=1, cex=2)
plot.cor(outp$lambda, outp$mn.phiA, x.lab="Nonbreeder Survival", yaxt="n", site=2)
plot.cor(outp$lambda, outp$mn.phiB, x.lab="Breeder Survival", site=2, yaxt="n")
plot.cor(outp$lambda, outp$mn.psiFYB, x.lab="First-year to Breeder", site=2)
plot.cor(outp$lambda, outp$mn.psiAB, x.lab="Nonbreeder to Breeder", yaxt="n", site=2)
plot.cor(outp$lambda, array( abs(constl$hacked.counts[,2]), c(13, 2, 2000)), x.lab="Number translocated", yaxt="n", site=2)
mtext("Population growth rate", side=2, outer=TRUE, line=2.5, cex=2)
# dev.off()
1.2.7 Transient Life Table Response Experiments (tLTREs)
Population Growth Rates and Contributions from Demographics
#********************
# Transient Life Table Response Experiment (tLTRE)
# Overall contributions
#********************
lambda <- expression( ( # numerator start
(nfy*phiFY*psiFYB*f*0.5 + # repro of first year breeders
nf*phiA*psiAB*f*0.5 + # repro of nonbreeders to breeders
nb*phiB*(1-psiBA)*f*0.5) + # repro of breeders remaining breeders
nfy*phiFY*(1-psiFYB) +
nf*phiA*(1-psiAB) +
nb*phiB*psiBA +
nfy*phiFY*psiFYB +
nb*phiB*(1-psiBA) ) / # numerator end
(nfy+nf+nb) # denominator
) # expression end
# Calculate proportional population sizes
nyr <- nrow(outp$NFY)-1
niter <- dim(outp$NFY)[[3]]
n1 <- outp$NFY[1:nyr, ,] / outp$Ntot[1:nyr, ,]
n2 <- outp$NF[1:nyr, ,] / outp$Ntot[1:nyr, ,]
n3 <- outp$NB[1:nyr, ,] / outp$Ntot[1:nyr, ,]
# Extract the mean demographic rates and population sizes and store them in a list
mu <- list(phiFY=outp$mus[1,,], phiA=outp$mus[2,,], phiB=outp$mus[3,,],
psiFYB=outp$mus[4,,], psiAB=outp$mus[5,,], psiBA=outp$mus[6,,],
f=exp(outp$lmu.prod),
nfy=apply(n1, c(2,3), mean), nf=apply(n2, c(2,3), mean), nb=apply(n3, c(2,3), mean))
# Calculate sensitivities
nms <- c('phiFY', 'phiA', 'phiB',
'psiFYB', 'psiAB', 'psiBA',
'f',
'nfy', 'nf', 'nb' )
sens.func <- function( name, expr=lambda, envir=mu) { eval(D(expr=expr, name=name), envir=envir) }
sensl <- lapply(nms, sens.func)
names(sensl) <- nms
sens <- abind(sensl, along=3)
# Define matrix to store results
cont <- array(NA, dim=c(2, niter, 10) )
dimnames(cont)[[3]] <- nms
# Calculate contributions for each demographic rate and stage-structured population size at each MCMC draw
for (i in 1:niter){
for (s in 1:2){
dp_stoch <- cbind(outp$mn.phiFY[1:12,s,i], outp$mn.phiF[1:12,s,i], outp$mn.phiB[1:12,s,i],
outp$mn.psiFYB[1:12,s,i], outp$mn.psiAB[1:12,s,i], outp$mn.psiBA[1:12,s,i],
outp$mn.prod[2:13,s,i],
n1[1:12,s,i], n2[1:12,s,i], n3[1:12,s,i])
# Derive process variance and covariance among demographic parameters using shrunk estimates of
# demographic rates and proportional pop. sizes
dp_varcov <- var(dp_stoch)
sensvec <- sens[s, i, ]
# Calculate demographic contributions
contmatrix <- dp_varcov * outer(sensvec, sensvec)
cont[s, i, ] <- rowSums(contmatrix)
}}
CRI <- apply(cont[,,], c(1,3), quantile, probs=c(0.025, 0.975))
names.arg <- c(expression(italic(phi)[italic(FY)]),
expression(italic(phi)[italic(A)]),
expression(italic(phi)[italic(B)]),
expression(italic(psi)[italic(FY:B)]),
expression(italic(psi)[italic(A:B)]),
expression(italic(psi)[italic(B:A)]),
expression(italic(f)),
expression(italic(n)[italic(FY)]),
expression(italic(n)[italic(A)]),
expression(italic(n)[italic(B)]))
# tiff(height=4, width=6, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//LTRE_overall.tiff")
par(mfrow=c(2,1), mar=c(3,5,1,1))
a <- barplot(colMeans(cont[1,,]), names.arg=names.arg,
ylab=expression( paste("Contribution to ", lambda) ), las=1,
ylim=range(CRI[,1,]), col=rep("grey65", 10),
border=rep("grey65", 10),
main="Los Haitises", font.main=1,
yaxt="n")
axis(2, at=c(0, 0.01, 0.02, 0.03), las=1)
segments(a, CRI[1,1,], a, CRI[2,1,], lwd=1.5)
b <- barplot(colMeans(cont[2,,]), names.arg=names.arg,
ylab=expression( paste("Contribution to ", lambda) ), las=1,
ylim=range(CRI[,2,]), col=rep("grey65", 10),
border=rep("grey65", 10),
main="Punta Cana", font.main=1,
yaxt="n")
axis(2, at=c(0, 0.01, 0.02, 0.03), las=1)
segments(a, CRI[1,2,], a, CRI[2,2,], lwd=1.5)
# dev.off()
#*******************
# Transient Life Table Response Experiment (tLTRE)
# Inter-annual Contributions
#*******************
sens <- diff <- array(NA, dim=c(constl$nyr-1, 2, 10, dim(outp$mn.phiA)[[3]]),
dimnames=list(Year=2012:2023,
site=c('Los Haitises', 'Punta Cana'),
Param=c('phiFY', 'phiA', 'phiB',
'psiFYB', 'psiAB', 'psiBA',
'f',
'nfy', 'nf', 'nb' ),
Iter=1:dim(outp$mn.phiA)[[3]]) )
getDiff <- function(x){
x[2:constl$nyr, , ] - x[1:(constl$nyr-1), , ]
}
diff[,,'phiFY',] <- getDiff(outp$mn.phiFY)
diff[,,'phiA',] <- getDiff(outp$mn.phiA)
diff[,,'phiB',] <- getDiff(outp$mn.phiB)
diff[,,'psiFYB',] <- getDiff(outp$mn.psiFYB)
diff[,,'psiAB',] <- getDiff(outp$mn.psiAB)
diff[,,'psiBA',] <- getDiff(outp$mn.psiBA)
diff[,,'f',] <- getDiff(outp$mn.prod)
diff[,,'nfy',] <- getDiff( outp$NFY/outp$Ntot )
diff[,,'nf',] <- getDiff( outp$NF/outp$Ntot )
diff[,,'nb',] <- getDiff( outp$NB/outp$Ntot )
getMn <- function(x) {
( x[2:constl$nyr, , ] + x[1:(constl$nyr-1), , ] ) / 2
}
means <- list('phiFY' = getMn(outp$mn.phiFY),
'phiA' = getMn(outp$mn.phiA),
'phiB' = getMn(outp$mn.phiB),
'psiFYB' = getMn(outp$mn.psiFYB),
'psiAB' = getMn(outp$mn.psiAB),
'psiBA' = getMn(outp$mn.psiBA),
'f' = getMn(outp$mn.prod),
'nfy' = getMn(outp$NFY),
'nf' = getMn(outp$NF),
'nb' = getMn(outp$NB) )
sens[,,'phiFY',] <- eval( D(lambda, "phiFY" ), envir = means )
sens[,,'phiA',] <- eval( D(lambda, "phiA" ), envir = means )
sens[,,'phiB',] <- eval( D(lambda, "phiB" ), envir = means )
sens[,,'psiFYB',] <- eval( D(lambda, "phiFYB" ), envir = means )
sens[,,'psiAB',] <- eval( D(lambda, "psiAB" ), envir = means )
sens[,,'psiBA',] <- eval( D(lambda, "psiBA" ), envir = means )
sens[,,'f',] <- eval( D(lambda, "f" ), envir = means )
sens[,,'nfy',] <- eval( D(lambda, "nfy" ), envir = means )
sens[,,'nf',] <- eval( D(lambda, "nf" ), envir = means )
sens[,,'nb',] <- eval( D(lambda, "nb" ), envir = means )
conts <- diff*sens
# ~~~~ figure 9.7 ~~~~
# Create a plot for the sequential differences in finite growth rates
niter <- dim(outp$mn.phiA)[[3]]
nyrs <- constl$nyr-1
diff.lam <- outp$lambda[2:nyrs, , ] - outp$lambda[1:(nyrs-1), ,]
differences <- cbind(2012:2022, apply(diff.lam, c(1,2), mean))
differences <- rbind(differences, c(2023, NA, NA))
# Mean contributions
V <- apply(conts, 3:1, mean)
V1 <- V2 <- V
V1[V1<0] <- 0
V2[V2>0] <- 0
# Make figure
colors <- c(colorRampPalette(c("blue4", "lightblue"))(3),
colorRampPalette(c("darkred", "yellow"))(3),
colorRampPalette(c("darkviolet"))(1),
colorRampPalette(c("lightgreen", "darkgreen"))(3))
legend.text <- c(expression(italic(phi)[italic(FY)]),
expression(italic(phi)[italic(A)]),
expression(italic(phi)[italic(B)]),
expression(italic(psi)[italic(FY:B)]),
expression(italic(psi)[italic(A:B)]),
expression(italic(psi)[italic(B:A)]),
expression(italic(f)),
expression(italic(n)[italic(FY)]),
expression(italic(n)[italic(A)]),
expression(italic(n)[italic(B)]))
# tiff(height=6.5, width=8, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//LTRE.tiff")
par(mfrow=c(2,2), mar=c(1, 2, 0, 1), oma=c(7,7,1,1))
txt.size <- 1.5
yl <- list(bquote( paste("Difference in population")),
bquote( paste("growth rate (", Delta, lambda, ")" )) )
barplot(differences[,2], ylim=c(-0.16, 0.16),
#ylab=expression( paste("Difference in population", "\ngrowth rate (", Delta, lambda, ")" )),
ylab= "",
xpd=NA,
axes=FALSE, border=NA,
cex.lab=txt.size, main="Los Haitises", cex.main=txt.size,
font.main=1 )
mtext(do.call(expression,yl), side=2, cex=txt.size,
line=c(5,3.5) )
axis(2, las=1)
abline(h=0)
barplot(differences[,3], ylim=c(-1.5, 1.5),
ylab="",
xpd=NA,
axes=FALSE, border=NA, main="Punta Cana",
cex.main=txt.size,
font.main=1)
axis(2, las=1)
abline(h=0)
barplot(V1[,1,], ylim=c(-0.19, 0.19),
col=colors, ylab="",
xlab="",
xpd=NA,
axes=FALSE,
border=NA,
cex.lab=txt.size)
mtext(expression( paste("Contribution to ", Delta, lambda) ), side=2, cex=txt.size,
line=4.25)
barplot(V2[,1,], add=TRUE, col=colors, axes=FALSE, border=NA)
axis(2, las=1)
abline(h=0)
barplot(V1[,2,], ylim=c(-0.12, 0.12),
col=colors, ylab="",
xlab="", axes=FALSE,
names.arg=2012:2023,
border=NA)
barplot(V2[,2,], add=TRUE, col=colors, axes=FALSE, border=NA)
axis(2, las=1)
abline(h=0)
mtext(side=1 , outer=TRUE, "Beginning year", line=2.5, cex=txt.size)
reset <- function() {
par(mfrow=c(1, 1), oma=rep(0, 4), mar=rep(0, 4), new=TRUE)
plot(0:1, 0:1, type="n", xlab="", ylab="", axes=FALSE)
}
reset()
legend(legend = legend.text, fill=colors,
bty="n", x="bottom", border=NA,
xpd=NA, x.intersp=0.005, text.width=NA,
cex=1.5, horiz=TRUE)
# dev.off()
#*******************2. Population Viability Analysis (PVA)
2.1 Model Code
Code for Population Viability Analysis.
library('nimble')
library('parallel')
library ('coda')
# load("/bsuscratch/brianrolek/riha_ipm/data.rdata")
# source("/bsuscratch/brianrolek/riha_ipm/MCMCvis.R")
# load("/bsuscratch/brianrolek/riha_ipm/outputs/ipm_shortrun.rdata")
# load IPM results to get good inits
load("C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//outputs//ipm_longrun.rdata")
load("data//data.rdata")
library ("MCMCvis")
#**********************
#* Set data to reflect PVA scenarios
#**********************
cpus <- 5 # number of processors
constl$K <- 50 # number of future years
constl$SC <- 45 # number of PVA scenarios
datl$constraint_data <- rbind(datl$constraint_data, array(1, dim=c(constl$K,2)) ) # to help constrain FYs born + hacked to be positive
constl$effort2 <- rbind(constl$effort2, array(0, dim=c(constl$K,2), dimnames=list(2024:(2024+constl$K-1), c("LH", "PC"))))
constl$num.treated <- rep( c(0, 15, 30, 45, 100, 0, 15, 30, 45, 100, 0, 15, 30, 45, 100), each=3 ) # 100s sub in for All but are over-ridden in model code
num.hacked <- rep( c(0, 0, 0, 0, 0, 5, 5, 5, 5, 5, 10, 10, 10, 10, 10), each=3 )
surv.diff <- array(NA, dim=c(constl$SC, 3, 2), dimnames=list(scenario=1:constl$SC, stage=c('FY', 'NB', 'B'), site=c('LH', 'PC')))
surv.diff[,1,] <- matrix(c(rep( c(0, 0.130, 0.260), 15),
rep( c(0, 0.146, 0.292), 15)), ncol=2)
surv.diff[,2,] <- matrix(c(rep( c(0, 0.006, 0.012), 15),
rep( c(0, 0.022, 0.044), 15)), ncol=2)
surv.diff[,3,] <- matrix(c(rep( c(0, 0.016, 0.032), 15),
rep( c(0, 0.026, 0.052), 15)), ncol=2)
constl$surv.diff <- surv.diff
hacked.counts <- array(0, dim=c(constl$SC, constl$nyr+constl$K, 2))
for (sc in 1:constl$SC){
hacked.counts[sc,1:constl$nyr, 1:2] <- constl$hacked.counts[, 1:2]
hacked.counts[sc,14:(constl$nyr+constl$K), 1] <- -num.hacked[sc]
hacked.counts[sc,14:(constl$nyr+constl$K), 2] <- num.hacked[sc]
}
constl$hacked.counts <- hacked.counts
#**********************
#* Specify initial values
#* from posterior of IPM
#**********************
#* This helps ensure that all chains run
# Identify chains with NAs that failed to initialize
out <- lapply(post, as.mcmc)
NAlist <- c()
for (i in 1:length(out)){
NAlist[i] <- any (is.na(out[[i]][,1:286]) | out[[i]][,1:286]<0)
}
out <- out[!NAlist]
outp <- MCMCpstr(out, type="chains")
!NAlist
Ni.func <- function (){
# take the means of the posterior for inits
Ni <- apply(outp$N, c(1,2,3), mean) |> round()
Ni.pva <- array(NA, dim=c(constl$SC, dim(Ni)+c(0,constl$K,0)))
for (sc in 1:constl$SC){
Ni.pva[sc, 1:7, 1:constl$nyr, 1:2] <- Ni
for (t in constl$nyr:(constl$nyr+constl$K)){
for (s in 1:2){
Ni.pva[sc, 1:7, t, s] <- Ni[1:7, 13, s]
}}} # t sc
return(Ni.pva)
} # function
u2 <- apply(outp$Ustar, c(1,2), mean)
inits.func.pva <- function (){
list(
# fecundity
lmu.prod = apply(outp$lmu.prod, 1, mean),
gamma = mean(outp$gamma),
rr = mean(outp$rr),
# survival
z = par_info[[1]]$inits$z,
mus = apply(outp$mus, c(1,2), mean),
betas = apply(outp$betas, 1, mean),
deltas = apply(outp$deltas, 1, mean),
sds = apply(outp$sds, 1, mean),
Ustar = u2,
r = mean(outp$r),
N = Ni.func()
)}
# Set seed then save for reproducibility
# then randomly draw for each chain
set.seed(1)
seeds <- sample(1:1000000, size=cpus, replace=FALSE)
par_info_pva <- list()
for (i in 1:cpus){
par_info_pva[[i]] <- list(seed=seeds[i], inits = inits.func.pva())
}
#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
# Observations (po) = y
# 1 seen first-year (age=0, just before 1st b-day)
# 2 seen nonbreeder
# 3 seen breeder
# 4 not seen
# States (ps)
# 1 alive first-year
# 2 alive nonbreeder
# 3 alive breeder
# 4 dead
# Groups
# 1 wild-born
# 2 translocated and hacked
###################################################
# PARAMETERS
# phiFY: survival probability first year
# phiA: survival probability nonbreeders
# phiB: survival probability breeders
# psiFYB: recruitment probability from first-year to breeder
# psiAB: recruitment probability from nonbreeders to breeder
# psiBA: recruitment probability from breeder to nonbreeders
# pFY: resight probability first-years
# pA: resight probability nonbreeders
# pB: resight probability breeders
# lmu.prod: mean productivity (males and females) per territory on the log scale
# sds: standard deviations for multivariate normal random effects for time and site
# R: correlation coefficients for multivariate normal random effects for time and site
# lambda: population growth rate (derived)
# extinct: binary indicator of extirpation at a site
# gamma: coefficient of effect from nest treatments
# betas: coefficient of effect from translocations
# deltas: coefficient of effect from survey effort
# mus: overall means for survival, recruitment, and detections
# r and rr: "r" parameter for negative binomial distribution
# also described as omega in manuscript
#**********************
#* Model code
#**********************
mycode <- nimbleCode(
{
###################################################
# Priors and constraints
###################################################
# survival, recruitment, and detection can be correlated
for (k in 1:8){
betas[k] ~ dnorm(0, sd=20) # prior for coefficients
deltas[k] ~ dnorm(0, sd=10)
} # k
for (j in 1:8){
for (s in 1:nsite){
lmus[j,s] <- logit(mus[j,s])
mus[j,s] ~ dbeta(1,1) # prior for means
}} # m population #s sex #h hacked
# Temporal random effects and correlations between sites
# Non-centered parameterization of the multivariate normal distribution to improve convergence
for (jj in 1:p){ sds[jj] ~ dexp(1) }# prior for temporal variation
R[1:p,1:p] <- t(Ustar[1:p,1:p]) %*% Ustar[1:p,1:p] # calculate rhos, correlation coefficients
Ustar[1:p,1:p] ~ dlkj_corr_cholesky(eta=1.1, p=p) # Ustar is the Cholesky decomposition of the correlation matrix
for (t in 1:(nyr+K)){ # survival params only have nyr-1, no problem to simulate from however
for (s in 1:nsite){
eta[1:p,s,t] <- diag(sds[1:p]) %*% t(Ustar[1:p,1:p]) %*% z.score[1:p,s,t]
for(j in 1:p){
z.score[j,s,t] ~ dnorm(0, sd=1) # z-scores
} # j
} } # s t
#######################
# Derived params
#######################
for(sc in 1:SC){
for (s in 1:nsite){
for (t in 1:(nyr-1+K)){
lambda[sc, t, s] <- Ntot[sc, t+1, s]/(Ntot[sc, t, s]) # population groewth rate
loglambda[sc, t, s] <- log(lambda[sc, t, s]) # r
}} #t
for (s in 1:nsite){
for (t in 1:K){
# quasi-extinction probabilities given population segments, N total, adults, and breeders
extinct[sc, t, s] <- equals(Ntot[sc, nyr+t, s], 0)
extinctAD[sc, t, s] <- equals(NAD[sc, nyr+t, s], 0)
extinctB[sc, t, s] <- equals(NB[sc, nyr+t, s], 0)
}}
} # sc
###############################
# Likelihood for productivity
###############################
# Priors
for (s in 1:nsite){
lmu.prod[s] ~ dnorm(0, sd=5)
} # s
gamma ~ dnorm(0, sd=10)
rr ~ dexp(0.05)
# Productivity likelihood
for (k in 1:npairsobs){
prod[k] ~ dnegbin(ppp[k], rr)
ppp[k] <- rr/(rr+mu.prod[k])
log(mu.prod[k]) <- lmu.prod[site.pair[k]] +
gamma*treat.pair[k] +
eta[9, site.pair[k], year.pair[k] ]
} # k
# Derive yearly productivity for population model
# need to reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.pair is a matrix of indices for each site
for (t in 1:nyr){
for (s in 1:nsite){
for (xxx in 1:pair.end[t,s]){
prodmat[t,s,xxx] <- mu.prod[ yrind.pair[xxx,t,s] ]
} # xxx
mn.prod[1,t,s] <- mean( prodmat[t,s,1:pair.end[t,s]] )
for(sc in 2:SC){
mn.prod[sc, t, s] <- mn.prod[1,t,s]
}}} # s t
# Future projections- fecundity
for(sc in 1:SC){
for (t in (nyr+1):(nyr+K) ){
for (s in 1:nsite){
numer.perc.treat[sc, t, s] <- # numerator
step( NB[sc, t, s]-(num.treated[sc]+1) ) * num.treated[sc] + # num.treated > NB, use num.treated
( 1-step( NB[sc, t, s]-(num.treated[sc]+1) )) * NB[sc, t, s] # num.treated < NB set to NB
denom.perc.treat[sc, t, s] <- # denominator, total nests available for treatment
step( NB[sc, t, s]-(num.treated[sc]+1) ) * NB[sc, t, s] + # num.treated > NB, use num.treated
( 1-step( NB[sc, t, s]-(num.treated[sc]+1) )) * numer.perc.treat[sc, t, s]
denom.perc.treat2[sc, t, s] <- equals(denom.perc.treat[sc, t, s], 0) +
(1- equals(denom.perc.treat[sc, t, s], 0)) *denom.perc.treat[sc, t, s]
perc.treat[sc, t, s] <- numer.perc.treat[sc, t, s] /
denom.perc.treat2[sc, t, s]
log(mn.prod[sc, t, s]) <- lmu.prod[s] +
gamma*perc.treat[sc, t, s] + # treat.pair set to one here.
eta[9, s, t]
}} } # s t sc
################################
# Likelihood for counts
################################
# Abundance for year=1
for (v in 1:7){
for (s in 1:nsite){
# Subtract one to allow dcat to include zero
N[1, v, 1, s] <- N2[1, v, 1, s] - 1
N2[1, v, 1, s] ~ dcat(pPrior[v, 1:s.end[v,s], s])
# Assign estimates for years with data to all scenarios
for (sc in 2:SC){
N[sc, v, 1, s] <- N[1, v, 1, s]
}
}} # s t
# Abundance for years > 1
for (t in 1:(nyr-1)){
for (s in 1:nsite){
# Number of wild born juvs
N[1, 1, t+1, s] ~ dpois( (NFY[1, t, s]*mn.phiFY[1, t, s]*mn.psiFYB[1, t, s] + # first year breeders
NF[1, t, s]*mn.phiA[1, t, s]*mn.psiAB[1, t, s] + # nonbreeders to breeders
NB[1, t, s]*mn.phiB[1, t, s]*(1-mn.psiBA[1, t, s])) # breeders remaining
*(mn.prod[1, t+1, s]/2) ) # end Poisson
# Abundance of nonbreeders
N[1, 2, t+1, s] ~ dbin(mn.phiFY[1, t, s]*(1-mn.psiFYB[1, t, s]), NFY[1, t, s]) # Nestlings to second year nonbreeders
N[1, 3, t+1, s] ~ dbin(mn.phiA[1, t, s]*(1-mn.psiAB[1, t, s]), NF[1, t, s]) # Nonbreeders to nonbreeders
N[1, 4, t+1, s] ~ dbin(mn.phiB[1, t, s]*mn.psiBA[1, t, s], NB[1, t, s]) # Breeders to nonbreeders
# Abundance of breeders
N[1, 5, t+1, s] ~ dbin(mn.phiFY[1, t, s]*mn.psiFYB[1, t, s], NFY[1, t, s]) # Nestlings to second year breeders
N[1, 6, t+1, s] ~ dbin(mn.phiA[1, t, s]*mn.psiAB[1, t, s], NF[1, t, s]) # Nonbreeder to breeder
N[1, 7, t+1, s] ~ dbin(mn.phiB[1, t, s]*(1-mn.psiBA[1, t, s]), NB[1, t, s]) # Breeder to breeder
# Assign estimates for years with data to
# all scenarios
for (v in 1:7){
for (sc in 2:SC){
N[sc, v, t+1, s] <- N[1, v, t+1, s]
} } # sc
}} # s t
# Future projections- population model
for (sc in 1:SC){
for (t in nyr:(nyr+K-1)){
for (s in 1:nsite){
# Number of wild born juvs
N[sc, 1, t+1, s] ~ dpois( (NFY[sc, t, s]*mn.phiFY[sc, t, s]*mn.psiFYB[sc, t, s] + # first year breeders
NF[sc, t, s]*mn.phiA[sc, t, s]*mn.psiAB[sc, t, s] + # nonbreeders to breeders
NB[sc, t, s]*mn.phiB[sc, t, s]*(1-mn.psiBA[sc, t, s])) # breeders remaining
*(mn.prod[sc, t+1, s]/2) ) # end Poisson
# Abundance of nonbreeders
## Second year nonbreeders
N[sc, 2, t+1, s] ~ dbin(mn.phiFY[sc, t, s]*(1-mn.psiFYB[sc, t, s]), NFY[sc, t, s]) # Nestlings to second year nonbreeders
## Adult nonbreeders
N[sc, 3, t+1, s] ~ dbin(mn.phiA[sc, t, s]*(1-mn.psiAB[sc, t, s]), NF[sc, t, s]) # Nonbreeders to nonbreeders
N[sc, 4, t+1, s] ~ dbin(mn.phiB[sc, t, s]*mn.psiBA[sc, t, s], NB[sc, t, s]) # Breeders to nonbreeders
# Abundance of breeders
## Second year breeders
N[sc, 5, t+1, s] ~ dbin(mn.phiFY[sc, t, s]*mn.psiFYB[sc, t, s], NFY[sc, t, s]) # Nestlings to second year breeders
## Adult breeders
N[sc, 6, t+1, s] ~ dbin(mn.phiA[sc, t, s]*mn.psiAB[sc, t, s], NF[sc, t, s]) # Nonbreeder to breeder
N[sc, 7, t+1, s] ~ dbin(mn.phiB[sc, t, s]*(1-mn.psiBA[sc, t, s]), NB[sc, t, s]) # Breeder to breeder
}} # s t
} # sc
# Count likelihoods, state-space model, and observation process
for (t in 1:nyr){
for (s in 1:nsite){
num.hacked[1, t, s] <- step( N[1, 1, t, 1] + hacked.counts[1, t, 1] ) * hacked.counts[1, t, s] +
(1- step( N[1, 1, t, 1] + hacked.counts[1, t, 1] )) * N[1, 1, t, 1] *
(-1*equals(s,1) + equals(s,2)) # change to
constraint_data[t, s] ~ dconstraint( (N[1, 1, t, s] + hacked.counts[1, t, s]) >= 0 ) # Transfers translocated first-year females
NFY[1, t, s] <- N[1, 1, t, s] + hacked.counts[1, t, s] # Transfers translocated first-year females
NF[1, t, s] <- sum(N[1, 2:4, t, s]) # number of adult nonbreeder females
NB[1, t, s] <- sum(N[1, 5:7, t, s]) # number of adult breeder females
NAD[1, t, s] <- NF[1, t, s] + NB[1, t, s] # number of adults
Ntot[1, t, s] <- sum(N[1, 1:7, t, s]) # total number of females
countsAdults[t, s] ~ dpois(lamAD[1, t, s]) # adult females
# constrain N1+hacked.counts to be >=0
log(lamAD[1, t, s]) <- log(NAD[1, t, s]) + deltas[1]*effort2[t, s] + deltas[2]*effort2[t, s]^2
log(lamFY[1, t, s]) <- log(N[1, 1, t, s]) + deltas[3]*effort2[t, s] + deltas[4]*effort2[t, s]^2
for (sc in 2:SC){
lamAD[sc, t, s] <- lamAD[1, t, s]
lamFY[sc, t, s] <- lamFY[1, t, s]
NFY[sc, t, s] <- NFY[1, t, s]
NF[sc, t, s] <- NF[1, t, s]
NB[sc, t, s] <- NB[1, t, s]
NAD[sc, t, s] <- NAD[1, t, s]
Ntot[sc, t, s] <- Ntot[1, t, s]
num.hacked[sc, t, s] <- num.hacked[1, t, s]
} # sc
}# s
# First-years at different sites have different distributions
# for better model fit
countsFY[t, 1] ~ dpois(lamFY[1, t, 1]) # doesn't have any zeroes so poisson
countsFY[t, 2] ~ dnegbin(pp[t], r) # first year females, includes translocated/hacked
pp[t] <- r/(r+(lamFY[1, t, 2] ))
} # t
r ~ dexp(0.05)
# Future projections of counts
for(sc in 1:SC){
for (t in (nyr+1):(nyr+K)){
for (s in 1:nsite){
# constrain N1+hacked.counts to be >=0, and allow it to shrink as the population shrinks
# If LH has fewer birds than translocations, none are translocated
num.hacked[sc, t, s] <- step( N[sc, 1, t, 1] + hacked.counts[sc, t, 1] ) * hacked.counts[sc, t, s] +
(1- step( N[sc, 1, t, 1] + hacked.counts[sc, t, 1] )) * N[sc, 1, t, 1] *
(-1*equals(s,1) + equals(s,2)) # change to negative value for LH when using
NFY[sc, t, s] <- N[sc, 1, t, s] + num.hacked[sc, t, s]
NF[sc, t, s] <- sum(N[sc, 2:4, t, s]) # number of adult nonbreeder females
NB[sc, t, s] <- sum(N[sc, 5:7, t, s]) # number of adult breeder females
NAD[sc, t, s] <- NF[sc, t, s] + NB[sc, t, s] # number of adults
Ntot[sc, t, s] <- sum(N[sc, 1:7, t, s]) # total number of females
}# s
} }# t sc
################################
# Likelihood for survival
################################
# Calculate yearly averages for integration
for (t in 1:(nyr-1)){
for (s in 1:nsite){
for (xxxx in 1:surv.end[t,s]){
# need to reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.surv is a matrix of indices for each site
phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
} # xxxx
mn.phiFY[1, t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] )
mn.phiA[1, t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
mn.phiB[1, t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
mn.psiFYB[1, t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
mn.psiAB[1, t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
mn.psiBA[1, t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
mn.pA[1, t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
mn.pB[1, t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
for (sc in 2:SC){
mn.phiFY[sc, t, s] <- mn.phiFY[1, t, s]
mn.phiA[sc, t, s] <- mn.phiA[1, t, s]
mn.phiB[sc, t, s] <- mn.phiB[1, t, s]
mn.psiFYB[sc, t, s] <- mn.psiFYB[1, t, s]
mn.psiAB[sc, t, s] <- mn.psiAB[1, t, s]
mn.psiBA[sc, t, s] <- mn.psiBA[1, t, s]
mn.pA[sc, t, s] <- mn.pA[1, t, s]
mn.pB[sc, t, s] <- mn.pB[1, t, s]
}
}}
for (i in 1:nind){
for (t in 1:(nyr-1)){
#Survival
logit(phiFY[i,t]) <- eta[1, site[i,t],t] +
lmus[1, site[i,t]] + betas[1]*hacked[i] # first year
logit(phiA[i,t]) <- eta[2, site[i,t],t] +
lmus[2, site[i,t]] #+ betas[2]*hacked[i] # nonbreeder
logit(phiB[i,t]) <- eta[3, site[i,t],t] +
lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
#Recruitment
logit(psiFYB[i,t]) <- eta[4, site[i,t],t] +
lmus[4, site[i,t]] #+ betas[4]*hacked[i] # first year to breeder
logit(psiAB[i,t]) <- eta[5, site[i,t],t] +
lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
logit(psiBA[i,t]) <-
lmus[6, site[i,t]] #+ betas[6]*hacked[i] # breeder to nonbreeder
#Re-encounter
logit(pA[i,t]) <- eta[7, site[i,t],t] +
lmus[7, site[i,t]] + betas[7]*hacked[i] +
deltas[5]*effort2[t, site[i,t]] + deltas[6]*effort2[t, site[i,t]]^2 # resight of nonbreeders
logit(pB[i,t]) <- eta[8, site[i,t],t] +
lmus[8, site[i,t]] + #betas[8]*hacked[i] +
deltas[7]*effort2[t, site[i,t]] + deltas[8]*effort2[t, site[i,t]]^2 # resight of breeders
}#t
}#i
# Future projections of survival, recruitment, and detection
for(sc in 1:SC){
for (t in 7:(nyr+K)){
for (s in 1:nsite){
# FYs we calculate the percent hacked for each year.
perc.hacked[sc, t, s] <- ( num.hacked[sc, t, s] / (NFY[sc, t, s]+0.0001) ) *
equals(s,2) # set LH to zero
# For adults we average over the previous 5 years as an approximation
perc.hacked.5yr[sc, t, s] <- ( sum( num.hacked[sc, (t-6):(t-1), s] ) /
(sum( NFY[sc, (t-6):(t-1), s] )+0.0001) )*
equals(s,2) # set LH to zero
}}}
cap <- 0.99 # cap on maximum survival
for(sc in 1:SC){
for (t in nyr:(nyr+K)){
for (s in 1:nsite){
# enforce the survival cap
mn.phiFY[sc, t, s] <- step(cap - mn.phiFY1[sc, t, s]) * mn.phiFY1[sc, t, s] +
(1 - step(cap - mn.phiFY1[sc, t, s])) * cap
mn.phiA[sc, t, s] <- step(cap - mn.phiA1[sc, t, s]) * mn.phiA1[sc, t, s] +
(1 - step(cap - mn.phiA1[sc, t, s])) * cap
mn.phiB[sc, t, s] <- step(cap - mn.phiB1[sc, t, s]) * mn.phiB1[sc, t, s] +
(1 - step(cap - mn.phiB1[sc, t, s])) * cap
# Survival
mn.phiFY1[sc, t, s] <- ilogit( eta[1, s, t] +
lmus[1, s] + betas[1]*perc.hacked[sc, t, s] ) + # change perc.hacked to zero at LH bc no birds translocated there
surv.diff[sc, 1, s]
mn.phiA1[sc, t, s] <- ilogit( eta[2, s, t] +
lmus[2, s] ) + surv.diff[sc, 2, s]
mn.phiB1[sc, t, s] <- ilogit( eta[3, s, t] +
lmus[3, s] + betas[3]*perc.hacked.5yr[sc, t, s] ) + # change perc.hacked to zero at LH bc no birds translocated there
surv.diff[sc, 3, s]
#Recruitment
logit( mn.psiFYB[sc, t, s] ) <- eta[4, s, t] +
lmus[4, s] # first year to breeder
logit( mn.psiAB[sc, t, s] ) <- eta[5, s, t] + # nonbreeder to breeder
lmus[5, s] + betas[5]*perc.hacked.5yr[sc, t, s] # change perc.hacked to zero at LH bc no birds translocated there
logit( mn.psiBA[sc, t, s] ) <-
lmus[6, s] # breeder to nonbreeder
#Re-encounter
logit( mn.pA[sc, t, s] ) <- eta[7, s, t] +
lmus[7, s] + betas[7]*perc.hacked.5yr[sc, t, s] # resight of nonbreeders
logit( mn.pB[sc, t, s] ) <- eta[8, s, t] +
lmus[8, s] # resight of breeders
} # s
} # t
} # sc
# Define state-transition and observation matrices
for (i in 1:nind){
# Define probabilities of state S(t+1) given S(t)
for (t in first[i]:(nyr-1)){
ps[1,i,t,1] <- 0
ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t])
ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
ps[1,i,t,4] <- (1-phiFY[i,t])
ps[2,i,t,1] <- 0
ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
ps[2,i,t,4] <- (1-phiA[i,t])
ps[3,i,t,1] <- 0
ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
ps[3,i,t,4] <- (1-phiB[i,t])
ps[4,i,t,1] <- 0
ps[4,i,t,2] <- 0
ps[4,i,t,3] <- 0
ps[4,i,t,4] <- 1
# Define probabilities of O(t) given S(t)
po[1,i,t,1] <- 1
po[1,i,t,2] <- 0
po[1,i,t,3] <- 0
po[1,i,t,4] <- 0
po[2,i,t,1] <- 0
po[2,i,t,2] <- pA[i,t]
po[2,i,t,3] <- 0
po[2,i,t,4] <- (1-pA[i,t])
po[3,i,t,1] <- 0
po[3,i,t,2] <- 0
po[3,i,t,3] <- pB[i,t]
po[3,i,t,4] <- (1-pB[i,t])
po[4,i,t,1] <- 0
po[4,i,t,2] <- 0
po[4,i,t,3] <- 0
po[4,i,t,4] <- 1
} #t
} #i
# Likelihood
for (i in 1:nind){
# Define latent state at first capture
z[i,first[i]] <- y[i,first[i]]
for (t in (first[i]+1):nyr){
# State process: draw S(t) given S(t-1)
z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
# Observation process: draw O(t) given S(t)
y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
} #t
} #i
} )
#**********************
#* Function to run model in NIMBLE
#**********************
run_pva <- function(info, datl, constl, code){
library('nimble')
library('coda')
params <- c(
# pop growth
"lambda",
# prod
"lmu.prod", "gamma", "rr", "mn.prod",
# survival
"mus", "lmus", "betas", "deltas",
# abundance
"NB", "NF", "NFY", "N", "NAD", "Ntot",
"r",
# error terms
"eta", "sds", "Ustar", "R",
# yearly summaries
'mn.phiFY', 'mn.phiA', 'mn.phiB',
'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
'mn.pA', 'mn.pB',
# pva
"extinct", "extinctAD", "extinctB"
)
#n.chains=1; n.thin=1; n.iter=50; n.burnin=25
n.chains=1; n.thin=50; n.iter=100000; n.burnin=50000
mod <- nimbleModel(code,
constants = constl,
data = datl,
inits = info$inits,
buildDerivs = FALSE,
calculate=T
)
cmod <- compileNimble(mod, showCompilerOutput = TRUE)
confhmc <- configureMCMC(mod)
confhmc$setMonitors(params)
hmc <- buildMCMC(confhmc)
chmc <- compileNimble(hmc, project = mod,
resetFunctions = TRUE,
showCompilerOutput = TRUE)
post <- runMCMC(chmc,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
thin = n.thin,
samplesAsCodaMCMC = T,
setSeed = info$seed)
return(post)
} # run_pva function end
#*****************
#* Run chains in parallel
#*****************
# WARNING: This takes 12 days to run parallel on an HPC
this_cluster <- makeCluster(cpus)
post <- parLapply(cl = this_cluster,
X = par_info_pva[1:cpus],
fun = run_pva,
datl = datl,
constl = constl,
code = mycode)
stopCluster(this_cluster)
# save(post, mycode, seeds, cpus,
# file="/bsuscratch/brianrolek/riha_ipm/outputs/pva_shortrun.rdata")2.2 PVA Projections
2.2.1 Local Extinction Probability
#* Plot PVA results
#*******************
library ('viridis')
library ("tidybayes")
library ('coda')
library ('MCMCvis')
library ('reshape2')
library('ggplot2')
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\pva_shortrun100k.rdata")
out <- lapply(post[-5], as.mcmc)
outp <- MCMCpstr(out, type="chains")
scen <- data.frame('Scenarios' = 1:45,
'Number translocated' = rep(c(0,5,10), each=15),
'Number of nests treated' = rep( rep(c(0, 15, 30, 45, 100), each=3), 3),
'Mortality reduction' = rep(c('100% Mortality', '80% Mortality', '60% Mortality'), 15) )
#*********************
#* Line plot of extinction probability over time
#*********************
pe <- apply(outp$extinct, c(1,2,3), mean)
lpe <- melt(pe)
colnames(lpe) <- c('Scenarios', 'Time', 'Site2', 'Extinct')
pe.df <- merge(lpe, scen, by='Scenarios', all.x=TRUE)
pe.df$Year <- pe.df$Time+2023
pe.df$Site <- ifelse(pe.df$Site2==1, "Los Haitises", "Punta Cana")
pe.df$Mortality.reduction <- factor(pe.df$Mortality.reduction,
levels= c('100% Mortality', '80% Mortality', '60% Mortality'))
p6 <- ggplot(pe.df, aes(x=Year, y=Extinct, group=Scenarios,
color = factor(Number.translocated),
linetype = factor(Number.of.nests.treated))) +
geom_line( linewidth = 1) +
facet_grid(rows= vars(Site, factor(Number.translocated)),
cols= vars(Mortality.reduction)) +
scale_color_viridis_d(option="viridis", direction=1,
name="Number translocated") +
scale_linetype(name="Number of nests treated") +
scale_x_continuous( name="Future year",
breaks=c(2020, 2030, 2040, 2050, 2060, 2070),
labels=c('', '2030', '', '2050', '', '2070'),
minor_breaks=NULL) +
scale_y_continuous( name="Extinction probability",
breaks=seq(0, 1, by=0.1),
labels=c('0', '', '', '', '', '0.5', '', '', '', '', '1'),
limits=c(0,1),
minor_breaks=NULL) +
theme_minimal()
p6
ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Extinction_lines.tiff",
plot=p6,
device="tiff",
width=6, height=6, dpi=300)#********************
#* PVA heatmap of extinction probability projected 50 years
#********************
pe.df2 <- pe.df[pe.df$Time==50, ]
p7 <- ggplot(pe.df2, aes(x = as.factor(Number.translocated),
y = as.factor(Number.of.nests.treated))) +
geom_tile( aes(fill = Extinct ) ) +
geom_text(aes(label=round(Extinct,2)), color="gray70") +
scale_fill_viridis_c(option="plasma", direction=1,
name = "Extinction\nprobability after\n50 years") +
facet_wrap(~ Site + factor(Mortality.reduction,
levels=c("100% Mortality", "80% Mortality", "60% Mortality"))) +
scale_y_discrete(breaks=c(0,15,30,45,100),
name="Number of territories treated") +
scale_x_discrete(name="Number of females translocated") +
theme_minimal( ) + theme(strip.background = element_rect(size = 0.5))
p7 
ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Extinction_50yrHeatmap.tiff",
plot=p7,
device="tiff",
width=6, height=4, dpi=300)
pe.df3 <- pe.df2[!is.na(pe.df2$Extinct),]
pe.df3 <- pe.df3[rev(order(pe.df3$Site2, pe.df3$Extinct)),]
pe.df32.2.3 Abundance
#***************************
#* Plot total Abundance
#***************************
lnt1 <- outp$Ntot[,1:13, , ] |> melt()
colnames(lnt1) <- c("Scenarios", "Time", "Site2", "Iter", "Abund")
lnt1$Year <- lnt1$Time + 2010
lnt1$Site <- factor(ifelse(lnt1$Site2==1, 'Los Haitises', 'Punta Cana'))
nt.df <- merge(lnt1, scen, by='Scenarios', all.x=TRUE)
nt.df$Mortality.reduction <- factor(nt.df$Mortality.reduction,
levels=c('100% Mortality','80% Mortality', '60% Mortality' ))
# calculate future means
lnt2 <- apply(outp$Ntot[,13:63, , ], c(1, 2, 3), mean) |> melt()
colnames(lnt2) <- c("Scenarios", "Time", "Site2", "Abund")
lnt2$Year <- lnt2$Time + 2022
lnt2$Site <- factor(ifelse(lnt2$Site2==1, 'Los Haitises', 'Punta Cana'))
nt.df2 <- merge(lnt2, scen, by='Scenarios', all.x=TRUE)
nt.df2$Mortality.reduction <- factor(nt.df2$Mortality.reduction,
levels=c('100% Mortality','80% Mortality', '60% Mortality' ))
p8 <- ggplot() +
stat_lineribbon(data= nt.df,
aes(x=Year, y=Abund, group=Scenarios),
fill='gray60',
linewidth=0.5,
point_interval = 'median_hdci', .width = 0.95,
na.rm=TRUE) +
geom_line(data= nt.df2,
aes(x=Year, y=Abund, #group=Scenarios,
color = factor(Number.translocated),
linetype = factor(Number.of.nests.treated) ),
linewidth=0.5) +
facet_grid(rows = vars(Site, Number.translocated), cols = vars(Mortality.reduction),
scales ='free_y') +
ylab('Number of females') +
scale_color_viridis_d(option="viridis", direction=1,
name="Number translocated") +
scale_linetype(name="Number of nests treated") +
theme_minimal()
p8
ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Abundance_projections.tiff",
plot=p8,
device="tiff",
width=8, height=6, dpi=300)2.2.4 Percent Change in Abundance
#******************************
#* Future percent change of population
#* between 2023 and 2073
#******************************
# calculate percent change and make a database for plotting
pcl <- ((( outp$Ntot[,63,,]/outp$Ntot[,13,,] )-1)*100) |> melt()
colnames(pcl) <- c("Scenario", "Site.num", "Iter", "PC")
pc.mns <- tapply(pcl$PC, list(pcl$Scenario, pcl$Site.num), mean) |> melt()
colnames(pc.mns) <- c("Scenario", "Site.num", "PC")
pc.df <- merge(pc.mns, scen, by.x="Scenario", by.y="Scenarios")
pc.df$Site <- ifelse(pc.df$Site.num==1, "Los Haitises", "Punta Cana")
# scale color so that percent change=0 is the middle color (gray),
# -100 is the extreme lower bound and max is the extreme upper bound
pc.df$pc.sc <- NA
for (i in 1:length(pc.df$PC)){
if(pc.df$PC[i]<=0){
pc.df$pc.sc[i] <- (pc.df$PC[i]-0)/sd(pc.df$PC[pc.df$PC<=0])}
else{ pc.df$pc.sc[i] <- (pc.df$PC[i]-0)/sd(pc.df$PC[pc.df$PC>0]) }
}
labs <- c(-90, 0, 700, 1400)
labs.sc <- c( (-90-0)/sd(pc.df$PC[pc.df$PC<=0]),
(0-0)/sd(pc.df$PC[pc.df$PC<=0]),
(700-0)/sd(pc.df$PC[pc.df$PC>0]),
(1400-0)/sd(pc.df$PC[pc.df$PC>0]))
# plot
p11 <- ggplot(pc.df, aes(x = as.factor(Number.translocated),
y = as.factor(Number.of.nests.treated))) +
geom_tile( aes(fill = pc.sc ) ) +
geom_text(aes(label=round(PC,0)), color="gray50") +
scale_fill_viridis_c(option="turbo", direction=-1,
name = "Percent change\nafter 50 years",
breaks = labs.sc, labels = labs,
begin=0, end=1) +
facet_wrap(~ Site + factor(Mortality.reduction,
levels=c("100% Mortality", "80% Mortality", "60% Mortality"))) +
scale_y_discrete(breaks=c(0, 15, 30, 45, 100),
name="Number of territories treated") +
scale_x_discrete(name="Number of females translocated") +
theme_minimal( ) + theme(strip.background = element_rect(size = 0.5))
p11 
ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Abundance_50yrHeatmap.tiff",
plot=p11,
device="tiff",
width=6, height=4, dpi=300)2.2.5 Probability of Decline in Abundance
#**********************************
# Probability of Future Population Decline
#**********************************
# p.decline <- apply( (outp$Ntot[,63,,]-outp$Ntot[,13,,])<0, c(1,2), mean)
# pdl <- melt(p.decline)
# colnames(pdl) <- c("Scenario", "Site.num", "pd")
# pdl$Site <- ifelse(pdl$Site.num==1, "Los Haitises", "Punta Cana")
# pdl <- merge(pdl, scen, by.x="Scenario", by.y="Scenarios")
p.decline <- tapply(pcl$PC<0, list(pcl$Scenario, pcl$Site.num), mean)
pdl <- melt(p.decline)
colnames(pdl) <- c("Scenario", "Site.num", "pd")
pdl$Site <- ifelse(pdl$Site.num==1, "Los Haitises", "Punta Cana")
pd <- merge(pdl, scen, by.x="Scenario", by.y="Scenarios")
p12 <- ggplot(pd, aes(x = as.factor(Number.translocated),
y = as.factor(Number.of.nests.treated))) +
geom_tile( aes(fill = pd ) ) +
geom_text(aes(label=round(pd,2)), color="gray50") +
scale_fill_viridis_c(option="magma", direction=1,
name = "Prob. of decline\nafter 50 years",
breaks = c(0,0.25,0.50,0.75,1.0),
labels = c(0,0.25,0.50,0.75,1.0),
begin=0, end=1) +
facet_wrap(~ Site + factor(Mortality.reduction,
levels=c("100% Mortality", "80% Mortality", "60% Mortality"))) +
scale_y_discrete(breaks=c(0, 15, 30, 45, 100),
name="Number of territories treated") +
scale_x_discrete(name="Number of females translocated") +
theme_minimal( ) + theme(strip.background = element_rect(size = 0.5))
p12 
ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//ProbDecline_50yrHeatmap.tiff",
plot=p12,
device="tiff",
width=6, height=4, dpi=300)2.2.6 Time to Extinction
#*****************************
# Time to extinction
#*****************************
# conditional on extinction threshold (D=X)
D <- 1
abund <- outp$Ntot[,14:(13+50),,]
dimnames(abund) <- list('Scenarios'=1:45,
'Time'= 2024:(2024+49),
'Site'=c("Los Haitises", "Punta Cana"),
'Iter'=1: dim(outp$extinct)[[4]])
qextinct <- abund == D
ext <- qextinct[,50,,] <= D
lte <- melt( qextinct )
colnames(lte)[5] <- "Extinct"
le <- melt( ext )
colnames(le)[4] <- "EverExtinct"
lte2 <- merge(lte, le, by=c("Scenarios", "Site", "Iter") )
lte3<- lte2[lte2$EverExtinct==TRUE, ]
lte4 <- lte3[lte3$Extinct==TRUE, ]
lte5 <- tapply(lte4$Time-2023, list(lte4$Scenarios, lte4$Site, lte4$Iter), min)
lte6 <- melt(lte5)
lte7 <- lte6[!is.na(lte6$value), ]
colnames(lte7) <- c('Scenarios', 'Site', 'Iter', 'TTE')
te <- table(lte7$TTE, lte7$Scenarios, lte7$Site) |> melt()
colnames(te) <- c('TE', 'Scen', 'Site', 'Count')
te$Scenarios <- factor(paste0("Scenario ", te$Scen),
levels=paste0('Scenario ', 1:45))
teLH <- te[te$Site=="Los Haitises",]
tePC <- te[te$Site=="Punta Cana",]
p9 <- ggplot() +
geom_col(data = teLH, aes(x= TE, y=Count)) +
facet_wrap(vars(Scenarios),
nrow=9, ncol=5, scales="free_y", shrink=FALSE)
p10 <- ggplot() +
geom_col(data = tePC, aes(x= TE, y=Count)) +
facet_wrap(vars(Scenarios),
nrow=9, ncol=5, scales="free_y", shrink=FALSE)
p9
p10
3. Additional IPM plots
3.1 Model Diagnostics
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_longrun.rdata")
load("data/data.rdata")
library ('MCMCvis')
library ('coda')
library ('ggplot2')
library('reshape2')
library('bayestestR')
out <- lapply(post, as.mcmc)
# Identify chains with NAs that
# failed to initialize
NAlist <- c()
for (i in 1:length(out)){
NAlist[i] <- any (is.na(out[[i]][,1:286]) | out[[i]][,1:286]<0)
}
# Subset chains to those with good initial values
out <- out[!NAlist]
out <- out[1:4] # omit chain 5 bc one parameter causes a convergence failure
outp <- MCMCpstr(out, type="chains")
!NAlist## [1] TRUE TRUE TRUE TRUE TRUE# default settings for plots
plt <- function(object, params,...) {
MCMCplot(object=out,
params=params,
guide_axis=TRUE,
HPD=TRUE, ci=c(80, 95), horiz=FALSE,
#ylim=c(-10,10),
...)
}3.1.1 Check Goodness-of-fit
Goodness-of-fit tests provide evidence that statistical distributions adequately describe the data. Here we test fit for brood size and counts. A Bayesian p-value nearest to 0.5 suggests good fitting statistical distributions, while values near 1 or 0 suggest poor fit.
# Goodness of fit check
fit.check <- function(out, ratio=FALSE,
name.rep="f.dmape.rep",
name.obs="f.dmape.obs",
jit=100,
ind=1,
lab=""){
par(mfrow=c(1,1))
# plot mean absolute percentage error
samps <- MCMCpstr(out, "all", type="chains")
rep <- samps[name.rep][[1]][ind,]
obs <- samps[name.obs][[1]][ind,]
mx <- max(c(rep, obs))
mn <- min(c(rep, obs))
plot(jitter(obs, amount=jit),
jitter(rep, amount=jit),
main=paste0("Mean absolute percentage error\n",lab),
ylab="Discrepancy replicate values",
xlab="Discrepancy observed values",
xlim=c(mn,mx), ylim=c(mn,mx),
pch=16, cex=0.5, col="gray10")
curve(1*x, from=mn, to=mx, add=T, lty=2, lwd=2, col="blue")
bp1 <- round(mean(rep > obs),2)
loc <- ifelse(bp1 < 0.5, "topleft", "bottomright")
legend(loc, legend=bquote(p[B]==.(bp1)), bty="n", cex=2)
if (ratio==TRUE){
t.rep <- samps["tvm.rep"][[1]][ind,]
t.obs <- samps["tvm.obs"][[1]][ind,]
# plot variance/mean ratio
hist(t.rep, nclass=50,
xlab="variance/mean ", main=NA, axes=FALSE)
abline(v=t.obs, col="red")
axis(1); axis(2)
}
return(list('Bayesian p-value'=bp1))
}
# Counts of adults (nonbreeders+breeders)
# tiff(height=4, width=5, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//GOF_adults.tiff")
fit.check(out, ratio=F,
name.rep="dmape.rep",
name.obs="dmape.obs",
ind=1,
lab="Adults(Breeder+Nonbreeder)- Poisson", jit=300)
## $`Bayesian p-value`
## [1] 0.39# dev.off()
# Counts of first-year fledglings, ind=2
# Poisson failed fit test bp=0
# So we assigned negative binomial only for Punta Cana.
# Los Haitises excluded from neg bin because no zeroes in counts.
# tiff(height=4, width=5, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//GOF_FY.tiff")
fit.check(out, ratio=F,
name.rep="dmape.rep",
name.obs="dmape.obs",
ind=2,
lab="First-year counts\nNeg binomial-Poisson", jit=300)
## $`Bayesian p-value`
## [1] 0.49# dev.off()
# Productivity
# tiff(height=4, width=5, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//GOF_prod.tiff")
fit.check(out, ratio=F,
name.rep="f.dmape.rep",
name.obs="f.dmape.obs",
ind=1,
lab="Productivity-Neg binomial", jit=300)
## $`Bayesian p-value`
## [1] 0.81# dev.off()3.1.2 Examine Traceplots, Compare Posteriors with Priors
Traceplots provide evidence that models adequately converged.
MCMCtrace(out, pdf=FALSE, params= "mus", Rhat=TRUE, priors=runif(20000, 0, 1), post_zm=FALSE)
MCMCtrace(out, pdf=FALSE, params= "lmu.prod", Rhat=TRUE, priors=rnorm(20000, 0, 5), post_zm=FALSE)
MCMCtrace(out, pdf=FALSE, params= "betas", Rhat=TRUE, priors=rnorm(20000, 0, 10), post_zm=FALSE)
MCMCtrace(out, pdf=FALSE, params= "deltas", Rhat=TRUE, priors=rnorm(20000, 0, 10), post_zm=FALSE)
MCMCtrace(out, pdf=FALSE, params= "gamma", Rhat=TRUE, priors=rnorm(20000, 0, 10), post_zm=FALSE)
MCMCtrace(out, pdf=FALSE, params= "sds", Rhat=TRUE, priors=rexp(20000, 1), post_zm=FALSE)
MCMCtrace(out, pdf=FALSE, params= "r", Rhat=TRUE, priors=rexp(20000, 0.05), post_zm=FALSE)
MCMCtrace(out, pdf=FALSE, params= "rr", Rhat=TRUE, priors=rexp(20000, 0.05), post_zm=FALSE)
3.2 Demographics- Annual Averages
# Demographics- Annual Averages
# Annual averages for integration into the population model
par(mfrow=c(1,1))
labs <- c(paste0("LH ",2011:2023), paste0("PC ",2011:2023))
plt(object=out, params="mn.phiFY", ylim=c(0,1),
main="First-year survival", labels = labs,
xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.phiA", ylim=c(0,1),
main="Adult nonbreeder", labels = labs,
xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.phiB", ylim=c(0,1),
main="Breeder", labels = labs,
xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.psiFYB", ylim=c(0,1),
main="First-year to breeder", labels = labs,
xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.psiAB", ylim=c(0,1),
main="Adult nonbreeder to breeder", labels = labs,
xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.psiBA", ylim=c(0,1),
main="Adult breeder to nonbreeder", labels = labs,
xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.pA", ylim=c(0,1),
main="Nonbreeder", labels = labs,
xlab = "Year", ylab= "Detection")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.pB", ylim=c(0,1),
main="Breeder", labels = labs,
xlab = "Year", ylab= "Detection")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.prod",
main="", labels=labs,
xlab = "Year", ylab= "Productivity")
abline(v=13.5, lwd=2)
3.3 Abundance of Finer Population Segments
Plot model estimates of demographic rates. Life Stages are abbreviated as B = breeder, NB = nonbreeder, FY = first year. First-year abundance accounts for translocated birds.
Population dynamics are determined by transitions, These transitions
between stages are abbreviated as the starting life stage to the final
life stage. For example a first-year recruiting to a breeder would be
abbreviated as “FY to B”. I’ll list them here for convenience:
* “FY to NB” is first-year to nonbreeder.
* “NB to NB” is nonbreeder adult to nonbreeder adult.
* “B to NB” is a breeding adult to a nonbreeder adult.
* “FY to B” is first-year to breeder.
* “NB to B” is nonbreeder adult to breeder adult.
* “B to B” is breeder adult to breeder adult.
# Finer population segments
par(mfrow=c(4,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 1, ", ", 1:13, ", 1]"),
main="Los Haitises\nFirst-years fledged",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 2, ", ", 1:13, ", 1]"),
main="Los Haitises\nFY to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 3, ", ", 1:13, ", 1]"),
main="Los Haitises\nNB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 4, ", ", 1:13, ", 1]"),
main="Los Haitises\nB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 5, ", ", 1:13, ", 1]"),
main="Los Haitises\nFY to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 6, ", ", 1:13, ", 1]"),
main="Los Haitises\nNB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 7, ", ", 1:13, ", 1]"),
main="Los Haitises\nB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
par(mfrow=c(4,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 1, ", ", 1:13, ", 2]"),
main="Punta Cana\nFY fledged",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 2, ", ", 1:13, ", 2]"),
main="Punta Cana\nFY to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 3, ", ", 1:13, ", 2]"),
main="Punta Cana\nNB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 4, ", ", 1:13, ", 2]"),
main="Punta Cana\nB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 5, ", ", 1:13, ", 2]"),
main="Punta Cana\nFY to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 6, ", ", 1:13, ", 2]"),
main="Punta Cana\nNB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 7, ", ", 1:13, ", 2]"),
main="Punta Cana\nB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")