飞跃：运行所有可能的线性模型答案

【问题标题】：Leaps: Running all possible linear models飞跃：运行所有可能的线性模型
【发布时间】：2015-08-25 20:34:07
【问题描述】：

我正在尝试使用 R 包 leaps 在单个因变量和超过 50 个可能的预测变量上运行所有可能的回归模型组合——所有可能的大小。每个变量是大约 50 个数值的向量。

这是我尝试过的，但没有成功：

vars <- c("var1", "var2", "var3", ... , "var60")
xs <- as.formula(paste("yvar ~ ", paste(vars, collapse="+")))
model1 <- regsubsets(x=xs, y=yvar, really.big=TRUE, 
    method="exhaustive", data=data)
summary(model1)

我收到的错误是：

Warning message:
In leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax, force.in = force.in,:
15  linear dependencies found

【问题讨论】：

(1) 这是警告，而不是错误。它不应该阻止您继续进行（但您应该尝试弄清楚它的含义）。 (2) 我看到regsubsets 实际上并不适合所有模型，这很好，因为有 2^50 个可能的模型。你知道这是否真的可行吗？ (3) 我们可以有一个可重现的例子吗？
@nongkrong，你能发布你的可复制示例吗（希望足够短以发表评论，否则作为答案发布并稍后删除）？
@BenBolker 当然，我做了dat <- setNames(data.frame(matrix(rnorm(61*50), 50)), c("yvar", paste0("var", 1:60))); model1 = regsubsets(yvar ~ ., y=yvar, really.big=T, method="exhaustive", data=dat)。运行需要一秒钟，我认为y=yvar 是多余的
设置 nvmax=50 以查看 all 子集的相当长的运行时间。默认最大尺寸为 8。

标签： r regression

【解决方案1】：

您正在尝试做的事情有几个问题，一些是计算的，一些是统计的。在 cmets 中已经提到了一些，但我认为最好将它们放在一个答案中。

首先，有 50 个可能的预测变量，有 $2^50$ 个可能的模型（以 10 为基数表示时为 16 位数字）：

> 2^50
[1] 1.1259e+15

虽然算法中有一些不错的捷径，但每个可能的模型仍然需要一些时间来计算（以及一些内存来存储结果）。如果我们假设您的计算机每秒可以计算 100,000 个模型，那么它只需要 356 年以上才能完成（假设您的计算机没有开始交换内存以保留所有结果，这确实会减慢速度）：

> 2^50/100000/60/60/24/365.24
[1] 356.7859

所以我建议坚持使用“最佳”模型，而不是尝试适应所有可能的模型。

另一个问题（导致警告消息的原因）是，如果您只有 50 个观察值（“每个变量是大约 50 个数值的向量”），那么您可以唯一估计的绝对最大斜率数（假设您还拟合截距）为 49，但您试图在最完整的模型中拟合 50 个斜率，因此所涉及的矩阵保证为奇异矩阵。即使限制为 49 个预测变量也不允许估计误差方差。一个常见的经验法则建议您应该尝试为每 10 个观测值拟合不超过 1 个预测变量，以便合理估计斜率和标准误差，这意味着您实际上应该最多只使用 5 个预测变量模型（如果您真的有50 个观察值，您描述中的“大约”意味着缺失值或其他可能使该值更少的问题）。

从最佳子集回归中学到的最好的东西不是“最佳”模型是什么，而是很少有“最佳”模型，并且多个模型可能具有类似的良好属性。更好的方法可能是使用 Lasso/LARS 样式拟合（取决于您真正想要回答的问题）。

未来的几个编程技巧：

了解用于指定模型的y ~ . 方法，这可能使您不必担心使用paste 和as.formula。

请参阅 fortune(77) 参考您的数据名称。

【讨论】：