R/ps_prioritize.R

In phylospatial: Spatial Phylogenetic Analysis

#' Calculate taxon conservation benefit
#'
#' Nonlinear function that converts proportion of range conserved into conservation "benefit."
#'
#' @param x Fraction of taxon range protected (value between 0 and 1).
#' @param lambda Shape parameter.
#'
#' @return Value between 0 and 1.
#' @export
benefit <- function(x, lambda = 1){
      lambda <- 2^lambda
      (1-(1-pmin(x, 1))^lambda)^(1/lambda)
      # (pmin needed in cases where numerical rounding results in values slightly >1)
}


#' Plot alternative lambda values
#'
#' Show a plot illustrating alternative values for the `lambda` parameter in \link{ps_prioritize}. Lambda determines the shape of
#' the "benefit" function that determines the conservation value of protecting a given proportion of the geographic range of a
#' species or clade. Positive values place a higher priority on protecting additional populations of largely unprotected taxa,
#' whereas negative values place a higher priority on protecting additional populations of relatively well-protected taxa. The
#' default value used by \link{ps_prioritize} is 1.
#'
#' @param lambda A vector of lambda values to plot
#' @return Plots a figure
#' @examples
#' plot_lambda()
#' plot_lambda(seq(0, 3, .1))
#'
#' @export
plot_lambda <- function(lambda = c(-1, -.5, 0, .5, 2, 1)){
      x <- seq(0, 1, .01)
      d <- sapply(lambda, function(l) benefit(x, lambda = l))
      graphics::matplot(x, d, type = "l", lty = 1, col = 1:ncol(d),
                        xlab = "proportion of taxon range protected",
                        ylab = "conservation benefit",
                        main = "Examples of how `labmda` affects\nthe conservation benefit function")
      graphics::legend("bottomright", legend = lambda,
                       col = 1:ncol(d), pch = 1, title = "lambda")
}


