Multiple bug fixes, added priors to engine_breg and allowed more th…

…an one offset (#7)

NAMESPACE

@@ -57,6 +57,7 @@ export(add_constrain_connectivity)
 export(add_constrain_dispersal)
 export(add_latent_spatial)
 export(add_log)
+export(add_offset)
 export(add_offset_range)
 export(add_predictor_range)
 export(add_predictors)
@@ -118,6 +119,7 @@ exportMethods(add_constrain_connectivity)
 exportMethods(add_constrain_dispersal)
 exportMethods(add_latent_spatial)
 exportMethods(add_log)
+exportMethods(add_offset)
 exportMethods(add_offset_range)
 exportMethods(add_predictor_range)
 exportMethods(add_predictors)

R/add_biodiversity.R

@@ -364,7 +364,7 @@ methods::setMethod(
 #'
 #' @param x A [`data.frame`], [`sf`] or [`Spatial`]) object of biodiversity information
 #' @param field_occurrence A [`numeric`] or [`character`] location of biodiversity records.
-#' @param field_space A [`vector`] on the column names (Default: \code{'X'}, \code{'Y'})
+#' @param field_space A [`vector`] on the column names (Default: \code{'x'}, \code{'y'})
 #' @param ... Other parameters passed down
 #'
 #' @import sf
@@ -373,7 +373,7 @@ methods::setMethod(
 #' @keywords internal
 #' @noRd
 
-format_biodiversity_data <- function(x, field_occurrence, field_space = c('X','Y'),...){
+format_biodiversity_data <- function(x, field_occurrence, field_space = c('x','y'),...){
   # Final checks
   if(inherits(x,'sf')) assertthat::assert_that(unique(sf::st_geometry_type(x)) %in% c('POINT','MULTIPOINT',
                                                                                       'POLYGON','MULTIPOLYGON'))
@@ -398,7 +398,7 @@ format_biodiversity_data <- function(x, field_occurrence, field_space = c('X','Y
       }
     }
     assertthat::assert_that( all(assertthat::has_name(x, field_space)),
-                             msg ='No spatial column found in the dataset. Specify manually or set to [X] and [Y].')
+                             msg ='No spatial column found in the dataset. Specify manually or set to [x] and [y].')
     # Select and format
     out <- subset(x, select = c(field_space, field_occurrence) ) %>%
       as_tibble()

R/add_offset.R

@@ -1,28 +1,111 @@
-#' Specify an offset object for a species
+#' Specify a spatial explicit offset
 #'
 #' @description
 #' Including offsets is another option to integrate spatial prior information
-#' in linear and additive regression models. Here, offset simply shift the intercept of
-#' the regression fit. Since offsets only make sense for linear regressions (and not for instance
+#' in linear and additive regression models. Offsets shift the intercept of
+#' the regression fit by a certain amount. Although only one offset can be added
+#' to a regression model, it is possible to combine several spatial-explicit estimates into
+#' one offset by calculating the sum of all spatial-explicit layers.
+#'
+#' @details
+#' This function allows to set any specific offset to a regression model. The offset
+#' has to be provided as spatial [`RasterLayer`] object. This function simply adds the layer to
+#' a [`distribution()`] object.
+#' **Note that any transformation of the offset (such as \code{log}) has do be done externally!**
+#'
+#' If the layer is range and requires additional formatting, consider using the
+#' function [`add_offset_range()`] which has additional functionalities such such distance
+#' transformations.
+#'
+#' @note
+#' Since offsets only make sense for linear regressions (and not for instance
 #' regression tree based methods such as [engine_bart]), they do not work for all engines.
+#' Offsets specified for non-supported engines are ignored during the estimation
+#'
+#' @param x [distribution()] (i.e. [`BiodiversityDistribution-class`]) object.
+#' @param layer A [`sf`] or [`RasterLayer`] object with the range for the target feature.
+#' @param add [`logical`] specifiying whether new offset is to be added. Setting
+#' this parameter to \code{FALSE} replaces the current offsets with the new one (Default: \code{TRUE}).
+#' @param ... Other parameters or arguments (currently not supported)
+#' @references
+#' * Merow, C., Allen, J.M., Aiello-Lammens, M., Silander, J.A., 2016. Improving niche and range estimates with Maxent and point process models by integrating spatially explicit information. Glob. Ecol. Biogeogr. 25, 1022–1036. https://doi.org/10.1111/geb.12453
+#' @family offset
+#' @keywords prior, offset
+#' @examples
+#' \dontrun{
+#'  x <- distribution(background) %>%
+#'    add_predictors(covariates) %>%
+#'    add_offset(samplingBias)
+#' }
+#' @name add_offset
+NULL
+
+#' @name add_offset
+#' @rdname add_offset
+#' @exportMethod add_offset
+#' @export
+methods::setGeneric(
+  "add_offset",
+  signature = methods::signature("x", "layer"),
+  function(x, layer, add = TRUE) standardGeneric("add_offset"))
+
+#' @name add_offset
+#' @rdname add_offset
+#' @usage \S4method{add_offset}{BiodiversityDistribution, raster}(x, layer)
+methods::setMethod(
+  "add_offset",
+  methods::signature(x = "BiodiversityDistribution", layer = "RasterLayer"),
+  function(x, layer, add = TRUE) {
+    assertthat::assert_that(inherits(x, "BiodiversityDistribution"),
+                            is.Raster(layer),
+                            is.logical(add)
+    )
+    # Messenger
+    if(getOption('ibis.setupmessages')) myLog('[Setup]','green','Adding spatial explicit offset...')
+    ori.name <- names(layer)
+
+    # Check that background and range align, otherwise raise error
+    if(compareRaster(layer, x$background,stopiffalse = FALSE)){
+      warning('Supplied layer does not align with background! Aligning them now...')
+      layer <- alignRasters(layer, x$background, method = 'bilinear', func = mean, cl = FALSE)
+      names(layer) <- ori.name
+    }
+    # Check whether an offset exists already
+    if(!is.Waiver(x$offset) && add){
+      # Add to current object
+      of <- x$offset
+      layer <- raster::resample(layer, of, method = 'bilinear', func = mean, cl = FALSE)
+      names(layer) <- ori.name # In case the layer name got lost
+      of <- raster::stack(of) |> raster::addLayer(layer)
+      x <- x$set_offset(of)
+    } else {
+      # Add as a new offset
+      x <- x$set_offset(layer)
+    }
+    return(x)
+  }
+)
+
+#' Specify a expert-based species range as offset
+#'
+#' @description
+#' This function has additional nuance compared to the more generic
+#' [`add_offset()`], allowing to specify transformations and modifications.
+#' It is customized specifically for expert-based ranges as offsets.
 #'
 #' @details
 #' As options allow specifying including the range either as binary or using distance
 #' as parameter following Merow et al. (2017). In this case the existing range needs
-#' to be transformed first, for instance using the [bossMaps] package.
-#' @note
-#' Only one offset per distribution object can currently be specified. This might
-#' change in the future.
-#' @param x [distribution()] (i.e. [`BiodiversityDistribution-class`]) object.
-#' @param range A [`sf`] or [`RasterLayer`]object with the range for the target feature.
+#' to be transformed first, for instance using the \pkg{bossMaps} package.
+#' @inheritParams add_offset
 #' @param method [`character`] describing how the range should be included (Options: \code{"binary"} | \code{"distance"}).
-#' @param name [`character`] A name of the species if the offset is not specified.
 #' @param distance_max Numeric threshold on the maximum distance (Default: \code{150000} [m]).
 #' @seealso [bossMaps]
 #' @references
 #' * Merow, C., Wilson, A.M., Jetz, W., 2017. Integrating occurrence data and expert maps for improved species range predictions. Glob. Ecol. Biogeogr. 26, 243–258. https://doi.org/10.1111/geb.12539
 #' * Merow, C., Allen, J.M., Aiello-Lammens, M., Silander, J.A., 2016. Improving niche and range estimates with Maxent and point process models by integrating spatially explicit information. Glob. Ecol. Biogeogr. 25, 1022–1036. https://doi.org/10.1111/geb.12453
-#' @keywords prior
+#' @keywords prior, offset
+#' @family offset
 #' @name add_offset_range
 NULL
 
@@ -32,55 +115,70 @@ NULL
 #' @export
 methods::setGeneric(
   "add_offset_range",
-  signature = methods::signature("x", "range", "method"),
-  function(x, range, method = 'distance', name = 'range_distance', distance_max = 150000) standardGeneric("add_offset_range"))
+  signature = methods::signature("x", "layer", "method"),
+  function(x, layer, method = 'distance', distance_max = 150000, add = TRUE) standardGeneric("add_offset_range"))
 
 #' Function for when raster is directly supplied (precomputed)
 #' @name add_offset_range
 #' @rdname add_offset_range
-#' @usage \S4method{add_offset_range}{BiodiversityDistribution, raster}(x, range)
+#' @usage \S4method{add_offset_range}{BiodiversityDistribution, raster}(x, layer)
 methods::setMethod(
   "add_offset_range",
-  methods::signature(x = "BiodiversityDistribution", range = "RasterLayer"),
-  function(x, range, method = 'range_distance') {
+  methods::signature(x = "BiodiversityDistribution", layer = "RasterLayer"),
+  function(x, layer, method = 'range_distance', add = TRUE) {
     assertthat::assert_that(inherits(x, "BiodiversityDistribution"),
-                            is.Raster(range),
-                            is.character(method)
+                            is.Raster(layer),
+                            is.character(method),
+                            is.logical(add)
     )
     # Messager
     if(getOption('ibis.setupmessages')) myLog('[Setup]','green','Adding range offset...')
 
+    # Save name
+    ori.name <- names(layer)
+
     # Check that background and range align, otherwise raise error
-    if(compareRaster(range, x$background,stopiffalse = FALSE)){
+    if(compareRaster(layer, x$background,stopiffalse = FALSE)){
       warning('Supplied range does not align with background! Aligning them now...')
-      range <- alignRasters(range, x$background, method = 'bilinear', func = mean, cl = FALSE)
+      layer <- alignRasters(layer, x$background, method = 'bilinear', func = mean, cl = FALSE)
+      names(layer) <- ori.name # In case the layer name got lost
     }
-    names(range) <- method
 
-    # Add as a new offset
-    x <- x$set_offset(range)
+    # Check whether an offset exists already
+    if(!is.Waiver(x$offset) && add){
+      # Add to current object
+      of <- x$offset
+      layer <- raster::resample(layer, of, method = 'bilinear', func = mean, cl = FALSE)
+      names(layer) <- ori.name # In case the layer name got lost
+      of <- raster::stack(of) |> raster::addLayer(layer)
+      x <- x$set_offset(of)
+    } else {
+      # Add as a new offset
+      x <- x$set_offset(layer)
+    }
     return(x)
   }
 )
 
 #' @name add_offset_range
 #' @rdname add_offset_range
-#' @usage \S4method{add_offset_range}{BiodiversityDistribution,sf, vector}(x, range, method)
+#' @usage \S4method{add_offset_range}{BiodiversityDistribution, sf, vector}(x, layer, method)
 methods::setMethod(
   "add_offset_range",
-  methods::signature(x = "BiodiversityDistribution", range = "sf", method = "character"),
-  function(x, range, method = 'distance', name = 'range_distance', distance_max = 150000 ) {
+  methods::signature(x = "BiodiversityDistribution", layer = "sf", method = "character"),
+  function(x, layer, method = 'distance', distance_max = 150000, add = TRUE ) {
     assertthat::assert_that(inherits(x, "BiodiversityDistribution"),
                             is.character(method),
-                            inherits(range, 'sf'),
+                            inherits(layer, 'sf'),
                             method %in% c('binary','distance'),
-                            is.null(distance_max) || is.numeric(distance_max)
+                            is.null(distance_max) || is.numeric(distance_max),
+                            is.logical(add)
     )
     # Messager
     if(getOption('ibis.setupmessages')) myLog('[Setup]','green','Adding range offset...')
 
     # Reproject if necessary
-    if(st_crs(range) != sf::st_crs(x$background)) range <- sf::st_transform(range, sf::st_crs(x$background))
+    if(st_crs(layer) != sf::st_crs(x$background)) layer <- sf::st_transform(layer, sf::st_crs(x$background))
 
     # Template raster for background
     if(!is.Waiver(x$predictors)){
@@ -92,41 +190,51 @@ methods::setMethod(
     }
 
     # Rasterize the range
-    range$id <- 1:nrow(range) # Assign an id if not already existing
+    layer$id <- 1:nrow(layer) # Assign an id if not already existing
     if( 'fasterize' %in% installed.packages()[,1] ){
-      ras_range <- fasterize::fasterize(range, temp, field = NULL)
+      ras_range <- fasterize::fasterize(layer, temp, field = NULL)
     } else {
-      ras_range <- raster::rasterize(range, temp,field = NULL)
+      ras_range <- raster::rasterize(layer, temp,field = NULL)
     }
 
     # -------------- #
     if(method == 'binary'){
       dis <- ras_range
-      names(dis) <- 'binary_range'
+      names(dis) <- "range_binary"
     } else if(method == 'distance'){
       # TODO: The below can be much more sophisticated.
       # - For instance adding a exponential decay
       # Calculate the linear distance
       dis <- raster::distance(ras_range)
       dis <- raster::mask(dis, x$background)
       # Set areas not intersecting with range to 0
-      suppressMessages( suppressWarnings( range <- sf::st_buffer(range, 0)) )
+      suppressMessages( suppressWarnings( layer <- sf::st_buffer(layer, 0)) )
 
       suppressMessages(
         dis <- raster::mask(dis,
-                        x$background[unlist( st_intersects(range, x$background) ),]
+                        x$background[unlist( sf::st_intersects(layer, x$background) ),]
         )
       )
       # If max distance is specified
       if(!is.null(distance_max)) dis[dis > distance_max] <- NA # Set values above threshold to NA
       # Convert to relative for better scaling
-      dis <- 1 - predictor_transform(dis,option = 'norm') #1 - (dis / cellStats(dis,'max') )
-      names(dis) <- name
+      dis <- 1 - predictor_transform(dis, option = 'norm') #1 - (dis / cellStats(dis,'max') )
+      names(dis) <- "range_distance"
     }
 
-    # Add as a new offset
-    x <- x$set_offset(dis)
+    # Check whether an offset exists already
+    if(!is.Waiver(x$offset) && add){
+      # Add to current object
+      of <- x$offset
+      ori.name <- names(dis)
+      dis <- raster::resample(dis, of, method = 'bilinear', func = mean, cl = FALSE)
+      names(dis) <- ori.name # In case the layer name got lost
+      of <- raster::stack(of) |> raster::addLayer(dis)
+      x <- x$set_offset(of)
+    } else {
+      # Add as a new offset
+      x <- x$set_offset(dis)
+    }
     return(x)
   }
 )
-

R/add_predictors.R

@@ -53,6 +53,9 @@ NULL
 #' Not every engine does actually need covariates. For instance it is perfectly legit
 #' to fit a model with only occurrence data and a spatial latent effect ([add_latent]).
 #' This correspondents to a spatial kernel density estimate.
+#'
+#' Certain names such \code{"offset"} are forbidden as predictor variable names. The function
+#' will return an error message if these are used.
 #' @aliases add_predictors
 #' @examples
 #' \dontrun{
@@ -134,6 +137,12 @@ methods::setMethod(
       names(env) <- names
     }
 
+    # Check that all names allowed
+    problematic_names <- grep("offset|w|weight|spatial_offset|Intercept|spatial.field", names(env),fixed = TRUE)
+    if( length(problematic_names)>0 ){
+      stop(paste0("Some predictor names are not allowed as they might interfere with model fitting:", paste0(names(env)[problematic_names],collapse = " | ")))
+    }
+
     # If priors have been set, save them in the distribution object
     if(!is.null(priors)) {
       # FIXME: Ideally attempt to match varnames against supplied predictors vis match.arg or similar

R/bdproto-biodiversitydataset.R

@@ -76,8 +76,13 @@ BiodiversityDatasetCollection <- bdproto(
     assertthat::assert_that(is.Id(id) || is.character(id))
     # Get data
     o <- self$get_data(id)
+    o <- guess_sf(o)
+    # Add lowercase coordinates for consistency
+    o$x <- sf::st_coordinates(o)[,1]
+    o$y <- sf::st_coordinates(o)[,2]
+
     # Return coordinates
-    if(hasName(o,'geom')) sf::st_coordinates(o$geom) else o[,c('X','Y')]
+    if(hasName(o,'geom')) sf::st_coordinates(o) else o[,c('x','y')]
   },
   # Remove a specific biodiversity dataset by id
   rm_data = function(self, id) {

R/bdproto-biodiversitydistribution.R

@@ -36,7 +36,7 @@ BiodiversityDistribution <- bdproto(
     # Query information from the distribution object
     ex =  self$show_background_info()
     pn = ifelse(is.Waiver(self$get_predictor_names()),'None',name_atomic(self$get_predictor_names(), "predictors"))
-    of = ifelse(is.Waiver(self$offset), '', paste0( "\n  offset:         <", self$get_offset(),">" ) )
+    of = ifelse(is.Waiver(self$offset), '', paste0( "\n  offset:         <", name_atomic(self$get_offset()),">" ) )
     pio = ifelse(is.Waiver(self$priors), '<Default>', paste0('Priors specified (',self$priors$length(), ')') )
 
     message(paste0('\033[1m','\033[36m','<',self$name(),'>','\033[39m','\033[22m',
@@ -164,15 +164,23 @@ BiodiversityDistribution <- bdproto(
   # Set offset
   # FIXME: For logical consistency could define a new bdproto object
   set_offset = function(self, x){
-    assertthat::assert_that(inherits(x, "Raster"))
-    if (!is.Waiver(self$offset)) warning("Overwriting previously defined offset.")
+    assertthat::assert_that(is.Raster(x))
     bdproto(NULL, self, offset = x )
   },
   # Get offset (print name)
   get_offset = function(self){
-    if(is.Waiver(self$offset)) return('None')
+    if(is.Waiver(self$offset)) return( self$offset() )
     names( self$offset )
   },
+  # Plot offset
+  plot_offset = function(self){
+    if(is.Waiver(self$offset)) return( self$offset() )
+    if(raster::nlayers(self$offset)>1){
+      of <- sum(self$offset,na.rm = TRUE)
+      of <- raster::mask(of, self$background)
+    } else {of <- self$offset}
+    raster::plot(of, col = ibis_colours$viridis_orig, main = "Combined offset")
+  },
   # Get log
   get_log = function(self){
     if(is.Waiver(self$log)){