r - 来自 gam predict$fit 和 predict$se.fit 的 ggplot 置信带

Question

我有以下变量：

prod: 正整数

tenure: 正数

cohort： 因素

这是一些具有这些规格的模拟数据。

set.seed(123)
my_data <- data.frame(prod   = rnbinom(10000, mu = 2.5, size = 1.5),
                      tenure = rexp(10000),
                      cohort = factor(sample(2011:2014, size = 10000, replace = TRUE,
                                             prob = c(0.17, 0.49, 0.26, 0.08))))

我使用以下模型拟合mgcv:gam：

library(mgcv)
mod <- gam(prod ~ s(tenure, by = cohort) + cohort, data = my_data, family = nb())

得到预测及其标准误：

preds   <- predict(mod, se.fit = TRUE)
my_data <- data.frame(my_data,
                      mu   = exp(preds$fit),
                      low  = exp(preds$fit - 1.96 * preds$se.fit),
                      high = exp(preds$fit + 1.96 * preds$se.fit))

package:ggplot2使用它来获取每个群组的平滑预测是相当简单的mu（同时也强制平滑器具有正值）：

library(magrittr)
library(ggplot2)
library(splines)
my_plot <-
  ggplot(my_data, aes(x = tenure, y = mu, color = cohort)) %>%
  + geom_smooth(method  = "glm",
                formula = y ~ ns(x, 3),
                family  = "quasipoisson",
                fill    = NA)

但我想平滑 GAM 的置信区间。我该如何添加这些？

不是答案

1. 删除fill = NA. 没有。这些置信区间将无限小，因为在一个队列中按任期进行的预测完全相同。
2. 添加对 的调用geom_ribbon(aes(x = tenure, ymin = low, ymax = high))。没有。这给了我一个超级摇摆不定、不平滑的信心带。
3. 使用package:ggvis！没有package:ggvis答案，请，除非没有办法做到这一点ggplot2。我当前的绘图框架是ggplot2，我现在坚持使用它，除非我必须切换才能执行此绘图。

Answer 1

这对我有用。

require(ggplot2)
require(mgcv)

set.seed(123)
my_data <- data.frame(prod   = rnbinom(10000, mu = 2.5, size = 1.5),
                      tenure = rexp(10000),
                      cohort = factor(sample(2011:2014, size = 10000, replace = TRUE,
                                             prob = c(0.17, 0.49, 0.26, 0.08))))
mod <- gam(prod ~ s(tenure, by = cohort) + cohort, data = my_data, family = nb())
preds   <- predict(mod, se.fit = TRUE)
my_data <- data.frame(my_data,
                      mu   = exp(preds$fit),
                      low  = exp(preds$fit - 1.96 * preds$se.fit),
                      high = exp(preds$fit + 1.96 * preds$se.fit))

ggplot(my_data, aes(x = tenure, y = prod, color = cohort)) +
  geom_point() + 
  geom_smooth(aes(ymin = low, ymax = high, y = mu), stat = "identity")