R/calibration_control.R

In proximetricsR: Spectral Preprocessing and Chemometric Calibration of NIR Sensors

#' @title A function that controls the calibration of models
#' @description
#'
#' \loadmathjax
#'
#' This function is used to further control some aspects of the calibration of
#' models (with the \code{\link{calibrate}} function) such as cross-validation
#' and outlier detection.
#' @usage
#' calibration_control(validation_type = c("lgo", "loo", "kfold", "none"),
#'                     number = ifelse(validation_type == "lgo", 100, 10),
#'                     p = 0.75,
#'                     folds = c("random", "sequential"),
#'                     tuning_parameter = c("rmse", "rsq", "none"),
#'                     learning_rates = c(maximum = 1.1, sequential = 1.05),
#'                     remove_outliers = 0,
#'                     cal_residual_limit = 2.5,
#'                     mahalanobis_limit = 5,
#'                     val_residual_limit = 3.5,
#'                     allow_parallel = TRUE,
#'                     fix_pls_factors = TRUE,
#'                     fixed_components = 0,
#'                     replacements = TRUE,
#'                     seed = NULL)
#' @param validation_type a character string indicating the type of
#' cross-validation (cv) to be conducted. Options are: \code{"lgo"} for
#' leave-group-out cv (default), \code{"loo"} for leave-one-out cv and
#' \code{"kfold"} for k-fold cv. See details.
#' @param number an integer indicating the number of sampling iterations or
#' sub-sample groups for the selected \code{validation_type} argument. Default
#' is \code{100} for leave-group out cv and \code{10} for k-fold cross-validation.
#' This parameter is ignored for leave-one-out cv.
#' @param p a numeric value indicating the percentage of calibration observations
#' to be retained at each sampling iteration at each local segment when
#' \code{"lgo"} is selected in the \code{validation_type} argument. Default is
#' 0.75 (i.e. 75 percent of the observations).
#' @param folds a character string indicating the way folds are created (valid
#' only when \code{validation_type = "kfold"}). Options are: \code{"random"}
#' (default) or \code{"sequential"}.
#' @param tuning_parameter a character string indicating which cross-validation
#' statistic to use for the optimization of the included number of components.
#' Options are: \code{"rmse"} (default, minimization of the root mean squared error),
#' \code{"rsq"} (maximization of the coefficient of determination) or
#' \code{"none"} (no tuning). Does not apply when \code{validation_type = "none"}.
#' @param learning_rates a vector of length 2 for additional control over the
#' selection of the optimal number of components. See details for its use. Defaults
#' to \code{c(1.1, 1.05)}.
#' @param remove_outliers an integer indicating the number of times the model should
#' automatically detect and remove outliers. Each time, a new model is fitted
#' with the outliers removed, until either no more outliers are found or the
#' \code{remove_outliers} has been reached. Outliers found and removed in each
#' step, as well as the first and last computed models are recorded. Outliers
#' are detected based on the limits set in the arguments \code{cal_residual_limit},
#' \code{mahalanobis_limit} and \code{val_residual_limit}. Setting
#' \code{remove_outliers} to \code{0} (default) disables automatic outlier removal,
#' whereas selecting it as \code{Inf} removes outliers until no more are found.
#' @param cal_residual_limit a numeric value which indicates the upper limit of
#' the standardized residuals for the fitted response variable. Observations with
#' absolute residuals above this limit are labeled as \code{"calibration outliers"}.
#' The standardized calibration residuals are calculated as the absolute
#' differences between the reference values and their corresponding fitted
#' values divided by the standard deviation of these absolute differences.
#' Default is 2.5 (as in NIRWise PLUS calibration software).
#' @param mahalanobis_limit a numeric value which indicates the upper limit of
#' the squared Mahalanobis distances of each sample in the score space to zero.
#' Observations with squared Mahalanobis distance above this limit are labeled as
#' \code{"Mahalanobis outliers"}. The squared Mahalanobis distances are calculated as
#' the squared Euclidean distance of the standardized scores to the origin.
#' Default is 5 (as in NIRWise PLUS calibration software).
#' @param val_residual_limit a numeric value which indicates the upper limit of the
#' standardized residuals for cross-validation predictions of the response
#' variable. This applies only to \code{"kfold"} or \code{"loo"} cross-validation.
#' Observations with absolute residuals above this limit are labeled as
#' \code{"validation outliers"}. The standardized validation residuals are
#' calculated as the absolute differences between the reference values
#' and their corresponding cross-validated predictions divided by the standard
#' deviation of these absolute differences. Default is 3.5 (as in NIRWise PLUS
#' calibration software).
#' @param allow_parallel a logical indicating if parallel execution is allowed.
#' If \code{TRUE}, parallelization is applied to the cross-validation procedure.
#' The parallelization of this for loop is implemented using the
#' \code{\link[foreach]{foreach}} function of the \code{foreach} package.
#' Default is \code{TRUE}.
#' @param fix_pls_factors a logical. This parameter only has an influence on the
#' produced application files, where it indicates whether the final number of
#' factors of the model should be fixed. Note that this has no influence on the
#' model in R itself, as the optimal number of components inside the model
#' remains the same (but it does influence the exported files). Default is
#' \code{TRUE}.
#' @param fixed_components a numerical value indicating a fixed number of
#' components to be used in the model (i.e. no optimization of the components).
#' The default value is \code{0}, which indicates that the number of components
#' is not fixed and it uses the one selected by the function.
#' @param replacements a logical. Only used in case \code{validation_type} is
#' selected as \code{"lgo"}. Specifies if the sampling for the calibration sets
#' must be done with replacements. See details for a more thorough explanation.
#' Defaults to \code{TRUE}.
#' @param seed an integer that can be used in any of the validation methods to
#' obtain reproducible results, using the \code{\link{set.seed}} function.
#' In case it is selected as \code{NULL}, no seed will be set. Note that the
#' seed will not be reset, and future random computations can be affected.
#' Furthermore, this parameter is meant as a way to provide reproducibility,
#' and should not be used to simply select the seed with the
#' best results. Default is \code{NULL}.
#' @details
#' This package extends the cross-validation methods implemented in the
#' NIRWise PLUS software, which is based only on k-fold cross validation.
#'
#' The validation methods available for assessing the predictive performance
#' of the models are:
#'  \itemize{
#'  \item \strong{Leave-group-out cross-validation (\code{"lgo"}):} The
#'  data is partitioned into different subsets of similar size. Each partition
#'  is based on a stratified random sampling using the distribution of the
#'  response variable. When \code{p} \mjeqn{\ge}{\geqslant} 0.5 (i.e.
#'  the number of calibration observations to retain is larger than 50% of the
#'  total samples), the sampling is conducted for selecting the validation
#'  samples, and when \code{p} is below 0.5 the sampling is conducted for
#'  selecting the calibration samples (samples used for model training). The
#'  model fitted with the selected calibration samples is used to predict the
#'  target response variable values of the validation samples. The accuracy
#'  and precision, indicated by the root mean square error (RMSE) and the
#'  coefficient of determination
#'  (\mjeqn{R^2}{R^2}) respectively, are computed. This process is repeated
#'  \mjeqn{m}{m} times (where \mjeqn{m}{m} is controlled by the \code{number}
#'  argument), and the final RMSE and \mjeqn{R^2}{R^2} are computed as the
#'  average over all respective results of the \mjeqn{m}{m} iterations. In case
#'  the parameter \code{replacements} is set to \code{TRUE}, the selection of the
#'  calibration sets is done by using sampling with replacement.
#'
#'  \item \strong{Leave-one-out cross-validation (\code{"loo"}):} The number of
#'  iterations is equal to the number of observations in the calibration set.
#'  In each iteration, one single observation is held out, while the remaining
#'  samples are used to fit a model, which is used to predict the response
#'  variable of the held out observation. The predictions are then compared
#'  to the reference ones and both the RMSE and the (\mjeqn{R^2}{R^2}) are
#'  computed.
#'
#'  \item \strong{k-fold cross-validation (\code{"kfold"}):} The data is split (either
#'  randomly or sequentially) into \mjeqn{k}{k} disjoint blocks of similar size,
#'  where \mjeqn{k}{k} is controlled by \code{number}. In the sequential splits,
#'  every block \mjeqn{B_i}{B_i} is selected as follows:
#'
#' \mjdeqn{B_i = \lbrace i + k(j - 1) | j \in N, \ i + k(j - 1) \leq n \rbrace}{B_i = \lbrace i + k(j - 1)| j \in \mathbb{N}, \ i + k(j - 1) \leq n \rbrace}
#'
#'  where \mjeqn{n}{n} is the total number of observations. In other words, the
#'  observations are put sequentially into the blocks until all observations
#'  have a block assigned.\cr
#'  A total of \mjeqn{k}{k} iterations is conducted. In each iteration, one block
#'  is considered as the validation set, while the remaining samples are used to
#'  fit a model, which is then used to predict the response variable of the
#'  held-out block.\cr
#'  The number observations in each block is given by the total number of
#'  observations divided by the number of blocks. Note that the maximum number
#'  of folds is limited to half of the number of observations. Note also that
#'  this implementation of k-fold cross-validation is an improved version of the
#'  one in the NIRWise PLUS software, where only the sequential sample selection
#'  is supported.
#'
#'  \item \strong{No validation (\code{"none"}):} No validation is carried out.
#'  }
#'
#'  For each validation type (except \code{"none"}), the optimal number of
#'  factors is not necessarily chosen to be the minimum of RMSE or the maximum of
#'  \mjeqn{R^2}{R^2} (depending on the \code{tuning_parameter}). Instead, since
#'  both are often monotonically decreasing respectively monotonically increasing
#'  as the number of components increases, an additional parameter \code{learning_rates}
#'  \mjeqn{\gamma}{\gamma} for fine-tuning of the determination of the number of
#'  factors is included:
#'
#'  For RMSE, consider the index where the minimum of all computed RMSE is attained:
#'
#'  \mjdeqn{n_{min} = arg\min_{n} \  RMSE_n}{n_{min} = arg\min_{n} \  RMSE_n},
#'
#'  Then, among all \mjteqn{1 < n < n_{min}}{1 \lt n \lt n_{min}}{1 \lt n \lt n_{min}} fulfilling
#'
#'  \mjtdeqn{RMSE_{n} < RMSE_{n_{min}} \cdot \gamma_{max}}{RMSE_{n} \lt RMSE_{n_{min}} \cdot \gamma_{max}}{RMSE_{n} \lt RMSE_{n_{min}} \cdot \gamma_{max}}
#'
#'  \mjtdeqn{RMSE_{n} < RMSE_{n+1} \cdot \gamma_{seq}}{RMSE_{n} \lt RMSE_{n+1} \cdot \gamma_{seq}}{RMSE_{n} \lt RMSE_{n+1} \cdot \gamma_{seq}}
#'
#'  we take the smallest \mjeqn{n}{n} as the optimal number of components.\cr
#'  For \mjeqn{R^2}{R^2}, a similar approach is taken, but with maxima instead of
#'  minima: \mjeqn{n_{max} = arg\max_{n} R^2_n}{n_{max} = arg\max_{n} R^2_n}
#'  Then, take the smallest \mjteqn{1 < n < n_{max}}{1 \lt n \lt n_{max}}{1 \lt n \lt n_{max}} still satisfying
#'
#'  \mjtdeqn{R^2_{n} > R^2_{n_{max}} \cdot \gamma_{max}^{-1}}{R^2_{n} \gt R^2_{n_{max}} \cdot \gamma_{max}^{-1}}{R^2_{n} \gt R^2_{n_{max}} \cdot \gamma_{max}^{-1}}
#'
#'  \mjtdeqn{R^2_{n} > R^2_{n+1} \cdot \gamma_{seq}^{-1}}{R^2_{n} \gt R^2_{n+1} \cdot \gamma_{seq}^{-1}}{R^2_{n} \gt R^2_{n+1} \cdot \gamma_{seq}^{-1}}
#'
#'  Note that in this case, we take the inverse of the learning rates. Furthermore,
#'  setting \code{learning_rates = c(1, 1)} retains the
#'  global minimum for RMSE, respectively maximum for \mjeqn{R^2}{R^2}.
#'
#' @return a list of class \code{calibration_control} mirroring
#' the specified parameters
#' @author Leonardo Ramirez-Lopez
#' @seealso \code{\link{calibrate}}, \code{\link{calibrate_models}}
#' @examples
#'
#' # 5-fold cross-validation with sequential sampling
#' calibration_control(
#'   validation_type = "kfold",
#'   number = 5,
#'   folds = "sequential"
#' )
#'
#' # leave-one-out cross_validation
#' calibration_control(validation_type = "loo")
#'
#' # 100 leave-group-out validations with 60% samples retained, with replacements
#' calibration_control(
#'   validation_type = "lgo",
#'   number = 100,
#'   p = 0.6,
#'   replacements = TRUE
#' )
#'
#' # 2-fold leave-group-out cross-validation with 75% samples retained, no replacements
#' calibration_control(
#'   validation_type = "lgo",
#'   number = 2,
#'   p = 0.75,
#'   replacements = FALSE
#' )
#'
#' # Same as before, but removing any outlier that is found
#' calibration_control(
#'   validation_type = "lgo",
#'   number = 2,
#'   p = 0.75,
#'   replacements = FALSE,
#'   remove_outliers = Inf
#' )
#' \donttest{
#' # no validation, gives warning
#' calibration_control(validation_type = "none")
#' }
#' @export
calibration_control <- function(
  validation_type = c("lgo", "loo", "kfold", "none"),
  number = ifelse(validation_type == "lgo", 100, 10),
  p = 0.75,
  folds = c("random", "sequential"),
  tuning_parameter = c("rmse", "rsq", "none"),
  learning_rates = c(maximum = 1.1, sequential = 1.05),
  remove_outliers = 0,
  cal_residual_limit = 2.5,
  mahalanobis_limit = 5,
  val_residual_limit = 3.5,
  allow_parallel = TRUE,
  fix_pls_factors = TRUE,
  fixed_components = 0,
  replacements = TRUE,
  seed = NULL
) {
  validation_type <- match.arg(validation_type)
  folds <- match.arg(folds)
  tuning_parameter <- match.arg(tuning_parameter)
  # Sanity checks
  if (!is.logical(allow_parallel)) {
    stop("allow_parallel must be a logical value")
  }
  if (!is.numeric(number)) {
    stop("'number' must be numeric")
  }
  if (validation_type == "none" & tuning_parameter %in% c("rmse", "rsq")) {
    warning("Cross-validation is required for model tuning, parameter tuning will not be conducted")
  }
  if (!is.numeric(remove_outliers)) {
    stop("'remove_outliers' must be an integer or 'Inf'.")
  }
  if (!is.numeric(p) || length(p) != 1 | p >= 1 | p <= 0) {
    stop("p must be a single numeric value larger than 0 and below than 1")
  }
  if (length(learning_rates) != 2) {
    stop("'learning_rates' must be of length 2.")
  }
  if (!is.numeric(learning_rates)) {
    stop("'learning_rates' must a vector of numericals.")
  }
  if (!is.logical(fix_pls_factors)) {
    stop("'fix_pls_factors' must be a logical.")
  }
  if (!is.numeric(fixed_components)) {
    stop("'fixed_components' must be numeric")
  }
  if (!is.logical(replacements)) {
    stop("'replacements' must be a logical.")
  }
  if (!is.null(seed) & !is.numeric(seed)) {
    stop("'seed' must be either NULL or an integer.")
  }
  if (!is.numeric(cal_residual_limit)) {
    stop("'cal_residual_limit' must be numerical.")
  }
  if (!is.numeric(mahalanobis_limit)) {
    stop("'mahalanobis_limit' must be numerical.")
  }
  if (!is.numeric(val_residual_limit)) {
    stop("'val_residual_limit' must be numerical.")
  }
  cntrl <- list(
    validation_type = validation_type,
    number = number,
    p = p,
    folds = folds,
    tuning_parameter = tuning_parameter,
    learning_rates = learning_rates,
    remove_outliers = remove_outliers,
    cal_residual_limit = cal_residual_limit,
    mahalanobis_limit = mahalanobis_limit,
    val_residual_limit = val_residual_limit,
    allow_parallel = allow_parallel,
    fix_pls_factors = fix_pls_factors,
    fixed_components = fixed_components,
    replacements = replacements,
    seed = seed
  )
  class(cntrl) <- c("calibration_control", "list")
  cntrl
}