R/MFreg.R

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

#' Calculate Multifunctionality Regularity (MFreg)
#'
#' @description
#' Calculates the multifunctionality regularity index, which measures how evenly
#' different ecosystem functions are distributed across the system. The function can
#' account for correlations between functions when specified.
#'
#' @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 in data). If NULL, equal weights are assigned.
#' @param cor Logical. If TRUE, function correlations are accounted for using redundancy correction. Default is FALSE.
#'
#' @details
#' Multifunctionality Regularity (MFreg) quantifies the evenness of function distribution
#' in an ecosystem. It is calculated as:
#'
#' \deqn{MFreg = \frac{-\sum_{i=1}^{n}\frac{w_i f_i}{\sum_{i=1}^{n}w_i f_i}\ln{\frac{w_i f_i}{\sum_{i=1}^{n}w_i f_i}}}{\ln(n)}}
#' where fi represents the normalized performance level of function i, wi is the weight of function i, and n is the total number of functions examined.
#'
#' 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 MFreg using these effective function values:
#'    \deqn{MFreg = \frac{-\sum_{i=1}^{n}\frac{w_i F_i}{\sum_{i=1}^{n}w_i F_i}\ln{\frac{w_i F_i}{\sum_{i=1}^{n}w_i F_i}}}{\ln(n)}}
#'
#' 5. The final result is the area under the curve (AUC) of MFreg values across different tau thresholds.
#'
#' @return A data frame with one column named "MFreg" containing the multifunctionality regularity
#' values for each observation (row) in the input data.
#'
#' @examples
#' data(forestfunctions)
#' head(forestfunctions)
#' MFreg(forestfunctions[,6:31], cor = FALSE)
#'
#' @export
MFreg <- function(data, weights = NULL, cor = FALSE) {
  # 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))
  MFreg_uncor <- function(wi,fi) {
    pi <- wi*fi / sum(wi*fi)
    wi <- wi[pi > 0]
    fi <- fi[pi > 0]
    pi <- pi[pi > 0]
    H <- -sum(pi * log(pi))
    S <- length(pi)
    e <- data.frame('MFreg' = H/log(S))
    return(e)
  }
  MFreg_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){
      pi <- wi*aivi$ai[,i] / sum(wi*aivi$ai[,i])
      H <- -sum(pi * log(pi))
      S <- length(pi)
      mfreg <- data.frame('MFreg' = H/log(S))
      data.frame(
        'MFreg' = mfreg,
        'tau' = tau[i]
      )
    }))
    tmp <- rbind(tmp,
                 data.frame(
                   'MFreg' = with(tmp, (sum(MFreg[seq_along(MFreg[-1])] * diff(tau)) +
                                          sum(MFreg[-1] * diff(tau))) / 2),
                   'tau' = 'AUC'
                 ))
    tmp_filtered <- subset(tmp, tau == 'AUC')
    tmp_selected <- tmp_filtered[, 1, drop = FALSE]
    return(tmp_selected)
  }
  if (!cor) {
    result <- do.call(rbind, lapply(1:nrow(data), function(i) MFreg_uncor(data[i,], weights)))
  }else{
    result <- do.call(rbind, lapply(1:nrow(data), function(i) {
      MFreg_cor(data[i,], weights, distM)
    }))
  }
  if (!is.null(rownames(data))) {
    rownames(result) <- rownames(data)
  }
  colnames(result) <- "MFreg"
  return(result)
}


# MFreg(forestfunctions[,6:31])