R/methods.R

In BiGaRx: Unified Classical and Canonical Biplots

.method_label <- function(x) {
  switch(x$method,
    hj = "HJ-Biplot (Galindo, 1986)",
    gh = "GH-Biplot (Gabriel, 1971)",
    jk = "JK-Biplot (Gabriel, 1971)",
    sqrt = "SQRT-Biplot (Gabriel, 1971)",
    alpha = sprintf("Biplot con alfa = %s (Gabriel, 1971)", format(x$alpha)))
}

.fmt <- function(x, digits = 1) formatC(x, format = "f", digits = digits)

#' @export
print.bigarx <- function(x, ...) {
  cat(.method_label(x), "\n", sep = "")
  cat(sprintf("%d filas x %d columnas | transformaci\u00F3n: %s | rango: %d\n",
              nrow(x$rows), nrow(x$cols), x$preprocessing$transform, x$rank))
  k <- min(3, nrow(x$inertia))
  cat("Inercia por eje (%): ",
      paste(.fmt(x$inertia$inercia[seq_len(k)]), collapse = ", "),
      sprintf("  | plano 1-2: %s%%\n", .fmt(x$inertia$acumulada[2])), sep = "")
  cat(sprintf("Bondad de ajuste plano 1-2 (%%): matriz %s | filas %s | columnas %s\n",
              .fmt(x$gof["matriz"]), .fmt(x$gof["filas"]), .fmt(x$gof["columnas"])))
  if (!is.null(x$imputation)) {
    cat(sprintf("Datos imputados: %d celdas (SVD de rango %d; ver $imputation)\n",
                nrow(x$imputation$missing), x$imputation$rank))
  }
  if (!is.null(x$sup_rows)) cat(sprintf("Filas suplementarias: %d\n", nrow(x$sup_rows)))
  if (!is.null(x$sup_cols)) cat(sprintf("Columnas suplementarias: %d\n", nrow(x$sup_cols)))
  cat("Usa summary() para ver la calidad de representaci\u00F3n.\n")
  invisible(x)
}

#' Resumen de un biplot
#'
#' Reúne la inercia por eje, la bondad de ajuste y la calidad de
#' representación de filas y columnas en los ejes elegidos.
#'
#' @param object Objeto `"bigarx"`.
#' @param axes Ejes para la calidad acumulada (por defecto `c(1, 2)`).
#' @param ... No se usa.
#' @return Objeto `"summary.bigarx"` con las tablas `inertia`, `gof`, `rows`
#'   y `cols`. En `rows` y `cols`, la columna `calidad` es la suma de las
#'   calidades en los ejes elegidos (QLR en MultBiplot).
#' @examples
#' fit <- bigarx(USArrests, method = "hj")
#' summary(fit)
#' @export
summary.bigarx <- function(object, axes = c(1, 2), ...) {
  .check_axes(object, axes)
  tab <- function(Q, C, coords) {
    data.frame(
      coords[, axes, drop = FALSE],
      calidad = rowSums(Q[, axes, drop = FALSE]),
      contribucion = C[, axes, drop = FALSE],
      check.names = FALSE
    )
  }
  out <- list(
    method = .method_label(object),
    axes = axes,
    inertia = object$inertia,
    gof = object$gof,
    rows = tab(object$quality_rows, object$contrib_rows, object$rows),
    cols = tab(object$quality_cols, object$contrib_cols, object$cols)
  )
  class(out) <- "summary.bigarx"
  out
}

#' @export
print.summary.bigarx <- function(x, digits = 2, ...) {
  cat(x$method, "\n\n", sep = "")
  cat("Inercia\n")
  inr <- x$inertia
  inr[-1] <- lapply(inr[-1], round, digits)
  print(utils::head(inr, 6), row.names = FALSE)
  cat("\nBondad de ajuste en el plano 1-2 (%)\n")
  print(round(x$gof, digits))
  lab <- paste(x$axes, collapse = "-")
  cat(sprintf("\nColumnas (calidad en el plano %s, %%)\n", lab))
  print(round(x$cols, digits))
  cat(sprintf("\nFilas (calidad en el plano %s, %%)\n", lab))
  rr <- round(x$rows, digits)
  if (nrow(rr) > 20) {
    print(utils::head(rr[order(-rr$calidad), ], 20))
    cat(sprintf("... %d filas m\u00E1s (ver summary(fit)$rows)\n", nrow(rr) - 20))
  } else {
    print(rr)
  }
  invisible(x)
}

