更改 draw（） plot 的正确 R 语法

问：

我需要帮助翻译“对我的数据中空间变化的长期趋势进行建模”，并将下面的图更改为正确的 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)

R 绘制 MGCV

你没有改变，等等，但我不认为这是问题的根源。您的模型会与？尝试对平滑进行重新排序以使算法更容易，因此在 {gratia} 无法直观地了解您想要的内容时，不可能强制执行此操作。这是{gratia}中已知的劣势：github.com/gavinsimpson/gratia/issues/191 由于我希望在几周内为我的课程解决这个问题，我希望能解决这个问题。同时，我将展示如何手动执行此操作（在我的会议之后）dbsti(CYR, Latitude, Longitude, d = c(1,2), bs = c("tp", "ds"), ....)bam()

0赞 Nate 3/3/2023

嗨，@GavinSimpson，是的，它用 ti（CYR， Latitude， Longitude， d = c（1,2）， bs = c（“tp”， “ds”）， ....）和 facets = Longitude 绘制。您是否也有该课程的链接，或者它会在 YouTube 上录制？我也希望看到它。

1赞 Gavin Simpson 3/6/2023

课程在这里：physalia-courses.org/courses-workshops/gams-in-r，但它是付费的，所以不，这次视频不会出现在 YouTube 上。材料将在（以前的版本已经存在）：github.com/gavinsimpson/physalia-gam-course

1赞 Gavin Simpson 3/17/2023

现在使用 GitHub 版本的 {gratia}（版本 0.8.1.24 或更高版本）应该会自动工作，因为我已经实现了此功能并关闭了 #191。

1赞 Nate 3/17/2023

太棒了，谢谢！如果你有时间，我还有一个关于何时在 plot.gam（）或 draw（）中使用“unconditional = FALSE”的问题 - 例如，如何/何时可以提前知道平滑参数而不是从数据中估计？stats.stackexchange.com/questions/609685/......

答：

1赞 Gavin Simpson 3/6/2023 #1

如果你知道一些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 轴上绘制和。x0x1

该代码产生

在包代码中执行此操作比较棘手，因此它会自动发生，如果模型中没有提示（没有 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)

这会产生

1赞 Gavin Simpson 3/17/2023 #2

从版本 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 样条，尽管在公式中定义项时项的排序方式不是这样。x0te()

与我对 OP 问题的评论相反，这个（行为）仅在 0.8.1.25 中添加; .24 仅添加了对及其方法的支持。draw.gam()smooth_estimates()draw()

我将其作为单独的答案发布，因为这种新行为仅适用于 2D 边缘平滑;如果我们有，代码目前没有确定我们可以通过在面板的 x 和 y 轴上生成和的数据来做得更好，然后 facet by 和。在这种情况下，您仍然需要我发布的原始答案中的想法来手动生成所需的输出。te(x0, x1, x2, x3, bs = c("cr", "ds"), d = c(1, 3))x1x2x3x1

与表示您的 dplyr 已过期相关的错误。最好将已安装的软件包设置为要使用的任意数量的 CPU 内核。希望您可以在没有 IT 支持的情况下更新软件包。但请注意，gratia 需要一个相当新的 R 版本（4.1.0 或更高版本）case_matchupdate.packages(ask = FALSE, checkBuilt = TRUE, Ncpus = X)X

1赞 Gavin Simpson 3/18/2023

将因子放入张量积的唯一方法是通过随机效应基，虽然这适用于，但我真的不确定张量积中的随机效应边际会如何处理。如果某些东西不能开箱即用，你最好按照前面的答案使用加上你自己的 ggplot 代码，但传递你想要评估平滑度的数据切片 - 这样你就可以修复你想要的值。smooth_estimates()draw()smooth_estimates()smooth_estimates()year

1赞 Gavin Simpson 3/18/2023

@Nate我在原始答案中添加了一个示例，展示了如何执行我上面提到的操作

上一个：无法找到签名“BAM”的函数“draw”的继承方法

下一个：扩大 ggplot 的边距

更改 draw（） plot 的正确 R 语法 - gratia package

Correct R syntax for changing draw() plot - gratia package

评论

评论