如何从 R 中的线性模型中获得 1000 个预测？答案

【问题标题】：How can I get 1000 predictions from a linear model in R?如何从 R 中的线性模型中获得 1000 个预测？
【发布时间】：2021-01-11 00:07:47
【问题描述】：

我想模拟 1000 次线性回归模型的预测值，并根据模型的自变量查看每种类型的汽车有多少次预测 mpg 最高。我使用测试集和训练集是因为我想在训练数据之外评估模型的拟合度。

data(mtcars)
library(caret)
trainingIndex <- createDataPartition(mtcars$mpg, p = 0.8, list = FALSE)
trainingset <- mtcars[trainingIndex,]
testingset <- mtcars[-trainingIndex,]

我创建了一个数据分区以包含一个训练集和一个测试集。现在我有一个测试集和一个训练集，我创建了线性模型

fit <- lm(mpg~., data = trainingset)

现在我有了线性模型，我尝试创建一个引导程序来从模拟中进行预测。我使用 boot_predict 但它给了我一个错误。

library(finalfit)
boot_predict(fit,testingset, type = "response", R = 1000, estimate_name = NULL, 
             confint_sep = "to", condense = TRUE, boot_compare = TRUE, compare_name = NULL,
             comparison = "difference", ref_symbol = "-", digits = c(2,3))

错误：mutate() 输入 term 有问题。 x 无效格式 '%.2f';对字符对象使用格式 %s i 输入term 是(function (x, digits) ...。运行rlang::last_error() 以查看错误发生的位置。另外：警告信息：在 predict.lm(object, newdata, se.fit, scale = 1, type = if (type == : 秩不足拟合的预测可能具有误导性

我不确定这是否是从引导程序获得 1000 预测的最佳方法

【问题讨论】：

该软件包很可能不是为此而设计的。首先看一下小插图表明它适用于分类类型的预测。我认为您仍然可以将其用于回归，但您的问题不是很清楚。
为什么需要训练集和测试集？而且，在引导程序中，如果训练数据中缺少汽车本身，会发生什么......你应该考虑一下这个问题并稍微改一下这个问题。可以编码
感谢您的回答。我想要的是模拟预测值，看看每辆车有多少次在数据集中具有最高预测值。例如，其中一辆汽车是马自达。我想根据自变量模拟 mpg 值的预测 1000 次，看看有多少次马自达在数据集中具有最高的预测值。 .
我编辑我的问题。我想模拟的是预测值，我想根据其他自变量预测每种汽车的 mpg。我想在训练数据之外测试模型。

标签： r linear-regression predict

【解决方案1】：

关于如何使用训练和测试的部分还不清楚，我建议你可以把它整理出来，作为另一个问题。这里面似乎包含了不止一个问题。

我可以尝试解决这个问题：

查看每种类型的汽车有多少次预测值最高 mpg 基于模型的自变量。

对于 1 个引导程序，适合的基本代码如下：

set.seed(111)
da = mtcars[sample(nrow(mtcars),replace=TRUE),]
fit = lm(mpg ~ .,data=da)

要获得排名我们可以做到：

rank(predict(fit,mtcars))

我们将其包装成一个函数并通过许多引导程序对其进行迭代：

bootpred = function(data){

da = da[sample(nrow(da),replace=TRUE),]
fit = lm(mpg ~ .,data=da)
rank(predict(fit,data))

}

predictions = replicate(1000,bootpred(mtcars))

结果是一个矩阵，每列1个bootstrap，每一行，汽车的预测排名：

head(predictions[,1:5],10)
                  [,1] [,2] [,3] [,4] [,5]
Mazda RX4           18   16   12   20   18
Mazda RX4 Wag       14   12   11   16   17
Datsun 710          24   24   27   26   23
Hornet 4 Drive      22   19   20   23   21
Hornet Sportabout   15   13   15   11   15
Valiant             16   11   18   21   20
Duster 360           7   29    5    6    9
Merc 240D           23   23   23   25   25
Merc 230            25   20   19   24   32
Merc 280            20   18   22   17   16

这会告诉你每辆车有多少次价值最高：

rowSums(predictions==1)
          Mazda RX4       Mazda RX4 Wag          Datsun 710      Hornet 4 Drive 
                 18                   0                   0                   0 
  Hornet Sportabout             Valiant          Duster 360           Merc 240D 
                  0                   0                   3                   1 
           Merc 230            Merc 280           Merc 280C          Merc 450SE 
                 12                   1                   0                   1 
         Merc 450SL         Merc 450SLC  Cadillac Fleetwood Lincoln Continental 
                  0                   0                  80                  72 
  Chrysler Imperial            Fiat 128         Honda Civic      Toyota Corolla 
                174                   0                   3                   0 
      Toyota Corona    Dodge Challenger         AMC Javelin          Camaro Z28 
                 18                   3                   0                  10 
   Pontiac Firebird           Fiat X1-9       Porsche 914-2        Lotus Europa 
                  3                   0                   0                   0 
     Ford Pantera L        Ferrari Dino       Maserati Bora          Volvo 142E 
                  0                  10                 591                   0

【讨论】：

【解决方案2】：

这里的问题是boot_predict() 调用broom::tidy()，如果newdata 的行名不等于seq_len(nrow(newdata))，它会添加一个term 列，这会破坏boot_predict() 中的格式化步骤（您'会发现debug() 是你的朋友）。您可能需要考虑向 finalfit here 的开发人员提交问题，并参考此 Stack Overflow 问题。

与此同时，您可以通过更改测试集的行名来解决此问题：

data(mtcars)
library(caret)
set.seed(42)
trainingIndex <- createDataPartition(mtcars$mpg, p = 0.8, list = FALSE)
trainingset <- mtcars[trainingIndex,]
testingset <- mtcars[-trainingIndex,]
rownames(testingset) <- 1:nrow(testingset) ## This is the new step that fixes it
fit <- lm(mpg~., data = trainingset)
library(finalfit)
boot_predict(fit, testingset)

   mpg cyl  disp  hp drat    wt  qsec vs am gear carb
1 21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
2 24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
3 17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
4 13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
                estimate                      Difference
1 21.02 (13.07 to 27.04)                               -
2  21.68 (6.25 to 40.14)  1.42 (-7.30 to 13.93, p=0.720)
3 15.39 (12.09 to 19.82) -4.82 (-10.83 to 2.68, p=0.180)
4  13.27 (3.26 to 24.34) -6.34 (-18.20 to 4.02, p=0.200)

【讨论】：