ordisurf vs mgcv:gam模型

Question

我想更新R vegan中ordisurf函数后面的GAM模型(添加随机因子)。为此，我首先构建了应该与ordisurf后面的模型相同的模型，如下所述：https://www.fromthebottomoftheheap.net/2011/06/10/what-is-ordisurf-doing/

然而，我在mgcv中创建的ordisurf和GAM模型给出了完全不同的结果。下面是沙丘数据的示例代码。有谁能解释一下吗？GAM模型应该是什么样子才能执行类似于ordisurf的操作？(在开始改进它之前，重要的是要知道…：)

ord <-metaMDS(dune, k=3) 

surf <- ordisurf(ord, dune.env$A1, choices = c(2,3))

scrs <- data.frame(scores(ord, display = "sites", choices = c(1,2,3)))

dat <- with(dune, cbind(scrs, dune.env$A1))

mod_23 <- gam(dune.env$A1 ~ s(NMDS2, NMDS3, k = 10), data = dat)

plot.gam(mod_12, se=FALSE, cex = 1, pch = 1, col="blue")

Answer 1

对这类问题的最好回答是阅读源代码，看看实际正在做什么。由于我写了那篇博文，我们决定调整ordisurf()的默认设置，以反映在评估GAM方面的良好实践，并更加严格地选择模型。

其中第一个意味着我们使用mgcv::gam()使用method = 'REML'来拟合模型，以利用REML平滑度选择。GCV平滑度选择是(目前) mgcv::gam()中的默认设置，但不建议这样做，mgcv的开发人员Simon Wood暗示，mgcv的某些未来版本可能会将默认设置从method = 'GCV.Cp'更改为method = 'REML'。原因是GCV平滑度选择在许多应用程序中被观察到不平滑。当GCV函数的轮廓在全局最优GCV得分附近变得非常平坦时，就会发生这种情况。数据中的细微差异可能会导致选择非常不同的平滑度参数，其中一些参数会导致严重的欠平滑。REML和ML平滑度选择已被证明受此问题的影响较小；REML得分的轮廓倾向于更加弯曲，并且具有明显的最小值。

第二点，也是我们现在使用select = TRUE的原因，是因为观察到GAM中的平滑度惩罚选择了估计曲面的摆动程度，只影响摆动基函数。基展开包含两个基函数，从惩罚的角度来看，它们是完全平滑的；有两个线性二维基函数代表两个线性平面，根据定义，这两个平面具有零曲率(零二阶导数)，因此根本不会影响摆动惩罚(默认情况下，对估计样条曲线的曲率使用惩罚)。最终结果是，平滑度惩罚可以一直惩罚到估计的平面表面，但不会进一步。这意味着在实践中，用户可以将ordisurf()模型的结果解释为线性表面，即使该线性表面在统计上不显着(许多用户只看图，而不是底层的GAM，这将告诉他们他们的平面是否显着)。select = TRUE所做的是在光滑度上增加第二个惩罚，这只影响完全光滑基函数。这具有GAM能够惩罚/收缩出我们正在估计的表面的效果。换句话说，它可以将表面收缩为平面，水平表面的零影响===排序配置和响应变量之间的零关系。

总而言之，我觉得这些选项最好地防止了用户获得误报。

如果您将gam()调用更改为：

mod_23 <- gam(dune.env$A1 ~ s(NMDS2, NMDS3, k = 10), data = dat,
              method = 'REML', select = TRUE)

然后，您应该得到与ordisurf()生成的输出相同的输出

r$> mod_23                                                                       

Family: gaussian 
Link function: identity 

Formula:
dune.env$A1 ~ s(NMDS2, NMDS3, k = 10)

Estimated degrees of freedom:
0.285  total = 1.28 

REML score: 43.25057     

r$> surf                                                                         

Family: gaussian 
Link function: identity 

Formula:
y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)

Estimated degrees of freedom:
0.285  total = 1.28 

REML score: 43.25057