R/MFric.R

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

#' Calculate Multifunctionality Richness (MFric)
#'
#' @description
#' This function calculates the Multifunctionality Richness (MFric) for each row in a dataset.
#' MFric represents the average level of multiple ecosystem function indicators, reflecting
#' the overall performance of an ecosystem across various functional metrics.
#'
#' @param data A numeric data frame or matrix where rows represent observations (e.g., sites, plots)
#'   and columns represent different ecosystem functions.
#' @param weights A numeric vector of weights for each function (column) in the data.
#'   If NULL (default), equal weights of 1 are assigned to all functions.
#' @param cor Logical. If FALSE (default), calculates uncorrected MFric.
#'   If TRUE, calculates correlation-corrected MFric accounting for redundancy among functions.
#'
#' @return A data frame with a single column named "MFric" containing the calculated
#'   Multifunctionality Richness values for each row in the input data.
#'   Row names are preserved from the input data if available.
#'
#' @details
#' The uncorrected MFric is calculated as:
#' \deqn{MFric = \frac{\sum_{i=1}^{n} w_i f_i}{\sum_{i=1}^{n} w_i}}
#' where fi represents the normalized performance level of function i, and wi denotes the weight assigned to function i.
#'
#' 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 MFric using these effective function values:
#'    \deqn{MFric = \frac{\sum_{i=1}^{n} w_i F_i}{\sum_{i=1}^{n} w_i}}
#'
#' 5. The final result is the area under the curve (AUC) of MFric values across different tau thresholds.
#'
#' @examples
#' data(forestfunctions)
#' head(forestfunctions)
#' MFric(forestfunctions[,6:31], cor = FALSE)
#'
#' @export
MFric <- 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))
  MFric_uncor <- function(fi, wi) {
    wi <- wi[fi > 0]
    fi <- fi[fi > 0]
    data.frame('MFric' = sum(wi * fi) / sum(wi))
  }
  MFric_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)
      data.frame(
        'MFric' = sum(wi * aivi$ai[,i]) / sum(wi),
        'tau' = tau[i]
      )
    ))
    tmp <- rbind(tmp,
                 data.frame(
                   'MFric' = with(tmp, (sum(MFric[seq_along(MFric[-1])] * diff(tau)) +
                                          sum(MFric[-1] * diff(tau))) / 2),
                   'tau' = 'AUC'
                 ))
    tmp <- subset(tmp, tau == 'AUC')
    tmp <- tmp[, 1, drop = FALSE]
    return(tmp)
  }
  if (!cor) {
    result <- do.call(rbind, lapply(1:nrow(data), function(i) MFric_uncor(data[i,], weights)))
  }else{
    result <- do.call(rbind, lapply(1:nrow(data), function(i) {
      MFric_cor(data[i,], weights, distM)
    }))
  }
  if (!is.null(rownames(data))) {
    rownames(result) <- rownames(data)
  }
  colnames(result) <- "MFric"
  return(result)
}

# data(forestfunctions)
# head(forestfunctions)
# MFric(forestfunctions[,6:31])