Merge pull request #7 from drnickisaac/master

Pull in Nicks changes

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: BRCindicators
 Type: Package
 Title: Creating multispecies biodiversity indicators
-Version: 1.3.5
-Date: 2020-02-20
+Version: 1.3.7
+Date: 2021-03-09
 Authors@R: c(person("Tom", "August", role = c("aut", "cre"), email = "[email protected]"),
              person("Gary", "Powney", role = c("aut")),
              person("Charlie", "Outhwaite", role = c("aut")),
@@ -27,7 +27,7 @@ VignetteBuilder: knitr
 License: GPL-3
 URL: https://github.com/BiologicalRecordsCentre/BRCindicators
 BugReports: https://github.com/BiologicalRecordsCentre/BRCindicators/issues
-RoxygenNote: 7.0.2
+RoxygenNote: 7.1.1
 Encoding: UTF-8
 Suggests: 
     testthat,

R/bayesian_meta_analysis.R

@@ -6,7 +6,8 @@
 #' @param plot Logical, should a trace plot be plotted to diagnose the model output?
 #' @param model The type of model to be used. See details.
 #' @param parallel if \code{TRUE} the model chains will be run in parallel using one fewer cores than
-#' are available on the computer. NOTE: this will typically not work for parallel use on cluster PCs.
+#' are available on the computer as default. NOTE: this will typically not work for parallel use on cluster PCs.
+#' @param n.cores if running the code in parallel this option specifies the number of cores to use. 
 #' @param incl.model if \code{TRUE} the model is added as an attribute of the object returned
 #' @param n.iter The number of iterations of the model to run. Defaults to 10,000 to avoid long run times
 #' though much longer runs are usually required to reach convergence
@@ -30,15 +31,16 @@
 #'  Includes three options: `seFromData` `Y1perfect` and `incl.2deriv`. See bayesian_meta_analysis for mode details. Using the default values `seFromData = FALSE` and `Y1perfect = TRUE` are the options used in Freeman  \emph{et al.} (2020).}
 #'  \item{\code{"smooth_det2"}}{ Equivalent to smooth with `seFromData = TRUE` and `Y1perfect = FALSE`. Retained for backwards compatability. Choosing this option will overwrite user-entered options for `seFromData` and `Y1perfect`.}
 #'  \item{\code{"smooth_det_sigtheta"}}{ Equivalent to smooth with `seFromData = FALSE` and `Y1perfect = FALSE`. Retained for backwards compatability. Choosing this option will overwrite user-entered options for `seFromData` and `Y1perfect`.}
-#'  \item{\code{"smooth_det"}}{ Specific variant of smooth_det2 - under review. Likely to be deprecated}
 #'  }
-#' @return Returns a dataframe with 4 columns: Year, Index, lower2.5, upper97.5. The last two columns are the credible intervals
+#' @return Returns a dataframe with 7 columns: Year, Index.Mprime, lowerCI.Mprime, upperCI.Mprime, Index.M, lowerCI.M and, upperCI.M. 
+#' Columns headed `M` and `Mprime` are means of the M and M' parameters  as defined in Freeman et al (2020). The 'upper' and 'lower' columns are the credible intervals, the width of which is defined by the `CI` argument.
+#' Note that M and M' are alternate ways of calculating the multispecies indicator: their means are nearly always virtually identical, but the uncertainty in M is usually much wider than in M'. See Freeman et al (2020) for more details.
 #' @import reshape2
 #' @importFrom boot inv.logit
 #' @importFrom coda mcmc.list as.mcmc
-#' @references Freeman, S.N., Isaac, N.J.B., Besbeas, P.T., Dennis, E.B. & Morgan, B.J.T. (2019) 
+#' @references Freeman, S.N., Isaac, N.J.B., Besbeas, P.T., Dennis, E.B. & Morgan, B.J.T. (2020) 
 #'             A generic method for estimating and smoothing multispecies biodiversity indices, robust to intermittent data. 
-#'             \emph{JABES}, in revision.
+#'             \emph{Journal of Agricultural Biological and Environmental Statistics}, in revision.
 #' @export
 #' @examples 
 #' # Create some example data in the format required
@@ -59,6 +61,7 @@ bma <- function (data,
                  plot = TRUE,
                  model = 'smooth',
                  parallel = FALSE,
+                 n.cores = parallel::detectCores()-1,
                  incl.model = TRUE,
                  n.iter = 1e4,
                  n.thin = 5,
@@ -108,10 +111,7 @@ bma <- function (data,
            model = "smooth"
            seFromData=FALSE
            Y1perfect = FALSE},
-         smooth_det = {
-           seFromData = TRUE
-           Y1perfect = FALSE
-         },
+         smooth_det = stop("smooth_det model has been deprecated"),
          random_walk = stop("Random walk model has been deprecated"),
          uniform = stop("Uniform model has been deprecated"),
          uniform_noeta = stop("Uniform model has been deprecated"),
@@ -152,7 +152,8 @@ bma <- function (data,
     # the user has specified that each species' time series should be scaled to start at a common value.
     # since the data are on the unbounded (log or logit) scale, we standardise them by addition/subtraction, rather than multiplication (as in rescale_species())
     # recall that rescale_indices has to be an integer
-    index <- t(apply(index, 1, function(x) rescale_indices + x - x[1]))
+    # without mising data this would be easy, but we have to identify the first year in each species' timeseries
+    index <- t(apply(index, 1, function(x) rescale_indices + x - x[min(which(!is.na(x)))]))
   }
 
 
@@ -182,22 +183,21 @@ bma <- function (data,
 
   if(model %in% c('smooth', 'smooth_stoch2', 'FNgr2')){
     # using row.names should ensure the same order in the bugs data
-    FY <- sapply(row.names(index), FUN = function(x){
-      min(data$year[!is.na(data$index) & data$species == x])
-    })
-    bugs_data[['FY']] <- FY - min(FY) + 1 # set lowest value to 1
-
+    #FY <- sapply(row.names(index), FUN = function(x){
+    #  min(data$year[!is.na(data$index) & data$species == x])
+    #})
+    #bugs_data[['FY']] <- FY - min(FY) + 1 # set lowest value to 1
+    bugs_data[['FY']] <- apply(index, 1, function(x) min(which(!is.na(x)))) # simpler alternative
   }
 
   # Setup parameters to monitor. NB Most of the model options have been deprecated, so much of this code is redundant
-  params = c("tau.spi")#, "sigma.obs")
+  params = c("tau.spi")
   if(model == 'smooth') params <- c(params, "Mprime")
-  #if(model %in% c('random_walk', 'uniform', 'uniform_noeta')) params <- c(params, "tau.eta")
-  #if(model %in% c('random_walk')) params <- c(params, "tau.I")
   if(model %in% c('smooth', 'smooth_stoch', 'smooth_det', 'FNgr',
                   'smooth_stoch2', 'FNgr2')) params <- c(params, "logLambda", "spgrowth", "M")
   if(model %in% c('smooth_stoch', 'smooth_det', 'FNgr')) params <- c(params, "tau.sg")
   if(model %in% c('smooth', 'smooth_stoch', 'smooth_det','smooth_stoch2')) params <- c(params, "beta", "taub")
+  if(incl.2deriv) params <- c(params, "t2dash")
   if(!seFromData) params <- c(params, "theta")
   if(save.sppars) {
     params <- c(params, "spindex")
@@ -211,7 +211,7 @@ bma <- function (data,
                             inits = NULL,
                             param = params,
                             parallel = parallel,
-                            n.cores = parallel::detectCores()-1,
+                            n.cores = n.cores,
                             model.file = bugs_path,
                             store.data = TRUE,
                             n.chains = 3,

R/bma_BUGS_models.R

@@ -42,7 +42,7 @@
 ################################################################################
 
 smoothing <- function() {'
-    ######################### Smoothing done here   #######################
+    ######################### Smoothing #######################
 
     beta[1] ~ dnorm(0, 0.000001)
     beta[2] ~ dnorm(0, 0.000001)
@@ -84,20 +84,20 @@ bma_model_Smooth <- function(incl.2deriv = FALSE,
     process_errors <- '
   # process errors
   tau.spi <- pow(sigma.spi,-2)
-  sigma.spi ~ dunif(0,30)'
+  sigma.spi ~ dunif(0.0001,30)'
 
     if(seFromData){
       obsErrors <- '
   # observation errors: one value per site-species
   for (s in 1:nsp){
    for (t in 1:nyears){
-    sigma.obs[s,t] ~ dunif(0, max_se) # for the missing values
+    sigma.obs[s,t] ~ dunif(0.0001, max_se) # for the missing values
   }}
     '
     } else {
       obsErrors <- '
   # observation error is constant
-  theta ~ dunif(0,10)'
+  theta ~ dunif(0.0001,10)'
     }
 
     return(paste(c(
@@ -175,6 +175,7 @@ bma_model_Smooth <- function(incl.2deriv = FALSE,
     #########################  second derivatives #######################
   
   I <- M
+  t2dash[1] <- 0
   t2dash[2] <- (I[2+1] - 2*I[2] + I[2-1])/1
   t2dash[nyears-1] <- (I[nyears] - 2*I[nyears-1] + I[nyears-2])/1
   t2dash[3] <- (-I[5]+16*I[4]-30*I[3]+16*I[2]-I[1] )/12
@@ -201,21 +202,25 @@ bma_model_Smooth <- function(incl.2deriv = FALSE,
 
 
 ################################################################################
-
-# OPTIONS UNDER REVIEW
+################################################################################
+# BEGIN DEPRECATED OPTIONS
+################################################################################
 
 bma_model_smooth_det <- function(){
   # Indicator defined by Growth rates, with Ruppert smoother (deterministic version)
   # Defined by equation 7 in Steve Freeman's document of 2/11/17
   # Also known as "smooth_indicator_1" in Steve's email of 2/11/17
   # this version takes standard errors on year 1 estimates. 
   # If actually zero then the model invents them (this is a fudge)
-
+  # THIS VERSION CONTAINS AN ERROR - to fix
+
   model <- '
   ###################  Define priors
   # process errors
-  tau.spi<-pow(sigma.spi,-2)
-  sigma.spi~dunif(0,30)
+  tau.spi <- pow(sigma.spi,-2)
+  sigma.spi ~ dunif(0,30)
+  tau.sg <- pow(sigma.sg,-2)
+  sigma.sg ~ dunif(0,1000)
 
   # observation errors
   # one value per site-species
@@ -266,12 +271,6 @@ bma_model_smooth_det <- function(){
   return(model)
 }
 
-
-
-
-################################################################################
-################################################################################
-# BEGIN DEPRECATED OPTIONS
 ################################################################################
 
 bma_model_FNgr2 <- function(){

R/butterflies_hs-data.R

@@ -1,6 +1,7 @@
 #' @name butterflies_hs
 #' @title Data - UK Butterfly Monitoring Scheme - Habitat Specialists
 #' @description This dataset from the UK Butterfly Monitoring Scheme has national abundance indices for 26 species regarded as habitat specialists from 1976-2017.
+#' There are 1000 rows of data. Ten species have a complete time-series (42 years of data); 15 species join late (9 in the 1970s, 3 in the 1980s and 3 in the 1990s); the Swallowtail (Papilio machaon) spans the full range of years but has no data for 1978.
 #' @docType data
 #' @format 
 #' There are four columns of data:

R/butterflies_wc-data.R

@@ -1,6 +1,7 @@
 #' @name butterflies_wc
 #' @title Data - UK Butterfly Monitoring Scheme - Wider Countryside Species
 #' @description This dataset from the UK Butterfly Monitoring Scheme has national abundance indices for 24 species regarded as wider countryside species from 1976-2017.
+#' There are 1005 lines of data, reflecting late entry into the time series for the Scotch Argus (Erebia aethops). Time-series for the other 23 species are complete.
 #' @docType data
 #' @usage data(butterflies_wc)
 #' @format 

R/dragonflies-data.R

@@ -1,7 +1,7 @@
 #' @name dragonflies
 #' @title Data - British Dragonfly Society, UK
 #' @description This dataset contains annual occupancy estimates for 45 species of dragonflies and damselflies in the UK, 1970-2015.
-#' There are 2039 rows: 44 species have an estimate for every year, but the Small red-eyed damselfly has no estimates prior to 2001 (the first year in which it was recorded in UK).
+#' There are 2039 rows: 44 species have an estimate for every year, but the Small red-eyed damselfly (Erythromma viridulum) has no estimates prior to 2001 (the first year in which it was recorded in UK).
 #' The models are based on biological records data from the British Dragonfly Society.
 #' The data are derived from Bayesian occupancy detection models, described in Outhwaite et al (2019).
 #' @usage data(dragonflies)

R/get_bmaBUGScode.R

@@ -13,7 +13,7 @@
 #'  Includes three options: `seFromData` `Y1perfect` and `incl.2deriv`. See bayesian_meta_analysis for mode details. Using the default values `seFromData = FALSE` and `Y1perfect = TRUE` are the options used in Freeman  \emph{et al.} (2020).}
 #'  \item{\code{"smooth_det2"}}{ Equivalent to smooth with `seFromData = TRUE` and `Y1perfect = FALSE`. Retained for backwards compatability. Choosing this option will overwrite user-entered options for `seFromData` and `Y1perfect`.}
 #'  \item{\code{"smooth_det_sigtheta"}}{ Equivalent to smooth with `seFromData = FALSE` and `Y1perfect = FALSE`. Retained for backwards compatability. Choosing this option will overwrite user-entered options for `seFromData` and `Y1perfect`.}
-#'  \item{\code{"smooth_det"}}{ Specific variant of smooth_det2 - under review. Likely to be deprecated}
+#'  \item{\code{"smooth_det"}}{ Specific variant of smooth_det2 - deprecated but retained here for backward compatibility}
 #'  }
 #' @export
 #' @examples
@@ -46,8 +46,8 @@ get_bmaBUGScode <- function(option="smooth",
                                                   incl.2deriv, 
                                                   seFromData=FALSE,
                                                   Y1perfect = FALSE)},
-         smooth_det = {model_code <- bma_model_smooth_det()},
          # Deprecated options can still be called directly from this function
+         smooth_det = {model_code <- bma_model_smooth_det()},
          random_walk = model_code <- bma_model_ranwalk(),
          uniform = model_code <- bma_model_uniform(),
          uniform_noeta = model_code <- bma_model_uniform_noeta(),