Предсказание.lm () с неизвестным уровнем фактора в тестовых данных

Question

Я подгоняю модель к факторам данных и прогнозированию.Если newdata в predict.lm() содержит один факторный уровень, который неизвестен модели, all из predict.lm() завершается ошибкой и возвращает ошибку.

Есть ли хороший способpredict.lm()

Пример кода:

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)

Мне бы очень хотелось, чтобы *1008* возвращал прогноз для тех уровней факторов, которые известны модели, и NA для неизвестных уровней факторов.последняя команда для возврата трех «реальных» предсказаний, соответствующих уровням факторов «A», «B» и «C», и NA, соответствующих неизвестному уровню «D».

Answer 1

Вы должны удалить дополнительные уровни перед любым вычислением, например:

> 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 

Это более общий способ, он установит все уровни, которые не встречаются в исходных данных, в NA,Как упомянул Хэдли в комментариях, они могли бы включить это в функцию predict(), но они не

Почему это нужно сделать, становится очевидным, если вы посмотрите на сам расчет.Внутренне прогнозы рассчитываются следующим образом:

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

Внизу у вас есть обе матрицы моделей.Вы видите, что у foo.new есть дополнительный столбец, поэтому вы больше не можете использовать матричный расчет.Если вы будете использовать новый набор данных для моделирования, вы также получите другую модель, представляющую собой дополнительную фиктивную переменную для дополнительного уровня.

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

Вы не можете просто удалить последний столбец изМатрица модели тоже, потому что, даже если вы сделаете это, оба других уровня все еще влияют.Код для уровня A будет (0,0).Для B это (1,0), для C это (0,1) ... и для D это снова (0,0)!Таким образом, ваша модель предполагает, что A и D - это один и тот же уровень, если она наивно опускает последнюю фиктивную переменную.

С более теоретической точки зрения: можно построить модель безуровни.Теперь, как я пытался объяснить ранее, эта модель только действительна для уровней, которые вы использовали при построении модели.Если вы сталкиваетесь с новыми уровнями, вы должны создать новую модель, включающую дополнительную информацию.Если вы этого не сделаете, единственное, что вы можете сделать, - это удалить дополнительные уровни из набора данных.Но тогда вы в основном теряете всю информацию, которая содержалась в нем, поэтому это обычно не считается хорошей практикой.

Answer 2

Убрал и расширил функцию на MorgenBall . Это теперь также реализовано в sperrorest .

Дополнительные функции

• сбрасывает неиспользуемые уровни факторов, а не просто устанавливает пропущенные значения на NA.
• выдает пользователю сообщение о том, что уровни факторов были сброшены
• проверяет наличие факторных переменных в test_data и возвращает исходный data.frame, если не присутствует
• работает не только для lm, glm, но и для glmmPQL

Примечание. Показанная здесь функция может со временем меняться (улучшаться).

#' @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)
}

Мы можем применить эту функцию к примеру в вопросе следующим образом:

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

Пытаясь улучшить эту функцию, я столкнулся с тем, что методы обучения SL, такие как lm, glm и т. Д., Требуют одинаковых уровней в обучении и тестировании, в то время как методы обучения ML (svm, randomForest) потерпеть неудачу, если уровни удалены. Эти методы требуют всех уровней в обучении и тестировании.

Общее решение довольно сложно достичь, поскольку каждая встроенная модель имеет свой способ хранения компонента уровня фактора (fit$xlevels для lm и fit$contrasts для glmmPQL). По крайней мере, кажется, что это согласованно для lm связанных моделей.

Answer 3

Если вы хотите разобраться с отсутствующими уровнями в ваших данных после создания вашей модели lm, но перед вызовом предиката (учитывая, что мы не знаем точно, какие уровни могут быть пропущены заранее), вот функция, которую я построил, чтобы установить все уровнине в модели для NA - прогноз также даст NA, и затем вы можете использовать альтернативный метод для прогнозирования этих значений.

object будет вашим выводом lm из lm (..., data = trainData)

data будет фреймом данных, для которого вы хотите создать прогнозы для

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 4

Одним из предположений о линейных / логистических регрессиях является небольшая или отсутствие мультиколлинеарности;поэтому, если предикторные переменные в идеале не зависят друг от друга, модели не нужно видеть все возможное разнообразие уровней факторов.Новый уровень фактора (D) является новым предиктором, и его можно установить равным NA, не влияя на способность прогнозирования остальных факторов A, B, C.Вот почему модель все еще должна быть в состоянии делать прогнозы.Но добавление нового уровня D сбрасывает ожидаемую схему.Вот и вся проблема.Настройка NA исправляет это.

Answer 5

Похоже, вам могут понравиться случайные эффекты. Посмотрите на что-то вроде glmer (пакет lme4). С байесовской моделью вы получите эффекты, которые приближаются к 0, когда для оценки их мало информации. Однако, предупреждаю, что вам придется делать прогнозирование самостоятельно, а не с помощью Foret ().

Кроме того, вы можете просто создать фиктивные переменные для уровней, которые вы хотите включить в модель, например, переменная 0/1 для понедельника, одна для вторника, одна для среды и т. д. Воскресенье будет автоматически удалено из модели, если она содержит все 0. Но наличие 1 в столбце Sunday в других данных не подведет шаг прогнозирования. Это будет просто предполагать, что воскресенье имеет эффект, который является средним в другие дни (что может быть или не быть правдой).

Answer 6

Пакет lme4 будет обрабатывать новые уровни, если вы установите флаг allow.new.levels=TRUE при вызове predict.

Пример: если ваш фактор дня недели находится в переменной dow и категориальномрезультат b_fail, вы можете запустить

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)

Это пример с логистической регрессией случайных эффектов.Конечно, вы можете выполнять регулярную регрессию ... или большинство моделей GLM.Если вы хотите идти дальше по Байесовскому пути, посмотрите на превосходную книгу Гельмана и Хилла и инфраструктуру Stan .

Answer 7

Быстрое и грязное решение для сплит-тестирования состоит в том, чтобы перекодировать редкие значения как «другие». Вот реализация:

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)
  }
}

Например, для data.table вызов будет выглядеть примерно так:

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

, где xcols - это любое подмножество colnames(dt).