4

我正在使用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 包于 2021-12-26 创建(v2.0.1)

4

2 回答 2

4

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

要为@JonManes 的答案添加更多技术细节，让我们看一下定义“model.frame”的 R 代码的第 545-553 行。

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

于 2021-12-27T00:45:03.017 回答
2

在您的第二次调用中，它看起来像是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
于 2021-12-27T00:24:40.543 回答