为什么在“gam(y ~ mgcv::s...)”中使用“mgcv::s”会导致错误？答案

【问题标题】：Why does using "mgcv::s" in "gam(y ~ mgcv::s...)" result in an error?为什么在“gam(y ~ mgcv::s...)”中使用“mgcv::s”会导致错误？
【发布时间】：2018-06-26 21:41:51
【问题描述】：

我想清楚并在行中使用:: 符号来拟合mgcv::gam。在mgcv::s 的模型调用中使用符号时，我偶然发现了一件事。带有可重现示例/错误的代码如下所示。

原因可能是因为我在模型公式中使用了这个符号，但我不知道为什么这不起作用/不允许。这可能是一些非常具体的语法（我猜可能不是 mgcv 特定的），但也许有人可以帮助我理解这一点以及我对 R 的理解。提前谢谢你。

library(mgcv)
dat <- data.frame(x = 1:10, y = 101:110)
# this results in an error: invalid type (list)...
mgcv::gam(y ~ mgcv::s(x, bs = "cs", k = -1), data = dat)
# after removing the mgcv:: in front of s everything works fine
mgcv::gam(y ~ s(x, bs = "cs", k = -1), data = dat)

# outside of the model call, both calls return the desired function
class(s)
# [1] "function"
class(mgcv::s)
# [1] "function"

【问题讨论】：

这是一个很好的问题，但它有点像“医生，我这样做的时候很痛”。 “好吧，那就别那么做了。”我很欣赏对清晰的渴望，但如果你想在gam() 公式中使用s() 来代替平滑术语，那将是非常疯狂的......
你是对的。我刚刚开始学习这种建模，不确定在其他建模包中是否有任何名为s 的函数，因此，我想弄清楚。我猜想，一旦调用了来自 mcgv 的gam，该调用可能会解析模型调用中的所有内容，这些内容被解释为来自mcgv 的语法。所以mgcv::前面的s是不必要的。尽管如此，使用这种表示法不应导致与一般表示法一致的错误。但既然我知道我不应该在这种特殊情况下使用这个符号，我会听从医生的建议;-)。
将公式视为模型的符号表示，而不是 R 代码。虽然您是对的，在搜索路径上潜伏在 mgcv:s 之前的任何 s() 都会干扰，但操作公式是您可以尝试使用 R 完成的非常困难的事情之一。
@GavinSimpson 正如您在下面所说的，我指出的问题非常具体，与处理这个问题相比，可以花费更多时间。不过，我在这里学到了一些关于 R 的东西，所以感谢您的解释和在这个小问题上花费的时间。

标签： r syntax mgcv

【解决方案1】：

说明

library(mgcv)
#Loading required package: nlme
#This is mgcv 1.8-24. For overview type 'help("mgcv-package")'.

f1 <- ~ s(x, bs = 'cr', k = -1)
f2 <- ~ mgcv::s(x, bs = 'cr', k = -1)

OK <- mgcv:::interpret.gam0(f1)$smooth.spec
FAIL <- mgcv:::interpret.gam0(f2)$smooth.spec

str(OK)
# $ :List of 10
#  ..$ term   : chr "x"
#  ..$ bs.dim : num -1
#  ..$ fixed  : logi FALSE
#  ..$ dim    : int 1
#  ..$ p.order: logi NA
#  ..$ by     : chr "NA"
#  ..$ label  : chr "s(x)"
#  ..$ xt     : NULL
#  ..$ id     : NULL
#  ..$ sp     : NULL
#  ..- attr(*, "class")= chr "cr.smooth.spec"

str(FAIL)
# list()

interpret.gam0的源代码第4行揭示了问题：

head(mgcv:::interpret.gam0)

