在测试数据中predive.lm（）具有未知因子水平
https://stackoverflow.com/questions/4285214
-
28-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian题
我将模型拟合以对数据和预测进行预测。如果是 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”级别。
解决方案
附加的功能
- 降低未使用的因子级别，而不仅仅是将丢失的值设置为
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 等等。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。正如哈德利（Hadley）在评论中提到的那样，他们本可以选择将其包括在 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
}
听起来您可能会喜欢随机效果。查看类似Glmer（LME4包）之类的东西。使用贝叶斯模型，当估算它们时几乎没有信息时，您将获得该方法0。但是，警告您必须自己进行预测，而不是使用Predivept（）。
另外，您可以简单地为要在模型中包括的级别进行虚拟变量，例如星期一的变量0/1，一个星期二，一个星期三等等。 0是。但是，在其他数据中的周日专栏中有1个不会使预测步骤失败。它只会假设周日的效果是其他几天的平均效果（这可能是或不正确的）。
线性/逻辑回归的假设之一是很少或没有多重共线性。因此，如果预测变量是理想情况下彼此独立的，那么该模型就不需要看到所有可能的因子水平。新的因子水平（D）是一个新的预测因子，可以设置为NA，而不会影响其余因素A，B，c的预测能力。这就是为什么模型仍然应该能够做出预测的原因。但是，增加新级别的D会抛出预期的模式。这就是整个问题。设置NA修复。
这 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的出色书籍和 斯坦 基础设施。
用于分裂测试的快速和划分解决方案是将稀有值重新描述为“其他”。这是一个实现：
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)
}
}
例如，使用数据。
dt[, (xcols) := mclapply(.SD, rare_to_other), .SDcol = xcols] # recode rare levels as other
在哪里 xcols 是任何子集 colnames(dt).
