r - 如何获得 glm 系数的阴影置信区间带？

Question

我想从 glm 模型(族二项式)绘制线条和阴影 95% 置信区间带(例如使用多边形)。对于线性模型(lm)，我之前已经能够根据预测绘制置信区间，因为它们包括拟合、较低和较高水平，请参见例如这个答案How to plot regression transformed back on original scale with colored confidence interval bands?但我不知道在这里该怎么做。提前感谢您的帮助。您可以在这里找到我使用的数据(它包含 3 个变量和 4582 个观察值):https://drive.google.com/file/d/1RbaN2vvczG0eiiqnJOKKFZE9GX_ufl7d/view?usp=sharing代码和图在这里:

# Models
hotglm=glm(hotspot~age+I(age^2),data = data, family = "binomial")
summary(hotglm)

coldglm=glm(coldspot~age+I(age^2),data = data, family = "binomial")
summary(coldglm)

# Plot
age = 1:200
lin=hotglm$coefficients[1]+hotglm$coefficients[2]*age+hotglm$coefficients[3]*age^2
pr = exp(lin)/(1+exp(lin))
par(mfrow=c(1,1))
plot(age, pr,type="l",col=2,lwd=2,ylim=c(0,.15))

lin=coldglm$coefficients[1]+coldglm$coefficients[2]*age+coldglm$coefficients[3]*age^2
pr = exp(lin)/(1+exp(lin))
lines(age, pr,type="l",col="blue", lwd=2)

Answer 1

结合@JamesCurran的回答，我相信这种方法可能适合您。

首先，您使用 purrr 中的 map2 将预测函数应用于两个模型并提取拟合值和标准误差。然后使用mutate加减1.96倍标准误差并变换。如果您不熟悉 purrr，了解 ~ 运算符会替换 function(x,y){} 并使特殊的对象 .x 和 .y 可用。

然后我们可以使用ggplot来绘制线条和置信区间。

library(tidyverse)
library(ggplot2)
hotglm <- glm(hotspot~age+I(age^2),data = data, family = "binomial")
coldglm <- glm(coldspot~age+I(age^2),data = data, family = "binomial")

plotdata <- map2(list(coldfit = coldglm,coldse = coldglm,hotfit = hotglm, hotse = hotglm),
     rep(c("fit","se.fit"),times=2),
     ~ predict(.x,data.frame(age=1:200),se.fit = TRUE)[[.y]]) %>%
        data.frame %>%
        mutate(age = 1:200,
         coldline = exp(coldfit)/(1+exp(coldfit)),       
         coldlower = exp(coldfit - (coldse * 1.96))/(1+exp(coldfit - (coldse * 1.96))),
         coldupper = exp(coldfit + (1.96 * coldse))/(1+exp(coldfit + (1.96 * coldse))),
         hotline = exp(hotfit)/(1+exp(hotfit)),       
         hotlower = exp(hotfit - (1.96 * hotse))/(1+exp(hotfit - (1.96 * hotse))),
         hotupper = exp(hotfit + (1.96 * hotse))/(1+exp(hotfit + (1.96 * hotse))))

ggplot(plotdata,aes(x=age,y=coldline)) +
  geom_line(color = "blue") +
  geom_line(aes(y=hotline),color="red")

ggplot(plotdata,aes(x=age,y=coldline)) +
  geom_line(color = "blue") +
  geom_ribbon(aes(ymin=coldlower, ymax=coldupper), alpha = 0.2,fill = "blue") +
  geom_line(aes(y=hotline),color="red") + geom_ribbon(aes(ymin=hotlower, ymax=hotupper), alpha = 0.2,fill = "red")