如何在 ggplot 中向 glm 模型添加置信区间？

Question

这是我的数据的示例：

DATA <- data.frame(
TotalAbund = sample(1:10),
TotalHab = sample(0:1),
TotalInv = sample(c("yes", "no"), 20, replace = TRUE)
)
DATA$TotalHab<-as.factor(DATA$TotalHab)
DATA

这是我的模型：

MOD.1<-glm(TotalAbund~TotalInv+TotalHab, family=quasipoisson, data=DATA)

这是我的情节：

NEWDATA <- with(DATA,
               expand.grid(TotalInv=unique(TotalInv),
                       TotalHab=unique(TotalHab)))

pred <- predict(MOD.1,newdata= NEWDATA,se.fit=TRUE)
gg1 <- ggplot(NEWDATA, aes(x=factor(TotalHab), y=TotalAbund,colour=TotalInv))

我收到以下错误...

Error in eval(expr, envir, enclos) : object 'TotalAbund' not found

...尝试运行最后一行代码时：

gg1 + geom_point(data=pframe,size=8,shape=17,alpha=0.7,
             position=position_dodge(width=0.75))

有人可以帮忙吗？另外，我如何将 95% 的置信区间添加到我的分数中？谢谢。

Answer 1

您需要自己计算 95% 置信区间。您使用predict 并要求se.fit 是正确的。我们将首先询问链接尺度上的预测，计算 95% 置信区间，然后将它们转换为真实尺度进行绘图。这是一个方便的函数，用于计算日志链接（您在模型中使用的）的 CI。

# get your prediction
pred <- predict(MOD.1,newdata= NEWDATA,se.fit=TRUE,
            type = "link")

# CI function
make_ci <- function(pred, data){

# fit, lower, and upper CI
fit <- pred$fit
lower <- fit - 1.96*pred$se.fit
upper <- fit + 1.96*pred$se.fit

return(data.frame(exp(fit), exp(lower), exp(upper), data))
}

my_pred <- make_ci(pred, NEWDATA)

# to be used in geom_errorbar
limits <- aes(x = factor(TotalHab), ymax = my_pred$exp.upper., ymin = my_pred$exp.lower.,
                     group = TotalInv)

然后我们把它画出来，我会把最后的调整留给你，让你弄清楚你想要的样子。

ggplot(my_pred, aes(x = factor(TotalHab), y = exp.fit., color = TotalInv))+
geom_errorbar(limits, position = position_dodge(width = 0.75),
            color = "black")+
geom_point(size = 8, position = position_dodge(width = 0.75), shape = 16)+
ylim(c(0,15))+
geom_point(data = DATA, aes(x = factor(TotalHab), y = TotalAbund, colour = TotalInv),
         size = 8, shape = 17, alpha = 0.7,
         position = position_dodge(width = 0.75))