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 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 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("nmds", "cmds", "pca"), k = 3, spatial = TRUE){

      method <- match.arg(method)
      enforce_ps(ps)
      stopifnot("Input data set contains no `dissim` data; add it first using `ps_add_dissim()`." = !is.null(ps$dissim))

      d <- as.matrix(ps$dissim)
      rownames(d) <- colnames(d) <- paste("cell", 1:ncol(d))
      a <- occupied(ps) # sites with taxa
      da <- d[a, a]

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

      # ordinate
      if(method == "cmds") ord <- stats::cmdscale(da, k = k)
      if(method == "nmds") ord <- vegan::metaMDS(da, k = k, trace = 0)$points
      if(method == "pca") ord <- stats::prcomp(da)$x[,1:k]

      # reinsert NA values
      b <- replicate(k, a)
      colnames(b) <- paste0("d", 1:k)
      b[a,] <- ord
      b[!a,] <- NA

      if(spatial & !is.null(ps$spatial)) b <- to_spatial(b, ps$spatial)
      return(b)
}


#' 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(moss())
#' 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){

      ord <- ps_ordinate(ps, method = method, k = 3, spatial = FALSE)

      a <- occupied(ps)
      rgb <- apply(ord, 2, function(x){
            x <- trans(x)
            x[!a] <- NA
            (x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))
      })

      colnames(rgb) <- c("r", "g", "b")
      if(spatial) rgb <- to_spatial(rgb, ps$spatial)
      rgb
}