lm() 和 predict.lm() 的奇怪行为取决于显式命名空间访问器的使用

Bizarre behaviour of lm() and predict.lm() depending on use of explicit namespace accessor

我对 R 中 lm 函数和关联的 predict.lm 函数的一些令人不安的行为很感兴趣。splines 基础包提供函数 bs 来生成b 样条展开,然后可以使用 lm 来拟合样条模型,这是一种通用的线性模型拟合函数。

lmpredict.lm 函数有很多利用公式和术语的内置便利。如果对 bs() 的调用嵌套在 lm 调用中,则用户可以向 predict 提供单变量数据,该数据将自动扩展为适当的 b 样条基。然后将像往常一样预测这个扩展的数据矩阵。

library(splines)

x <- sort(runif(50, 0, 10))
y <- x^2

splineModel <- lm(y ~ bs(x, y, degree = 3, knots = c(3, 6)))

newData <- data.frame(x = 4)
prediction <- predict(splineModel, newData) # 16

plot(x, y)
lines(x, splineModel$fitted.values, col = 'blue3')
points(newData$x, prediction, pch = 3, cex = 3, col = 'red3')
legend("topleft", legend = c("Data", "Fitted Values", "Predicted Value"),
       pch = c(1, NA, 3), col = c('black', 'blue3', 'red3'), lty = c(NA, 1, NA))

如我们所见,这非常有效:

当使用 :: 运算符明确指示 bs 函数是从 splines 包的名称空间导出时,就会出现奇怪的情况。以下代码片段除此更改外是相同的:

library(splines)

x <- sort(runif(50, 0, 10))
y <- x^2

splineModel <- lm(y ~ splines::bs(x, y, degree = 3, knots = c(3, 6)))

newData <- data.frame(x = 4) 
prediction <- predict(splineModel, newData) # 6.40171

plot(x, y)
lines(x, splineModel$fitted.values, col = 'blue3')
points(newData$x, prediction, pch = 3, cex = 3, col = 'red3')
legend("topleft", legend = c("Data", "Fitted Values", "Predicted Value"),
       pch = c(1, NA, 3), col = c('black', 'blue3', 'red3'), lty = c(NA, 1, NA))

如果一开始就没有使用 library 附加 splines 包,则在第二个片段中会产生完全相同的结果。我想不出另一种情况,即在已加载的包上使用 :: 运算符会改变程序行为。

使用 splines 中的其他函数也会出现相同的行为,例如自然样条基础实现 ns。有趣的是,在这两种情况下 "y hat" 或拟合值都是合理的并且彼此匹配。据我所知,除了属性名称外,拟合模型对象是相同的。

我无法确定此行为的来源。虽然这可能读起来像错误报告,但我的问题

  1. 为什么会这样?我一直在努力跟进 predict.lm 但无法确定分歧发生的位置。
  2. 这是某种有意为之的行为吗?如果是,我在哪里可以了解更多相关信息?

所以问题是模型需要跟踪使用原始数据计算的节点,并在预测新数据时使用这些值。这通常发生在 lm() 调用中的 model.frame() 调用中。 bs() 函数 returns 是 "bs" 的 class 并且在创建 model.frame 时,该列被分派到 splines:::makepredictcall.bs 以尝试捕获边界结。 (您可以在 model.frame.default 函数中看到 makepredictcall 调用。)

但是如果我们比较结果

splineModel1 <- lm(y ~ bs(x, y, degree = 3, knots = c(3, 6)))
attr(terms(splineModel1), "predvar")
# list(y, bs(x, degree = 3L, knots = c(3, 6), Boundary.knots =  c(0.275912734214216, 
# 9.14309860439971), intercept = FALSE))

splineModel2 <- lm(y ~ splines::bs(x, y, degree = 3, knots = c(3, 6)))
attr(terms(splineModel2), "predvar")
# list(y, splines::bs(x, y, degree = 3, knots = c(3, 6)))

注意第二个没有捕获 Boundary.knots。这是因为 splines:::makepredictcall.bs 函数实际上查看调用的名称

function (var, call) {
    if (as.character(call)[1L] != "bs") 
        return(call)
    ...
}

当您在公式中使用 splines::bs 时,则 as.character(call)[1L] returns "splines::bs""bs" 不匹配,因此不会发生任何事情。我不清楚为什么会有这张支票。似乎方法调度应该足以假设它是一个 bs 对象。

在我看来,这似乎不是期望的行为,可能应该修复。但是函数 bs() 不应该在没有加载包的情况下真正被调用,因为像 makepredictcall.bs 这样的函数也不会被导入,所以对这些对象的自定义调度会被破坏。

似乎与样条模型'terms'部分的'predvars'属性中的边界节点值有关。

如果我们称它们为 splineModel_1 和 splineModel_2

predict(splineModel_1, newData)
16
predict(splineModel_2, newData)
6.969746

attr(splineModel_2[["terms"]], "predvars") <- attr(splineModel_1[["terms"]], "predvars")

predict(splineModel_1, newData)
16
predict(splineModel_2, newData)
16

attr(splineModel_1[["terms"]], "predvars")
list(y, bs(x, degree = 3L, knots = c(3, 6), Boundary.knots = c(0.323248628992587, 9.84225275926292), intercept = FALSE))

attr(splineModel_2[["terms"]], "predvars")
list(y, splines::bs(x, y, degree = 3, knots = c(3, 6)))

如您所见,区别在于 Boundary.knots。唯一的区别是拦截默认为 FALSE,因此这可能不相关。 Boundary.knots 取自 x 的最小值和最大值。至于它是由一个版本的 bs 而不是另一个版本设置的,我只能假设这是 lm 代码中的一个遗留物,它寻找 'bs' 而不是 'splines::bs' 来设置 Boundary.knots正确。