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

Question

我想弄清楚，并在行中使用::符号来拟合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"

Answer 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确实允许您通过通过参数extra.special传递"mgcv::s"来解决这个问题。

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

我同意本·博尔克的观点。这是一个很好的练习，以挖掘内部发生的事情，但这是一个过激反应，认为这是一个错误，并修复它。

更有洞察力：

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

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

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

1. mgcv::interpret.gam的解析/解释，它生成一个平滑的概要文件；
2. mgcv::smooth.construct的初步构建，建立基础/设计矩阵和惩罚矩阵(主要是C级的)；
3. mgcv::smoothCon的二次构造，它通过“变量(复制平滑因子”)、线性函数项、空空间惩罚(如果使用select = TRUE)、惩罚重标度、定心约束等进行二次构造；
4. 最后由mgcv:::gam.setup进行集成，它将所有平滑器组合在一起，返回一个模型矩阵等。

因此，这是一个复杂得多的过程。

Answer 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这样的计算结果相同的表达式。