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

Posted 2023-02-16

【中文标题】predict.lm() 在测试数据中具有未知因子水平

【英文标题】：predict.lm() with an unknown factor level in test data

【发布时间】：2011-05-16 04:09:59

【问题描述】：

我正在拟合一个模型来分解数据并进行预测。如果predict.lm() 中的newdata 包含模型未知的单个因子水平，则predict.lm() 的all 将失败并返回错误。

有没有一种好方法可以让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”的三个“真实”预测以及对应于未知水平“D”的NA。

【参考方案1】：

拆分测试的快速而简单的解决方案是将稀有值重新编码为“其他”。这是一个实现：

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) 的任意子集。

【参考方案2】：

如果您在调用predict 时设置标志allow.new.levels=TRUE，lme4 包将处理新级别。

示例：如果您的星期几因素在变量 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 基础架构。

【讨论】：

这听起来很有帮助。您能否编辑您的答案以包含可运行代码？如果我只是将lm 更改为lmer，R 会抱怨我没有指定任何随机效果。

【参考方案3】：

通过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://***.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))

在尝试改进此功能时，我发现lm、glm 等 SL 学习方法在训练和测试中需要相同的水平，而 ML 学习方法（svm、randomForest）如果级别被删除，则失败。这些方法需要所有级别的训练和测试。

一般的解决方案很难实现，因为每个拟合模型都有不同的存储因子水平分量的方式（fit$xlevels 用于 lm，fit$contrasts 用于 glmmPQL）。至少在lm 相关模型中似乎是一致的。

【讨论】：

虽然您编写了一个非常方便的函数，但我注意到此代码不适用于变量名称以数字结尾的数据集。

sperrorest 现在已被 mlr 纳入。这个方法在 mlr 的什么地方？

@Muno 在makeLearner() 中使用fix.factor.prediction，例如makeLearner("regr.lm", fix.factors.prediction = TRUE)

【参考方案4】：

如果您想在创建 lm 模型之后但在调用 predict 之前处理数据中缺失的级别（假设我们事先并不确切知道哪些级别可能会丢失）这里是我为设置所有级别而构建的函数不在模型中为 NA - 预测也会给出 NA，然后您可以使用替代方法来预测这些值。

object 将是 lm(...,data=trainData)

的 lm 输出

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

【讨论】：

感谢您提供此功能。我认为 predict() 应该在内部执行此操作，并发送警告，而不是完全失败。

【参考方案5】：

您必须在任何计算之前删除额外的级别，例如：

> 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 在 cmets 中提到的，他们本可以选择将其包含在 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 是同一级别，如果它会天真地删除最后一个虚拟变量。

在更理论的部分：可以在没有所有级别的情况下构建模型。现在，正如我之前试图解释的那样，该模型仅对您在构建模型时使用的关卡有效。如果你遇到新的关卡，你必须建立一个新的模型来包含额外的信息。如果你不这样做，你唯一能做的就是从数据集中删除额外的级别。但是，您基本上会丢失其中包含的所有信息，因此通常不被认为是好的做法。

【讨论】：

我不完全确定为什么这在理论上是不可能的......如果（如果！我应该提前指定这个）我使用 contr.treatment 模型矩阵，其他因素水平应该不受影响，他们应该吗？

非常感谢您的解释，但我仍然不明白...是的，当然，3 级因子和 4 级因子不携带相同的信息。但是为什么不应该对已经看到的因子水平做出预测呢？是的，4 水平因子的模型矩阵不适合 3 水平因子的系数，但可以简单地删除与未知水平对应的列。我的应用程序根据星期几预测销售额 - 即使商店在周日从未营业，难道不能预测周一的销售额（我们已经看到）吗？

@Stephan ：当然。但是，如果您在周日有销售数据，但您没有将其带入原始模型中，则不会。因为在星期天销售的商店在星期一的销售量与在星期天没有营业的商店的销售量不同。因此，模型和新数据是不兼容的，因为它们没有完全谈论同一件事。统计学就是这样：它是数学，不是一般理论。

@Stephan：添加了另一个角度来看待它，也许这样可以解决问题。

我认为你在这里偏离了基础——在很多情况下，你可能事先不知道所有可能的值，当遇到新值时返回缺失值是一个明智的选择。模型矩阵具有不同表示的事实是一个红鲱鱼。

【参考方案6】：

线性/逻辑回归的假设之一是很少或没有多重共线性；因此，如果预测变量理想地相互独立，则模型不需要查看所有可能的因子水平变化。一个新的因子水平 (D) 是一个新的预测因子，可以设置为 NA，而不影响其余因子 A、B、C 的预测能力。这就是模型应该仍然能够做出预测的原因。但是添加新的 D 级会抛出预期的模式。这就是整个问题。设置 NA 可以解决这个问题。

【参考方案7】：

听起来您可能喜欢随机效果。研究一下 glmer （lme4 包）之类的东西。使用贝叶斯模型，当估计它们时使用的信息很少时，您将获得接近 0 的效果。但是，警告您必须自己进行预测，而不是使用 predict()。

或者，您可以简单地为要包含在模型中的级别创建虚拟变量，例如变量 0/1 代表星期一，一个代表星期二，一个代表星期三，等等。如果模型中包含全 0，星期日将自动从模型中删除。但是在其他数据的星期日列中具有 1 不会使预测步骤失败。它只是假设周日的影响与其他日子的平均值相同（这可能是真的，也可能不是）。

【讨论】：

谢谢，我发现这个答案对我的简历问题很有帮助：stats.stackexchange.com/questions/172696/…