1 function (gf, textra = NULL, extra.special = NULL)              
2 {                                                               
3     p.env <- environment(gf)                                    
4     tf <- terms.formula(gf, specials = c("s", "te", "ti", "t2", 
5         extra.special))                                         
6     terms <- attr(tf, "term.labels")

由于"mgcv::s" 不匹配，所以您遇到了问题。但是mgcv 确实允许您通过参数"mgcv::s" 传递extra.special 来解决这个问题：

FIX <- mgcv:::interpret.gam0(f, extra.special = "mgcv::s")$smooth.spec
all.equal(FIX, OK)
# [1] TRUE

只是这在高级程序中不是用户可控的：

head(mgcv::gam, n = 10)

#1  function (formula, family = gaussian(), data = list(), weights = NULL, 
#2      subset = NULL, na.action, offset = NULL, method = "GCV.Cp",        
#3      optimizer = c("outer", "newton"), control = list(), scale = 0,     
#4      select = FALSE, knots = NULL, sp = NULL, min.sp = NULL, H = NULL,  
#5      gamma = 1, fit = TRUE, paraPen = NULL, G = NULL, in.out = NULL,    
#6      drop.unused.levels = TRUE, drop.intercept = NULL, ...)             
#7  {                                                                      
#8      control <- do.call("gam.control", control)                         
#9      if (is.null(G)) {                                                  
#10         gp <- interpret.gam(formula)  ## <- default to extra.special = NULL

我同意本·博克的观点。挖掘内部发生的事情是一个很好的练习，但将其视为错误并修复它是一种过度反应。

更多洞察：

mgcv 中的s、te 等与stats::poly 和splines::bs 的逻辑不同。

例如，当您执行X <- splines::bs(x, df = 10, degree = 3) 时，它会评估 x 并直接创建设计矩阵X。
当你做s(x, bs = 'cr', k = 10)时，不做评价；它被解析。

mgcv 的顺利构建需要几个阶段：

由mgcv::interpret.gam 解析/解释，生成更平滑的配置文件；
由mgcv::smooth.construct 进行初始构建，建立基础/设计矩阵和惩罚矩阵（主要在 C 级完成）；
mgcv::smoothCon 的二次构造，它选取“by”变量（例如，为因子“by”复制平滑）、线性函数项、零空间惩罚（如果您使用 select = TRUE）、惩罚重新缩放、居中约束等；
mgcv:::gam.setup 的最终集成，将所有平滑器组合在一起，返回模型矩阵等。

所以，这是一个复杂得多的过程。

【讨论】：

感谢您对 mgcv 中处理步骤的详细解释和见解。我同意这不是必须解决的问题。仍然与 :: 的典型用法不完全一致

【解决方案2】：

这看起来像是mgcv 问题。例如，lm() 函数接受 poly() 或 stats::poly() 并给出相同的结果（除了事物的名称）：

> x <- 1:100
> y <- rnorm(100)
> lm(y ~ poly(x, 3))

Call:
lm(formula = y ~ poly(x, 3))

Coefficients:
(Intercept)  poly(x, 3)1  poly(x, 3)2  poly(x, 3)3  
    0.07074      0.13631     -1.52845     -0.93285  

> lm(y ~ stats::poly(x, 3))

Call:
lm(formula = y ~ stats::poly(x, 3))

Coefficients:
       (Intercept)  stats::poly(x, 3)1  stats::poly(x, 3)2  stats::poly(x, 3)3  
           0.07074             0.13631            -1.52845            -0.93285

它也适用于 splines::bs 函数，因此这不是特定于 poly() 的。

您应该联系mgcv 维护者并指出该包中的这个错误。我猜它是专门寻找s，而不是寻找像mgcv::s 这样的表达式，其计算结果相同。

【讨论】：

感谢您的交叉检查以及您与lm 的示例。至少，我似乎不必更新我对使用:: 的理解。我会按照建议报告这个问题。
我认为这行得通（你的两个例子）因为poly() 和bs() 只是将基扩展作为基函数矩阵返回。所以公式解析可以只评估模型框架中公式的 RHS，然后你会得到一些合理的东西。 s() 所做的不仅仅是返回在 x 处评估的基函数矩阵，因此在处理 gam 公式时，比您提到的 lm 案例中发生的事情要多得多。跨度>
什么是错误，或者至少是一个错误，在你的全局环境中定义的任何s() 都会干扰公式的解析。理想情况下，gam() 在评估平滑时只会使用它的s()。然而，我怀疑在这个解析东西的公式中存在龙，与其试图让这个对这些问题变得健壮，不如让它保持现状更好？
嗯，它是 - 这些规则可能不是我们所期望的（即我们希望它使用它自己的 s() 而不是其他东西）。这与解析公式不同——正如我所提到的，如果你看的话，lm 通过评估模型框架环境中的 rhs 来工作。我认为mgcv::s() 不可能做到这一点，因为它不会返回任何符合 R 的评估公式愿景的东西。 mgcv 随后创建正确的模型矩阵，但 it 必须这样做；与mgcv::s 的橙子相比，你的例子是苹果。
我还没有深入研究代码，但似乎发生了这样的事情：gam 查看名称 s 的公式，并通过替换其他内容来特别处理它，然后在结果上调用函数s()。如果它找到一个名为 s() 的不同函数（如您所提到的），或者如果该函数以不同的方式指定（如提到的原始问题），则会失败。如果它做两件事不同会更好：1，如果公式中有函数调用，则使用标准的R语义来查找函数，以及2，查看函数是否需要特殊处理。