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

Question

我想从 glm 模型（二项式）中绘制直线和阴影 95% 置信区间带（例如使用多边形）。对于线性模型 (lm)，我之前已经能够根据预测绘制置信区间，因为它们包括拟合、下限和上限，请参见例如这个答案 but I do not know how to do it here. Thanks for help in advance. You can find here the data that I used (it contains 3 variables and 4582 observations): 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

predict.glm 有一个可选参数 se.fit，通常设置为 FALSE。将其设置为 TRUE，然后您可以使用预测 +/- 1.96 * std.error 来计算您的 Wald 置信区间。

绘制它们 - 取决于您是否需要线条或阴影区域，但是 lines 或 polygon 应该可以满足您的要求。

Answer 2

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

首先，您使用 purrr 中的 map2 将预测函数应用于两个模型并提取拟合和标准误差。然后用mutate加减1.96倍的标准误和transform。如果您不熟悉 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")

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

How to get shaded confidence interval bands for glm coefficients?

r

polygon

confidence-interval

glm