How to get shaded confidence interval bands for glm coefficients?

Question

I would like to plot the line and the shaded 95% confidence interval bands (for example using polygon)from a glm model (family binomial). For linear models (lm), I have previously been able to plot the confidence intervals from the predictions as they included the fit, lower and upper level see e.g. this answer How to plot regression transformed back on original scale with colored confidence interval bands? 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 Code and figure here:

# 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)

Ian Campbell Ian Campbell · Accepted Answer · 2020-04-02T20:31:21

Incorporating on @JamesCurran's answer, I believe this approach may work for you.

First, you use map2 from purrr to apply the prediction function to both models and extract the fit and the standard error. Then use mutate to add and subtract 1.96 times the standard error and transform. If you're unfamiliar with purrr, it's helpful to know that the ~ operator replaces function(x,y){} and makes the special objects .x and .y available.

Then we can use ggplot to plot the lines and confidence intervals.

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")

How to get shaded confidence interval bands for glm coefficients?

2 Answers