diff --git a/R/occDetFunc.r b/R/occDetFunc.r index f806a06..bdb26d4 100644 --- a/R/occDetFunc.r +++ b/R/occDetFunc.r @@ -1,11 +1,11 @@ #' Occupancy detection Function -#' +#' #' Run occupancy detection models using the output from \code{formatOccData} -#' +#' #' This function requires both the R package R2jags and the program JAGS. #' These are not installed by default when sparta is loaded and so should be #' installed by the user. More details can be found in teh vignette. -#' +#' #' @param taxa_name A character giving the name of the species to be modelled. #' @param occDetdata The 2nd element of the object returned by formatOccData. #' @param spp_vis The 1st element of the object returned by formatOccData. @@ -24,13 +24,13 @@ #' @param modeltype A character string or vector of strings that specifies the model to use. See details. If #' used then model.function is ignored. #' @param regional_codes A data.frame object detailing which site is associated with which region. -#' each row desginates a site and each column represents a region. The first column represents the +#' each row desginates a site and each column represents a region. The first column represents the #' site name (as in \code{site}). Subsequent columns are named for each regions with 1 representing #' the site is in that region and 0 that it is not. NOTE a site should only be in one region #' @param region_aggs A named list giving aggregations of regions that you want trend #' estimates for. For example \code{region_aggs = list(GB = c('england', 'scotland', 'wales'))} #' will produced a trend for GB (Great Britain) as well as its constituent nations. Note that -#' 'england', scotland' and 'wales' must appear as names of columns in \code{regional_codes}. +#' 'england', scotland' and 'wales' must appear as names of columns in \code{regional_codes}. #' More than one aggregate can be given, eg \code{region_aggs = list(GB = c('england', 'scotland', #' 'wales'), UK = c('england', 'scotland', 'wales', 'northern_ireland'))}. #' @param max_year numeric, final year to which analysis will be run, this can be set if it is beyond @@ -40,45 +40,45 @@ #' @param model.function optionally a user defined BUGS model coded as a function (see \code{?jags}, #' including the example there, for how this is done) #' @param additional.parameters A character vector of additional parameters to monitor -#' @param additional.BUGS.elements A named list giving additioanl bugs elements passed +#' @param additional.BUGS.elements A named list giving additioanl bugs elements passed #' to \code{R2jags::jags} 'data' argument -#' @param additional.init.values A named list giving user specified initial values to +#' @param additional.init.values A named list giving user specified initial values to #' be added to the defaults. #' @param return_data Logical, if \code{TRUE} (default) the BUGS data object is returned with the data #' @param saveMatrix Logical, if \code{FALSE} (default) the sims.matrix element of the jags object is omitted, in order to reduce the filesize. #' @param criterion Determines whether the model should be run. If an integer then this defines the threshold number of records (50 in Outhwaite et al 2019). -#' Other options are `EqualWt` or `HighSpec`, which define the application of "rules of thumb" defined in Pocock et al 2019. +#' Other options are `EqualWt` or `HighSpec`, which define the application of "rules of thumb" defined in Pocock et al 2019. #' Defaults to 1, in which case the model is applied for so long there is a single record of the focal species. #' @param provenance An optional text string allowing the user to identify the dataset. #' @param rem_aggs_with_missing_regions An option which if TRUE will remove all aggregates which contain at least one region with no data. #' If `FALSE`, only aggregates where ALL regions in that aggregate contain no data, are dropped. Defaults to TRUE #' @param allowSitesMultiRegions An option that permits sites to be included in more than one region if `TRUE`. If `FALSE` then these sites are dropped. Defaults to `FALSE` -#' +#' #' @details \code{modeltype} is used to choose the model as well as the associated initial values, #' and parameters to monitor. Elements to choose from can be separated into the following components: -#' +#' #' A. Prior type: this has 3 options, each of which was tested in Outhwaite et al (in review): #' 1. sparta - This uses the same model as in Isaac et al (2014). #' 2. indran - This is the adaptive stationary model. #' 3. ranwalk - This is the random walk model. -#' +#' #' B. Hyperprior type: This has 3 options, each of these are discussed in Outhwaite et al (in review): #' 1. halfuniform - the original formulation in Isaac et al (2014). #' 2. halfcauchy - preferred form, tested in Outhwaite et al (2018). #' 3. inversegamma - alternative form presented in the literature. -#' +#' #' C. List length specification: This has 3 options: #' 1. catlistlength - list length as a categorical variable. #' 2. contlistlength - list length as a continuous variable. #' 3. nolistlength - no list length variable. -#' +#' #' D. Julian date: this is an additional option for including Julian date within the detection model: #' 1. jul_date. -#' +#' #' Not all combinations are available in sparta. You will get an error if you try and use #' a combination that is not supported. There is usually a good reason why that #' combination is not a good idea. Here are the model elements available: -#' +#' #' \itemize{ #' \item{\code{"sparta"}}{ - This uses the same model as in Isaac et al (2014)} #' \item{\code{"indran"}}{ - Here the prior for the year effect of the state model is modelled as a random effect. This allows the model to adapt to interannual variability.} @@ -93,7 +93,7 @@ #' \item{\code{"centering"}}{ - No longer available. Includes hierarchical centering of the model parameters. Centring does not change the model explicitly but writes it in a way that allows parameter estimates to be updated simultaneously.} #' } #' These options are provided as a vector of characters, e.g. \code{modeltype = c('indran', 'halfcauchy', 'catlistlength')} -#' +#' #' @return A list including the model, JAGS model output, the path of the model file used and information on the number of iterations, first year, last year, etc. #' Key aspects of the model output include: #' \itemize{ @@ -131,33 +131,33 @@ #' Statistics for citizen science: extracting signals of change from noisy ecological data. #' \emph{Methods in Ecology and Evolution}, 5: 1052-1060. #' @references Outhwaite, C.L., Chandler, R.E., Powney, G.D., Collen, B., Gregory, R.D. & Isaac, N.J.B. (2018). -#' Prior specification in Bayesian occupancy modelling improves analysis of species occurrence data. +#' Prior specification in Bayesian occupancy modelling improves analysis of species occurrence data. #' \emph{Ecological Indicators}, 93: 333-343. #' @references Pocock, Logie, Isaac, Outhwaite & August. Rapid assessment of the suitability of multi-species citizen science datasets for occupancy trend analysis. \emph{bioRxiv} 813626 (2019) doi:10.1101/813626. -#' +#' #' @examples #' \dontrun{ -#' +#' #' set.seed(123) -#' +#' #' # Create data -#' n <- 15000 #size of dataset +#' n <- 15000 # size of dataset #' nyear <- 20 # number of years in data #' nSamples <- 100 # set number of dates #' nSites <- 50 # set number of sites -#' +#' #' # Create somes dates -#' first <- as.Date(strptime("2010/01/01", format="%Y/%m/%d")) -#' last <- as.Date(strptime(paste(2010+(nyear-1),"/12/31", sep=''), format="%Y/%m/%d")) -#' dt <- last-first -#' rDates <- first + (runif(nSamples)*dt) -#' +#' first <- as.Date(strptime("2010/01/01", format = "%Y/%m/%d")) +#' last <- as.Date(strptime(paste(2010 + (nyear - 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) -#' +#' #' # sites are visited randomly -#' site <- sample(paste('A', 1:nSites, sep=''), size = n, TRUE) -#' +#' 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) #' @@ -165,473 +165,522 @@ #' visitData <- formatOccData(taxa = taxa, site = site, survey = survey) #' #' # run the model with these data for one species (very small number of iterations) -#' results <- occDetFunc(taxa_name = taxa[1], -#' n_iterations = 50, -#' burnin = 15, -#' occDetdata = visitData$occDetdata, -#' spp_vis = visitData$spp_vis, -#' write_results = FALSE, -#' provenance = "sparta test dataset") +#' results <- occDetFunc( +#' taxa_name = taxa[1], +#' n_iterations = 50, +#' burnin = 15, +#' occDetdata = visitData$occDetdata, +#' spp_vis = visitData$spp_vis, +#' write_results = FALSE, +#' provenance = "sparta test dataset" +#' ) #' } #' @export #' @importFrom reshape2 acast #' @import LearnBayes -occDetFunc <- function (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, criterion = 1, provenance = NULL, saveMatrix = FALSE, - rem_aggs_with_missing_regions = TRUE, allowSitesMultiRegions = FALSE){ - +occDetFunc <- function(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, criterion = 1, provenance = NULL, saveMatrix = FALSE, + rem_aggs_with_missing_regions = TRUE, allowSitesMultiRegions = FALSE) { + ################## BASIC CHECKS # first run the error checks - errorChecks(n_iterations = n_iterations, burnin = burnin, - thinning = thinning, n_chains = n_chains, seed = seed) - + errorChecks( + n_iterations = n_iterations, burnin = burnin, + thinning = thinning, n_chains = n_chains, seed = seed + ) + # Set seed for repeatability - if(!is.null(seed)) set.seed(seed) - + if (!is.null(seed)) set.seed(seed) + # Check the taxa_name is one of my species - if(!taxa_name %in% colnames(spp_vis)) stop('taxa_name is not the name of a taxa in spp_vis') - ################## + if (!taxa_name %in% colnames(spp_vis)) stop("taxa_name is not the name of a taxa in spp_vis") + ################## min_year_original <- min_year <- min(occDetdata$TP) - + # only include sites which have more than nyr of records - yps <- rowSums(acast(occDetdata, site ~ TP, length, value.var = 'L') > 0) + yps <- rowSums(acast(occDetdata, site ~ TP, length, value.var = "L") > 0) sites_to_include <- names(yps[yps >= nyr]) - + # strip out the visits to sites that were visited in just one year i <- occDetdata$site %in% sites_to_include - - if(sum(i) > 0){ - occDetdata <- occDetdata[i,] - spp_vis <- spp_vis[i,] - } else stop(paste0("There are no sites visited in at least ", nyr, " years.")) - # calcluate a set of data metrics for this species - data_Metrics <- dataMetrics(sp = taxa_name, - formattedData = list(occDetdata=occDetdata, spp_vis=spp_vis)) + if (sum(i) > 0) { + occDetdata <- occDetdata[i, ] + spp_vis <- spp_vis[i, ] + } else { + stop(paste0("There are no sites visited in at least ", nyr, " years.")) + } + + # calculate a set of data metrics for this species + data_Metrics <- dataMetrics( + sp = taxa_name, + formattedData = list(occDetdata = occDetdata, spp_vis = spp_vis) + ) is.wholenumber <- - function(x, tol = .Machine$double.eps^0.5) { - if(is.numeric(x)) abs(x - round(x)) < tol - else FALSE + function(x, tol = .Machine$double.eps^0.5) { + if (is.numeric(x)) { + abs(x - round(x)) < tol + } else { + FALSE } + } # check there is enough data to run a model. If so, proceed with the main event - if(is.wholenumber(criterion)) { - # the criterion is a whole number. this defines the number of records - # check whether the number of records meets this value - proceed <- sum(spp_vis[,taxa_name]) >= criterion - } else if(criterion == "EqualWt") { - proceed <- applyRuleOfThumb(data_Metrics, "EqualWt") - } else if(criterion == "HighSpec") { - proceed <- applyRuleOfThumb(data_Metrics, "HighSpec") - } else - stop("Criterion must be either an integer, `EqualWt` or `HighSpec`") + if (is.wholenumber(criterion)) { + # the criterion is a whole number. this defines the number of records + # check whether the number of records meets this value + proceed <- sum(spp_vis[, taxa_name]) >= criterion + } else if (criterion == "EqualWt") { + proceed <- applyRuleOfThumb(data_Metrics, "EqualWt") + } else if (criterion == "HighSpec") { + proceed <- applyRuleOfThumb(data_Metrics, "HighSpec") + } else { + stop("Criterion must be either an integer, `EqualWt` or `HighSpec`") + } # 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)]) + occDetdata <- merge(occDetdata, spp_vis[, c("visit", taxa_name)]) names(occDetdata)[names(occDetdata) == taxa_name] <- "focal" - if(nrow1 != nrow(occDetdata)) stop('some visits have been lost') - - if(!proceed){ + if (nrow1 != nrow(occDetdata)) stop("some visits have been lost") + + if (!proceed) { # there is not enough data: set the outputs accordingly - bugs_data <- list(y=0,nsite=0,nvisit=0) + bugs_data <- list(y = 0, nsite = 0, nvisit = 0) BD_MD <- error_status <- site_match <- modelcode <- NA - warning(paste(taxa_name, - "has insufficient data after site filtering. Either decrease nyr or change the criterion")) - out <- list(message = "No model run: insufficient data") + warning(paste( + taxa_name, + "has insufficient data after site filtering. Either decrease nyr or change the criterion" + )) + out <- list(message = "No model run: insufficient data") } else { # There is enough data: we can proceed with the main event - + # Check if R2jags is installed if (!requireNamespace("R2jags", quietly = TRUE)) { stop("Package 'R2jags' is needed for the 'occDetModel' function to work. Please insatll this from CRAN. You will also be required to install JAGS, which you can download from https://sourceforge.net/projects/mcmc-jags/files/JAGS/", - call. = FALSE) + call. = FALSE + ) } - + # If doing regional we take control of model specification - if(!is.null(regional_codes)){ + if (!is.null(regional_codes)) { oldmodeltype <- modeltype - modeltype <-c('ranwalk', 'halfcauchy', - c('jul_date', 'catlistlength')[c('jul_date', 'catlistlength') %in% modeltype]) - if(!all(oldmodeltype %in% modeltype)) - message('When using regional data the model specification will be set to ranwalk, halfcauchy. jul_date and catlistlength can still be specified by the user') - } - - if(!'catlistlength' %in% modeltype & !'nolistlength' %in% modeltype){ - modeltype <- c(modeltype, 'contlistlength') - } - if(!any(c('catlistlength', 'nolistlength', 'contlistlength') %in% modeltype)){ + modeltype <- c( + "ranwalk", "halfcauchy", + c("jul_date", "catlistlength")[c("jul_date", "catlistlength") %in% modeltype] + ) + if (!all(oldmodeltype %in% modeltype)) { + message("When using regional data the model specification will be set to ranwalk, halfcauchy. jul_date and catlistlength can still be specified by the user") + } + } + + if (!"catlistlength" %in% modeltype & !"nolistlength" %in% modeltype) { + modeltype <- c(modeltype, "contlistlength") + } + if (!any(c("catlistlength", "nolistlength", "contlistlength") %in% modeltype)) { stop('modeltype should contain one of "catlistlength", "nolistlength", "contlistlength", which specify the list-length effect to be included') } - + # Do we have JAGS installed - this works only on windows - if(.Platform$OS.type == "windows"){ - JAGS_test <- Sys.which(names = 'jags-terminal.exe') - if(JAGS_test[[1]] == '') stop('R cannot find jags-terminal.exe, check that you have installed JAGS') + if (.Platform$OS.type == "windows") { + JAGS_test <- Sys.which(names = "jags-terminal.exe") + if (JAGS_test[[1]] == "") stop("R cannot find jags-terminal.exe, check that you have installed JAGS") } - + # If we are using regional codes do some checks - if(!is.null(regional_codes)){ - if(!inherits(regional_codes, 'data.frame')) stop("regional_codes should be a data.frame") - - # check whether there is a column called "site". - #If not, let's assume that the site column is the first in the dataframe - #NB previous behaviour was to assume *both* that it was column 1 and named 'site' - if(!"site" %in% names(regional_codes)) { + if (!is.null(regional_codes)) { + if (!inherits(regional_codes, "data.frame")) stop("regional_codes should be a data.frame") + + # check whether there is a column called "site". + # If not, let's assume that the site column is the first in the dataframe + # NB previous behaviour was to assume *both* that it was column 1 and named 'site' + if (!"site" %in% names(regional_codes)) { warning(paste0("renaming ", names(regional_codes)[1], " as 'site'")) - names(regional_codes)[1] <- "site" + names(regional_codes)[1] <- "site" } siteCol <- grep("site", names(regional_codes)) - + # remove locations that are not in the data abs_sites <- as.character(regional_codes$site)[!as.character(regional_codes$site) %in% as.character(occDetdata$site)] - if(length(abs_sites) > 0){ - warning(paste(length(abs_sites), 'sites are in regional_codes but not in occurrence data')) + if (length(abs_sites) > 0) { + warning(paste(length(abs_sites), "sites are in regional_codes but not in occurrence data")) } - if(any(is.na(regional_codes))){ - warning(paste(sum(is.na(regional_codes)), - "NAs are present in regional_codes, these will be replaced with 0s")) + if (any(is.na(regional_codes))) { + warning(paste( + sum(is.na(regional_codes)), + "NAs are present in regional_codes, these will be replaced with 0s" + )) regional_codes[is.na(regional_codes)] <- 0 } - - sites_no_region <- as.character(regional_codes$site[rowSums(regional_codes[,-siteCol]) == 0]) - sites_multi_region <- as.character(regional_codes$site[rowSums(regional_codes[,-siteCol]) > 1]) + + sites_no_region <- as.character(regional_codes$site[rowSums(regional_codes[, -siteCol]) == 0]) + sites_multi_region <- as.character(regional_codes$site[rowSums(regional_codes[, -siteCol]) > 1]) site_counts <- table(regional_codes$site) sites_multi_row <- names(site_counts[site_counts > 1]) - - if(length(sites_no_region) > 0) - warning(paste(length(sites_no_region), 'sites are not assigned to a region in regional_codes and will be removed')) - if(length(sites_multi_row) > 0) - warning(paste(length(sites_multi_row), 'sites appear more than once in regional_codes. These will be removed')) - if(length(sites_multi_region) > 0) - if(allowSitesMultiRegions) - warning(paste(length(sites_multi_region), 'sites are assigned to multiple region in regional_codes: you have specified this is allowed')) - else - warning(paste(length(sites_multi_region), 'sites are assigned to more than one region in regional_codes and will be removed')) - + + if (length(sites_no_region) > 0) { + warning(paste(length(sites_no_region), "sites are not assigned to a region in regional_codes and will be removed")) + } + if (length(sites_multi_row) > 0) { + warning(paste(length(sites_multi_row), "sites appear more than once in regional_codes. These will be removed")) + } + if (length(sites_multi_region) > 0) { + if (allowSitesMultiRegions) { + warning(paste(length(sites_multi_region), "sites are assigned to multiple region in regional_codes: you have specified this is allowed")) + } else { + warning(paste(length(sites_multi_region), "sites are assigned to more than one region in regional_codes and will be removed")) + } + } + # finally check that every site with species data also has a region sites_no_region2 <- setdiff(sites_to_include, as.character(regional_codes$site)) - if(length(sites_no_region2) >= 1) - warning(paste(length(sites_no_region2), 'sites are in occurrence data but not in regional_codes and will be removed')) - + if (length(sites_no_region2) >= 1) { + warning(paste(length(sites_no_region2), "sites are in occurrence data but not in regional_codes and will be removed")) + } + # strip these same sites out of the occDetdata & the regional codes bad_sites <- c(abs_sites, sites_multi_row, sites_no_region, sites_no_region2) - if(!allowSitesMultiRegions) bad_sites <- c(bad_sites, sites_multi_region) + if (!allowSitesMultiRegions) bad_sites <- c(bad_sites, sites_multi_region) bad_sites <- unique(bad_sites) regional_codes <- subset(regional_codes, !site %in% bad_sites) occDetdata <- subset(occDetdata, !site %in% bad_sites) } - + # If we are using regional aggregates do some checks - if(!is.null(region_aggs)){ - if(is.null(regional_codes)) stop('Cannot use regional aggregates if regional_codes is not supplied') - stopifnot(inherits(region_aggs, 'list')) - if(!all(unique(unlist(region_aggs)) %in% colnames(regional_codes)[-siteCol])){ - stop(paste0('Aggregate members [', - paste(unique(unlist(region_aggs))[!unique(unlist(region_aggs)) %in% colnames(regional_codes)[-siteCol]], - collapse = ', '), - '] not in regional_codes column names [', - paste(colnames(regional_codes)[-siteCol], - collapse = ', '), - ']')) + if (!is.null(region_aggs)) { + if (is.null(regional_codes)) stop("Cannot use regional aggregates if regional_codes is not supplied") + stopifnot(inherits(region_aggs, "list")) + if (!all(unique(unlist(region_aggs)) %in% colnames(regional_codes)[-siteCol])) { + stop(paste0( + "Aggregate members [", + paste(unique(unlist(region_aggs))[!unique(unlist(region_aggs)) %in% colnames(regional_codes)[-siteCol]], + collapse = ", " + ), + "] not in regional_codes column names [", + paste(colnames(regional_codes)[-siteCol], + collapse = ", " + ), + "]" + )) } } - + # 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) { + years <- (max(occDetdata$TP) - min(occDetdata$TP)) + 1 + 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 means there is no data, for any species in this year. BUGS cannot run if there is a year with no data. Quitting...') - else - error_msg <- paste0('There are ', length(missing_yrs),' years with no visits, including ', missing_yrs[1],'. This means there is no data, for any species in these years. BUGS cannot run if any year has no data. Quitting...') + 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 means there is no data, for any species in this year. BUGS cannot run if there is a year with no data. Quitting...") + } else { + error_msg <- paste0("There are ", length(missing_yrs), " years with no visits, including ", missing_yrs[1], ". This means there is no data, for any species in these years. BUGS cannot run if any year has no data. Quitting...") + } stop(error_msg) } # year and site need to be numeric starting from 1 to length of them. This is due to the way the bugs code is written nsite <- length(unique(occDetdata$site)) site_match <- data.frame(name = unique(occDetdata$site), id = 1:nsite) - occDetdata <- merge(occDetdata, site_match, by.x='site', by.y="name") - + occDetdata <- merge(occDetdata, site_match, by.x = "site", by.y = "name") + # need to get a measure of whether the species was on that site in that year, unequivocally, in zst - zst <- acast(occDetdata, id ~ factor(TP), value.var = 'focal', max, fill = 0) # initial values for the latent state = observed state + zst <- acast(occDetdata, id ~ factor(TP), value.var = "focal", max, fill = 0) # initial values for the latent state = observed state # Calculate min year. We're doing it now as it's fine if we've dropped the first year(s) of data and nothing in the middle min_year <- min(occDetdata$TP) - if(min_year != min_year_original){ - warning(paste0('\nThe first year of the data has changed, as the first years of data were dropped.\n', - 'The original first year was ',min_year_original,'. It is now ',min_year,'\n')) + if (min_year != min_year_original) { + warning(paste0( + "\nThe first year of the data has changed, as the first years of data were dropped.\n", + "The original first year was ", min_year_original, ". It is now ", min_year, "\n" + )) } - + # if the max_year is not null, edit the zst table to add the additional years required - if(!is.null(max_year)){ - + if (!is.null(max_year)) { # check that max_year is a numeric value - if(!is.numeric(max_year)) stop('max_year should be a numeric value') - + if (!is.numeric(max_year)) stop("max_year should be a numeric value") + # check that max_year is greater than the final year of the dataset - if(max_year <= max(occDetdata$TP)) stop('max_year should be greater than the final year of available data') - + if (max_year <= max(occDetdata$TP)) stop("max_year should be greater than the final year of available data") + nTP <- max_year - min_year + 1 - - # if nyear is greater than the number of years due to different specification of max_year, + + # if nyear is greater than the number of years due to different specification of max_year, # add on additional columns to zst so that inital values can be create for these years. - - if(nTP > ncol(zst)){ + + if (nTP > ncol(zst)) { # work out how many columns need to be added to_add <- nTP - ncol(zst) zst <- cbind(zst, matrix(0, ncol = to_add, nrow = nrow(zst))) # add column names - colnames(zst) <- 1:nTP + colnames(zst) <- 1:nTP } - + # if a value has not been selected for max_year then continue as before - }else{ + } else { # record the max year max_year <- max(occDetdata$TP) nTP <- max_year - min_year + 1 } - + # look for time periods with no data - if(length(unique(occDetdata$TP)) != nTP) stop('It looks like you have time periods with no data. This will crash BUGS') - + if (length(unique(occDetdata$TP)) != nTP) stop("It looks like you have time periods with no data. This will crash BUGS") + # TP and site need to be numeric starting from 1 to length of them. This is due to the way the bugs code is written occDetdata$TP <- occDetdata$TP - min(occDetdata$TP) + 1 - + # Parameter you wish to monitor, shown in the output parameters <- c("psi.fs", "tau2", "tau.lp", "alpha.p", "a") - - # If ranwalk + halfcauchy monitor mu.lp - if(all(c('ranwalk', 'halfcauchy') %in% modeltype)){ - if(!'centering' %in% tolower(modeltype) & !'intercept' %in% tolower(modeltype)){ + + # If ranwalk + halfcauchy monitor mu.lp + if (all(c("ranwalk", "halfcauchy") %in% modeltype)) { + if (!"centering" %in% tolower(modeltype) & !"intercept" %in% tolower(modeltype)) { parameters <- c(parameters, "mu.lp") } } - - # If sparta monitor mu.lp - if('sparta' %in% tolower(modeltype)) { + + # If sparta monitor mu.lp + if ("sparta" %in% tolower(modeltype)) { parameters <- c(parameters, "mu.lp") } - + # Add user specified paramters if given - if(!is.null(additional.parameters)) parameters <- c(parameters, additional.parameters) - + if (!is.null(additional.parameters)) parameters <- c(parameters, additional.parameters) + # Add parameters for each of the model types - for(ptype in modeltype){ + for (ptype in modeltype) { parameters <- getParameters(parameters, modeltype = ptype) } - + # add parameters for regions - if(!is.null(regional_codes)){ + if (!is.null(regional_codes)) { # remove spaces from region names, then extract them - colnames(regional_codes) <- gsub(' ', '_', colnames(regional_codes)) - region_names <- setdiff(colnames(regional_codes), "site") - - parameters <- c(parameters, - paste0("psi.fs.r_", region_names), - paste0("a_", region_names)) + colnames(regional_codes) <- gsub(" ", "_", colnames(regional_codes)) + region_names <- setdiff(colnames(regional_codes), "site") + + parameters <- c( + parameters, + paste0("psi.fs.r_", region_names), + paste0("a_", region_names) + ) # ignore some parameters not used in regions model - parameters <- parameters[!parameters %in% c('a')] + parameters <- parameters[!parameters %in% c("a")] } - + # add parameters for regional aggregates - if(!is.null(region_aggs)){ - parameters <- c(parameters, paste0('psi.fs.r_', names(region_aggs))) + if (!is.null(region_aggs)) { + parameters <- c(parameters, paste0("psi.fs.r_", names(region_aggs))) } - + # now assemble the bugs_data and related objects # HERE IS THE BUGS DATA - bugs_data <- with(occDetdata, - list(y = as.numeric(focal), Year = TP, Site = id, - nyear = nTP, nsite = nsite, nvisit = nrow(occDetdata))) + bugs_data <- with( + occDetdata, + list( + y = as.numeric(focal), Year = TP, Site = id, + nyear = nTP, nsite = nsite, nvisit = nrow(occDetdata) + ) + ) # temporary test - if(max(occDetdata$id) != bugs_data$nsite) stop(paste0("Site idenitifier exceeds nsite (", - max(occDetdata$id), nsite,")")) - - - for(btype in modeltype){ # one call per element of modeltype: each adds a section - bugs_data <- getBugsData(bugs_data, modeltype = btype, - occDetData = occDetdata) + if (max(occDetdata$id) != bugs_data$nsite) { + stop(paste0( + "Site idenitifier exceeds nsite (", + max(occDetdata$id), nsite, ")" + )) + } + + + for (btype in modeltype) { # one call per element of modeltype: each adds a section + bugs_data <- getBugsData(bugs_data, + modeltype = btype, + occDetData = occDetdata + ) } # Add additional elements if specified - if(!is.null(additional.BUGS.elements)){ - if(!is.list(additional.BUGS.elements)){ + if (!is.null(additional.BUGS.elements)) { + if (!is.list(additional.BUGS.elements)) { stop("additional.BUGS.elements must be a list") } else { bugs_data <- c(bugs_data, additional.BUGS.elements) } } - - # Add regional elements to bugs data - if(!is.null(regional_codes)){ + # Add regional elements to bugs data + if (!is.null(regional_codes)) { # removed unwanted bugs elements - bugs_data <- bugs_data[!names(bugs_data) %in% c('psi0.a', 'psi0.b')] - + bugs_data <- bugs_data[!names(bugs_data) %in% c("psi0.a", "psi0.b")] + # expand the lookup table to include regions - regional_lookup <- merge(regional_codes, site_match, by.y="name", by.x="site") - + regional_lookup <- merge(regional_codes, site_match, by.y = "name", by.x = "site") + zero_sites <- NULL - for(region in region_names){ - if(sum(regional_codes[ , region]) != 0){ - bugs_data[paste0('r_', region)] <- list(regional_lookup[order(regional_lookup$id),region]) - bugs_data[paste0('nsite_r_', region)] <- sum(regional_codes[, region]) + for (region in region_names) { + if (sum(regional_codes[, region]) != 0) { + bugs_data[paste0("r_", region)] <- list(regional_lookup[order(regional_lookup$id), region]) + bugs_data[paste0("nsite_r_", region)] <- sum(regional_codes[, region]) } else { zero_sites <- c(zero_sites, region) } } - - if(!is.null(zero_sites)){ - warning(paste('The following regions have no data and should not be modelled:', - paste(zero_sites, collapse = ', '), - '- These regions will not be included in the model')) + + if (!is.null(zero_sites)) { + warning(paste( + "The following regions have no data and should not be modelled:", + paste(zero_sites, collapse = ", "), + "- These regions will not be included in the model" + )) # remove parameters - parameters <- parameters[!parameters %in% c(paste0("psi.fs.r_", zero_sites), - paste0("a_", zero_sites))] + parameters <- parameters[!parameters %in% c( + paste0("psi.fs.r_", zero_sites), + paste0("a_", zero_sites) + )] # remove regions for regions_codes - regional_codes <- regional_codes[ ,!colnames(regional_codes) %in% zero_sites] + regional_codes <- regional_codes[, !colnames(regional_codes) %in% zero_sites] region_names <- setdiff(region_names, zero_sites) # remove region aggregates - if(rem_aggs_with_missing_regions){ + if (rem_aggs_with_missing_regions) { rem_aggs <- unlist(lapply(region_aggs, FUN = function(x) any(zero_sites %in% x))) rem_aggs_names <- names(region_aggs)[rem_aggs] - + # remove aggs if you need to - if(length(rem_aggs_names) > 0){ - warning(paste('The following region aggregates have to be removed as they contain a region with no data:', - paste(rem_aggs_names, collapse = ', '), - '- These region aggregates will not be included in the model\n', - 'If you want to keep aggregates with one or more missing regions,', - 'set rem_aggs_with_missing_regions=FALSE')) + if (length(rem_aggs_names) > 0) { + warning(paste( + "The following region aggregates have to be removed as they contain a region with no data:", + paste(rem_aggs_names, collapse = ", "), + "- These region aggregates will not be included in the model\n", + "If you want to keep aggregates with one or more missing regions,", + "set rem_aggs_with_missing_regions=FALSE" + )) region_aggs <- region_aggs[!names(region_aggs) %in% rem_aggs_names] - parameters <- parameters[!parameters %in% paste0('psi.fs.r_', rem_aggs_names)] + parameters <- parameters[!parameters %in% paste0("psi.fs.r_", rem_aggs_names)] } } else { rem_aggs <- unlist(lapply(region_aggs, FUN = function(x) all(x %in% zero_sites))) rem_aggs_names <- names(region_aggs)[rem_aggs] edit_aggs <- unlist(lapply(region_aggs, FUN = function(x) any(zero_sites %in% x))) edit_aggs_names <- names(region_aggs)[edit_aggs & !(rem_aggs)] - + # edit aggs if you need to - if(length(edit_aggs_names) > 0){ - warning(paste('The following region aggregates have to be edited as they contain regions with no data:', - paste(edit_aggs_names, collapse = ', '), - '\n- These region aggregates will still be included in the model\n')) + if (length(edit_aggs_names) > 0) { + warning(paste( + "The following region aggregates have to be edited as they contain regions with no data:", + paste(edit_aggs_names, collapse = ", "), + "\n- These region aggregates will still be included in the model\n" + )) # Recreate aggs, removing regions without data - region_aggs_new <- lapply(region_aggs, FUN = function(agg){ + region_aggs_new <- lapply(region_aggs, FUN = function(agg) { agg[!(agg %in% zero_sites)] }) names(region_aggs_new) <- names(region_aggs) region_aggs <- region_aggs_new } - + # remove aggs completely if you need to - if(length(rem_aggs_names) > 0){ - warning(paste('The following region aggregates have to be removed as all regions within them have no data:', - paste(rem_aggs_names, collapse = ', '), - '- These region aggregates will not be included in the model')) + if (length(rem_aggs_names) > 0) { + warning(paste( + "The following region aggregates have to be removed as all regions within them have no data:", + paste(rem_aggs_names, collapse = ", "), + "- These region aggregates will not be included in the model" + )) region_aggs <- region_aggs[!names(region_aggs) %in% rem_aggs_names] - parameters <- parameters[!parameters %in% paste0('psi.fs.r_', rem_aggs_names)] + parameters <- parameters[!parameters %in% paste0("psi.fs.r_", rem_aggs_names)] } } - } + } regions_years <- list() temp <- bugs_data[grepl("^r_", names(bugs_data))] regions_nobs <- sapply(temp, function(x) sum(x[bugs_data$Site] * bugs_data$y)) regions_sites <- unlist(bugs_data[grepl("nsite_r", names(bugs_data))]) - + # make a copy of the bugs_data to calculate metadata from bugs_data_copy <- with(occDetdata, data.frame(y = as.numeric(focal), year = TP, site = site)) bugs_data_copy <- merge(bugs_data_copy, regional_codes, all.x = TRUE) - + # add regional codes to this copy and get max and min years and year gaps for each region - for(region_name in region_names){ - - #regions_nobs[paste0('n_obs_','r_', region_name)] <- sum(bugs_data_copy$y * bugs_data_copy[,region_name]) - #regions_sites[paste0('n_sites_','r_', region_name)] <- sum(bugs_data_copy[,region_name]) - current_r <- bugs_data_copy$y * bugs_data_copy[,region_name] * bugs_data_copy$year - current_r <- subset(current_r,current_r !=0) - - if(length(current_r) > 2){ - - current_rmin <- (min_year-1) + min(current_r) - current_rmax <- (min_year-1) + max(current_r) - regions_years[paste0('min_year_data_','r_', region_name)] <- current_rmin - regions_years[paste0('max_year_data_','r_', region_name)] <- current_rmax + for (region_name in region_names) { + # regions_nobs[paste0('n_obs_','r_', region_name)] <- sum(bugs_data_copy$y * bugs_data_copy[,region_name]) + # regions_sites[paste0('n_sites_','r_', region_name)] <- sum(bugs_data_copy[,region_name]) + current_r <- bugs_data_copy$y * bugs_data_copy[, region_name] * bugs_data_copy$year + current_r <- subset(current_r, current_r != 0) + + if (length(current_r) > 2) { + current_rmin <- (min_year - 1) + min(current_r) + current_rmax <- (min_year - 1) + max(current_r) + regions_years[paste0("min_year_data_", "r_", region_name)] <- current_rmin + regions_years[paste0("max_year_data_", "r_", region_name)] <- current_rmax current_datagaps <- dataGaps(current_r, min_year, max_year, current_rmin, current_rmax) - regions_years[paste0('gap_start_','r_', region_name)] <- current_datagaps$gap_start - regions_years[paste0('gap_end_','r_', region_name)] <- current_datagaps$gap_end - regions_years[paste0('gap_middle_','r_', region_name)] <- current_datagaps$gap_middle - - } else if(length(current_r) == 1) { - - current_rmin <- (min_year-1) + min(current_r) - current_rmax <- (min_year-1) + max(current_r) - regions_years[paste0('min_year_data_','r_', region_name)] <- current_rmin - regions_years[paste0('max_year_data_','r_', region_name)] <- current_rmax + regions_years[paste0("gap_start_", "r_", region_name)] <- current_datagaps$gap_start + regions_years[paste0("gap_end_", "r_", region_name)] <- current_datagaps$gap_end + regions_years[paste0("gap_middle_", "r_", region_name)] <- current_datagaps$gap_middle + } else if (length(current_r) == 1) { + current_rmin <- (min_year - 1) + min(current_r) + current_rmax <- (min_year - 1) + max(current_r) + regions_years[paste0("min_year_data_", "r_", region_name)] <- current_rmin + regions_years[paste0("max_year_data_", "r_", region_name)] <- current_rmax current_datagaps <- dataGaps(current_r, min_year, max_year, current_rmin, current_rmax) - regions_years[paste0('gap_start_','r_', region_name)] <- current_datagaps$gap_start - regions_years[paste0('gap_end_','r_', region_name)] <- current_datagaps$gap_end - regions_years[paste0('gap_middle_','r_', region_name)] <- NA - - } else if(length(current_r) < 1){ - - regions_years[paste0('min_year_data_','r_', region_name)] <- NA - regions_years[paste0('max_year_data_','r_', region_name)] <- NA - regions_years[paste0('gap_start_','r_', region_name)] <- NA - regions_years[paste0('gap_end_','r_', region_name)] <- NA - regions_years[paste0('gap_middle_','r_', region_name)] <- NA - + regions_years[paste0("gap_start_", "r_", region_name)] <- current_datagaps$gap_start + regions_years[paste0("gap_end_", "r_", region_name)] <- current_datagaps$gap_end + regions_years[paste0("gap_middle_", "r_", region_name)] <- NA + } else if (length(current_r) < 1) { + regions_years[paste0("min_year_data_", "r_", region_name)] <- NA + regions_years[paste0("max_year_data_", "r_", region_name)] <- NA + regions_years[paste0("gap_start_", "r_", region_name)] <- NA + regions_years[paste0("gap_end_", "r_", region_name)] <- NA + regions_years[paste0("gap_middle_", "r_", region_name)] <- NA } } } - + # start preparing metadata BD_MD <- list() - + # add max and min data years for the whole dataset all_years_data <- bugs_data$y * bugs_data$Year - all_years_data <- subset(all_years_data, all_years_data !=0) - BD_MD$min_year_data <- (min_year-1) + min(all_years_data) - BD_MD$max_year_data <- (min_year-1) + max(all_years_data) - + all_years_data <- subset(all_years_data, all_years_data != 0) + BD_MD$min_year_data <- (min_year - 1) + min(all_years_data) + BD_MD$max_year_data <- (min_year - 1) + max(all_years_data) + # use these to find year gap data - BD_MD$yeargaps<-dataGaps(all_years_data, min_year, max_year, BD_MD$min_year_data, BD_MD$max_year_data) - - + BD_MD$yeargaps <- dataGaps(all_years_data, min_year, max_year, BD_MD$min_year_data, BD_MD$max_year_data) + + initiate <- function(z, nTP, bugs_data) { - init <- list (z = z, - alpha.p = rep(runif(1, -2, 2), - nTP), - a = rep(runif(1, -2, 2), nTP), - eta = rep(runif(1, -2, 2), bugs_data$nsite)) - + init <- list( + z = z, + alpha.p = rep( + runif(1, -2, 2), + nTP + ), + a = rep(runif(1, -2, 2), nTP), + eta = rep(runif(1, -2, 2), bugs_data$nsite) + ) + # add extra init values if needed - for(itype in modeltype){ + for (itype in modeltype) { init <- getInitValues(init, modeltype = itype) } - + # if ranwalk + centreing a -> aa - if(all(c('ranwalk', 'centering') %in% modeltype)){ - names(init)[names(init) == 'a'] <- 'aa' + if (all(c("ranwalk", "centering") %in% modeltype)) { + names(init)[names(init) == "a"] <- "aa" } - + # add user specified values - if(!is.null(additional.init.values)){ - if(!is.list(additional.init.values)){ + if (!is.null(additional.init.values)) { + if (!is.list(additional.init.values)) { stop("additional.BUGS.elements must be a list") } else { init <- c(init, additional.init.values) @@ -639,157 +688,169 @@ occDetFunc <- function (taxa_name, occDetdata, spp_vis, n_iterations = 5000, nyr } return(init) } - - # set the initial values... - init.vals <- replicate(n_chains, initiate(z = zst, - nTP = nTP, - bugs_data = bugs_data), - simplify = F) - + + # set the initial values... + init.vals <- replicate(n_chains, initiate( + z = zst, + nTP = nTP, + bugs_data = bugs_data + ), + simplify = F + ) + # modify initial values for regional model - if(!is.null(regional_codes)){ + if (!is.null(regional_codes)) { # remove initial values for a and psi - init.vals <- lapply(init.vals, FUN = function(x){ - x[!names(x) %in% c('psi0', 'a')] + init.vals <- lapply(init.vals, FUN = function(x) { + x[!names(x) %in% c("psi0", "a")] }) } - + # Select the correct model file - if(is.null(model.function)){ + if (is.null(model.function)) { model.file <- getModelFile(modeltype, - regional_codes = regional_codes, - region_aggs = region_aggs) + regional_codes = regional_codes, + region_aggs = region_aggs + ) } else { - cat('Using user model.function') + cat("Using user model.function") model.file <- model.function } - - modelcode <- paste(readLines(model.file), collapse = '\n') - + + modelcode <- paste(readLines(model.file), collapse = "\n") + ### REVIEW CODE - cat('#### PLEASE REVIEW THE BELOW ####\n\n') - cat('Your model settings:', paste(modeltype, collapse = ', ')) - cat('\n\nModel File:\n\n') + cat("#### PLEASE REVIEW THE BELOW ####\n\n") + cat("Your model settings:", paste(modeltype, collapse = ", ")) + cat("\n\nModel File:\n\n") cat(modelcode) - cat('\n\nbugs_data:\n\n') + cat("\n\nbugs_data:\n\n") cat(str(bugs_data)) - cat('\n\ninit.vals:\n\n') + cat("\n\ninit.vals:\n\n") cat(str(init.vals)) - cat('\n\nparameters:\n\n') + cat("\n\nparameters:\n\n") cat(parameters) ### - - if(identical(model.file, occDetBUGScode)){ - warning("You have selected a formulation with potentially informative priors that are subject to boundary effects (See Outhwaite et al 2018 for details). - This option is retained within sparta for backwards compatibility only: we strongly recommend that you do + + if (identical(model.file, occDetBUGScode)) { + warning("You have selected a formulation with potentially informative priors that are subject to boundary effects (See Outhwaite et al 2018 for details). + This option is retained within sparta for backwards compatibility only: we strongly recommend that you do not use this option for inferring changes in species' distributions") } - - error_status <- try( - out <- R2jags::jags(bugs_data, init.vals, parameters, model.file = model.file, - n.chains = n_chains, n.iter = n_iterations, n.thin = thinning, - n.burnin = burnin, DIC = TRUE) + + error_status <- try( + out <- R2jags::jags(bugs_data, init.vals, parameters, + model.file = model.file, + n.chains = n_chains, n.iter = n_iterations, n.thin = thinning, + n.burnin = burnin, DIC = TRUE + ) ) - + dir.create(path = output_dir, showWarnings = FALSE) # create the top results folder - if (class(error_status) == "try-error" ){ - - warning('JAGS returned an error when modelling', taxa_name, 'error:', error_status[1]) + if (class(error_status) == "try-error") { + warning("JAGS returned an error when modelling", taxa_name, "error:", error_status[1]) return(NULL) - - } - } # end of "if(proceed)" - - ########################################## Add metadata - - # calcluate number of site:year combinations with repeat visits (across the whole dataset) - temp <- as.data.frame(with(occDetdata, table(site, TP)))$Freq - prop_visits_repeated <- mean(temp[temp>0] > 1) - - if(is.null(attributes(occDetdata))){ - metadata <- list() - } else if('metadata' %in% names(attributes(occDetdata))){ - metadata <- attr(occDetdata, 'metadata') - } else { - metadata <- list() } - - # get the sessionInfo and coerce into a useable format - session.info <- sessionInfo() - - packages <- c(sapply(session.info$otherPkgs, function(x) x$Version), - sapply(session.info$loadedOnly, function(x) x$Version)) - - MD <- list(method = 'sparta::occDetFunc', - call = call <- match.call(), - date = Sys.Date(), - user = Sys.info()['user'], - summary = list(species = taxa_name, - n_sites = length(unique(occDetdata$site)), - n_years = length(unique(occDetdata$TP)), - n_visits = nrow(occDetdata), - n_obs = sum(occDetdata$focal), - n_species_sites = length(unique(subset(occDetdata, focal==TRUE)$site)), - min_year_model = min_year, - max_year_model = max_year), - gaps = ifelse(is.na(BD_MD), NA, list( - min_year_data = BD_MD$min_year_data, - max_year_data = BD_MD$max_year_data, - gap_start = BD_MD$yeargaps$gap_start, - gap_end = BD_MD$yeargaps$gap_end, - gap_middle = BD_MD$yeargaps$gap_middle)), - spp_Metrics = as.list(data_Metrics), - dataset_Metrics = list(# dataset properties - totalObservations = sum(occDetdata$L), - propRepeats = prop_visits_repeated), - provenance = provenance, - output_path = ifelse(test = write_results, - file.path(getwd(), output_dir, paste(taxa_name, ".rdata", sep = "")), - NA), - session_info = list(session.info[-c(7:8)], - packages) - ) - - # add regional elements if applicable - if(!is.null(regional_codes) & proceed){ - MD$summary$region_nobs <- regions_nobs - MD$summary$region_years <- unlist(regions_years) - MD$summary$region_nsite <- regions_sites - }else{ - MD$summary$region_nobs <- NA - MD$summary$region_years <- NA - MD$summary$region_nsite <- NA - } - - # If the data coming in is the result of analysis we want to - # append this data - name <- 'analysis' - i = 1 - while(name %in% names(metadata)){ - name <- paste0(name, i) - i = i + 1 - } - - metadata[name] <- list(MD) - attr(out, 'metadata') <- metadata - - if(!saveMatrix) out$BUGSoutput$sims.matrix <- NULL - - out$SPP_NAME <- taxa_name - out$min_year <- min_year - out$max_year <- max_year - out$sites_included <- ifelse(test = proceed, yes = site_match, no = NA) - out$nsites <- bugs_data$nsite - out$nvisits <- bugs_data$nvisit - out$species_observations <- sum(bugs_data$y) - out$sparta_version <- packages["sparta"] - if(!is.null(regional_codes)) out$regions <- region_names - if(!is.null(region_aggs)) out$region_aggs <- region_aggs - if(return_data) out$bugs_data <- bugs_data - attr(out, 'modeltype') <- modeltype - attr(out, 'modelcode') <- modelcode - class(out) <- 'occDet' - if(write_results) save(out, file = file.path(output_dir, paste(taxa_name, ".rdata", sep = ""))) - return(out) - } - \ No newline at end of file + } # end of "if(proceed)" + + ########################################## Add metadata + + # calcluate number of site:year combinations with repeat visits (across the whole dataset) + temp <- as.data.frame(with(occDetdata, table(site, TP)))$Freq + prop_visits_repeated <- mean(temp[temp > 0] > 1) + + if (is.null(attributes(occDetdata))) { + metadata <- list() + } else if ("metadata" %in% names(attributes(occDetdata))) { + metadata <- attr(occDetdata, "metadata") + } else { + metadata <- list() + } + + # get the sessionInfo and coerce into a useable format + session.info <- sessionInfo() + packages <- c( + sapply(session.info$basePkgs, packageVersion, simplify = FALSE), + sapply(session.info$otherPkgs, function(x) x$Version) + ) + + MD <- list( + method = "sparta::occDetFunc", + call = call <- match.call(), + date = Sys.Date(), + user = Sys.info()["user"], + summary = list( + species = taxa_name, + n_sites = length(unique(occDetdata$site)), + n_years = length(unique(occDetdata$TP)), + n_visits = nrow(occDetdata), + n_obs = sum(occDetdata$focal), + n_species_sites = length(unique(subset(occDetdata, focal == TRUE)$site)), + min_year_model = min_year, + max_year_model = max_year + ), + gaps = ifelse(is.na(BD_MD), NA, list( + min_year_data = BD_MD$min_year_data, + max_year_data = BD_MD$max_year_data, + gap_start = BD_MD$yeargaps$gap_start, + gap_end = BD_MD$yeargaps$gap_end, + gap_middle = BD_MD$yeargaps$gap_middle + )), + spp_Metrics = as.list(data_Metrics), + dataset_Metrics = list( # dataset properties + totalObservations = sum(occDetdata$L), + propRepeats = prop_visits_repeated + ), + provenance = provenance, + output_path = ifelse(test = write_results, + file.path(getwd(), output_dir, paste(taxa_name, ".rdata", sep = "")), + NA + ), + session_info = list( + session.info[1:2], + packages + ) + ) + + # add regional elements if applicable + if (!is.null(regional_codes) & proceed) { + MD$summary$region_nobs <- regions_nobs + MD$summary$region_years <- unlist(regions_years) + MD$summary$region_nsite <- regions_sites + } else { + MD$summary$region_nobs <- NA + MD$summary$region_years <- NA + MD$summary$region_nsite <- NA + } + + # If the data coming in is the result of analysis we want to + # append this data + name <- "analysis" + i <- 1 + while (name %in% names(metadata)) { + name <- paste0(name, i) + i <- i + 1 + } + + metadata[name] <- list(MD) + attr(out, "metadata") <- metadata + + if (!saveMatrix) out$BUGSoutput$sims.matrix <- NULL + + out$SPP_NAME <- taxa_name + out$min_year <- min_year + out$max_year <- max_year + out$sites_included <- ifelse(test = proceed, yes = site_match, no = NA) + out$nsites <- bugs_data$nsite + out$nvisits <- bugs_data$nvisit + out$species_observations <- sum(bugs_data$y) + out$sparta_version <- packages["sparta"] + if (!is.null(regional_codes)) out$regions <- region_names + if (!is.null(region_aggs)) out$region_aggs <- region_aggs + if (return_data) out$bugs_data <- bugs_data + attr(out, "modeltype") <- modeltype + attr(out, "modelcode") <- modelcode + class(out) <- "occDet" + if (write_results) save(out, file = file.path(output_dir, paste(taxa_name, ".rdata", sep = ""))) + return(out) +} diff --git a/tests/testthat/helper-errorChecks.R b/tests/testthat/helper-errorChecks.R new file mode 100644 index 0000000..8f99d65 --- /dev/null +++ b/tests/testthat/helper-errorChecks.R @@ -0,0 +1,49 @@ +# Create data +n <- 3000 # size of dataset +nyr <- 10 # number of years in data +nSamples <- 30 # set number of dates +nSites <- 15 # set number of sites + +# Create somes dates +first <- as.POSIXct(strptime("2010/01/01", "%Y/%m/%d")) +last <- as.POSIXct(strptime(paste(2010 + (nyr - 1), "/12/31", sep = ""), "%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 +time_period <- sample(rDates, size = n, TRUE) +time_period_missing <- sample(c(rDates, NA), size = n, TRUE) + +dist_sub <- rnorm(n, 10, 1) +sim_sub <- rnorm(n, 10, 1) + +dist_sub_fac <- as.factor(rnorm(n, 10, 1)) +sim_sub_fac <- as.factor(rnorm(n, 10, 1)) + + +dist_sub_chr <- as.character(rnorm(n, 10, 1)) +sim_sub_chr <- as.character(rnorm(n, 10, 1)) + +dist_sub_fac <- as.factor(rnorm(n, 10, 1)) +sim_sub_fac <- as.factor(rnorm(n, 10, 1)) + + +# combine this to a dataframe +df <- data.frame( + taxa = taxa, + site = site, + time_period = as.character(time_period), + time_period_missing = time_period_missing, + dist_sub = dist_sub, + sim_sub = sim_sub, + dist_sub_chr = dist_sub_chr, + sim_sub_chr = sim_sub_chr, + dist_sub_fac = dist_sub_fac, + sim_sub_fac = sim_sub_fac +) diff --git a/tests/testthat/test-errorChecks.R b/tests/testthat/test-errorChecks.R new file mode 100644 index 0000000..320b5fd --- /dev/null +++ b/tests/testthat/test-errorChecks.R @@ -0,0 +1,272 @@ +# Test error for startDate not in date format +test_that("Does it detect when startDate is not in a date format?", { + expect_error( + errorChecks(startDate = df$time_period), + 'startDate is not in a date format. This should be of class "Date" or "POSIXct"' + ) +}) + +# Test error for startDate containing NAs +test_that("Does it detect when startDate contains NAs?", { + expect_error( + errorChecks(startDate = df$time_period_missing), + "startDate must not contain NAs" + ) +}) + +# Test error for Date not being a data.frame or date vector +test_that("Does it detect when date is not a dataframe or date vector?", { + expect_error( + errorChecks(Date = df$time_period), + "Date must be a data.frame or date vector" + ) +}) + +# Test error for Date containing NAs +test_that("Does it detect when time_period_missing contains NAs?", { + expect_error( + errorChecks(Date = df$time_period_missing), + "Date must not contain NAs" + ) +}) + +# Test error for endDate not in date format +test_that("Does it detect when time_period is not in a date format?", { + expect_error( + errorChecks(endDate = df$time_period), + 'endDate is not in a date format. This should be of class "Date" or "POSIXct"' + ) +}) + +# Test error for endDate containing NAs +test_that("Does it detect when endDate contains NAs?", { + expect_error( + errorChecks(endDate = df$time_period_missing), + "endDate must not contain NAs" + ) +}) + +# Test error for dist_sub must be integer or numeric +test_that("Does it detect when dist_sub is not an integer or numeric value?", { + expect_error( + errorChecks(dist_sub = df$dist_sub_chr, sim_sub = df$sim_sub), + "dist_sub must be integer or numeric" + ) + expect_error( + errorChecks(dist_sub = df$dist_sub_fac, sim_sub = df$sim_sub), + "dist_sub must be integer or numeric" + ) +}) + +test_that("Does it detect when sim_sub is not an integer or numeric value?", { + expect_error( + errorChecks(dist_sub = df$dist_sub, sim_sub = df$sim_sub_fac), + "sim_sub must be integer or numeric" + ) + + expect_error( + errorChecks(dist_sub = df$dist_sub, sim_sub = df$sim_sub_chr), + "sim_sub must be integer or numeric" + ) +}) + +# Test error for useIterations must be logical +test_that("Does it detect when iterations are not logical?", { + useIterations_num <- 1 + useIterations_chr <- "TRUE" + + + expect_error( + errorChecks(useIterations = useIterations_chr), + "useIterations must be logical" + ) + expect_error( + errorChecks(useIterations = useIterations_num), + "useIterations must be logical" + ) +}) + +# Test error for iterations, family, thinning, and numeric parameters +test_that("Does it detect when iterations are not an integer?", { + expect_error( + errorChecks(iterations = "1000"), + "iterations must be numeric or integer" + ) +}) + +test_that("Does it detect when Family is not Binomial or Bernoulli?", { + expect_error( + errorChecks(family = "Poisson"), + "family must be either Binomial or Bernoulli" + ) +}) + +test_that("Does it detect when thinning is larger than the number of iterations?", { + expect_error( + errorChecks(n_iterations = 1000, burnin = 500, thinning = 1100, n_chains = 3), + "thinning must not be larger that the number of iteration (n_iterations)", + fixed = TRUE + ) +}) + +test_that("Does it detect when n_iterations is not an integer", { + expect_error( + errorChecks(n_iterations = "1000", burnin = 500, thinning = 5, n_chains = 3), + "n_iterations should be numeric" + ) +}) + +test_that("Does it detect when burnin is not an integer?", { + expect_error( + errorChecks(n_iterations = 1000, burnin = "500", thinning = 5, n_chains = 3), + "burnin should be numeric" + ) +}) + +test_that("Does it detect when thinning is not an integer?", { + expect_error( + errorChecks(n_iterations = 1000, burnin = 500, thinning = "5", n_chains = 3), + "thinning should be numeric" + ) +}) + +test_that("Does it detect when n_chains is not an integer?", { + expect_error( + errorChecks(n_iterations = 1000, burnin = 500, thinning = 5, n_chains = "3"), + "n_chains should be numeric" + ) +}) + +test_that("Does it detect when The seed is not an integer?", { + expect_error( + errorChecks(seed = "1"), + "seed muct be numeric" + ) +}) + +# Test for detecting missing year_col when start_col and end_col are NA +test_that("can it detect when year_col is missing with start_col and end_col as NA?", { + expect_error( + errorChecks(year_col = NA, start_col = NA, end_col = df$time_period[1]), + "year_col or start_col and end_col must be given" + ) +}) + +# Test for detecting missing end_col when year_col is NA +test_that("can it detect when end_col is missing with year_col as NA?", { + expect_error( + errorChecks(year_col = NA, start_col = df$time_period[1], end_col = NA), + "year_col or start_col and end_col must be given" + ) +}) + +# Test for detecting simultaneous use of year_col with start_col and end_col +test_that("can it detect the simultaneous use of year_col with start_col and end_col?", { + expect_error( + errorChecks(year_col = NA, start_col = df$time_period[1], end_col = df$time_period[1]), + "year_col cannot be used at the same time as start_col and end_col" + ) +}) + +# Test for detecting phi value below the permitted range +test_that("can it detect when phi is below the permitted range?", { + expect_error( + errorChecks(phi = 0.1), + "phi is outside permitted range of 0.50 to 0.95" + ) +}) + +# Test for detecting phi value above the permitted range +test_that("can it detect when phi is above the permitted range?", { + expect_error( + errorChecks(phi = 0.99), + "phi is outside permitted range of 0.50 to 0.95" + ) +}) + +# Test for detecting alpha value below the permitted range +test_that("can it detect when alpha is below the permitted range?", { + expect_error( + errorChecks(alpha = 0.05), + "alpha is outside permitted range of 0.08 to 0.50" + ) +}) + +# Test for detecting alpha value above the permitted range +test_that("can it detect when alpha is above the permitted range?", { + expect_error( + errorChecks(alpha = 0.99), + "alpha is outside permitted range of 0.08 to 0.50" + ) +}) + +# Test for detecting fres_site_filter not being a character vector +test_that("can it detect when fres_site_filter is not a character vector?", { + expect_error( + errorChecks(fres_site_filter = c(1, 5, 10, 12)), + "fres_site_filter must be a character vector" + ) +}) + +# Test for detecting non_benchmark_sp not being a character vector +test_that("can it detect when non_benchmark_sp is not a character vector?", { + expect_error( + errorChecks(non_benchmark_sp = c(1, 5, 10, 12)), + "non_benchmark_sp must be a character vector" + ) +}) + +# Test for detecting fres_site_filter not being a character vector +test_that("can it detect when fres_site_filter is not a character vector and is a single number?", { + expect_error( + errorChecks(fres_site_filter = 1), + "fres_site_filter must be a character vector" + ) +}) + +# Test for detecting fres_site_filter not being a character vector +test_that("can it detect when fres_site_filter is not a character vector and is a single number?", { + expect_error( + errorChecks(non_benchmark_sp = 1), + "non_benchmark_sp must be a character vector" + ) +}) + +# Test for detecting incorrect time_periods format +test_that("can it detect when time_periods format is incorrect?", { + expect_error( + errorChecks(time_periods = as.matrix(df$time_period)), + 'time_periods should be a data.frame. e.g. "data.frame(start=c(1980,1990),end=c(1989,1999))"', + fixed = TRUE + ) +}) + +# Test for detecting incorrect filepath for a .exe file +test_that("can it detect when the filepath is not to a '.exe' file?", { + temp <- tempfile() + + dir.create(temp) + + + expect_error( + errorChecks(frespath = temp), + "filepath is not the path to a '.exe' file" + ) +}) + +# Test for detecting non-existent .exe file +test_that("can it detect when a .exe file does not exist?", { + expect_error( + errorChecks(frespath = "file.exe"), + "file.exe does not exist" + ) +}) + +# Test error for site must not contain empty values +test_that("Does it detect empty string values?", { + expect_error( + errorChecks(site = c("a", "b", "")), + "site must not contain empty values (i.e. '')", + fixed = TRUE + ) +}) diff --git a/tests/testthat/testoccDetFunc.r b/tests/testthat/testoccDetFunc.r index ff0a200..228e011 100755 --- a/tests/testthat/testoccDetFunc.r +++ b/tests/testthat/testoccDetFunc.r @@ -1,39 +1,45 @@ context("Test occDetFunc") # Create data -n <- 15000 #size of dataset +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", "%Y/%m/%d")) -last <- as.Date(strptime(paste(2010+(nyr-1),"/12/31", sep=''), "%Y/%m/%d")) -dt <- last-first -rDates <- first + (runif(nSamples)*dt) +first <- as.Date(strptime("2010/01/01", "%Y/%m/%d")) +last <- as.Date(strptime(paste(2010 + (nyr - 1), "/12/31", sep = ""), "%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) +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) # set up regions for some testing -regionsDF <- data.frame(site = unique(site), - region1 = c(rep(1, 20), rep(0, 30)), - region2 = c(rep(0, 20), rep(1, 15), rep(0, 15)), - region3 = c(rep(0, 20), rep(0, 15), rep(1, 15))) +regionsDF <- data.frame( + site = unique(site), + region1 = c(rep(1, 20), rep(0, 30)), + region2 = c(rep(0, 20), rep(1, 15), rep(0, 15)), + region3 = c(rep(0, 20), rep(0, 15), rep(1, 15)) +) # format data -suppressWarnings({visitData <- formatOccData(taxa = taxa, site = site, - survey = survey)}) +suppressWarnings({ + visitData <- formatOccData( + taxa = taxa, site = site, + survey = survey + ) +}) ## additional data for testing missing years -# remove one years +# remove one years time_period_missing <- sub("2018-", "2019-", survey) time_period_missing <- as.Date(time_period_missing) @@ -42,415 +48,617 @@ 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)}) +suppressWarnings({ + visitData_missing <- formatOccData( + taxa = taxa, site = site, + survey = time_period_missing + ) +}) +suppressWarnings({ + visitData_missing2 <- formatOccData( + taxa = taxa, site = site, + survey = time_period_missing2 + ) +}) # additional data for testing species lost after nyr filter -taxa_filterError <- c(taxa,"aa") -site_filterError <- c(site,"B1") -survey_filterError <- c(survey,sample(rDates, size = 1, TRUE)) +taxa_filterError <- c(taxa, "aa") +site_filterError <- c(site, "B1") +survey_filterError <- c(survey, sample(rDates, size = 1, TRUE)) -suppressWarnings({visitData_filterError <- formatOccData(taxa = taxa_filterError, site = site_filterError, - survey = survey_filterError)}) +suppressWarnings({ + visitData_filterError <- formatOccData( + taxa = taxa_filterError, site = site_filterError, + survey = survey_filterError + ) +}) test_that("Test occDetFunc errors", { - - expect_error(results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 500, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111), - 'must not be larger that the number of iteration') - - expect_error(results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - criterion = "OnlyGoodData", - seed = 111), - "Criterion must be either an integer, `EqualWt` or `HighSpec`") - - - expect_error(results <- occDetFunc(taxa_name = 'apple', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - criterion = "EqualWt", - seed = 111), - 'taxa_name is not the name of a taxa in spp_vis') - - expect_error(results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData_missing$occDetdata, - spp_vis = visitData_missing$spp_vis, - write_results = FALSE, - criterion = "HighSpec", - seed = 111), - 'There are no visits in year 2018. This means there is no data, for any species in this year. BUGS cannot run if there is a year with no data. Quitting...') - - 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 means there is no data, for any species in these years. BUGS cannot run if any year has no data. Quitting...') - - - expect_error(results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - max_year = "2028", - write_results = FALSE, - seed = 111), - 'max_year should be a numeric value') - - expect_error(results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - max_year = 2028, - write_results = FALSE, - seed = 111), - 'max_year should be greater than the final year of available data') - + # skip the test if jagsUI is not installed + if ((!requireNamespace("jagsUI", quietly = TRUE))) { + skip("jagsUI software not available") + } + + expect_error( + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 500, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111 + ), + "must not be larger that the number of iteration" + ) + + expect_error( + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + criterion = "OnlyGoodData", + seed = 111 + ), + "Criterion must be either an integer, `EqualWt` or `HighSpec`" + ) + + + expect_error( + results <- occDetFunc( + taxa_name = "apple", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + criterion = "EqualWt", + seed = 111 + ), + "taxa_name is not the name of a taxa in spp_vis" + ) + + expect_error( + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData_missing$occDetdata, + spp_vis = visitData_missing$spp_vis, + write_results = FALSE, + criterion = "HighSpec", + seed = 111 + ), + "There are no visits in year 2018. This means there is no data, for any species in this year. BUGS cannot run if there is a year with no data. Quitting..." + ) + + 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 means there is no data, for any species in these years. BUGS cannot run if any year has no data. Quitting..." + ) + + + expect_error( + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + max_year = "2028", + write_results = FALSE, + seed = 111 + ), + "max_year should be a numeric value" + ) + + expect_error( + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + max_year = 2028, + write_results = FALSE, + seed = 111 + ), + "max_year should be greater than the final year of available data" + ) }) test_that("Test occDetFunc with defaults", { - - sink(file=ifelse(Sys.info()["sysname"] == "Windows", - "NUL", - "/dev/null")) - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111) + # skip the test if jagsUI is not installed + if ((!requireNamespace("jagsUI", quietly = TRUE))) { + skip("jagsUI software not available") + } + + sink(file = ifelse(Sys.info()["sysname"] == "Windows", + "NUL", + "/dev/null" + )) + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111 + ) sink() - expect_identical(results$SPP_NAME, 'a') + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations","sparta_version")) - + 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_observations", "sparta_version" + ) + ) }) test_that("Test occDetFunc with model types", { - - sink(file=ifelse(Sys.info()["sysname"] == "Windows", - "NUL", - "/dev/null")) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('indran', 'centering', - 'halfcauchy')) - expect_identical(results$SPP_NAME, 'a') + # skip the test if jagsUI is not installed + if ((!requireNamespace("jagsUI", quietly = TRUE))) { + skip("jagsUI software not available") + } + + sink(file = ifelse(Sys.info()["sysname"] == "Windows", + "NUL", + "/dev/null" + )) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c( + "indran", "centering", + "halfcauchy" + ) + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations","sparta_version")) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('indran', 'centering', - 'inversegamma')) - expect_identical(results$SPP_NAME, 'a') + 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_observations", "sparta_version" + ) + ) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c( + "indran", "centering", + "inversegamma" + ) + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('indran', 'intercept')) - expect_identical(results$SPP_NAME, 'a') + 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_observations", "sparta_version" + ) + ) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c("indran", "intercept") + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('indran', 'intercept', - 'halfcauchy')) - expect_identical(results$SPP_NAME, 'a') + 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_observations", "sparta_version" + ) + ) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c( + "indran", "intercept", + "halfcauchy" + ) + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('indran', 'intercept', - 'inversegamma')) - expect_identical(results$SPP_NAME, 'a') + 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_observations", "sparta_version" + ) + ) + + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c( + "indran", "intercept", + "inversegamma" + ) + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('ranwalk', 'centering', - 'halfcauchy')) - expect_identical(results$SPP_NAME, 'a') + 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_observations", "sparta_version" + ) + ) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c( + "ranwalk", "centering", + "halfcauchy" + ) + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('ranwalk', 'centering', - 'inversegamma')) - expect_identical(results$SPP_NAME, 'a') + 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_observations", "sparta_version" + ) + ) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c( + "ranwalk", "centering", + "inversegamma" + ) + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('ranwalk', 'intercept')) - expect_identical(results$SPP_NAME, 'a') + 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_observations", "sparta_version" + ) + ) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c("ranwalk", "intercept") + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('ranwalk', 'intercept', - 'halfcauchy')) - expect_identical(results$SPP_NAME, 'a') + 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_observations", "sparta_version" + ) + ) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c( + "ranwalk", "intercept", + "halfcauchy" + ) + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('ranwalk', 'intercept', - 'inversegamma')) - expect_identical(results$SPP_NAME, 'a') + 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_observations", "sparta_version" + ) + ) + + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c( + "ranwalk", "intercept", + "inversegamma" + ) + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('ranwalk', - 'halfcauchy')) - expect_identical(results$SPP_NAME, 'a') + 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_observations", "sparta_version" + ) + ) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c( + "ranwalk", + "halfcauchy" + ) + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - + 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_observations", "sparta_version" + ) + ) + sink() - }) test_that("Test occDetFunc with julian date", { - - sink(file=ifelse(Sys.info()["sysname"] == "Windows", - "NUL", - "/dev/null")) - - suppressWarnings({visitData <- formatOccData(taxa = taxa, site = site, - survey = survey, - includeJDay = TRUE)}) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('ranwalk', 'intercept', - 'inversegamma', 'jul_date')) - expect_identical(results$SPP_NAME, 'a') + # skip the test if jagsUI is not installed + if ((!requireNamespace("jagsUI", quietly = TRUE))) { + skip("jagsUI software not available") + } + + sink(file = ifelse(Sys.info()["sysname"] == "Windows", + "NUL", + "/dev/null" + )) + + suppressWarnings({ + visitData <- formatOccData( + taxa = taxa, site = site, + survey = survey, + includeJDay = TRUE + ) + }) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c( + "ranwalk", "intercept", + "inversegamma", "jul_date" + ) + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - 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)) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('indran', 'centering', - 'halfcauchy', 'jul_date')) - expect_identical(results$SPP_NAME, 'a') + 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_observations", "sparta_version" + ) + ) + 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)) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c( + "indran", "centering", + "halfcauchy", "jul_date" + ) + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - 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)) - + 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_observations", "sparta_version" + ) + ) + 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)) + sink() - -}) +}) test_that("Test occDetFunc with catagorical list length", { - - sink(file=ifelse(Sys.info()["sysname"] == "Windows", - "NUL", - "/dev/null")) - - suppressWarnings({visitData <- formatOccData(taxa = taxa, site = site, - survey = survey, - includeJDay = TRUE)}) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('ranwalk', 'intercept', - 'inversegamma', 'catlistlength')) - expect_identical(results$SPP_NAME, 'a') + # skip the test if jagsUI is not installed + if ((!requireNamespace("jagsUI", quietly = TRUE))) { + skip("jagsUI software not available") + } + + sink(file = ifelse(Sys.info()["sysname"] == "Windows", + "NUL", + "/dev/null" + )) + + suppressWarnings({ + visitData <- formatOccData( + taxa = taxa, site = site, + survey = survey, + includeJDay = TRUE + ) + }) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c( + "ranwalk", "intercept", + "inversegamma", "catlistlength" + ) + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - expect_true('dtype2.p' %in% row.names(results$BUGSoutput$summary)) - expect_true('dtype3.p' %in% row.names(results$BUGSoutput$summary)) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, modeltype = c('indran', 'centering', - 'halfcauchy', 'catlistlength')) - expect_identical(results$SPP_NAME, 'a') + 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_observations", "sparta_version" + ) + ) + expect_true("dtype2.p" %in% row.names(results$BUGSoutput$summary)) + expect_true("dtype3.p" %in% row.names(results$BUGSoutput$summary)) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, modeltype = c( + "indran", "centering", + "halfcauchy", "catlistlength" + ) + ) + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations", "sparta_version")) - expect_true('dtype2.p' %in% row.names(results$BUGSoutput$summary)) - expect_true('dtype3.p' %in% row.names(results$BUGSoutput$summary)) - + 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_observations", "sparta_version" + ) + ) + expect_true("dtype2.p" %in% row.names(results$BUGSoutput$summary)) + expect_true("dtype3.p" %in% row.names(results$BUGSoutput$summary)) + sink() - -}) +}) test_that("Test occDetFunc using regions and region aggregates", { - - sink(file=ifelse(Sys.info()["sysname"] == "Windows", - "NUL", - "/dev/null")) - - results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, - modeltype = c("ranwalk", "halfcauchy"), - regional_codes = regionsDF, - region_aggs = list(agg1 = c('region1', 'region2'))) - - expect_identical(results$SPP_NAME, 'a') + # skip the test if jagsUI is not installed + if ((!requireNamespace("jagsUI", quietly = TRUE))) { + skip("jagsUI software not available") + } + + sink(file = ifelse(Sys.info()["sysname"] == "Windows", + "NUL", + "/dev/null" + )) + + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, + modeltype = c("ranwalk", "halfcauchy"), + regional_codes = regionsDF, + region_aggs = list(agg1 = c("region1", "region2")) + ) + + expect_identical(results$SPP_NAME, "a") expect_identical(results$n.iter, 50) - 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_observations","sparta_version", - "regions", "region_aggs")) - expect_identical(results$regions, - c("region1", "region2", "region3")) + 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_observations", "sparta_version", + "regions", "region_aggs" + ) + ) + expect_identical( + results$regions, + c("region1", "region2", "region3") + ) expect_identical(names(results$region_aggs), "agg1") RNs <- row.names(results$BUGSoutput$summary) expect_true("a_region1[1]" %in% RNs) @@ -460,143 +668,183 @@ test_that("Test occDetFunc using regions and region aggregates", { expect_true("psi.fs.r_region2[1]" %in% RNs) expect_true("psi.fs.r_region3[1]" %in% RNs) expect_true("psi.fs.r_agg1[1]" %in% RNs) - + # test with a region with no data - regionsempty<-regionsDF + regionsempty <- regionsDF regionsempty$region4 <- 0 - - expect_warning(results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, - modeltype = c("ranwalk", "halfcauchy"), - regional_codes = regionsempty, - region_aggs = list(agg1 = c('region1', 'region2'))), - 'The following regions have no data and should not be modelled: region4 - These regions will not be included in the model') - + + expect_warning( + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, + modeltype = c("ranwalk", "halfcauchy"), + regional_codes = regionsempty, + region_aggs = list(agg1 = c("region1", "region2")) + ), + "The following regions have no data and should not be modelled: region4 - These regions will not be included in the model" + ) + # test with regional_codes not as a dataframe - expect_error(results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, - modeltype = c("ranwalk", "halfcauchy"), - regional_codes = as.matrix(regionsDF), - region_aggs = list(agg1 = c('region1', 'region2'))), - 'regional_codes should be a data.frame') - + expect_error( + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, + modeltype = c("ranwalk", "halfcauchy"), + regional_codes = as.matrix(regionsDF), + region_aggs = list(agg1 = c("region1", "region2")) + ), + "regional_codes should be a data.frame" + ) + # test with NAs in regional_codes regionsNA <- regionsDF - regionsNA[1,3] <- NA - - expect_warning(results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, - modeltype = c("ranwalk", "halfcauchy"), - regional_codes = regionsNA, - region_aggs = list(agg1 = c('region1', 'region2'))), - "NAs are present in regional_codes, these will be replaced with 0s") - + regionsNA[1, 3] <- NA + + expect_warning( + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, + modeltype = c("ranwalk", "halfcauchy"), + regional_codes = regionsNA, + region_aggs = list(agg1 = c("region1", "region2")) + ), + "NAs are present in regional_codes, these will be replaced with 0s" + ) + # test with sites in multiple regions regionsmulti <- regionsDF - regionsmulti[1,3] <- 1 - - expect_warning(results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, - modeltype = c("ranwalk", "halfcauchy"), - regional_codes = regionsmulti, - region_aggs = list(agg1 = c('region1', 'region2'))), - '1 sites are assigned to more than one region in regional_codes and will be removed') - + regionsmulti[1, 3] <- 1 + + expect_warning( + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, + modeltype = c("ranwalk", "halfcauchy"), + regional_codes = regionsmulti, + region_aggs = list(agg1 = c("region1", "region2")) + ), + "1 sites are assigned to more than one region in regional_codes and will be removed" + ) + # test with sites in no regions regionsmulti <- regionsDF - regionsmulti[1,2] <- 0 - - expect_warning(results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, - modeltype = c("ranwalk", "halfcauchy"), - regional_codes = regionsmulti, - region_aggs = list(agg1 = c('region1', 'region2'))), - '1 sites are not assigned to a region in regional_codes and will be removed') - + regionsmulti[1, 2] <- 0 + + expect_warning( + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, + modeltype = c("ranwalk", "halfcauchy"), + regional_codes = regionsmulti, + region_aggs = list(agg1 = c("region1", "region2")) + ), + "1 sites are not assigned to a region in regional_codes and will be removed" + ) + # test for sites in occurence data but not regions - regionsmissing <- regionsDF[2:50,] - - expect_warning(results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, - modeltype = c("ranwalk", "halfcauchy"), - regional_codes = regionsmissing, - region_aggs = list(agg1 = c('region1', 'region2'))), - '1 sites are in occurrence data but not in regional_codes and will be removed') - + regionsmissing <- regionsDF[2:50, ] + + expect_warning( + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, + modeltype = c("ranwalk", "halfcauchy"), + regional_codes = regionsmissing, + region_aggs = list(agg1 = c("region1", "region2")) + ), + "1 sites are in occurrence data but not in regional_codes and will be removed" + ) + # test for regional aggregate containing unknown region - expect_error(results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, - modeltype = c("ranwalk", "halfcauchy"), - regional_codes = regionsDF, - region_aggs = list(agg1 = c('region1', 'region2','region4'))), - 'Aggregate members [region4] not in regional_codes column names [region1, region2, region3]',fixed=TRUE) - + expect_error( + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, + modeltype = c("ranwalk", "halfcauchy"), + regional_codes = regionsDF, + region_aggs = list(agg1 = c("region1", "region2", "region4")) + ), + "Aggregate members [region4] not in regional_codes column names [region1, region2, region3]", + fixed = TRUE + ) + # test for regional aggregate without regional_codes - expect_error(results <- occDetFunc(taxa_name = 'a', - n_iterations = 50, - burnin = 15, - occDetdata = visitData$occDetdata, - spp_vis = visitData$spp_vis, - write_results = FALSE, - seed = 111, - modeltype = c("ranwalk", "halfcauchy"), - regional_codes = NULL, - region_aggs = list(agg1 = c('region1', 'region2'))), - 'Cannot use regional aggregates if regional_codes is not supplied') + expect_error( + results <- occDetFunc( + taxa_name = "a", + n_iterations = 50, + burnin = 15, + occDetdata = visitData$occDetdata, + spp_vis = visitData$spp_vis, + write_results = FALSE, + seed = 111, + modeltype = c("ranwalk", "halfcauchy"), + regional_codes = NULL, + region_aggs = list(agg1 = c("region1", "region2")) + ), + "Cannot use regional aggregates if regional_codes is not supplied" + ) sink() - -}) +}) test_that("Test occDetFunc with empty species post nyr filter", { - - sink(file=ifelse(Sys.info()["sysname"] == "Windows", - "NUL", - "/dev/null")) - - expect_warning(results <- occDetFunc(taxa_name = 'aa', - n_iterations = 50, - burnin = 15, - occDetdata = visitData_filterError$occDetdata, - spp_vis = visitData_filterError$spp_vis, - write_results = FALSE, - seed = 111, - modeltype = c("ranwalk", "halfcauchy")), - 'aa has insufficient data after site filtering. Either decrease nyr or change the criterion') + # skip the test if jagsUI is not installed + if ((!requireNamespace("jagsUI", quietly = TRUE))) { + skip("jagsUI software not available") + } + + sink(file = ifelse(Sys.info()["sysname"] == "Windows", + "NUL", + "/dev/null" + )) + + expect_warning( + results <- occDetFunc( + taxa_name = "aa", + n_iterations = 50, + burnin = 15, + occDetdata = visitData_filterError$occDetdata, + spp_vis = visitData_filterError$spp_vis, + write_results = FALSE, + seed = 111, + modeltype = c("ranwalk", "halfcauchy") + ), + "aa has insufficient data after site filtering. Either decrease nyr or change the criterion" + ) sink() - }) diff --git a/tests/testthat/testoccDetModel.r b/tests/testthat/testoccDetModel.r index 99c39b8..7885ddc 100644 --- a/tests/testthat/testoccDetModel.r +++ b/tests/testthat/testoccDetModel.r @@ -1,23 +1,23 @@ context("Test occDetModel") # Create data -n <- 15000 #size of dataset +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", "%Y/%m/%d")) -last <- as.Date(strptime(paste(2010+(nyr-1),"/12/31", sep=''), "%Y/%m/%d")) -dt <- last-first -rDates <- first + (runif(nSamples)*dt) +first <- as.Date(strptime("2010/01/01", "%Y/%m/%d")) +last <- as.Date(strptime(paste(2010 + (nyr - 1), "/12/31", sep = ""), "%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) +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) @@ -25,25 +25,35 @@ survey <- sample(rDates, size = n, TRUE) df <- data.frame(taxa, site, survey) test_that("Test occDetModel", { - - sink(file=ifelse(Sys.info()["sysname"] == "Windows", - "NUL", - "/dev/null")) - suppressWarnings(results <- occDetModel(taxa = df$taxa, - write_results = FALSE, - site = df$site, - survey = df$survey, - species_list = c('a','m','g'), - n_iterations = 100, - burnin = 10, - thinning = 2, - seed = 111)) + # skip the test if jagsUI is not installed + if ((!requireNamespace("jagsUI", quietly = TRUE))) { + skip("jagsUI software not available") + } + + sink(file = ifelse(Sys.info()["sysname"] == "Windows", + "NUL", + "/dev/null" + )) + suppressWarnings(results <- occDetModel( + taxa = df$taxa, + write_results = FALSE, + site = df$site, + survey = df$survey, + species_list = c("a", "m", "g"), + n_iterations = 100, + burnin = 10, + thinning = 2, + seed = 111 + )) sink() - expect_identical(names(results), c('a','m','g')) - expect_identical(names(results[[1]]), - c("model", "BUGSoutput", "parameters.to.save", "model.file", - "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", - "nsites", "nvisits", "species_observations","sparta_version")) - -}) \ No newline at end of file + expect_identical(names(results), c("a", "m", "g")) + expect_identical( + names(results[[1]]), + c( + "model", "BUGSoutput", "parameters.to.save", "model.file", + "n.iter", "DIC", "SPP_NAME", "min_year", "max_year", "sites_included", + "nsites", "nvisits", "species_observations", "sparta_version" + ) + ) +})