R/bigarx.R

In BiGaRx: Unified Classical and Canonical Biplots

#' Biplots clásicos GH, JK y HJ con un motor único
#'
#' Calcula un biplot a partir de la descomposición en valores singulares (SVD)
#' de la matriz preprocesada \eqn{Z = U \Sigma V^T}. Los métodos se diferencian
#' solo en cómo se reparten los valores singulares entre filas y columnas:
#'
#' | Método | Filas | Columnas | Calidad máxima para | Referencia |
#' |---|---|---|---|---|
#' | `"jk"` | \eqn{U\Sigma} | \eqn{V} | filas | Gabriel (1971) |
#' | `"gh"` | \eqn{U} | \eqn{V\Sigma} | columnas | Gabriel (1971) |
#' | `"sqrt"` | \eqn{U\Sigma^{1/2}} | \eqn{V\Sigma^{1/2}} | ninguna | Gabriel (1971) |
#' | `"hj"` | \eqn{U\Sigma} | \eqn{V\Sigma} | filas y columnas | Galindo (1986) |
#'
#' Con `alpha` numérico se obtiene la familia general
#' \eqn{G = U\Sigma^{\alpha}}, \eqn{H = V\Sigma^{1-\alpha}}; `alpha = 1` es el
#' JK, `alpha = 0` el GH y `alpha = 0.5` el SQRT. El HJ-Biplot no pertenece a
#' esa familia: asigna \eqn{\Sigma} a ambos lados, por lo que el producto de
#' marcadores no reproduce los datos, pero filas y columnas logran a la vez la
#' máxima calidad de representación (Galindo, 1986).
#'
#' @param X Matriz o data.frame numérico (individuos en filas, variables en
#'   columnas).
#' @param method `"hj"` (por defecto), `"gh"`, `"jk"` o `"sqrt"`. Se ignora si se
#'   da `alpha`.
#' @param alpha Número entre 0 y 1 para la familia general de Gabriel. `NULL`
#'   usa `method`.
#' @param transform Preprocesamiento: `"standardize"` (por defecto),
#'   `"center"`, `"double"`, `"grand"` o `"none"`. Ver [transformaciones].
#' @param ndim Número de ejes que se guardan en las coordenadas (por defecto 3,
#'   o el rango si es menor). Las medidas de calidad siempre usan todos los
#'   ejes.
#' @param pca_scaling Si es `TRUE`, multiplica las filas por \eqn{\sqrt{n-1}} y
#'   divide las columnas por \eqn{\sqrt{n-1}} (Gabriel, 1971, ecuación 52). Con
#'   `method = "gh"` y `transform = "center"` la longitud de cada flecha
#'   aproxima la desviación típica de la variable y el producto escalar de dos
#'   flechas, su covarianza.
#' @param sup_rows Matriz opcional de filas suplementarias (mismas columnas que
#'   `X`). Se proyectan sin intervenir en el cálculo de los ejes.
#' @param sup_cols Matriz opcional de columnas suplementarias (mismas filas que
#'   `X`).
#' @param groups Factor opcional de longitud `nrow(X)`, solo para colorear o
#'   agrupar en los gráficos. No interviene en el cálculo.
#' @param na_action Qué hacer con los `NA`: `"fail"` (por defecto) detiene el
#'   análisis; `"impute"` los rellena con [impute_svd()] antes de calcular el
#'   biplot. Los valores imputados quedan en `$imputation`.
#' @param impute_rank Rango de la SVD usada para imputar (por defecto 2).
#' @param fix_signs Si es `TRUE` (por defecto), orienta cada eje para que la
#'   carga de mayor valor absoluto en `V` sea positiva. Hace los resultados
#'   reproducibles entre plataformas.
#'
#' @return Un objeto de clase `"bigarx"`: lista con
#' \describe{
#'   \item{`rows`, `cols`}{Coordenadas de filas y columnas en los primeros
#'     `ndim` ejes.}
#'   \item{`sup_rows`, `sup_cols`}{Coordenadas de los elementos suplementarios,
#'     si se dieron.}
#'   \item{`d`}{Todos los valores singulares no nulos.}
#'   \item{`eigen`}{Valores propios \eqn{\sigma_k^2}.}
#'   \item{`inertia`}{Tabla con valor propio, porcentaje de inercia y
#'     porcentaje acumulado por eje.}
#'   \item{`quality_rows`, `quality_cols`}{Calidad de representación (en %) de
#'     cada elemento en cada eje: \eqn{100\, f_{ik}^2 / \sum_h f_{ih}^2}. Suman
#'     100 sobre todos los ejes y se acumulan sumando columnas.}
#'   \item{`contrib_rows`, `contrib_cols`}{Contribución (en %) de cada elemento
#'     a la construcción de cada eje: \eqn{100\, f_{ik}^2 / \sum_i f_{ik}^2}.}
#'   \item{`gof`}{Bondad de ajuste en el plano 1-2 para la matriz, las filas y
#'     las columnas (ver Detalles).}
#'   \item{`method`, `alpha`, `preprocessing`, `groups`, `call`}{Información
#'     del ajuste.}
#' }
#'
#' @details
#' **Bondad de ajuste.** Siguiendo a Gabriel (1971, ecuaciones 35, 66 y 67) y
#' Galindo (1986, tabla 1), para un plano de dimensión 2:
#' * matriz: \eqn{(\sigma_1^2+\sigma_2^2)/\sum\sigma_k^2};
#' * filas con \eqn{G = U\Sigma^{\alpha}}:
#'   \eqn{\sum_{k\le 2}\sigma_k^{4\alpha}/\sum_k\sigma_k^{4\alpha}};
#' * columnas con \eqn{H = V\Sigma^{1-\alpha}}:
#'   \eqn{\sum_{k\le 2}\sigma_k^{4(1-\alpha)}/\sum_k\sigma_k^{4(1-\alpha)}}.
#'
#' Así, en el GH las columnas tienen \eqn{\sum\sigma^4} y las filas solo
#' \eqn{2/r}; en el JK es al revés; y en el HJ ambas tienen la medida en
#' \eqn{\sigma^4}, que es la propiedad que lo define.
#'
#' @references
#' Gabriel, K. R. (1971). The biplot graphic display of matrices with
#' application to principal component analysis. *Biometrika*, 58(3), 453-467.
#'
#' Galindo, M. P. (1986). Una alternativa de representación simultánea:
#' HJ-Biplot. *Qüestiió*, 10(1), 13-23.
#'
#' @examples
#' fit <- bigarx(USArrests, method = "hj")
#' fit
#' summary(fit)
#' plot(fit)
#'
#' # La familia general de Gabriel
#' bigarx(USArrests, alpha = 0.5)
#'
#' # Filas suplementarias
#' bigarx(USArrests[-(1:5), ], sup_rows = USArrests[1:5, ])
#' @export
bigarx <- function(X, method = c("hj", "gh", "jk", "sqrt"), alpha = NULL,
                   transform = c("standardize", "center", "double", "grand", "none"),
                   ndim = 3, pca_scaling = FALSE, sup_rows = NULL,
                   sup_cols = NULL, groups = NULL, fix_signs = TRUE,
                   na_action = c("fail", "impute"), impute_rank = 2) {
  cl <- match.call()
  na_action <- match.arg(na_action)
  X <- .as_data_matrix(X, allow_na = na_action == "impute")
  transform <- match.arg(transform)
  imputation <- NULL
  if (anyNA(X)) {
    imputation <- impute_svd(X, rank = impute_rank, transform = transform)
    X <- imputation$X
  }

  if (!is.null(alpha)) {
    if (!is.numeric(alpha) || length(alpha) != 1 || alpha < 0 || alpha > 1) {
      stop("`alpha` debe ser un n\u00FAmero entre 0 y 1.", call. = FALSE)
    }
    method <- "alpha"
  } else {
    method <- match.arg(method)
    alpha <- switch(method, jk = 1, gh = 0, sqrt = 0.5, hj = NA_real_)
  }

  if (!is.null(groups)) {
    if (length(groups) != nrow(X)) {
      stop("`groups` debe tener una entrada por fila de `X`.", call. = FALSE)
    }
    groups <- as.factor(groups)
  }

  pp <- .preprocess(X, transform)
  Z <- pp$Z
  s <- svd(Z)
  tol <- max(dim(Z)) * max(s$d) * .Machine$double.eps
  r <- sum(s$d > tol)
  if (r < 2) stop("La matriz preprocesada tiene rango menor que 2.", call. = FALSE)
  d <- s$d[seq_len(r)]
  U <- s$u[, seq_len(r), drop = FALSE]
  V <- s$v[, seq_len(r), drop = FALSE]

  if (fix_signs) {
    sg <- apply(V, 2, function(v) sign(v[which.max(abs(v))]))
    U <- sweep(U, 2, sg, "*")
    V <- sweep(V, 2, sg, "*")
  }

  ex <- .exponents(method, alpha)
  Fr <- sweep(U, 2, d^ex$row, "*")
  Gc <- sweep(V, 2, d^ex$col, "*")
  k_pca <- if (pca_scaling) sqrt(nrow(Z) - 1) else 1
  Fr <- Fr * k_pca
  Gc <- Gc / k_pca

  axn <- paste0("Eje", seq_len(r))
  dimnames(Fr) <- list(rownames(X), axn)
  dimnames(Gc) <- list(colnames(X), axn)

  ndim <- max(2L, min(as.integer(ndim), r))
  keep <- seq_len(ndim)

  out <- list(
    rows = Fr[, keep, drop = FALSE],
    cols = Gc[, keep, drop = FALSE],
    sup_rows = NULL,
    sup_cols = NULL,
    d = stats::setNames(d, axn),
    eigen = stats::setNames(d^2, axn),
    inertia = .inertia_table(d),
    quality_rows = .quality(Fr),
    quality_cols = .quality(Gc),
    contrib_rows = .contrib(Fr),
    contrib_cols = .contrib(Gc),
    gof = .gof(d, ex),
    method = method,
    alpha = alpha,
    pca_scaling = pca_scaling,
    preprocessing = pp$info,
    groups = groups,
    rank = r,
    ndim = ndim,
    imputation = imputation,
    .U = U, .V = V,
    call = cl
  )
  class(out) <- "bigarx"

  if (!is.null(sup_rows)) out$sup_rows <- predict.bigarx(out, sup_rows, type = "rows")
  if (!is.null(sup_cols)) out$sup_cols <- predict.bigarx(out, sup_cols, type = "cols")
  out
}

