R/boot.R

In BiGaRx: Unified Classical and Canonical Biplots

#' Biplot con bootstrap: intervalos de confianza
#'
#' Calcula intervalos de confianza bootstrap para los parámetros de un
#' biplot, siguiendo la metodología de Nieto, Galindo, Leiva y
#' Vicente-Galindo (2014), implementada originalmente en `biplotbootGUI`.
#'
#' **Algoritmo.** Se extraen `B` muestras de individuos (filas) con
#' reemplazo. En cada una se repite el preprocesamiento (por ejemplo, se
#' vuelve a estandarizar) y el biplot, y se guardan los parámetros:
#' * `inercia`: porcentaje de inercia de cada eje;
#' * `bondad`: bondad de ajuste de la matriz en el plano;
#' * `coordenada`: coordenadas de cada variable en cada eje del plano;
#' * `calidad`: calidad de representación de cada variable en el plano;
#' * `contribucion`: contribución de cada variable a cada eje;
#' * `angulo_variables`: ángulo entre cada par de variables (en grados);
#' * `angulo_eje`: ángulo entre cada variable y cada eje (en grados).
#'
#' **Alineación de ejes.** La SVD de cada réplica puede devolver un eje con
#' el signo invertido, lo que inflaría artificialmente los intervalos de las
#' coordenadas y de los ángulos con los ejes. Por eso, en cada réplica se
#' invierte el signo del eje si sus cargas apuntan en sentido contrario a las
#' del análisis original. Los ángulos entre variables no dependen del signo.
#' Si dos valores propios son muy parecidos, los ejes pueden intercambiarse
#' entre réplicas; conviene revisar `inercia` antes de interpretar.
#'
#' **Intervalos** (`ci`):
#' * `"percentile"` (boot-q): cuantiles \eqn{\alpha/2} y \eqn{1-\alpha/2} de
#'   las réplicas (Algoritmo 3 del artículo);
#' * `"normal"` (boot-t): estimación \eqn{\pm z_{1-\alpha/2}} por el error
#'   estándar bootstrap, con \eqn{t_{n-1}} si \eqn{n < 25} (Algoritmo 2).
#'
#' @param X Matriz o data.frame numérico.
#' @param B Número de réplicas bootstrap (por defecto 1000).
#' @param method,transform,alpha Como en [bigarx()].
#' @param axes Plano en el que se calculan coordenadas, calidades y ángulos.
#' @param conf Nivel de confianza (por defecto 0,95).
#' @param ci `"percentile"` (por defecto) o `"normal"`.
#' @param seed Semilla opcional para reproducir los resultados.
#' @param groups Factor opcional para colorear el biplot.
#' @param progress Mostrar una barra de progreso.
#'
#' @return Un objeto `"bigarx_boot"` con
#' \describe{
#'   \item{`fit`}{El biplot original ([bigarx()]).}
#'   \item{`table`}{data.frame con una fila por parámetro: tipo, elemento,
#'     eje, estimación, media y error estándar bootstrap, sesgo, límites del
#'     intervalo y método.}
#'   \item{`replicates`}{Matriz de réplicas (una fila por réplica, una
#'     columna por parámetro).}
#'   \item{`B`, `conf`, `ci`, `axes`}{Configuración.}
#' }
#'
#' @references
#' Nieto, A. B., Galindo, M. P., Leiva, V. y Vicente-Galindo, P. (2014). A
#' methodology for biplots based on bootstrapping with R. *Revista
#' Colombiana de Estadística*, 37(2), 367-397. \doi{10.15446/rce.v37n2spe.47944}
#'
#' Efron, B. y Tibshirani, R. J. (1993). *An Introduction to the Bootstrap*.
#' Chapman & Hall.
#'
#' @examples
#' bt <- bigarx_boot(USArrests, B = 200, seed = 1)
#' bt
#' summary(bt, type = "angulo_variables")
#' plot(bt)
#' plot(bt, type = "angulo_variables")
#' @export
bigarx_boot <- function(X, B = 1000, method = c("hj", "gh", "jk", "sqrt"),
                        transform = c("standardize", "center", "double", "grand", "none"),
                        alpha = NULL, axes = c(1, 2), conf = 0.95,
                        ci = c("percentile", "normal"), seed = NULL,
                        groups = NULL, progress = interactive()) {
  cl <- match.call()
  method <- match.arg(method)
  transform <- match.arg(transform)
  ci <- match.arg(ci)
  X <- .as_data_matrix(X)
  if (!is.numeric(B) || B < 50) stop("`B` debe ser al menos 50 (se recomiendan 1000 o m\u00E1s).", call. = FALSE)
  if (!is.numeric(conf) || conf <= 0 || conf >= 1) stop("`conf` debe estar entre 0 y 1.", call. = FALSE)
  if (length(axes) != 2) stop("`axes` debe tener dos ejes.", call. = FALSE)
  nd <- max(axes)
  fit <- bigarx(X, method = method, alpha = alpha, transform = transform,
                ndim = max(3, nd), groups = groups)
  if (nd > fit$ndim) stop("El plano pedido supera el rango de los datos.", call. = FALSE)
  V0 <- fit$.V[, axes, drop = FALSE]

  stat_fun <- function(Xb, align = TRUE) {
    f <- bigarx(Xb, method = method, alpha = alpha, transform = transform,
                ndim = max(3, nd), fix_signs = FALSE)
    if (f$rank < nd) return(NULL)
    H <- f$cols[, axes, drop = FALSE]
    if (align) {
      sg <- sign(colSums(f$.V[, axes, drop = FALSE] * V0))
      sg[sg == 0] <- 1
      H <- sweep(H, 2, sg, "*")
    }
    .boot_stats(f, H, axes)
  }

  est <- stat_fun(X, align = TRUE)
  if (!is.null(seed)) set.seed(seed)
  n <- nrow(X)
  reps <- matrix(NA_real_, B, length(est), dimnames = list(NULL, names(est)))
  pb <- if (progress) utils::txtProgressBar(min = 0, max = B, style = 3) else NULL
  fallos <- 0L
  for (b in seq_len(B)) {
    idx <- sample.int(n, n, replace = TRUE)
    Xb <- X[idx, , drop = FALSE]
    s <- tryCatch(stat_fun(Xb), error = function(e) NULL)
    if (is.null(s)) fallos <- fallos + 1L else reps[b, ] <- s
    if (!is.null(pb)) utils::setTxtProgressBar(pb, b)
  }
  if (!is.null(pb)) close(pb)
  reps <- reps[stats::complete.cases(reps), , drop = FALSE]
  if (fallos > 0) {
    warning(sprintf("%d r\u00E9plicas no se pudieron calcular (por ejemplo, columnas constantes en la muestra) y se descartaron.",
                    fallos), call. = FALSE)
  }

  a <- 1 - conf
  mu <- colMeans(reps)
  se <- apply(reps, 2, stats::sd)
  if (ci == "percentile") {
    lo <- apply(reps, 2, stats::quantile, probs = a / 2, names = FALSE)
    hi <- apply(reps, 2, stats::quantile, probs = 1 - a / 2, names = FALSE)
  } else {
    q <- if (n >= 25) stats::qnorm(1 - a / 2) else stats::qt(1 - a / 2, n - 1)
    lo <- est - q * se
    hi <- est + q * se
  }
  info <- do.call(rbind, strsplit(names(est), "|", fixed = TRUE))
  tab <- data.frame(tipo = info[, 1], elemento = info[, 2], eje = info[, 3],
                    estimacion = unname(est), media = unname(mu), ee = unname(se),
                    sesgo = unname(mu - est), inferior = unname(lo), superior = unname(hi),
                    row.names = NULL, stringsAsFactors = FALSE)
  out <- list(fit = fit, table = tab, replicates = reps, B = nrow(reps), conf = conf,
              ci = ci, axes = axes, call = cl)
  class(out) <- "bigarx_boot"
  out
}

