提取用于在 mgcv 中绘制平滑图的数据

Question

几年前的This thread 描述了如何提取用于绘制拟合 gam 模型的平滑分量的数据。它有效，但仅当存在一个平滑变量时。我有不止一个平滑变量，不幸的是我只能从系列的最后一个中提取平滑。这是一个例子：

library(mgcv)
a = rnorm(100)
b = runif(100)
y = a*b/(a+b)

mod = gam(y~s(a)+s(b))
summary(mod)

plotData <- list()
trace(mgcv:::plot.gam, at=list(c(25,3,3,3)), 
        #this gets you to the location where plot.gam calls plot.mgcv.smooth (see ?trace)
        #plot.mgcv.smooth is the function that does the actual plotting and
        #we simply assign its main argument into the global workspace
        #so we can work with it later.....
        quote({
                    #browser()
                    plotData <<- c(plotData, pd[[i]])
                }))
plot(mod,pages=1)
plotData

我正在尝试获取a 和b 的估计平滑函数，但列表plotData 只给出了b 的估计值。我研究了plot.gam 函数的内容，但我很难理解它是如何工作的。如果有人已经解决了这个问题，我将不胜感激。

Answer 1

mgcv >= 1.8-6 的更新答案

从 mgcv 的 1.8-6 版开始，plot.gam() 现在以不可见的方式返回绘图数据（来自 ChangeLog）：

• plot.gam 现在静默返回绘图数据列表，以帮助高级 用户 (Fabian Scheipl) 制作自定义情节。

因此，在原始答案中使用下面示例中的mod，可以这样做

> plotdata <- plot(mod, pages = 1)
> str(plotdata)
List of 2
 $ :List of 11
  ..$ x      : num [1:100] -2.45 -2.41 -2.36 -2.31 -2.27 ...
  ..$ scale  : logi TRUE
  ..$ se     : num [1:100] 4.23 3.8 3.4 3.05 2.74 ...
  ..$ raw    : num [1:100] -0.8969 0.1848 1.5878 -1.1304 -0.0803 ...
  ..$ xlab   : chr "a"
  ..$ ylab   : chr "s(a,7.21)"
  ..$ main   : NULL
  ..$ se.mult: num 2
  ..$ xlim   : num [1:2] -2.45 2.09
  ..$ fit    : num [1:100, 1] -0.251 -0.242 -0.234 -0.228 -0.224 ...
  ..$ plot.me: logi TRUE
 $ :List of 11
  ..$ x      : num [1:100] 0.0126 0.0225 0.0324 0.0422 0.0521 ...
  ..$ scale  : logi TRUE
  ..$ se     : num [1:100] 1.25 1.22 1.18 1.15 1.11 ...
  ..$ raw    : num [1:100] 0.859 0.645 0.603 0.972 0.377 ...
  ..$ xlab   : chr "b"
  ..$ ylab   : chr "s(b,1.25)"
  ..$ main   : NULL
  ..$ se.mult: num 2
  ..$ xlim   : num [1:2] 0.0126 0.9906
  ..$ fit    : num [1:100, 1] -0.83 -0.818 -0.806 -0.794 -0.782 ...
  ..$ plot.me: logi TRUE

其中的数据可用于自定义绘图等。

下面的原始答案仍然包含有用的代码，用于生成用于生成这些图的相同类型的数据。

---

原答案

有几种方法可以轻松做到这一点，并且都涉及在协变量范围内从模型进行预测。然而，诀窍是将一个变量保持在某个值（比如它的样本平均值），同时在其范围内改变另一个变量。

这两种方法涉及：

1. 预测数据的拟合响应，包括截距和所有模型项（其他协变量保持固定值），或
2. 根据上述模型进行预测，但返回每个项的贡献

其中第二个更接近（如果不完全是）plot.gam 所做的。

这是一些适用于您的示例并实现上述想法的代码。

library("mgcv")
set.seed(2)
a <- rnorm(100)
b <- runif(100)
y <- a*b/(a+b)
dat <- data.frame(y = y, a = a, b = b)

mod <- gam(y~s(a)+s(b), data = dat)

现在生成预测数据

pdat <- with(dat,
             data.frame(a = c(seq(min(a), max(a), length = 100),
                              rep(mean(a), 100)),
                        b = c(rep(mean(b), 100),
                              seq(min(b), max(b), length = 100))))

预测模型对新数据的拟合响应

这会从上方执行项目符号 1

pred <- predict(mod, pdat, type = "response", se.fit = TRUE)

> lapply(pred, head)
$fit
        1         2         3         4         5         6 
0.5842966 0.5929591 0.6008068 0.6070248 0.6108644 0.6118970 

$se.fit
       1        2        3        4        5        6 
2.158220 1.947661 1.753051 1.579777 1.433241 1.318022

然后您可以针对pdat 中的协变量绘制$fit - 尽管请记住我的预测保持b 不变然后保持a 不变，因此在针对@ 绘制拟合时只需要前100 行987654334@ 或针对b 的第二个 100 行。例如，首先将fitted和upper和lower置信区间数据添加到预测数据的数据框

pdat <- transform(pdat, fitted = pred$fit)
pdat <- transform(pdat, upper = fitted + (1.96 * pred$se.fit),
                        lower = fitted - (1.96 * pred$se.fit))