# Exponentes de los valores singulares para filas y columnas.
.exponents <- function(method, alpha) {
  if (identical(method, "hj")) list(row = 1, col = 1) else list(row = alpha, col = 1 - alpha)
}

.inertia_table <- function(d) {
  ev <- d^2
  pct <- 100 * ev / sum(ev)
  data.frame(
    eje = paste0("Eje", seq_along(d)),
    valor_singular = d,
    valor_propio = ev,
    inercia = pct,
    acumulada = cumsum(pct),
    row.names = NULL
  )
}

.quality <- function(M) {
  ss <- rowSums(M^2)
  ss[ss == 0] <- NA_real_
  100 * M^2 / ss
}

.contrib <- function(M) {
  100 * sweep(M^2, 2, colSums(M^2), "/")
}

.gof <- function(d, ex, plane = 2) {
  ratio <- function(p) {
    w <- d^p
    sum(w[seq_len(plane)]) / sum(w)
  }
  c(matriz = 100 * ratio(2),
    filas = 100 * ratio(4 * ex$row),
    columnas = 100 * ratio(4 * ex$col))
}

#' Proyecta filas o columnas suplementarias
#'
#' Sitúa en el biplot nuevas filas (individuos) o columnas (variables) sin
#' recalcular los ejes. Las filas se transforman con los parámetros de la
#' matriz activa; las columnas, con sus propias medias y desviaciones.
#'
#' @param object Objeto `"bigarx"`.
#' @param newdata Matriz o data.frame numérico.
#' @param type `"rows"` (por defecto) o `"cols"`.
#' @param ... No se usa.
#' @return Matriz de coordenadas en los mismos ejes que `object`.
#' @examples
#' fit <- bigarx(USArrests[-(1:5), ])
#' predict(fit, USArrests[1:5, ])
#' @export
predict.bigarx <- function(object, newdata, type = c("rows", "cols"), ...) {
  type <- match.arg(type)
  Xs <- .as_data_matrix(newdata, "newdata")
  ex <- .exponents(object$method, object$alpha)
  d <- object$d
  keep <- seq_len(object$ndim)
  k_pca <- if (isTRUE(object$pca_scaling)) sqrt(nrow(object$.U) - 1) else 1
  if (type == "rows") {
    if (ncol(Xs) != nrow(object$.V)) {
      stop("`newdata` debe tener las mismas columnas que los datos activos.", call. = FALSE)
    }
    Zs <- .transform_sup_rows(Xs, object$preprocessing)
    # Filas: U Sigma^a = Z V Sigma^(a-1)
    out <- sweep(Zs %*% object$.V, 2, d^(ex$row - 1), "*") * k_pca
  } else {
    if (nrow(Xs) != nrow(object$.U)) {
      stop("`newdata` debe tener las mismas filas que los datos activos.", call. = FALSE)
    }
    Zs <- .transform_sup_cols(Xs, object$preprocessing)
    # Columnas: V Sigma^b = Z' U Sigma^(b-1)
    out <- sweep(t(Zs) %*% object$.U, 2, d^(ex$col - 1), "*") / k_pca
  }
  colnames(out) <- names(d)
  out[, keep, drop = FALSE]
}