# Parametros de un biplot para el bootstrap. Devuelve un vector con nombres
# "tipo|elemento|eje".
.boot_stats <- function(f, H, axes) {
  vars <- rownames(H)
  ejes <- paste0("Eje", axes)
  out <- c()
  k <- min(3, length(f$inertia$inercia))
  ine <- f$inertia$inercia[seq_len(k)]
  names(ine) <- paste0("inercia|", paste0("Eje", seq_len(k)), "|", paste0("Eje", seq_len(k)))
  d <- f$d
  bond <- 100 * sum(d[axes]^2) / sum(d^2)
  names(bond) <- paste0("bondad|matriz|", paste(axes, collapse = "-"))
  coo <- as.vector(H)
  names(coo) <- paste0("coordenada|", rep(vars, 2), "|", rep(ejes, each = length(vars)))
  cal <- rowSums(f$quality_cols[, axes, drop = FALSE])
  names(cal) <- paste0("calidad|", vars, "|", paste(axes, collapse = "-"))
  ctr <- as.vector(f$contrib_cols[, axes, drop = FALSE])
  names(ctr) <- paste0("contribucion|", rep(vars, 2), "|", rep(ejes, each = length(vars)))
  nr <- sqrt(rowSums(H^2))
  cosm <- (H %*% t(H)) / outer(nr, nr)
  pares <- utils::combn(length(vars), 2)
  angv <- acos(pmin(1, pmax(-1, cosm[t(pares)]))) * 180 / pi
  names(angv) <- paste0("angulo_variables|", vars[pares[1, ]], "-", vars[pares[2, ]], "|",
                        paste(axes, collapse = "-"))
  ange <- as.vector(acos(pmin(1, pmax(-1, H / nr))) * 180 / pi)
  names(ange) <- paste0("angulo_eje|", rep(vars, 2), "|", rep(ejes, each = length(vars)))
  c(ine, bond, coo, cal, ctr, angv, ange)
}

