R/frequentist_unconstrained.R

In OceanTracer: Marine Source Unmixing Toolkit

#' Source apportionment via unconstrained least-squares with classical variability propagation
#'
#' Estimates the fractional contribution of each sediment source to a mixture
#' by solving an overdetermined mass-balance system with least-squares
#' minimisation. Source variability is propagated by Monte Carlo resampling:
#' at each iteration, synthetic source signatures are drawn from a normal
#' distribution defined by the measured mean and standard deviation, and the
#' system is solved independently for each draw.
#'
#' This function is the computational back-end called by \code{\link{unmixing}}.
#' For most users, calling \code{unmixing} directly is recommended.
#'
#' @param data Data frame in the standard OceanTracer format (sources + mixture).
#' @param variability Character. \code{"SEM"} for Standard Error of the Mean,
#'   \code{"SD"} for Standard Deviation.
#' @param iter Integer or character. Number of Monte Carlo iterations (e.g., \code{1000}, or character keywords like \code{"short"}, \code{"medium"}, \code{"long"}). Default is 1000.
#' @param means Logical. \code{TRUE} when the input is in averaged format
#'   (\code{Mean_}/\code{SD_} columns); \code{FALSE} for raw observations.
#' @param seed Integer. Random seed for reproducibility. Default is 123456L.
#' @param verbose Logical. Whether to show progress bars and summary messages in the console. Default is TRUE.
#'
#' @return A data frame with one row per iteration containing the mixture ID,
#'   goodness-of-fit (GOF), and one column per source with the estimated
#'   fractional contribution.
#'
#' @importFrom stats lm.fit coef rt sd
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @importFrom crayon cyan green bold
#'
#' @export
frequentist_unconstrained <- function(data, variability = "SEM", iter = 1000, means = FALSE, seed = 123456L, verbose = TRUE)
{
  set.seed(seed)
  
  sr_names <- unique(data[, 2])
  
  if (is.numeric(iter)) { 
    iter <- iter
  } else {
    if (iter == "very short") iter <- 500
    if (iter == "short")      iter <- 1000
    if (iter == "medium")     iter <- 5000
    if (iter == "long")       iter <- 20000
    if (iter == "bananas")    iter <- 100000
  }
  
  if (means == TRUE) {
    source <- data[c(1:(nrow(data) - 1)), c(2:ncol(data))]
    x <- round((ncol(data) - 3) / 2 + 2)
    mixture <- inputMixture(data[nrow(data), 1:x])
    colnames(source)[1] <- "id"
    colnames(mixture) <- colnames(source[, c(1:(ncol(source) / 2))])
    
  } else {
    source <- inputSource(data)
    mixture <- inputMixture(data)
  }
  
  if (nrow(source) - 1 >= ncol(mixture)) { 
    warning("As a minimum, n - 1 properties are required to discriminate rigorously between n sources. Additional properties are frequently required to increase the reliability of the results")
  }
  
  snames <- source[, 1]
  
  source <- data.matrix(source[-1])
  mixture <- data.matrix(mixture[-1])
  
  # normalize
  cols <- (ncol(source) - 1) / 2
  for (col in c(1:cols)) {
    mx <- max(source[, col] + source[, cols + col])
    mn <- min(source[, col] - source[, cols + col])
    source[, col] <- (source[, col] - mn) / (mx - mn)
    source[, cols + col] <- source[, cols + col] / (mx - mn)
    mixture[, col] <- (mixture[, col] - mn) / (mx - mn)
  }
  
  nsource <- nrow(source)
  ntracer <- (ncol(source) - 1) / 2
  
  r <- matrix(nrow = iter, ncol = nsource + 2)
  csource <- matrix(nrow = nsource, ncol = ntracer)
  
  if (verbose) {
    message("Summary of the model inputs:")
    message(paste(ntracer, "variables from", nsource, "sources"))
  }
  
  # Introduce the progress bar
  if (verbose) {
    pb <- utils::txtProgressBar(min = 0, max = iter, style = 3, width = 50, char = "=")
  }
  
  for (i in c(1:iter)) {
    y <- c()
    x <- matrix(nrow = ntracer, ncol = nsource - 1)
    
    if (i == 2 || i == 3) { # central solution
      for (j in c(1:ntracer)) {
        y <- c(y, mixture[1, j][[1]] - source[nsource, j])
        for (k in c(1:(nsource - 1))) {
          x[j, k] <- source[k, j] - source[nsource, j]
        }
      }
      for (j in c(1:ntracer)) {
        for (k in c(1:nsource)) {
          csource[k, j] <- source[k, j]
        }
      }
    } else {
      ls <- c()
      for (j in c(1:ntracer)) {
        if (variability == "SEM") {
          x1 <- stats::rt(1, source[nsource, ntracer * 2 + 1]) / sqrt(source[nsource, ntracer * 2 + 1])
        } else if (variability == "SD") {
          x1 <- stats::rt(1, source[nsource, ntracer * 2 + 1])
        } else {
          stop("Variability method not implemented. Please use 'SEM' or 'SD'.")
        }
        
        csource[nsource, j] <- source[nsource, j] + source[nsource, ntracer + j] * x1
        ls <- c(ls, csource[nsource, j])
      }
      
      for (j in c(1:ntracer)) {
        y <- c(y, mixture[1, j][[1]] - ls[j])
        for (k in c(1:(nsource - 1))) {
          if (variability == "SEM") {
            x1 <- stats::rt(1, source[k, ntracer * 2 + 1]) / sqrt(source[k, ntracer * 2 + 1])
          } else if (variability == "SD") {
            x1 <- stats::rt(1, source[k, ntracer * 2 + 1])
          } else {
            stop("Variability method not implemented. Please use 'SEM' or 'SD'.")
          }
          csource[k, j] <- source[k, j] + source[k, ntracer + j] * x1
          x[j, k] <- csource[k, j] - ls[j]
        }
      }
    }
    
    # least squares method
    model <- stats::lm.fit(x, y)
    w <- as.vector(stats::coef(model))
    w <- c(w, 1 - sum(w))
    
    gof <- c()
    for (j in c(1:ntracer)) {
      x1 <- 0.0
      for (k in c(1:nsource)) {
        x1 <- x1 + w[k] * csource[k, j]
      }
      gof <- c(gof, (mixture[1, j][[1]] - x1)^2)
    }
    
    gof <- 1.0 - mean(gof)
    gof <- max(gof, 0)
    
    w <- c(1, gof, w)
    r[i, ] <- w
    
    if (verbose) {
      utils::setTxtProgressBar(pb, min(i))
    }
  }
  
  if (verbose) {
    close(pb)
  }
  
  r <- as.data.frame(r)
  colnames(r) <- c('id', 'gof', paste(sr_names[c(1:nsource)], sep = ""))
  
  if (verbose) {
    message("\nSummary of the model outputs:")
    message(paste("Unmixing completed over", iter, "iterations incorporating source variability."))
  }
  
  invisible(return(r))
}