Workflow for Snowy Owl Population Trends
Brian W. Rolek
2024-05-23
- 1 Simulations: Can regression models detect population declines when counts have cycles?
- 2 Snowy Owl Analyses
Supplemental materials for: McCabe, R.A., Aarvak, T., Aebischer, A., Bates, K., Bêty, J., Bollache, L., Brinker, D., Driscoll, C., Elliott, K.H., Fitzgerald, G., Fuller, M., Gauthier, G., Gilg, O., Gousy-Leblanc, M., Holt, D., Jacobsen, K.O., Johnson, D., Kulikova, O., Lang, J., Lecomte, N., McClure, C., McDonald, T., Menyushina, I., Miller, E., Morozov, V.V., Øien, I.J., Robillard, A., Rolek, B., Sittler, B., Smith, N., Sokolov, A., Sokolova, N., Solheim, R., Soloviev, M., Stoffel, M., Weidensaul, S., Wiebe, K.L.,  Zazelenchuck, D., and Therrien, J.F.. 2024. Status assessment and conservation priorities for a circumpolar raptor: the snowy owl. Bird Conservation International.
Contact information: rolek.brian@peregrinefund.org
Metadata, data, and scripts used in analyses can be found at https://github.com/The-Peregrine-Fund/Snowy-Owl-Population-Trends.
The full workflow below is visible as a html website at: https://the-peregrine-fund.github.io/Snowy-Owl-Population-Trends/.
A permanent archive and DOI is available at: https://zenodo.org/doi/10.5281/zenodo.11186858
1 Simulations: Can regression models detect population declines when counts have cycles?
Here, we simulate data of an irruptive species using a zero-inflated poisson to induce cycles. Then we analyze it using a negative binomial.
library(nimble)
library(MCMCvis)
library (MASS)
library (pscl)1.1 A function to simulate cyclical time series
sim.ts <- function( period = 4,
ntime = 30,
mn.psi = -1.38, # about 1/5
mn.abund = 10,
beta0 = -0.7,
beta1 = 2,
beta2 = 2,
plot.ts = TRUE){
yr2 <- 1:ntime
yr <- (yr2-(ntime/2+0.5))/(ntime/2-0.5)
period2 <- (yr[2]-yr[1])*period
# decreasing abundance over time
lambda <- exp( log(mn.abund) + beta0*yr )
# irruptive cycles as a zero-inflation parameter
psi <- plogis(mn.psi +
beta1 * sin(2*pi*(yr/period2)) +
beta2 * cos(2*pi*(yr/period2)))
# simulate abundance
z <- rbinom(length(psi), size=1, prob=psi)
mu <- lambda*z
abund <- rpois(length(mu), mu)
if(plot.ts==TRUE){
plot(abund, type="l",
xlab="Year", ylab="Abundance")
}
return(list(datl = list(y = abund),
constl = list(ntime = ntime,
yr = yr,
period=period2*5),
truth = list(mn.psi = mn.psi,
mn.abund = mn.abund,
beta0 = beta0,
beta1 = beta1,
beta2 = beta2)
))
}1.2 Simulate data and implement negative binomial regression
#*******************
#* Frequentist analysis
#* using a
#* negative binomial distribution
#*******************
set.seed(123456)
# Population trends vary
df <- data.frame(est=NA, truth=NA, bias=NA, cov=NA, lci=NA, uci=NA)
war <- rep(FALSE, 5000)
ind <- 1
beta.scenarios <- c(-0.7, -0.35, 0, 0.35, 0.7)
for (j in 1:length(beta.scenarios)){
for (i in 1:1000){
dat <- sim.ts(beta0=beta.scenarios[j], plot.ts=FALSE)
dfm <- data.frame(y=dat$datl$y, yr=dat$constl$yr)
fit <- tryCatch(glm(y ~ 1 + yr, family = negative.binomial(2),
data = dfm),
warning= function(w) { return(TRUE) } )
if (is.logical(fit)){ war[ind] <- fit
ind <- ind+1
next} else{
out <- summary(fit)
uci <- out$coefficients[2,1] + 1.96*out$coefficients[2,2]
lci <- out$coefficients[2,1] - 1.96*out$coefficients[2,2]
estimate <- out$coefficients[2,1]
truth <- dat$truth$beta0
relbias <- ifelse(truth==0,
((estimate+1)-(truth+1))/(truth+1),
(estimate-truth)/truth)
coverage <- ifelse(lci>=truth | uci<=truth, 0, 1)
df[ind, 1:6] <- c(estimate, truth, relbias, coverage, lci, uci)
ind <- ind + 1
}
}}
df1 <- df[war!=TRUE,]
#* Abundance varies
set.seed(123456)
df <- data.frame(est=NA, truth=NA, bias=NA, cov=NA, lci=NA, uci=NA, abund=NA)
war <- rep(FALSE, 4000)
ind <- 1
abund.scenarios <- c(0.5, 1, 2, 3)
for (j in 1:length(abund.scenarios)){
for (i in 1:1000){
dat <- sim.ts(mn.abund=abund.scenarios[j], plot.ts=FALSE)
dfm <- data.frame(y=dat$datl$y, yr=dat$constl$yr)
fit <- tryCatch(glm(y ~ 1 + yr, family = negative.binomial(2),
data = dfm),
warning= function(w) { return(TRUE) } )
if (is.logical(fit)){ war[ind] <- fit
ind <- ind+1
next} else{
out <- summary(fit)
uci <- out$coefficients[2,1] + 1.96*out$coefficients[2,2]
lci <- out$coefficients[2,1] - 1.96*out$coefficients[2,2]
estimate <- out$coefficients[2,1]
truth <- dat$truth$beta0
relbias <- ifelse(truth==0,
((estimate+1)-(truth+1))/(truth+1),
(estimate-truth)/truth)
coverage <- ifelse(lci>=truth |
uci<=truth, 0, 1)
df[ind, 1:7] <- c(estimate, truth, relbias, coverage, lci, uci, abund.scenarios[j])
ind <- ind + 1
}
}}
# need to remove mods that didn't converge.
df2 <- df[war!=TRUE,]
#* Lengths of time vary
set.seed(123456)
df <- data.frame(est=NA, truth=NA, bias=NA, cov=NA, lci=NA, uci=NA, time=NA)
war <- rep(FALSE, 3000)
ind <- 1
time.scenarios <- c(5, 10, 15)
for (j in 1:length(time.scenarios)){
for (i in 1:1000){
dat <- sim.ts(ntime=time.scenarios[j], plot.ts=FALSE)
dfm <- data.frame(y=dat$datl$y, yr=dat$constl$yr)
fit <- tryCatch(glm(y ~ 1 + yr,
family = negative.binomial(2),
data = dfm),
warning= function(w) { return(TRUE) } )
if (is.logical(fit)){ war[ind] <- fit
ind <- ind+1
next} else{
out <- summary(fit)
uci <- out$coefficients[2,1] + 1.96*out$coefficients[2,2]
lci <- out$coefficients[2,1] - 1.96*out$coefficients[2,2]
estimate <- out$coefficients[2,1]
truth <- dat$truth$beta0
relbias <- ifelse(truth==0,
((estimate+1)-(truth+1))/(truth+1),
(estimate-truth)/truth)
coverage <- ifelse(lci>=truth |
uci<=truth, 0, 1)
df[ind, 1:7] <- c(estimate, truth, relbias, coverage, lci, uci, time.scenarios[j])
ind <- ind + 1
}
}}
# need to remove mods that didn't converge
df3 <- df[war!=TRUE,]1.3 Plot simulation results
# plot a sample of time series
dat <- list()
for (i in 1:1000){
dat[[i]] <- sim.ts(beta0=beta.scenarios[1], plot.ts=FALSE)[[1]]$y
}
samps <- sample(1:1000, size=6, replace = FALSE)
df <- do.call(rbind, dat[samps])
library (reshape2)
library(ggplot2)
rownames(df) <- samps
ldf <- melt(df)
ggplot(ldf)+
geom_line(aes(x=Var2, y=value)) +
facet_wrap('Var1') +
ylab("Simulated count") +
xlab("Time")
Fig. S1. Six randomly selected simulations among a total of 1000 for one scenario (mean abundance = 10, duration = 30-years and temporal trend of abundance = -0.7) designed to resemble counts from an irruptive species.
# Plots
par(mfrow=c(1,3))
boxplot(df1$bias~df1$truth,
xlab=expression(paste("Trend (", delta[1],")")),
ylab="Relative bias of trend")
abline(h=0, lty=2, lwd=2)
mtext("A", cex=1.5, side=3, outer=F)
boxplot(df2$bias~df2$abund,
xlab=expression(paste("Mean abundance (", delta[0],")")),
ylab="")
abline(h=0, lty=2, lwd=2)
mtext("B", cex=1.5, side=3, outer=F)
boxplot(df3$bias~df3$time,
xlab="Duration of monitoring",
ylab="")
abline(h=0, lty=2, lwd=2)
mtext("C", cex=1.5, side=3, outer=F)
Fig. S2. Relative bias of population trends estimated by generalized linear models with a negative binomial distribution according to: (A) magnitude of 30-year temporal trend of abundance, (B) mean abundance, and (C) duration of monitoring. The dashed line depicts unbiased estimates. Solid horizontal lines within each box depict medians (50th percentile); boxes depict the 25th and 75th percentiles; whiskers depict the largest and smallest values within 1.5 times the interquartile range; outlier data are depicted as open circles; and a dashed line at no bias (zero) is provided for reference. For each scenario, we simulated 1000 time series data sets.
1.4 Create a table of simulation results
# tables
tab1 <- data.frame(
# different population trends
varied.by= "population trends",
mn.abund=10,
trend= c(-0.7, -0.35, 0, 0.35, 0.7),
duration=30,
mean.rel.bias=tapply(df1$bias, df1$truth, mean),
sd.rel.bias=tapply(df1$bias, df1$truth, sd),
coverage=tapply(df1$cov, df1$truth, mean)
)
tab2 <- data.frame(
#* test different levels of abundance
varied.by= "abundance",
mn.abund=c(0.5, 1, 2, 3),
trend=-0.7,
duration=30,
mean.rel.bias=tapply(df2$bias, df2$abund, mean),
sd.rel.bias=tapply(df2$bias, df2$abund, sd),
coverage=tapply(df2$cov, df2$abund, mean)
)
tab3 <- data.frame(
#* test different lengths of time
varied.by= "duration",
mn.abund=10,
trend=-0.7,
duration=c(5, 10, 15),
mean.rel.bias=tapply(df3$bias, df3$time, mean),
sd.rel.bias=tapply(df3$bias, df3$time, sd),
coverage=tapply(df3$cov, df3$time, mean)
)
tab <- rbind(tab1, tab2, tab3)
rownames(tab) <- NULL
knitr::kable(tab, digits=c(0, 2, 2, 2, 2),
col.names = c("Scenario", "Abundance", "Trend", "Duration",
"Mean relative bias", "SD relative bias", "Coverage"),
row.names=FALSE,
caption="Table S1. Simulation scenarios and results.")| Scenario | Abundance | Trend | Duration | Mean relative bias | SD relative bias | Coverage |
|---|---|---|---|---|---|---|
| population trends | 10.0 | -0.70 | 30 | -0.11 | 1 | 0.98 |
| population trends | 10.0 | -0.35 | 30 | -0.18 | 1 | 0.98 |
| population trends | 10.0 | 0.00 | 30 | 0.03 | 0 | 0.97 |
| population trends | 10.0 | 0.35 | 30 | 0.07 | 1 | 0.98 |
| population trends | 10.0 | 0.70 | 30 | 0.02 | 1 | 0.98 |
| abundance | 0.5 | -0.70 | 30 | 0.02 | 2 | 0.93 |
| abundance | 1.0 | -0.70 | 30 | -0.10 | 1 | 0.96 |
| abundance | 2.0 | -0.70 | 30 | -0.09 | 1 | 0.97 |
| abundance | 3.0 | -0.70 | 30 | -0.15 | 1 | 0.98 |
| duration | 10.0 | -0.70 | 5 | -1.92 | 5 | 0.95 |
| duration | 10.0 | -0.70 | 10 | -0.02 | 2 | 0.96 |
| duration | 10.0 | -0.70 | 15 | 0.25 | 1 | 0.96 |
2 Snowy Owl Analyses
2.1 Posterior predictive checks to assess goodness-of-fit
We tested the goodness-of-fit using posterior predictive checks for three distributions: poisson, negative binomial, and zero-inflated Poisson. First we load the required packages and snowy owl data.
library ('nimble')
library ('nimbleHMC')
library ('MCMCvis')
library ('coda')
library('parallel')
load("data/data.Rdata")
set.seed(5757575)Next, we implemented goodness-of-fit checks for a Poisson, negative binomial, and zero-inflated Poisson distribution. The model using a negative binomial had much greater fit; therefore, we retained this for inference in the manuscript.
2.1.1 Poisson distribution
#***********************
#* Poisson model
#***********************
run <- function(seed, datl, constl){
library('nimble')
library('coda')
library ('nimbleHMC')
code <- nimbleCode(
{
# priors
for (j in 1:nsite){
mu[j] ~ dnorm(0, sd=5) # constrain to reasonable values <exp(5) or <148
beta[j] ~ dnorm(0, sd=5)
}
sigma.time ~ dexp(2)
# likelihood
for (t in 1:ntime){
for (j in 1:nsite){
y[t,j] ~ dpois(lambda[t,j])
# abundance
log(lambda[t,j]) <- mu[j] +
beta[j]*time[t] +
eps.time[t]
} # j
eps.time[t] ~ dnorm(0, sd=sigma.time)
} # t
###################
# Assess GOF of the 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
###################
# GOF for model
for (t in 1:ntime){
for (j in 1:nsite){
c.obs[t,j] <- y[t,j] # observed counts
c.exp[t,j] <- lambda[t,j] # expected counts adult breeder
c.rep[t,j] ~ dpois(lambda[t,j]) # expected counts
# Compute fit statistics, Mean absolute error
dssm.obs[t,j] <- abs( ( (c.obs[t,j]) - (c.exp[t,j]) ) / (c.obs[t,j]+0.001) )
dssm.rep[t,j] <- abs( ( (c.rep[t,j]) - (c.exp[t,j]) ) / (c.rep[t,j]+0.001) )
} } # j # t
dmape.obs <- sum(dssm.obs[1:ntime,1:nsite])
dmape.rep <- sum(dssm.rep[1:ntime,1:nsite])
# variance-mean ratio
tvm.rep <- sd(c.rep[1:ntime,1:nsite])^2/mean(c.rep[1:ntime,1:nsite])
tvm.obs <- sd(y[1:ntime,1:nsite])^2/mean(y[1:ntime,1:nsite])
}
) # end model
params <- c( "sigma.time",
"mu", "beta",
"dssm.obs", "dssm.rep",
"dmape.obs", "dmape.rep",
"tvm.rep", "tvm.obs"
)
# create initial values for missing y data
y.inits <- array(NA, dim(datl$y))
y.inits[is.na(datl$y)] <- rpois(sum(is.na(datl$y)), 2)
inits <- function(){ list (y = y.inits,
sigma.time = rexp(1),
mu = runif(constl$nsite, -2, 2),
beta = runif(constl$nsite, -2, 2)
)}
n.chains=1; n.thin=10; n.iter=20000; n.burnin=10000
mod <- nimbleModel(code, calculate=T, constants = constl,
data = datl, inits = inits(),
buildDerivs = TRUE)
mod$calculate()
cmod <- compileNimble(mod )
confhmc <- configureHMC(mod)
confhmc$setMonitors(params)
hmc <- buildMCMC(confhmc)
chmc <- compileNimble(hmc, project=mod, resetFunctions = TRUE)
post <- runMCMC(chmc,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
thin = n.thin,
samplesAsCodaMCMC = T)
return(post)
}
# approximately 15 minute runtime
this_cluster <- makeCluster(4)
post <- parLapply(cl = this_cluster,
X = 1:4,
fun = run,
dat=datl,
const=constl)
stopCluster(this_cluster)
pois <- list(as.mcmc(post[[1]]),
as.mcmc(post[[2]]),
as.mcmc(post[[3]]),
as.mcmc(post[[4]]))
# save(pois, post, run,
# file="C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\SnowyOwl_HawkMountain\\pois.Rdata")
# Check for convergence
params_pois <-c( "sigma.time",
"mu", "beta")
MCMCtrace(pois, params_pois, pdf=F,
ind = TRUE, Rhat = TRUE, n.eff = TRUE)
par(mfrow=c(1,1))
MCMCplot(object = pois, params = params_pois)2.1.2 Negative binomial distribution
#***********************
#* Negative binomial model
#* Best fitting model
#***********************
run <- function(seed, datl, constl){
library('nimble')
library('coda')
library ('nimbleHMC')
code <- nimbleCode(
{
# priors
# priors
for (j in 1:nsite){
mu[j] ~ dnorm(0, sd=5) # constrain to reasonable values <exp(5) or <148
beta[j] ~ dnorm(0, sd=5)
r[j] ~ dexp(0.1)
} # j
sigma.time ~ dexp(2)
# likelihood
for (t in 1:ntime){
for (j in 1:nsite){
p[t,j] <- r[j]/(r[j]+lambda[t,j])
y[t,j] ~ dnegbin(p[t,j], r[j])
# abundance
log(lambda[t,j]) <- mu[j] +
beta[j]*time[t] +
eps.time[t]
} # j
eps.time[t] ~ dnorm(0, sd=sigma.time)
} # t
###################
# Assess GOF of the 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
###################
# GOF for model
for (t in 1:ntime){
for (j in 1:nsite){
c.obs[t,j] <- y[t,j] # observed counts
c.exp[t,j] <- lambda[t,j] # expected counts adult breeder
c.rep[t,j] ~ dnegbin(p[t,j],r[j]) # expected counts
# Compute fit statistics, Mean absolute error
dssm.obs[t,j] <- abs( ( (c.obs[t,j]) - (c.exp[t,j]) ) / (c.obs[t,j]+0.001) )
dssm.rep[t,j] <- abs( ( (c.rep[t,j]) - (c.exp[t,j]) ) / (c.rep[t,j]+0.001) )
} } # j # t
dmape.obs <- sum(dssm.obs[1:ntime,1:nsite])
dmape.rep <- sum(dssm.rep[1:ntime,1:nsite])
# variance-mean ratio
tvm.rep <- sd(c.rep[1:ntime,1:nsite])^2/mean(c.rep[1:ntime,1:nsite])
tvm.obs <- sd(y[1:ntime,1:nsite])^2/mean(y[1:ntime,1:nsite])
}
) # end model
params <- c( "sigma.time",
"mu", "beta", "r", "p",
"dssm.obs", "dssm.rep",
"dmape.obs", "dmape.rep",
"tvm.rep", "tvm.obs"
)
# create initial values for missing y data
y.inits <- array(NA, dim(datl$y))
y.inits[is.na(datl$y)] <- rpois(sum(is.na(datl$y)), 2)
inits <- function(){ list (y = y.inits,
sigma.time = rexp(1),
mu = runif(constl$nsite, -2, 2),
beta = runif(constl$nsite, -2, 2),
r = runif(constl$nsite)
)}
n.chains=1; n.thin=10; n.iter=20000; n.burnin=10000
mod <- nimbleModel(code, calculate=T, constants = constl,
data = datl, inits = inits(),
buildDerivs = TRUE)
mod$calculate()
cmod <- compileNimble(mod )
confhmc <- configureHMC(mod)
confhmc$setMonitors(params)
hmc <- buildMCMC(confhmc)
chmc <- compileNimble(hmc, project=mod, resetFunctions = TRUE)
post <- runMCMC(chmc,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
thin = n.thin,
samplesAsCodaMCMC = T)
return(post)
}
# approximately 15 minute runtime
this_cluster <- makeCluster(4)
post <- parLapply(cl = this_cluster,
X = 1:4,
fun = run,
dat=datl,
const=constl)
stopCluster(this_cluster)
nb <- list(as.mcmc(post[[1]]),
as.mcmc(post[[2]]),
as.mcmc(post[[3]]),
as.mcmc(post[[4]]))
params_nb <-c( "sigma.time",
"mu", "beta", "r")
# Check for convergence
MCMCtrace(nb, params_nb, pdf=F,
ind = TRUE, Rhat = TRUE, n.eff = TRUE)
par(mfrow=c(1,1))
MCMCplot(object = nb, params = params_nb)
save(out=nb, post=post, run,
file="C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\SnowyOwl_HawkMountain\\nb.Rdata")2.1.3 Zero-inflated Poisson distribution
#***********************
#* Zero-inflated Poisson model
#***********************
run <- function(seed, datl, constl){
library('nimble')
library('coda')
library ('nimbleHMC')
code <- nimbleCode(
{
# priors
for (j in 1:nsite){
mu[j] ~ dnorm(0, sd=5) # constrain to reasonable values <exp(5) or <148
beta[j] ~ dnorm(0, sd=5)
psi[j] ~ dunif(0, 1)
}
sigma.time ~ dexp(2)
# likelihood
for (t in 1:ntime){
for (j in 1:nsite){
y[t,j] ~ dpois(lambda.star[t,j])
lambda.star[t,j] <- lambda[t,j]*z[t,j]
z[t,j] ~ dbern(psi[j])
# abundance
log(lambda[t,j]) <- mu[j] +
beta[j]*time[t] +
eps.time[t]
} # j
eps.time[t] ~ dnorm(0, sd=sigma.time)
} # t
###################
# Assess GOF of the 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
###################
# GOF for model
for (t in 1:ntime){
for (j in 1:nsite){
c.obs[t,j] <- y[t,j] # observed counts
c.exp[t,j] <- z[t,j]*lambda[t,j] # expected counts adult breeder
c.rep[t,j] ~ dpois(lambda[t,j]*psi[j]) # expected counts
# Compute fit statistics, Mean absolute error
dssm.obs[t,j] <- abs( ( (c.obs[t,j]) - (c.exp[t,j]) ) / (c.obs[t,j]+0.001) )
dssm.rep[t,j] <- abs( ( (c.rep[t,j]) - (c.exp[t,j]) ) / (c.rep[t,j]+0.001) )
} } # j # t
dmape.obs <- sum(dssm.obs[1:ntime,1:nsite])
dmape.rep <- sum(dssm.rep[1:ntime,1:nsite])
# variance-mean ratio
tvm.rep <- sd(c.rep[1:ntime,1:nsite])^2/mean(c.rep[1:ntime,1:nsite])
tvm.obs <- sd(y[1:ntime,1:nsite])^2/mean(y[1:ntime,1:nsite])
}
) # end model
params <- c( "sigma.time",
"mu", "beta", "psi",
"dssm.obs", "dssm.rep",
"dmape.obs", "dmape.rep",
"tvm.rep", "tvm.obs"
)
# create initial values for missing y data
y.inits <- array(NA, dim(datl$y))
y.inits[is.na(datl$y)] <- rpois(sum(is.na(datl$y)), 2)
inits <- function(){ list (y = y.inits,
sigma.time = rexp(1),
mu = runif(constl$nsite, -2, 2),
beta = runif(constl$nsite, -2, 2),
psi = runif(constl$nsite)
)}
n.chains=1; n.thin=10; n.iter=20000; n.burnin=10000
mod <- nimbleModel(code, calculate=T, constants = constl,
data = datl, inits = inits(),
buildDerivs = TRUE)
mod$calculate()
cmod <- compileNimble(mod )
# Use MCMC rather than HMC
# because HMC samplers are stuck
confhmc <- configureMCMC(mod)
confhmc$setMonitors(params)
hmc <- buildMCMC(confhmc)
chmc <- compileNimble(hmc, project=mod, resetFunctions = TRUE)
post <- runMCMC(chmc,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
thin = n.thin,
samplesAsCodaMCMC = T)
return(post)
}
# approximately 15 minute runtime
this_cluster <- makeCluster(4)
post <- parLapply(cl = this_cluster,
X = 1:4,
fun = run,
dat=datl,
const=constl)
stopCluster(this_cluster)
zip <- list(as.mcmc(post[[1]]),
as.mcmc(post[[2]]),
as.mcmc(post[[3]]),
as.mcmc(post[[4]]))
# save(zip, post, run,
# file="C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\SnowyOwl_HawkMountain\\zip.Rdata")
params_zip <- c( "sigma.time", "psi",
"mu", "beta")
MCMCtrace(zip, params_zip, pdf=F,
ind = TRUE, Rhat = TRUE, n.eff = TRUE)
par(mfrow=c(1,1))
MCMCplot(object = zip, params = params_zip)2.1.4 A function to assess goodness-of-fit
# Function for posterior predictive checks
# to assess goodness-of-fit
plot.diag <- function(out, ratio=FALSE, lab=""){
par(mfrow=c(1,1))
# plot mean absolute percentage error
samps <- MCMCpstr(out, "all", type="chains")
mx <- max(c(samps$dmape.rep, samps$dmape.obs))
mn <- min(c(samps$dmape.rep, samps$dmape.obs))
plot(jitter(samps$dmape.obs, amount=300),
jitter(samps$dmape.rep, amount=300),
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(samps$dmape.rep > samps$dmape.obs),2)
loc <- ifelse(bp1 < 0.5, "topleft", "bottomright")
legend(loc, legend=bquote(p[B]==.(bp1)), bty="n", cex=2)
if (ratio==TRUE){
# plot variance/mean ratio
hist(samps$tvm.rep, nclass=50,
xlab="variance/mean ", main=NA, axes=FALSE)
abline(v=samps$tvm.obs, col="red")
axis(1); axis(2)
}
return(list('Bayesian p-value'=bp1))
}2.1.5 Check goodness-of-fit
params_pois <-c( "sigma.time",
"mu", "beta")
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\SnowyOwl_HawkMountain\\pois.Rdata")
plot.diag(pois, lab="Poisson") # posterior predictive check
Fig. S3. Posterior predictive check for the Poisson model.
## $`Bayesian p-value`
## [1] 0params_nb <-c( "sigma.time",
"mu", "beta", "r")
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\SnowyOwl_HawkMountain\\nb.Rdata")
plot.diag(nb, lab="Negative binomial") # posterior predictive check
Fig. S4. Posterior predictive check for the negative binomial model.
## $`Bayesian p-value`
## [1] 0.43# Check for convergence
params_zip <- c( "sigma.time", "psi",
"mu", "beta")
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\SnowyOwl_HawkMountain\\zip.Rdata")
plot.diag(zip, lab="Zero-inflated Poisson") # posterior predictive check
Fig. S5. Posterior predictive check for the zero-inflated Poisson model.
## $`Bayesian p-value`
## [1] 0.972.2 Postprocess the best-fitting model (negative binomial)
Load the required libraries, data, and model output.
library("zoo")
library ("HDInterval")
library ("MCMCvis")
library ("ggplot2")
library ("reshape2")
library("gridBase")
library("grid")
library ("gridExtra")
library ("cowplot")
library("bayestestR")
library ("ggdist")
library ("viridis")
options(scipen=999)
load("data/data.Rdata")
# load your negative binomial model here
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\SnowyOwl_HawkMountain\\nb.Rdata")
# Function for transparency in base R plots
makeTransparent<-function(someColor, alpha=100){
newColor<-col2rgb(someColor)
apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
}2.2.1 Print table of parameter estimates
#***********************
#* table with parameter estimates
#***********************
library ("knitr")
nms <- c("Utqiagvik", "Bylot Island Core",
"Karupelv Valley",
"Fennoscandia", "Wrangel")
alll <- MCMCpstr(nb, params = c("mu", "beta", "r","sigma.time"),
type="chains")
allm <- do.call(rbind, alll)
tab <- data.frame(
Median = apply(allm, 1, median),
LHDI95 = apply(allm, 1, HDInterval::hdi, credMass=0.95)[1,],
UHDI95 = apply(allm, 1, HDInterval::hdi, credMass=0.95)[2,],
LHDI80 = apply(allm, 1, HDInterval::hdi, credMass=0.80)[1,],
UHDI80 = apply(allm, 1, HDInterval::hdi, credMass=0.80)[2,],
pd = apply(allm, 1, pd)
)
tab$sig95 <- ifelse(tab$LHDI95>0 | tab$UHDI95<0, TRUE, FALSE )
tab$sig80 <- ifelse(tab$LHDI80>0 | tab$UHDI80<0, TRUE, FALSE )
tab$Site <- c(rep(nms, 3), "All sites")
tab$Parameter <- c(rep("mu", 5), rep("beta", 5),
rep("r", 5), "sigma.time")
tab <- tab[, c(9, 10, 1:8)]
tab <- tab[order(tab$Site), ]
tab[,c(3:8)] <- round(tab[,c(3:8)],3)
knitr::kable(tab[,1:8], digits=2,
caption="Table S2. Coefficient estimates from the negative binomial model.",
row.names=FALSE,
col.names = c("Site", "Parameter", "Median",
"95% Lower HDI", "95% Upper HDI",
"80% Lower HDI", "80% Upper HDI",
"Prob. direction"))| Site | Parameter | Median | 95% Lower HDI | 95% Upper HDI | 80% Lower HDI | 80% Upper HDI | Prob. direction |
|---|---|---|---|---|---|---|---|
| All sites | sigma.time | 0.17 | 0.01 | 0.53 | 0.01 | 0.34 | 1.00 |
| Bylot Island Core | mu | 1.43 | 0.45 | 2.78 | 0.70 | 2.17 | 1.00 |
| Bylot Island Core | beta | 0.12 | -1.77 | 2.25 | -1.18 | 1.32 | 0.55 |
| Bylot Island Core | r | 0.15 | 0.04 | 0.28 | 0.07 | 0.22 | 1.00 |
| Fennoscandia | mu | 2.80 | 2.03 | 3.75 | 2.26 | 3.30 | 1.00 |
| Fennoscandia | beta | 0.69 | -0.48 | 1.86 | -0.06 | 1.37 | 0.88 |
| Fennoscandia | r | 0.61 | 0.23 | 1.20 | 0.30 | 0.88 | 1.00 |
| Karupelv Valley | mu | 0.57 | -0.52 | 1.85 | -0.15 | 1.28 | 0.87 |
| Karupelv Valley | beta | -0.57 | -2.18 | 0.95 | -1.52 | 0.37 | 0.78 |
| Karupelv Valley | r | 0.14 | 0.03 | 0.28 | 0.06 | 0.21 | 1.00 |
| Utqiagvik | mu | 2.51 | 1.90 | 3.13 | 2.12 | 2.88 | 1.00 |
| Utqiagvik | beta | -0.73 | -1.86 | 0.43 | -1.43 | 0.04 | 0.90 |
| Utqiagvik | r | 0.46 | 0.22 | 0.78 | 0.28 | 0.64 | 1.00 |
| Wrangel | mu | 4.10 | 3.71 | 4.56 | 3.81 | 4.36 | 1.00 |
| Wrangel | beta | -0.04 | -1.03 | 0.86 | -0.57 | 0.59 | 0.54 |
| Wrangel | r | 1.92 | 0.81 | 3.39 | 1.09 | 2.69 | 1.00 |
2.2.2 Derive a weighted population growth rate among sites
#***********************
#* Plot an overall trend
#* weighted by abundance
#***********************
# weight by modeled abundance for each year
# weight by data, mean abundance across years
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\SnowyOwl_HawkMountain\\nb.Rdata")
beta <- MCMCpstr(nb, "beta", type="chains")[[1]]
mu <- MCMCpstr(nb, "mu", type="chains")[[1]]
pred.lam <- array(NA, dim=c(length(constl$time), nrow(mu), ncol(mu)))
for (i in 1:nrow(mu)){
for (t in 1:length(constl$time)){
pred.lam[t,i,] <- mu[i,] +
beta[i,]*constl$time[t]
}}
# set NAs for unsurveyed years
wna <- which(is.na(datl$y), arr.ind=TRUE)
for (i in 1:nrow(wna)){
pred.lam[wna[i,1],wna[i,2], ] <- NA
}
wmns <- as.matrix(datl$y+1)
wmns <- ifelse(is.na(wmns), 0, wmns)
zo <- ifelse(is.na(datl$y), 0, 1)
wts2 <- array(NA, dim(zo))
for (t in 1:nrow(zo)){
mns <- wmns[t,]*zo[t,] # impute zeroes for unmonitored years
wts2[t,1:5] <- unlist(mns/sum(mns, na.rm=T)) # calculate proportions
}
lam.post <- array (NA, dim(pred.lam))
# calculate lambda (population growth rate)
for (t in 2:nrow(pred.lam)){
for (j in 1:ncol(pred.lam)){
lam.post[t,j,] <- exp(pred.lam[t,j,]) / exp(pred.lam[t-1,j,])
}}
# calculate a weighted mean (wm)
# of population growth rate
wm.post.lam <- array (NA, c(nrow(pred.lam),dim(pred.lam)[[3]]), dimnames=list(Year=1988:2020, Iter=1:4000))
for (t in 1:nrow(pred.lam)){
for (k in 1:dim(pred.lam)[[3]]){
wm.post.lam[t,k] <- weighted.mean(x=lam.post[t,,k], w=wts2[t,], na.rm=T)
}}
lam.post.long <- melt(wm.post.lam)
# calculate overall average pop growth rate
av.post <- apply(wm.post.lam, 2, mean, na.rm=T)
median(av.post, na.rm=T)## [1] 0.9839614HDInterval::hdi(av.post, credMass = 0.95, na.rm=T)## lower upper
## 0.9480951 1.0230132
## attr(,"credMass")
## [1] 0.95HDInterval::hdi(av.post, credMass = 0.80, na.rm=T)## lower upper
## 0.9598772 1.0086814
## attr(,"credMass")
## [1] 0.8pd(av.post, null=1, na.rm=T)## Probability of Direction: 0.79# calculate pop growth rates for each year
wm.df <- data.frame(
year= 1988:2020,
m= apply(wm.post.lam, 1, median, na.rm=T),
lci95= apply(wm.post.lam, 1, HDInterval::hdi)[1,],
uci95= apply(wm.post.lam, 1, HDInterval::hdi)[2,],
lci80 = apply(wm.post.lam, 1, HDInterval::hdi, credMass=.8)[1,],
uci80 = apply(wm.post.lam, 1, HDInterval::hdi, credMass=.8)[2,],
pd = apply(wm.post.lam, 1, pd, null=1)
)
knitr::kable(wm.df, digits=c(0,3,3,3,3,3,2),
row.names=FALSE,
col.names = c("Year", "Median", "95% Lower HDI", "95% Upper HDI",
"80% Lower HDI", "80% Upper HDI", "Prob. direction"),
caption="Table S3. Population growth rates for each year.")| Year | Median | 95% Lower HDI | 95% Upper HDI | 80% Lower HDI | 80% Upper HDI | Prob. direction |
|---|---|---|---|---|---|---|
| 1988 | NA | NA | NA | NA | NA | 1.00 |
| 1989 | 1.020 | 0.962 | 1.084 | 0.980 | 1.055 | 0.76 |
| 1990 | 0.965 | 0.871 | 1.060 | 0.909 | 1.023 | 0.78 |
| 1991 | 0.997 | 0.939 | 1.055 | 0.960 | 1.031 | 0.54 |
| 1992 | 0.997 | 0.939 | 1.055 | 0.960 | 1.031 | 0.54 |
| 1993 | 0.989 | 0.942 | 1.039 | 0.960 | 1.019 | 0.69 |
| 1994 | 0.993 | 0.942 | 1.044 | 0.962 | 1.025 | 0.62 |
| 1995 | 0.981 | 0.936 | 1.023 | 0.953 | 1.009 | 0.80 |
| 1996 | 0.971 | 0.918 | 1.025 | 0.939 | 1.008 | 0.86 |
| 1997 | 1.004 | 0.953 | 1.048 | 0.975 | 1.035 | 0.56 |
| 1998 | 0.967 | 0.878 | 1.045 | 0.914 | 1.014 | 0.79 |
| 1999 | 0.988 | 0.940 | 1.033 | 0.961 | 1.018 | 0.71 |
| 2000 | 0.989 | 0.947 | 1.034 | 0.965 | 1.020 | 0.68 |
| 2001 | 0.997 | 0.944 | 1.058 | 0.960 | 1.030 | 0.55 |
| 2002 | 0.993 | 0.941 | 1.045 | 0.960 | 1.024 | 0.60 |
| 2003 | 0.992 | 0.946 | 1.037 | 0.964 | 1.021 | 0.64 |
| 2004 | 1.001 | 0.946 | 1.056 | 0.968 | 1.038 | 0.52 |
| 2005 | 0.972 | 0.926 | 1.016 | 0.942 | 0.998 | 0.89 |
| 2006 | 0.978 | 0.934 | 1.020 | 0.949 | 1.005 | 0.84 |
| 2007 | 1.006 | 0.961 | 1.054 | 0.978 | 1.036 | 0.61 |
| 2008 | 0.979 | 0.938 | 1.023 | 0.949 | 1.004 | 0.83 |
| 2009 | 0.997 | 0.940 | 1.055 | 0.960 | 1.031 | 0.55 |
| 2010 | 0.998 | 0.945 | 1.051 | 0.965 | 1.031 | 0.54 |
| 2011 | 0.962 | 0.909 | 1.016 | 0.929 | 0.997 | 0.91 |
| 2012 | 0.961 | 0.911 | 1.019 | 0.928 | 0.996 | 0.92 |
| 2013 | 0.976 | 0.921 | 1.036 | 0.936 | 1.011 | 0.78 |
| 2014 | 0.979 | 0.916 | 1.054 | 0.933 | 1.021 | 0.73 |
| 2015 | 0.964 | 0.906 | 1.017 | 0.927 | 0.998 | 0.90 |
| 2016 | 0.962 | 0.904 | 1.019 | 0.922 | 0.996 | 0.90 |
| 2017 | 0.976 | 0.921 | 1.036 | 0.936 | 1.011 | 0.78 |
| 2018 | 0.962 | 0.904 | 1.019 | 0.922 | 0.996 | 0.90 |
| 2019 | 0.979 | 0.915 | 1.049 | 0.937 | 1.022 | 0.74 |
| 2020 | 0.985 | 0.917 | 1.060 | 0.939 | 1.030 | 0.67 |
c3 <- makeTransparent("gray40", alpha=10)
c4 <- makeTransparent(4, alpha=180)
rownames(wts2) <- rownames(datl$y) <- 1988:2020
colnames(wts2) <- colnames(datl$y)
names(dimnames(datl$y)) <- c("Year", "Site")
all.df <- melt(as.matrix(datl$y), value.name="Count")
all.df$weights <- melt(as.matrix(wts2), value.name="weights")$weights
colnames(all.df) <- c("Year", "Site", "Count", "Weights")
labels <- c("Utqiagvik", "Bylot Island", "Karupelv Valley",
"Fennoscandia", "Wrangel Island")
sl <- all.df[!is.na(all.df$Count),]
p1 <- ggplot() + theme_minimal() +
coord_cartesian(clip = "off") +
geom_line(data=all.df,aes(x=Year,y=Site,color=Weights), linewidth=2)+
geom_point(data=sl,aes(x=Year,y=Site),shape="I", size=4)+
xlab("Year")+ylab("Site")+scale_color_viridis_c()+
scale_y_discrete(labels=labels)+
labs(color='Weight') +
xlim(1988, 2020) + theme(legend.position = "top") +
annotate(geom = "text", x = 1990, y = 6, label = "A", size=6)
p2 <- ggplot() + theme_minimal() +
geom_line(data=lam.post.long, aes(x=Year, y=value, group=Iter),
color="gray40", linewidth=0.5, alpha=0.05 ) +
geom_hline(yintercept=1, lwd=2, color="black", linetype="dashed") +
geom_line(data=wm.df, aes(x=year, y=lci95), color="deepskyblue3", linewidth=0.5 ) +
geom_line(data=wm.df, aes(x=year, y=uci95 ), color="deepskyblue3", linewidth=0.5) +
geom_line(data=wm.df, aes(x=year, y=lci80), color="deepskyblue3", linewidth=1 ) +
geom_line(data=wm.df, aes(x=year, y=uci80), color="deepskyblue3", linewidth=1 ) +
geom_line(data=wm.df, aes(x=year, y=m), color="deepskyblue3", linewidth=2 ) +
xlab("Year") +
ylab( expression(paste("Population growth (", lambda,")")) )+
xlim(1988, 2020) +
annotate(geom = "text", x = 1990, y = 1.1, label = "B", size=6)+
coord_cartesian(xlim=c(1988, 2020), ylim=c(0.85, 1.15))
ap12 <- align_plots(p1, p2, align="v", axis="l")
p12 <- plot_grid(ap12[[1]], ap12[[2]], nrow = 2, align="v")p12
Fig. S6. Population growth rates of Snowy Owl.
2.2.3 Derive percent change for IUCN Red List Criteria
2.2.3.1 Calculate percent change with 8-year and 10.7-year generation time
#***********
#* Calculate proportion of base year
#* or percent change
#* Short generation time
#* 8-year generations
#*************
# calculate lambda from 1996
startyr <- 2020-(3*8) # 3 generations of 8 years
dnames <- list(year=1988:2020,
site=nms,
Iter=1:4000)
iucn <- array(NA, dim=dim(wm.post.lam), dimnames=dimnames(wm.post.lam))
iucn[9,] <- 1
for (t in 10:nrow(wm.post.lam)){
iucn[t,] <- iucn[(t-1),] * wm.post.lam[t,]
}
for (t in 8:1){
iucn[t,] <- iucn[(t+1),] / wm.post.lam[(t+1),]
}
# convert to percent change
iucn <- (iucn-1) * 100
iucn.post <- melt(iucn)
i.df <- data.frame(
year = 1988:2020,
m = apply(iucn, 1, median, na.rm=T),
lci95 = apply(iucn, 1, HDInterval::hdi)[1,],
uci95 = apply(iucn, 1, HDInterval::hdi)[2,],
lci80 = apply(iucn, 1, HDInterval::hdi, credMass=0.8)[1,],
uci80 = apply(iucn, 1, HDInterval::hdi, credMass=0.8)[2,],
pd = apply(iucn, 1, pd)
)
knitr::kable(i.df, digits=c(0,1,1,1,1,1,2),
row.names=FALSE,
col.names = c("Year", "Median", "95% Lower HDI", "95% Upper HDI",
"80% Lower HDI", "80% Upper HDI", "Prob. direction"),
caption="Table S4. Percent change since 1996 using 8-year generations.")| Year | Median | 95% Lower HDI | 95% Upper HDI | 80% Lower HDI | 80% Upper HDI | Prob. direction |
|---|---|---|---|---|---|---|
| 1988 | 9.1 | -24.2 | 44.3 | -12.4 | 30.0 | 0.71 |
| 1989 | 11.4 | -18.3 | 47.7 | -10.1 | 30.6 | 0.76 |
| 1990 | 7.5 | -22.4 | 36.6 | -10.6 | 25.6 | 0.71 |
| 1991 | 7.3 | -16.2 | 31.0 | -7.9 | 21.4 | 0.74 |
| 1992 | 6.9 | -11.3 | 25.4 | -4.7 | 18.3 | 0.78 |
| 1993 | 5.8 | -6.8 | 19.7 | -3.7 | 13.4 | 0.81 |
| 1994 | 4.9 | -4.4 | 14.4 | -1.3 | 10.8 | 0.86 |
| 1995 | 2.9 | -2.7 | 8.6 | -0.8 | 6.5 | 0.86 |
| 1996 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.00 |
| 1997 | 0.4 | -4.7 | 4.8 | -2.5 | 3.5 | 0.56 |
| 1998 | -3.0 | -13.5 | 8.7 | -10.1 | 4.2 | 0.69 |
| 1999 | -4.1 | -15.3 | 9.4 | -11.4 | 4.4 | 0.73 |
| 2000 | -5.1 | -18.2 | 11.6 | -14.6 | 4.5 | 0.74 |
| 2001 | -5.3 | -21.0 | 15.5 | -16.3 | 7.0 | 0.71 |
| 2002 | -6.0 | -27.3 | 17.0 | -20.0 | 7.8 | 0.69 |
| 2003 | -6.8 | -30.4 | 21.5 | -23.1 | 9.2 | 0.69 |
| 2004 | -6.5 | -32.4 | 27.8 | -24.4 | 12.7 | 0.66 |
| 2005 | -9.2 | -37.5 | 26.6 | -28.3 | 11.4 | 0.71 |
| 2006 | -11.5 | -41.2 | 27.6 | -32.7 | 10.3 | 0.73 |
| 2007 | -10.8 | -44.1 | 32.7 | -33.5 | 14.3 | 0.70 |
| 2008 | -12.8 | -46.7 | 34.4 | -37.6 | 13.1 | 0.72 |
| 2009 | -12.9 | -52.4 | 37.6 | -41.4 | 14.4 | 0.70 |
| 2010 | -13.2 | -55.9 | 42.5 | -43.8 | 16.7 | 0.69 |
| 2011 | -16.4 | -56.0 | 43.4 | -45.2 | 16.2 | 0.73 |
| 2012 | -19.7 | -59.3 | 41.1 | -48.8 | 13.3 | 0.76 |
| 2013 | -21.5 | -64.5 | 40.0 | -52.2 | 11.7 | 0.77 |
| 2014 | -23.2 | -66.0 | 42.8 | -55.6 | 11.0 | 0.77 |
| 2015 | -26.3 | -69.3 | 41.9 | -60.6 | 7.2 | 0.79 |
| 2016 | -28.8 | -71.6 | 42.3 | -63.8 | 4.9 | 0.81 |
| 2017 | -30.5 | -74.3 | 44.3 | -65.1 | 5.7 | 0.80 |
| 2018 | -33.1 | -77.7 | 42.8 | -70.0 | 1.3 | 0.82 |
| 2019 | -34.7 | -78.4 | 47.4 | -73.4 | 0.4 | 0.82 |
| 2020 | -35.6 | -81.7 | 50.1 | -74.9 | 1.1 | 0.81 |
p3 <- ggplot() + theme_minimal() +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-20, ymax=300), color="green4", fill="green4") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-30, ymax=-20), color="green3", fill="green3") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-50, ymax=-30), color="yellow", fill="yellow") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-80, ymax=-50), color="orange", fill="orange") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-100, ymax=-80), color="red", fill="red") +
geom_line(data=iucn.post, aes(x=Year, y=value, group=Iter),
color="gray40", size=0.5, alpha=0.05 ) +
geom_hline(yintercept=1, lwd=2, color="black", linetype="dashed") +
geom_line(data=i.df, aes(x=year, y=lci95), color="black", size=0.5 ) +
geom_line(data=i.df, aes(x=year, y=uci95 ), color="black", size=0.5) +
geom_line(data=i.df, aes(x=year, y=lci80), color="black", size=1 ) +
geom_line(data=i.df, aes(x=year, y=uci80), color="black", size=1 ) +
geom_line(data=i.df, aes(x=year, y=m), color="black", size=2 ) +
annotate(geom = "text", x = 1995, y = 25, label = "Least concern", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -25, label = "Near threatened", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -40, label = "Vulnerable", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -70, label = "Endangered", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -90, label = "Critically endangered", size=4, color="black") +
xlab("") +
ylab( "Percent change")+
coord_cartesian(xlim=c(1988, 2020), ylim=c(-100, 60)) +
annotate(geom = "text", x = 1995, y = 55, label = "A", size=8) +
theme(axis.text= element_text(size=12),
axis.title=element_text(size=14),
axis.ticks = element_line(color = "black"),
plot.title = element_text(size=22)) +
labs(title="8-year generation time")
iucn2020 <- iucn.post[iucn.post$Year==2020, ]
i.df2020 <- i.df[i.df$year==2020,]
p4 <- ggplot() + theme_minimal() +
geom_rect(aes(ymin=0, ymax=2, xmin=-20, xmax=300), color="green4", fill="green4") +
geom_rect(aes(ymin=0, ymax=2, xmin=-30, xmax=-20), color="green3", fill="green3") +
geom_rect(aes(ymin=0, ymax=2, xmin=-50, xmax=-30), color="yellow", fill="yellow") +
geom_rect(aes(ymin=0, ymax=2, xmin=-80, xmax=-50), color="orange", fill="orange") +
geom_rect(aes(ymin=0, ymax=2, xmin=-100, xmax=-80), color="red", fill="red") +
stat_halfeye(data=iucn2020, aes(x = value, y=0), alpha=0.65,
slab_color="gray20", slab_fill="gray40",
point_interval="median_hdi", .width = c(0.80, 0.95),
point_size=5) +
scale_size_continuous(range = c(7, 15)) +
geom_vline(xintercept=0, lwd=2, color="black", linetype="dashed") +
xlab("") + ylab("") +
annotate(geom = "text", x = 55, y = 0.25, label = "C", size=8) +
theme(axis.text= element_text(size=12),
axis.title=element_text(size=14),
axis.ticks = element_line(color = "black")) +
coord_flip(xlim=c(-100, 60), ylim=c(0, 1), clip="on")
# calculate proportion of distribution in each category and mode
cuttab <- table(cut(iucn2020$value, breaks=c(-100,-80, -50, -30, -20, 300)))
prop <- cuttab/sum(cuttab)
df.iucn <- data.frame(IUCN.criteria=rev( c('Least concern', 'Near threatened', 'Vulnerable', 'Endangered', 'Critically endangered')),
Proportion.within=prop,
Proportion.within.and.worse= cumsum(prop)
)
df.iucn <- df.iucn[nrow(df.iucn):1,]
knitr::kable(df.iucn, digits=c(0, 0, 2, 2),
row.names=FALSE,
col.names = c("IUCN Category", "A2 Criteria", "Prop within", "Prop within and worse"),
caption="Table S5. Percent change over three generations using 8-year generations.")| IUCN Category | A2 Criteria | Prop within | Prop within and worse |
|---|---|---|---|
| Least concern | (-20,300] | 0.33 | 1.00 |
| Near threatened | (-30,-20] | 0.11 | 0.67 |
| Vulnerable | (-50,-30] | 0.28 | 0.57 |
| Endangered | (-80,-50] | 0.28 | 0.29 |
| Critically endangered | (-100,-80] | 0.01 | 0.01 |
# calculate the mode
dens <- density(iucn2020$value)
# mode
dens$x[dens$y==max(dens$y)]## [1] -45.02996# median
median(iucn2020$value)## [1] -35.60979# mean
mean(iucn2020$value)## [1] -26.587# 80% HDI
HDInterval::hdi(iucn2020$value, credMass = 0.80)## lower upper
## -74.898768 1.149189
## attr(,"credMass")
## [1] 0.8#***********
#* Calculate proportion of base year
#* or percent change
#* Longer generation time
#* 10.7 yr generations
#*************
# calculate lambda from 1998
max.gen.time <- (2020-1988)/3
startyr <- 1988 # 3 generations of 8 years
dnames <- list(year=1988:2020,
site=nms,
Iter=1:4000)
iucn <- array(NA, dim=dim(wm.post.lam), dimnames=dimnames(wm.post.lam))
iucn[1,] <- 1
for (t in 2:nrow(wm.post.lam)){
iucn[t,] <- iucn[(t-1),] * wm.post.lam[t,]
}
# convert to percent change
iucn <- (iucn-1) * 100
iucn.post <- melt(iucn)
i.df <- data.frame(
year = 1988:2020,
m = apply(iucn, 1, median, na.rm=T),
lci95 = apply(iucn, 1, HDInterval::hdi)[1,],
uci95 = apply(iucn, 1, HDInterval::hdi)[2,],
lci80 = apply(iucn, 1, HDInterval::hdi, credMass=0.8)[1,],
uci80 = apply(iucn, 1, HDInterval::hdi, credMass=0.8)[2,],
pd = apply(iucn, 1, pd)
)
knitr::kable(i.df, digits=c(0,1,1,1,1,1,2),
row.names=FALSE,
col.names = c("Year", "Median", "95% Lower HDI", "95% Upper HDI",
"80% Lower HDI", "80% Upper HDI", "Prob. direction"),
caption="Table S5. Percent change since 1996 using 10.7-year generations.")| Year | Median | 95% Lower HDI | 95% Upper HDI | 80% Lower HDI | 80% Upper HDI | Prob. direction |
|---|---|---|---|---|---|---|
| 1988 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.00 |
| 1989 | 2.0 | -3.8 | 8.4 | -2.0 | 5.5 | 0.76 |
| 1990 | -1.6 | -15.5 | 11.9 | -9.6 | 6.9 | 0.60 |
| 1991 | -1.7 | -16.2 | 13.6 | -10.9 | 7.2 | 0.60 |
| 1992 | -2.0 | -18.8 | 17.1 | -12.4 | 10.1 | 0.59 |
| 1993 | -3.1 | -21.9 | 19.7 | -16.3 | 10.0 | 0.62 |
| 1994 | -3.6 | -27.6 | 22.4 | -19.0 | 12.2 | 0.62 |
| 1995 | -5.5 | -31.6 | 23.9 | -22.9 | 11.6 | 0.65 |
| 1996 | -8.3 | -34.4 | 22.9 | -26.5 | 9.4 | 0.71 |
| 1997 | -7.9 | -35.4 | 26.4 | -28.1 | 10.9 | 0.69 |
| 1998 | -10.8 | -40.2 | 29.1 | -32.6 | 10.6 | 0.73 |
| 1999 | -11.9 | -44.9 | 30.1 | -37.0 | 9.9 | 0.73 |
| 2000 | -12.8 | -48.2 | 32.2 | -37.6 | 12.6 | 0.73 |
| 2001 | -13.1 | -50.8 | 37.3 | -38.8 | 15.7 | 0.71 |
| 2002 | -13.8 | -53.3 | 42.3 | -42.2 | 16.2 | 0.71 |
| 2003 | -14.5 | -56.1 | 46.2 | -44.9 | 17.3 | 0.70 |
| 2004 | -14.5 | -62.4 | 46.5 | -49.4 | 16.6 | 0.69 |
| 2005 | -17.0 | -63.7 | 47.9 | -48.8 | 18.8 | 0.71 |
| 2006 | -18.8 | -64.7 | 50.2 | -55.1 | 14.7 | 0.72 |
| 2007 | -18.1 | -66.6 | 56.0 | -54.4 | 19.5 | 0.71 |
| 2008 | -19.9 | -70.1 | 55.9 | -57.5 | 18.2 | 0.72 |
| 2009 | -20.2 | -70.4 | 63.5 | -60.0 | 20.2 | 0.71 |
| 2010 | -20.5 | -72.8 | 69.5 | -62.0 | 22.6 | 0.70 |
| 2011 | -23.5 | -75.1 | 66.2 | -67.5 | 16.4 | 0.72 |
| 2012 | -26.4 | -76.8 | 63.6 | -67.6 | 15.8 | 0.75 |
| 2013 | -28.1 | -77.9 | 66.9 | -70.8 | 13.0 | 0.75 |
| 2014 | -29.6 | -80.0 | 67.6 | -71.6 | 13.2 | 0.76 |
| 2015 | -32.0 | -82.4 | 65.4 | -77.0 | 7.5 | 0.77 |
| 2016 | -34.6 | -82.6 | 66.5 | -78.0 | 6.0 | 0.79 |
| 2017 | -36.2 | -84.4 | 68.6 | -80.1 | 4.7 | 0.79 |
| 2018 | -38.9 | -85.8 | 66.8 | -81.7 | 2.7 | 0.80 |
| 2019 | -40.2 | -87.6 | 68.7 | -82.1 | 3.8 | 0.80 |
| 2020 | -41.0 | -91.9 | 69.7 | -84.2 | 3.5 | 0.80 |
p5 <- ggplot() + theme_minimal() +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-20, ymax=300), color="green4", fill="green4") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-30, ymax=-20), color="green3", fill="green3") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-50, ymax=-30), color="yellow", fill="yellow") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-80, ymax=-50), color="orange", fill="orange") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-100, ymax=-80), color="red", fill="red") +
geom_line(data=iucn.post, aes(x=Year, y=value, group=Iter),
color="gray40", size=0.5, alpha=0.05 ) +
geom_hline(yintercept=1, lwd=2, color="black", linetype="dashed") +
geom_line(data=i.df, aes(x=year, y=lci95), color="black", size=0.5 ) +
geom_line(data=i.df, aes(x=year, y=uci95 ), color="black", size=0.5) +
geom_line(data=i.df, aes(x=year, y=lci80), color="black", size=1 ) +
geom_line(data=i.df, aes(x=year, y=uci80), color="black", size=1 ) +
geom_line(data=i.df, aes(x=year, y=m), color="black", size=2 ) +
annotate(geom = "text", x = 1995, y = 25, label = "Least concern", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -25, label = "Near threatened", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -40, label = "Vulnerable", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -70, label = "Endangered", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -90, label = "Critically endangered", size=4, color="black") +
xlab("Year") +
ylab( "Percent change")+
coord_cartesian(xlim=c(1988, 2020), ylim=c(-100, 60)) +
annotate(geom = "text", x = 1995, y = 55, label = "B", size=8) +
theme(axis.text= element_text(size=12),
axis.title=element_text(size=14),
axis.ticks = element_line(color = "black"),
plot.title = element_text(size=22)) +
labs(title="10.7-year generation time")
iucn2020 <- iucn.post[iucn.post$Year==2020, ]
i.df2020 <- i.df[i.df$year==2020,]
p6 <- ggplot() + theme_minimal() +
geom_rect(aes(ymin=0, ymax=2, xmin=-20, xmax=300), color="green4", fill="green4") +
geom_rect(aes(ymin=0, ymax=2, xmin=-30, xmax=-20), color="green3", fill="green3") +
geom_rect(aes(ymin=0, ymax=2, xmin=-50, xmax=-30), color="yellow", fill="yellow") +
geom_rect(aes(ymin=0, ymax=2, xmin=-80, xmax=-50), color="orange", fill="orange") +
geom_rect(aes(ymin=0, ymax=2, xmin=-100, xmax=-80), color="red", fill="red") +
stat_halfeye(data=iucn2020, aes(x = value, y=0), alpha=0.65,
slab_color="gray20", slab_fill="gray40",
point_interval="median_hdi", .width = c(0.80, 0.95),
point_size=5) +
scale_size_continuous(range = c(7, 15)) +
geom_vline(xintercept=0, lwd=2, color="black", linetype="dashed") +
xlab("") + ylab("Density (scaled)\nof percent change\nover three generations") +
annotate(geom = "text", x = 55, y = 0.25, label = "D", size=8) +
theme(axis.text= element_text(size=12),
axis.title=element_text(size=14),
axis.ticks = element_line(color = "black")) +
coord_flip(xlim=c(-100, 60), ylim=c(0, 1), clip="on")
ap56 <- align_plots(p5, p6, align="h", axis="l")
p56 <- plot_grid(ap56[[1]], ap56[[2]], nrow = 1, align="h", rel_widths = c(2, 1))
ap56 <- align_plots(p3, p4, p5, p6, align="h", axis="l")
p56 <- plot_grid(ap56[[1]], ap56[[2]],
ap56[[3]], ap56[[4]],
nrow = 2, align="h", rel_widths = c(2, 1)) 2.2.3.2 Plot percent change
p56
Fig. S7. Percent change in abundance of Snowy Owl.
# calculate proportion of distribution in each category and mode
cuttab <- table(cut(iucn2020$value, breaks=c(-100,-80, -50, -30, -20, 300)))
prop <- cuttab/sum(cuttab)
df.iucn <- data.frame(IUCN.criteria=rev( c('Least concern', 'Near threatened', 'Vulnerable', 'Endangered', 'Critically endangered')),
Proportion.within=prop,
Proportion.within.and.worse= cumsum(prop)
)
df.iucn <- df.iucn[nrow(df.iucn):1,]
knitr::kable(df.iucn, digits=c(0, 0, 2, 2),
row.names=FALSE,
col.names = c("IUCN Category", "A2 Criteria", "Prop within", "Prop within and worse"),
caption="Table S6. Percent change over three generations using 10.7-year generations.")| IUCN Category | A2 Criteria | Prop within | Prop within and worse |
|---|---|---|---|
| Least concern | (-20,300] | 0.31 | 1.00 |
| Near threatened | (-30,-20] | 0.08 | 0.69 |
| Vulnerable | (-50,-30] | 0.23 | 0.61 |
| Endangered | (-80,-50] | 0.34 | 0.38 |
| Critically endangered | (-100,-80] | 0.04 | 0.04 |
# calculate the mode
dens <- density(iucn2020$value)
# mode
dens$x[dens$y==max(dens$y)]## [1] -53.66316# median
median(iucn2020$value)## [1] -40.95308# mean
mean(iucn2020$value)## [1] -27.24941# 80% HDI
HDInterval::hdi(iucn2020$value, credMass = 0.80)## lower upper
## -84.193152 3.495978
## attr(,"credMass")
## [1] 0.82.3 Replicate analysis but include sites with shorter time series and site as a random effect
We tested whether a model including a random factor for site (N=9). This analysis includes sites with less abundance and shorter time series that were found to be less reliable in simulations.
2.3.1 Load packages and manipulate data
library(readxl)
library(reshape2)
library(MCMCvis)
library(HDInterval)
library(ggplot2)
library(tidybayes)
library(bayestestR)
library(cowplot)
library(parallel)
dat <- read.csv("data\\data.csv", na.strings="NA")
# change Bylot Island into
# 2 separate sites
# one of 100 km^2 and
# one of 300 km^2
dat$BylotIsland_300 <- dat$BylotIsland_400-dat$BylotIsland_100
colSums(!is.na(dat))## year Utqiagvik_213 BylotIsland_100
## 35 29 28
## BylotIsland_400 Igloolik.Island_114 Karupelv.Valley_75
## 20 10 33
## Hochstetter.Forland_100 Fennoscandia_x Alaska
## 11 13 32
## Wrangel BylotIsland_300
## 20 20# Bundle data
y <- as.data.frame(dat[-c(1,2), -c(4)])
# create data lists for nimble
y.nim <- y[,-1]
datl <- list(y = y.nim) # exclude 8-fennoscania because no area information
constl <- list(nsite = ncol(y.nim),
ntime = nrow(y.nim),
time = (y$year-2004)/16
)2.3.2 Run the model
This model includes data from all 9 sites and has a random effect for both site and time.
#***********************
#* Negative binomial model
#* Best fitting model
#* having site and time random effects
#* in Appendix 3
#***********************
run <- function(seed, datl, constl){
library('nimble')
library('coda')
library ('nimbleHMC')
code <- nimbleCode(
{
# priors
mu ~ dnorm(0, sd=5) # constrain to reasonable values <exp(5) or <148
beta ~ dnorm(0, sd=5)
for (j in 1:nsite){
r[j] ~ dexp(0.2)
}
sigma.time ~ dexp(2)
sigma.site ~ dexp(2)
# likelihood
for (t in 1:ntime){
for (j in 1:nsite){
p[t,j] <- r[j]/(r[j]+lambda[t,j])
y[t,j] ~ dnegbin(p[t,j], r[j])
# abundance
log(lambda[t,j]) <- mu +
beta*time[t] +
eps.time[t] + eps.site[j]
} # j
eps.time[t] ~ dnorm(0, sd=sigma.time)
} # t
for (j in 1:nsite){ eps.site[j] ~ dnorm(0, sigma.site) }
###################
# Assess GOF of the 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
###################
# GOF for model
for (t in 1:ntime){
for (j in 1:nsite){
c.obs[t,j] <- y[t,j] # observed counts
c.exp[t,j] <- lambda[t,j] # expected counts adult breeder
c.rep[t,j] ~ dnegbin(p[t,j],r[j]) # expected counts
# Compute fit statistics, Mean absolute error
dssm.obs[t,j] <- abs( ( (c.obs[t,j]) - (c.exp[t,j]) ) / (c.obs[t,j]+0.001) )
dssm.rep[t,j] <- abs( ( (c.rep[t,j]) - (c.exp[t,j]) ) / (c.rep[t,j]+0.001) )
} } # j # t
dmape.obs <- sum(dssm.obs[1:ntime,1:nsite])
dmape.rep <- sum(dssm.rep[1:ntime,1:nsite])
# variance-mean ratio
tvm.rep <- sd(c.rep[1:ntime,1:nsite])^2/mean(c.rep[1:ntime,1:nsite])
tvm.obs <- sd(y[1:ntime,1:nsite])^2/mean(y[1:ntime,1:nsite])
}
) # end model
params <- c( "sigma.time", "sigma.site",
"mu", "beta", "r", "p",
"dssm.obs", "dssm.rep",
"dmape.obs", "dmape.rep",
"tvm.rep", "tvm.obs"
)
# create initial values for missing y data
y.inits <- array(NA, dim(datl$y))
y.inits[is.na(datl$y)] <- rpois(sum(is.na(datl$y)), 2)
inits <- function(){ list (y = y.inits,
sigma.time = rexp(1),
sigma.site = rexp(1),
mu = runif(1, -2, 2),
beta = runif(1, -2, 2),
r = runif(constl$nsite)
)}
n.chains=1; n.thin=10; n.iter=20000; n.burnin=10000
mod <- nimbleModel(code, calculate=T, constants = constl,
data = datl, inits = inits(),
buildDerivs = TRUE)
mod$calculate()
cmod <- compileNimble(mod )
confhmc <- configureHMC(mod)
confhmc$setMonitors(params)
hmc <- buildMCMC(confhmc)
chmc <- compileNimble(hmc, project=mod, resetFunctions = TRUE)
post <- runMCMC(chmc,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
thin = n.thin,
samplesAsCodaMCMC = T)
return(post)
}
# approximately 15 minute runtime
this_cluster <- makeCluster(4)
post <- parLapply(cl = this_cluster,
X = 1:4,
fun = run,
dat=datl,
const=constl)
stopCluster(this_cluster)
nb <- list(as.mcmc(post[[1]]),
as.mcmc(post[[2]]),
as.mcmc(post[[3]]),
as.mcmc(post[[4]]))2.3.3 Postprocess and plot population growth
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\SnowyOwl_HawkMountain\\nb_site-re_allsites.Rdata")
# Check for convergence
params_nb <-c( "sigma.time", "sigma.site",
"mu", "beta", "r")
MCMCtrace(nb, params_nb, pdf=F,
ind = TRUE, Rhat = TRUE, n.eff = TRUE)
par(mfrow=c(1,1))
MCMCplot(object = nb, params = params_nb)
coeftab95 <- MCMCsummary(nb, params_nb,
HPD=TRUE, hpd_prob = 0.95,
round=2, pg0 = TRUE, func = median,
func_name = "Median")
coeftab80 <- MCMCsummary(nb, params_nb,
HPD=TRUE, hpd_prob = 0.80,
round=2, pg0 = TRUE, func = median,
func_name = "Median")
mu <- MCMCpstr(nb, "mu", type="chains")[[1]]
beta <- MCMCpstr(nb, "beta", type="chains")[[1]]
pred.lam <- array(NA, dim=c(length(constl$time), ncol(mu)))
for (t in 1:length(constl$time)){
pred.lam[t,] <- mu +
beta[1,]*constl$time[t]
}
lam.post <- array (NA, dim(pred.lam), dimnames=list(Year=1988:2020, Iter=1:4000))
# calculate lambda (population growth rate)
for (t in 2:nrow(pred.lam)){
lam.post[t,] <- exp(pred.lam[t,]) / exp(pred.lam[t-1,])
}
lam.post.long <- melt(lam.post)
# calculate overall average pop growth rate
av.post <- apply(lam.post, 2, mean, na.rm=T)
median(av.post, na.rm=T)## [1] 0.9883551HDInterval::hdi(av.post, credMass = 0.95, na.rm=T)## lower upper
## 0.9567679 1.0178284
## attr(,"credMass")
## [1] 0.95HDInterval::hdi(av.post, credMass = 0.80, na.rm=T)## lower upper
## 0.9695829 1.0084350
## attr(,"credMass")
## [1] 0.8pd(av.post, null=1, na.rm=T)## Probability of Direction: 0.79# calculate pop growth rates for each year
lam.df <- data.frame(
year= 1988:2020,
m= apply(lam.post, 1, median, na.rm=T),
lci= apply(lam.post, 1, HDInterval::hdi)[1,],
uci= apply(lam.post, 1, HDInterval::hdi)[2,],
lci80 = apply(lam.post, 1, HDInterval::hdi, credMass=.8)[1,],
uci80 = apply(lam.post, 1, HDInterval::hdi, credMass=.8)[2,],
pd = apply(lam.post, 1, pd, null=1)
)
makeTransparent<-function(someColor, alpha=100){
newColor<-col2rgb(someColor)
apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
}
c3 <- makeTransparent("gray40", alpha=10)
c4 <- makeTransparent(4, alpha=180)
rownames(datl$y) <- 1988:2020
names(dimnames(datl$y)) <- c("Year", "Site")
all.df <- melt(as.matrix(datl$y), value.name="Count")
colnames(all.df) <- c("Year", "Site", "Count")
labels <- unique(all.df$Site)
sl <- all.df[!is.na(all.df$Count),]
p1 <- ggplot() + theme_minimal() +
coord_cartesian(clip = "off") +
geom_line(data=all.df,aes(x=Year,y=Site),linewidth=2, color='gray60')+
geom_point(data=sl,aes(x=Year,y=Site),shape="I",size=4)+
xlab("Year")+ylab("Site")+scale_color_viridis_c()+
scale_y_discrete(labels=labels)+
labs(color='black') +
xlim(1988, 2020) + #theme(legend.position = "top") +
annotate(geom = "text", x = 1990, y = 9.5, label = "A", size=6)
p2 <- ggplot() + theme_minimal() +
geom_line(data=lam.post.long, aes(x=Year, y=value, group=Iter),
color="gray40", linewidth=0.5, alpha=0.05 ) +
geom_hline(yintercept=1, lwd=2, color="black", linetype="dashed") +
geom_line(data=lam.df, aes(x=year, y=lci), color="deepskyblue3", linewidth=0.5 ) +
geom_line(data=lam.df, aes(x=year, y=uci ), color="deepskyblue3", linewidth=0.5) +
geom_line(data=lam.df, aes(x=year, y=lci80), color="deepskyblue3", linewidth=1 ) +
geom_line(data=lam.df, aes(x=year, y=uci80), color="deepskyblue3", linewidth=1 ) +
geom_line(data=lam.df, aes(x=year, y=m), color="deepskyblue3", linewidth=2 ) +
xlab("Year") +
ylab( expression(paste("Population growth (", lambda,")")) )+
xlim(1988, 2020) +
annotate(geom = "text", x = 1990, y = 1.1, label = "B", size=6)+
coord_cartesian(xlim=c(1988, 2020), ylim=c(0.85, 1.15))
ap12 <- align_plots(p1, p2, align="v", axis="l")
p12 <- plot_grid(ap12[[1]], ap12[[2]], nrow = 2, align="v")p12
Fig. S8. Population growth rates of Snowy Owl using all data and a random effect for site.
2.3.4 Postprocess and plot IUCN Criteria as percent change
# plot percent change
perc.post <- array(NA, dim=dim(pred.lam), dimnames=list(Year=1988:2020, Iter=1:4000))
# calc lambda for each site relative to 1996
for (t in 1:nrow(pred.lam)){
perc.post[t,] <- exp(pred.lam[t,]) / exp(pred.lam[9,])
}
perc.post.long <- melt(perc.post)
iucn <- (perc.post-1) * 100
iucn.post <- melt(iucn)
i.df <- data.frame(
year = 1988:2020,
m = apply(iucn, 1, median, na.rm=T),
lci95 = apply(iucn, 1, HDInterval::hdi)[1,],
uci95 = apply(iucn, 1, HDInterval::hdi)[2,],
lci80 = apply(iucn, 1, HDInterval::hdi, credMass=0.8)[1,],
uci80 = apply(iucn, 1, HDInterval::hdi, credMass=0.8)[2,],
pd = apply(iucn, 1, pd)
)
knitr::kable(i.df, digits=c(0,1,1,1,1,1,2),
row.names=FALSE,
col.names = c("Year", "Median", "95% Lower HDI", "95% Upper HDI",
"80% Lower HDI", "80% Upper HDI", "Prob. direction"),
caption="Table S7. Percent change since 1996 using 8-year generations.")| Year | Median | 95% Lower HDI | 95% Upper HDI | 80% Lower HDI | 80% Upper HDI | Prob. direction |
|---|---|---|---|---|---|---|
| 1988 | 9.8 | -15.3 | 39.3 | -7.7 | 26.8 | 0.79 |
| 1989 | 8.5 | -13.5 | 33.7 | -5.7 | 24.1 | 0.79 |
| 1990 | 7.3 | -11.7 | 28.2 | -4.9 | 20.4 | 0.79 |
| 1991 | 6.0 | -9.9 | 23.0 | -4.1 | 16.7 | 0.79 |
| 1992 | 4.8 | -8.0 | 18.0 | -3.3 | 13.2 | 0.79 |
| 1993 | 3.6 | -5.9 | 13.4 | -2.5 | 9.7 | 0.79 |
| 1994 | 2.4 | -3.9 | 8.8 | -1.7 | 6.4 | 0.79 |
| 1995 | 1.2 | -1.8 | 4.5 | -0.8 | 3.1 | 0.79 |
| 1996 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.00 |
| 1997 | -1.2 | -4.3 | 1.8 | -3.0 | 0.8 | 0.79 |
| 1998 | -2.3 | -8.5 | 3.6 | -6.0 | 1.7 | 0.79 |
| 1999 | -3.5 | -12.4 | 5.4 | -9.1 | 2.3 | 0.79 |
| 2000 | -4.6 | -16.2 | 7.3 | -12.0 | 3.0 | 0.79 |
| 2001 | -5.7 | -19.8 | 9.2 | -14.7 | 3.8 | 0.79 |
| 2002 | -6.8 | -23.3 | 11.2 | -18.1 | 3.8 | 0.79 |
| 2003 | -7.9 | -26.6 | 13.2 | -20.8 | 4.5 | 0.79 |
| 2004 | -8.9 | -30.6 | 14.4 | -23.4 | 5.1 | 0.79 |
| 2005 | -10.0 | -33.7 | 16.3 | -26.9 | 4.8 | 0.79 |
| 2006 | -11.1 | -36.6 | 18.3 | -29.4 | 5.4 | 0.79 |
| 2007 | -12.1 | -39.4 | 20.3 | -31.8 | 5.9 | 0.79 |
| 2008 | -13.1 | -42.1 | 22.3 | -34.1 | 6.5 | 0.79 |
| 2009 | -14.1 | -44.7 | 24.4 | -36.6 | 6.8 | 0.79 |
| 2010 | -15.1 | -47.2 | 26.5 | -40.5 | 5.7 | 0.79 |
| 2011 | -16.1 | -49.5 | 28.7 | -42.6 | 6.1 | 0.79 |
| 2012 | -17.1 | -51.8 | 30.8 | -44.7 | 6.5 | 0.79 |
| 2013 | -18.1 | -57.8 | 29.1 | -46.7 | 6.9 | 0.79 |
| 2014 | -19.0 | -59.9 | 31.1 | -48.7 | 7.3 | 0.79 |
| 2015 | -20.0 | -61.8 | 33.0 | -50.5 | 7.8 | 0.79 |
| 2016 | -20.9 | -63.7 | 35.1 | -53.3 | 7.2 | 0.79 |
| 2017 | -21.8 | -65.5 | 37.1 | -55.1 | 7.6 | 0.79 |
| 2018 | -22.7 | -67.2 | 39.2 | -56.7 | 7.9 | 0.79 |
| 2019 | -23.6 | -68.8 | 41.3 | -58.4 | 8.3 | 0.79 |
| 2020 | -24.5 | -70.4 | 43.4 | -59.9 | 8.7 | 0.79 |
p3 <- ggplot() + theme_minimal() +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-20, ymax=300), color="green4", fill="green4") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-30, ymax=-20), color="green3", fill="green3") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-50, ymax=-30), color="yellow", fill="yellow") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-80, ymax=-50), color="orange", fill="orange") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-100, ymax=-80), color="red", fill="red") +
geom_line(data=iucn.post, aes(x=Year, y=value, group=Iter),
color="gray40", size=0.5, alpha=0.05 ) +
geom_hline(yintercept=1, lwd=2, color="black", linetype="dashed") +
geom_line(data=i.df, aes(x=year, y=lci95), color="black", size=0.5 ) +
geom_line(data=i.df, aes(x=year, y=uci95 ), color="black", size=0.5) +
geom_line(data=i.df, aes(x=year, y=lci80), color="black", size=1 ) +
geom_line(data=i.df, aes(x=year, y=uci80), color="black", size=1 ) +
geom_line(data=i.df, aes(x=year, y=m), color="black", size=2 ) +
annotate(geom = "text", x = 1995, y = 25, label = "Least concern", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -25, label = "Near threatened", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -40, label = "Vulnerable", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -70, label = "Endangered", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -90, label = "Critically endangered", size=4, color="black") +
xlab("Year") +
ylab( "Percent change")+
coord_cartesian(xlim=c(1988, 2020), ylim=c(-100, 60)) +
annotate(geom = "text", x = 1995, y = 55, label = "A", size=8) +
theme(axis.text= element_text(size=12),
axis.title=element_text(size=14),
axis.ticks = element_line(color = "black"))
iucn2020 <- iucn.post[iucn.post$Year==2020, ]
i.df2020 <- i.df[i.df$year==2020,]
p4 <- ggplot() + theme_minimal() +
geom_rect(aes(ymin=0, ymax=2, xmin=-20, xmax=300), color="green4", fill="green4") +
geom_rect(aes(ymin=0, ymax=2, xmin=-30, xmax=-20), color="green3", fill="green3") +
geom_rect(aes(ymin=0, ymax=2, xmin=-50, xmax=-30), color="yellow", fill="yellow") +
geom_rect(aes(ymin=0, ymax=2, xmin=-80, xmax=-50), color="orange", fill="orange") +
geom_rect(aes(ymin=0, ymax=2, xmin=-100, xmax=-80), color="red", fill="red") +
stat_halfeye(data=iucn2020, aes(x = value, y=0), alpha=0.65,
slab_color="gray20", slab_fill="gray40",
point_interval="median_hdi", .width = c(0.80, 0.95),
point_size=5) +
scale_size_continuous(range = c(7, 15)) +
geom_vline(xintercept=0, lwd=2, color="black", linetype="dashed") +
xlab("") + ylab("Density (scaled)\nof percent change\nover three generations") +
annotate(geom = "text", x = 55, y = 0.25, label = "B", size=8) +
theme(axis.text= element_text(size=12),
axis.title=element_text(size=14),
axis.ticks = element_line(color = "black")) +
coord_flip(xlim=c(-100, 60), ylim=c(0, 1), clip="on")
ap45 <- align_plots(p3, p4, align="h", axis="l")
p45 <- plot_grid(ap45[[1]], ap45[[2]], nrow = 1, align="h", rel_widths = c(2, 1))p45
Fig. S9. Percent change in abundance of Snowy Owl since 1996 using all data and a random effect for site.
# calculate the mode
dens <- density(iucn2020$value)
# mode
dens$x[dens$y==max(dens$y)]## [1] -35.76325# median
median(iucn2020$value)## [1] -24.50609# mean
mean(iucn2020$value)## [1] -19.36453# 80% HDI
HDInterval::hdi(iucn2020$value, credMass = 0.80)## lower upper
## -59.907153 8.675562
## attr(,"credMass")
## [1] 0.82.4 Replicate analysis but exclude the only site with significant declines (Utqiagvik)
Load the required libraries, data, and model output.
library("zoo")
library ("HDInterval")
library ("MCMCvis")
library ("ggplot2")
library ("reshape2")
library("gridBase")
library("grid")
library ("gridExtra")
library ("cowplot")
library("bayestestR")
library ("ggdist")
library ("viridis")
library ('nimble')
library ('nimbleHMC')
library ('MCMCvis')
library ('coda')
library('parallel')
set.seed(5757575)
options(scipen=999)
load("data/data.Rdata")
# remove Utqiagvik
datl$y <- datl$y[, -1]
# update number of sites
constl$nsite <- ncol(datl$y)Use the same model as before but without Utqiagvik data
#***********************
#* Negative binomial model
#* Best fitting model
#***********************
run <- function(seed, datl, constl){
library('nimble')
library('coda')
library ('nimbleHMC')
code <- nimbleCode(
{
# priors
# priors
for (j in 1:nsite){
mu[j] ~ dnorm(0, sd=5) # constrain to reasonable values <exp(5) or <148
beta[j] ~ dnorm(0, sd=5)
r[j] ~ dexp(0.1)
} # j
sigma.time ~ dexp(2)
# likelihood
for (t in 1:ntime){
for (j in 1:nsite){
p[t,j] <- r[j]/(r[j]+lambda[t,j])
y[t,j] ~ dnegbin(p[t,j], r[j])
# abundance
log(lambda[t,j]) <- mu[j] +
beta[j]*time[t] +
eps.time[t]
} # j
eps.time[t] ~ dnorm(0, sd=sigma.time)
} # t
###################
# Assess GOF of the 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
###################
# GOF for model
for (t in 1:ntime){
for (j in 1:nsite){
c.obs[t,j] <- y[t,j] # observed counts
c.exp[t,j] <- lambda[t,j] # expected counts adult breeder
c.rep[t,j] ~ dnegbin(p[t,j],r[j]) # expected counts
# Compute fit statistics, Mean absolute error
dssm.obs[t,j] <- abs( ( (c.obs[t,j]) - (c.exp[t,j]) ) / (c.obs[t,j]+0.001) )
dssm.rep[t,j] <- abs( ( (c.rep[t,j]) - (c.exp[t,j]) ) / (c.rep[t,j]+0.001) )
} } # j # t
dmape.obs <- sum(dssm.obs[1:ntime,1:nsite])
dmape.rep <- sum(dssm.rep[1:ntime,1:nsite])
# variance-mean ratio
tvm.rep <- sd(c.rep[1:ntime,1:nsite])^2/mean(c.rep[1:ntime,1:nsite])
tvm.obs <- sd(y[1:ntime,1:nsite])^2/mean(y[1:ntime,1:nsite])
}
) # end model
params <- c( "sigma.time",
"mu", "beta", "r", "p",
"dssm.obs", "dssm.rep",
"dmape.obs", "dmape.rep",
"tvm.rep", "tvm.obs"
)
# create initial values for missing y data
y.inits <- array(NA, dim(datl$y))
y.inits[is.na(datl$y)] <- rpois(sum(is.na(datl$y)), 2)
inits <- function(){ list (y = y.inits,
sigma.time = rexp(1),
mu = runif(constl$nsite, -2, 2),
beta = runif(constl$nsite, -2, 2),
r = runif(constl$nsite)
)}
n.chains=1; n.thin=10; n.iter=20000; n.burnin=10000
mod <- nimbleModel(code, calculate=T, constants = constl,
data = datl, inits = inits(),
buildDerivs = TRUE)
mod$calculate()
cmod <- compileNimble(mod )
confhmc <- configureHMC(mod)
confhmc$setMonitors(params)
hmc <- buildMCMC(confhmc)
chmc <- compileNimble(hmc, project=mod, resetFunctions = TRUE)
post <- runMCMC(chmc,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
thin = n.thin,
samplesAsCodaMCMC = T)
return(post)
}
# approximately 15 minute runtime
this_cluster <- makeCluster(4)
post <- parLapply(cl = this_cluster,
X = 1:4,
fun = run,
dat=datl,
const=constl)
stopCluster(this_cluster)
nb <- list(as.mcmc(post[[1]]),
as.mcmc(post[[2]]),
as.mcmc(post[[3]]),
as.mcmc(post[[4]]))
params_nb <-c( "sigma.time",
"mu", "beta", "r")
# Check for convergence
MCMCtrace(nb, params_nb, pdf=F,
ind = TRUE, Rhat = TRUE, n.eff = TRUE)
par(mfrow=c(1,1))
MCMCplot(object = nb, params = params_nb)
# save(out=nb, post=post, run,
# file="C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\SnowyOwl_HawkMountain\\nb_noUtqiagvik.Rdata")2.4.1 Print table of parameter estimates
#***********************
#* table with parameter estimates
#***********************
# load your negative binomial model here
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\SnowyOwl_HawkMountain\\nb_noUtqiagvik.Rdata")
# Function for transparency in base R plots
makeTransparent<-function(someColor, alpha=100){
newColor<-col2rgb(someColor)
apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
}
library ("knitr")
nms <- c("Bylot Island Core",
"Karupelv Valley",
"Fennoscandia", "Wrangel")
alll <- MCMCpstr(nb, params = c("mu", "beta", "r","sigma.time"),
type="chains")
allm <- do.call(rbind, alll)
tab <- data.frame(
Median = apply(allm, 1, median),
LHDI95 = apply(allm, 1, HDInterval::hdi, credMass=0.95)[1,],
UHDI95 = apply(allm, 1, HDInterval::hdi, credMass=0.95)[2,],
LHDI80 = apply(allm, 1, HDInterval::hdi, credMass=0.80)[1,],
UHDI80 = apply(allm, 1, HDInterval::hdi, credMass=0.80)[2,],
pd = apply(allm, 1, pd)
)
tab$sig95 <- ifelse(tab$LHDI95>0 | tab$UHDI95<0, TRUE, FALSE )
tab$sig80 <- ifelse(tab$LHDI80>0 | tab$UHDI80<0, TRUE, FALSE )
tab$Site <- c(rep(nms, 3), "All sites")
tab$Parameter <- c(rep("mu", 4), rep("beta", 4),
rep("r", 4), "sigma.time")
tab <- tab[, c(9, 10, 1:8)]
tab <- tab[order(tab$Site), ]
tab[,c(3:8)] <- round(tab[,c(3:8)],3)
knitr::kable(tab[,1:8], digits=2,
caption="Table S8. Coefficient estimates from the negative binomial model.",
row.names=FALSE,
col.names = c("Site", "Parameter", "Median",
"95% Lower HDI", "95% Upper HDI",
"80% Lower HDI", "80% Upper HDI",
"Prob. direction"))| Site | Parameter | Median | 95% Lower HDI | 95% Upper HDI | 80% Lower HDI | 80% Upper HDI | Prob. direction |
|---|---|---|---|---|---|---|---|
| All sites | sigma.time | 0.32 | 0.02 | 1.09 | 0.02 | 0.73 | 1.00 |
| Bylot Island Core | mu | 1.56 | 0.41 | 3.03 | 0.73 | 2.37 | 1.00 |
| Bylot Island Core | beta | 0.03 | -2.26 | 2.23 | -1.29 | 1.48 | 0.51 |
| Bylot Island Core | r | 0.14 | 0.05 | 0.28 | 0.07 | 0.21 | 1.00 |
| Fennoscandia | mu | 2.67 | 1.80 | 3.66 | 2.09 | 3.23 | 1.00 |
| Fennoscandia | beta | 0.63 | -0.61 | 1.84 | -0.08 | 1.47 | 0.85 |
| Fennoscandia | r | 0.71 | 0.20 | 2.15 | 0.28 | 1.11 | 1.00 |
| Karupelv Valley | mu | 0.51 | -0.55 | 1.86 | -0.21 | 1.30 | 0.83 |
| Karupelv Valley | beta | -0.64 | -2.28 | 1.07 | -1.67 | 0.33 | 0.79 |
| Karupelv Valley | r | 0.14 | 0.04 | 0.32 | 0.06 | 0.23 | 1.00 |
| Wrangel | mu | 4.10 | 3.64 | 4.60 | 3.81 | 4.39 | 1.00 |
| Wrangel | beta | 0.01 | -0.96 | 1.06 | -0.59 | 0.66 | 0.51 |
| Wrangel | r | 2.24 | 0.67 | 12.42 | 0.94 | 3.92 | 1.00 |
2.4.2 Derive a weighted population growth rate among sites
#***********************
#* Plot an overall trend
#* weighted by abundance
#***********************
# weight by modeled abundance for each year
# weight by data, mean abundance across years
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\SnowyOwl_HawkMountain\\nb_noUtqiagvik.Rdata")
beta <- MCMCpstr(nb, "beta", type="chains")[[1]]
mu <- MCMCpstr(nb, "mu", type="chains")[[1]]
pred.lam <- array(NA, dim=c(length(constl$time), nrow(mu), ncol(mu)))
for (i in 1:nrow(mu)){
for (t in 1:length(constl$time)){
pred.lam[t,i,] <- mu[i,] +
beta[i,]*constl$time[t]
}}
# set NAs for unsurveyed years
wna <- which(is.na(datl$y), arr.ind=TRUE)
for (i in 1:nrow(wna)){
pred.lam[wna[i,1],wna[i,2], ] <- NA
}
wmns <- as.matrix(datl$y+1)
wmns <- ifelse(is.na(wmns), 0, wmns)
zo <- ifelse(is.na(datl$y), 0, 1)
wts2 <- array(NA, dim(zo))
for (t in 1:nrow(zo)){
mns <- wmns[t,]*zo[t,] # impute zeroes for unmonitored years
wts2[t,1:4] <- unlist(mns/sum(mns, na.rm=T)) # calculate proportions
}
lam.post <- array (NA, dim(pred.lam))
# calculate lambda (population growth rate)
for (t in 2:nrow(pred.lam)){
for (j in 1:ncol(pred.lam)){
lam.post[t,j,] <- exp(pred.lam[t,j,]) / exp(pred.lam[t-1,j,])
}}
# calculate a weighted mean (wm)
# of population growth rate
wm.post.lam <- array (NA, c(nrow(pred.lam),dim(pred.lam)[[3]]), dimnames=list(Year=1988:2020, Iter=1:4000))
for (t in 1:nrow(pred.lam)){
for (k in 1:dim(pred.lam)[[3]]){
wm.post.lam[t,k] <- weighted.mean(x=lam.post[t,,k], w=wts2[t,], na.rm=T)
}}
lam.post.long <- melt(wm.post.lam)
# calculate overall average pop growth rate
av.post <- apply(wm.post.lam, 2, mean, na.rm=T)
median(av.post, na.rm=T)## [1] 0.9940119HDInterval::hdi(av.post, credMass = 0.95, na.rm=T)## lower upper
## 0.9404984 1.0441790
## attr(,"credMass")
## [1] 0.95HDInterval::hdi(av.post, credMass = 0.80, na.rm=T)## lower upper
## 0.9591881 1.0245977
## attr(,"credMass")
## [1] 0.8pd(av.post, null=1, na.rm=T)## Probability of Direction: 0.60# calculate pop growth rates for each year
wm.df <- data.frame(
year= 1988:2020,
m= apply(wm.post.lam, 1, median, na.rm=T),
lci95= apply(wm.post.lam, 1, HDInterval::hdi)[1,],
uci95= apply(wm.post.lam, 1, HDInterval::hdi)[2,],
lci80 = apply(wm.post.lam, 1, HDInterval::hdi, credMass=.8)[1,],
uci80 = apply(wm.post.lam, 1, HDInterval::hdi, credMass=.8)[2,],
pd = apply(wm.post.lam, 1, pd, null=1)
)
knitr::kable(wm.df, digits=c(0,3,3,3,3,3,2),
row.names=FALSE,
col.names = c("Year", "Median", "95% Lower HDI", "95% Upper HDI",
"80% Lower HDI", "80% Upper HDI", "Prob. direction"),
caption="Table S9. Population growth rates for each year.")| Year | Median | 95% Lower HDI | 95% Upper HDI | 80% Lower HDI | 80% Upper HDI | Prob. direction |
|---|---|---|---|---|---|---|
| 1988 | NA | NA | NA | NA | NA | 1.00 |
| 1989 | 1.017 | 0.952 | 1.080 | 0.976 | 1.056 | 0.71 |
| 1990 | 0.961 | 0.865 | 1.067 | 0.901 | 1.021 | 0.79 |
| 1991 | 1.000 | 0.939 | 1.064 | 0.962 | 1.039 | 0.50 |
| 1992 | 1.000 | 0.939 | 1.064 | 0.962 | 1.039 | 0.50 |
| 1993 | 1.000 | 0.940 | 1.065 | 0.963 | 1.040 | 0.50 |
| 1994 | 0.999 | 0.940 | 1.059 | 0.962 | 1.036 | 0.52 |
| 1995 | 1.000 | 0.943 | 1.067 | 0.963 | 1.040 | 0.50 |
| 1996 | 0.998 | 0.910 | 1.089 | 0.945 | 1.054 | 0.52 |
| 1997 | 1.011 | 0.950 | 1.072 | 0.975 | 1.049 | 0.65 |
| 1998 | 0.963 | 0.868 | 1.060 | 0.908 | 1.020 | 0.78 |
| 1999 | 0.999 | 0.939 | 1.061 | 0.962 | 1.036 | 0.51 |
| 2000 | 1.001 | 0.946 | 1.056 | 0.966 | 1.037 | 0.51 |
| 2001 | 1.000 | 0.943 | 1.067 | 0.963 | 1.040 | 0.50 |
| 2002 | 1.000 | 0.941 | 1.064 | 0.963 | 1.039 | 0.50 |
| 2003 | 1.002 | 0.945 | 1.061 | 0.965 | 1.036 | 0.52 |
| 2004 | 1.002 | 0.934 | 1.059 | 0.963 | 1.041 | 0.52 |
| 2005 | 0.995 | 0.945 | 1.047 | 0.962 | 1.027 | 0.58 |
| 2006 | 1.000 | 0.942 | 1.062 | 0.962 | 1.037 | 0.50 |
| 2007 | 1.008 | 0.959 | 1.061 | 0.979 | 1.042 | 0.63 |
| 2008 | 1.000 | 0.948 | 1.060 | 0.965 | 1.037 | 0.51 |
| 2009 | 1.000 | 0.943 | 1.067 | 0.963 | 1.039 | 0.50 |
| 2010 | 1.001 | 0.943 | 1.058 | 0.965 | 1.038 | 0.51 |
| 2011 | 0.971 | 0.889 | 1.059 | 0.921 | 1.024 | 0.76 |
| 2012 | 0.967 | 0.878 | 1.057 | 0.915 | 1.021 | 0.77 |
| 2013 | 0.983 | 0.895 | 1.077 | 0.927 | 1.036 | 0.66 |
| 2014 | 1.000 | 0.859 | 1.127 | 0.912 | 1.078 | 0.50 |
| 2015 | 0.983 | 0.895 | 1.077 | 0.927 | 1.036 | 0.66 |
| 2016 | 0.983 | 0.895 | 1.077 | 0.927 | 1.036 | 0.66 |
| 2017 | 0.983 | 0.895 | 1.077 | 0.927 | 1.036 | 0.66 |
| 2018 | 0.983 | 0.895 | 1.077 | 0.927 | 1.036 | 0.66 |
| 2019 | 0.998 | 0.870 | 1.120 | 0.919 | 1.075 | 0.51 |
| 2020 | 0.985 | 0.895 | 1.084 | 0.930 | 1.043 | 0.64 |
c3 <- makeTransparent("gray40", alpha=10)
c4 <- makeTransparent(4, alpha=180)
rownames(wts2) <- rownames(datl$y) <- 1988:2020
colnames(wts2) <- colnames(datl$y)
names(dimnames(datl$y)) <- c("Year", "Site")
all.df <- melt(as.matrix(datl$y), value.name="Count")
all.df$weights <- melt(as.matrix(wts2), value.name="weights")$weights
colnames(all.df) <- c("Year", "Site", "Count", "Weights")
labels <- c( "Bylot Island", "Karupelv Valley",
"Fennoscandia", "Wrangel Island")
sl <- all.df[!is.na(all.df$Count),]
p1 <- ggplot() + theme_minimal() +
coord_cartesian(clip = "off") +
geom_line(data=all.df,aes(x=Year,y=Site,color=Weights), linewidth=2)+
geom_point(data=sl,aes(x=Year,y=Site),shape="I", size=4)+
xlab("Year")+ylab("Site")+scale_color_viridis_c()+
scale_y_discrete(labels=labels)+
labs(color='Weight') +
xlim(1988, 2020) + theme(legend.position = "top") +
annotate(geom = "text", x = 1990, y = 6, label = "A", size=6)
p2 <- ggplot() + theme_minimal() +
geom_line(data=lam.post.long, aes(x=Year, y=value, group=Iter),
color="gray40", linewidth=0.5, alpha=0.05 ) +
geom_hline(yintercept=1, lwd=2, color="black", linetype="dashed") +
geom_line(data=wm.df, aes(x=year, y=lci95), color="deepskyblue3", linewidth=0.5 ) +
geom_line(data=wm.df, aes(x=year, y=uci95 ), color="deepskyblue3", linewidth=0.5) +
geom_line(data=wm.df, aes(x=year, y=lci80), color="deepskyblue3", linewidth=1 ) +
geom_line(data=wm.df, aes(x=year, y=uci80), color="deepskyblue3", linewidth=1 ) +
geom_line(data=wm.df, aes(x=year, y=m), color="deepskyblue3", linewidth=2 ) +
xlab("Year") +
ylab( expression(paste("Population growth (", lambda,")")) )+
xlim(1988, 2020) +
annotate(geom = "text", x = 1990, y = 1.1, label = "B", size=6)+
coord_cartesian(xlim=c(1988, 2020), ylim=c(0.85, 1.15))
ap12 <- align_plots(p1, p2, align="v", axis="l")
p12 <- plot_grid(ap12[[1]], ap12[[2]], nrow = 2, align="v")2.4.3 Derive percent change for IUCN Red List Criteria
2.4.3.1 Calculate percent change with 8-year generation time
#***********
#* Calculate proportion of base year
#* or percent change
#* Short generation time
#* 8-year generations
#*************
# calculate lambda from 1996
startyr <- 2020-(3*8) # 3 generations of 8 years
dnames <- list(year=1988:2020,
site=nms,
Iter=1:4000)
iucn <- array(NA, dim=dim(wm.post.lam), dimnames=dimnames(wm.post.lam))
iucn[9,] <- 1
for (t in 10:nrow(wm.post.lam)){
iucn[t,] <- iucn[(t-1),] * wm.post.lam[t,]
}
for (t in 8:1){
iucn[t,] <- iucn[(t+1),] / wm.post.lam[(t+1),]
}
# convert to percent change
iucn <- (iucn-1) * 100
iucn.post <- melt(iucn)
i.df <- data.frame(
year = 1988:2020,
m = apply(iucn, 1, median, na.rm=T),
lci95 = apply(iucn, 1, HDInterval::hdi)[1,],
uci95 = apply(iucn, 1, HDInterval::hdi)[2,],
lci80 = apply(iucn, 1, HDInterval::hdi, credMass=0.8)[1,],
uci80 = apply(iucn, 1, HDInterval::hdi, credMass=0.8)[2,],
pd = apply(iucn, 1, pd)
)
knitr::kable(i.df, digits=c(0,1,1,1,1,1,2),
row.names=FALSE,
col.names = c("Year", "Median", "95% Lower HDI", "95% Upper HDI",
"80% Lower HDI", "80% Upper HDI", "Prob. direction"),
caption="Table S10. Percent change since 1996.")| Year | Median | 95% Lower HDI | 95% Upper HDI | 80% Lower HDI | 80% Upper HDI | Prob. direction |
|---|---|---|---|---|---|---|
| 1988 | 2.0 | -30.2 | 49.4 | -25.6 | 23.3 | 0.54 |
| 1989 | 3.8 | -27.7 | 48.0 | -22.9 | 24.6 | 0.58 |
| 1990 | 0.0 | -30.7 | 35.8 | -21.5 | 20.6 | 0.50 |
| 1991 | 0.0 | -24.9 | 29.9 | -17.5 | 17.3 | 0.50 |
| 1992 | 0.2 | -19.5 | 24.0 | -13.8 | 13.7 | 0.51 |
| 1993 | 0.2 | -15.5 | 17.0 | -10.1 | 10.7 | 0.51 |
| 1994 | 0.2 | -10.7 | 13.2 | -7.2 | 7.6 | 0.51 |
| 1995 | 0.2 | -8.3 | 9.7 | -5.8 | 5.2 | 0.52 |
| 1996 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.00 |
| 1997 | 1.1 | -5.0 | 7.2 | -2.5 | 4.9 | 0.65 |
| 1998 | -2.6 | -15.6 | 11.5 | -10.7 | 5.8 | 0.66 |
| 1999 | -2.4 | -18.4 | 12.3 | -12.5 | 6.7 | 0.64 |
| 2000 | -2.2 | -20.2 | 17.6 | -13.2 | 10.6 | 0.60 |
| 2001 | -2.2 | -24.8 | 21.2 | -18.2 | 11.0 | 0.58 |
| 2002 | -2.3 | -26.1 | 29.8 | -18.8 | 16.5 | 0.56 |
| 2003 | -1.8 | -32.7 | 32.7 | -22.0 | 20.0 | 0.54 |
| 2004 | -1.8 | -36.3 | 38.7 | -25.0 | 23.0 | 0.54 |
| 2005 | -2.3 | -38.9 | 45.6 | -30.2 | 24.1 | 0.54 |
| 2006 | -2.1 | -45.4 | 49.0 | -32.0 | 28.8 | 0.54 |
| 2007 | -1.1 | -49.0 | 56.0 | -31.5 | 35.1 | 0.52 |
| 2008 | -1.2 | -50.8 | 64.2 | -40.9 | 32.1 | 0.52 |
| 2009 | -1.0 | -57.7 | 68.8 | -43.9 | 35.9 | 0.51 |
| 2010 | -1.1 | -60.9 | 77.4 | -48.9 | 37.3 | 0.51 |
| 2011 | -3.7 | -63.0 | 78.1 | -51.8 | 36.1 | 0.54 |
| 2012 | -7.0 | -65.9 | 81.0 | -55.9 | 33.3 | 0.57 |
| 2013 | -8.4 | -69.0 | 86.4 | -56.3 | 37.3 | 0.59 |
| 2014 | -8.8 | -76.0 | 92.0 | -58.1 | 43.1 | 0.58 |
| 2015 | -10.7 | -80.1 | 99.1 | -64.7 | 41.1 | 0.59 |
| 2016 | -12.0 | -81.8 | 108.5 | -68.9 | 41.5 | 0.60 |
| 2017 | -13.5 | -84.5 | 115.6 | -69.7 | 45.6 | 0.61 |
| 2018 | -15.3 | -87.9 | 122.2 | -72.8 | 46.5 | 0.61 |
| 2019 | -14.9 | -92.0 | 138.8 | -79.0 | 47.7 | 0.61 |
| 2020 | -16.3 | -92.5 | 151.4 | -81.2 | 50.4 | 0.61 |
p3 <- ggplot() + theme_minimal() +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-20, ymax=300), color="green4", fill="green4") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-30, ymax=-20), color="green3", fill="green3") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-50, ymax=-30), color="yellow", fill="yellow") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-80, ymax=-50), color="orange", fill="orange") +
geom_rect(aes(xmin=1986, xmax=2022, ymin=-100, ymax=-80), color="red", fill="red") +
geom_line(data=iucn.post, aes(x=Year, y=value, group=Iter),
color="gray40", size=0.5, alpha=0.05 ) +
geom_hline(yintercept=1, lwd=2, color="black", linetype="dashed") +
geom_line(data=i.df, aes(x=year, y=lci95), color="black", size=0.5 ) +
geom_line(data=i.df, aes(x=year, y=uci95 ), color="black", size=0.5) +
geom_line(data=i.df, aes(x=year, y=lci80), color="black", size=1 ) +
geom_line(data=i.df, aes(x=year, y=uci80), color="black", size=1 ) +
geom_line(data=i.df, aes(x=year, y=m), color="black", size=2 ) +
annotate(geom = "text", x = 1995, y = 25, label = "Least concern", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -25, label = "Near threatened", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -40, label = "Vulnerable", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -70, label = "Endangered", size=4, color="black") +
annotate(geom = "text", x = 1995, y = -90, label = "Critically endangered", size=4, color="black") +
xlab("") +
ylab( "Percent change")+
coord_cartesian(xlim=c(1988, 2020), ylim=c(-100, 60)) +
annotate(geom = "text", x = 1995, y = 55, label = "A", size=8) +
theme(axis.text= element_text(size=12),
axis.title=element_text(size=14),
axis.ticks = element_line(color = "black"),
plot.title = element_text(size=22)) +
labs(title="8-year generation time")
iucn2020 <- iucn.post[iucn.post$Year==2020, ]
i.df2020 <- i.df[i.df$year==2020,]
p4 <- ggplot() + theme_minimal() +
geom_rect(aes(ymin=0, ymax=2, xmin=-20, xmax=300), color="green4", fill="green4") +
geom_rect(aes(ymin=0, ymax=2, xmin=-30, xmax=-20), color="green3", fill="green3") +
geom_rect(aes(ymin=0, ymax=2, xmin=-50, xmax=-30), color="yellow", fill="yellow") +
geom_rect(aes(ymin=0, ymax=2, xmin=-80, xmax=-50), color="orange", fill="orange") +
geom_rect(aes(ymin=0, ymax=2, xmin=-100, xmax=-80), color="red", fill="red") +
stat_halfeye(data=iucn2020, aes(x = value, y=0), alpha=0.65,
slab_color="gray20", slab_fill="gray40",
point_interval="median_hdi", .width = c(0.80, 0.95),
point_size=5) +
scale_size_continuous(range = c(7, 15)) +
geom_vline(xintercept=0, lwd=2, color="black", linetype="dashed") +
xlab("") + ylab("") +
annotate(geom = "text", x = 55, y = 0.25, label = "B", size=8) +
theme(axis.text= element_text(size=12),
axis.title=element_text(size=14),
axis.ticks = element_line(color = "black")) +
coord_flip(xlim=c(-100, 60), ylim=c(0, 1), clip="on")
# calculate proportion of distribution in each category and mode
cuttab <- table(cut(iucn2020$value, breaks=c(-100,-80, -50, -30, -20, 300)))
prop <- cuttab/sum(cuttab)
df.iucn <- data.frame(IUCN.criteria=rev( c('Least concern', 'Near threatened', 'Vulnerable', 'Endangered', 'Critically endangered')),
Proportion.within=prop,
Proportion.within.and.worse= cumsum(prop)
)
df.iucn <- df.iucn[nrow(df.iucn):1,]
knitr::kable(df.iucn, digits=c(0, 0, 2, 2),
row.names=FALSE,
col.names = c("IUCN Category", "A2 Criteria", "Prop within", "Prop within and worse"),
caption="Table S11. Percent change over three generations.")| IUCN Category | A2 Criteria | Prop within | Prop within and worse |
|---|---|---|---|
| Least concern | (-20,300] | 0.52 | 1.00 |
| Near threatened | (-30,-20] | 0.08 | 0.48 |
| Vulnerable | (-50,-30] | 0.18 | 0.40 |
| Endangered | (-80,-50] | 0.20 | 0.22 |
| Critically endangered | (-100,-80] | 0.02 | 0.02 |
# calculate the mode
dens <- density(iucn2020$value)
# mode
dens$x[dens$y==max(dens$y)]## [1] -40.94626# median
median(iucn2020$value)## [1] -16.33671# mean
mean(iucn2020$value)## [1] 4.868512# 80% HDI
HDInterval::hdi(iucn2020$value, credMass = 0.80)## lower upper
## -81.19603 50.44685
## attr(,"credMass")
## [1] 0.8ap34 <- align_plots(p3, p4, align="h", axis="l")
p34 <- plot_grid(ap34[[1]], ap34[[2]],
nrow = 1, align="h", rel_widths = c(2, 1)) 2.4.3.2 Plot percent change
p34
Fig. S11. Percent change in abundance of Snowy Owl excluding data from Utqiagvik.