.boot_types <- c("inercia", "bondad", "coordenada", "calidad", "contribucion",
                 "angulo_variables", "angulo_eje")

#' @export
print.bigarx_boot <- function(x, ...) {
  cat(.method_label(x$fit), "con bootstrap (Nieto et al., 2014)\n")
  cat(sprintf("%d r\u00E9plicas | intervalos %s al %s%% | plano %s\n", x$B,
              if (x$ci == "percentile") "percentiles (boot-q)" else "normales (boot-t)",
              format(100 * x$conf), paste(x$axes, collapse = "-")))
  t <- x$table[x$table$tipo %in% c("inercia", "bondad"), ]
  t$parametro <- ifelse(t$tipo == "bondad", "Bondad de ajuste (%)",
                        paste0("Inercia ", t$elemento, " (%)"))
  cat("\n")
  print(data.frame(parametro = t$parametro, estimacion = round(t$estimacion, 2),
                   ee = round(t$ee, 2), inferior = round(t$inferior, 2),
                   superior = round(t$superior, 2)), row.names = FALSE)
  cat("\nUsa summary(x, type = ...) con: ", paste(.boot_types, collapse = ", "), "\n", sep = "")
  invisible(x)
}

#' Resumen de un biplot con bootstrap
#'
#' @param object Objeto `"bigarx_boot"`.
#' @param type Parámetros a mostrar: uno o varios de `"inercia"`, `"bondad"`,
#'   `"coordenada"`, `"calidad"`, `"contribucion"`, `"angulo_variables"` y
#'   `"angulo_eje"`. Por defecto, todos.
#' @param digits Decimales.
#' @param ... No se usa.
#' @return Un data.frame (de forma invisible) con las filas mostradas.
#' @examples
#' \donttest{
#' b <- bigarx_boot(USArrests, B = 200, method = "hj", seed = 1)
#' summary(b, type = "inercia")
#' }
#' @export
summary.bigarx_boot <- function(object, type = .boot_types, digits = 3, ...) {
  type <- match.arg(type, .boot_types, several.ok = TRUE)
  t <- object$table[object$table$tipo %in% type, ]
  cat(sprintf("Intervalos bootstrap al %s%% (%s, %d r\u00E9plicas)\n", format(100 * object$conf),
              object$ci, object$B))
  for (ty in type) {
    s <- t[t$tipo == ty, -1]
    if (!nrow(s)) next
    cat("\n", ty, "\n", sep = "")
    num <- vapply(s, is.numeric, logical(1))
    s[num] <- lapply(s[num], round, digits)
    print(s, row.names = FALSE)
  }
  invisible(t)
}

#' Gráficos de un biplot con bootstrap
#'
#' `type = "biplot"` dibuja el biplot original con la nube de réplicas de la
#' punta de cada flecha y la elipse que contiene el `conf` de ellas. Los
#' demás tipos dibujan los intervalos de confianza como segmentos.
#'
#' @param x,object Objeto `"bigarx_boot"`.
#' @param type `"biplot"` (por defecto), `"inercia"`, `"calidad"`,
#'   `"contribucion"`, `"angulo_variables"`, `"angulo_eje"` o `"coordenada"`.
#' @param cloud Dibujar las réplicas individuales en el biplot.
#' @param ... Argumentos adicionales para [plot.bigarx()] (en `"biplot"`).
#' @return `x` de forma invisible, o un objeto `ggplot`.
#' @examples
#' \donttest{
#' b <- bigarx_boot(USArrests, B = 200, method = "hj", seed = 1)
#' plot(b)
#' }
#' @export
plot.bigarx_boot <- function(x, type = c("biplot", .boot_types), cloud = TRUE, ...) {
  type <- match.arg(type)
  if (type == "biplot") {
    f <- x$fit
    k <- .col_rescale(f, x$axes, "auto")
    plot(f, axes = x$axes, ...)
    tips <- .boot_tips(x)
    for (v in names(tips)) {
      P <- tips[[v]] * k
      if (cloud) graphics::points(P, pch = ".", col = grDevices::adjustcolor("#C2410C", 0.35))
      el <- .ellipse_pts(P, x$conf)
      graphics::lines(el, col = "#C2410C", lwd = 1.2)
    }
    return(invisible(x))
  }
  t <- x$table[x$table$tipo == type, ]
  lab <- if (type %in% c("coordenada", "contribucion", "angulo_eje")) paste(t$elemento, t$eje) else t$elemento
  ord <- order(t$estimacion)
  t <- t[ord, ]; lab <- lab[ord]
  m <- nrow(t)
  op <- graphics::par(mar = c(4, max(6, max(nchar(lab)) * 0.55), 3, 1))
  on.exit(graphics::par(op))
  graphics::plot(t$estimacion, seq_len(m), xlim = range(c(t$inferior, t$superior)),
                 yaxt = "n", pch = 19, col = "#0B7A72", ylab = "",
                 xlab = .boot_xlab(type),
                 main = sprintf("Intervalos bootstrap al %s%%", format(100 * x$conf)))
  graphics::segments(t$inferior, seq_len(m), t$superior, seq_len(m), col = "#0B7A72", lwd = 2)
  graphics::axis(2, at = seq_len(m), labels = lab, las = 1, cex.axis = 0.75)
  if (type %in% c("angulo_variables", "angulo_eje")) graphics::abline(v = 90, lty = 3, col = "grey50")
  invisible(x)
}

