Извлечение данных, используемых для создания гладкого участка в mgcv

Этот поток из нескольких лет назад описывает, как извлекать данные, используемые для построения гладких компонентов установленной модели гаммы. Он работает, но только тогда, когда есть одна плавная переменная. У меня есть более чем одна плавная переменная, и, к сожалению, я могу только извлечь сглаживания из последней серии. Вот пример:

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, и мне сложно понять, как это работает. Если кто-то уже решил эту проблему, я был бы благодарен.

Ответ 1

Обновлен ответ для mgcv >= 1,8-6

Начиная с версии 1.8-6 из mgcv, 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

Данные в нем могут использоваться для пользовательских графиков и т.д.

В исходном ответе ниже содержится полезный код для генерации тех же данных, которые используются для создания этих графиков.

Оригинальный ответ

Существует несколько способов сделать это легко, и оба предполагают предсказание модели по диапазону ковариаций. Однако трюк состоит в том, чтобы удерживать одну переменную при некотором значении (например, ее среднее значение), изменяя другую по ее диапазону.

Два метода включают:

Предсказание встроенных ответов для данных, включая перехват и все модельные термины (с другими ковариатами, хранящимися при фиксированных значениях), или
Предскажите модель, как указано выше, но верните вклады каждого термина

Вторая из них ближе к (если не совсем) 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.

Ответ 2

Гэвин дал отличный ответ, но я хотел представить его с точки зрения исходной ссылки (поскольку я просто потратил много времени на выяснение, как это работает:).

Я использовал код непосредственно из https://stats.stackexchange.com/info/7795/how-to-obtain-the-values-used-in-plot-gam-in-mgcv, а также обнаружил, что я вернул только последнюю модель. Причина этого заключается в том, что фрагмент кода трассировки помещается в функцию mgcv:: plot.gam. Вы должны убедиться, что код помещен внутри цикла for, который выполняет итерацию над m, и вы управляете этим аргументом 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))