拟合由 `leap::regsubsets` 选择的多项式回归模型

Question

我使用leaps::regsubsets 执行了线性回归模型的最佳子集选择。然后我选择了具有 14 个预测变量的模型，并使用 coef(model, 14) 给了我以下输出：

structure(c(16.1303774392893, -0.0787496652705482, -0.104929454314886, 
-1.22322411065346, 1.14718778105312, 0.75468065020279, 0.455617836039703, 
0.521951041899427, 0.0124590834643436, -0.0002293804247409, 
1.26667965342874e-07, 1.4002805624594e-06, -9.90560347112683e-07, 
1.8809273394337e-06, 5.48249071436573e-07), .Names = c("(Intercept)", "X1", 
"X2", "poly(X4, 2)1", "poly(X5, 2)1", "poly(X6, 2)2", "poly(X7, 2)2", 
"poly(X9, 2)1", "X10", "X12", "X13", "X14", "X16", "X17", "X18"))

要获得此模型，我需要将其与lm 匹配。因为poly(X, 2)1 是线性的，poly(X, 2)2 是二次的，所以我做到了：

lm(X20 ~ X1 + X2 + X4 + X5 + I(X6 ^ 2) + I(X7 ^ 2) +
         X9 + X10 + X12 + X13 + X14 + X16 + X17 + X18, df)

我想我知道为什么系数不同（请参阅poly() in lm(): difference between raw vs. orthogonal），但为什么它们不给出相同的拟合值和调整后的 R2？

当然，在公式中使用poly(X, 2)[,2] 与regsubsets 输出完全一致。但是只使用第二项正交多项式并指定模型如下是否有效？

lm(X20 ~ X1 + X2 + X4 + X5 + poly(X6, 2)[,2] + poly(X7, 2)[,2] +
   X9 + X10 + X12 + X13 + X14 + X16 + X17 + X18, df) 

有没有比手动指定模型更直接的方法从regsubsets 输出中检索单个模型？

Answer 1

但他们为什么不给出相同的拟合值和调整后的 R2？

如果您不使用 poly 中的所有列，拟合值不一定相同。

set.seed(0)
y <- runif(100)
x <- runif(100)
X <- poly(x, 3)

all.equal(lm(y ~ X)$fitted, lm(y ~ x + I(x ^ 2) + I(x ^ 3))$fitted)
#[1] TRUE

all.equal(lm(y ~ X[, 1:2])$fitted, lm(y ~ x + I(x ^ 2))$fitted)
#[1] TRUE

all.equal(lm(y ~ X - 1)$fitted, lm(y ~ x + I(x ^ 2) + I(x ^ 3) - 1)$fitted)  ## no intercept
#[1] "Mean relative difference: 33.023"

all.equal(lm(y ~ X[, c(1, 3)])$fitted, lm(y ~ x + I(x ^ 3))$fitted)
#[1] "Mean relative difference: 0.03008166"

all.equal(lm(y ~ X[, c(2, 3)])$fitted, lm(y ~ I(x ^ 2) + I(x ^ 3))$fitted)
#[1] "Mean relative difference: 0.03297488"

对于任何k <= degree，我们只有~ 1 + poly(x, degree)[, 1:k] 等效于~ 1 + x + I(x ^ 2) + ... + I(x ^ k)。 （我明确写出了截距，以强调我们必须从 0 次多项式开始。）

(原因与如何生成正交多项式有关。有关详细信息，请参阅How `poly()` generates orthogonal polynomials? How to understand the "coefs" returned?。请注意，在进行 QR 分解时X = QR，因为R 是上三角矩阵（不是对角矩阵），对于任意子集ind，Q[, ind] 与X[, ind] 的列空间不同，除非ind = 1:k。）

所以，I(x ^ 2) 不等于 ploy(x, 2)[, 2]，因此您将得到不同的拟合值（调整后的）R2。

仅使用第二项正交多项式并指定模型如下是否有效？

leaps（或通常任何建模者）从正交多项式中删除列确实是个坏主意。正交多项式是一个类似因子的项，其重要性由 F 统计量（即，将所有列视为一个整体）而不是单个列的 t 统计量确定。

事实上，即使对于原始多项式，省略任何低阶项也不是一个好主意。例如，y ~ 1 + I(x ^ 2) 省略线性项不是一个好主意。这里的一个基本问题是线性移位不是不变的。例如，如果我们将x 转换为x1：

shift <- runif(1)  ## an arbitrary value; can be `mean(x)`
x1 <- x - shift

那么y ~ 1 + I(x ^ 2) 不等同于y ~ 1 + I(x1 ^ 2)，但y ~ 1 + x + I(x ^ 2) 仍然等同于y ~ 1 + x1 + I(x1 ^ 2)。

all.equal(lm(y ~ 1 + I(x ^ 2))$fitted, lm(y ~ 1 + I(x1 ^ 2))$fitted)
#[1] "Mean relative difference: 0.02020984"

all.equal(lm(y ~ 1 + x + I(x ^ 2))$fitted, lm(y ~ 1 + x1 + I(x1 ^ 2))$fitted)
#[1] TRUE

我在R: How to or should I drop an insignificant orthogonal polynomial basis in a linear model? 简要提到了删除列的问题，但我在这里的示例可以让您更深入地了解。

有没有比手动指定模型更直接的方法从regsubsets 输出中检索单个模型？

我不知道；至少两年前我在回答这个帖子时没有弄清楚：Get all models from leaps regsubsets。

---

还有一个问题。假设leaps 返回poly(X, 2)1，我绝对应该在我的模型中保留poly(X, 2)1。但是如果leaps 只返回poly(X, 2)1 怎么办？那么高阶项可以去掉吗？

删除高阶项没有问题（在这种情况下，您最初拟合了一个二次多项式）。正如我所说，我们等价于ind = 1:j，其中j <= degree。但请确保您理解这一点。举以下两个例子。

• 如果leaps 删除poly(x, 5)3 和poly(x, 5)5。您可以安全地删除poly(x, 5)5，但仍建议保留poly(x, 5)3。也就是说，不是拟合 5 阶多项式，而是拟合 4 阶多项式。
• 如果leaps 丢弃poly(x, 6)3 和poly(x, 6)5。由于poly(x, 6)6 没有被删除，因此建议您不要删除任何条款。