diff --git a/.Rbuildignore b/.Rbuildignore index cbc8813..c5b1d46 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -8,4 +8,5 @@ pre_vignette logo.png frescalo.exe ^docs$ -vignettes\sparta_vignette.Rmd \ No newline at end of file +vignettes\sparta_vignette.Rmd +misc/ \ No newline at end of file diff --git a/.gitignore b/.gitignore index 0bf2a9c..da8b20b 100644 --- a/.gitignore +++ b/.gitignore @@ -8,4 +8,6 @@ # etc frescalo.exe -sparta.Rproj \ No newline at end of file +sparta.Rproj + +misc/ \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION old mode 100755 new mode 100644 index 3e84b22..647aca7 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: sparta Type: Package Title: Trend Analysis for Unstructured Data -Version: 0.1.48 -Date: 2019-01-15 +Version: 0.1.49 +Date: 2019-03-07 Authors@R: c(person("Tom", "August", role = c("aut", "cre"), email = "tomaug@ceh.ac.uk"), person("Gary", "Powney", role = c("aut")), person("Charlie", "Outhwaite", role = c("aut")), @@ -26,7 +26,8 @@ Imports: gdata, R2jags, LearnBayes, - rmarkdown + rmarkdown, + boot Suggests: testthat (>= 0.8.0), roxygen2, diff --git a/NAMESPACE b/NAMESPACE index b945f49..5a816c3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -9,6 +9,7 @@ export(WSS) export(createWeights) export(dataDiagnostics) export(date2timeperiod) +export(detection_phenology) export(formatOccData) export(frescalo) export(getBugsData) @@ -24,6 +25,7 @@ export(occDetModel) export(occurrenceChange) export(recsOverTime) export(reportingRateModel) +export(simOccData) export(siteSelection) export(siteSelectionMinL) export(siteSelectionMinTP) @@ -35,6 +37,7 @@ import(ggplot2) import(lme4) import(rmarkdown) import(sp) +importFrom(boot,inv.logit) importFrom(dplyr,distinct) importFrom(dplyr,group_by) importFrom(dplyr,summarise) diff --git a/Occ_viz_html.html b/Occ_viz_html.html deleted file mode 100644 index 7cb0be1..0000000 --- a/Occ_viz_html.html +++ /dev/null @@ -1,269 +0,0 @@ - - - - - - - - - - - - - - -Occupancy model summary - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - -

This is a visualiation of an occupancy model produced using the r-package sparta. For more information of sparta visit https://github.com/biologicalrecordscentre/sparta

-
-

Table of contents

-
    -
  1. Basic Information
    2. -
  2. Species trend
    3. -
  3. Traceplots - Annnual occupancy estimates
    4. -
  4. Traceplots - Annnual detectability estimates
    5. -
  5. Traceplots - Other parameters
    6. -
  6. Rhat values
    7. -
-
-

Basic Information

-
## Species: a
-
## Year range: 1980 - 1999
-
## Iterations: 1000
-
## Chains: 3
-
## Burn in: 10
-
## Thinning: 2
-
## Number of sites: 50
-
## Number of visits: 4700
-
## Number of sites with records of a: 50
-
## Number of observations of a: 511
-
-
-

Species trend

-

-
-
-

Traceplots - Annnual occupancy estimates

-

-
-
-

Traceplots - Annnual detectability estimates

-

-
-
-

Traceplots - Other parameters

-

-
-
-

Rhat values

-

