拟合由“leaps::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)

我想我知道为什么系数不同（参见），但为什么它们不给出相同的拟合值和调整后的 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

but why don't they give the same fitted values and adjusted 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次多项式开始。）

(原因与正交多项式的生成方式有关。有关详细信息，请参阅。请注意，在进行 QR 因式分解时 X = QR，因为 R 是上三角矩阵（不是对角矩阵），Q[, ind] 不会有与 X[, ind] 相同的列 space 对于任意子集 ind，除非 ind = 1:k.)

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

is it valid to use only second term orthogonal polynomial and specify the model as follows?

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

我在简要提到了删除列的问题，但我这里的示例可以让您更深入地了解。

Is there more direct way to retrieve single model from regsubsets output than specifying the model by hand?

我不知道；至少我在大约 2 年前回答这个话题时没有弄清楚：.

One remaining question though. Assuming that leaps returns poly(X, 2)1 I should definitely retain poly(X, 2)1 in my model. But what if only poly(X, 2)1 is returned by leaps? Can higher order term can be dropped then?

删除高阶项没有问题（在这种情况下，您最初拟合的是二次多项式）。正如我所说，我们对 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 未被删除，因此建议您完全不删除任何条款。

拟合由“leaps::regsubsets”选择的多项式回归模型

Fitting a polynomial regression model selected by `leaps::regsubsets`

regression

r

linear-regression

lm

polynomials