#' Convierte un objeto `prcomp` en un biplot de BiGaRx
#'
#' Reutiliza la SVD ya calculada por [stats::prcomp()] para obtener cualquiera
#' de los biplots de [bigarx()].
#'
#' @param x Objeto de clase `"prcomp"` con `retx = TRUE`.
#' @param method,alpha,ndim,fix_signs Como en [bigarx()].
#' @param ... No se usa.
#' @return Objeto `"bigarx"`. No incluye elementos suplementarios.
#' @examples
#' p <- prcomp(USArrests, scale. = TRUE)
#' as_bigarx(p, method = "hj")
#' @export
as_bigarx <- function(x, ...) UseMethod("as_bigarx")

#' @rdname as_bigarx
#' @export
as_bigarx.prcomp <- function(x, method = c("hj", "gh", "jk", "sqrt"),
                             alpha = NULL, ndim = 3, fix_signs = TRUE, ...) {
  if (is.null(x$x)) stop("El objeto prcomp debe calcularse con retx = TRUE.", call. = FALSE)
  n <- nrow(x$x)
  Z <- x$x %*% t(x$rotation)
  transform <- if (!isFALSE(x$scale)) "standardize" else if (!isFALSE(x$center)) "center" else "none"
  # La matriz reconstruida ya está preprocesada; se analiza sin transformar
  # y se registra la transformación original.
  fit <- bigarx(Z, method = method, alpha = alpha, transform = "none",
                ndim = ndim, fix_signs = fix_signs)
  fit$preprocessing <- list(transform = transform,
                            col_center = if (isFALSE(x$center)) NULL else x$center,
                            col_scale = if (isFALSE(x$scale)) NULL else x$scale,
                            row_center = NULL, grand = NULL)
  fit$call <- match.call()
  fit
}