-
-
- - - - -
- - - - - - - - diff --git a/R/detection_phenology.R b/R/detection_phenology.R new file mode 100644 index 0000000..471fbe5 --- /dev/null +++ b/R/detection_phenology.R @@ -0,0 +1,84 @@ +#' 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 +#' +#' @details +#' Takes a object of \code{OccDet} fitted with the \code{jul_date} option +#' +#' Calculate the phenology of detection and produces a plot of detectability over time for the reference data type. +#' +#' @return Some numbers. +#' +#' @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. +#' @export + + +########################################### +# calculate some stats +# test +# correct for starting on 1 July +########################################### + + +detection_phenology <- function(model, spname=NULL, bins=12){ + + require(sparta) + require(reshape2) + require(plyr) + require(boot) + require(ggplot2) + + data <- model$BUGSoutput$sims.list + + if(!"beta1" %in% names(data)) + stop("no phenological effect was modelled!") + + # 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) + + 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] + })) + + 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, + mean_pDet = mean(pDet), + lower95CI = quantile(pDet, 0.025), + upper95CI = quantile(pDet, 0.975)) + + # now convert the jds back to their equivalent Julian Dates + 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) + + ylab("Detection probability") + + ggtitle(spname) + + theme_bw() + gp + + # calculate some simple stats + +} + diff --git a/R/errorChecks.r b/R/errorChecks.r index 111da03..5981273 100644 --- a/R/errorChecks.r +++ b/R/errorChecks.r @@ -1,6 +1,6 @@ #' @importFrom dplyr distinct -errorChecks <- function(taxa = NULL, site = NULL, survey = NULL, closure_period = NULL, time_period = NULL, +errorChecks <- function(taxa = NULL, site = NULL, survey = NULL, replicate = NULL, closure_period = NULL, time_period = NULL, startDate = NULL, endDate = NULL, Date = NULL, time_periodsDF = NULL, dist = NULL, sim = NULL, dist_sub = NULL, sim_sub = NULL, minSite = NULL, useIterations = NULL, @@ -17,6 +17,7 @@ errorChecks <- function(taxa = NULL, site = NULL, survey = NULL, closure_period site=site, survey=survey, closure_period=closure_period, + replicate=replicate, time_period=time_period, startDate=startDate, endDate=endDate) @@ -51,21 +52,30 @@ errorChecks <- function(taxa = NULL, site = NULL, survey = NULL, closure_period } + ###### Make sure there are no NAs + ### Checks for taxa ### if(!is.null(taxa)){ - # Make sure there are no NAs if(!all(!is.na(taxa))) stop('taxa must not contain NAs') } ### Checks for site ### if(!is.null(site)){ - # Make sure there are no NAs if(!all(!is.na(site))) stop('site must not contain NAs') } + ### Checks for closure period ### + if(!is.null(closure_period)){ + if(!all(!is.na(closure_period))) stop('closure_period must not contain NAs') + } + + ### Checks for replicate ### + if(!is.null(replicate)){ + if(!all(!is.na(replicate))) stop('replicate must not contain NAs') + } + ### Checks for time_period ### if(!is.null(time_period)){ - # Make sure there are no NAs if(!all(!is.na(time_period))) stop('time_period must not contain NAs') } diff --git a/R/formatOccData.r b/R/formatOccData.r index 3d480c4..ade150d 100644 --- a/R/formatOccData.r +++ b/R/formatOccData.r @@ -8,6 +8,7 @@ #' @param site A character vector of site names, as long as the number of observations. #' @param survey A vector as long as the number of observations. #' This must be a Date if either closure_period is not supplied or if includeJDay = \code{TRUE} +#' @param replicate An optional vector to identify replicate samples (visits) per survey. Need not be globally unique (e.g can be 1, 2, .. n within surveys) #' @param closure_period An optional vector of integers specifying the closure period. #' If \code{FALSE} then closure_period will be extracted as the year from the survey. #' @param includeJDay Logical. If \code{TRUE} a Julian day column is returned in the @@ -85,17 +86,23 @@ #' site = site, #' survey = survey, #' closure_period = closure_period) -#' } +#' +#' # format the unicorns data +#' formatted_data <- formatOccData(taxa = unicorns$CONCEPT, +#' survey = unicorns$Date, +#' site = unicorns$kmsq) +#'} +#' #' @export #' @importFrom dplyr distinct #' @importFrom dplyr summarise #' @importFrom dplyr group_by #' @importFrom reshape2 dcast -formatOccData <- function(taxa, site, survey, closure_period = NULL, includeJDay = FALSE){ +formatOccData <- function(taxa, site, survey, replicate = NULL, closure_period = NULL, includeJDay = FALSE){ # Do error checks - errorChecks(taxa = taxa, site = site, survey = survey, closure_period = closure_period) + errorChecks(taxa = taxa, site = site, survey = survey, replicate = replicate, closure_period = closure_period) # Additional error check whether survey is a date if(!'POSIXct' %in% class(survey) & !'Date' %in% class(survey)){ @@ -109,8 +116,10 @@ formatOccData <- function(taxa, site, survey, closure_period = NULL, includeJDay } } # otherwise survey can be anything (including a character vector) + if(is.null(replicate)) replicate <- rep(1, length(survey)) # all + # Create dataframe from vectors - taxa_data <- data.frame(taxa, site, survey) + taxa_data <- data.frame(taxa, site, survey, replicate) # now add the time period (TP). It's defined by closure_period if supplied # TP was formerly referred to as year @@ -122,9 +131,10 @@ formatOccData <- function(taxa, site, survey, closure_period = NULL, includeJDay taxa_data$TP <- lookup$TP[match(closure_period, lookup$cp)] #check that surveys are nested within closure_periods - temp <- taxa_data - temp <- group_by(temp, survey) - temp <- summarise(temp, ns = length(unique(temp$TP))) + temp <- taxa_data %>% + group_by(survey) %>% + group_by(replicate) %>% + summarise(ns = length(unique(TP))) if(any(temp$ns) > 1) { # return a warning but assume they know what they're doing warning(paste(as.numeric(table(temp$ns >1)[2]), 'survey identities appear in multiple closure periods')) @@ -139,7 +149,7 @@ formatOccData <- function(taxa, site, survey, closure_period = NULL, includeJDay taxa_data <- distinct(taxa_data) # add list length column - taxa_data$visit <- paste(taxa_data$site, taxa_data$survey, sep="") # create a factor which combines date and site + taxa_data$visit <- paste(taxa_data$site, taxa_data$survey, taxa_data$replicate, sep="") # create a factor which combines site, survey & replicate visit_Ls <- as.data.frame(table(taxa_data$visit)) names(visit_Ls) <- c("visit", "L") taxa_data <- merge(taxa_data, visit_Ls) diff --git a/R/getBugsData.R b/R/getBugsData.R index 2b52c02..c1c2af7 100644 --- a/R/getBugsData.R +++ b/R/getBugsData.R @@ -55,7 +55,9 @@ getBugsData <- function(bugs_data, modeltype, verbose = FALSE, jul_date = { if(verbose) cat('Adding bugs_data elements for Julian Date\n') - JulDate <- occDetData$Jul_date + + JulDate <- occDetData$Jul_date #- 182 # centering has already been done in formatOccData! + bugs_data <- c(bugs_data, JulDate = list(as.numeric(JulDate))) return(bugs_data) diff --git a/R/occDetFunc.r b/R/occDetFunc.r index ec51975..b08b83a 100644 --- a/R/occDetFunc.r +++ b/R/occDetFunc.r @@ -168,7 +168,7 @@ occDetFunc <- function (taxa_name, occDetdata, spp_vis, n_iterations = 5000, nyr seed = NULL, model.function = NULL, regional_codes = NULL, region_aggs = NULL, additional.parameters = NULL, additional.BUGS.elements = NULL, additional.init.values = NULL, - return_data = FALSE){ + return_data = TRUE){ J_success <- requireNamespace("R2jags", quietly = TRUE) @@ -205,7 +205,11 @@ occDetFunc <- function (taxa_name, occDetdata, spp_vis, n_iterations = 5000, nyr if(!taxa_name %in% colnames(spp_vis)) stop('taxa_name is not the name of a taxa in spp_vis') # Add the focal column (was the species recorded on the visit?). Use the spp_vis dataframe to extract this info + nrow1 <- nrow(occDetdata) occDetdata <- merge(occDetdata, spp_vis[,c("visit", taxa_name)]) + + if(nrow1 != nrow(occDetdata)) stop('some visits have been lost') + names(occDetdata)[names(occDetdata) == taxa_name] <- "focal" # If we are using regional codes do some checks @@ -257,8 +261,15 @@ occDetFunc <- function (taxa_name, occDetdata, spp_vis, n_iterations = 5000, nyr # look for missing years before time frame can be extended using max_year parameter years <- (max(occDetdata$TP) - min(occDetdata$TP))+1 - if(length(unique(occDetdata$TP)) != years) stop('It looks like you have years with no data. This will crash BUGS') - + if(length(unique(occDetdata$TP)) != years) { + # find out which years have no data + missing_yrs <- with(occDetdata, setdiff(min(TP):max(TP), unique(TP))) + if(length(missing_yrs) ==1) + error_msg <- paste0('There are no visits in year ', missing_yrs,'. This will crash BUGS') + else + error_msg <- paste0('There are ', length(missing_yrs),' years with no visits, including ', missing_yrs[1],'. This will crash BUGS') + stop(error_msg) + } # Record the max and min values of TP min_year <- min(occDetdata$TP) diff --git a/R/simOccData.R b/R/simOccData.R new file mode 100644 index 0000000..9b96de1 --- /dev/null +++ b/R/simOccData.R @@ -0,0 +1,99 @@ +#' Simulate for Occupancy detection models +#' +#' Simulates some data suitable for use in sparta +#' The user defines the parameters for the data generation +#' At present it works with just one species and generates the list length probabalistically +#' +#' @return A list, the first two elements of which ('spp_vis' & 'occDetData') mimic the output of occDetFunc. +#' The third element ('Z') is the presence-absence state variable and the fourth ('p') is the true probability of detection. +#' +#' @examples +#' \dontrun{ +#' +#' # set the sparta options +#'sparta_options <- c('ranwalk', # prior on occupancy is set by last year's posterior +#' 'jul_date', # use the Julian date as a covariate on the detection probability +#' 'catlistlength', # categorises the visits into three sets of 'qualities' +#' 'halfcauchy') # prior on the precisions +#' +#' # simulate some data +#'mydata <- simOccData(nvisit=200, nsite=10, nTP=5, psi=0.5, beta1=0.1, beta2=-2e-3) +#'with(mydata, plot(occDetdata$Jul_date, p)) +#' +#'# run the occupancy model model +#'out <- occDetFunc('mysp', mydata$occDetdata, mydata$spp_vis, n_iter = 1e4, +#' modeltype = sparta_options, return_data=TRUE) +#' +#'out$BUGSoutput +#'detection_phenology(out) +#' +#'qplot(data=melt(out$BUGSoutput$sims.array), geom='line', +#' x=Var1, col=factor(Var2), y=value) + +#' facet_wrap(~Var3, ncol=4, scales='free') +#' +#' } +#' @importFrom reshape2 melt +#' @importFrom boot inv.logit +#' @export + +simOccData <- function( + nsites = 20, + nvisits = 100, + nTP = 10, + psi = 0.5, + trend = -0.01, # this proportion of sites changes state each TP + mu.lp = -1, + tau.lp = 10, + beta1=0.1, beta2=-2e-3, + dtype2.p = 3, + dtype3.p = 10 + ){ + + #-------------------------- State variable + Z <- matrix(NA, nrow=nsites, ncol=nTP) + + # set the initial distribution + Z[,1] <- rbinom(n=nsites, size=1, prob=psi) + + for(t in 2:nTP){ + transition <- rbinom(n=nsites, size=1, prob=abs(trend)) + if(trend > 0) Zt <- apply(cbind(Z[,(t-1)], transition), 1, max) + else if(trend < 0) Zt <- apply(cbind(Z[,(t-1)], transition), 1, min) + Z[,t] <- Zt + } + + mZ <- melt(Z) + names(mZ) <- c('site', 'TP', "Z") + + #-------------------------- observations + + alpha.p <- rnorm(n = nTP, mean = mu.lp, sd = sqrt(1/tau.lp)) + + # visits are randomly allocated to sites and TPs + # there is a nonzero probability that some TPs will have no visits, which could be a problem + occDetdata <- data.frame(visit = 1:nvisits, + site = sample.int(n=nsites, size=nvisits, repl=TRUE), + L = sample(c(1,2,4), size=nvisits, repl=TRUE), + TP = sample.int(n=nTP, size=nvisits, repl=TRUE), + Jul_date = sample.int(n=365, size=nvisits, repl=TRUE) - 182 + ) + # probability of detection + p <- inv.logit( alpha.p[occDetdata$TP] + + beta1 * (occDetdata$Jul_date) + + beta2 * (occDetdata$Jul_date)^2 + + dtype2.p * (occDetdata$L %in% 2:3) + + dtype3.p * (occDetdata$L == 4) + ) + + # true occupancy + mZ <- merge(occDetdata, mZ, all.x=TRUE)$Z + # Detection history + Y <- rbinom(n=nvisits, size=1, p= p*mZ) + spp_vis <- data.frame(visit = 1: nvisits, mysp = Y) + + return(list(spp_vis = spp_vis, + occDetdata = occDetdata, + Z=Z, p=p)) + +} + diff --git a/man/createWeights.Rd b/man/createWeights.Rd index 1a7f8c1..a8b27aa 100644 --- a/man/createWeights.Rd +++ b/man/createWeights.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/createWeights.r +% Please edit documentation in R/createWeights.R \name{createWeights} \alias{createWeights} \title{Create frescalo weights file} diff --git a/man/dataDiagnostics.Rd b/man/dataDiagnostics.Rd index ea2af82..e7af8e7 100644 --- a/man/dataDiagnostics.Rd +++ b/man/dataDiagnostics.Rd @@ -4,8 +4,7 @@ \alias{dataDiagnostics} \title{Data Diagnostics} \usage{ -dataDiagnostics(taxa, site, time_period, plot = TRUE, - progress_bar = TRUE) +dataDiagnostics(taxa, site, time_period, plot = TRUE, progress_bar = TRUE) } \arguments{ \item{taxa}{A character vector of taxon names, as long as the number of observations.} diff --git a/man/detection_phenology.Rd b/man/detection_phenology.Rd new file mode 100644 index 0000000..13f753b --- /dev/null +++ b/man/detection_phenology.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/detection_phenology.R +\name{detection_phenology} +\alias{detection_phenology} +\title{Diagnostics for the detectability with respect to Julian Date} +\usage{ +detection_phenology(model, spname = NULL, bins = 12) +} +\arguments{ +\item{model}{a fitted sparta model of class \code{OccDet}.} + +\item{spname}{optional name of the species (used for plotting)} + +\item{bins}{number of points to estimate across the year. Defaults to 12} +} +\value{ +Some numbers. +} +\description{ +Creates a plot of detectability over the season and calculates some simple statistics +} +\details{ +Takes a object of \code{OccDet} fitted with the \code{jul_date} option + +Calculate the phenology of detection and produces a plot of detectability over time for the reference data type. +} +\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. +} diff --git a/man/formatOccData.Rd b/man/formatOccData.Rd index 9c25e9e..30dc0c8 100644 --- a/man/formatOccData.Rd +++ b/man/formatOccData.Rd @@ -4,7 +4,7 @@ \alias{formatOccData} \title{Format data for Occupancy detection models} \usage{ -formatOccData(taxa, site, survey, closure_period = NULL, +formatOccData(taxa, site, survey, replicate = NULL, closure_period = NULL, includeJDay = FALSE) } \arguments{ @@ -15,6 +15,8 @@ formatOccData(taxa, site, survey, closure_period = NULL, \item{survey}{A vector as long as the number of observations. This must be a Date if either closure_period is not supplied or if includeJDay = \code{TRUE}} +\item{replicate}{An optional vector to identify replicate samples (visits) per survey. Need not be globally unique (e.g can be 1, 2, .. n within surveys)} + \item{closure_period}{An optional vector of integers specifying the closure period. If \code{FALSE} then closure_period will be extracted as the year from the survey.} @@ -91,7 +93,13 @@ formatted_data <- formatOccData(taxa = taxa, site = site, survey = survey, closure_period = closure_period) + +# format the unicorns data +formatted_data <- formatOccData(taxa = unicorns$CONCEPT, + survey = unicorns$Date, + site = unicorns$kmsq) } + } \references{ Isaac, N.J.B., van Strien, A.J., August, T.A., de Zeeuw, M.P. and Roy, D.B. (2014). diff --git a/man/htmlSummary.Rd b/man/htmlSummary.Rd index 0fa183b..1330589 100644 --- a/man/htmlSummary.Rd +++ b/man/htmlSummary.Rd @@ -4,8 +4,8 @@ \alias{htmlSummary} \title{Create HTML Report} \usage{ -htmlSummary(occDet, open = TRUE, output_dir = getwd(), - output_file = NULL, ...) +htmlSummary(occDet, open = TRUE, output_dir = getwd(), output_file = NULL, + ...) } \arguments{ \item{occDet}{An object of class occDet} diff --git a/man/occDetFunc.Rd b/man/occDetFunc.Rd index 870b197..71b3b16 100644 --- a/man/occDetFunc.Rd +++ b/man/occDetFunc.Rd @@ -4,13 +4,13 @@ \alias{occDetFunc} \title{Occupancy detection Function} \usage{ -occDetFunc(taxa_name, occDetdata, spp_vis, n_iterations = 5000, - nyr = 2, burnin = 1500, thinning = 3, n_chains = 3, - write_results = TRUE, output_dir = getwd(), modeltype = "sparta", - max_year = NULL, seed = NULL, model.function = NULL, - regional_codes = NULL, region_aggs = NULL, - additional.parameters = NULL, additional.BUGS.elements = NULL, - additional.init.values = NULL, return_data = FALSE) +occDetFunc(taxa_name, occDetdata, spp_vis, n_iterations = 5000, nyr = 2, + burnin = 1500, thinning = 3, n_chains = 3, write_results = TRUE, + output_dir = getwd(), modeltype = "sparta", max_year = NULL, + seed = NULL, model.function = NULL, regional_codes = NULL, + region_aggs = NULL, additional.parameters = NULL, + additional.BUGS.elements = NULL, additional.init.values = NULL, + return_data = TRUE) } \arguments{ \item{taxa_name}{A character giving the name of the species to be modelled.} diff --git a/man/plot.occDet.Rd b/man/plot.occDet.Rd index 7e93f92..3310c4c 100644 --- a/man/plot.occDet.Rd +++ b/man/plot.occDet.Rd @@ -4,8 +4,7 @@ \alias{plot.occDet} \title{Plot occDet Objects} \usage{ -\method{plot}{occDet}(x, y = NULL, main = x$SPP_NAME, reg_agg = "", - ...) +\method{plot}{occDet}(x, y = NULL, main = x$SPP_NAME, reg_agg = "", ...) } \arguments{ \item{x}{An object of class occDet} diff --git a/man/plot_GIS.Rd b/man/plot_GIS.Rd index 2197b5c..e79a915 100644 --- a/man/plot_GIS.Rd +++ b/man/plot_GIS.Rd @@ -6,12 +6,11 @@ \usage{ plot_GIS(gis_data = NULL, main = "", xlab = "", ylab = "", xlim = NULL, ylim = NULL, show.axis = TRUE, show.grid = TRUE, - grid.div = 1, round.grid = FALSE, grid.col = "grey", - fill.col = NA, line.col = NULL, bg.col = "white", box.col = NA, - new.window = TRUE, no.margin = FALSE, set.margin = TRUE, - max.dimen = 13, cex.main = 1.2, cex.lab = 1, cex.axis = 0.8, - blank.plot = FALSE, plot.shape = TRUE, additions = FALSE, - return.dimen = TRUE) + grid.div = 1, round.grid = FALSE, grid.col = "grey", fill.col = NA, + line.col = NULL, bg.col = "white", box.col = NA, new.window = TRUE, + no.margin = FALSE, set.margin = TRUE, max.dimen = 13, cex.main = 1.2, + cex.lab = 1, cex.axis = 0.8, blank.plot = FALSE, plot.shape = TRUE, + additions = FALSE, return.dimen = TRUE) } \arguments{ \item{gis_data}{GIS object, or list of objects, to be plotted.} diff --git a/man/simOccData.Rd b/man/simOccData.Rd new file mode 100644 index 0000000..9a1f36c --- /dev/null +++ b/man/simOccData.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simOccData.R +\name{simOccData} +\alias{simOccData} +\title{Simulate for Occupancy detection models} +\usage{ +simOccData(nsites = 20, nvisits = 100, nTP = 10, psi = 0.5, + trend = -0.01, mu.lp = -1, tau.lp = 10, beta1 = 0.1, beta2 = -0.002, + dtype2.p = 3, dtype3.p = 10) +} +\value{ +A list, the first two elements of which ('spp_vis' & 'occDetData') mimic the output of occDetFunc. +The third element ('Z') is the presence-absence state variable and the fourth ('p') is the true probability of detection. +} +\description{ +Simulates some data suitable for use in sparta +The user defines the parameters for the data generation +At present it works with just one species and generates the list length probabalistically +} +\examples{ +\dontrun{ + +# set the sparta options +sparta_options <- c('ranwalk', # prior on occupancy is set by last year's posterior + 'jul_date', # use the Julian date as a covariate on the detection probability + 'catlistlength', # categorises the visits into three sets of 'qualities' + 'halfcauchy') # prior on the precisions + +# simulate some data +mydata <- simOccData(nvisit=200, nsite=10, nTP=5, psi=0.5, beta1=0.1, beta2=-2e-3) +with(mydata, plot(occDetdata$Jul_date, p)) + +# run the occupancy model model +out <- occDetFunc('mysp', mydata$occDetdata, mydata$spp_vis, n_iter = 1e4, + modeltype = sparta_options, return_data=TRUE) + +out$BUGSoutput +detection_phenology(out) + +qplot(data=melt(out$BUGSoutput$sims.array), geom='line', + x=Var1, col=factor(Var2), y=value) + + facet_wrap(~Var3, ncol=4, scales='free') + +} +} diff --git a/tests/testthat/testformatOccData.r b/tests/testthat/testformatOccData.r index 1c29f8e..c21834b 100644 --- a/tests/testthat/testformatOccData.r +++ b/tests/testthat/testformatOccData.r @@ -23,15 +23,21 @@ site <- sample(paste('A', 1:nSites, sep=''), size = n, TRUE) survey <- sample(rDates, size = n, TRUE) # set the closure period to be in 2 year bins -closure_period <- ceiling((as.numeric(format(rDates,'%Y')) - 2009)/2) +closure_period <- ceiling((as.numeric(format(survey,'%Y')) - 2009)/2) + +# create survey variable that is not a date and a closure period to match +survey_numbered <- as.integer(as.factor(survey)) + + +replicate <- rep(1, n) test_that("Test formatOccData", { expect_warning(visitData <- formatOccData(taxa = taxa, site = site, survey = survey), '854 out of 15000 observations will be removed as duplicates') - head_spp_vis <- structure(list(visit = c("A102010-02-17", "A102010-04-14", "A102010-04-22", - "A102010-08-29", "A102010-11-04", "A102011-02-09"), a = c(TRUE, + head_spp_vis <- structure(list(visit = c("A102010-02-171", "A102010-04-141", "A102010-04-221", + "A102010-08-291", "A102010-11-041", "A102011-02-091"), a = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE), b = c(FALSE, FALSE, FALSE, FALSE, FALSE, FALSE), c = c(FALSE, FALSE, FALSE, FALSE, TRUE, FALSE), d = c(FALSE, FALSE, FALSE, FALSE, TRUE, FALSE), e = c(FALSE, @@ -57,8 +63,8 @@ test_that("Test formatOccData", { 6L), class = "data.frame") - head_occDetdata <- structure(list(visit = c("A102010-02-17", "A102010-04-14", "A102010-04-22", - "A102010-08-29", "A102010-11-04", "A102011-02-09"), + head_occDetdata <- structure(list(visit = c("A102010-02-171", "A102010-04-141", "A102010-04-221", + "A102010-08-291", "A102010-11-041", "A102011-02-091"), site = structure(c(2L, 2L, 2L, 2L, 2L, 2L), .Label = c("A1", "A10", "A11", "A12", "A13", "A14", "A15", "A16", "A17", "A18", "A19", "A2", "A20", "A21", "A22", "A23", "A24", "A25", "A26", "A27", "A28", "A29", "A3", @@ -78,16 +84,47 @@ test_that("Test formatOccData", { test_that("Test formatOccData errors", { - expect_error(visitData <- formatOccData(taxa = head(taxa), site = site, survey = survey), - 'The following arguements are not of equal length: taxa, site, survey') - expect_error(visitData <- formatOccData(taxa = taxa, site = head(site), survey = survey, closure_period = closure_period), - 'The following arguements are not of equal length: taxa, site, survey, closure_period') - expect_error(visitData <- formatOccData(taxa = taxa, site = site, survey = head(survey), closure_period = closure_period), - 'The following arguements are not of equal length: taxa, site, survey, closure_period') - expect_error(visitData <- formatOccData(taxa = taxa, site = site, survey = survey, closure_period=head(closure_period)), - 'The following arguements are not of equal length: taxa, site, survey, closure_period') + expect_error(visitData <- formatOccData(taxa = head(taxa), site = site, survey = survey, closure_period = closure_period, replicate = replicate), + 'The following arguements are not of equal length: taxa, site, survey, closure_period, replicate') + expect_error(visitData <- formatOccData(taxa = taxa, site = head(site), survey = survey, closure_period = closure_period, replicate = replicate), + 'The following arguements are not of equal length: taxa, site, survey, closure_period, replicate') + expect_error(visitData <- formatOccData(taxa = taxa, site = site, survey = head(survey), closure_period = closure_period, replicate = replicate), + 'The following arguements are not of equal length: taxa, site, survey, closure_period, replicate') + expect_error(visitData <- formatOccData(taxa = taxa, site = site, survey = survey, closure_period=head(closure_period), replicate = replicate), + 'The following arguements are not of equal length: taxa, site, survey, closure_period, replicate') + expect_error(visitData <- formatOccData(taxa = taxa, site = site, survey = survey, closure_period=closure_period, replicate = head(replicate)), + 'The following arguements are not of equal length: taxa, site, survey, closure_period, replicate') }) +test_that("Test formatOccData date requirement errors", { + + expect_warning(visitData <- formatOccData(taxa = taxa, site = site, survey = survey, closure_period = closure_period), + '854 out of 15000 observations will be removed as duplicates') + expect_error(visitData <- formatOccData(taxa = taxa, site = site, survey = survey_numbered), + 'survey must be a date if closure_period not supplied') + expect_error(visitData <- formatOccData(taxa =taxa, site = site, survey = survey_numbered, includeJDay = TRUE, closure_period = closure_period), + 'survey must be a date if Julian Date is to be included') +}) - +test_that("Test formatOccData specified closure period", { + + expect_warning(visitData <- formatOccData(taxa = taxa, site = site, survey = survey, closure_period = closure_period), + '854 out of 15000 observations will be removed as duplicates') + + head_occDetdata_cp <- structure(list(visit = c("A102010-02-17", "A102010-04-14", "A102010-04-22", + "A102010-08-29", "A102010-11-04", "A102011-02-09"), + site = structure(c(2L, 2L, 2L, 2L, 2L, 2L), + .Label = c("A1", "A10", "A11", "A12", "A13", "A14", "A15", "A16", "A17", "A18", "A19", "A2", "A20", "A21", + "A22", "A23", "A24", "A25", "A26", "A27", "A28", "A29", "A3", + "A30", "A31", "A32", "A33", "A34", "A35", "A36", "A37", "A38", + "A39", "A4", "A40", "A41", "A42", "A43", "A44", "A45", "A46", + "A47", "A48", "A49", "A5", "A50", "A6", "A7", "A8", "A9"), class = "factor"), + L = c(5L, 2L, 1L, 5L, 5L, 1L), + TP = c(1L, 1L, 1L, 1L, 1L, 1L)), + .Names = c("visit", "site", "L", "TP"), + row.names = c(1L, 6L, 8L, 9L, 14L, 19L), class = "data.frame") + + expect_identical(head(visitData$occDetdata), head_occDetdata_cp) + +}) \ No newline at end of file diff --git a/tests/testthat/testfrescalo.r b/tests/testthat/testfrescalo.r index 086e693..0b25cc9 100644 --- a/tests/testthat/testfrescalo.r +++ b/tests/testthat/testfrescalo.r @@ -41,8 +41,8 @@ frespath <- file.path(tempdir(), 'fres.exe') test_that("Test errors", { - if (!capabilities('libcurl') | .Platform$OS.type != "windows") skip('skipping as libcurl not supported') - skip('Carbon black blocks Frescalo') + if (!capabilities('libcurl')) skip('skipping as libcurl not supported') + if (.Platform$OS.type == "windows") skip('Carbon black blocks Frescalo') if(.Platform$OS.type == "windows"){ download.file(url = 'https://github.com/BiologicalRecordsCentre/frescalo/raw/master/Frescalo_3a_windows.exe', @@ -54,6 +54,7 @@ test_that("Test errors", { destfile = frespath, method = "libcurl", quiet = TRUE) + system(command = paste('chmod', '+x', normalizePath(frespath))) } else{ stop(paste('frescalo is not supported on', .Platform$OS.type)) } @@ -100,10 +101,10 @@ test_that("Test errors", { test_that("Runs without error", { - - if (!capabilities('libcurl') | .Platform$OS.type != "windows") skip('skipping as libcurl not supported') - skip('Carbon black blocks Frescalo') + if (!capabilities('libcurl')) skip('skipping as libcurl not supported') + + if (.Platform$OS.type == "windows") skip('Carbon black blocks Frescalo') # This first run is done using years temp <- tempfile(pattern = 'dir') diff --git a/tests/testthat/testoccDetFunc.r b/tests/testthat/testoccDetFunc.r index 4ab093a..336b86f 100644 --- a/tests/testthat/testoccDetFunc.r +++ b/tests/testthat/testoccDetFunc.r @@ -33,13 +33,20 @@ suppressWarnings({visitData <- formatOccData(taxa = taxa, site = site, survey = survey)}) ## additional data for testing missing years -# remove one year +# remove one years time_period_missing <- sub("2018-", "2019-", survey) time_period_missing <- as.Date(time_period_missing) +# remove a second year +time_period_missing2 <- sub("2017-", "2019-", time_period_missing) +time_period_missing2 <- as.Date(time_period_missing2) + # format data suppressWarnings({visitData_missing <- formatOccData(taxa = taxa, site = site, survey = time_period_missing)}) +suppressWarnings({visitData_missing2 <- formatOccData(taxa = taxa, site = site, + survey = time_period_missing2)}) + test_that("Test occDetFunc errors", { @@ -69,8 +76,18 @@ test_that("Test occDetFunc errors", { spp_vis = visitData_missing$spp_vis, write_results = FALSE, seed = 111), - 'It looks like you have years with no data. This will crash BUGS') + 'There are no visits in year 2018. This will crash BUGS') + + expect_error(results <- occDetFunc(taxa_name = 'a', + n_iterations = 50, + burnin = 15, + occDetdata = visitData_missing2$occDetdata, + spp_vis = visitData_missing2$spp_vis, + write_results = FALSE, + seed = 111), + 'There are 2 years with no visits, including 2017. This will crash BUGS') + expect_error(results <- occDetFunc(taxa_name = 'a', n_iterations = 50, burnin = 15, @@ -112,7 +129,7 @@ test_that("Test occDetFunc with defaults", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) }) @@ -135,7 +152,7 @@ test_that("Test occDetFunc with model types", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) results <- occDetFunc(taxa_name = 'a', n_iterations = 50, @@ -150,7 +167,7 @@ test_that("Test occDetFunc with model types", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) results <- occDetFunc(taxa_name = 'a', n_iterations = 50, @@ -164,7 +181,7 @@ test_that("Test occDetFunc with model types", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) results <- occDetFunc(taxa_name = 'a', n_iterations = 50, @@ -179,7 +196,7 @@ test_that("Test occDetFunc with model types", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) results <- occDetFunc(taxa_name = 'a', @@ -195,7 +212,7 @@ test_that("Test occDetFunc with model types", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) results <- occDetFunc(taxa_name = 'a', n_iterations = 50, @@ -210,7 +227,7 @@ test_that("Test occDetFunc with model types", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) results <- occDetFunc(taxa_name = 'a', n_iterations = 50, @@ -225,7 +242,7 @@ test_that("Test occDetFunc with model types", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) results <- occDetFunc(taxa_name = 'a', n_iterations = 50, @@ -239,7 +256,7 @@ test_that("Test occDetFunc with model types", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) results <- occDetFunc(taxa_name = 'a', n_iterations = 50, @@ -254,7 +271,7 @@ test_that("Test occDetFunc with model types", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) results <- occDetFunc(taxa_name = 'a', @@ -270,7 +287,7 @@ test_that("Test occDetFunc with model types", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) results <- occDetFunc(taxa_name = 'a', n_iterations = 50, @@ -285,7 +302,7 @@ test_that("Test occDetFunc with model types", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) sink() @@ -314,7 +331,7 @@ test_that("Test occDetFunc with julian date", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) expect_true('beta1' %in% row.names(results$BUGSoutput$summary)) expect_true('beta2' %in% row.names(results$BUGSoutput$summary)) expect_true('beta3' %in% row.names(results$BUGSoutput$summary)) @@ -332,7 +349,7 @@ test_that("Test occDetFunc with julian date", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) expect_true('beta1' %in% row.names(results$BUGSoutput$summary)) expect_true('beta2' %in% row.names(results$BUGSoutput$summary)) expect_true('beta3' %in% row.names(results$BUGSoutput$summary)) @@ -364,7 +381,7 @@ test_that("Test occDetFunc with catagorical list length", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) expect_true('dtype2.p' %in% row.names(results$BUGSoutput$summary)) expect_true('dtype3.p' %in% row.names(results$BUGSoutput$summary)) @@ -381,7 +398,7 @@ test_that("Test occDetFunc with catagorical list length", { expect_identical(names(results), c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_sites", "species_observations")) + "nsites", "nvisits", "species_sites", "species_observations", "bugs_data")) expect_true('dtype2.p' %in% row.names(results$BUGSoutput$summary)) expect_true('dtype3.p' %in% row.names(results$BUGSoutput$summary)) @@ -412,7 +429,7 @@ test_that("Test occDetFunc using regions and region aggregates", { c("model", "BUGSoutput", "parameters.to.save", "model.file", "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", "nsites", "nvisits", "species_sites", "species_observations", - "regions", "region_aggs")) + "regions", "region_aggs", "bugs_data")) expect_identical(results$regions, c("region1", "region2", "region3")) expect_identical(names(results$region_aggs), "agg1")