Pregunta

Soy ajustar un modelo a los datos de los factores y la predicción. Si el newdata en predict.lm() contiene un único nivel de factor que es desconocido para el modelo, todos de predict.lm() falla y devuelve un error.

¿Hay una buena manera de tener predict.lm() devuelve una predicción para los niveles de los factores del modelo conoce y NA para niveles de los factores desconocidos, en lugar de solamente un error?

Ejemplo código:

foo <- data.frame(response=rnorm(3),predictor=as.factor(c("A","B","C")))
model <- lm(response~predictor,foo)
foo.new <- data.frame(predictor=as.factor(c("A","B","C","D")))
predict(model,newdata=foo.new)

Me gustaría que el último comando para volver tres predicciones "reales" que corresponde a los niveles del factor "A", "B" y "C" y un NA correspondiente al nivel desconocido "D".

¿Fue útil?

Solución

ordenaban y extendió la función por MorgenBall . También se implementa en sperrorest ahora.

Las características adicionales

  • cae niveles de los factores utilizados en lugar de sólo la creación de los valores que faltan para NA.
  • emite un mensaje al usuario de que los niveles de factor han sido retirados
  • comprueba la existencia de variables de factor en test_data y devuelve data.frame original si no están presentes
  • funciona no sólo para lm, y glm sino también para glmmPQL

Nota: La función que se muestra aquí puede cambiar (mejorar) con el tiempo.

#' @title remove_missing_levels
#' @description Accounts for missing factor levels present only in test data
#' but not in train data by setting values to NA
#'
#' @import magrittr
#' @importFrom gdata unmatrix
#' @importFrom stringr str_split
#'
#' @param fit fitted model on training data
#'
#' @param test_data data to make predictions for
#'
#' @return data.frame with matching factor levels to fitted model
#'
#' @keywords internal
#'
#' @export
remove_missing_levels <- function(fit, test_data) {

  # https://stackoverflow.com/a/39495480/4185785

  # drop empty factor levels in test data
  test_data %>%
    droplevels() %>%
    as.data.frame() -> test_data

  # 'fit' object structure of 'lm' and 'glmmPQL' is different so we need to
  # account for it
  if (any(class(fit) == "glmmPQL")) {
    # Obtain factor predictors in the model and their levels
    factors <- (gsub("[-^0-9]|as.factor|\\(|\\)", "",
                     names(unlist(fit$contrasts))))
    # do nothing if no factors are present
    if (length(factors) == 0) {
      return(test_data)
    }

    map(fit$contrasts, function(x) names(unmatrix(x))) %>%
      unlist() -> factor_levels
    factor_levels %>% str_split(":", simplify = TRUE) %>%
      extract(, 1) -> factor_levels

    model_factors <- as.data.frame(cbind(factors, factor_levels))
  } else {
    # Obtain factor predictors in the model and their levels
    factors <- (gsub("[-^0-9]|as.factor|\\(|\\)", "",
                     names(unlist(fit$xlevels))))
    # do nothing if no factors are present
    if (length(factors) == 0) {
      return(test_data)
    }

    factor_levels <- unname(unlist(fit$xlevels))
    model_factors <- as.data.frame(cbind(factors, factor_levels))
  }

  # Select column names in test data that are factor predictors in
  # trained model

  predictors <- names(test_data[names(test_data) %in% factors])

  # For each factor predictor in your data, if the level is not in the model,
  # set the value to NA

  for (i in 1:length(predictors)) {
    found <- test_data[, predictors[i]] %in% model_factors[
      model_factors$factors == predictors[i], ]$factor_levels
    if (any(!found)) {
      # track which variable
      var <- predictors[i]
      # set to NA
      test_data[!found, predictors[i]] <- NA
      # drop empty factor levels in test data
      test_data %>%
        droplevels() -> test_data
      # issue warning to console
      message(sprintf(paste0("Setting missing levels in '%s', only present",
                             " in test data but missing in train data,",
                             " to 'NA'."),
                      var))
    }
  }
  return(test_data)
}

Podemos aplicar esta función para el ejemplo en la cuestión de la siguiente manera:

