R/indexing.R

In SQIpro: Comprehensive Soil Quality Index Computation and Visualization

# ============================================================
# SQIpro: Indexing Functions
# Six methods: Linear, Regression, PCA, Fuzzy, Entropy, TOPSIS
# ============================================================

# Internal helper ─────────────────────────────────────────────
.min_max <- function(x) {
  r <- range(x, na.rm = TRUE)
  if (!is.finite(r[1]) || !is.finite(r[2])) return(rep(0.5, length(x)))
  if (r[1] == r[2]) return(rep(0.5, length(x)))
  (x - r[1]) / (r[2] - r[1])
}

.group_means <- function(scored, group_cols, vars) {
  syms <- rlang::syms(group_cols)
  scored %>%
    dplyr::group_by(!!!syms) %>%
    dplyr::summarise(dplyr::across(dplyr::all_of(vars), mean, na.rm = TRUE),
                     .groups = "drop")
}


#' Soil Quality Index: Linear Scoring Method
#'
#' @description
#' Computes the Soil Quality Index (SQI) using the linear additive scoring
#' method of Doran & Parkin (1994) and Andrews et al. (2004).  Each variable
#' score (0--1, from \code{\link{score_all}}) is averaged across replicates
#' within each group, optionally weighted, and then min-max normalised to
#' produce the final index.
#'
#' \deqn{SQI_g = \frac{\sum_{j=1}^{p} w_j \bar{S}_{gj}}{\sum_{j=1}^{p} w_j}}
#'
#' where \eqn{\bar{S}_{gj}} is the mean score of variable \eqn{j} in group
#' \eqn{g} and \eqn{w_j} is the weight of variable \eqn{j}.
#'
#' @param scored A scored data frame from \code{\link{score_all}}.
#' @param config A \code{sqi_config} object (see \code{\link{make_config}}).
#' @param group_cols Character vector of grouping column names.
#' @param mds_vars Character vector. If supplied, only these variables are
#'   used.  Otherwise all numeric variables in \code{config} are used.
#' @param weights Named numeric vector of variable weights.  Defaults to
#'   equal weights (1 for all).  Names must match variable names.
#'
#' @return A data frame with group columns plus:
#'   \describe{
#'     \item{SQI_linear}{Final normalised Soil Quality Index (0--1).}
#'     \item{Raw_score}{Weighted mean score before normalisation.}
#'   }
#'
#' @references
#' Doran, J.W., & Parkin, T.B. (1994). Defining and assessing soil quality.
#' In J.W. Doran et al. (Eds.), \emph{Defining Soil Quality for a
#' Sustainable Environment}, pp. 1--21. SSSA Special Publication 35.
#' \doi{10.2136/sssaspecpub35.c1}
#'
#' Andrews, S.S., Karlen, D.L., & Cambardella, C.A. (2004). The soil
#' management assessment framework. \emph{Soil Science Society of America
#' Journal}, 68(6), 1945--1962. \doi{10.2136/sssaj2004.1945}
#'
#' @examples
#' data(soil_data)
#' cfg <- make_config(
#'   variable = c("pH","EC","BD","OC","MBC","Clay"),
#'   type     = c("opt","less","less","more","more","opt"),
#'   opt_low  = c(6.0, NA, NA, NA, NA, 20),
#'   opt_high = c(7.0, NA, NA, NA, NA, 35)
#' )
#' scored <- score_all(soil_data, cfg, group_cols = c("LandUse","Depth"))
#' result <- sqi_linear(scored, cfg, group_cols = c("LandUse","Depth"))
#' print(result)
#'
#' @export
sqi_linear <- function(scored, config, group_cols = "LandUse",
                       mds_vars = NULL, weights = NULL) {

  vars <- if (!is.null(mds_vars)) mds_vars else config$variable
  vars <- intersect(vars, names(scored))

  gm <- .group_means(scored, group_cols, vars)

  # Weights
  w <- if (!is.null(weights)) {
    weights[vars]
  } else if ("weight" %in% names(config)) {
    idx <- match(vars, config$variable)
    stats::setNames(ifelse(is.na(idx), 1, config$weight[idx]), vars)
  } else {
    stats::setNames(rep(1, length(vars)), vars)
  }
  w[is.na(w)] <- 1

  score_mat  <- as.matrix(gm[, vars])
  raw_score  <- as.numeric(score_mat %*% w / sum(w))
  sqi_index  <- .min_max(raw_score)

  result <- gm[, group_cols, drop = FALSE]
  result$Raw_score  <- round(raw_score, 4)
  result$SQI_linear <- round(sqi_index, 4)
  result[order(-result$SQI_linear), ]
}


