你可以使用我的{gratia}和
compare_smooths()
函数来完成这个操作:
library("gratia")
library("mgcv")
data("swiss")
fit1 <- gam(Fertility ~ s(Examination) + s(Education),
data = swiss, method = "REML")
fit2 <- gam(Agriculture ~ s(Examination) + s(Education),
data = swiss, method = "REML")
comp <- compare_smooths(fit1, fit2)
draw(comp)
这会产生
![enter image description here](https://istack.dev59.com/Gj0iY.webp)
compare_smooth()
的输出是一个嵌套数据框(tibble)。
r$> comp
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)
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
为了创建自己的图表,我们需要从非嵌套的数据框开始,并添加置信区间。
ucomp <- unnest(comp, data) %>%
add_confint()
依次绘制每个面板
library("ggplot2")
library("dplyr")
p_edu <- ucomp |>
filter(smooth == "s(Education)") |>
ggplot(aes(x = Education, y = est)) +
geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci, fill = model),
alpha = 0.2) +
geom_line(aes(colour = model)) +
scale_fill_brewer(palette = "Set1") +
scale_colour_brewer(palette = "Set1") +
geom_rug(data = swiss,
mapping = aes(x = Education, y = NULL),
sides = "b", alpha = 0.4) +
labs(title = "s(Education)", y = "Estimate",
colour = "Model", fill = "Model")
p_exam <- ucomp |>
filter(smooth == "s(Examination)") |>
ggplot(aes(x = Examination, y = est)) +
geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci, fill = model),
alpha = 0.2) +
geom_line(aes(colour = model)) +
scale_fill_brewer(palette = "Set1") +
scale_colour_brewer(palette = "Set1") +
geom_rug(data = swiss,
mapping = aes(x = Examination, y = NULL),
sides = "b", alpha = 0.4) +
labs(title = "s(Examination)", y = "Estimate",
colour = "Model", fill = "Model")
现在使用{patchwork}包将图形放在一起。
library("patchwork")
p_edu + p_exam + plot_layout(guides = "collect")
这会使用{ggplot2},因此如果您想对颜色有更多控制,需要查看其他比例尺,例如?scale_fill_manual
,或者提供其他现成的离散比例尺。
我可以在{gratia}中使其中一些操作更加简单 - 我可以允许用户提供用于颜色和填充的比例尺,如果他们提供原始数据,我还可以画出地毯图。
![enter image description here](https://istack.dev59.com/NPlJ1.webp)