Для этого вы можете просто вернуть значения сплайна на шкале ссылок (без пересечения), а затем возвести в степень значения, чтобы получить вещи по шкале шансов
Если вы используете mgcv::gam(), то выможет сделать это следующим образом:
library('mgcv')
set.seed(1)
dat <- gamSim(1, dist = "binary")
m1 <- gam(y ~ s(x2), data = dat, method = "REML", family = binomial())
pdat <- with(dat, data.frame(x2 = seq(min(x2), max(x2), length = 500)))
pred <- predict(m1, newdata = pdat, se.fit = TRUE, type = "iterms")
pred <- data.frame(x2 = pdat$x2, fit = pred$fit[,1], se.fit = pred$se.fit[,1])
## compute CI on the logit (log-odds) scale
pred <- transform(pred,
upper = fit + (2 * se.fit),
lower = fit - (2 * se.fit))
## transform fitted values + CI to odds scale
pred <- transform(pred,
odds = exp(fit),
oupper = exp(upper),
olower = exp(lower))
## plot
library("ggplot2")
library("cowplot")
theme_set(theme_bw())
## plot on the logit-scale
p1 <- ggplot(pred, aes(x = x2, y = fit)) +
geom_ribbon(aes(x= x2, ymin = lower, ymax = upper),
inherit.aes = FALSE, alpha = 0.1) +
geom_line()
## plot on the odds scale
p2 <- ggplot(pred, aes(x = x2, y = odds)) +
geom_ribbon(aes(x= x2, ymin = olower, ymax = oupper),
inherit.aes = FALSE, alpha = 0.1) +
geom_line()
plot_grid(p1, p2, ncol = 1)
, который производит это:
Верхняя панель является просто представлением ggplotсюжет вы показали.Нижняя панель преобразуется в шкалу шансов.
Вам нужно будет немного ее изменить, если в модели несколько сглаживаний.Линия
pred <- data.frame(....)
должна будет выбрать другие столбцы из компонентов $fit и $se.fit.
Один быстрый способ сделать это, если вы не хотите делать все самостоятельноэто захватить вывод из plot(model)
m2 <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat,
method = "REML", family = binomial())
plt_data <- plot(m2, pages = 1, seWithMean = TRUE)
Теперь plt_data - это список с компонентом для каждого сглаживания.Чтобы воссоздать график, полученный при выполнении plot(m2), нам нужно использовать:
x - это координаты x для сглаживания. fit и* Компоненты se содержат данные о координатах Y (установленные значения) и их стандартные ошибки
Мы напишем функцию для добавления доверительного интервала и, возможно, применим преобразование:
add_ci <- function(df, trans = function(eta) { eta }) {
df <- transform(df, yhat = trans(fit),
upper = trans(fit + (2 * se)),
lower = trans(fit - (2 * se)))
df
}
И примените его к каждому из объектов данных в списке plt_data:
p1dat <- add_ci(as.data.frame(plt_data[[1]][c('x','se','fit')]))
p2dat <- add_ci(as.data.frame(plt_data[[2]][c('x','se','fit')]))
p3dat <- add_ci(as.data.frame(plt_data[[3]][c('x','se','fit')]))
p4dat <- add_ci(as.data.frame(plt_data[[4]][c('x','se','fit')]))
Теперь мы можем построить
p1 <- ggplot(data = p1dat,
aes(x = x, y = yhat)) +
geom_ribbon(aes(x = x, ymin = lower, ymax = upper),
inherit.aes = FALSE, alpha = 0.2) +
geom_line() + labs(y = 's(x0)', x = 'x0')
p2 <- p1 %+% p2dat + labs(y = 's(x1)', x = 'x1')
p3 <- p1 %+% p3dat + labs(y = 's(x2)', x = 'x2')
p4 <- p1 %+% p4dat + labs(y = 's(x3)', x = 'x3')
plot_grid(p1, p2, p3, p4, ncol = 2)
, давая
Далее мы можем применить преобразование
p1dat <- add_ci(as.data.frame(plt_data[[1]][c('x','se','fit')]), trans = exp)
p2dat <- add_ci(as.data.frame(plt_data[[2]][c('x','se','fit')]), trans = exp)
p3dat <- add_ci(as.data.frame(plt_data[[3]][c('x','se','fit')]), trans = exp)
p4dat <- add_ci(as.data.frame(plt_data[[4]][c('x','se','fit')]), trans = exp)
pt1 <- p1 %+% p1dat + labs(y = 's(x0)', x = 'x0') + coord_cartesian(ylim = c(0, 100))
pt2 <- p1 %+% p2dat + labs(y = 's(x1)', x = 'x1') + coord_cartesian(ylim = c(0, 4000))
pt3 <- p1 %+% p3dat + labs(y = 's(x2)', x = 'x2') + coord_cartesian(ylim = c(0, 250))
pt4 <- p1 %+% p4dat + labs(y = 's(x3)', x = 'x3') + coord_cartesian(ylim = c(0, 5))
plot_grid(pt1, pt2, pt3, pt4, ncol = 2)
, которое производит
Как видите, вам нужно много играть с ограничениями осей, когда CI взрывается.