#' Soil Quality Index: Regression-Based Method
#'
#' @description
#' Computes the SQI using stepwise multiple linear regression to identify
#' and weight the most predictive soil variables.  The dependent variable
#' (e.g., crop yield, total biomass) determines which variables enter the
#' model.  Regression coefficients serve as weights in the index.
#'
#' This follows the method described by Masto et al. (2008) and
#' Mukherjee & Lal (2014).
#'
#' @param scored A scored data frame from \code{\link{score_all}}.
#' @param config A \code{sqi_config} object.
#' @param dep_var Character. Name of the dependent variable column in
#'   \code{scored}.
#' @param group_cols Character vector of grouping column names.
#' @param mds_vars Character vector of candidate predictor variable names.
#'   If \code{NULL}, all variables in \code{config} are used.
#' @param direction Character. Direction for stepwise selection:
#'   \code{"both"} (default), \code{"forward"}, or \code{"backward"}.
#'
#' @return A data frame with group columns plus:
#'   \describe{
#'     \item{SQI_regression}{Normalised SQI (0--1).}
#'     \item{selected_vars}{(attribute) Character vector of selected
#'       predictors.}
#'   }
#'
#' @references
#' Masto, R.E., Chhonkar, P.K., Singh, D., & Patra, A.K. (2008).
#' Alternative soil quality indices. \emph{Environmental Monitoring and
#' Assessment}, 136, 419--435. \doi{10.1007/s10661-007-9697-z}
#'
#' Mukherjee, A., & Lal, R. (2014). Comparison of soil quality index using
#' three methods. \emph{PLOS ONE}, 9(8), e105981.
#' \doi{10.1371/journal.pone.0105981}
#'
#' @examples
#' data(soil_data)
#' cfg <- make_config(
#'   variable = c("pH","EC","BD","OC","MBC","Clay"),
#'   type     = c("opt","less","less","more","more","opt"),
#'   opt_low  = c(6.0, NA, NA, NA, NA, 20),
#'   opt_high = c(7.0, NA, NA, NA, NA, 35)
#' )
#' scored <- score_all(soil_data, cfg, group_cols = c("LandUse","Depth"))
#' # OC used as surrogate dependent variable
#' result <- sqi_regression(scored, cfg, dep_var = "OC",
#'                          group_cols = c("LandUse","Depth"))
#' print(result)
#'
#' @export
sqi_regression <- function(scored, config, dep_var,
                            group_cols = "LandUse",
                            mds_vars   = NULL,
                            direction  = "both") {

  if (!dep_var %in% names(scored))
    stop("`dep_var` '", dep_var, "' not found in scored data.", call. = FALSE)

  vars <- if (!is.null(mds_vars)) mds_vars else config$variable
  vars <- intersect(vars, names(scored))
  vars <- setdiff(vars, dep_var)

  num_data           <- scored[, c(dep_var, vars)]
  names(num_data)[1] <- "Dep"
  num_data           <- stats::na.omit(num_data)

  full_model <- stats::lm(Dep ~ ., data = num_data)
  best_model <- stats::step(full_model, direction = direction, trace = 0)
  coefs      <- stats::coef(best_model)
  sel_vars   <- names(coefs)[names(coefs) != "(Intercept)"]

  gm      <- .group_means(scored, group_cols, sel_vars)
  score_m <- as.matrix(gm[, sel_vars, drop = FALSE])
  w       <- coefs[sel_vars]
  raw     <- as.numeric(score_m %*% w / length(sel_vars))

  result               <- gm[, group_cols, drop = FALSE]
  result$SQI_regression <- round(.min_max(raw), 4)
  attr(result, "selected_vars") <- sel_vars
  result[order(-result$SQI_regression), ]
}


