r - 提取用于在 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 的答案

从 的 1.8-6 版开始mgcv , plot.gam()现在不可见地返回绘图数据(来自 ChangeLog):

• plot.gam now silently returns a list of plotting data, to help advanced users (Fabian Scheipl) to produce custimized plot.

因此，并使用 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

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

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

原答案

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

这两种方法涉及:

预测数据的拟合响应，包括截距和所有模型项(其他协变量保持固定值)，或

从上面的模型预测，但返回每一项的贡献

其中的第二个更接近(如果不完全是) 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

然后您可以绘制 $fit反对 pdat 中的协变量- 虽然记得我有预测 b恒则持 a常数，所以在绘制拟合时只需要前 100 行 a或针对 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 的平滑器的值显示。