Rao's quadratic entropy

README.md

Computing Rao's Quadratic entropy

Takes a sample-by-species abundance matrix and species-by-species functional distance matrix and returns values of Rao's quadratic entropy. Values are optionally converted into effective numbers of species through the transformation: 1/(1 - D).

EXAMPLE

# Create sample-by-species abundance matrix
abund <- matrix(sample(1:100, 100, replace = T), 10, 10)

colnames(abund) <- paste0("species", 1:10)

# Create species-by-species functional distance matrix
traits <- FD::gowdis(matrix(runif(100, 0, 1), 10, 10, dimnames = list(colnames(abund), colnames(abund))))

# Compute diversity
raoQ(abund, traits)

raoQ.R

#' `raoQ` an R function to compute Rao's Quadratic entropy

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

#' @param abund a sample-by-species abundance matrix
#' @param trait a species-by-trait matrix
#' @param Hill = TRUE, whether Hill numbers (effective numbers of species) should be returned
#' @param scale = FALSE, whether values of Q should be scaled by their maximum value
#' @param method whether the "default" formulation of Rao's Q should be used or the "divc" formulation used 
#' @return a vector of sample diversities
#' 
#' The "divc" formulation uses the formula in:
#' Champely, S. and Chessel, D. (2002) Measuring biological diversity using Euclidean metrics. Environmental and Ecological Statistics, 9, 167–177.

raoQ <- function(abund, trait, Hill = TRUE, scale = FALSE, method = "default") {
  
  abund <- as.matrix(abund)
  
  anames <- colnames(abund)[order(colnames(abund))]
  
  abund <- abund[, anames, drop = FALSE]
  
  trait <- as.matrix(trait)
  
  trait <- trait[anames, anames]
  
  if(ncol(abund) != ncol(trait))
    
    stop("Not all species in the abundance matrix appear in the trait matrix!") 
  
  abund <- abund / rowSums(abund)
  
  if(method == "default") 
    
    Q <- apply(abund, 1, function(x) crossprod(x, trait %*% x))
  
  if(method == "divc")
    
    Q <- apply(abund, 1, function(x) x %*% trait^2 %*% (x/2/sum(x)^2))
  
  if(Hill == TRUE) Q <- 1/(1 - Q)
  
  if(scale == TRUE) Q <- Q / max(Q)
  
  names(Q) <- rownames(abund)
  
  return(Q)
  
}