提取用于以mgcv形成平滑图的数据

Question

几年前的这个主题描述了如何提取用于绘制拟合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)生成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

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

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

原始答案

有两种方法可以轻松地完成这项工作,并且都涉及从模型中预测协变量的范围.然而,技巧是将一个变量保持在某个值(比如它的样本均值),同时在其范围内改变另一个变量.

这两种方法涉及:

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

然后你可以$fit对协变量进行绘图pdat- 虽然记得我有预测保持b常数然后保持a常数,所以你只需要在绘制拟合时的前100行a或者反对的第二行100 b.例如,首先将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显示平滑器的值.