#' Soil Quality Index: PCA-Based Method
#'
#' @description
#' Computes SQI using Principal Component Analysis, weighting selected MDS
#' variables by the proportion of variance their component explains.
#' This is the most widely cited data-driven approach in soil quality
#' research (Andrews et al., 2004; Bastida et al., 2008).
#'
#' \deqn{SQI_{PCA} = \sum_{k=1}^{m} \frac{V_k}{\sum V} \bar{S}_{g, j_k}}
#'
#' where \eqn{V_k} is the variance explained by component \eqn{k},
#' \eqn{j_k} is the MDS variable selected from component \eqn{k}, and
#' \eqn{\bar{S}_{g, j_k}} is the group mean score of that variable.
#'
#' @param scored A scored data frame from \code{\link{score_all}}.
#' @param config A \code{sqi_config} object.
#' @param group_cols Character vector of grouping column names.
#' @param mds Object returned by \code{\link{select_mds}}.  If \code{NULL},
#'   \code{select_mds} is run internally with default parameters.
#'
#' @return A data frame with group columns and \code{SQI_pca} (0--1).
#'
#' @references
#' Andrews, S.S., Karlen, D.L., & Cambardella, C.A. (2004). The soil
#' management assessment framework. \emph{Soil Science Society of America
#' Journal}, 68(6), 1945--1962. \doi{10.2136/sssaj2004.1945}
#'
#' Bastida, F., Zsolnay, A., Hernandez, T., & Garcia, C. (2008). Past,
#' present and future of soil quality indices: A biological perspective.
#' \emph{Geoderma}, 147(3--4), 159--171.
#' \doi{10.1016/j.geoderma.2008.08.007}
#'
#' @examples
#' data(soil_data)
#' cfg <- make_config(
#'   variable = c("pH","EC","BD","OC","MBC","Clay"),
#'   type     = c("opt","less","less","more","more","opt"),
#'   opt_low  = c(6.0, NA, NA, NA, NA, 20),
#'   opt_high = c(7.0, NA, NA, NA, NA, 35)
#' )
#' scored <- score_all(soil_data, cfg, group_cols = c("LandUse","Depth"))
#' result <- sqi_pca(scored, cfg, group_cols = c("LandUse","Depth"))
#' print(result)
#'
#' @export
sqi_pca <- function(scored, config, group_cols = "LandUse", mds = NULL) {

  if (is.null(mds))
    mds <- select_mds(scored, group_cols = group_cols, verbose = FALSE)

  sel_vars  <- mds$mds_vars
  n_comp    <- ncol(mds$loadings)
  var_expl  <- mds$var_explained[seq_len(n_comp)]
  weights   <- var_expl / sum(var_expl)

  # Pair each component weight with its selected variable
  comp_vars <- character(n_comp)
  for (pc in seq_len(n_comp)) {
    ld     <- abs(mds$loadings[, pc])
    cands  <- names(ld[ld >= 0.6])
    cands  <- intersect(cands, sel_vars)
    if (length(cands) == 0) cands <- names(which.max(ld))
    comp_vars[pc] <- cands[1]
  }

  gm       <- .group_means(scored, group_cols, unique(comp_vars))
  raw      <- numeric(nrow(gm))
  for (pc in seq_len(n_comp)) {
    raw <- raw + weights[pc] * gm[[comp_vars[pc]]]
  }

  result          <- gm[, group_cols, drop = FALSE]
  result$SQI_pca  <- round(.min_max(raw / n_comp), 4)
  result[order(-result$SQI_pca), ]
}


