为什么带有子集参数的 lm() 给出的答案与提前子集化给出的答案不同？

Question

我正在使用lm()包含多项式的训练数据集。当我提前进行子集化时，与在函数调用中[ ]使用参数相比，我得到了不同的系数。为什么？subsetlm()

library(ISLR2)

set.seed (1)
train <- sample(392, 196)

auto_train <- Auto[train,]


lm.fit.data <- lm(mpg ~ poly(horsepower, 2), data = auto_train)
summary(lm.fit.data)
#> 
#> Call:
#> lm(formula = mpg ~ poly(horsepower, 2), data = auto_train)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -12.8711  -2.6655  -0.0096   2.0806  16.1063 
#> 
#> Coefficients:
#>                      Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)           23.8745     0.3171  75.298  < 2e-16 ***
#> poly(horsepower, 2)1 -89.3337     4.4389 -20.125  < 2e-16 ***
#> poly(horsepower, 2)2  33.2985     4.4389   7.501 2.25e-12 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 4.439 on 193 degrees of freedom
#> Multiple R-squared:  0.705,  Adjusted R-squared:  0.702 
#> F-statistic: 230.6 on 2 and 193 DF,  p-value: < 2.2e-16


lm.fit.subset <- lm(mpg ~ poly(horsepower, 2), data = Auto, subset = train)
summary(lm.fit.subset)
#> 
#> Call:
#> lm(formula = mpg ~ poly(horsepower, 2), data = Auto, subset = train)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -12.8711  -2.6655  -0.0096   2.0806  16.1063 
#> 
#> Coefficients:
#>                       Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)            23.5496     0.3175  74.182  < 2e-16 ***
#> poly(horsepower, 2)1 -123.5881     6.4587 -19.135  < 2e-16 ***
#> poly(horsepower, 2)2   47.7189     6.3613   7.501 2.25e-12 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 4.439 on 193 degrees of freedom
#> Multiple R-squared:  0.705,  Adjusted R-squared:  0.702 
#> F-statistic: 230.6 on 2 and 193 DF,  p-value: < 2.2e-16

^{由reprex 包(v2.0.1)于 2021 年 12 月 26 日创建}

Answer 1

tl;dr正如其他评论和答案中所建议的，正交多项式基的特征是在考虑子集化之前计算的。

为了向 @JonManes 的答案添加更多技术细节，让我们看看R 代码的第 545-553 行，其中定义了 'model.frame'。

首先我们有（第 545-549 行）

 if(is.null(attr(formula, "predvars"))) {
        for (i in seq_along(varnames))
            predvars[[i+1L]] <- makepredictcall(variables[[i]], vars[[i+1L]])
        attr(formula, "predvars") <- predvars
    }

• 此时，在代码中，formula将不是一个实际的公式（这太简单了！），而是一个terms对象，其中包含有关模型结构的各种对开发人员有用的信息......
• predvars是定义正确重建数据相关基础（如正交多项式和样条曲线）所需信息的属性（请参阅 参考资料 获取?makepredictcall更多信息，或此处，尽管一般来说，这些内容的记录确实很差；我希望它能被记录下来在这里，但它不是......）。例如，

attr(terms(model.frame(mpg ~ poly(horsepower, 2), data = auto_train)),  "predvars")

给出

list(mpg, poly(horsepower, 2, coefs = list(alpha = c(102.612244897959, 
142.498828460405), norm2 = c(1, 196, 277254.530612245, 625100662.205702
))))

这些是多项式的系数，取决于输入数据的分布。

只有在该信息建立之后，在第 553 行，我们才能得到

subset <- eval(substitute(subset), data, env)

换句话说，在确定多项式特征之前，子集参数甚至不会被求值（所有这些信息然后都会传递给内部函数C_modelframe，你真的不想看它......）

请注意，此问题不会导致统计学习环境中的训练集和测试集之间的信息泄漏：多项式的参数化根本不会影响模型的预测（理论上，尽管与通常的浮点结果一样）不太可能完全相同）。在最坏的情况下（如果训练和全套非常不同），它可能会稍微降低数值稳定性。

FWIW（对我来说）这一切都令人惊讶，并且似乎值得在r-devel@r-project.org邮件列表上提出（至少文档中的注释似乎是正确的）。