Description: This program uses a Bayesian species richness model and data augmentation to analyze species richness of each route and year.
model.richness <- function(){
# Priors for determining the probability of detection based on broadcast species
beta.prebroad ~ dt(0, pow(2.5, -2), 1) #p(detection) before any broadcast
beta.pacific ~ dt(0, pow(2.5, -2), 1)
beta.mottled ~ dt(0, pow(2.5, -2), 1)
beta.crested ~ dt(0, pow(2.5, -2), 1)
beta.bw ~ dt(0, pow(2.5, -2), 1)
beta.spectacled ~ dt(0, pow(2.5, -2), 1)
beta.whiskered ~ dt(0, pow(2.5, -2), 1)
beta.gbarred ~ dt(0, pow(2.5, -2), 1)
beta.stygian ~ dt(0, pow(2.5, -2), 1)
beta.ghorned ~ dt(0, pow(2.5, -2), 1)
# Hyperpriors for psi
lnrho.psi ~ dnorm(5,1)%_%T(0.01,10)
rho.psi <- exp(lnrho.psi)
for(hh in 1:n.route){
# Prior for omega, which is prob(species membership)
omega[hh] ~ dbeta(0.001, 1)
for(tt in 1:n.year){
#Hyperpriors for Psi ensure a flexible and uninformed prior
logit.mu.psi[hh,tt] ~ dnorm(0, 0.2)
mu.psi[hh,tt] <- exp(logit.mu.psi[hh,tt])/(1+exp(logit.mu.psi[hh,tt]))
for(ss in 1:n.species.aug){
# Hyperpriors based off route/year random effects, shared mean but
# flexible to relate to each species' prob(occupancy) psi
a.psi[hh,tt,ss] <- mu.psi[hh,tt]*rho.psi
b.psi[hh,tt,ss] <- rho.psi-(mu.psi[hh,tt]*rho.psi)
# Prior for Psi, which will vary by route (hh) and year (tt) and species (ss),
# but based on shared mean probability of occupancy for each route/year (mu.psi)
psi[hh,tt,ss] ~ dbeta(a.psi[hh,tt,ss], b.psi[hh,tt,ss])%_%T(0.0001,0.99)
for(ii in 1:n.survey){ # 1 to 3 surveys per year
for(jj in 1:10){ # 10 stations per route
for(kk in 1:n.broadcast){ # before or after broadcast
# Function that creates logistic equation for p(detection)
# p = p(detection), which varies by route, year, survey, station, and
# broadcast period (pre- or post-broadcast)
p[hh,tt,ss,ii,jj,kk] <- exp(logit.p[hh,tt,ss,ii,jj,kk])/
(1+exp(logit.p[hh,tt,ss,ii,jj,kk]))
# Logistic regression equation
logit.p[hh,tt,ss,ii,jj,kk] <-
beta.prebroad*ks.prebroad[hh,jj,kk]+
beta.pacific*ks.pacific[hh,jj,kk]+
beta.mottled*ks.mottled[hh,jj,kk]+
beta.crested*ks.crested[hh,jj,kk]+
beta.bw*ks.bw[hh,jj,kk]+
beta.spectacled*ks.spectacled[hh,jj,kk]+
beta.whiskered*ks.whiskered[hh,jj,kk]+
beta.gbarred*ks.gbarred[hh,jj,kk]+
beta.stygian*ks.stygian[hh,jj,kk]+
beta.ghorned*ks.ghorned[hh,jj,kk]
}
}
}
}
}
} # end priors
# Likelihood
for(hh in 1:n.route){ # 6 routes
for(tt in 1:n.year){ # all years
for(ss in 1:n.species.aug){ # all augmented species
# Is each species present in that route and year?
w[hh,tt,ss] ~ dbern(omega[hh])
for(ii in 1:n.survey){ # 1 to 3 surveys per year
# Occupancy by route and year and survey based on
# psi, probability of occupancy for each route/year
z[hh,tt,ss,ii] ~ dbern(psi[hh,tt,ss]*w[hh,tt,ss])
for(jj in 1:n.station){ # 10 stations per route
for(kk in 1:n.broadcast){ # before or after broadcast
# Detection by route, year, survey station, and broadcast period
# Binary observations by route, year, survey, station, pre/post broadcast
y[hh,tt,ss,ii,jj,kk] ~
dbern(eff.p[hh,tt,ss,ii,jj,kk])
# Effective p(detection), which depends on occupany (z) = 1 for
# that route/year/survey/station
eff.p[hh,tt,ss,ii,jj,kk] <- p[hh,tt,ss,ii,jj,kk]*z[hh,tt,ss,ii]
}
}
}
}
}
} # end likelihood
# Derived quantities
for(ss in 1:n.species.aug){
spp.occ[ss] <- sum(z[,,ss,]) # number of occupied routes/years among sampled
species.present[ss] <- ifelse(spp.occ[ss]>0, 1, 0) # if each species exists
}
Nsmall <- sum(species.present) # small estimate of number of species present
for(hh in 1:n.route){
for(tt in 1:n.year){
richness[hh,tt] <- sum(z[hh,tt,,]) # number of species at each route/year
}
}
}
Session Info
devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────────────────────
## setting value
## version R version 4.0.2 (2020-06-22)
## os macOS Mojave 10.14.6
## system x86_64, darwin17.0
## ui RStudio
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/Chicago
## date 2020-11-10
##
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────
## package * version date lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 4.0.2)
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.2)
## backports 1.1.8 2020-06-17 [1] CRAN (R 4.0.2)
## boot 1.3-25 2020-04-26 [1] CRAN (R 4.0.2)
## callr 3.4.3 2020-03-28 [1] CRAN (R 4.0.2)
## cli 2.0.2 2020-02-28 [1] CRAN (R 4.0.2)
## coda * 0.19-3 2019-07-05 [1] CRAN (R 4.0.2)
## colorspace 1.4-1 2019-03-18 [1] CRAN (R 4.0.2)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 4.0.2)
## desc 1.2.0 2018-05-01 [1] CRAN (R 4.0.2)
## devtools * 2.3.1 2020-07-21 [1] CRAN (R 4.0.2)
## digest 0.6.25 2020-02-23 [1] CRAN (R 4.0.2)
## dplyr * 1.0.1 2020-07-31 [1] CRAN (R 4.0.2)
## ellipsis 0.3.1 2020-05-15 [1] CRAN (R 4.0.2)
## evaluate 0.14 2019-05-28 [1] CRAN (R 4.0.1)
## ezknitr * 0.6 2016-09-16 [1] CRAN (R 4.0.2)
## fansi 0.4.1 2020-01-08 [1] CRAN (R 4.0.2)
## farver 2.0.3 2020-01-16 [1] CRAN (R 4.0.2)
## fs 1.5.0 2020-07-31 [1] CRAN (R 4.0.2)
## generics 0.0.2 2018-11-29 [1] CRAN (R 4.0.2)
## ggplot2 * 3.3.2 2020-06-19 [1] CRAN (R 4.0.2)
## ggthemes * 4.2.0 2019-05-13 [1] CRAN (R 4.0.2)
## glue 1.4.1 2020-05-13 [1] CRAN (R 4.0.2)
## gtable 0.3.0 2019-03-25 [1] CRAN (R 4.0.2)
## highr 0.8 2019-03-20 [1] CRAN (R 4.0.2)
## knitr * 1.29 2020-06-23 [1] CRAN (R 4.0.2)
## labeling 0.3 2014-08-23 [1] CRAN (R 4.0.2)
## lattice 0.20-41 2020-04-02 [1] CRAN (R 4.0.2)
## lifecycle 0.2.0 2020-03-06 [1] CRAN (R 4.0.2)
## magrittr 1.5 2014-11-22 [1] CRAN (R 4.0.2)
## markdown 1.1 2019-08-07 [1] CRAN (R 4.0.2)
## MCMCvis * 0.14.0 2020-03-25 [1] CRAN (R 4.0.2)
## memoise 1.1.0 2017-04-21 [1] CRAN (R 4.0.2)
## mime 0.9 2020-02-04 [1] CRAN (R 4.0.2)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.0.2)
## pillar 1.4.6 2020-07-10 [1] CRAN (R 4.0.2)
## pkgbuild 1.1.0 2020-07-13 [1] CRAN (R 4.0.2)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.0.2)
## pkgload 1.1.0 2020-05-29 [1] CRAN (R 4.0.2)
## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.0.2)
## processx 3.4.3 2020-07-05 [1] CRAN (R 4.0.2)
## ps 1.3.3 2020-05-08 [1] CRAN (R 4.0.2)
## purrr 0.3.4 2020-04-17 [1] CRAN (R 4.0.2)
## R.methodsS3 1.8.0 2020-02-14 [1] CRAN (R 4.0.2)
## R.oo 1.23.0 2019-11-03 [1] CRAN (R 4.0.2)
## R.utils 2.9.2 2019-12-08 [1] CRAN (R 4.0.2)
## R2jags * 0.6-1 2020-04-27 [1] CRAN (R 4.0.2)
## R2WinBUGS 2.1-21 2015-07-30 [1] CRAN (R 4.0.2)
## R6 2.4.1 2019-11-12 [1] CRAN (R 4.0.2)
## remotes 2.2.0 2020-07-21 [1] CRAN (R 4.0.2)
## rjags * 4-10 2019-11-06 [1] CRAN (R 4.0.2)
## rlang 0.4.7 2020-07-09 [1] CRAN (R 4.0.2)
## rprojroot 1.3-2 2018-01-03 [1] CRAN (R 4.0.2)
## rstudioapi 0.11 2020-02-07 [1] CRAN (R 4.0.2)
## scales 1.1.1 2020-05-11 [1] CRAN (R 4.0.2)
## sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 4.0.2)
## stringi 1.4.6 2020-02-17 [1] CRAN (R 4.0.2)
## stringr 1.4.0 2019-02-10 [1] CRAN (R 4.0.2)
## testthat 2.3.2 2020-03-02 [1] CRAN (R 4.0.2)
## tibble 3.0.3 2020-07-10 [1] CRAN (R 4.0.2)
## tidyr * 1.1.1 2020-07-31 [1] CRAN (R 4.0.2)
## tidyselect 1.1.0 2020-05-11 [1] CRAN (R 4.0.2)
## usethis * 1.6.1 2020-04-29 [1] CRAN (R 4.0.2)
## vctrs 0.3.2 2020-07-15 [1] CRAN (R 4.0.2)
## withr 2.2.0 2020-04-20 [1] CRAN (R 4.0.2)
## xfun 0.16 2020-07-24 [1] CRAN (R 4.0.2)
##
## [1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library
This document was “spun” with:
ezspin(file = “programs/b_richness_model_global.R”, out_dir = “output”, fig_dir = “figures”, keep_md = F)