Description: This program graphs occupancy analysis results for owls of El Salvador. Plots show the posteriors of Psi and P(detection).
Load libraries
library(knitr) # documentation-related
library(ezknitr) # documentation-related
library(devtools) # documentation-related
# analysis-related
library(ggplot2)
library(ggthemes)
library(maptools)
library(plyr)
library(rgdal)
library(raster)
library(mapproj)
library(stringr)
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(587453)
Psi Posteriors by year and route
load(file = "data/plotting_data/specd_psi_posteriors_RtYr.Rdata")
load(file = "data/plotting_data/ferpy_psi_posteriors_RtYr.Rdata")
load(file = "data/plotting_data/mottd_psi_posteriors_RtYr.Rdata")
load(file = "data/plotting_data/richness_psi_posteriors_RtYr.Rdata")
Psi posteriors across years by species and route
load(file = "data/plotting_data/specd_psi_posteriors_RtSpp.Rdata")
load(file = "data/plotting_data/ferpy_psi_posteriors_RtSpp.Rdata")
load(file = "data/plotting_data/mottd_psi_posteriors_RtSpp.Rdata")
Probability of detection by broadcast species and species of analysis
load(file = "data/plotting_data/specd_p_detection_posteriors.Rdata")
load(file = "data/plotting_data/ferpy_p_detection_posteriors.Rdata")
load(file = "data/plotting_data/mottd_p_detection_posteriors.Rdata")
load(file = "data/plotting_data/richness_p_detection_posteriors.Rdata")
Species accounts
load(file = "data/output_data/richness_species_accounts.Rdata")
Add in blanks for M1 for FerPy and Specd
psi.means.specd <- rbind(psi.means.specd, 0)
psi.means.specd$Route[max(nrow(psi.means.specd))] <- "M1"
psi.means.specd$Region[max(nrow(psi.means.specd))] <- "Montecristo"
psi.means.specd$Species[max(nrow(psi.means.specd))] <- "Specd"
psi.means.ferpy <- rbind(psi.means.ferpy, 0)
psi.means.ferpy$Route[max(nrow(psi.means.ferpy))] <- "M1"
psi.means.ferpy$Region[max(nrow(psi.means.ferpy))] <- "Montecristo"
psi.means.ferpy$Species[max(nrow(psi.means.ferpy))] <- "FerPy"
Merge data
psi.means <- rbind(psi.means.ferpy, psi.means.specd, psi.means.mottd)
psi.means$code <- ifelse(psi.means$Species == "FerPy", "Ferruginous Pygmy-Owl",
ifelse(psi.means$Species == "Mottd", "Mottled Owl", "Spectacled Owl"))
Last minute fix to route names
table(psi.means$Route)
##
## EI1 EI2 M1 M2 N1 N2
## 3 3 3 3 3 3
psi.means$Route <- ifelse(psi.means$Route == "EI1", "EI-1",
ifelse(psi.means$Route == "EI2", "EI-2",
ifelse(psi.means$Route == "M1", "M-1",
ifelse(psi.means$Route == "M2", "M-2",
ifelse(psi.means$Route == "N1", "N-1",
"N-2")))))
ggplot(data = psi.means, aes(x = Route, y = Psi.median, group = str_wrap(code,10)))+
geom_bar(stat = "identity", aes(fill = str_wrap(code,10)), position= position_dodge())+
geom_linerange(aes(ymin = Psi.LL05, ymax = Psi.UL95), position = position_dodge(width = 0.9))+
facet_wrap(~Route, nrow = 3,scales = "free_x")+
ylim(c(0,1))+
scale_fill_manual(values = c("#cccccc", "#969696", "#525252"))+
theme_minimal()+
ylab("Probability of occupancy")+
xlab("Route")+
theme(legend.position = "bottom")+
theme(legend.text=element_text(size=6),
legend.title = element_blank(),
axis.text.y = element_text(size=6),
axis.text.x = element_blank(),
axis.title.y = element_text(size=8),
axis.title.x = element_blank(),
strip.text = element_text(size = 6),
legend.key.size = unit(1, "line"))
Merge psi posts
psi.post.all <- rbind(psi.post.ferpy, psi.post.mottd, psi.post.specd)
psi.post.all <- rbind(psi.post.all, 0)
psi.post.all$Year[psi.post.all$Psi.mean == 0] <- 2000
psi.post.all$Species[psi.post.all$Psi.mean == 0] <- "FerPy"
psi.post.all$Route[psi.post.all$Psi.mean == 0] <- "M.1"
psi.post.all <- rbind(psi.post.all, 0.00001)
psi.post.all$Year[psi.post.all$Psi.mean == 0.00001] <- 2000
psi.post.all$Species[psi.post.all$Psi.mean == 0.00001] <- "Specd"
psi.post.all$Route[psi.post.all$Psi.mean == 0.00001] <- "M.1"
psi.post.all$code <- ifelse(psi.post.all$Species == "FerPy", "Ferruginous Pygmy-Owl",
ifelse(psi.post.all$Species == "Mottd", "Mottled Owl", "Spectacled Owl"))
Last minute fix for Route names
psi.post.all$Route <- ifelse(psi.post.all$Route == "EI.1", "EI-1",
ifelse(psi.post.all$Route == "EI.2", "EI-2",
ifelse(psi.post.all$Route == "M.1", "M-1",
ifelse(psi.post.all$Route == "M.2", "M-2",
ifelse(psi.post.all$Route == "N.1", "N-1",
"N-2")))))
table(psi.post.all$Route)
##
## EI-1 EI-2 M-1 M-2 N-1 N-2
## 30 24 6 23 30 30
ggplot(data = psi.post.all,
aes(x = Year, y = Psi.median))+
geom_linerange(aes(ymin = Psi.LL05, ymax = Psi.UL95), color = "#b5b5b5")+
geom_line()+
geom_point()+
#facet_grid(code~Route)+
facet_grid(Route~code)+
theme_minimal()+
xlab("Year")+
ylab("Probability of occupancy")+
scale_x_continuous(breaks = 2003:2013, limits = c(2003,2013))+
theme(axis.text.x = element_text(size=6, angle = 90),
axis.text.y = element_text(size=6),
axis.title = element_text(size=8),
#axis.title.x = element_blank(),
panel.border = element_rect(linetype = "solid", fill = NA, color = "#969696"),
strip.text = element_text(size =6, margin = margin(0,0.1,0.1,0.1, "cm")),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_segment).
## geom_path: Each group consists of only one observation. Do you need to adjust the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust the group aesthetic?
## Warning: Removed 2 rows containing missing values (geom_point).
Merge data
p.det.post.richness$Species <- "All Species"
p.det.post <- rbind(p.det.post.ferpy, p.det.post.mottd,
p.det.post.specd,p.det.post.richness)
table(p.det.post$Broadcast)
##
## Pre-broadcast Mottled Owl Pacific Screech-owl Crested Owl
## 4 4 4 4
## Black-and-white Owl Spectacled Owl Whiskered Screech-owl Fulvous Owl
## 4 4 4 4
## Stygian Owl Great Horned Owl
## 4 4
table(p.det.post$Species)
##
## All Species FerPy Mottd Specd
## 10 10 10 10
p.det.post$Species[p.det.post$Species == "FerPy"] <- "Ferruginous Pygmy-Owl"
p.det.post$Species[p.det.post$Species == "Mottd"] <- "Mottled Owl"
p.det.post$Species[p.det.post$Species == "Specd"] <- "Spectacled Owl"
Probability of detection was a function of what broadcast species was used, with a constant probability of detection for all pre-broadcast time periods.
ggplot(data = p.det.post, aes(y = median.plogis, x = Broadcast, group = Species))+
scale_fill_manual(values = c("darkgrey", "#7fc97f", "#beaed4", "#fdc086"))+
geom_hline(data = p.det.post[p.det.post$broadcast.param == "beta.prebroad",],
aes(yintercept = LL05.plogis), color = "#999999")+
geom_hline(data = p.det.post[p.det.post$broadcast.param == "beta.prebroad",],
aes(yintercept = UL95.plogis), color = "#999999")+
geom_point(size = 0.5)+
geom_linerange(aes(ymin = LL05.plogis, ymax = UL95.plogis),
position = position_dodge(0.9))+
facet_wrap(~Species, nrow = 2)+
ylab("Probability of detection")+
xlab("Broadcast species")+
theme_minimal()+
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 5),
legend.position = "none",
axis.text.y = element_text(size=6),
axis.title = element_text(size=6),
panel.border = element_rect(linetype = "solid", fill = NA, color = "#969696"),
strip.text = element_text(size =6, margin = margin(0,0.1,0.1,0.1, "cm")),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())
Last minute fix for Route names
table(richness.RtYr.post$Route)
##
## EI1 EI2 M1 M2 N1 N2
## 10 8 4 7 10 10
richness.RtYr.post$Route <- ifelse(richness.RtYr.post$Route == "EI1", "EI-1",
ifelse(richness.RtYr.post$Route == "EI2", "EI-2",
ifelse(richness.RtYr.post$Route == "M1", "M-1",
ifelse(richness.RtYr.post$Route == "M2", "M-2",
ifelse(richness.RtYr.post$Route == "N1", "N-1",
"N-2")))))
table(richness.RtYr.post$Route)
##
## EI-1 EI-2 M-1 M-2 N-1 N-2
## 10 8 4 7 10 10
All 6 routes
ggplot(data = richness.RtYr.post,
aes(x = Year, y = Richness.median, group = Route))+
geom_ribbon(aes(ymax = richness.detected, ymin = 0), fill = "#cccccc")+
geom_ribbon(aes(ymin = Richness.median,
ymax = Richness.UL95), fill = "#999999")+
geom_ribbon(aes(ymin = richness.detected,
ymax = Richness.median), fill = "#999999")+
geom_pointrange(aes(ymin = Richness.LL05, ymax = Richness.UL95))+
geom_line()+
facet_wrap(~Route, nrow = 3)+
theme_minimal()+
xlab("Year")+
ylab("Species richness")+
scale_x_continuous(breaks = 2003:2013)+
scale_y_continuous(breaks = seq(0,14,by=2), limits = c(0,6))+
theme(axis.text.x = element_text(size=6, angle = 90),
axis.text.y = element_text(size=6),
axis.title = element_text(size=8),
panel.border = element_rect(linetype = "solid",
fill = NA, color = "#969696"),
strip.text = element_text(size =6, margin = margin(0,0.1,0.1,0.1, "cm")),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())
Extract country data
ElSalvador <- getData("GADM", country = "SV", level = 0)
Guatemala <- getData("GADM", country = "GT", level = 0)
Honduras <- getData("GADM", country = "HN", level = 0)
# Process into df for ggplot
ElSalvador <- fortify(ElSalvador)
## Regions defined for each Polygons
Guatemala <- fortify(Guatemala)
## Regions defined for each Polygons
Honduras <- fortify(Honduras)
## Regions defined for each Polygons
# Range for mapping
range(ElSalvador$long)
## [1] -90.12486 -87.68375
range(ElSalvador$lat)
## [1] 13.15264 14.45055
Map of protected areas
ggplot(data = ElSalvador, aes(x = long, y = lat, group = group))+
geom_polygon(fill = "grey", color = "#848484")+
geom_polygon(data = Honduras, aes(x = long, y = lat, group = group),
fill = "lightgrey", color = "#848484")+
geom_polygon(data = Guatemala, aes(x = long, y = lat, group = group),
fill = "lightgrey", color = "#848484")+
coord_map(xlim = c(-90.2, -87.2), ylim = c(13.0, 14.5))+
# El Imposible
annotate(geom = "text", y = 13.74, x = -89.8, label = "EINP", size = 4)+
annotate(geom = "point", y = 13.83, x = -89.97, size = 5, shape = 21, fill = "lightgrey")+
# Montecristo
annotate(geom = "text", y = 14.23, x = -89.2, label = "MNP", size = 4)+
annotate(geom = "point", y = 14.33, x = -89.36, size = 5, shape = 21, fill = "lightgrey")+
# Nancuchiname
annotate(geom = "text", y = 13.45, x = -88.64, label = "NF", size = 4)+
annotate(geom = "point", y = 13.33, x = -88.72, size = 5, shape = 21, fill = "lightgrey")+
# Country & Ocean names
annotate(geom = "text", y = 13.8, x = -89.0, label = "EL SALVADOR", size = 3)+
annotate(geom = "text", y = 14.2, x = -90.0, label = "GUATEMALA", size = 3)+
annotate(geom = "text", y = 14.2, x = -87.8, label = "HONDURAS", size = 3)+
annotate(geom = "text", y = 13.3, x = -89.5, label = "Pacific Ocean", size = 2.5,
fontface = "italic")+
ylab("Latitude")+
xlab("Longitude")+
theme_bw()
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-02-28
##
## ─ Packages ──────────────────────────────────────────────────────────────────────────────────────────────
## package * version date lib source
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.2)
## callr 3.5.1 2020-10-13 [1] CRAN (R 4.0.2)
## cli 2.2.0 2020-11-20 [1] CRAN (R 4.0.2)
## codetools 0.2-16 2018-12-24 [1] CRAN (R 4.0.3)
## 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)
## 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)
## foreign 0.8-80 2020-05-24 [1] CRAN (R 4.0.3)
## fs 1.5.0 2020-07-31 [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)
## 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)
## htmltools 0.5.0 2020-06-16 [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)
## lifecycle 0.2.0 2020-03-06 [1] CRAN (R 4.0.2)
## magrittr 2.0.1 2020-11-17 [1] CRAN (R 4.0.2)
## mapproj * 1.2.7 2020-02-03 [1] CRAN (R 4.0.2)
## maps * 3.3.0 2018-04-03 [1] CRAN (R 4.0.2)
## maptools * 1.0-2 2020-08-24 [1] CRAN (R 4.0.2)
## markdown 1.1 2019-08-07 [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.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)
## plyr * 1.8.6 2020-03-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)
## R6 2.5.0 2020-10-28 [1] CRAN (R 4.0.2)
## raster * 3.4-5 2020-11-14 [1] CRAN (R 4.0.2)
## Rcpp 1.0.5 2020-07-06 [1] CRAN (R 4.0.2)
## remotes 2.2.0 2020-07-21 [1] CRAN (R 4.0.2)
## rgdal * 1.5-18 2020-10-13 [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)
## 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)
## sp * 1.4-4 2020-10-07 [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)
## 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)
## 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)
## 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)
## 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:
ezknitr::ezspin(file = “programs/g_graphing_results_global.R”, out_dir = “output”, fig_dir = “figures”, keep_md = F)