R/ps_ordinate.R

In phylospatial: Spatial Phylogenetic Analysis

#' Community phylogenetic ordination
#'
#' Perform an ordination that reduces a spatial phylogenetic data set into `k` dimensions, using one of
#' several alternative ordination algorithms.
#'
#' @param ps A `phylospatial` object. Unless `method = "pca"`, `ps` must have a non-null `dissim` component,
#'   generated by \link{ps_add_dissim}.
#' @param method Ordination method, either `"cmds"` (the default, classical MDS, implemented via `stats::cmdscale()`,
#'   `"nmds"` (nonmetric MDS, implemented via `vegan::metaMDS()`; this is slower but often preferred), or `"pca"`
#'   (principal component analysis, implemented via `stats::prcomp()`),.
#' @param k Positive integer giving the desired number of output dimensions; default is `3`.
#' @param spatial Logical indicating whether a spatial object (inherited from `ps`) should be returned.
#'    Default is TRUE.
#' @seealso For visualization using ordination onto RGB color space, see [ps_rgb()].
#' @return A matrix or spatial object with `k` variables.
#' @examples
#' ps <- ps_add_dissim(ps_simulate(50, 5, 5))
#' ord <- ps_ordinate(ps, method = "cmds", k = 4)
#' terra::plot(ord)
#'
#' @export
ps_ordinate <- function(ps, method = c("cmds", "nmds", "pca"), k = 3, spatial = TRUE){

      method <- match.arg(method)
      enforce_ps(ps)

      ord <- ps_ordinate_internal(ps, method = method, k = k)

      # expand to full extent
      ps_expand(ps, ord, spatial = spatial && !is.null(ps$spatial))
}

# Internal version returning occupied-only ordination matrix
ps_ordinate_internal <- function(ps, method, k) {
      if(method == "pca") {
            comm <- t(t(ps$comm) * ps$tree$edge.length)
            # Remove zero-variance columns to avoid prcomp errors
            col_var <- apply(comm, 2, var, na.rm = TRUE)
            comm <- comm[, col_var > 0, drop = FALSE]
            ord <- stats::prcomp(comm, scale. = TRUE)$x[, 1:k]
      } else {
            stopifnot("Input data set contains no `dissim` data, which is required unless `method = 'pca'`; add it first using `ps_add_dissim()`." = !is.null(ps$dissim))

            # dissim is already occupied-only, use directly
            d <- as.matrix(ps$dissim)

            # sites fully segregated by the 2 basal clades have Inf distance;
            # set distance to 2x max observed distance
            d[is.infinite(d)] <- max(d[!is.infinite(d)]) * 2

            # ordinate
            if(method == "cmds") ord <- stats::cmdscale(d, k = k)
            if(method == "nmds") ord <- vegan::metaMDS(stats::as.dist(d), k = k, trace = 0)$points
      }

      colnames(ord) <- paste0("d", 1:k)
      ord
}


#' Map phylospatial data onto RGB color bands
#'
#' Perform an ordination that reduces a spatial phylogenetic data set into three dimensions that can be
#' plotted as the RGB bands of color space to visualize spatial patterns of community phylogenetic composition.
#' This function is a wrapper around `ps_ordinate()`.
#'
#' @param ps A `phylospatial` object with a non-null `dissim` component, generated by \link{ps_add_dissim}.
#' @param method Ordination method, either "pca" (principal component analysis implemented via `stats::prcomp()`),
#'    "cmds" (classical MDS, implemented via `stats::cmdscale()`), or "nmds" (the default, nonmetric MDS,
#'    implemented via `vegan::metaMDS()`; this is slower but often preferred).
#' @param trans A function giving a transformation to apply to each dimension of the ordinated data.
#'    The default is the identity function. Specifying \code{rank} generates a more uniform color distribution.
#' @param spatial Logical indicating whether a spatial object (inherited from `ps`) should be returned.
#'    Default is TRUE.
#'
#' @return A matrix or spatial object with three variables containing RGB color values in the range 0-1.
#' @examples
#' ps <- ps_add_dissim(ps_simulate(50, 20, 20))
#' RGB <- ps_rgb(ps, method = "cmds")
#' terra::plotRGB(RGB * 255, smooth = FALSE)
#'
#' @export
ps_rgb <- function(ps, method = c("nmds", "cmds", "pca"), trans = identity, spatial = TRUE){

      method <- match.arg(method)
      enforce_ps(ps)

      # get occupied-only ordination directly
      ord <- ps_ordinate_internal(ps, method = method, k = 3)

      rgb <- apply(ord, 2, function(x){
            x <- trans(x)
            (x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))
      })

      colnames(rgb) <- c("r", "g", "b")
      ps_expand(ps, rgb, spatial = spatial && !is.null(ps$spatial))
}