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

【问题标题】：How to get shaded confidence interval bands for glm coefficients?如何获得 glm 系数的阴影置信区间带？
【发布时间】：2020-07-14 22:06:27
【问题描述】：

我想绘制来自 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)

【问题讨论】：

标签： r polygon glm confidence-interval

【解决方案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")

【讨论】：

【解决方案2】：

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

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

【讨论】：