R/projectModel.R

In julienvollering/MIAmaxent: A Modular, Integrated Approach to Maximum Entropy Distribution Modeling

#' Project model across explanatory data.
#'
#' \code{projectModel} calculates model predictions for any points where values
#' of the explanatory variables in the model are known. It can be used to get
#' model predictions for the training data, or to project the model to a new
#' space or time.
#'
#' Missing data (NA) for a continuous variable will result in NA output for that
#' point. Missing data for a categorical variable is treated as belonging to
#' none of the categories.
#'
#' When \code{rescale = FALSE} the scale of the maxent model output (PRO or raw)
#' returned by this function is dependent on the data used to train the model.
#' For example, a location with PRO = 2 can be interpreted as having a
#' probability of presence twice as high as an average site in the
#' \emph{training} data (Halvorsen, 2013, Halvorsen et al., 2015). When
#' \code{rescale = TRUE}, the output is linearly rescaled with respect to the
#' data onto which the model is projected. In this case, a location with PRO = 2
#' can be interpreted as having a probability of presence twice as high as an
#' average site in the \emph{projection} data. Similarly, raw values are on a
#' scale which is dependent on the size of either the training data extent
#' (\code{rescale = FALSE}) or projection data extent (\code{rescale = TRUE}).
#'
#' @param model The model to be projected. This may be the object returned by
#'   \code{\link{chooseModel}}, or the 'selectedmodel' returned by
#'   \code{\link{selectEV}}.
#' @param transformations Transformation functions used to create the derived
#'   variables in the model. I.e. the 'transformations' returned by
#'   \code{\link{deriveVars}}. Equivalently, the full file pathway of the
#'   'transformations.Rdata' file saved as a result of \code{\link{deriveVars}}.
#' @param data Data frame of all the explanatory variables (EVs) included in the
#'   model (see \code{\link{readData}}). Alternatively, an object of class
#'   'SpatRaster' containing rasters for all EVs included in the model. Column
#'   or raster names must match EV names.
#' @param clamping Logical. Do clamping \emph{sensu} Phillips et al. (2006).
#'   Default is \code{FALSE}.
#' @param raw Logical. Return raw maxent output instead of probability ratio
#'   output (PRO)? Default is FALSE. Irrelevant for "LR"-type models.
#' @param rescale Logical. Linearly rescale model output (PRO or raw) with
#'   respect to the projection \code{data}? This has implications for the
#'   interpretation of output values with respect to reference values (e.g. PRO
#'   = 1). See details. Irrelevant for "LR"-type models.
#' @param filename Full file pathway to write raster model predictions if
#'   \code{data} is an object of class 'SpatRaster'. File format is inferred
#'   from the filename extension as in \code{terra::writeRaster}.

#'
#'@return List of 2: \enumerate{ \item output: A data frame with the model
#'  output in column 1 and the corresponding explanatory data in subsequent
#'  columns, or a raster containing predictions if \code{data} is a SpatRaster.
#'  \item ranges: A list showing the range of \code{data} compared to the
#'  training data, on a 0-1 scale.} If \code{data} is a SpatRaster, the output
#'  is also plotted.
#'
#'
#'@references Halvorsen, R. (2013) A strict maximum likelihood explanation of
#'  MaxEnt, and some implications for distribution modelling. Sommerfeltia, 36,
#'  1-132.
#'@references Halvorsen, R., Mazzoni, S., Bryn, A. & Bakkestuen, V. (2015)
#'  Opportunities for improved distribution modelling practice via a strict
#'  maximum likelihood interpretation of MaxEnt. Ecography, 38, 172-183.
#'@references Phillips, S.J., Anderson, R.P. & Schapire, R.E. (2006) Maximum
#'  entropy modeling of species geographic distributions. Ecological Modelling,
#'  190, 231-259.
#'
#' @examples
#' \dontrun{
#' # From vignette:
#' EVfiles <- c(list.files(system.file("extdata", "EV_continuous", package="MIAmaxent"),
#'                         full.names=TRUE),
#'              list.files(system.file("extdata", "EV_categorical", package="MIAmaxent"),
#'                         full.names=TRUE))
#' EVstack <- rast(EVfiles)
#' names(EVstack) <- gsub(".asc", "", basename(EVfiles))
#' grasslandPreds <- projectModel(model = grasslandmodel,
#'                                transformations = grasslandDVs$transformations,
#'                                data = EVstack)
#' grasslandPreds
#' }
#'
#'@export


