R/missRanger.R

In domus11move/IMPMUL2: Imputacion Con Arboles De Decisión

#' Imputacion rapida de valores faltantes por "chained random forest"
#'
#' @importFrom stats var reformulate terms.formula predict setNames
#' @importFrom ranger ranger
#'
#' @description Usa el paquete "ranger" [1] para hacer una imputacion de valores faltantes por "chained random forest", ver [2] y [3].
#' Entre el ajuste del modelo iterativo, ofrece la opción de "predictive mean matching".
#' Esto hace que, en primer lugar, evite la imputación con valores que no estén presentes en la base original
#' (como el valor 0.558 en una variable binaria). Ademas, el ppm intenta elevar las varianzas en las distribuciones condicionales
#' resultantes a un nivel realistico y, como tal, permite una imputacion multiple al repetir la llamada al comando.
#' La encadenacion iterativa se detiene cuando \code{maxiter}
#' es alcanzado or si el promedio de la estimacion out-of-bag de rendimiento deja de mejorar. En el ultimo caso, excepto para la primera iteracion, la segunda ultima (o mejor) data imputada es retornada.
#'
#' @param data Un \code{data.frame} o \code{tibble} con valores faltantes para imputar.
#' @param formula Una formula de dos lados especificando las variables a ser imputadas (lado izquierdo) y variables usadas para imputar (lado derecho). Por defecto es . ~ ., en otras palabras, usar todas las variables para imputar todas las variables.
#' Si por ejemplo todas las variables (con datos faltantes) deberian ser imputadas por todas las variables excepto la variable "ID", usar . ~ . - ID. Notar que un "." es evaluado separadamente por cada lado de la fórmula. Tener en cuenta que las variables
#' con datos faltantes deben aparecer en el lado izquierdo si ellos deberian ser usados en el lado derecho.
#' @param pmm.k Numero de candidatos con valores no faltantes para muestrear desde la perspectiva del ppm. 0 para evitar este paso.
#' @param maxiter Maximo numero de iteraciones en cadena.
#' @param seed Semilla.
#' @param verbose Controla cuanta informacion es mostrada en la pantalla. 0 para mostrar nada. 1 (por defecto) muestra un "." por iteracion y variable, 2 para imprimir el error de prediccion OOB por iteracion y variable (1 menos R-cuadrado para regresion).
#' Ademas, si \code{verbose} es positivo, las variables usadas para imputacion son listadas como las variables a ser imputadas (en el orden de imputacion). Esto va a ser util para detectar si algunas variables son inesperadamente salteadas.
#' @param returnOOB Bandera logica. Si TRUE, el error promedio de la estimacion out-of-bag es adicionado al resultado como atributo "oob". Esto no funciona en el caso especial que las variables se imputen univariadamente.
#' @param case.weights Vector con caso de pesos no negativos.
#' @param ... Argumentos pasados a \code{ranger}. Si la base de datos es larga, mejor usar menos arboles (e.g. \code{num.trees = 20}) y/o un bajo valor \code{sample.fraction}.
#' Los siguiente argumentos son, por ejemplo, incompatbles con \code{ranger}: \code{mtry}, \code{write.forest}, \code{probability}, \code{split.select.weights}, \code{dependent.variable.name}, y \code{classification}.
#'
#' @return Un imputado \code{data.frame}.
#'
#' @title missRanger
#' @author ecoteam2019
#'
#' @references
#' [1] Wright, M. N. & Ziegler, A. (2016). ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R. Journal of Statistical Software, in press.
#' http://arxiv.org/abs/1508.04409.
#'
#' [2] Stekhoven, D.J. and Buehlmann, P. (2012). 'MissForest - nonparametric missing value imputation for mixed-type data', Bioinformatics, 28(1) 2012, 112-118.
#' https://doi.org/10.1093/bioinformatics/btr597.
#'
#' [3] Van Buuren, S., Groothuis-Oudshoorn, K. (2011). mice: Multivariate Imputation by Chained Equations in R. Journal of Statistical Software, 45(3), 1-67.
#' http://www.jstatsoft.org/v45/i03/
#' @export
#'
#' @examples
#' irisWithNA <- genNAs(iris, seed = 34)
#' irisImputed <- missRanger(irisWithNA, pmm.k = 3, num.trees = 100)
#' head(irisImputed)
#' head(irisWithNA)
#'
#' \dontrun{
#' # Algoritmo para arboles extra
#' irisImputed_et <- missRanger(irisWithNA, pmm.k = 3, num.trees = 100, splitrule = "extratrees")
#' head(irisImputed_et)
#'
#' # No imputar "Species". Notar: Como esta variable contiene valores faltantes, no será usada
#' # para imputar otras variables.
#' head(irisImputed <- missRanger(irisWithNA, . - Species ~ ., pmm.k = 3, num.trees = 100))
#'
#' # Imputar solo univariadamente.
#' head(irisImputed <- missRanger(irisWithNA, . ~ 1))
#'
#' # Usar "Species" y "Petal.Length" para imputar "Species" y "Petal.Length".
#' head(irisImputed <- missRanger(irisWithNA, Species + Petal.Length ~ Species + Petal.Length,
#'                                pmm.k = 3, num.trees = 100))
#'
#' # Imputacion multiple: Llena los datos 20 veces, corre 20 analisis y agrupa sus resultados.
#' require(mice)
#' filled <- replicate(20, missRanger(irisWithNA, verbose = 0, num.trees = 100, pmm.k = 5),
#'                     simplify = FALSE)
#' models <- lapply(filled, function(x) lm(Sepal.Length ~ ., x))
#' summary(pooled_fit <- pool(models)) # Realistically inflated standard errors and p values
#'
#' # Un data set con valores logicos, numericos, caracteres and factores.
#' n <- 100
#' X <- data.frame(x1 = seq_len(n),
#'                 x2 = log(seq_len(n)),
#'                 x3 = sample(LETTERS[1:3], n, replace = TRUE),
#'                 x4 = factor(sample(LETTERS[1:3], n, replace = TRUE)),
#'                 x5 = seq_len(n) > 50)
#' head(X)
#' X_NA <- generateNA(X, p = seq(0, 0.8, by = .2))
#' head(X_NA)
#'
#' head(X_imp <- missRanger(X_NA))
#' head(X_imp <- missRanger(X_NA, pmm = 3))
#' head(X_imp <- missRanger(X_NA, pmm = 3, verbose = 0))
#' head(X_imp <- missRanger(X_NA, pmm = 3, verbose = 2, returnOOB = TRUE))
#' attr(X_imp, "oob") # OOB prediction errors per column.
#'
#' # El interfaz de la formula
#' head(X_imp <- missRanger(X_NA, x2 ~ x2 + x3, pmm = 3)) # Does not use x3 because of NAs
#' head(X_imp <- missRanger(X_NA, x2 + x3 ~ x2 + x3, pmm = 3))
#' head(X_imp <- missRanger(X_NA, x2 + x3 ~ 1, pmm = 3)) # Univariate imputation
#' }
missRanger <- function(data, formula = . ~ ., pmm.k = 0L, maxiter = 10L, seed = NULL,
                       verbose = 1, returnOOB = FALSE, case.weights = NULL, ...) {
  if (verbose) {
    cat("\nMissing value imputation by random forests\n")
  }

  # 1) CHECKS INICIALES

  stopifnot(is.data.frame(data), dim(data) >= 1L,
            inherits(formula, "formula"),
            length(formula <- as.character(formula)) == 3L,
            is.numeric(pmm.k), length(pmm.k) == 1L, pmm.k >= 0L,
            is.numeric(maxiter), length(maxiter) == 1L, maxiter >= 1L,
            !(c("write.forest", "probability", "split.select.weights", "mtry",
                "dependent.variable.name", "classification") %in% names(list(...))))

  if (!is.null(case.weights)) {
    stopifnot(length(case.weights) == nrow(data), !anyNA(case.weights))
  }

  if (!is.null(seed)) {
    set.seed(seed)
  }

  # 2) SELECCIONA Y CONVIERTE VARIABLES A INPUTAR

  # Extrae lhs y rhs de formula
  relevantVars <- lapply(formula[2:3], function(z) attr(terms.formula(
    reformulate(z), data = data[1, ]), "term.labels"))

  # Elige variables desde lhs con algunos pero no con todos los valores faltantes.
  toImpute <- relevantVars[[1]][vapply(data[, relevantVars[[1]], drop = FALSE],
                                       FUN.VALUE = TRUE, function(z) anyNA(z) && !all(is.na(z)))]

  # Intenta convertir variables especiales a valores numericos/factores para que sea predecido seguramente por ranger.
  converted <- convert(data[, toImpute, drop = FALSE], check = TRUE)
  data[, toImpute] <- converted$X

  # Remueve variables que no pueden ser convertidas con seguridad.
  visitSeq <- setdiff(toImpute, converted$bad)

  if (verbose) {
    cat("\n  Variables to impute:\t\t")
    cat(visitSeq, sep = ", ")
  }

  if (!length(visitSeq)) {
    if (verbose) {
      cat("\n")
    }
    return(data)
  }

  # Obten indicadores de datos faltantes y ordena variables por el numero de valores faltantes
  dataNA <- is.na(data[, visitSeq, drop = FALSE])
  visitSeq <- names(sort(colSums(dataNA)))

  # 3) SELECCIONA VARIABLES USADAS PARA IMPUTAR

  # Variables en las rhs deberian o aparecer en "visitSeq" o no contener ningun valor faltante.
  imputeBy <- relevantVars[[2]][relevantVars[[2]] %in% visitSeq |
                                  !vapply(data[, relevantVars[[2]], drop = FALSE], anyNA, TRUE)]
  completed <- setdiff(imputeBy, visitSeq)

  if (verbose) {
    cat("\n  Variables used to impute:\t")
    cat(imputeBy, sep = ", ")
  }

  # 4) IMPUTACION

  # Inicializacion
  j <- 1L             # iterador
  crit <- TRUE        # criterior en la prediccion del error OOB a seguir iterando
  verboseDigits <- 4L # formateo de los errores de prediccion OOB (si verbose = 2)
  predError <- setNames(rep(1, length(visitSeq)), visitSeq)

  if (verbose >= 2) {
    cat("\n", abbreviate(visitSeq, minlength = verboseDigits + 2L), sep = "\t")
  }

  # Recorriendo iteraciones y variables para imputar
  while (crit && j <= maxiter) {
    if (verbose) {
      cat("\niter ", j, ":\t", sep = "")
    }
    dataLast <- data
    predErrorLast <- predError

    for (v in visitSeq) {
      v.na <- dataNA[, v]

      if (length(completed) == 0L) {
        data[[v]] <- impUniv(data[[v]])
      } else {
        fit <- ranger(formula = reformulate(completed, response = v),
                      data = data[!v.na, union(v, completed), drop = FALSE],
                      case.weights = case.weights[!v.na],
                      ...)
        pred <- predict(fit, data[v.na, completed, drop = FALSE])$predictions
        data[v.na, v] <- if (pmm.k) pmm(xtrain = fit$predictions,
                                        xtest = pred,
                                        ytrain = data[[v]][!v.na],
                                        k = pmm.k) else pred
        predError[[v]] <- fit$prediction.error / (if (fit$treetype == "Regression") var(data[[v]][!v.na]) else 1)

        if (is.nan(predError[[v]])) {
          predError[[v]] <- 0
        }
      }

      if (j == 1L && (v %in% imputeBy)) {
        completed <- union(completed, v)
      }

      if (verbose == 1) {
        cat(".")
      } else if (verbose >= 2) {
        cat(format(round(predError[[v]], verboseDigits), nsmall = verboseDigits), "\t")
      }
    }

    j <- j + 1L
    crit <- mean(predError) < mean(predErrorLast)
  }

  if (verbose) {
    cat("\n")
  }

  if (j == 2L || (j == maxiter && crit)) {
    dataLast <- data
    predErrorLast <- predError
  }

  if (returnOOB) {
    attr(dataLast, "oob") <- predErrorLast
  }

  # Revierte las conversiones
  revert(converted, X = dataLast)
}

