R/MFdiv.R

In emf: Ecosystem Multifunctionality: Richness, Divergence, and Regularity

#' Calculate Multifunctionality Divergence (MFdiv)
#'
#' @description
#' Multifunctionality divergence (MFdiv) was calculated to quantify the degree of difference between different functions of an ecosystem.
#' MFdiv was calculated by the mean pairwise distance method.
#'
#' @param data A data frame or matrix where rows represent observations and columns represent functions.
#' @param weights A numeric vector of weights for each function (column). If NULL, equal weights are assigned.
#' @param cor Logical. If TRUE, calculates MFdiv with redundancy correction based on correlation between functions.
#'            If FALSE, calculates uncorrected MFdiv.
#'
#' @details
#' To measure MFdiv quantitatively, we employ the mean pairwise distance method (Webb et al., 2002).
#'
#' For uncorrected MFdiv, the formula is:
#' \deqn{MFdiv = \frac{\sum_{i=1}^{n}\sum_{j=i+1}^{n}{w_i w_j D_{ij}}}{\sum_{i=1}^{n}\sum_{j=i+1}^{n}{w_i w_j}}}{MFdiv = (sum of weighted distances)/(sum of weights)}
#' where \eqn{D_{ij} = |f_i - f_j|}
#'
#' When redundancy correction is applied (`cor = TRUE`), the function accounts for correlations
#' between ecosystem functions. The correction process involves:
#'
#' 1. Calculating a distance matrix based on correlations: \eqn{d_{ij} = \sqrt{1 - |r_{ij}|}}
#'
#' 2. Applying threshold-based correction: \eqn{d_{ij}(\tau) = \min(d_{ij}, \tau)}
#'
#' 3. Computing effective function values:
#'    \eqn{F_i(\tau) = \sum_{j=1}^{L}(1 - \frac{d_{ij}(\tau)}{\tau})f_j}
#'
#' 4. Calculating the corrected MFdiv using these effective function values:
#'    \deqn{D_{ij}(\tau) = |F_i - F_j|}
#'    \deqn{MFdiv = \frac{\sum_{i=1}^{n}\sum_{j=i+1}^{n}{w_i w_j D_{ij}(\tau)}}{\sum_{i=1}^{n}\sum_{j=i+1}^{n}{w_i w_j}}}{MFdiv = (sum of weighted distances)/(sum of weights)}
#'
#' 5. The final result is the area under the curve (AUC) of MFdiv values across different tau thresholds.
#'
#' @return A data frame with MFdiv values for each observation (row) in the input data.
#'
#' @examples
#' data(forestfunctions)
#' head(forestfunctions)
#' MFdiv(forestfunctions[,6:31], cor = FALSE)
#'
#' @export
MFdiv <- function(data, weights = NULL, cor = FALSE) {
  if (ncol(data) <= 2) {
    stop("The number of functions must be greater than 2!")
  }
  # If no weights are provided, create a weight vector with all 1's
  if (is.null(weights)) {
    weights <- rep(1, ncol(data))
  }
  if (length(weights) != ncol(data)) {
    stop("The length of the weight vector must be equal to the number of columns in the data frame")
  }
  correlation <- cor(data)
  distM <- sqrt(1 - abs(correlation))
  MFdiv_uncor <- function(fi, wi) {
    wi <- wi[fi > 0]
    fi <- fi[fi > 0]
    n <- length(fi)
    numerator <- 0
    denominator <- 0
    for (i in 1:(n-1)) {
      for (j in (i+1):n) {
        d_ij <- abs(fi[i] - fi[j])
        numerator <- numerator + wi[i] * wi[j] * d_ij
        denominator <- denominator + wi[i] * wi[j]
      }
    }
    data.frame('MFdiv' = numerator / denominator)
  }
  MFdiv_cor <- function(fi, wi, distM, tau = seq(0, 1, 0.01)) {
    distM <- distM[fi > 0, fi > 0]
    wi <- wi[fi > 0]
    fi <- fi[fi > 0]
    data_transform <- function (fi, dij, tau) {
      out <- lapply(tau, function(tau_) {
        dij_ <- dij
        if (tau_ == 0) {
          dij_[dij_ > 0] <- 1
          a <- as.vector((1 - dij_/1) %*% fi)
        } else {
          dij_[which(dij_ > tau_, arr.ind = T)] <- tau_
          a <- as.vector((1 - dij_/tau_) %*% fi)
        }
        v <- fi/a
        v[a == 0] = 1
        cbind(a, v)
      })
      out_a <- matrix(sapply(out, function(x) x[, 1]), ncol = length(tau))
      out_v <- matrix(sapply(out, function(x) x[, 2]), ncol = length(tau))
      colnames(out_a) <- colnames(out_v) <- paste0("tau_", round(tau, 3))
      list(ai = out_a, vi = out_v)
    }
    aivi <- data_transform(fi, distM, tau)
    tmp <- do.call(rbind, lapply(1:length(tau), function(i){
      Fi <- aivi$ai[, i]
      n <- length(Fi)
      numerator <- 0
      denominator <- 0
      for (j in 1:(n-1)) {
        for (k in (j+1):n) {
          d_ij <- abs(Fi[j] - Fi[k])
          numerator <- numerator + wi[j] * wi[k] * d_ij
          denominator <- denominator + wi[j] * wi[k]
        }
      }
      mfdiv <- numerator / denominator
      data.frame(
        'MFdiv' = mfdiv,
        'tau' = tau[i]
      )}
    ))
    tmp <- rbind(tmp,
                 data.frame(
                   'MFdiv' = with(tmp, (sum(MFdiv[seq_along(MFdiv[-1])] * diff(tau)) +
                                          sum(MFdiv[-1] * diff(tau))) / 2),
                   'tau' = 'AUC'
                 ))
    tmp <- subset(tmp, tau == 'AUC', select = 1)
    return(tmp)
  }
  if (!cor) {
    result <- do.call(rbind, lapply(1:nrow(data), function(i) MFdiv_uncor(data[i,], weights)))
  }else{
    result <- do.call(rbind, lapply(1:nrow(data), function(i) {
      MFdiv_cor(data[i,], weights, distM)
    }))
  }
  if (!is.null(rownames(data))) {
    rownames(result) <- rownames(data)
  }
  colnames(result) <- "MFdiv"
  return(result)
}


# MFdiv(forestfunctions[,6:31])