r - 绘图时在 Y 轴上变换变量

Question

我已经使用 gam 来拟合包括样条项的广义加法模型。它带有附加的情节。我想在 Y 轴上显示优势比 (OR)，而不是当前显示的图表。我感谢您的帮助。

标签： rplotgam

Answer 1

为此，您可以只返回链接刻度上的样条值（没有截距），然后对这些值取幂以获得赔率刻度

如果您正在使用，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 的爆发，您将需要大量调整轴限制。