.check_axes <- function(x, axes) {
  if (length(axes) != 2 || any(axes > x$ndim) || any(axes < 1)) {
    stop(sprintf("`axes` debe tener dos ejes entre 1 y %d (usa `ndim` mayor en bigarx()).",
                 x$ndim), call. = FALSE)
  }
}

# Factor para que las flechas ocupen el mismo rango que los puntos cuando
# las dos nubes tienen escalas distintas (GH, JK, alfa). En el HJ ambas
# nubes comparten sistema de referencia y no se reescala.
.col_rescale <- function(x, axes, col_scale) {
  if (is.numeric(col_scale)) return(col_scale)
  if (identical(x$method, "hj") || identical(col_scale, "none")) return(1)
  r <- max(abs(x$rows[, axes]))
  c <- max(sqrt(rowSums(x$cols[, axes, drop = FALSE]^2)))
  0.9 * r / c
}

#' Gráfico de un biplot con gráficos base
#'
#' @param x Objeto `"bigarx"`.
#' @param axes Ejes a representar.
#' @param col_scale `"auto"` (por defecto) reescala las flechas en GH, JK y
#'   alfa para que se vean junto a los puntos; `"none"` las deja como están;
#'   un número las multiplica por ese valor. En el HJ-Biplot `"auto"` no
#'   reescala, porque filas y columnas comparten sistema de referencia.
#' @param min_quality Calidad mínima (en %) en el plano para mostrar los
#'   elementos. Un número se aplica a filas y columnas; un vector de dos,
#'   `c(filas, columnas)`, fija un umbral para cada una. Por ejemplo,
#'   `c(60, 0)` filtra solo las filas. Por defecto 0 (todo). En el JK-Biplot
#'   las columnas tienen baja calidad por construcción (Galindo, 1986), así
#'   que un umbral común puede ocultar variables.
#' @param labels Si es `TRUE`, rotula filas y columnas.
#' @param row_col,arrow_col Colores de los puntos y de las flechas. Si hay
#'   `groups`, los puntos se colorean por grupo.
#' @param main,xlab,ylab Título y nombres de los ejes. Por defecto, el nombre
#'   del método y el porcentaje de inercia de cada eje.
#' @param ... Argumentos adicionales para [graphics::plot()], como `cex`,
#'   `cex.main` o `las`.
#' @return `x`, de forma invisible.
#' @examples
#' fit <- bigarx(USArrests, method = "hj")
#' plot(fit, min_quality = 50)
#' @export
plot.bigarx <- function(x, axes = c(1, 2), col_scale = "auto", min_quality = 0,
                        labels = TRUE, row_col = "#0B7A72", arrow_col = "#C2410C",
                        main = .method_label(x), xlab = NULL, ylab = NULL, ...) {
  .check_axes(x, axes)
  k <- .col_rescale(x, axes, col_scale)
  R <- x$rows[, axes, drop = FALSE]
  C <- x$cols[, axes, drop = FALSE] * k
  mq <- rep_len(min_quality, 2)
  qr <- rowSums(x$quality_rows[, axes, drop = FALSE]) >= mq[1]
  qc <- rowSums(x$quality_cols[, axes, drop = FALSE]) >= mq[2]
  lim <- range(c(R[qr, ], C[qc, ], 0)) * 1.12
  ax <- x$inertia[axes, ]
  plot(R[qr, , drop = FALSE], type = "n", asp = 1, xlim = lim, ylim = lim,
       xlab = if (is.null(xlab)) sprintf("Eje %d (%s%%)", axes[1], .fmt(ax$inercia[1])) else xlab,
       ylab = if (is.null(ylab)) sprintf("Eje %d (%s%%)", axes[2], .fmt(ax$inercia[2])) else ylab,
       main = main, ...)
  graphics::abline(h = 0, v = 0, lty = 3, col = "grey60")
  pc <- if (is.null(x$groups)) row_col else grDevices::hcl.colors(nlevels(x$groups), "Dark 3")[as.integer(x$groups)]
  if (length(pc) > 1) pc <- pc[qr]
  graphics::points(R[qr, , drop = FALSE], pch = 19, cex = 0.7, col = pc)
  graphics::arrows(0, 0, C[qc, 1], C[qc, 2], length = 0.08, col = arrow_col, lwd = 1.5)
  if (labels) {
    graphics::text(R[qr, , drop = FALSE], labels = rownames(R)[qr], pos = 3,
                   cex = 0.65, col = pc)
    graphics::text(C[qc, , drop = FALSE] * 1.07, labels = rownames(C)[qc],
                   cex = 0.75, col = arrow_col, font = 2)
  }
  if (!is.null(x$sup_rows)) {
    S <- x$sup_rows[, axes, drop = FALSE]
    graphics::points(S, pch = 17, col = "grey30")
    if (labels) graphics::text(S, labels = rownames(S), pos = 1, cex = 0.65, col = "grey30")
  }
  if (!is.null(x$sup_cols)) {
    S <- x$sup_cols[, axes, drop = FALSE] * k
    graphics::arrows(0, 0, S[, 1], S[, 2], length = 0.08, col = "grey30", lty = 2)
    if (labels) graphics::text(S * 1.07, labels = rownames(S), cex = 0.7, col = "grey30")
  }
  if (!is.null(x$groups)) {
    graphics::legend("topright", legend = levels(x$groups), pch = 19, bty = "n", cex = 0.75,
                     col = grDevices::hcl.colors(nlevels(x$groups), "Dark 3"))
  }
  invisible(x)
}

