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

Posted 2023-02-16

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

【英文标题】：Bizarre behaviour of lm() and predict.lm() depending on use of explicit namespace accessor

【发布时间】：2017-09-16 05:36:18

【问题描述】：

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

lm 和 predict.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 包的命名空间导出时，就会发生奇怪的事情。以下代码 sn-p 除了该更改之外是相同的：

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))

如果splines 包从未使用library 附加在第二个sn-p 中，则会产生完全相同的结果。我想不出另一种情况，在已加载的包上使用 :: 运算符会改变程序行为。

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

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

为什么会发生这种情况？我一直在尝试跟进predict.lm，但无法确定分歧发生在哪里。 这是某种预期的行为，如果是，我可以从哪里了解更多信息？

【问题讨论】：

另一个奇怪的事情是，如果你查看每个模型的系数，它们是相同的，但预测是不同的。顺便说一句，您不应该两次创建数据，因为每次都会不同（除非您每次都设置相同的种子）。在这里没有什么区别，因为无论哪种方式，数据都是完全确定的，导致相同的模型输出，但最好设置种子并只创建一次数据。

你说得对，设置种子或重用数据会更好。但我想强调的是，第二个 sn-p 是自包含的并且自相矛盾的，与第一个无关——第二个图中的预测值不可能与原始数据拟合的值相差如此之远/

是的，两个模型对象中的系数以及所有数值内容都是相同的。问题出现在预测步骤中的某个地方，该步骤使用拟合模型对象的“调用”和“术语”元素的组合来自动将新的 x 值扩展为 b 样条向量。

【参考方案1】：

所以问题在于模型需要跟踪使用原始数据计算的节点，并在预测新数据时使用这些值。这通常发生在model.frame() 调用中的lm() 调用中。 bs() 函数返回一个 "bs" 类，并在制作 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] 返回的 "splines::bs" 与 "bs" 不匹配，因此不会发生任何事情。我不清楚为什么要进行此检查。似乎方法分派应该足以假设它是一个bs 对象。

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

【讨论】：

啊，太好了。那里的“调度”是通过字符串比较完成的，这似乎很奇怪。考虑到我新预测的点不仅在边界结内，而且在内部结内，我在 x = 4 处的特定预测会受到影响似乎仍然很奇怪。边界结的位置不应影响这些点的估计。

干得好。也就是说，我认为在断定它是不可取的并且应该修复之前，了解代码作者为何包含该初始检查非常重要。

@mb7744 这并不奇怪——当您深入研究 R 的模型拟合机制时，您会发现很多与 formula 和 terms 相关的操作都使用这种类型的操作字符向量经过解析，最终从公式中推导出来。

@JoshO'Brien 对。我可以看到用户将bs 对象直接传递给函数，因此我可以看到它存在的原因。我猜“问题”实际上是它不检查“splines::bs”或“splines:::bs”

当然可以，但是为什么不直接检查对象的类呢？或者是否有必要进行字符串比较，因为在那个阶段可用的只是调用的字符串文字？ （似乎更好的系统是将调用中对象的类/属性与​​字符串call分开记录。例如，可以分配B <- bs(..)，然后调用fit <- lm(y ~ B)并能够使用@ 987654345@... 但在游戏中为时已晚。）

【参考方案2】：

似乎与 splineModel 的 '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。

【讨论】：

很好地注意到了不同的边界结，但关于“lm 代码中寻找'bs'而不是'splines::bs'以正确设置Boundary.knots 的遗物” ，请注意，在基本 R 中没有其他 bs 函数。在没有附加 splines 库的情况下调用我的第一个片段会导致错误。

@mb7744 确实如此。我推测过去可能有一个名为 bs 的函数（在基本 R 或其他常见包中），这使得当时拥有 splines::bs 是个好主意。