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

  • 本文关键字:子集 答案 设置 参数 lm r lm poly
  • 更新时间 :
  • 英文 :


我在包含多项式的数据训练集上使用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

最新更新