然后使用行1:100 为变量a 和101:200 为变量b 绘制平滑

layout(matrix(1:2, ncol = 2))
## plot 1
want <- 1:100
ylim <- with(pdat, range(fitted[want], upper[want], lower[want]))
plot(fitted ~ a, data = pdat, subset = want, type = "l", ylim = ylim)
lines(upper ~ a, data = pdat, subset = want, lty = "dashed")
lines(lower ~ a, data = pdat, subset = want, lty = "dashed")
## plot 2
want <- 101:200
ylim <- with(pdat, range(fitted[want], upper[want], lower[want]))
plot(fitted ~ b, data = pdat, subset = want, type = "l", ylim = ylim)
lines(upper ~ b, data = pdat, subset = want, lty = "dashed")
lines(lower ~ b, data = pdat, subset = want, lty = "dashed")
layout(1)

这会产生

如果您想要一个通用的 y 轴刻度，请删除上面的两个 ylim 行，将第一行替换为：

ylim <- with(pdat, range(fitted, upper, lower))

预测各个平滑项对拟合值的贡献

上面2中的想法几乎是用同样的方式完成的，但是我们要求type = "terms"。

pred2 <- predict(mod, pdat, type = "terms", se.fit = TRUE)

这将返回$fit 和$se.fit 的矩阵

> lapply(pred2, head)
$fit
        s(a)       s(b)
1 -0.2509313 -0.1058385
2 -0.2422688 -0.1058385
3 -0.2344211 -0.1058385
4 -0.2282031 -0.1058385
5 -0.2243635 -0.1058385
6 -0.2233309 -0.1058385

$se.fit
      s(a)      s(b)
1 2.115990 0.1880968
2 1.901272 0.1880968
3 1.701945 0.1880968
4 1.523536 0.1880968
5 1.371776 0.1880968
6 1.251803 0.1880968

只需将$fit 矩阵中的相关列与pdat 中的相同协变量绘制出来，同样只使用第一组或第二组 100 行。再举个例子

pdat <- transform(pdat, fitted = c(pred2$fit[1:100, 1], 
                                   pred2$fit[101:200, 2]))
pdat <- transform(pdat, upper = fitted + (1.96 * c(pred2$se.fit[1:100, 1], 
                                                   pred2$se.fit[101:200, 2])),
                        lower = fitted - (1.96 * c(pred2$se.fit[1:100, 1], 
                                                   pred2$se.fit[101:200, 2])))

然后使用行1:100 为变量a 和101:200 为变量b 绘制平滑

layout(matrix(1:2, ncol = 2))
## plot 1
want <- 1:100
ylim <- with(pdat, range(fitted[want], upper[want], lower[want]))
plot(fitted ~ a, data = pdat, subset = want, type = "l", ylim = ylim)
lines(upper ~ a, data = pdat, subset = want, lty = "dashed")
lines(lower ~ a, data = pdat, subset = want, lty = "dashed")
## plot 2
want <- 101:200
ylim <- with(pdat, range(fitted[want], upper[want], lower[want]))
plot(fitted ~ b, data = pdat, subset = want, type = "l", ylim = ylim)
lines(upper ~ b, data = pdat, subset = want, lty = "dashed")
lines(lower ~ b, data = pdat, subset = want, lty = "dashed")
layout(1)

这会产生

请注意此情节与之前制作的情节之间的细微差别。第一个图包括截距项的影响和来自b 平均值的贡献。在第二个图中，仅显示了 a 的平滑器值。

Answer 2

Gavin 给出了一个很好的答案，但我想根据原始引用的帖子提供一个（因为我只是花了很多时间弄清楚它是如何工作的 :)。

我直接使用了https://stats.stackexchange.com/questions/7795/how-to-obtain-the-values-used-in-plot-gam-in-mgcv的代码，还发现我只得到了最后一个返回的模型。原因是跟踪代码 sn-p 放置在 mgcv::plot.gam 函数中的位置。您需要确保将代码放在迭代 m 的 for 循环中，并通过 at 参数进行控制。

以下跟踪非常适合我的 mgcv:::plot.gam 版本

plotData <<- list()
trace(mgcv:::plot.gam, at=list(c(26,3,4,3)), 
quote({
       plotData[[i]] <<- pd[[i]]
  })
)

它在 mgcv:::plot.gam 函数中的这个块之后插入跟踪调用：

if (m > 0) 
    for (i in 1:m) if (pd[[i]]$plot.me && (is.null(select) || 
        i == select)) {

现在 plotData 的元素将对应于绘制的不同变量。我发现两个函数对于找出插入此跟踪调用的正确位置非常有帮助

edit(mgcv:::plot.gam)
as.list(body(mgcv::::plot.gam))

Answer 3

除了 Gavin Simpson 的出色回答之外，现在还有一个名为 itsadug 的 R 包，它提供了几个函数来可视化与 mgcv 匹配的 GAM。

其中有 plot_smooth（根据帮助“绘制总和效果图并可选择移除随机效果”）。如果我正确理解文档，这接近 Gavin Simpson 提到的选项 1。

还有 get_modelterm，它返回一个列表（或可选的 data.frame），其中包含所选平滑项的估计值。这似乎等同于选项 2（或从 plot.gam 返回的值，但没有绘图）。