#' Datos ordenados de un biplot
#'
#' Devuelve las coordenadas, calidades y contribuciones en formato largo
#' (una fila por elemento), listo para 'ggplot2' u otros paquetes.
#'
#' @param x Objeto `"bigarx"`.
#' @param matrix `"rows"`, `"cols"`, `"both"` (por defecto) o `"eigen"`.
#' @param axes Ejes para la columna `calidad`.
#' @param ... No se usa.
#' @return Un data.frame.
#' @examples
#' fit <- bigarx(USArrests)
#' head(bigarx_tidy(fit))
#' @export
bigarx_tidy <- function(x, ...) UseMethod("bigarx_tidy")

#' @rdname bigarx_tidy
#' @export
bigarx_tidy.default <- function(x, matrix = c("both", "rows", "cols", "eigen"), axes = c(1, 2), ...) {
  matrix <- match.arg(matrix)
  if (matrix == "eigen") return(x$inertia)
  one <- function(coords, Q, C, tipo, sup = FALSE) {
    df <- data.frame(elemento = rownames(coords), tipo = tipo,
                     suplementario = sup, coords, row.names = NULL,
                     check.names = FALSE)
    df$calidad <- if (is.null(Q)) NA_real_ else rowSums(Q[, axes, drop = FALSE])
    df
  }
  parts <- list()
  if (matrix %in% c("both", "rows")) {
    r <- one(x$rows, x$quality_rows, x$contrib_rows, "fila")
    r$grupo <- if (is.null(x$groups)) NA else as.character(x$groups)
    parts$r <- r
    if (!is.null(x$sup_rows)) {
      s <- one(x$sup_rows, NULL, NULL, "fila", TRUE)
      s$grupo <- NA
      parts$sr <- s
    }
  }
  if (matrix %in% c("both", "cols")) {
    cc <- one(x$cols, x$quality_cols, x$contrib_cols, "columna")
    cc$grupo <- NA
    parts$c <- cc
    if (!is.null(x$sup_cols)) {
      s <- one(x$sup_cols, NULL, NULL, "columna", TRUE)
      s$grupo <- NA
      parts$sc <- s
    }
  }
  out <- do.call(rbind, unname(parts))
  rownames(out) <- NULL
  out
}