.boot_xlab <- function(type) {
  switch(type, inercia = "Inercia (%)", bondad = "Bondad de ajuste (%)",
         coordenada = "Coordenada", calidad = "Calidad en el plano (%)",
         contribucion = "Contribuci\u00F3n al eje (%)",
         angulo_variables = "\u00C1ngulo entre variables (grados)",
         angulo_eje = "\u00C1ngulo con el eje (grados)")
}

# Puntas de flecha de cada replica: lista con una matriz B x 2 por variable.
.boot_tips <- function(x) {
  vars <- rownames(x$fit$cols)
  ejes <- paste0("Eje", x$axes)
  stats::setNames(lapply(vars, function(v) {
    cbind(x$replicates[, paste0("coordenada|", v, "|", ejes[1])],
          x$replicates[, paste0("coordenada|", v, "|", ejes[2])])
  }), vars)
}

.ellipse_pts <- function(P, conf, npts = 100) {
  mu <- colMeans(P)
  S <- stats::cov(P)
  e <- eigen(S, symmetric = TRUE)
  t <- seq(0, 2 * pi, length.out = npts)
  r <- sqrt(stats::qchisq(conf, 2))
  pts <- e$vectors %*% (sqrt(pmax(e$values, 0)) * r * rbind(cos(t), sin(t)))
  cbind(mu[1] + pts[1, ], mu[2] + pts[2, ])
}

#' @rdname plot.bigarx_boot
#' @exportS3Method ggplot2::autoplot
autoplot.bigarx_boot <- function(object, type = c("biplot", .boot_types), cloud = TRUE, ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("Instala 'ggplot2' para usar autoplot(): install.packages(\"ggplot2\")", call. = FALSE)
  }
  type <- match.arg(type)
  x <- object
  if (type == "biplot") {
    k <- .col_rescale(x$fit, x$axes, "auto")
    g <- autoplot.bigarx(x$fit, axes = x$axes, ...)
    tips <- .boot_tips(x)
    nube <- do.call(rbind, lapply(names(tips), function(v) {
      data.frame(x = tips[[v]][, 1] * k, y = tips[[v]][, 2] * k, variable = v)
    }))
    el <- do.call(rbind, lapply(names(tips), function(v) {
      e <- .ellipse_pts(tips[[v]] * k, x$conf)
      data.frame(x = e[, 1], y = e[, 2], variable = v)
    }))
    if (cloud) {
      g <- g + ggplot2::geom_point(data = nube, ggplot2::aes(x = .data$x, y = .data$y),
                                   colour = "#C2410C", alpha = 0.12, size = 0.5)
    }
    g <- g + ggplot2::geom_path(data = el, ggplot2::aes(x = .data$x, y = .data$y,
                                                        group = .data$variable),
                                colour = "#C2410C", linewidth = 0.6) +
      ggplot2::labs(subtitle = sprintf("Elipses bootstrap al %s%% (%d r\u00E9plicas)",
                                       format(100 * x$conf), x$B))
    return(g)
  }
  t <- x$table[x$table$tipo == type, ]
  t$etiqueta <- if (type %in% c("coordenada", "contribucion", "angulo_eje")) {
    paste(t$elemento, t$eje)
  } else {
    t$elemento
  }
  t$etiqueta <- stats::reorder(t$etiqueta, t$estimacion)
  g <- ggplot2::ggplot(t, ggplot2::aes(y = .data$etiqueta, x = .data$estimacion,
                                       xmin = .data$inferior, xmax = .data$superior)) +
    ggplot2::geom_pointrange(colour = "#0B7A72") +
    ggplot2::labs(x = .boot_xlab(type), y = NULL,
                  title = sprintf("Intervalos bootstrap al %s%%", format(100 * x$conf)))
  if (type %in% c("angulo_variables", "angulo_eje")) {
    g <- g + ggplot2::geom_vline(xintercept = 90, linetype = 3, colour = "grey50")
  }
  g
}