Compute ultrametric trait matrix for functional traits

README.md

Computing ultrametric distance matrices from species' functional traits

fdist takes a functional trait matrix, and returns an ultrametric distance matrix using the method that best preserves the original (non-ultrametric) distances.

References:

Mouchet, M., Guilhaumon, F., Villéger, S., Mason, N. W., Tomasini, J. A., & Mouillot, D. (2008). Towards a consensus for calculating dendrogram‐based functional diversity indices. Oikos, 117(5), 794-800.

Mérigot, B., Durbec, J. P., & Gaertner, J. C. (2010). On goodness-of-fit measure for dendrogram-based analyses. Ecology, 91(6), 1850-1859.

Example

# Create example dataset
example.data <- matrix(runif(100, 0, 1), nrow = 10, ncol = 10)

# Get distance matrix from example data
example.dist <- dist(example.data)

# Load required libraries
library(FD)
library(clue)

# Run function gow2dis
fdist(example.dist, consensus = TRUE)

# Run on raw trait matrix
fdist(example.data)

# Replace continuous with categorical trait
example.data[, 8] <- sample(letters[1:4], replace = T)

fdist(example.data)

fdist.R

#' `fdist` an R function to compute the ultrametric distance based on species' functional traits

#' Author: Jon Lefcheck
#' Last updated: 13 September 2017

#' REQUIRES: `FD` and `clue`

#' @param d a species-by-trait (distance) matrix
#' @param w an optional vector of weights
#' @param consensus = TRUE whether the consensus tree should be evaluated
#' @param ... additional arguments to `gowdis`
#' @return a species-by-species distance matrix of class 'dist'

fdist <- function(d, w, consensus = TRUE, ...) {
  
  if(any(class(d) != "dist")) {
    
    if(missing(w)) w <- rep(1, ncol(d))
    
    if(!any(apply(d, 1, is.numeric))) {
      
      d <- FD::gowdis(as.data.frame(d), w = w, ...) 
      
      class(d) <- c(class(d), "gowdis")
      
      message("Some traits are not numeric: computing Gower distances.")
      
    } else {
      
      d <- sweep(d, 2, w, function(x, y) scale(x * sqrt(y)))
      
      d <- dist(d)
      
    }
    
  }
  
  tree_methods <- c("single", "complete", "average", "mcquitty", "ward.D")
  
  trees <- lapply(tree_methods, function(j) hclust(d, method = j))
  
  trees.ultra <- lapply(trees,  function(j) clue::cl_ultrametric(as.hclust(j)))
  
  names(trees.ultra) <- tree_methods
  
  if(consensus == TRUE) {
    
    ensemble.trees <- clue::cl_ensemble(list = trees)
    
    consensus.tree <- clue::cl_consensus(ensemble.trees) 
    
    trees.ultra <- c(trees.ultra, consensus.tree[1])
    
    names(trees.ultra)[length(trees.ultra)] <- "consensus"
    
  }
  
  trees.dissim <- lapply(trees.ultra,  function(j) clue::cl_dissimilarity(j, d, method = "spectral"))
  
  trees.dissim <- do.call(rbind, trees.dissim)
  
  min.tree <- which.min(trees.dissim)
  
  message(paste("The best tree was:", names(trees.ultra)[min.tree]))
  
  finald <- trees.ultra[names(trees.ultra) == names(trees.ultra)[min.tree]][[1]]
  
  return(finald)
  
}