In USFWS/SppDistMonProj: Species Distribution Monitoring Project Data Workflow Tools

# Load packages
library(raster)
library(unmarked)
library(knitr)
library(rmarkdown)

knitr::opts_chunk$set(echo = FALSE, results = "asis", warning = TRUE)

dir.create("figure")
figdir.path <- "./figure/"

# Get user parameters
data <- params$data
year <- params$year

x <- load(data)
data <- get(x)
rm(x)

years <- data$years
nyears <- length(years)

nsites <- length(data$habcovData)

\newpage

Project description

The Species Occurrence Project is an annual refuge Inventory and Monitoring survey began in 2011 to monitor species occurrence across a National Wildlife Refuge. The survey consists of conducting multiple aerial snow track surveys for a species during winter months. The survey design as described in the Species Occurrence Project Protocol requires the surveys be conducted at each site every day for a 5-day period, each day (i.e., occasion) constituting an independent sample of occurrence data. Sites are equal-area polygons (i.e., grid cells) as defined by a lattice overlaid on the Refuge boundary (Fig. \@ref(fig:surveymap)) which divides the refuge into r nsites unique sites. Pilot observers are instructed to search each cell to which he/she is assigned with the same intensity to standardize detection by survey effort. Detection is known to be affected by visibility conditions and, therefore, a categorical visibility covariate (i.e., good vs poor) is collected for each survey occasion in each year to be used in data analysis. Habitat quality is a known factor affecting the probability of species occurrence. Because habitat quality across the the refuge does not substantially vary over time (i.e., year to year), a static spatial data set measuring habitat quality at the level of the survey site (Fig. \@ref(fig:habitatmap)) is used in the analysis to account for variation in occurrence due to habitat quality.

verts = as.vector(sp::bbox(data$habcovData))
rasterpts = data.frame(raster::rasterToPoints(data$habcovData))
raster::plot(raster::rasterToPolygons(data$habcovData), axes = TRUE, xlab = "X", ylab= " Y", asp = 1)
text(rasterpts$x, rasterpts$y, label = 1:nrow(rasterpts), cex = 0.75)
polygon(x = verts, y = verts[c(1,3,4,2)], lwd = 2, border = "blue")

rasterpts = data.frame(rasterToPoints(data$habcovData))
raster::plot(data$habcovData, xlab = "X", ylab = "Y", asp = 1)
text(rasterpts$x, rasterpts$y, label = 1:nrow(rasterpts), cex = 0.75)

\newpage

Analysis methods

The Species Occurrence Project began in 2011 and surveys have been conducted each year up to this year, r max(years). Four single-species occupancy models are fit to each year of data and model parameters are estimated for each model. The 4 models are a null model (Null) with constant occupancy probability ($\Psi$) and constant detection probability (\textit{p}), a detection-covariate model (DetCov) with constant $\Psi$ and \textit{p} varying by an occasion-specific covariate (i.e., good and poor visibility), a occupancy-covariate model (OccuCov) with $\Psi$ varying by a site-specific covariate (i.e., habitat quality measure) and constant \textit{p}, and a global model (Global) with the same covariates in the DetCov and OccuCov models. Year-specific model sets are ranked using AIC and model-averaged parameter estimates are derived.

\newpage

Analysis results

The survey area corresponds to the entire refuge and is divided into r nsites equal-area polygons, each constituting a single survey site with a unique site ID (Fig. \@ref(fig:surveymap)). Each site is characterized by a measure of habitat quality measured on the the real number scale (Fig. \@ref(fig:habitatmap)).

models <- c(Null = as.formula("~ 1 ~ 1"), DetCov = as.formula("~ viscov ~ 1"),
            OccuCov = as.formula("~ 1 ~ habcov"), Global = as.formula("~ viscov ~ habcov"))
nmods <- length(models)
modnames <- c("Null","DetCov","OccuCov","Global")
# fit all candidate models to each year
allfits <- lapply(data$unmarkedDataList,
                  function(data)fitList(fits = lapply(models, function(x, data)occu(x, data), data = data)))
modseltables <- lapply(allfits, function(x)as(modSel(x), "data.frame"))
modranks <- as.data.frame(lapply(modseltables, function(x)match(modnames, x$model)),
                          row.names = modnames, col.names = paste0("Y", years))
newdata1 <- data.frame(habcov = seq(-3,3,0.1) )
Epsi <- lapply(allfits, predict, type = "state", newdata = newdata1, appendData = TRUE)
newdata2 <- data.frame(viscov = factor(c("good","poor")))
Ep <- lapply(allfits, predict, type = "det", newdata = newdata2, appendData = TRUE)
topmods <- mapply(get.topmod, allfits, modranks)
n.sample.sites <- sapply(data$siteDataList, function(x)nrow(x))
PAO <- as.data.frame(do.call("rbind", mapply(calc.pao, topmod = topmods, n = n.sample.sites,
                                            MoreArgs <- list(stat = "mean"), SIMPLIFY = FALSE)))
PAO$Year <- years
class(PAO) <- c("pao", "data.frame")
betas <- lapply(topmods, get.betas)
# save(allfits, file = paste0("./output/raw_analysis/all_model_fits_", year, ".gzip"))

\captionsetup{justification = raggedright, singlelinecheck = false}

\blandscape

dummy <- mapply(annual.tables, modseltables, years, SIMPLIFY = FALSE)

\elandscape

pao.plot(PAO, years, modranks)

betas.plot(betas, years, modranks)

\newpage

Individual year predicted values

This section contains year-specific estimates ...

\newpage

\blandscape

dummy <- mapply(annual.estimates, Epsi, Ep, years, MoreArgs = list(figdir.path = figdir.path), SIMPLIFY = FALSE)

\elandscape

\newpage

Appendix

\newpage

R session info

sessionInfo()

---

title: | | Species Distribution Monitoring Project | Annual data analysis report for r max(years)

---

USFWS/SppDistMonProj documentation built on May 27, 2020, 9:44 a.m.