summary 函数对 regsubsets 的输出做了什么？

Question

首先，我认为这个问题是编码问题，而不是统计问题。它几乎肯定会在 Stats.SE.

结束

R 中的 leaps 包对模型 selection 有一个有用的函数，称为 regsubsets，对于任何给定的模型大小，它可以找到产生最小残差的变量平方和。现在我正在阅读 Julian Faraway 的书 Linear Models with R，第二版，。在第 154-5 页，他有一个使用 AIC 模型 selection 的例子。重现示例的完整代码如下运行：

data(state)
statedata = data.frame(state.x77, row.names=state.abb)
require(leaps)
b = regsubsets(Life.Exp~.,data=statedata)
rs = summary(b)
rs$which
AIC = 50*log(rs$rss/50) + (2:8)*2
plot(AIC ~ I(1:7), ylab="AIC", xlab="Number of Predictors")

rs$which 命令生成 regsubsets 函数的输出，并允许您在绘制 AIC 并找到最小化参数数量后 select 模型工商会。但这就是问题所在：虽然键入的示例工作正常，但当我尝试使用此代码并使其适应其他数据时，我遇到了数组中元素数量错误的问题。例如：

require(faraway)
data(odor, package='faraway')
b=regsubsets(odor~temp+gas+pack+
              I(temp^2)+I(gas^2)+I(pack^2)+
              I(temp*gas)+I(temp*pack)+I(gas*pack),data=odor)
rs=summary(b)
rs$which
AIC=50*log(rs$rss/50) + (2:10)*2

产生一条警告消息：

Warning message:
In 50 * log(rs$rss/50) + (2:10) * 2 :
  longer object length is not a multiple of shorter object length

果然length(rs$rss)=8，但是length(2:10)=9。现在我需要做的是模型 selection，这意味着我真的应该为每个模型尺寸设置一个 RSS 值。但是如果我在 AIC 公式中选择 b$rss，它不适用于原始示例！

所以这是我的问题： summary() 对 regsubsets() 函数的输出做了什么？ RSS值的个数不仅不一样，而且值本身也不一样

Answer 1

regsubsets 有一个 nvmax 参数来控制 “要检查的子集的最大大小”。默认情况下，这是 8。如果将其增加到 9 或更高，您的代码就可以工作。

不过请注意，AIC 公式中的 50 是样本量（即 statedata 中的 50 个州）。所以对于你的第二个例子，这应该是 nrow(odor)，所以 15.

Answer 2

好的，所以你知道 regsubsets 的帮助页面说

regsubsets returns an object of class "regsubsets" containing no user-serviceable parts. It is designed to be processed by summary.regsubsets.

你马上就会知道原因了。

regsubsets中的代码调用 Alan Miller 的 Fortran 77 代码进行子集选择。也就是说，它不是我写的，它是用 Fortran 77 写的。我确实理解算法。 1996 年我写 leaps 时（2017 年我又写了一个重要的 modification）我花了足够的时间阅读代码以了解变量在做什么，但 regsubsets 主要遵循代码附带的 Fortran 驱动程序的结构。

regsubsets 对象的 rss 字段具有该名称，因为它在 Fortran 代码中存储了一个名为 RSS 的变量。这个变量是不是最佳模型的残差平方和。 RSS 在设置阶段计算，在任何子集选择完成之前，由子路径 SSLEAPS 注释 'Calculates partial residual sums of squares from an orthogonal reduction from AS75.1.' 也就是说，RSS 描述了模型的 RSS在设计矩阵中没有从左到右拟合的选择：模型只有最左边的变量，然后是最左边的两个变量，依此类推。如果任何人不打算阅读 Fortran，就没有理由需要知道这一点，因此没有记录在案。

summary.regsubsets中的代码从对象的$ress分量中提取输出中的残差平方和，该分量来自Fortran代码中的RESS变量。这是一个数组，其 [i,j] 元素是大小为 i.

的第 j 个最佳模型的残差平方和

所有模型标准都是在 summary.regsubsets 的同一循环中从 $ress 计算得出的，可以编辑为：

    for (i in ll$first:min(ll$last, ll$nvmax)) {
        for (j in 1:nshow) {
            vr <- ll$ress[i, j]/ll$nullrss
            rssvec <- c(rssvec, ll$ress[i, j])
            rsqvec <- c(rsqvec, 1 - vr)
            adjr2vec <- c(adjr2vec, 1 - vr * n1/(n1 + ll$intercept - 
                i))
            cpvec <- c(cpvec, ll$ress[i, j]/sigma2 - (n1 + ll$intercept - 
                2 * i))
            bicvec <- c(bicvec, (n1 + ll$intercept) * log(vr) + 
                i * log(n1 + ll$intercept))
        }
    }

cpvec 为您提供与 AIC 相同的信息，但如果您需要 AIC，则可以直接执行相同的循环并计算它。

summary 函数对 regsubsets 的输出做了什么？

What does the summary function do to the output of regsubsets?

r

package