predict.lm()在测试数据中具有未知因子级别

Question

我正在拟合一个模型来分析数据和预测.如果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".

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$xlevelsfor lm和fit$contrastsfor glmmPQL).至少它似乎在lm相关模型中是一致的.

Answer 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

}