当因子水平只有一个水平时，将 predict() 与 RcppArmadillo/RcppEigen 结合使用

Question

我有一个关于 predict() 函数与 RcppArmadillo 和 RcppEigen 包一起使用的问题，当因子变量只有一个水平时。我在下面使用 iris 数据集构建了一个 MWE。

我想先使用 RcppArmadillo 估计线性回归模型，然后用它来预测值。我用于估计的数据包含因子变量（具有多个级别且没有 NA）。我想做的预测在一个方面有点不寻常：我想为所有观察使用相同的因子水平来预测值（这个水平在估计中使用的水平）。在下面的示例中，这意味着我想预测 Sepal.Length 就好像所有观察结果都来自 "versicolor" 物种。

这在我使用 lm() 函数估计模型时效果很好，但不适用于 RcppArmadillo::fastLm() 或 RcppEigen::fastLm() 函数。我收到以下错误：Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : contrasts can be applied only to factors with 2 or more levels。如果其中一个因子水平缺失，同样的错误会再次发生。我很清楚为什么至少需要两个级别来进行估计，但我不明白为什么只有一个级别在模型被正确估计后的预测是一个问题。

显而易见的解决方案是使用 lm() 而不是 fastLm()，但遗憾的是这是不可能的，因为我的数据非常大。经过反复试验，我发现了这个肮脏的解决方法：

我堆叠了两个版本的数据：第一个是原始数据（所有因子水平），第二个是修改后的数据（所有观测值具有相同的因子水平）；
我预测此数据集的值（之所以有效，是因为此数据集中存在所有因子水平）；
我只保留修改后的数据子集。

有没有人有比这更好的解决方案，或者至少解释为什么会出现此错误？

library(data.table)

# Loading iris data
iris <- as.data.table(iris)

# Estimating the model
model <-
  RcppArmadillo::fastLm(Sepal.Length ~ 
                 factor(Species)
               + Sepal.Width 
               + Petal.Length 
               + Petal.Width,
               data=iris)

summary(model)

#### 
#### Here is the error I don't understand
#### 

# This is the standard use of the predict function
iris2 <- copy(iris)
iris2[, predict := predict(model, iris2)]

# This is the way I want to use the predict function
# This does not work for some reason
iris2 <- copy(iris)
iris2[, Species := "versicolor"]
iris2[, predict2 := predict(model, iris2)]

#### 
#### This is a dirty work-around
#### 

# Creating a modified dataframe
iris3 <- copy(iris)
iris3[, `:=`(Species = "versicolor",
             data = "Modified data")]

# copying the original dataframe
iris4 <- copy(iris)
iris4[, data := "Original data"]

# Stacking the original data and the modified data
iris5 <- rbind(iris3, iris4)
iris5[, predict := predict(model, iris5)]

# Keeping only the modified data
iris_final <- iris5[data == "Modified data"]

Answer 1

不是解决方案，而是解释为什么会这样。

如果我们检查RcppAramdillo:::predict.fastLm()的源代码，我们会发现它通过

构建预测点的设计矩阵

x <- model.matrix(object$formula, newdata)

另一方面，如果我们检查 stats::predict.lm() 的来源，我们会发现

tt <- terms(object)
## Some source omitted here
Terms <- delete.response(tt)
m <- model.frame(Terms, newdata, na.action = na.action, xlev = object$xlevels)
if (!is.null(cl <- attr(Terms, "dataClasses")))  .checkMFClasses(cl, m)
X <- model.matrix(Terms, m, contrasts.arg = object$contrasts)

这表明 lm() 在其结果中存储了有关预测变量的因子水平和对比的信息，而 fastLm() 在 predict() 调用中重建了该信息：

names(model)
# [1] "coefficients"  "stderr"        "df.residual"   "fitted.values"
# [5] "residuals"     "call"          "intercept"     "formula"      
names(lm_mod) ## Constructed with `lm()` call with same formula
#  [1] "coefficients"  "residuals"     "effects"       "rank"         
#  [5] "fitted.values" "assign"        "qr"            "df.residual"  
#  [9] "contrasts"     "xlevels"       "call"          "terms"        
# [13] "model"

请注意 lm 对象中的 "xlevels" 和 "contrasts" 元素在 fastLm 对象中不存在。不过，从 help("fastLM"):

Linear models should be estimated using the lm function. In some cases, lm.fit may be appropriate.

如果我错了，Dirk 可以纠正我，但我认为 fastLm() 的重点不是提供一个丰富的 OLS 实现来涵盖 stats::lm() 所做的所有用例；我认为它更能说明问题。

如果你的问题是大数据，这就是你不想使用 stats::lm() 的原因，我可以建议像 biglm::biglm() 这样的东西吗？（例如，参见 here）。如果你真的打算使用 RcppArmadillo::fastLm()，你可以做一个较小版本的解决方法；无需复制整个数据，只需为每个未使用的因子水平在预测集中附加一行。

当因子水平只有一个水平时，将 predict() 与 RcppArmadillo/RcppEigen 结合使用

Using predict() with RcppArmadillo/RcppEigen when a factor level has only one level

r

predict

rcpp