我在包含多项式的数据训练集上使用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
创建于2021-12-26由reprex包(v2.0.1(
tl;dr正如其他评论和答案中所建议的,在考虑子集之前计算正交多项式基的特征。
为了给@JonManes的答案添加更多技术细节,让我们看看R代码中定义了"model.frame"的第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
了解little位的更多信息,或者在这里,尽管通常情况下,这些东西的文档记录得很差;我希望它在这里记录下来,但它不是…(。例如
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
邮件列表中提出(至少文档中的一条注释似乎是有序的(。
在第二次调用中,看起来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