如何在同一面板中绘制不同 GAM 的平滑组件?

How do you plot smooth components of different GAMs in same panel?

我有两个 GAM,它们具有相同的预测变量但不同的自变量。我想将两个 GAM 组合成一组图,其中每个预测变量的平滑分量(部分残差)在同一面板中(例如用颜色区分)。可重现的例子:

# Required packages
require(mgcv)
require(mgcViz)

# Dataset
data("swiss")

# GAM models
fit1 <- mgcv::gam(Fertility ~ s(Examination) + s(Education), data = swiss)
fit2 <- mgcv::gam(Agriculture ~ s(Examination) + s(Education), data = swiss)

# Converting GAM objects to a gamViz objects
viz_fit1 <- mgcViz::getViz(fit1)
viz_fit2 <- mgcViz::getViz(fit2)

# Make plotGAM objects
trt_fit1 <- plot(viz_fit1, allTerms = T) + l_fitLine()
trt_fit2 <- plot(viz_fit2, allTerms = T) + l_fitLine()

# Print plots
print(trt_fit1, pages = 1)
print(trt_fit2, pages = 1)

fit1 的绘图如下所示:

而 fit2 是这样的:

所以我想将两个考试合并到一个面板中,将两个教育合并到另一个面板中,显示具有不同 color/linetype 的自变量(来自不同的 GAM)。

如果您希望它们在同一个图中,您可以使用 trt_fit1[["plots"]][[1]]$data$fit 从您的拟合中提取数据并自己绘制它们。我从mgcViz github看剧情风格。您可以根据需要添加第二个轴或刻度。

library(tidyverse)
exam_dat <- 
  bind_rows(trt_fit1[["plots"]][[1]]$data$fit %>% mutate(fit = "Fit 1"), 
            trt_fit2[["plots"]][[1]]$data$fit %>% mutate(fit = "Fit 2"))
  

ggplot(data = exam_dat, aes(x = x, y = y, colour = fit)) +
  geom_line() +
  labs(x = "Examination", y = "s(Examination)") + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

为了简单地将它们放在同一个面板上,您可以使用 gridExtra 因为 fit1fit2 有一个 ggplot 对象。

gridExtra::grid.arrange(
  trt_fit1[["plots"]][[2]]$ggObj, 
  trt_fit2[["plots"]][[2]]$ggObj, 
  nrow = 1)

reprex package (v2.0.1)

于 2022-02-18 创建

您也可以使用我的 {gratia} 和 compare_smooths() 功能来做到这一点:

library("gratia")
library("mgcv")

# Dataset 
data("swiss") 

# GAM models 
fit1 <- gam(Fertility ~ s(Examination) + s(Education),
            data = swiss, method = "REML")
fit2 <- gam(Agriculture ~ s(Examination) + s(Education),
            data = swiss, method = "REML")

# create and object that contains the info to compare smooths
comp <- compare_smooths(fit1, fit2)

# plot
draw(comp)

这会产生

compare_smooth() 的输出是嵌套数据框 (tibble)

r$> comp                                                                        
# A tibble: 4 × 5
  model smooth         type  by    data              
  <chr> <chr>          <chr> <chr> <list>            
1 fit1  s(Education)   TPRS  NA    <tibble [100 × 3]>
2 fit2  s(Education)   TPRS  NA    <tibble [100 × 3]>
3 fit1  s(Examination) TPRS  NA    <tibble [100 × 3]>
4 fit2  s(Examination) TPRS  NA    <tibble [100 × 3]>

因此,如果您想对绘图等进行自定义,您需要知道如何使用嵌套数据框或只做

library("tidyr")
unnest(comp, data)

这让你:

r$> unnest(comp, data)                                                          
# A tibble: 400 × 8
   model smooth       type  by      est    se Education Examination
   <chr> <chr>        <chr> <chr> <dbl> <dbl>     <dbl>       <dbl>
 1 fit1  s(Education) TPRS  NA     1.19  3.48      1             NA
 2 fit1  s(Education) TPRS  NA     1.37  3.20      1.53          NA
 3 fit1  s(Education) TPRS  NA     1.56  2.94      2.05          NA
 4 fit1  s(Education) TPRS  NA     1.75  2.70      2.58          NA
 5 fit1  s(Education) TPRS  NA     1.93  2.49      3.10          NA
 6 fit1  s(Education) TPRS  NA     2.11  2.29      3.63          NA
 7 fit1  s(Education) TPRS  NA     2.28  2.11      4.15          NA
 8 fit1  s(Education) TPRS  NA     2.44  1.95      4.68          NA
 9 fit1  s(Education) TPRS  NA     2.59  1.82      5.20          NA
10 fit1  s(Education) TPRS  NA     2.72  1.71      5.73          NA
# … with 390 more rows