predict(model,newdata=remove_missing_levels (fit=model, test_data=foo.new))

Al tratar de mejorar esta función, me encontré con el hecho de que SL métodos como lm aprendizaje, etc necesidad glm los mismos niveles en tren y prueba mientras ML métodos de aprendizaje (svm, randomForest) fallará si se eliminan los niveles. Estos métodos necesitan todos los niveles en tren y prueba.

Una solución general es bastante difícil de lograr ya que cada modelo ajustado tiene una forma diferente de almacenar su componente nivel de factor (fit$xlevels para lm y fit$contrasts para glmmPQL). Al menos parece ser consistente a través de los modelos relacionados lm.

Otros consejos

Usted tiene que quitar los niveles extra antes de cualquier cálculo, como:

> id <- which(!(foo.new$predictor %in% levels(foo$predictor)))
> foo.new$predictor[id] <- NA
> predict(model,newdata=foo.new)
         1          2          3          4 
-0.1676941 -0.6454521  0.4524391         NA 

Esta es una manera más general, de hacerlo, se establece todos los niveles que no se producen en los datos originales a NA. Como se mencionó Hadley en los comentarios, que podrían haber optado por incluir esto en la función predict(), pero no lo hicieron

¿Por qué usted tiene que hacer que se vuelve evidente si nos fijamos en el propio cálculo. Internamente, las predicciones se calculan como:

model.matrix(~predictor,data=foo) %*% coef(model)
        [,1]
1 -0.1676941
2 -0.6454521
3  0.4524391

En la parte inferior tiene dos matrices modelo. Se ve que el de foo.new tiene una columna adicional, por lo que no puede utilizar el cálculo matricial más. Si desea utilizar el nuevo conjunto de datos al modelo, también se obtendría un modelo diferente, siendo uno con una variable ficticia extra para el nivel adicional.

> model.matrix(~predictor,data=foo)
  (Intercept) predictorB predictorC
1           1          0          0
2           1          1          0
3           1          0          1
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$predictor
[1] "contr.treatment"

> model.matrix(~predictor,data=foo.new)
  (Intercept) predictorB predictorC predictorD
1           1          0          0          0
2           1          1          0          0
3           1          0          1          0
4           1          0          0          1
attr(,"assign")
[1] 0 1 1 1
attr(,"contrasts")
attr(,"contrasts")$predictor
[1] "contr.treatment"

Puede no sólo eliminar la última columna de la matriz de modelo o bien, porque incluso si lo hace, los otros dos niveles están siendo influenciados. El código para A nivel sería (0,0). Para B esto es (1,0), por este C (0,1) ... y para D es nuevo (0,0)! Así que su modelo podría suponer que A y D están al mismo nivel si sería ingenuo dejar caer la última variable ficticia.

En una parte más teórica: Es posible construir un modelo sin tener todos los niveles. Ahora, como he tratado de explicar antes, ese modelo es Sólo válido para los niveles que se utiliza para crear el modelo. Si se encuentra con nuevos niveles, usted tiene que construir un nuevo modelo para incluir la información adicional. Si no se hace eso, lo único que puede hacer es eliminar las dosis adicionales del conjunto de datos. Pero entonces, básicamente, pierde toda la información que estaba contenida en el mismo, por lo que generalmente no se considera una buena práctica.

Si usted quiere tratar con los niveles que faltan en los datos después de la creación de su modelo lm pero antes de llamar a predecir (dado que no sabemos exactamente qué niveles podrían faltar antemano) aquí es la función que he construido para configurar todos los niveles no en el modelo a NA - la predicción también dará entonces NA y, a continuación, se puede utilizar un método alternativo para predecir estos valores

.

objetivo será su salida de película de película (..., data = trainData)

datos será la trama de datos que desea crear predicciones para

missingLevelsToNA<-function(object,data){

  #Obtain factor predictors in the model and their levels ------------------

  factors<-(gsub("[-^0-9]|as.factor|\\(|\\)", "",names(unlist(object$xlevels))))
  factorLevels<-unname(unlist(object$xlevels))
  modelFactors<-as.data.frame(cbind(factors,factorLevels))


  #Select column names in your data that are factor predictors in your model -----

  predictors<-names(data[names(data) %in% factors])


  #For each factor predictor in your data if the level is not in the model set the value to NA --------------

  for (i in 1:length(predictors)){
    found<-data[,predictors[i]] %in% modelFactors[modelFactors$factors==predictors[i],]$factorLevels
    if (any(!found)) data[!found,predictors[i]]<-NA
  }

  data

}

