在R中合并mgcv包GAM模型s(x)+s(y)+ti(x,y)的平滑效应可视化

问题背景 #

我在R中常用mgcv包的te()张量平滑函数构建二维GAM模型，代码如下：

library(mgcv)
library(mgcViz)

m = gam(z~te(x,y),data=df)

gv = getViz(m)
plot(gv)

这段代码会生成x和y对z的联合效应二维图。

但我发现采用单变量平滑项加张量交互项的模型s(x) + s(y) + ti(x,y)对数据拟合效果更好：

library(mgcv)
library(mgcViz)

m = gam(z~s(x) + s(y) + ti(x,y),data=df)

gv = getViz(m)
plot(gv)

该模型会分别生成x、y的边际效应图以及二者的交互效应图，虽便于模型诊断，但我需要将这些效应合并成类似te()生成的单张图用于发表级绘图。现请教两个问题：

• 是否有简便的实现方法？
• 是否有充分理由不合并这些效应进行可视化？

解决方案 #

基于Gavin Simpson的建议，我使用mgcv::vis.gam开发了适配复杂场景的解决方案。我的实际模型涉及地理经纬度坐标，且按季度分层，还需在效应图上叠加世界地图。

简化后的模型代码如下：

m = gam(z~s(lon,by=quarter) + s(lat,by=quarter) + 
      ti(lon,lat,by=quarter),data=df)

我编写了以下函数来绘制指定季度的经纬度联合效应：

my_plot_gam_geo_quarter = 
  function(model,quarter,lon.lat=c("lon","lat"),contour.col="black",
           too.far=0.1,asp=1,
           world.col="lightgray",world.border="darkgray",
           cond=list(quarter=quarter)) {
    # 从tmap包获取世界地图数据
    suppressWarnings(data("World",package="tmap",overwrite=FALSE))
    
    mgcv::vis.gam(model,view=lon.lat,cond=cond,plot.type="contour",
                  contour.col=contour.col,too.far=too.far,asp=asp)
    plot(World[1],add=TRUE,col=world.col,border=world.border) 
  }

调用函数绘制第一季度的经纬度效应图：

my_plot_gam_geo_quarter(m,1)

内容的提问来源于stack exchange，提问作者David M. Kaplan