Ответ 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.