Предсказание.lm() с неизвестным уровнем фактора в тестовых данных
Я подгоняю модель к факторам данных и прогнозированию. Если newdata в 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", и NA соответствующий неизвестному уровню "D".
7 ответов
Приведены в порядок и расширены функции MorganBall. В настоящее время это также реализовано в сперроресте.
Дополнительные возможности
- сбрасывает неиспользуемые уровни факторов, а не просто устанавливает пропущенные значения в
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://stackru.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 и т.д. требуют одинаковых уровней в обучении и тестировании во время обученияsvm, randomForest) не удастся, если уровни будут удалены. Эти методы требуют всех уровней в обучении и тестировании.
Общее решение довольно сложно достичь, так как каждая подобранная модель имеет свой способ хранения компонента факторного уровня (fit$xlevels за lm а также fit$contrasts за glmmPQL). По крайней мере, это кажется последовательным через lm связанные модели.
Вы должны удалить дополнительные уровни перед любым расчетом, например:
> 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 один и тот же уровень, если бы он наивно отбрасывал последнюю фиктивную переменную.
На более теоретической части: возможно построить модель, не имея всех уровней. Теперь, как я пытался объяснить ранее, эта модель действительна только для уровней, которые вы использовали при построении модели. Если вы сталкиваетесь с новыми уровнями, вы должны создать новую модель, включающую дополнительную информацию. Если вы этого не сделаете, единственное, что вы можете сделать, - это удалить дополнительные уровни из набора данных. Но тогда вы в основном теряете всю информацию, которая содержалась в нем, поэтому это обычно не считается хорошей практикой.
Если вы хотите разобраться с отсутствующими уровнями в ваших данных после создания вашей модели 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
}
Одно из предположений о линейных / логистических регрессиях заключается в малой или отсутствии мультиколлинеарности; поэтому, если предикторные переменные в идеале не зависят друг от друга, тогда модели не нужно видеть все возможное разнообразие уровней факторов. Новый уровень фактора (D) является новым предиктором, и его можно установить равным NA, не влияя на способность прогнозирования оставшихся факторов A,B,C. Вот почему модель все еще должна быть в состоянии делать прогнозы. Но добавление нового уровня D сбрасывает ожидаемую схему. Вот и вся проблема. Настройка NA исправляет это.
Похоже, вам могут понравиться случайные эффекты. Посмотрите на что-то вроде glmer (пакет lme4). С байесовской моделью вы получите эффекты, которые приближаются к 0, когда для оценки их недостаточно информации. Однако, предупреждаю, что вам придется делать прогнозирование самостоятельно, а не с помощью Foret().
Кроме того, вы можете просто создать фиктивные переменные для уровней, которые вы хотите включить в модель, например, переменную 0/1 для понедельника, одну для вторника, одну для среды и т. Д. Воскресенье будет автоматически удалено из модели, если она содержит все 0'. Но наличие 1 в столбце Sunday в других данных не подведет шаг прогнозирования. Это будет просто предполагать, что воскресенье имеет эффект, который является средним в другие дни (что может быть или не быть правдой).
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 & Hill и инфраструктуру Stan.
Быстрое и грязное решение для сплит-тестирования состоит в том, чтобы перекодировать редкие значения как "другие". Вот реализация:
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),