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

Question

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

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

当我使用 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()，您可以使用较小版本的解决方法；无需复制整个数据，只需将一行添加到每个未使用的因子水平的预测集中即可。