Skip to content

Commit

Permalink
minor test fix to plot_DetectionPhenology
Browse files Browse the repository at this point in the history
  • Loading branch information
DylanCarbone committed Apr 3, 2024
1 parent 4c15cbd commit e3bc003
Show file tree
Hide file tree
Showing 3 changed files with 68 additions and 70 deletions.
87 changes: 44 additions & 43 deletions R/plot_DetectionPhenology.R
Original file line number Diff line number Diff line change
@@ -1,26 +1,26 @@
#' Diagnostics for the detectability with respect to Julian Date
#'
#' Diagnostics for the detectability with respect to Julian Date
#'
#' Creates a plot of detectability over the season and calculates some simple statistics
#'
#' @param model a fitted sparta model of class \code{OccDet}.
#' @param spname optional name of the species (used for plotting)
#' @param bins number of points to estimate across the year. Defaults to 12
#' @param density_function whether the model used a density function to fit Julian date. This form was implemented from version 0.1.48 onwards. For models ran using earlier versions of the package this should be set to FALSE
#'
#' @details
#'
#' @details
#' Takes a object of \code{OccDet} fitted with the \code{jul_date} option
#'
#'
#' Calculates the phenology of detection and produces a plot of detectability over time for the reference data type.
#'
#' @return This function returns plot showing the detection probability on the y axis and Julian day on the x.
#' The data within the output list shows the Julian day for each point estimated (equal to the number of bins)
#' and the mean detection probability with 95% credible intervals.
#'
#' @references van Strien, A.J., Termaat, T., Groenendijk, D., Mensing, V. & Kéry, M. (2010)
#' Site-occupancy models may offer new opportunities for dragonfly monitoring based on daily species lists.
#'
#' @references van Strien, A.J., Termaat, T., Groenendijk, D., Mensing, V. & Kéry, M. (2010)
#' Site-occupancy models may offer new opportunities for dragonfly monitoring based on daily species lists.
#' \emph{Basic and Applied Ecology}, 11, 495–503.
#'
#' @importFrom reshape2 melt
#' @importFrom reshape2 melt
#' @importFrom plyr ddply
#' @importFrom ggplot2 ggplot
#' @importFrom boot inv.logit
Expand All @@ -34,64 +34,65 @@
###########################################


