Skip to content

Instantly share code, notes, and snippets.

@jslefche
Last active May 17, 2021 05:58
Show Gist options
  • Save jslefche/09756ff84afc7b6a82ea0582e663d098 to your computer and use it in GitHub Desktop.
Save jslefche/09756ff84afc7b6a82ea0582e663d098 to your computer and use it in GitHub Desktop.
Rao's quadratic entropy

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` 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)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment