如何将基函数叠加到 GAM 图上

Question

我正在尝试找到一种方法来复制与此类似的图，其中样条曲线和构成这些样条曲线的基函数都绘制在同一窗口中。

我在下面分别成功地完成了这两项工作：

#### Load Libraries ####
library(mgcv)
library(tidyverse)
library(gratia)
library(gamair)
library(ggpubr)

#### Set Theme ####
theme_set(theme_bw())

#### Add Data ####
data("wesdr")
wes <- as_tibble(wesdr)
wes

#### Fit GAM ####
fit <- gam(
  ret ~ s(dur, bs = "cr"),
  method = "REML",
  family = binomial,
  data = wes
)

#### Plot Basis Functions ####
b <- draw(basis(fit))
s <- draw(fit)
ggarrange(b,s)

但是，我不确定如何将它们混合在一起。简单地把它们放在一起显然是行不通的：

#### Attempt at Plotting BF and Spline ####
wes %>% 
  ggplot(aes(x=dur,
             y=ret))+
  stat_smooth(method = "gam",
              method.args = list(family = binomial),
              formula = y ~ s(x, bs = "cr"),
              se = T,
              color = "steelblue")+
  geom_line(data = basis(fit),
            aes(x=dur,
                y=value,
                color=bf))

如何才能做到这一点？

Answer 1

图中并没有真正显示或使用任何响应数据，只有样条协变量的值，它并不真正需要它，除非你想要漂亮、平滑的基函数。如果你想为估计的样条绘制基础，那就是另一回事了。假设您可能需要两者（第一个用于教授或解释样条曲线的工作原理，第二个用于根据特定拟合进行解释），下面我将展示如何生成这两种图形。

选项 1，使用基础和用户指定的权重

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

df <- data.frame(x = seq(0, 1, length = 100))
bs <- basis(s(x, bs = "bs", k = 10), data = df)

# let's weight the basis functions (simulating model coefs)
set.seed(1)
betas <- data.frame(bf = factor(1:10), beta = rnorm(10))

# we need to merge the weights for each basis function with the basis object
bs <- bs |>
    left_join(betas, by = join_by("bf" == "bf")) |>
    mutate(value_w = value * beta)

# now we want to sum the weighted basis functions for each value of `x`
spl <- bs |>
    group_by(x) |>
    summarise(spline = sum(value_w))

# now plot
bs |> 
    ggplot(aes(x = x, y = value_w, colour = bf, group = bf)) +
    geom_line(show.legend = FALSE) +
    geom_line(aes(x = x, y = spline), data = spl, linewidth = 1.5,
              inherit.aes = FALSE) +
    labs(y = expression(f(x)), x = "x")

这会产生：

选项 2，使用估计模型

如果你想为实际模型拟合做这个，你可以按照上面的例子，但是你需要在样条中包含可识别性约束（见

?basis

）并从向量中提取基函数的正确权重

coef(m)

.

返回的模型系数

{gratia} 的

basis()

有一种拟合模型的方法，可以自动执行此过程。

dat <- data_sim("eg1", seed = 4)
m <- gam(y ~ s(x0) + s(x1) + s(x2, bs = "bs") + s(x3), data = dat, method = "REML")

# data to evaluate the basis at
ds <- data_slice(dat, x2 = evenly(x2, n = 200))

# generate a tidy representation of the fitted basis functions
x2_bs <- basis(m, term = "s(x2)", data = ds)

# compute values of the spline by summing basis functions at each x2
x2_spl <- x2_bs |>
    group_by(x2) |>
    summarise(spline = sum(value))

# now plot
x2_bs |> 
    ggplot(aes(x = x2, y = value, colour = bf, group = bf)) +
    geom_line(show.legend = FALSE) +
    geom_line(aes(x = x2, y = spline), data = x2_spl, linewidth = 1.5,
              inherit.aes = FALSE) +
    labs(y = expression(f(x2)), x = "x2")

这产生

要获得您想要的最终版本（具有可信区间），请使用

smooth_estimates()

在相同的协变量值下评估样条，而不是手动对基函数求和：

# evaluate the spline at the same values as we evaluated the basis functions
x2_sm <- smooth_estimates(m, "s(x2)", data = ds) |>
    add_confint()

# now plot
x2_bs |> 
    ggplot(aes(x = x2, y = value, colour = bf, group = bf)) +
    geom_line(show.legend = FALSE) +
    geom_ribbon(aes(x = x2, ymin = lower_ci, ymax = upper_ci),
                data = x2_sm, # <---- new !
                inherit.aes = FALSE, alpha = 0.2) +
    geom_line(aes(x = x2, y = est), data = x2_sm, # <---- new !
              linewidth = 1.5, inherit.aes = FALSE) +
    labs(y = expression(f(x2)), x = "x2")

产生

你哪里错了？

我认为您的方法不是出于几个原因。

1. draw()

   方法不返回基础数据。通过设计（因为

   ggplot()

   的工作方式），他们返回 ggplot 对象。最好使用函数（就像您使用

   basis()

   所做的那样）来获得您想要的输出，然后使用

   ggplot()

   自己绘制它们，正如我在选项 2 的最后一个示例中展示的那样。

2. 永远不要使用

   geom_smooth()

   或

   stat_smooth()

   来适应GAM。很容易犯错误；在这里你忘了要求

   method = "REML"

   ，你需要在

   method_args = list(method = "REML")

   电话中通过

   stat_smooth()

   完成。

不过你的方法并没有错；请注意，图中左侧的许多基函数都是负的，因此它们将拟合样条拉低，即使其他一些基函数在拟合样条上方达到峰值。

一个最后的评论；使用 {patchwork} 包来安排

draw()

返回的对象，因为你会得到更好的对齐方式。

library("patchwork")
b + s + plot_layout(ncol = 2)

draw.gam()

和 {gratia} 中的许多其他

draw()

方法已经返回 patchworks，而不是简单的 ggplot 对象，因此如果您使用 {patchwork} 的布局工具，您将获得最佳兼容性。