更改 draw() 图的正确 R 语法 - 感谢包

Question

我需要帮助翻译“建模数据中空间变化的长期趋势”并将下面的图更改为正确的 R 语法。

我想在 x 轴上有经度，在 y 轴上有纬度，每个日历年都有一个面板（“CYR”）。更改平滑函数中变量的顺序会给我这个错误：

     ...ti(Latitude, CYR, Longitude, d = c(2,1), bs = c('ds','tp'), k = c(25, 14)), ..
     ...ti(CYR, Latitude, Longitude, d = c(2,1), bs = c('ds','tp'), k = c(25, 14)), ..

Error in check.term(termi, rec) : 
  bam can not discretize with this nesting structure

它也以某种方式连接到我的其他协变量（删除它们会停止错误，但我无法正确“绘制”绘图）。

library(mgcv)
library(gratia)

m <- bam(occur ~ s(temp) + 
           sal + 
           s(DO) + 
           sed_depth + 
           water_depth + 

           s(fCYR, bs = "re") + # Long-term trend
           
           # Spatial variation
           s(Longitude, Latitude, k = 100, bs = 'tp') + # Try bs = tp (default)
           
           s(fSite, bs = "re") + # Repeated measures design
         
         # long-term trend varies spatially (code I need to change, I think)
         ti(Latitude, CYR, Longitude, d = c(2,1), bs = c('ds','tp'), k = c(25, 14)),
         
         data = toad, 
         method = 'fREML', 
         # knots = knots,
         nthreads = 4, 
         discrete = TRUE,
         family = binomial(link = "logit"), 
         # select = TRUE,
         gamma = 1.5)

draw(m, select = 6)

Answer 1

从版本 0.8.1.25（撰写本文时软件包的开发版本）开始，{gratia} 现在应该能够自动处理此问题。

library("gratia")
library("mgcv")
library("ggplot2")
library("dplyr")

df <- data_sim("eg1", n = 1000,  dist = "normal", scale = 2, seed = 1)
m <- gam(y ~ te(x0, x1, x2, k = c(25, 10), d = c(1,2), bs = c("cr", "ds")),
         data = df, method = "REML")
draw(m, n = 25)

生产

可以看出，它在面板中绘制了 2D Duchon 样条，在小平面上绘制了

x0

，尽管在公式中定义

te()

项时，这并不是这些项的排序方式。

与我对 OP 问题的评论相反，此（

draw.gam()

行为）仅在 0.8.1.25 中添加； .24 添加了对

smooth_estimates()

及其

draw()

方法的支持。

我将其作为单独的答案发布，因为这种新行为仅适用于 2D 边缘平滑；如果我们有

te(x0, x1, x2, x3, bs = c("cr", "ds"), d = c(1, 3))

，代码目前无法确定我们可以通过为

x1

和

x2

生成数据在面板的 x 和 y 轴上进行，然后通过

x3

和

x1

。在那种情况下，您仍然需要我发布的原始答案中的想法来手动生成所需的输出。

Answer 2

如果您了解一些 ggplot，则手工执行此操作相对容易：

library("gratia")
library("mgcv")
library("ggplot2")
library("dplyr")

df <- data_sim("eg1", n = 1000,  dist = "normal", scale = 2, seed = 1)
m <- gam(y ~ te(x2, x0, x1, k = c(5, 20), d = c(1,2), bs = c("ds", "cr")),
         data = df, method = "REML")

sms <- smooth_estimates(m, n = 50, n_3d = 10)

est_lim <- c(-1, 1) * max(abs(sms[["est"]]), na.rm = TRUE)

sms |>
    mutate(fx2 = factor(x2)) |>
    ggplot(aes(x = x0, y = x1, fill = est, group = fx2)) +
    geom_raster(aes(x = x0, y = x1, fill = est, group = fx2)) +
    geom_contour(aes(z = est, group = fx2, fill = NULL), colour = "black") +
    facet_wrap(~ fx2) +
    scale_fill_distiller(palette = "RdBu", type = "div") +
    expand_limits(fill = est_lim)

这里我手动识别示例中的第二个平滑是 2d 样条，所以我知道在绘图的 x 和 y 轴上绘制

x0

和

x1

。

代码生成

在包代码中执行此操作比较棘手，因此它会自动发生，并且如果模型中没有提示（没有 2d 平滑）并且您希望能够指定绘制的内容，那么它会编写一个接口仍然。我想我很快就会在 {gratia} 中提供前者（确定一个 2d 平滑并将其绘制在 x 和 y 轴上，在其他变量上分面），但更难的部分需要更多思考并且可能不会可以直接从

draw(m)

之类的命令执行。

为了完全控制（例如，为面的

x2

指定特定[漂亮]值），生成要评估平滑值的数据切片：

ds <- data_slice(m, x0 = evenly(x0, n = 50), x1 = evenly(x1, n = 50),
                 x2 = seq(0, 1, by = 0.25))
sms <- smooth_estimates(m, data = ds)
est_lim <- c(-1, 1) * max(abs(sms[["est"]]), na.rm = TRUE)

sms |>
    ggplot(aes(x = x0, y = x1, fill = est, group = x2)) +
    geom_raster(aes(x = x0, y = x1, fill = est, group = x2)) +
    geom_contour(aes(z = est, group = x2, fill = NULL), colour = "black") +
    facet_wrap(~ x2) +
    scale_fill_distiller(palette = "RdBu", type = "div") +
    expand_limits(fill = est_lim)

产生