stepwise - stepaic r

predict.lm() con un nivel de factor desconocido en los datos de prueba

Question

Estoy ajustando un modelo para factorizar datos y predecir. Si newdata en predict.lm() contiene un único nivel de factor desconocido para el modelo, todo predict.lm() falla y devuelve un error.

¿Hay una buena manera de que predict.lm() devuelva una predicción para los niveles de factores que conoce el modelo y NA para los niveles de factores desconocidos, en lugar de solo un error?

Código de ejemplo:

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 devuelva tres predicciones "reales" correspondientes a los niveles de factor "A", "B" y "C" y una NA correspondiente al nivel desconocido "D".

Answer 1

Debe eliminar los niveles adicionales 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 forma más general de hacerlo, establecerá todos los niveles que no ocurren en los datos originales a NA. Como Hadley mencionó en los comentarios, podrían haber elegido incluir esto en la función de predict() , pero no lo hicieron

Por qué tienes que hacer eso se vuelve obvio si miras el cálculo en sí mismo. 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 tienes ambas matrices modelo. Verá que el de foo.new tiene una columna adicional, por lo que ya no puede usar el cálculo de la matriz. Si utilizara el nuevo conjunto de datos para modelar, también obtendría un modelo diferente, siendo uno con una variable ficticia adicional para el nivel extra.

> 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"

Tampoco puede borrar la última columna de la matriz del modelo, porque incluso si lo hace, los otros dos niveles seguirán influidos. El código para el nivel A sería (0,0). Para B esto es (1,0), para C esto (0,1) ... ¡y para D es de nuevo (0,0)! Entonces su modelo supondría que A y D tienen el mismo nivel si cayeran ingenuamente la última variable ficticia.

En una parte más teórica: es posible construir un modelo sin tener todos los niveles. Ahora, como traté de explicar antes, ese modelo solo es válido para los niveles que usó al construir el modelo. Si te encuentras con nuevos niveles, debes construir un nuevo modelo para incluir la información adicional. Si no lo hace, lo único que puede hacer es eliminar los niveles adicionales del conjunto de datos. Pero básicamente pierde toda la información que contenía, por lo que generalmente no se considera una buena práctica.

Answer 2

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

Ejemplo: si su factor de día de la semana está en una variable dow y un resultado categórico b_fail , podría 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, puede realizar una regresión regular ... o la mayoría de los modelos GLM. Si desea seguir por el camino de Bayes, mire el excelente libro de Gelman & Hill y la infraestructura de Stan .

Answer 3

Parece que te gustaría efectos aleatorios. Mira algo como glmer (paquete lme4). Con un modelo bayesiano, obtendrás efectos que se aproximan a 0 cuando hay poca información para usar al estimarlos. Sin embargo, advierta que tendrá que hacer la predicción usted mismo, en lugar de usar predicción ().

Alternativamente, puede simplemente hacer variables ficticias para los niveles que desea incluir en el modelo, por ejemplo, una variable 0/1 para el lunes, una para el martes, una para el miércoles, etc. El domingo se eliminará automáticamente del modelo si contiene todo 0''s. Pero tener un 1 en la columna del domingo en los otros datos no fallará en el paso de predicción. Simplemente asumirá que el domingo tiene un efecto que es promedio los otros días (que puede ser cierto o no).

Answer 4

Si quiere lidiar con los niveles faltantes en sus datos después de crear su modelo de lm pero antes de llamar a predecir (dado que no sabemos exactamente qué niveles pueden faltar de antemano), aquí está la función que he creado para configurar todos los niveles, no en el modelo a NA - la predicción también dará NA y luego puede utilizar un método alternativo para predecir estos valores.

objeto será tu salida de lm de lm (..., data = trainData)

los datos serán el marco 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

}

Answer 5

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

Características adicionales

• elimina los niveles de factor no utilizados en lugar de simplemente establecer los valores faltantes en NA .
• emite un mensaje al usuario que indica que los niveles de factor se han descartado
• comprueba la existencia de variables de factor en test_data y devuelve data.frame original si no hay presente
• funciona no solo para lm , glm y sino también para glmmPQL

Nota: La función que se muestra aquí puede cambiar (mejorar) a lo largo del 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://.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 al ejemplo en la pregunta de la siguiente manera:

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

Al tratar de mejorar esta función, encontré el hecho de que los métodos de aprendizaje de SL como lm , glm , etc. necesitan los mismos niveles en tren y prueba, mientras que los métodos de aprendizaje ML ( svm , randomForest ) fallan 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 de nivel de factor ( fit$xlevels para lm y fit$contrasts para glmmPQL ). Al menos parece ser consistente en todos los modelos relacionados con lm .

Answer 6

Una de las suposiciones de Regresiones lineales / logísticas es la colinealidad múltiple o la colinealidad mínima; Entonces, si las variables predictoras son idealmente independientes entre sí, entonces el modelo no necesita ver toda la variedad posible de niveles de factores. Un nuevo nivel de factor (D) es un nuevo predictor, y puede establecerse en NA sin afectar la capacidad de predicción de los factores A, B, C restantes. Esta es la razón por la cual el modelo aún debería ser capaz de hacer predicciones. Pero la adición del nuevo nivel D arroja el esquema esperado. Ese es todo el problema. Establecer arreglos de NA eso.

Answer 7

Una solución rápida y sucia para las pruebas de división, es recodificar valores raros como "otros". Aquí hay una implementació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 así como:

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

donde xcols es cualquier subconjunto de colnames(dt) .