#' Una version de \code{typeof} internamente usado por \code{missRanger}.
#'
#' @description Retorna o el valor "numeric" (double o integer), "factor", "character", "logical", "special" (modo numerico, pero cualquier double o integer) o "" (cualquier otro).
#' \code{missRanger} requiere esta informacion para tratar con respuestas tipo no originalmente soportado por \code{ranger}.
#'
#' @author ecoteam2019
#'
#' @param object Algun objeto.
#'
#' @return Un string.
typeof2 <- function(object) {
  if (is.numeric(object)) "numeric" else
    if (is.factor(object)) "factor" else
      if (is.character(object)) "character" else
        if (is.logical(object)) "logical" else
          if (mode(object) == "numeric") "special" else ""
}

#' Conversion de variables no-factores/no-numericas.
#'
#' @description Convierte variables no-factores/no-numericas en un data frame de factores/numericas. Guarda la informacion para su reversion.
#'
#' @author ecoteam2019
#'
#' @param X Un data frame.
#' @param check Si \code{TRUE}, la funcion chequea si las columnas convertidas pueden ser revertidas sin cambios.
#'
#' @return Una lista con los siguientes elementos: \code{X} es la data frame convertida, \code{vars}, \code{types}, \code{classes} son los nombres, tipos y clases de las variables convertidas. Finalmente, \code{bad} los nombres de las variables en \code{X} que deberian ser convertidas pero que no lo fueron.
convert <- function(X, check = FALSE) {
  stopifnot(is.data.frame(X))

  if (!ncol(X)) {
    return(list(X = X, bad = character(0), vars = character(0),
                types = character(0), classes = character(0)))
  }

  types <- vapply(X, typeof2, FUN.VALUE = "")
  bad <- types == "" | if (check) mapply(function(a, b)
    isFALSE(all.equal(a, b)), X, revert(convert(X))) else FALSE
  types <- types[!(types %in% c("numeric", "factor") | bad)]
  vars <- names(types)
  classes <- lapply(X[, vars, drop = FALSE], class)

  X[, vars] <- lapply(X[, vars, drop = FALSE], function(v)
    if (is.character(v) || is.logical(v)) as.factor(v) else as.numeric(v))

  list(X = X, bad = names(X)[bad], vars = vars, types = types, classes = classes)
}

#' Revierte conversion
#'
#' @description Revierte la conversion hecha por \code{convert}.
#'
#' @author ecoteam2019
#'
#' @param con Una lista retornada con \code{convert}.
#' @param X Un data frame con algunas columnas a ser convertidas de vuelta de acuerdo con la information almacenada en \code{converted}.
#'
#' @return Un data frame.
revert <- function(con, X = con$X) {
  stopifnot(c("vars", "types", "classes") %in% names(con), is.data.frame(X))

  if (!length(con$vars)) {
    return(X)
  }

  f <- function(v, ty, cl) {
    switch(ty, logical = as.logical(v), character = as.character(v),
           special = {class(v) <- cl; v}, v)
  }
  X[, con$vars] <- Map(f, X[, con$vars, drop = FALSE], con$types, con$classes)
  X
}