ggplot2 geom_ribbon 来自 mgcv::gamm

Question

我正在尝试根据 gamm 模型的预测添加功能区，这似乎比预期的要难一些，因为 gamm 与 gam 有所不同。

我第一次尝试直接使用 geom_stat，但这不起作用（并且不会使用我的整个模型，其中还包括其他几个协变量）

library(tidyverse); library(mgcv)

dt = cbind(V1=scale(sample(1000)), 
    Age=rnorm(n = 1000, mean = 40, sd = 10), 
    ID=rep(seq(1:500),each=2) %>% as.data.frame()

# Works fine ----
dt %>% ggplot(aes(x=Age, y=V1)) + 
   stat_smooth(method="gam", formula= y~s(x,bs="cr")) 

# Fails horribly :P 
dt %>% ggplot(aes(x=Age, y=V1)) + 
    stat_smooth(method="gamm", formula= y~s(x,bs="cr"))

Maximum number of PQL iterations:  20   
iteration 1  
Warning message:  
Computation failed in `stat_smooth()`:  
no applicable method for 'predict' applied to an object of class "c('gamm', 'list')"   

我已经尝试在 model$gamm 上使用 predict 函数，但我不确定如何使用它，以及如何制作 CI 功能区

dt.model = gamm(V1 ~ s(Age, bs="cr") + s(ID, bs = 're'), data=dt, family="gaussian", discrete=T)

dt$pred = predict(dt.model$gam)

dt %>% ggplot(aes(x = Age, y = V1)) +
   geom_line(aes(group=ID), alpha=.3) +
   geom_point(alpha=.2) +
   geom_smooth(aes(y=pred))

我认识到这是一个糟糕的示例数据，因为这是一个愚蠢的形状。 但我希望能够沿着 model.fit 预测的线添加带有 CI 的功能区。而且我更喜欢在 ggplot 中执行此操作，特别是因为我想要在后台绘制意大利面条图。

Answer 1

在predict 中使用se.fit=TRUE：

library(tidyverse)
library(mgcv)

dt <- cbind(V1=scale(sample(1000)), 
    Age=rnorm(n = 1000, mean = 40, sd = 10), 
    ID=rep(seq(1:500),each=2)) %>% as.data.frame()

dt.model <- gamm(V1 ~ s(Age, bs="cr") + s(ID, bs = "re"), 
           data=dt, family="gaussian", discrete=T)

pred <- predict(dt.model$gam, se.fit=T)

dt %>% ggplot(aes(x = Age, y = V1)) +
   geom_line(aes(group=ID), alpha=.3) +
   geom_point(alpha=.2) +
   geom_ribbon(aes(ymin=pred$fit-1.96*pred$se.fit,
                   ymax=pred$fit+1.96*pred$se.fit), alpha=0.2, fill="red")+
   geom_line(aes(y=pred$fit), col="blue", lwd=1)