R/cosi.R

In sergioventurini/SAuR: Sensitivity Analysis Using R

#' Calculation of Sensitivity Indexes From Given Data.
#'
#' Returns the sensitivity indexes for given input and output arguments
#    based upon a discrete cosine transformation.
#'
#' @param data An object of class \code{SAuR_data} containing the model's data
#' @param M A length-one numeric vector indicating the number of coefficients
#' @param gfx A length-one character vector providing the conditional regression
#'   line plot main title
#' @param outcol A length-one numeric vector indicating the output column
#'    to focus on
#' @return A \code{list} object containing the sensitivity measures described above.
#' @author Sergio Venturini \email{sergio.venturini@unito.it}
#' @seealso \code{\link{betaKS3_R}} for computing some sensitivity measures.
#' @references
#'   Venturini, S., Borgonovo, E. (2020), "Sensitivity Analysis Using 
#'   \code{R}: the \pkg{SAuR} Package", Technical report.
#' @examples
#' \dontrun{
#' data(simdata_sub, package = "SAuR")
#'
#' }
#'
#' @export
cosi <- function(data, M = NULL, gfx = NULL, outcol = NULL) {
  check_rstudio <- try(rstudioapi::isAvailable(), silent = TRUE)

  x <- data@X
  y <- data@Y
  if (!is.null(outcol))
    y <- y[, outcol, drop = FALSE]
  n <- data@n
  k <- data@p
  if (!is.matrix(y)) y <- as.matrix(y)
  nn <- nrow(y)
  kk <- ncol(y)

  if (nn != n) {
    stop("input/output size do not match.")
  }
  index <- apply(x, 2, order)
  xr <- apply(x, 2, sort)

  yr <- matrix(0, nrow = n, ncol = (k*kk))
  for (i in 1:kk) {
    z <- y[, i]
    for (j in 1:k) {
      yr[, (i - 1)*k + j] <- z[index[, j]]
    }
  }

  # frequency selection
  if (is.null(M)) {
    M <- max(c(ceiling(sqrt(n)), 3))
    print(paste0("cosi() function: using ", as.integer(M), " coefficients"))
  }

  if (M < 0) {
    M <- -M
    unscaled <- TRUE
  } else {
    unscaled <- FALSE
  }

  # compute transformation
  allcoeff <- matrix(0, nrow = nrow(yr), ncol = ncol(yr))
  for (i in 1:ncol(yr)) {
    allcoeff[, i] <- dct_fft(yr[, i], norm = TRUE)
  }

  # transformation is orthogonal, so by Parseval's theorem
  V <- colSums(allcoeff[-1, ]^2)
  if (M == 0) {
    # estimate approximation error or mean
    Si <- 1 - median(n*cumsum(t(allcoeff[seq(nrow(yr), 2), ]^2))/as.numeric(outer(0:(n - 2), V)))
  } else {
    Vi <- colSums(allcoeff[2:M, ]^2)
    if (!unscaled) {
      Si <- Vi/V
    } else {
      Si <- Vi/(n - 1)
      V <- V/(n - 1)
    }
  }

  # return predictions
  yhats <- coeffs <- list(kk)
  for (j in seq(kk)) {
    yhats_j <- matrix(0, nrow = n, ncol = k)
    coeffs_j <- allcoeff[1:M, (k*(j - 1) + 1):(k*j)]
    for (i in 1:k) {
      yhat <- numeric(n)
      yhat[1:(M + 1)] <- allcoeff[1:(M + 1), k*(j - 1) + i]
      yhats_j[, i] <- idct_fft(yhat, norm = TRUE)
    }
    yhats[[j]] <- yhats_j
    coeffs[[j]] <- coeffs_j
  }
  if (is.null(colnames(data@Y))) {
    res_names <- paste0("Y", seq(kk))
  } else {
    res_names <- colnames(data@Y)
  }
  if (!is.null(outcol))
    res_names <- res_names[outcol]
  names(yhats) <- names(coeffs) <- res_names

  if (!is.null(gfx)) {
    for (j in seq(kk)) {
      if ((class(check_rstudio) != "try-error") & check_rstudio) {
        # do nothing --> the graph is automatically provided in a new plot window
      } else {
        dev.new()
      }
      gfxrows <- ifelse(k < 4, k, 3)
      if (k > 1)
        par(mfrow = c(ceiling(k/gfxrows), gfxrows))
      for (i in seq(k)) {
        yhat <- numeric(n)
        yhat[1:(M + 1)] <- allcoeff[1:(M + 1), k*(j - 1) + i]
        plot(x[index[, i], i], y[index[, i], j], col = "#0072B2", pch = 19, cex = .5,
          xlab = paste0("x_", i), ylab = "y")
        lines(x[index[, i], i], idct_fft(yhat, norm = TRUE), col = "#E69F00", lwd = 2)
        title(paste0(as.character(gfx), " - ", res_names[j]))
      }
      par(mfrow = c(1, 1))
    }
  }

  if (kk > 1) {
    dim(Si) <- c(k, kk)
    dim(V) <- c(k, kk)
  }

  return(list(Si = Si, V = V, yhats = yhats, coeffs = coeffs))
}