从r中的重复多元回归模型预测和计算可靠性测试统计量答案

【问题标题】：predicting and calculating reliability test statistics from repeated multiple regression model in r从r中的重复多元回归模型预测和计算可靠性测试统计量
【发布时间】：2015-04-24 10:32:56
【问题描述】：

我想使用 R 中的 lm 函数对我的数据运行 MLR。但是，我正在使用数据拆分交叉验证方法来访问模型的可靠性。我打算使用“样本”函数以 80:20 的比例将数据随机拆分为校准和验证数据集。我想重复说 100 次。在不设置种子的情况下，我相信来自不同样本的模型会有所不同。我在这里遇到了上一篇文章中的函数，它解决了第一部分；

lst <- lapply(1:100, function(repetition) {
mod <- lm(...)   
# Replace this with the code you need to train your model
return(mod)
})
save(lst, file="myfile.RData")

现在的问题是我如何验证这 100 个模型中的每一个，并获得每个模型的可靠性测试统计数据，例如 RSME、ME、Rsquare，并希望获得置信区间。

如果我可以得到包含所有 100 个模型的预测值的数据框形式的输出，那么我应该从那里开始。

有什么帮助吗？

谢谢

【问题讨论】：

平均拟合值并获得回归的 RMSE。采取多数投票（如果 T，T，F，然后 T），如果分类，则获得平均误分类率。 |我是什么？ | R 平方将不再相关 - 哪个是正确的 R 平方值？ R 平方只是单个回归线的拟合度量，您有 100，并且没有一个将是您的最终拟合度量。|如果您有 100 个系数集，则可以凭经验获得它们的置信区间。

标签： r validation lm

【解决方案1】：

快速回顾一下您的问题：您似乎希望将 MLR 模型拟合到大型训练集，然后使用该模型对剩余的验证集进行预测。您希望将此过程重复 100 次，然后您希望能够分析各个模型的特征和预测。

要实现这一点，您可以在模型生成和预测过程中将临时模型信息存储在数据结构中。然后，您可以在之后重新获取和处理所有信息。你没有在描述中提供你自己的数据集，所以我将使用 R 的内置数据集之一来演示它是如何工作的：

> library(car)
> Prestige <- Prestige[,c("prestige","education","income","women")]
> Prestige[,c("income")] <- log2(Prestige[,c("income")])
> head(Prestige,n=5)
                    prestige education      income women
gov.administrators      68.8     13.11 -0.09620212 11.16
general.managers        69.1     12.26 -0.04955335  4.02
accountants             63.4     12.77 -0.11643822 15.70
purchasing.officers     56.8     11.42 -0.11972061  9.11
chemists                73.5     14.62 -0.12368966 11.68

我们首先初始化一些变量。假设您要创建 100 个模型并将 80% 的数据用于训练目的：

nrIterations=100
totalSize <- nrow(Prestige)
trainingSize <- floor(0.80*totalSize)

我们还想创建用于保存中间模型信息的数据结构。在这方面，R 是一种相当通用的高级语言，所以我们将只创建一个列表列表。这意味着每个listentry 可以自己再次保存另一个信息列表。这使我们可以灵活地添加我们需要的任何内容：

trainTestTuple <- list(mode="list",length=nrIterations)

我们现在已准备好创建模型和预测。在每次循环过程中，都会创建一个不同的随机训练子集，同时将剩余数据用于测试目的。接下来，我们将我们的模型拟合到训练数据，然后我们使用这个获得的模型对测试数据进行预测。请注意，我们明确使用自变量来预测因变量：

for(i in 1:nrIterations)
{
  trainIndices <- sample(seq_len(totalSize),size = trainingSize)
  trainSet <- Prestige[trainIndices,]
  testSet <- Prestige[-trainIndices,]

  trainingFit <- lm(prestige ~ education + income + women, data=trainSet)

  # Perform predictions on the testdata
  testingForecast <- predict(trainingFit,newdata=data.frame(education=testSet$education,income=testSet$income,women=testSet$women),interval="confidence",level=0.95)
  # Do whatever else you want to do (compare with actual values, calculate other stuff/metrics ...)
  # ...

  # add your training and testData to a tuple and add it to a list
  tuple <- list(trainingFit,testingForecast) # Add whatever else you need ..
  trainTestTuple[[i]] <- tuple # Add this list to the "list of lists"
}

