Description: This program takes posteriors from a Bayesian occupancy model to analyze occupancy of each route and year for each species of owl and prepares it for plotting.
Load libraries
library(knitr) # documentation-related
library(ezknitr) # documentation-related
library(devtools) # documentation-related
# analysis-related
library(R2jags)
library(tidyr)
library(dplyr)
library(MCMCvis)
Clear environment and set seed
Note: for reproducibility, it is important to include these. Clearing the environment ensures that you have specified all pertinent files that need to be loaded, and setting the seed ensures that your analysis is repeatable
remove(list = ls())
set.seed(258854)
Jagsout FerPy
load(file = "data/output_data/ferpy_jagsout.Rdata")
Survey list
(route.names <- c("EI1", "EI2", "M2", "N1", "N2") )
## [1] "EI1" "EI2" "M2" "N1" "N2"
route.index <- 1:length(route.names)
(year.names <- 2003:2013)
## [1] 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013
year.index <- 1:length(year.names)
Create traceplots of parameters (creates PDF)
Ferpy
# MCMCtrace(ferpy.jagsout,
# pdf = TRUE,
# open_pdf = FALSE,
# filename = "ferpyTrace",
# wd = "output")
and year
For all analysis, remove years and routes that were not surveyed
(exclude.byrow <- c(3,16:20,33,37,47))
## [1] 3 16 17 18 19 20 33 37 47
all.rows <- 1:(length(route.names)*length(year.index))
(include.byrow <- all.rows[!all.rows %in% exclude.byrow])
## [1] 1 2 4 5 6 7 8 9 10 11 12 13 14 15 21 22 23 24 25 26 27 28 29 30 31 32
## [27] 34 35 36 38 39 40 41 42 43 44 45 46 48 49 50 51 52 53 54 55
Ferpy
psi.post.ferpy <- MCMCsummary(ferpy.jagsout,
params = "psi",
Rhat = TRUE,
n.eff = TRUE,
probs = c(0.05, 0.5, 0.95))
psi.post.ferpy$Year <- rep(year.names, each = length(route.index))
psi.post.ferpy$Route <- rep(c("EI.1", "EI.2", "M.2", "N.1", "N.2"), length(year.index))
psi.post.ferpy$Species <- "FerPy"
psi.post.ferpy <- psi.post.ferpy[include.byrow,]
Populate Region
psi.post.ferpy$Region[psi.post.ferpy$Route == "EI.1" |
psi.post.ferpy$Route == "EI.2"] <- "El Imposible"
psi.post.ferpy$Region[psi.post.ferpy$Route == "N.1" |
psi.post.ferpy$Route == "N.2"] <- "Nancuchiname"
psi.post.ferpy$Region[psi.post.ferpy$Route == "M.1" |
psi.post.ferpy$Route == "M.2"] <- "Montecristo"
colnames(psi.post.ferpy)
## [1] "mean" "sd" "5%" "50%" "95%" "Rhat" "n.eff"
## [8] "Year" "Route" "Species" "Region"
colnames(psi.post.ferpy) <- c("Psi.mean", "Psi.sd", "Psi.LL05", "Psi.median",
"Psi.UL95", "Psi.Rhat", "Psi.neff", "Year", "Route",
"Species", "Region")
Need to delete a few years by route:
all.years <- 1:length(year.index)
exc.EI.1 <- 4
inc.EI.1 <- all.years[!all.years %in% exc.EI.1]
exc.EI.2 <- c(4,8,10)
inc.EI.2 <- all.years[!all.years %in% exc.EI.2]
exc.M.2 <- c(1,4,7,11)
inc.M.2 <- all.years[!all.years %in% exc.M.2]
exc.N.1 <- 4
inc.N.1 <- all.years[!all.years %in% exc.N.1]
exc.N.2 <- 4
inc.N.2 <- all.years[!all.years %in% exc.N.2]
Psi posteriors for Ferpy
ferpy.chains <- MCMCpstr(ferpy.jagsout,
params = "psi",
type = "chains")
Calculate by species and route
psi.means.ferpy <- as.data.frame(matrix(NA, ncol = 5, nrow = length(route.names)))
colnames(psi.means.ferpy) <- c("Species", "Route", "Psi.LL05", "Psi.median", "Psi.UL95")
psi.means.ferpy$Species <- "FerPy"
psi.means.ferpy$Route <- route.names
Populate quantile results
psi.means.ferpy[psi.means.ferpy$Route == "EI1",3:5] <-
quantile(ferpy.chains$psi[1,inc.EI.1,], probs = c(0.05, 0.5, 0.95))
psi.means.ferpy[psi.means.ferpy$Route == "EI2",3:5] <-
quantile(ferpy.chains$psi[2,inc.EI.2,], probs = c(0.05, 0.5, 0.95))
psi.means.ferpy[psi.means.ferpy$Route == "M2",3:5] <-
quantile(ferpy.chains$psi[3,inc.M.2,], probs = c(0.05, 0.5, 0.95))
psi.means.ferpy[psi.means.ferpy$Route == "N1",3:5] <-
quantile(ferpy.chains$psi[4,inc.N.1,], probs = c(0.05, 0.5, 0.95))
psi.means.ferpy[psi.means.ferpy$Route == "N2",3:5] <-
quantile(ferpy.chains$psi[5,inc.N.2,], probs = c(0.05, 0.5, 0.95))
Populate Region
psi.means.ferpy$Region[psi.means.ferpy$Route == "EI1" |
psi.means.ferpy$Route == "EI2"] <- "El Imposible"
psi.means.ferpy$Region[psi.means.ferpy$Route == "N1" |
psi.means.ferpy$Route == "N2"] <- "Nancuchiname"
psi.means.ferpy$Region[psi.means.ferpy$Route == "M1" |
psi.means.ferpy$Route == "M2"] <- "Montecristo"
Probability of detection was a function of broadcast species
p.det.post.ferpy <- MCMCsummary(ferpy.jagsout,
params = grep("^beta",
ferpy.jagsout$parameters.to.save,
value = T),
Rhat = TRUE,
n.eff = TRUE,
probs = c(0.05, 0.5, 0.95))
p.det.post.ferpy$Species <- "FerPy"
p.det.post.ferpy$broadcast.param <- grep("^beta",
ferpy.jagsout$parameters.to.save,
value = T)
p.det.post.ferpy
## mean sd 5% 50% 95% Rhat
## beta.prebroad -1.5368169 0.1383735 -1.7668055 -1.53441128 -1.31240155 1.00
## beta.pacific -0.2920913 0.2345886 -0.6818964 -0.29062581 0.09067835 1.00
## beta.mottled -0.8136654 0.2563250 -1.2372612 -0.80979594 -0.40180457 1.00
## beta.crested -1.4264304 0.2990514 -1.9337012 -1.41654697 -0.95257623 1.00
## beta.bw -1.6110965 0.3176010 -2.1485495 -1.59902930 -1.10844683 1.00
## beta.spectacled -1.1792589 0.2785869 -1.6506382 -1.17136430 -0.72815511 1.00
## beta.whiskered -3.0351227 17.3980611 -23.1605623 -0.76724541 10.22943412 1.12
## beta.gbarred -1.1561691 22.0679651 -16.2930367 -0.10981085 13.89062990 1.10
## beta.stygian 1.1614821 22.8662784 -14.5532606 -0.08019815 16.55408854 1.04
## beta.ghorned 11.2242188 72.3037207 -15.0974908 0.11283999 38.00314765 1.38
## n.eff Species broadcast.param
## beta.prebroad 16382 FerPy beta.prebroad
## beta.pacific 16540 FerPy beta.pacific
## beta.mottled 15910 FerPy beta.mottled
## beta.crested 16290 FerPy beta.crested
## beta.bw 16095 FerPy beta.bw
## beta.spectacled 16535 FerPy beta.spectacled
## beta.whiskered 819 FerPy beta.whiskered
## beta.gbarred 571 FerPy beta.gbarred
## beta.stygian 619 FerPy beta.stygian
## beta.ghorned 777 FerPy beta.ghorned
Merge data together and change column names
colnames(p.det.post.ferpy)
## [1] "mean" "sd" "5%" "50%"
## [5] "95%" "Rhat" "n.eff" "Species"
## [9] "broadcast.param"
colnames(p.det.post.ferpy) <- c("mean", "sd", "LL05", "median", "UL95", "Rhat",
"n.eff", "Species", "broadcast.param")
Transform to probability scale
p.det.post.ferpy$LL05.plogis <- plogis(p.det.post.ferpy$LL05)
p.det.post.ferpy$median.plogis <- plogis(p.det.post.ferpy$median)
p.det.post.ferpy$UL95.plogis <- plogis(p.det.post.ferpy$UL95)
Save broadcast species as a factor, ordered by confidence
p.det.post.ferpy$Broadcast <- factor(p.det.post.ferpy$broadcast.param,
levels = c("beta.prebroad", "beta.mottled",
"beta.pacific", "beta.crested",
"beta.bw", "beta.spectacled",
"beta.whiskered", "beta.gbarred",
"beta.stygian", "beta.ghorned"),
labels = c("Pre-broadcast", "Mottled Owl",
"Pacific Screech-owl", "Crested Owl",
"Black-and-white Owl",
"Spectacled Owl",
"Whiskered Screech-owl", "Fulvous Owl",
"Stygian Owl", "Great Horned Owl"))
Psi Posteriors by year and route
save(psi.post.ferpy, file = "data/plotting_data/ferpy_psi_posteriors_RtYr.Rdata")
Psi posteriors across years by species and route
save(psi.means.ferpy, file = "data/plotting_data/ferpy_psi_posteriors_RtSpp.Rdata")
Probability of detection by broadcast species and species of analysis
save(p.det.post.ferpy, file = "data/plotting_data/ferpy_p_detection_posteriors.Rdata")
Session Info
devtools::session_info()
## ─ Session info ─────────────────────────────────────────────────────────────────
## setting value
## version R version 4.0.3 (2020-10-10)
## os macOS Big Sur 10.16
## system x86_64, darwin17.0
## ui RStudio
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/Chicago
## date 2021-01-14
##
## ─ 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.2.1 2020-12-09 [1] CRAN (R 4.0.2)
## base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.0.2)
## boot 1.3-25 2020-04-26 [1] CRAN (R 4.0.3)
## callr 3.5.1 2020-10-13 [1] CRAN (R 4.0.2)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.0.2)
## checkmate 2.0.0 2020-02-06 [1] CRAN (R 4.0.2)
## chron * 2.3-56 2020-08-18 [1] CRAN (R 4.0.2)
## cli 2.2.0 2020-11-20 [1] CRAN (R 4.0.2)
## cluster 2.1.0 2019-06-19 [1] CRAN (R 4.0.3)
## coda * 0.19-4 2020-09-30 [1] CRAN (R 4.0.2)
## colorspace 2.0-0 2020-11-11 [1] CRAN (R 4.0.2)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 4.0.2)
## data.table 1.13.6 2020-12-30 [1] CRAN (R 4.0.2)
## desc 1.2.0 2018-05-01 [1] CRAN (R 4.0.2)
## devtools * 2.3.2 2020-09-18 [1] CRAN (R 4.0.2)
## digest 0.6.27 2020-10-24 [1] CRAN (R 4.0.2)
## dplyr * 1.0.2 2020-08-18 [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)
## forcats 0.5.0 2020-03-01 [1] CRAN (R 4.0.2)
## foreign 0.8-80 2020-05-24 [1] CRAN (R 4.0.3)
## Formula * 1.2-4 2020-10-16 [1] CRAN (R 4.0.2)
## fs 1.5.0 2020-07-31 [1] CRAN (R 4.0.2)
## gdata * 2.18.0 2017-06-06 [1] CRAN (R 4.0.2)
## generics 0.1.0 2020-10-31 [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.2 2020-08-27 [1] CRAN (R 4.0.2)
## gridExtra 2.3 2017-09-09 [1] CRAN (R 4.0.2)
## gtable 0.3.0 2019-03-25 [1] CRAN (R 4.0.2)
## gtools 3.8.2 2020-03-31 [1] CRAN (R 4.0.2)
## haven 2.3.1 2020-06-01 [1] CRAN (R 4.0.2)
## highr 0.8 2019-03-20 [1] CRAN (R 4.0.2)
## Hmisc * 4.4-2 2020-11-29 [1] CRAN (R 4.0.2)
## hms 0.5.3 2020-01-08 [1] CRAN (R 4.0.2)
## htmlTable 2.1.0 2020-09-16 [1] CRAN (R 4.0.2)
## htmltools 0.5.0 2020-06-16 [1] CRAN (R 4.0.2)
## htmlwidgets 1.5.3 2020-12-10 [1] CRAN (R 4.0.2)
## ImportExport * 1.3 2020-09-21 [1] CRAN (R 4.0.2)
## jpeg 0.1-8.1 2019-10-24 [1] CRAN (R 4.0.2)
## knitr * 1.30 2020-09-22 [1] CRAN (R 4.0.2)
## labeling 0.4.2 2020-10-20 [1] CRAN (R 4.0.2)
## lattice * 0.20-41 2020-04-02 [1] CRAN (R 4.0.3)
## latticeExtra 0.6-29 2019-12-19 [1] CRAN (R 4.0.2)
## lifecycle 0.2.0 2020-03-06 [1] CRAN (R 4.0.2)
## lubridate * 1.7.9.2 2020-11-13 [1] CRAN (R 4.0.2)
## magrittr 2.0.1 2020-11-17 [1] CRAN (R 4.0.2)
## markdown 1.1 2019-08-07 [1] CRAN (R 4.0.2)
## Matrix 1.2-18 2019-11-27 [1] CRAN (R 4.0.3)
## 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)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.0.2)
## nnet 7.3-14 2020-04-26 [1] CRAN (R 4.0.3)
## pillar 1.4.7 2020-11-20 [1] CRAN (R 4.0.2)
## pkgbuild 1.2.0 2020-12-15 [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)
## png 0.1-7 2013-12-03 [1] CRAN (R 4.0.2)
## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.0.2)
## processx 3.4.5 2020-11-30 [1] CRAN (R 4.0.2)
## ps 1.5.0 2020-12-05 [1] CRAN (R 4.0.2)
## purrr 0.3.4 2020-04-17 [1] CRAN (R 4.0.2)
## R.methodsS3 1.8.1 2020-08-26 [1] CRAN (R 4.0.2)
## R.oo 1.24.0 2020-08-26 [1] CRAN (R 4.0.2)
## R.utils 2.10.1 2020-08-26 [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.5.0 2020-10-28 [1] CRAN (R 4.0.2)
## RColorBrewer 1.1-2 2014-12-07 [1] CRAN (R 4.0.2)
## Rcpp 1.0.5 2020-07-06 [1] CRAN (R 4.0.2)
## readxl * 1.3.1 2019-03-13 [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.9 2020-11-26 [1] CRAN (R 4.0.2)
## rmarkdown 2.6 2020-12-14 [1] CRAN (R 4.0.2)
## RODBC * 1.3-17 2020-05-11 [1] CRAN (R 4.0.2)
## rpart 4.1-15 2019-04-12 [1] CRAN (R 4.0.3)
## rprojroot 2.0.2 2020-11-15 [1] CRAN (R 4.0.2)
## rstudioapi 0.13 2020-11-12 [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.5.3 2020-09-09 [1] CRAN (R 4.0.2)
## stringr 1.4.0 2019-02-10 [1] CRAN (R 4.0.2)
## survival * 3.2-7 2020-09-28 [1] CRAN (R 4.0.3)
## testthat 3.0.1 2020-12-17 [1] CRAN (R 4.0.2)
## tibble 3.0.4 2020-10-12 [1] CRAN (R 4.0.2)
## tidyr * 1.1.2 2020-08-27 [1] CRAN (R 4.0.2)
## tidyselect 1.1.0 2020-05-11 [1] CRAN (R 4.0.2)
## usethis * 2.0.0 2020-12-10 [1] CRAN (R 4.0.2)
## utf8 1.1.4 2018-05-24 [1] CRAN (R 4.0.2)
## vctrs 0.3.6 2020-12-17 [1] CRAN (R 4.0.2)
## withr 2.3.0 2020-09-22 [1] CRAN (R 4.0.2)
## writexl 1.3.1 2020-08-26 [1] CRAN (R 4.0.2)
## xfun 0.19 2020-10-30 [1] CRAN (R 4.0.2)
## yaml 2.2.1 2020-02-01 [1] CRAN (R 4.0.2)
##
## [1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library
This document was “spun” with:
ezspin(file = “programs/e01_processing_results_ferpy.R”, out_dir = “output”, fig_dir = “figures”, keep_md = F)