Подтвердить что ты не робот

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

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

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

Пример кода:

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)

Я бы хотел, чтобы самая последняя команда возвращала три "реальных" прогноза, соответствующих уровням факторов "A", "B" и "C" и a NA, соответствующим неизвестному уровню "D".

4b9b3361

Ответ 1

Подчеркнул и расширил функцию 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.

Ответ 2

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

> 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 являются одинаковым уровнем, если бы наивно отбросить последнюю фиктивную переменную.

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

Ответ 3

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

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

данные - это кадр данных, который вы хотите создать для

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

}

Ответ 4

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

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

Ответ 5

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

Ответ 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. Если вы хотите отправиться дальше по байесовскому пути, посмотрите на отличную книгу Gelman and Hill и Stan.