Parece que usted pueden gustar efectos aleatorios. Mira en algo así como glmer (paquete lme4). Con un modelo bayesiano, obtendrá efectos que se acercan a 0 cuando hay poca información para su uso en la estimación de ellos. Advertencia, sin embargo, que usted tiene que hacer la predicción a sí mismo, en lugar de utilizar predicen ().

Si lo prefiere, puede simplemente hacer las variables ficticias para los niveles que desee incluir en el modelo, por ejemplo, una variable de 0/1 para el lunes, uno para el martes, uno para el miércoles, etc. domingo se elimina automáticamente del modelo si contiene los 0 de. Sin embargo, tener un 1 en la columna de domingo en los otros datos no dejará la etapa de predicción. Se acaba de asumir que el domingo tiene un efecto que de un promedio de los otros días (que pueden o no ser verdad).

Uno de los supuestos de lineales regresiones logísticas / es de poca o ninguna multicolinealidad; por lo que si las variables predictoras son idealmente independientes entre sí, entonces el modelo no tiene que ver toda la variedad posible de los niveles del factor. Un nuevo nivel de factor (D) es un nuevo predictor, y se puede configurar para NA sin afectar la capacidad de predicción de los factores restantes A, B, C. Es por esto que el modelo aún debe ser capaz de hacer predicciones. Pero además de la nueva categoría de D arroja el esquema esperado. Eso es todo el tema. Configuración de correcciones de NA eso.

El paquete lme4 manejará nuevos niveles si se establece el indicador allow.new.levels=TRUE al llamar predict.

Ejemplo: si su día de la semana es el factor en un dow variable y un b_fail resultado categórico, puede ejecutar

M0 <- lmer(b_fail ~ x + (1 | dow), data=df.your.data, family=binomial(link='logit')) M0.preds <- predict(M0, df.new.data, allow.new.levels=TRUE)

Este es un ejemplo con una regresión logística de efectos aleatorios. Por supuesto, se puede realizar la regresión normal ... o la mayoría de los modelos GLM. Si quieres ir más abajo en el camino bayesiano, vistazo a excelente libro Gelman y de la colina de la infraestructura y la Stan .

Una solución rápida y sucia-para pruebas de división, es recodificar valores raras como "otros". Aquí es una aplicación:

rare_to_other <- function(x, fault_factor = 1e6) {
  # dirty dealing with rare levels:
  # recode small cells as "other" before splitting to train/test,
  # assuring that lopsided split occurs with prob < 1/fault_factor
  # (N.b. not fully kosher, but useful for quick and dirty exploratory).

  if (is.factor(x) | is.character(x)) {
    min.cell.size = log(fault_factor, 2) + 1
    xfreq <- sort(table(x), dec = T)
    rare_levels <- names(which(xfreq < min.cell.size))
    if (length(rare_levels) == length(unique(x))) {
      warning("all levels are rare and recorded as other. make sure this is desirable")
    }
    if (length(rare_levels) > 0) {
      message("recoding rare levels")
      if (is.factor(x)) {
        altx <- as.character(x)
        altx[altx %in% rare_levels] <- "other"
        x <- as.factor(altx)
        return(x)
      } else {
        # is.character(x)
        x[x %in% rare_levels] <- "other"
        return(x)
      }
    } else {
      message("no rare levels encountered")
      return(x)
    }
  } else {
    message("x is neither a factor nor a character, doing nothing")
    return(x)
  }
}

Por ejemplo, con data.table, la llamada sería algo como:

dt[, (xcols) := mclapply(.SD, rare_to_other), .SDcol = xcols] # recode rare levels as other

donde xcols es un cualquier subconjunto de colnames(dt).

Licenciado bajo: CC-BY-SA con atribución
No afiliado a StackOverflow
scroll top