R/Fav.R

In fuzzySim: Fuzzy Similarity in Species Distributions

Fav <- function(model = NULL, obs = NULL, pred = NULL, n1n0 = NULL, sample.preval = NULL, method = "RBV", true.preval = NULL, inv = FALSE, verbosity = 2) {
  # version 2.0 (4 Mar 2025)

  if (!is.null(model)) {

    if (verbosity > 0) {
      if (!is.null(pred)) message("Argument 'pred' ignored in favour of 'model'.")
      if (!is.null(obs)) message("Argument 'obs' ignored in favour of 'model'.")
      if (!is.null(n1n0)) message("Argument 'n1n0' ignored in favour of 'model'.")
      if (!is.null(sample.preval)) message("Argument 'sample.preval' ignored in favour of 'model'.")
    }  # end if verbosity

    if (inv)
      stop("'inv=TRUE' requires 'pred' instead of 'model'\n(plus either one of 'obs', 'n1n0' or 'sample.preval'),\nand 'pred' must consist of favourability values.")

    #   if (is(model, "glm") || is(model, "gam")) {
    #     obs <- model$y
    #     pred <- model$fitted.values
    #   } else if (is(model, "gbm")) {
    #     obs <- model$data$y
    #     exp_y <- exp(model$fit)  # gbm provides 'fit' on the predictors' scale
    #     pred <- exp_y / (1 + exp_y)  # logit to probability
    #   } else if (is(model, "randomForest")) {
    #     obs <- as.integer(as.character(model$y))
    #     pred <- predict(model, type = "prob")[ , "1"]
    #   } else if (is(model, "bart")) {
    #     if (is.null(model$fit$data)) stop("'$fit$data' section not available in 'model'. Try computing the 'bart' model with 'keeptrees=TRUE', or providing other arguments here instead of 'model', i.e. 'pred' plus one of 'obs', 'n1n0' or 'sample.preval'.")
    #     obs <- model$fit$data@y  # requires model ran with keeptrees=TRUE
    #     pred <- stats::fitted(model, type = "response")
    #   } else stop("'model' is of a non-implemented class.")

    obspred <- modEvA::mod2obspred(model)
    obs <- obspred[ , "obs"]
    pred <- obspred[ , "pred"]
  }  # end if model

  # if(!is.null(obs) & !is.null(pred) & !class(pred) == "RasterLayer" & length(obs) != length(pred)) {
  #   stop("'obs' and 'pred' must have the same length (and be in the same order).")
  # }

  #if (!is.vector(obs)) stop("'obs' must be of class 'vector'.")
  #if (!(is.vector(pred) || inherits(pred, "SpatRaster") || inherits(rast(pred), "SpatRaster"))) stop("'pred' must be of class 'vector', 'Raster' or 'SpatRaster'.")
  obs <- unlist(obs)  # in case input was tibble
  pred <- unlist(pred)  # in case input was tibble

  if (any(na.omit(pred[]) < 0) || any(na.omit(pred[]) > 1))  # '[]' in case 'pred' is raster
    stop ("'pred' contains values that are not between 0 and 1.")

  vals <- na.omit(obs)
  if (is(model, "randomForest"))  vals <- as.integer(as.character(vals))
  if (!all(vals %in% c(0, 1))) stop("Favourability is only applicable when the response variable is binary, taking only values of 0 or 1.")

  if (!is.null(obs)) {
    n1 <- sum(obs == 1, na.rm = TRUE)
    n0 <- sum(obs == 0, na.rm = TRUE)
  }  else if (!is.null(n1n0)) {
    if (verbosity > 0) {
      if (!is.null(obs)) message("Argument 'n1n0' ignored in favour of 'obs'")
    }
    n1 <- n1n0[1]
    n0 <- n1n0[2]
    #if(n1 + n0 != length(pred)) stop("n1+n0 must equal the length of 'pred'.")
  } else if (!is.null(sample.preval)) {
    if (verbosity > 0) {
      if (!is.null(obs)) message("Argument 'sample.preval' ignored in favour of 'obs'")
      else if (!is.null(n1n0)) message("Argument 'sample.preval' ignored in favour of 'n1n0'")
    }
    n1 <- sample.preval * 100
    n0 <- 100 - n1
  } else stop("You need to provide either 'model', or 'pred' plus either one of 'obs', 'n1n0' or 'sample.preval'.")

  # slightly reduce probabilities of exactly 1, which would cause division by zero (resulting favourability is still 1):
  pred[pred == 1] <- 1 - 2.2e-16

  if (method == "RBV") {  # original Real, Barbosa & Vargas (2006)
    if (inv)
      out <- (pred * n1 / n0) / (pred * n1 / n0 + 1 - pred)
    else
      out <- (pred / (1 - pred)) / ((n1 / n0) + (pred / (1 - pred)))

  } else if (method == "AT") {  # Albert & Thuiller (2008); but see Acevedo & Real (2012)!
    sample.preval <- n1 / (n1 + n0)
    if (inv)
      out <- (pred * sample.preval) / (sample.preval + true.preval + pred * sample.preval)
    else
      out <- (pred / (1 - pred)) / ((sample.preval / true.preval) + (pred / (1 - pred)))

  } else stop("method must be either 'RBV' or 'AT'")

  if (is.vector(out) && !is.null(names(out)))  names(out) <- NULL  # in case input was tibble

  return(out)
}