在同一张图中组合两个 gamm 输出?

Combine two gamm outputs in same graph?

轴和变量相同,但原始数据框不同

mod_a <- gamm(Response ~ s(variable1) + s(variable2) + s(variable3), data=df1)
mod_b <- gamm(Response ~ s(variable1) + s(variable2) + s(variable3), data=df2)

如何将它们组合成一个图并为每个图进行颜色编码?所以它看起来像这样(下图)?以便绘图同时显示 mod_amod_b,即使它们最初来自不同的数据框?

示例数据集:

df1 <- data.frame (Response = c(00, 17, 03, 23, 02, 21, 24, 21, 16, 24, 15, 28, 07, 30, 11, 07, 21, 14, 10, 05, 14, 17, 02, 03, 18, 28, 05, 16, 14, 02, 18, 26, 30, 06, 11, 06, 25, 03, 20, 19, 30, 16, 24, 12, 22, 20, 23, 20, 14, 26),
                   variable1 = c(26, 00, 26, 03, 29, 25, 18, 24, 22, 17, 18, 15, 20, 23, 29, 17, 02, 21, 25, 05, 28, 17, 13, 03, 29, 01, 12, 06, 05, 09, 04, 17, 12, 27, 25, 14, 06, 05, 05, 06, 01, 26, 26, 08, 19, 25, 30, 29, 18, 07),
                   variable2 = c(08, 03, 22, 09, 10, 00, 06, 22, 23, 02, 06, 08, 19, 06, 29, 27, 14, 24, 01, 08, 15, 10, 24, 04, 27, 09, 19, 20, 16, 04, 00, 02, 26, 21, 09, 26, 29, 19, 03, 19, 30, 14, 26, 28, 28, 15, 11, 19, 08, 07),
                   variable3 = c(12, 07, 15, 21, 23, 19, 02, 00, 28, 27, 08, 22, 04, 18, 14, 18, 15, 20, 27, 19, 24, 07, 05, 26, 05, 28, 21, 26, 22, 30, 18, 01, 19, 05, 24, 18, 29, 15, 06, 11, 19, 13, 16, 07, 22, 08, 27, 17, 21, 25),
                   variable4 = c(07, 21, 24, 16, 30, 14, 27, 14, 24, 13, 28, 15, 11, 24, 19, 12, 02, 30, 19, 27, 03, 12, 23, 16, 17, 12, 04, 17, 01, 07, 29, 12, 03, 20, 04, 27, 19, 10, 18, 08, 15, 29, 11, 03, 16, 08, 11, 19, 25, 13),
                   variable5 = c("sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5"))

df2 <- data.frame (Response  = c(24, 29, 16, 03, 01, 04, 08, 03, 17, 09, 27, 11, 28, 02, 11, 15, 26, 12, 05, 03, 06, 06, 11, 24, 19, 25, 07, 14, 29, 02, 04, 27, 15, 06, 18, 10, 30, 16, 17, 22, 07, 24, 02, 24, 17, 09, 00, 20, 06, 27),
                   variable1 = c(22, 11, 19, 08, 03, 16, 04, 20, 12, 25, 08, 21, 04, 07, 09, 28, 25, 04, 27, 17, 00, 22, 29, 08, 17, 06, 12, 16, 08, 00, 16, 24, 20, 09, 10, 10, 04, 24, 11, 00, 07, 21, 15, 11, 05, 00, 07, 05, 25, 03),
                   variable2 = c(11, 21, 01, 06, 18, 22, 10, 19, 26, 16, 12, 08, 18, 11, 25, 16, 16, 25, 02, 29, 22, 02, 01, 03, 10, 08, 16, 19, 07, 10, 05, 17, 04, 24, 20, 29, 23, 00, 01, 18, 10, 24, 15, 09, 14, 26, 30, 30, 04, 29),
                   variable3 = c(15, 06, 24, 29, 04, 07, 26, 14, 21, 15, 18, 02, 27, 09, 09, 24, 09, 15, 23, 15, 09, 13, 08, 07, 14, 03, 03, 07, 27, 21, 06, 30, 03, 03, 27, 11, 01, 05, 03, 14, 10, 20, 30, 10, 22, 23, 03, 30, 30, 25),
                   variable4 = c(03, 22, 10, 07, 23, 08, 12, 06, 25, 17, 12, 28, 21, 28, 18, 21, 15, 17, 23, 10, 11, 21, 12, 10, 26, 04, 18, 18, 26, 25, 20, 02, 15, 28, 17, 04, 14, 28, 01, 13, 16, 05, 14, 02, 06, 15, 16, 26, 29, 07),
                   variable5 = c("sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5","sq1", "sq2", "sq3", "sq4", "sq5"))

library(mgcv)
mod_a <- gamm(Response ~ s(variable1) + s(variable2) + s(variable3), data=df1)
mod_b <- gamm(Response ~ s(variable1) + s(variable2) + s(variable3), data=df2)

plot(mod_a$gam, pages = 1, shade = T, shade.col = 'gray', residuals = T)
plot(mod_b$gam, pages = 1, shade = T, shade.col = 'gray', residuals = T)

一个选项是我的{gratia}套餐:

library('dplyr')
library('gratia')
# can't handle gamm objects just yet so extract the $gam compoents
ma <- mod_a$gam
mb <- mod_b$gam

然后使用 compare_smooths(),它有 gam 个对象的方法

compare_smooths(ma, mb)

这是一个 returns 嵌套的 tibble

r$> compare_smooths(ma, mb)
# A tibble: 6 × 5
  model smooth       type  by    data
  <chr> <chr>        <chr> <chr> <list>
1 ma    s(variable1) TPRS  NA    <tibble [100 × 3]>
2 mb    s(variable1) TPRS  NA    <tibble [100 × 3]>
3 ma    s(variable2) TPRS  NA    <tibble [100 × 3]>
4 mb    s(variable2) TPRS  NA    <tibble [100 × 3]>
5 ma    s(variable3) TPRS  NA    <tibble [100 × 3]>
6 mb    s(variable3) TPRS  NA    <tibble [100 × 3]>

其中有一个 draw() 方法:

compare_smooths(ma, mb) %>%
  draw()

产生

如果你想做一个特定的平滑使用smooths参数

r$> compare_smooths(ma, mb, smooths = "s(variable1)")
# A tibble: 2 × 5
  model smooth       type  by    data
  <chr> <chr>        <chr> <chr> <list>
1 ma    s(variable1) TPRS  NA    <tibble [100 × 3]>
2 mb    s(variable1) TPRS  NA    <tibble [100 × 3]>

我将为 gamm 个对象添加一个方法,这样你以后就可以做到

compare_smooths(mod_a, mod_b)