plot_DetectionPhenology <- function(model, spname=NULL, bins=12, density_function = TRUE){

plot_DetectionPhenology <- function(model, spname = NULL, bins = 12, density_function = TRUE) {
data <- model$BUGSoutput$sims.list
if(!"beta1" %in% names(data))

if (!"beta1" %in% names(data)) {
stop("no phenological effect was modelled!")
}

# the base: alpha.p is common to all models:
# the base: alpha.p is common to all models:
# it's the logit probability of detection on a single species list
# use the average value across years
pDet1 <- rowMeans(data$alpha.p)
# pDet1 is an array of dims equal to (niter, nyr)

# Julian dates are 1:
jul_dates <- seq(from=1, to=365, length.out=bins)
if(density_function == TRUE){
if(!"beta3" %in% names(data))
jul_dates <- seq(from = 1, to = 365, length.out = bins)

if (density_function == TRUE) {
if (!"beta3" %in% names(data)) {
stop("beta3 not found. Please check that the density function method was used")
}
# we ran the Julian Date option
# So let's calculate the detection probabilities
pDet <- melt(sapply(jul_dates, function(jd){
pDet1 + data$beta3[,1] * (1/((2*pi)^0.5 * data$beta2[,1]) * exp(-((jd - data$beta1[,1])^2 / (2* data$beta2[,1]^2))))

pDet <- melt(sapply(jul_dates, function(jd) {
pDet1 + data$beta3[, 1] * (1 / ((2 * pi)^0.5 * data$beta2[, 1]) * exp(-((jd - data$beta1[, 1])^2 / (2 * data$beta2[, 1]^2))))
}))
}else{
} else {
cjd <- jul_dates - 182
# we ran the Julian Date option
# So let's scale the detection probabilities to end June (day 180)
pDet <- melt(sapply(cjd, function(jd){
pDet1 + jd * data$beta1[,1] + jd^2 * data$beta2[,1]
}))
pDet <- melt(sapply(cjd, function(jd) {
pDet1 + jd * data$beta1[, 1] + jd^2 * data$beta2[, 1]
}))
}
names(pDet) <- c("it", "bin","lgt_pDet")

names(pDet) <- c("it", "bin", "lgt_pDet")
pDet$pDet <- inv.logit(pDet$lgt_pDet)

# now summarize these posterior distributions
pDet_summary <-ddply(
pDet, .(bin), summarise,
pDet_summary <- ddply(
pDet, .(bin), summarise,
mean_pDet = mean(pDet),
lower95CI = quantile(pDet, 0.025),
upper95CI = quantile(pDet, 0.975))

# now convert the jds back to their equivalent Julian Dates
if(density_function == FALSE){pDet_summary$cJulDate <- cjd[pDet_summary$bin]}
pDet_summary$JulianDay <- jul_dates[pDet_summary$bin]
upper95CI = quantile(pDet, 0.975)
)

# now convert the jds back to their equivalent Julian Dates
if (density_function == FALSE) {
pDet_summary$cJulDate <- cjd[pDet_summary$bin]
}
pDet_summary$JulianDay <- jul_dates[pDet_summary$bin]

# now plot the detection over time
gp <- ggplot(data=pDet_summary, x=JulianDay, y=mean_pDet) +
geom_line(aes(x=JulianDay, y=mean_pDet)) +
geom_ribbon(aes(x=JulianDay, ymin=lower95CI, ymax=upper95CI), alpha=0.2) +
gp <- ggplot(data = pDet_summary, x = JulianDay, y = mean_pDet) +
geom_line(aes(x = JulianDay, y = mean_pDet)) +
geom_ribbon(aes(x = JulianDay, ymin = lower95CI, ymax = upper95CI), alpha = 0.2) +
ylab("Detection probability") +
ggtitle(spname) +
theme_bw()
gp

# calculate some simple stats

return(gp)
}

23 changes: 23 additions & 0 deletions tests/testthat/helper-plot_DetectionPhenology.r
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
# Create data
n <- 15000 # size of dataset
nyr <- 20 # number of years in data
nSamples <- 100 # set number of dates
nSites <- 50 # set number of sites
set.seed(125)

# Create somes dates
first <- as.Date(strptime("2010/01/01", format = "%Y/%m/%d"))
last <- as.Date(strptime(paste(2010 + (nyr - 1), "/12/31", sep = ""),
format = "%Y/%m/%d"
))
dt <- last - first
rDates <- first + (runif(nSamples) * dt)

# taxa are set as random letters
taxa <- sample(letters, size = n, TRUE)

# three sites are visited randomly
site <- sample(paste("A", 1:nSites, sep = ""), size = n, TRUE)

# the date of visit is selected at random from those created earlier
survey <- sample(rDates, size = n, TRUE)
Original file line number Diff line number Diff line change
@@ -1,28 +1,3 @@
context("Test plot_DetectionPhenology")

# Create data
n <- 15000 #size of dataset
nyr <- 20 # number of years in data
nSamples <- 100 # set number of dates
nSites <- 50 # set number of sites
set.seed(125)

# Create somes dates
first <- as.Date(strptime("2010/01/01", format = "%Y/%m/%d"))
last <- as.Date(strptime(paste(2010+(nyr-1), "/12/31", sep=''),
format = "%Y/%m/%d"))
dt <- last-first
rDates <- first + (runif(nSamples)*dt)

# taxa are set as random letters
taxa <- sample(letters, size = n, TRUE)

# three sites are visited randomly
site <- sample(paste('A', 1:nSites, sep = ''), size = n, TRUE)

# the date of visit is selected at random from those created earlier
survey <- sample(rDates, size = n, TRUE)

test_that("Test plot_DetectionPhenology", {

sink(file = ifelse(Sys.info()["sysname"] == "Windows",
Expand Down Expand Up @@ -79,8 +54,7 @@ test_that("Test plot_DetectionPhenology", {
head_data_JulianDay <- c(1, 34.0909090909091, 67.1818181818182, 100.272727272727, 133.363636363636,
166.454545454545)


expect_identical(names(results), c("data", "layers", "scales", "mapping", "theme", "coordinates", "facet", "plot_env", "labels"))
expect_identical(names(results), c("data", "layers", "scales", "guides", "mapping", "theme", "coordinates", "facet", "plot_env", "layout", "labels"))
expect_identical(names(results$data), c("bin", "mean_pDet", "lower95CI", "upper95CI", "JulianDay"))
expect_equal(head(results$data$bin), head_data_bin, tolerance = 1e-3)
expect_equal(head(results$data$mean_pDet), head_data_mean_pDet, tolerance = 1e-3)
Expand Down

0 comments on commit e3bc003

Please sign in to comment.