为什么带子集参数的 lm() 给出的答案与预先设置子集的答案不同？

Question

我在包含多项式的训练数据集上使用 lm()。当我使用 [ ] 预先设置子集时，与在 lm() 函数调用中使用 subset 参数相比，我得到了不同的系数。为什么？

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 package (v2.0.1)}

于 2021-12-26 创建

Answer 1

在您的第二次调用中，看起来 poly() 在子集化之前首先计算。比较下面 model.frame() 的输出：

# first call
model.frame(mpg ~ poly(horsepower, 2), data = auto_train)[1:5,]
#>      mpg poly(horsepower, 2).1 poly(horsepower, 2).2
#> 326 44.3         -0.1037171808          0.1498371034
#> 169 23.0         -0.0372467155         -0.0099055358
#> 131 26.0         -0.0429441840         -0.0000530004
#> 301 23.9         -0.0239526225         -0.0300950106
#> 272 23.2          0.0045347198         -0.0601592336

# second call
model.frame(mpg ~ poly(horsepower, 2), data = Auto, subset = train)[1:5,]
#>      mpg poly(horsepower, 2).1 poly(horsepower, 2).2
#> 326 44.3         -0.0741931315          0.1133792778
#> 169 23.0         -0.0282078693         -0.0034299423
#> 131 26.0         -0.0321494632          0.0039029222
#> 301 23.9         -0.0190108168         -0.0185862638
#> 272 23.2          0.0006971527         -0.0418538153

# same as
model.frame(mpg ~ poly(horsepower, 2), data = Auto)[train,][1:5,]
#>      mpg poly(horsepower, 2).1 poly(horsepower, 2).2
#> 326 44.3         -0.0741931315          0.1133792778
#> 169 23.0         -0.0282078693         -0.0034299423
#> 131 26.0         -0.0321494632          0.0039029222
#> 301 23.9         -0.0190108168         -0.0185862638
#> 272 23.2          0.0006971527         -0.0418538153

Answer 2

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

要为@JonManes 的回答添加更多技术细节，让我们看看 lines 545-553 of the R code where 'model.frame' is defined。

首先我们有（第 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 little更多信息，或 , although in general this stuff is really poorly documented; I'd expect it to be documented here 但它不是 ...）。例如，

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 邮件列表中提出（至少文档中的注释似乎是有序的）。

为什么带子集参数的 lm() 给出的答案与预先设置子集的答案不同？

Why does lm() with the subset argument give a different answer than subsetting in advance?

r

lm

poly