projectModel <- function(model, transformations, data, clamping = FALSE,
                         raw = FALSE, rescale = FALSE, filename = NULL) {

  dvnamesni <- names(model$betas)[grep(":", names(model$betas), invert = TRUE)]
  evnames <- unique(sub("_.*", "", dvnamesni))

  map <- FALSE
  if (inherits(data, c("RasterLayer","RasterStack","RasterBrick"))) {
    data <- terra::rast(data)
  }
  if (inherits(data, "SpatRaster")) {
    map <- TRUE
    names(data) <- make.names(names(data), allow_ = FALSE)
    evstack <- data[[evnames]]
    data <- terra::as.data.frame(evstack, na.rm = TRUE)
    cells <- as.numeric(row.names(data))
  }

  for (i in evnames) {
    if (sum(colnames(data) == i) != 1) {
      stop(paste(i, "must be represented in 'data' (exactly once)"),
           call. = FALSE) }
  }

  alltransf <- .load.transf(transformations)
  .check.dvs.in.transf(dvnamesni, alltransf)

  Ranges <- lapply(evnames, function(x) {
    evdata <- data[, x]
    anevtransf <- alltransf[grepl(paste0(x, "_"), names(alltransf))][[1]]
    xnull <- environment(anevtransf)$xnull
    if (inherits(xnull, c("numeric", "integer"))) {
      L <- (evdata - range(xnull)[1]) / diff(range(xnull))
      return(range(L, na.rm = TRUE))
    }
    if (inherits(xnull, c("factor", "character"))) {
      if (all(evdata %in% xnull)) {return("inside")} else {return("outside")}
    }
  })
  names(Ranges) <- evnames

  dvdatani <- lapply(dvnamesni, function(x) {
    evdata <- data[, sub("_.*", "", x)]
    y <- alltransf[[paste0(x, "_transf")]](evdata)
    if (clamping == TRUE) {
      y[y > 1] <- 1
      y[y < 0] <- 0
    }
    return(y)
  })
  newdata <- as.data.frame(do.call(cbind, dvdatani))
  names(newdata) <- dvnamesni
  type <- if (inherits(model, "MIAmaxent_iwlr")) {
    ifelse(raw == TRUE, "raw", "PRO")
  } else { "response" }

  if (any(is.nan(unlist(newdata)))) {
    nans <- which(apply(newdata, 1, function(x) {any(is.nan(x))}))
    warning("Transformed 'data' has ", length(nans), " rows with NaN. Predictions for these will be NA.",
            call. = FALSE)
    preds <- rep(NA, nrow(newdata))
    preds[-nans] <- stats::predict(model, newdata[-nans,], type)
      } else {
    preds <- stats::predict(model, newdata, type)
  }

  if (inherits(model, "MIAmaxent_iwlr") && rescale == TRUE) {
    if (raw == TRUE) { preds <- preds/sum(preds, na.rm = TRUE)
    } else { preds <- (preds/sum(preds, na.rm = TRUE)) * sum(!is.na(preds)) }
  }

  Output <- data.frame(preds, data)
  colnames(Output)[1] <- type

  if (map == TRUE) {
    values <- rep(NA, terra::ncell(evstack))
    values[cells] <- Output[,1]
    outraster <- evstack[[1]]
    outraster <- terra::setValues(outraster, values)
    names(outraster) <- type
    if (!is.null(filename)) {
      terra::writeRaster(outraster, filename)
    }
    terra::plot(outraster)
    return(list(output = outraster, ranges = Ranges))
  } else {
    return(list(output = Output, ranges = Ranges))
  }

}