#' Phylogenetic conservation prioritization
#'
#' Create a ranking of conservation priorities using optimal or probabilistic forward stepwise selection. Prioritization accounts for
#' the occurrence quantities for all lineages present in the site, including terminal taxa and larger clades; the evolutionary branch
#' lengths of these lineages on the phylogeny, which represent their unique evolutionary heritage; the impact that protecting the site
#' would have on these lineages' range-wide protection levels; the compositional complementarity between the site, other high-priority
#' sites, and existing protected areas; the site's initial protection level; the relative cost of protecting the site; and a free
#' parameter "lambda" determining the shape of the conservation benefit function.
#'
#' @param ps phylospatial object.
#' @param init Optional numeric vector or spatial object giving the starting protection status of each site across the study area.
#'    Values should be between 0 and 1 and represent the existing level of conservation effectiveness in each site. If this argument
#'    is not specified, it is assumed that no existing reserves are present.
#' @param cost Optional numeric vector or spatial object giving the relative cost of protecting each site. Values should be positive,
#'    with greater values indicating higher cost of conserving a site. If this argument is not specified, cost is assumed to be uniform
#'    across sites.
#' @param lambda Shape parameter for taxon conservation benefit function. This can be any real number. Positive values, such as the default
#'    value \code{1}, place higher priority on conserving the first part of the range of a given species or clade, while negative values
#'    (which are not typically used) place higher priority on fully protecting the most important taxa (those with small ranges and long branches)
#'    rather than partially protecting all taxa. See the function \link{plot_lambda} for an illustration of alternative `lambda` values.
#' @param protection Degree of protection of proposed new reserves (number between 0 and 1, with same meaning as \code{init}).
#' @param max_iter Integer giving max number of iterations to perform before stopping, i.e. max number of sites to rank.
#' @param method Procedure for selecting which site to add to the reserve network at each iteration:
#'  \itemize{
#'    \item "optimal": The default, this selects the site with the highest marginal value at each iteration. This is a
#'    optimal approach that gives the same result each time.
#'    \item "probable": This option selects a site randomly, with selection probabilities calculated as a function of sites' marginal values. This
#'    approach gives a different prioritization ranking each time an optimization is performed, so \code{n_reps} optimizations are performed,
#'    and ranks for each site are summarized across repetitions.
#' }
#' @param trans A function that transforms marginal values into relative selection probabilities; only used if `method = "probable"`.
#'    The function should take a vector of positive numbers representing marginal values and return an equal-length vector of positive numbers
#'    representing a site's relative likelihood of being selected. The default function returns the marginal value if a site is in the top 25
#'    highest-value sites, and zero otherwise.
#' @param n_reps Number of random repetitions to do; only used if `method = "probable"`. Depending on the data set, a large number of reps
#'    (more than the default of 100) may be needed in order to achieve a stable result. This may be a computational barrier for large data
#'    sets; multicore processing via \code{n_cores} can help.
#' @param n_cores Number of compute cores to use for parallel processing; only used if `method = "probable"`.
#' @param summarize Logical: should summary statistics across reps (TRUE, default) or the reps themselves (FALSE) be returned? Only relevant
#'    if `method = "probable"`.
#' @param spatial Logical: should the function return a spatial object (TRUE, default) or a matrix (FALSE)?
#' @param progress Logical: should a progress bar be displayed?
#' @details This function uses the forward stepwise selection algorithm of Kling et al. (2019) to generate a ranked conservation prioritization.
#'    Prioritization begins with the starting protected lands network identified in `init`, if provided. At each iteration, the marginal
#'    conservation value of fully protecting each site is calculated, and a site is selected to be added to the reserve network. Selection can
#'    happen either in an "optimal" or "probable" fashion as described under the `method` argument. This process is repeated until all sites
#'    are fully protected or until \code{max_iter} has been reached, with sites selected early in the process considered higher conservation
#'    priorities.
#'
#'    The benefit of the probabilistic approach is that it relaxes the potentially unrealistic assumption that protected land will actually be
#'    added in the optimal order. Since the algorithm avoids compositional redundancy between high-priority sites, the optimal approach will
#'    never place high priority on a site that has high marginal value but is redundant with a slightly higher-value site, whereas the
#'    probabilistic approach will select them at similar frequencies (though never in the same randomized run).
#'
#'    Every time a new site is protected as the algorithm progresses, it changes the marginal conservation value of the other sites. Marginal
#'    value is the increase in conservation benefit that would arise from fully protecting a given site, divided by the cost of protecting the
#'    site. This is calculated as a function of the site's current protection level, the quantitative presence probability or abundance of all
#'    terminal taxa and larger clades present in the site, their evolutionary branch lengths on the phylogeny, the impact that protecting the
#'    site would have on their range-wide protection levels, and the free parameter `lambda`. `lambda` determines the relative importance of
#'    protecting a small portion of every taxon's range, versus fully protecting the ranges of more valuable taxa (those with longer
#'    evolutionary branches and smaller geographic ranges).
#' @seealso [benefit()], [plot_lambda()]
#' @references Kling, M. M., Mishler, B. D., Thornhill, A. H., Baldwin, B. G., & Ackerly, D. D. (2019). Facets of phylodiversity: evolutionary
#'    diversification, divergence and survival as conservation targets. Philosophical Transactions of the Royal Society B, 374(1763), 20170397.
#' @return Matrix or spatial object containing a ranking of conservation priorities. Lower rank values represent higher
#'    conservation priorities. All sites with a lower priority than \code{max_iter} have a rank value equal to the number
#'    of sites in the input data set (i.e. the lowest possible priority).
#'  \describe{
#'    \item{If `method = "optimal"`. }{the result contains a single variable "priority" containing the ranking.}
#'    \item{If `method = "probable"` and `summarize = TRUE`, }{the "priority" variable gives the average rank across reps,
#'    variables labeled "pctX" give the Xth percentile of the rank distribution for each site, variables labeled "topX"
#'    give the proportion of reps in which a site was in the top X highest-priority sites, and variables labeled "treX" give
#'    a ratio representing "topX" relative to the null expectation of how often "topX" should occur by chance alone.}
#'    \item{If `method = "probable"` and `summarize = FALSE`, }{the result contains the full set of \code{n_rep} solutions,
#'    each representing the the ranking, with low values representing higher priorities.. }
#' }
#' @examples
#' \donttest{
#' # simulate a toy `phylospatial` data set
#' set.seed(123)
#' ps <- ps_simulate()
#'
#' # basic prioritization
#' p <- ps_prioritize(ps)
#'
#' # specifying locations of initial protected areas
#' # (can be binary, or can be continuous values between 0 and 1)
#' # here we'll create an `init` raster with arbitrary values ranging from 0-1,
#' # using the reference raster layer that's part of our `phylospatial` object
#' protected <- terra::setValues(ps$spatial, seq(0, 1, length.out = 400))
#' cost <- terra::setValues(ps$spatial, rep(seq(100, 20, length.out = 20), 20))
#' p <- ps_prioritize(ps, init = protected, cost = cost)
#'
#' # using probabilistic prioritization
#' p <- ps_prioritize(ps, init = protected, cost = cost,
#'       method = "prob", n_reps = 1000, max_iter = 10)
#' terra::plot(p$top10)
#' }
#' @export
ps_prioritize <- function(ps,
                          init = NULL,
                          cost = NULL,
                          lambda = 1,
                          protection = 1,
                          max_iter = NULL,
                          method = c("optimal", "probable"),
                          trans = function(x) replace(x, which(rank(-x) > 25), 0),
                          n_reps = 100, n_cores = 1, summarize = TRUE,
                          spatial = TRUE,
                          progress = interactive()){

      method <- match.arg(method)
      enforce_ps(ps)
      if(lambda < 0) warning("choosing a negative value for `lambda` is not generally recommended.")

      e <- ps$tree$edge.length / sum(ps$tree$edge.length) # edges evolutionary value

      occ <- ps$occupied
      n_occ <- nrow(ps$comm)

      # extract init/cost values for occupied sites only
      if(is.null(init)){
            p <- rep(0, n_occ)
      }else{
            p_full <- init[]
            p <- p_full[occ]
            stopifnot("`init` may not contain NA values, or values outside the 0-1 range, for sites that contain taxa." =
                            all(is.finite(p)) & min(p) >= 0 & max(p) <= 1)
      }

      if(is.null(cost)){
            cost <- rep(1, n_occ)
      }else{
            cost_full <- cost[]
            cost <- cost_full[occ]
            stopifnot("`cost` may only contain finite, nonnegative values for sites that contain taxa." =
                            all(is.finite(cost)) & min(cost) >= 0)
      }

      n_ranks <- n_occ
      n_poss <- sum(p < 1)

      y <- rep(NA_real_, n_occ)  # prioritization rankings (occupied sites only)
      r <- rep(n_ranks, n_ranks)

      m <- apply(ps$comm, 2, function(x) x / sum(x, na.rm = TRUE)) # normalize to fraction of range

      n_iter <- ifelse(is.null(max_iter), n_ranks, min(max_iter, n_ranks))

      # Pre-compute lambda transformation for benefit function
      lambda_t <- 2^lambda
      inv_lambda_t <- 1 / lambda_t

      # Inline benefit function for speed
      benefit_fast <- function(x) {
            (1 - (1 - pmin(x, 1))^lambda_t)^inv_lambda_t
      }

      ranks <- function(y, progress = TRUE){
            p_local <- p
            r_local <- r

            if(progress) pb <- utils::txtProgressBar(min = 0, max = n_iter,
                                                     initial = 0, style = 3)
            for(i in 1:n_iter){
                  if(progress) utils::setTxtProgressBar(pb, i)

                  # Range protection: vectorized with colSums
                  b <- colSums(m * p_local)

                  # Current network value
                  ben_b <- benefit_fast(b)
                  v <- sum(e * ben_b)

                  # Marginal protection boost
                  delta_p <- protection - p_local
                  delta_p[delta_p <= 0] <- NA

                  # Marginal value per site: the gain from protecting that site
                  mv <- sapply(1:n_occ, function(j){
                        if(is.na(delta_p[j])) return(NA)
                        b_new <- b + m[j,] * delta_p[j]
                        (sum(e * benefit_fast(b_new)) - v) / cost[j]
                  })

                  if(all(is.na(mv))) break

                  if(method == "optimal"){
                        sel <- which.max(mv)
                  }else{
                        probs <- trans(mv)
                        probs[is.na(probs)] <- 0
                        if(sum(probs) == 0) break
                        sel <- sample(length(probs), 1, prob = probs)
                  }

                  p_local[sel] <- protection
                  y[sel] <- i
                  r_local[sel] <- i
            }
            if(progress) close(pb)
            return(y)
      }

      if(method == "optimal"){
            y <- ranks(y, progress = progress)
            result <- matrix(y, ncol = 1)
            colnames(result) <- "priority"
      }else{
            # probabilistic: multiple reps
            all_ranks <- matrix(NA, n_occ, n_reps)
            if(n_cores == 1){
                  for(rep in 1:n_reps){
                        all_ranks[, rep] <- ranks(y, progress = FALSE)
                  }
            }else{
                  if (!requireNamespace("furrr", quietly = TRUE)) {
                        stop("To use `n_cores` greater than 1, package `furrr` must be installed.", call. = FALSE)
                  }
                  plan <- future::plan()
                  future::plan(future::multisession, workers = n_cores)
                  rnd <- furrr::future_map(
                        1:n_reps,
                        function(i) ranks(y, progress = FALSE),
                        .progress = progress,
                        .options = furrr::furrr_options(seed = TRUE)
                  )
                  future::plan(plan)
                  for(i in 1:n_reps) all_ranks[, i] <- rnd[[i]]
            }

            if(summarize){
                  avg_rank <- rowMeans(all_ranks, na.rm = TRUE)
                  result <- matrix(avg_rank, ncol = 1)
                  colnames(result) <- "priority"

                  # Add summary columns
                  for(pct in c(5, 25, 50, 75, 95)){
                        q <- apply(all_ranks, 1, stats::quantile, probs = pct/100, na.rm = TRUE)
                        result <- cbind(result, q)
                        colnames(result)[ncol(result)] <- paste0("pct", pct)
                  }
                  for(top_n in c(10, 25, 50)){
                        top_prop <- rowMeans(all_ranks <= top_n, na.rm = TRUE)
                        tre <- top_prop / (top_n / n_occ)
                        result <- cbind(result, top_prop, tre)
                        colnames(result)[(ncol(result)-1):ncol(result)] <- c(paste0("top", top_n), paste0("tre", top_n))
                  }
            }else{
                  result <- all_ranks
            }
      }

      # expand to full extent
      if(spatial & !is.null(ps$spatial)){
            result <- ps_expand(ps, result, spatial = TRUE)
      } else {
            result <- ps_expand(ps, result, spatial = FALSE)
      }
      result
}