#' Soil Quality Index: Fuzzy Logic Method
#'
#' @description
#' Computes SQI using a fuzzy membership aggregation approach.  Each scored
#' variable (already 0--1) is treated as a fuzzy membership value, and
#' groups are aggregated using either the arithmetic mean (equivalent to the
#' linear method) or the fuzzy weighted average operator.
#'
#' This approach is appropriate when variable importance is uncertain or
#' when expert-elicited weights are available
#' (Zhu et al., 2006; Torbert & Wood, 1992).
#'
#' @param scored A scored data frame from \code{\link{score_all}}.
#' @param config A \code{sqi_config} object.
#' @param group_cols Character vector of grouping column names.
#' @param mds_vars Character vector of MDS variable names.
#' @param fuzzy_weights Named numeric vector of fuzzy importance weights
#'   (sum need not equal 1; they are normalised internally).
#'   Defaults to equal weights.
#' @param operator Character. Aggregation operator: \code{"mean"} (default)
#'   or \code{"geometric"} (product-based, penalises low scores on any
#'   variable).
#'
#' @return A data frame with group columns and \code{SQI_fuzzy} (0--1).
#'
#' @references
#' Zhu, A.X., Liu, F., Li, B., Pei, T., Qin, C., Liu, G., Wang, Y.,
#' Chen, Y., Ma, X., Qi, F., & Li, R. (2010). Differentiation of soil
#' conditions over flat areas using land surface feedback dynamic patterns
#' extracted from MODIS. \emph{Soil Science Society of America Journal},
#' 74(1), 861--869.
#'
#' Torbert, H.A., & Wood, C.W. (1992). Effects of soil compaction and
#' water-filled pore space on soil microbial activity and N losses.
#' \emph{Communications in Soil Science and Plant Analysis}, 23,
#' 1321--1331. \doi{10.1080/00103629209368668}
#'
#' @examples
#' data(soil_data)
#' cfg <- make_config(
#'   variable = c("pH","EC","BD","OC","MBC","Clay"),
#'   type     = c("opt","less","less","more","more","opt"),
#'   opt_low  = c(6.0, NA, NA, NA, NA, 20),
#'   opt_high = c(7.0, NA, NA, NA, NA, 35)
#' )
#' scored <- score_all(soil_data, cfg, group_cols = c("LandUse","Depth"))
#' result <- sqi_fuzzy(scored, cfg, group_cols = c("LandUse","Depth"))
#' print(result)
#'
#' @export
sqi_fuzzy <- function(scored, config, group_cols = "LandUse",
                      mds_vars = NULL, fuzzy_weights = NULL,
                      operator = c("mean", "geometric")) {

  operator <- match.arg(operator)
  vars <- if (!is.null(mds_vars)) mds_vars else config$variable
  vars <- intersect(vars, names(scored))

  gm <- .group_means(scored, group_cols, vars)

  w <- if (!is.null(fuzzy_weights)) {
    fw <- fuzzy_weights[vars]
    fw[is.na(fw)] <- 1
    fw / sum(fw)
  } else {
    rep(1 / length(vars), length(vars))
  }

  score_mat <- as.matrix(gm[, vars])
  raw <- if (operator == "mean") {
    as.numeric(score_mat %*% w)
  } else {
    # Geometric mean (penalises any near-zero score)
    apply(score_mat, 1, function(row) {
      prod(pmax(row, 1e-6)^w)
    })
  }

  result           <- gm[, group_cols, drop = FALSE]
  result$SQI_fuzzy <- round(.min_max(raw), 4)
  result[order(-result$SQI_fuzzy), ]
}


#' Soil Quality Index: Entropy Weighting Method
#'
#' @description
#' Computes SQI using Shannon entropy to derive objective weights for each
#' variable.  Variables with higher information entropy (greater
#' discriminating power among groups) receive higher weights.  This removes
#' subjectivity from weight assignment.
#'
#' The entropy weight for variable \eqn{j} is:
#' \deqn{e_j = -\frac{1}{\ln n} \sum_{i=1}^{n} p_{ij} \ln(p_{ij})}
#' \deqn{w_j = \frac{1 - e_j}{\sum_k (1 - e_k)}}
#'
#' where \eqn{p_{ij} = \bar{S}_{ij} / \sum_i \bar{S}_{ij}}.
#'
#' @param scored A scored data frame from \code{\link{score_all}}.
#' @param config A \code{sqi_config} object.
#' @param group_cols Character vector of grouping column names.
#' @param mds_vars Character vector of MDS variable names.
#'
#' @return A data frame with group columns, \code{SQI_entropy},
#'   and attribute \code{entropy_weights} (named numeric vector).
#'
#' @references
#' Shannon, C.E. (1948). A mathematical theory of communication.
#' \emph{Bell System Technical Journal}, 27(3), 379--423.
#' \doi{10.1002/j.1538-7305.1948.tb01338.x}
#'
#' Li, P., Qian, H., & Wu, J. (2010). Groundwater quality assessment based
#' on improved water quality index in Pengyang County, Ningxia, Northwest
#' China. \emph{E-Journal of Chemistry}, 7, 209--216.
#' \doi{10.1155/2010/451304}
#'
#' @examples
#' data(soil_data)
#' cfg <- make_config(
#'   variable = c("pH","EC","BD","OC","MBC","Clay"),
#'   type     = c("opt","less","less","more","more","opt"),
#'   opt_low  = c(6.0, NA, NA, NA, NA, 20),
#'   opt_high = c(7.0, NA, NA, NA, NA, 35)
#' )
#' scored <- score_all(soil_data, cfg, group_cols = c("LandUse","Depth"))
#' result <- sqi_entropy(scored, cfg, group_cols = c("LandUse","Depth"))
#' attr(result, "entropy_weights")
#'
#' @export
sqi_entropy <- function(scored, config, group_cols = "LandUse",
                        mds_vars = NULL) {

  vars <- if (!is.null(mds_vars)) mds_vars else config$variable
  vars <- intersect(vars, names(scored))

  gm       <- .group_means(scored, group_cols, vars)
  score_m  <- as.matrix(gm[, vars])
  n        <- nrow(score_m)

  # Normalise columns to get proportions p_ij
  col_sums <- colSums(score_m)
  col_sums[col_sums == 0] <- 1e-9
  P <- sweep(score_m, 2, col_sums, "/")
  P[P == 0] <- 1e-9   # avoid log(0)

  # Entropy per variable
  e_j <- -colSums(P * log(P)) / log(n)
  e_j <- pmin(e_j, 1)
  w_j <- (1 - e_j) / sum(1 - e_j)

  raw <- as.numeric(score_m %*% w_j)

  result              <- gm[, group_cols, drop = FALSE]
  result$SQI_entropy  <- round(.min_max(raw), 4)
  attr(result, "entropy_weights") <- round(w_j, 4)
  result[order(-result$SQI_entropy), ]
}