现在，相关部分：在迭代结束时，我们将拟合模型和样本外预测结果放在一个列表中。此列表包含我们要为当前迭代保存的所有中间信息。最后，我们将此列表放入列表列表中。

现在我们已经完成了建模，我们仍然可以访问我们需要的所有信息，并且可以以任何我们想要的方式处理和分析它。我们来看看模型 50 的建模和预测结果。首先，我们从列表列表中提取模型和预测结果：

> tuple_50 <- trainTestTuple[[50]]
> trainingFit_50 <- tuple_50[[1]]
> testingForecast_50 <- tuple_50[[2]]

我们来看看模型摘要：

> summary(trainingFit_50)

Call:
lm(formula = prestige ~ education + log2(income) + women, data = trainSet)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.9552  -4.6461   0.5016   4.3196  18.4882 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -287.96143   70.39697  -4.091 0.000105 ***
education       4.23426    0.43418   9.752  4.3e-15 ***
log2(income)  155.16246   38.94176   3.984 0.000152 ***
women           0.02506    0.03942   0.636 0.526875    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.308 on 77 degrees of freedom
Multiple R-squared:  0.8072,    Adjusted R-squared:  0.7997 
F-statistic: 107.5 on 3 and 77 DF,  p-value: < 2.2e-16

然后我们显式地获得模型 R-squared 和 RMSE：

> summary(trainingFit_50)$r.squared
[1] 0.8072008
> summary(trainingFit_50)$sigma
[1] 7.308057

我们来看看样本外预测：

> testingForecast_50
        fit       lwr      upr
1  67.38159 63.848326 70.91485
2  74.10724 70.075823 78.13865
3  64.15322 61.284077 67.02236
4  79.61595 75.513602 83.71830
5  63.88237 60.078095 67.68664
6  71.76869 68.388457 75.14893
7  60.99983 57.052282 64.94738
8  82.84507 78.145035 87.54510
9  72.25896 68.874070 75.64384
10 49.19994 45.033546 53.36633
11 48.00888 46.134464 49.88329
12 20.14195  8.196699 32.08720
13 33.76505 27.439318 40.09079
14 24.31853 18.058742 30.57832
15 40.79585 38.329835 43.26187
16 40.35038 37.970858 42.72990
17 38.38186 35.818814 40.94491
18 40.09030 37.739428 42.44117
19 35.81084 33.139461 38.48223
20 43.43717 40.799715 46.07463
21 29.73700 26.317428 33.15657

最后，我们得到了一些关于第二个预测值和相应置信区间的更详细的结果：

> testingPredicted_2ndprediction <- testingForecast_50[2,1]
> testingLowerConfidence_2ndprediction <- testingForecast_50[2,2]
> testingUpperConfidence_2ndprediction <- testingForecast_50[2,3]

编辑重读后，我突然想到您显然不是每次都拆分相同的确切数据集。您在每次迭代期间使用完全不同的数据分区，它们应该以 80/20 的方式进行拆分。但是，只需稍作修改，仍然可以应用相同的解决方案。

另外：出于交叉验证的目的，您可能应该查看cv.lm()

R 帮助中的描述：此函数为多元线性回归提供预测准确性的内部和交叉验证测量。（对于二元逻辑回归，使用 CVbinary 函数。）数据被随机分配到多个“折叠”。依次删除每个折叠，而剩余的数据用于重新拟合回归模型并在删除的观察值处进行预测。

编辑：回复评论。您可以使用您保存的相关性能指标。例如，您可以在trainTestTuple 上使用sapply，以便从每个子列表中提取相关元素。 sapply 会将这些元素作为向量返回，您可以从中计算 mean。这应该有效：

mean_ME <- mean(sapply(trainTestTuple,"[[",2))
mean_MAD <- mean(sapply(trainTestTuple,"[[",3))
mean_MSE <- mean(sapply(trainTestTuple,"[[",4))
mean_RMSE <- mean(sapply(trainTestTuple,"[[",5))
mean_adjRsq <- mean(sapply(trainTestTuple,"[[",6))

另一个小编辑：你的 MAD 的计算看起来很奇怪。仔细检查这是否正是您想要的，这可能是一件好事。

【讨论】：

感谢 Jellen Vermeir 的精彩回应。下面是我的数据dropbox.com/s/l1rul8cbxvygvr6/data.csv?dl=0 结果变量是 Bud。我为需要测试添加了以下几行： valResiduals
出于可读性目的，我在答案底部回复了您的评论。