#' @rdname bigarx_tidy
#' @exportS3Method generics::tidy
tidy.bigarx <- function(x, matrix = c("both", "rows", "cols", "eigen"), axes = c(1, 2), ...) {
  bigarx_tidy(x, matrix = matrix, axes = axes, ...)
}

#' Gráfico de un biplot con 'ggplot2'
#'
#' Devuelve un objeto `ggplot` que se puede modificar con capas, temas y
#' escalas de 'ggplot2' como cualquier otro gráfico.
#'
#' @inheritParams plot.bigarx
#' @param object Objeto `"bigarx"`.
#' @param ... No se usa.
#' @return Un objeto `ggplot`.
#' @examples
#' if (requireNamespace("ggplot2", quietly = TRUE)) {
#'   fit <- bigarx(USArrests)
#'   ggplot2::autoplot(fit) + ggplot2::theme_minimal()
#' }
#' @exportS3Method ggplot2::autoplot
autoplot.bigarx <- function(object, axes = c(1, 2), col_scale = "auto",
                            min_quality = 0, labels = TRUE,
                            row_col = "#0B7A72", arrow_col = "#C2410C", ...) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("Instala 'ggplot2' para usar autoplot(): install.packages(\"ggplot2\")", call. = FALSE)
  }
  .check_axes(object, axes)
  k <- .col_rescale(object, axes, col_scale)
  df <- bigarx_tidy(object, axes = axes)
  ax <- paste0("Eje", axes)
  df$x <- df[[ax[1]]]
  df$y <- df[[ax[2]]]
  isc <- df$tipo == "columna"
  df$x[isc] <- df$x[isc] * k
  df$y[isc] <- df$y[isc] * k
  mq <- rep_len(min_quality, 2)
  umbral <- ifelse(df$tipo == "fila", mq[1], mq[2])
  df <- df[is.na(df$calidad) | df$calidad >= umbral, ]
  rows <- df[df$tipo == "fila", ]
  cols <- df[df$tipo == "columna", ]
  inr <- object$inertia[axes, ]
  g <- ggplot2::ggplot() +
    ggplot2::geom_hline(yintercept = 0, linetype = 3, colour = "grey60") +
    ggplot2::geom_vline(xintercept = 0, linetype = 3, colour = "grey60")
  if (!is.null(object$groups)) {
    g <- g + ggplot2::geom_point(data = rows, ggplot2::aes(x = .data$x, y = .data$y,
                                                           colour = .data$grupo,
                                                           shape = .data$suplementario))
  } else {
    g <- g + ggplot2::geom_point(data = rows, ggplot2::aes(x = .data$x, y = .data$y,
                                                           shape = .data$suplementario),
                                 colour = row_col)
  }
  g <- g +
    ggplot2::geom_segment(data = cols, ggplot2::aes(x = 0, y = 0, xend = .data$x, yend = .data$y,
                                                    linetype = .data$suplementario),
                          colour = arrow_col,
                          arrow = ggplot2::arrow(length = ggplot2::unit(0.18, "cm"))) +
    ggplot2::scale_shape_manual(values = c(`FALSE` = 19, `TRUE` = 17), guide = "none") +
    ggplot2::scale_linetype_manual(values = c(`FALSE` = 1, `TRUE` = 2), guide = "none") +
    ggplot2::coord_equal(clip = "off") +
    ggplot2::expand_limits(x = range(c(df$x, 0)) * 1.15, y = range(c(df$y, 0)) * 1.15) +
    ggplot2::labs(x = sprintf("Eje %d (%s%%)", axes[1], .fmt(inr$inercia[1])),
                  y = sprintf("Eje %d (%s%%)", axes[2], .fmt(inr$inercia[2])),
                  title = .method_label(object), colour = NULL)
  if (labels) {
    g <- g +
      ggplot2::geom_text(data = rows, ggplot2::aes(x = .data$x, y = .data$y, label = .data$elemento),
                         size = 2.6, vjust = -0.7, colour = if (is.null(object$groups)) row_col else "grey25") +
      ggplot2::geom_text(data = cols, ggplot2::aes(x = .data$x * 1.08, y = .data$y * 1.08,
                                                   label = .data$elemento),
                         size = 3, fontface = "bold", colour = arrow_col)
  }
  g
}