#' Soil Quality Index: TOPSIS Method
#'
#' @description
#' Computes SQI using the Technique for Order of Preference by Similarity
#' to Ideal Solution (TOPSIS), a multi-criteria decision analysis method.
#' Each group is ranked by its Euclidean distance to the positive ideal
#' solution (all scores = 1) and negative ideal solution (all scores = 0).
#'
#' \deqn{C_i^* = \frac{d_i^-}{d_i^+ + d_i^-}}
#'
#' where \eqn{d_i^+} and \eqn{d_i^-} are distances to the positive and
#' negative ideal solutions.  \eqn{C_i^* \in [0, 1]} with higher values
#' indicating better soil quality.
#'
#' @param scored A scored data frame from \code{\link{score_all}}.
#' @param config A \code{sqi_config} object.
#' @param group_cols Character vector of grouping column names.
#' @param mds_vars Character vector of MDS variable names.
#' @param weights Named numeric vector of criteria weights.  Defaults to
#'   equal weights.
#'
#' @return A data frame with group columns and \code{SQI_topsis} (0--1).
#'
#' @references
#' Hwang, C.L., & Yoon, K. (1981). \emph{Multiple Attribute Decision
#' Making: Methods and Applications}. Springer, Berlin.
#' \doi{10.1007/978-3-642-48318-9}
#'
#' Yoon, K. (1987). A reconciliation among discrete compromise solutions.
#' \emph{Journal of the Operational Research Society}, 38, 277--286.
#' \doi{10.1057/jors.1987.44}
#'
#' @examples
#' data(soil_data)
#' cfg <- make_config(
#'   variable = c("pH","EC","BD","OC","MBC","Clay"),
#'   type     = c("opt","less","less","more","more","opt"),
#'   opt_low  = c(6.0, NA, NA, NA, NA, 20),
#'   opt_high = c(7.0, NA, NA, NA, NA, 35)
#' )
#' scored <- score_all(soil_data, cfg, group_cols = c("LandUse","Depth"))
#' result <- sqi_topsis(scored, cfg, group_cols = c("LandUse","Depth"))
#' print(result)
#'
#' @export
sqi_topsis <- function(scored, config, group_cols = "LandUse",
                       mds_vars = NULL, weights = NULL) {

  vars <- if (!is.null(mds_vars)) mds_vars else config$variable
  vars <- intersect(vars, names(scored))

  gm      <- .group_means(scored, group_cols, vars)
  score_m <- as.matrix(gm[, vars])

  # Normalise by vector norm
  norms      <- sqrt(colSums(score_m^2))
  norms[norms == 0] <- 1e-9
  norm_m     <- sweep(score_m, 2, norms, "/")

  # Weights
  w <- if (!is.null(weights)) {
    weights[vars]
  } else {
    rep(1 / length(vars), length(vars))
  }
  w[is.na(w)] <- 1 / length(vars)
  w <- w / sum(w)

  weighted_m <- sweep(norm_m, 2, w, "*")

  # Ideal solutions
  pis <- apply(weighted_m, 2, max)   # positive ideal
  nis <- apply(weighted_m, 2, min)   # negative ideal

  d_pos <- sqrt(rowSums(sweep(weighted_m, 2, pis, "-")^2))
  d_neg <- sqrt(rowSums(sweep(weighted_m, 2, nis, "-")^2))

  C <- d_neg / (d_pos + d_neg)

  result            <- gm[, group_cols, drop = FALSE]
  result$SQI_topsis <- round(C, 4)
  result[order(-result$SQI_topsis), ]
}


