R/impute.R

In BiGaRx: Unified Classical and Canonical Biplots

#' Imputa datos faltantes con una SVD de rango reducido
#'
#' Rellena las celdas `NA` con la aproximación de rango bajo de la matriz,
#' iterando hasta que los valores imputados se estabilizan. Es el ajuste por
#' mínimos cuadrados alternados descrito por Gabriel y Zamir (1979) y
#' revisado por Galindo (1986, sección 3.1).
#'
#' Algoritmo:
#' 1. Se rellenan los `NA` con la media de su columna.
#' 2. Se aplica la transformación elegida (por ejemplo, estandarizar).
#' 3. Se calcula la SVD y la aproximación de rango `rank`.
#' 4. Se vuelve a la escala original y se reemplazan solo las celdas
#'    faltantes por los valores aproximados.
#' 5. Se repiten los pasos 2 a 4 hasta que el cambio relativo de los valores
#'    imputados sea menor que `tol`.
#'
#' Los valores observados nunca se modifican.
#'
#' @param X Matriz o data.frame numérico con `NA`.
#' @param rank Rango de la aproximación (por defecto 2, el plano del biplot).
#'   Debe ser menor que el número de filas y de columnas.
#' @param transform Transformación usada en cada iteración, como en
#'   [bigarx()]. Conviene usar la misma del análisis posterior.
#' @param maxit Número máximo de iteraciones.
#' @param tol Tolerancia para el cambio relativo de los valores imputados.
#'
#' @return Un objeto de clase `"bigarx_imputation"`: lista con
#' \describe{
#'   \item{`X`}{La matriz completa.}
#'   \item{`missing`}{Matriz de posiciones (fila, columna) imputadas.}
#'   \item{`values`}{data.frame con fila, columna y valor imputado.}
#'   \item{`iterations`, `converged`, `rank`, `transform`}{Detalles del ajuste.}
#' }
#'
#' @references
#' Gabriel, K. R. y Zamir, S. (1979). Lower rank approximation of matrices
#' by least squares with any choice of weights. *Technometrics*, 21(4),
#' 489-498.
#'
#' Galindo, M. P. (1986). Una alternativa de representación simultánea:
#' HJ-Biplot. *Qüestiió*, 10(1), 13-23.
#'
#' @examples
#' X <- as.matrix(USArrests)
#' X[c(3, 10, 25), 2] <- NA
#' imp <- impute_svd(X, rank = 2)
#' imp
#' imp$values
#' @export
impute_svd <- function(X, rank = 2,
                       transform = c("standardize", "center", "double", "grand", "none"),
                       maxit = 500, tol = 1e-6) {
  X <- .as_data_matrix(X, allow_na = TRUE)
  transform <- match.arg(transform)
  miss <- is.na(X)
  if (!is.numeric(rank) || length(rank) != 1 || rank < 1 ||
      rank >= min(dim(X))) {
    stop(sprintf("`rank` debe ser un entero entre 1 y %d.", min(dim(X)) - 1),
         call. = FALSE)
  }
  rank <- as.integer(rank)
  vacias_col <- colSums(!miss) == 0
  if (any(vacias_col)) {
    stop("Columnas sin ning\u00FAn dato: ", paste(colnames(X)[vacias_col], collapse = ", "),
         call. = FALSE)
  }
  vacias_fil <- rowSums(!miss) == 0
  if (any(vacias_fil)) {
    stop("Filas sin ning\u00FAn dato: ", paste(rownames(X)[vacias_fil], collapse = ", "),
         ". Elim\u00EDnalas antes de imputar.", call. = FALSE)
  }

  Xc <- X
  it <- 0L
  converged <- TRUE
  if (any(miss)) {
    cm <- colMeans(X, na.rm = TRUE)
    Xc[miss] <- cm[col(X)[miss]]
    converged <- FALSE
    for (it in seq_len(maxit)) {
      pp <- .preprocess(Xc, transform)
      s <- svd(pp$Z, nu = rank, nv = rank)
      Zh <- s$u %*% (s$d[seq_len(rank)] * t(s$v))
      Xh <- .unpreprocess(Zh, pp$info)
      old <- Xc[miss]
      new <- Xh[miss]
      Xc[miss] <- new
      change <- sqrt(sum((new - old)^2)) / max(sqrt(sum(old^2)), .Machine$double.eps)
      if (change < tol) {
        converged <- TRUE
        break
      }
    }
    if (!converged) {
      warning(sprintf("La imputaci\u00F3n no convergi\u00F3 en %d iteraciones; aumenta `maxit`.", maxit),
              call. = FALSE)
    }
  }
  pos <- which(miss, arr.ind = TRUE)
  out <- list(
    X = Xc,
    missing = pos,
    values = data.frame(fila = rownames(X)[pos[, 1]],
                        columna = colnames(X)[pos[, 2]],
                        valor = Xc[pos],
                        row.names = NULL),
    iterations = it,
    converged = converged,
    rank = rank,
    transform = transform
  )
  class(out) <- "bigarx_imputation"
  out
}

#' @export
print.bigarx_imputation <- function(x, ...) {
  cat(sprintf("Imputaci\u00F3n por SVD de rango %d (transformaci\u00F3n: %s)\n", x$rank, x$transform))
  cat(sprintf("Celdas imputadas: %d | iteraciones: %d | %s\n", nrow(x$missing), x$iterations,
              if (x$converged) "convergi\u00F3" else "NO convergi\u00F3"))
  if (nrow(x$values) > 0) {
    print(utils::head(x$values, 10), row.names = FALSE)
    if (nrow(x$values) > 10) cat(sprintf("... %d m\u00E1s (ver $values)\n", nrow(x$values) - 10))
  }
  invisible(x)
}