#' Compare All SQI Methods
#'
#' @description
#' Runs all six SQI methods (\code{sqi_linear}, \code{sqi_regression},
#' \code{sqi_pca}, \code{sqi_fuzzy}, \code{sqi_entropy}, \code{sqi_topsis})
#' on the same scored dataset and returns a combined results table for
#' method comparison.
#'
#' @param scored A scored data frame from \code{\link{score_all}}.
#' @param config A \code{sqi_config} object.
#' @param group_cols Character vector of grouping column names.
#' @param dep_var Character. Dependent variable for \code{sqi_regression}.
#'   If \code{NULL}, the regression method is skipped.
#' @param mds Object from \code{\link{select_mds}}, or \code{NULL} to
#'   compute automatically.
#'
#' @return A data frame with one row per group and columns for each SQI
#'   method.  Also includes \code{Mean_SQI} and \code{Rank} columns.
#'
#' @examples
#' data(soil_data)
#' cfg <- make_config(
#'   variable = c("pH","EC","BD","OC","MBC","Clay"),
#'   type     = c("opt","less","less","more","more","opt"),
#'   opt_low  = c(6.0, NA, NA, NA, NA, 20),
#'   opt_high = c(7.0, NA, NA, NA, NA, 35)
#' )
#' scored  <- score_all(soil_data, cfg, group_cols = c("LandUse","Depth"))
#' results <- sqi_compare(scored, cfg, group_cols = c("LandUse","Depth"),
#'                         dep_var = "OC")
#' print(results)
#'
#' @export
sqi_compare <- function(scored, config, group_cols = "LandUse",
                        dep_var = NULL, mds = NULL) {

  if (is.null(mds))
    mds <- select_mds(scored, group_cols = group_cols, verbose = FALSE)
  mds_v <- mds$mds_vars

  res_lin  <- sqi_linear(scored, config, group_cols, mds_vars = mds_v)
  res_pca  <- sqi_pca(scored, config, group_cols, mds = mds)
  res_fuz  <- sqi_fuzzy(scored, config, group_cols, mds_vars = mds_v)
  res_ent  <- sqi_entropy(scored, config, group_cols, mds_vars = mds_v)
  res_top  <- sqi_topsis(scored, config, group_cols, mds_vars = mds_v)

  # Use merge for reliable multi-column joining
  .safe_merge <- function(base_df, res_df, sqi_col) {
    merged <- merge(base_df, res_df[, c(group_cols, sqi_col)],
                    by = group_cols, all.x = TRUE, sort = FALSE)
    merged[[sqi_col]]
  }

  base <- res_lin[, group_cols, drop = FALSE]
  base$SQI_linear  <- .safe_merge(base, res_lin, "SQI_linear")
  base$SQI_pca     <- .safe_merge(base, res_pca, "SQI_pca")
  base$SQI_fuzzy   <- .safe_merge(base, res_fuz, "SQI_fuzzy")
  base$SQI_entropy <- .safe_merge(base, res_ent, "SQI_entropy")
  base$SQI_topsis  <- .safe_merge(base, res_top, "SQI_topsis")

  if (!is.null(dep_var)) {
    res_reg <- sqi_regression(scored, config, dep_var, group_cols,
                              mds_vars = mds_v)
    base$SQI_regression <- .safe_merge(base, res_reg, "SQI_regression")
  }

  sqi_cols      <- grep("^SQI_", names(base), value = TRUE)
  base$Mean_SQI <- round(rowMeans(base[, sqi_cols], na.rm = TRUE), 4)
  base$Rank     <- rank(-base$Mean_SQI, ties.method = "min")
  base[order(base$Rank), ]
}