Как получить закрашенные полосы доверительного интервала для коэффициентов glm?

Question

Я хотел бы построить линию и заштрихованные полосы 95% доверительного интервала (например, используя многоугольник) из модели glm (семейный бином). Для линейных моделей (лм) ранее я был в состоянии построить доверительные интервалы из прогнозов, поскольку они включали подгонку, нижний и верхний уровни, см., Например, этот ответ Как построить регрессию, преобразованную обратно в исходную шкалу с цветным доверительным интервалом группы? но я не знаю, как это сделать здесь. Спасибо за помощь заранее. Здесь вы можете найти данные, которые я использовал (он содержит 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, я считаю, что этот подход может сработать для вас.

Сначала вы используете map2 из purrr, чтобы применить функцию прогнозирования к обеим моделям и извлечь подгонку и стандартную ошибку. Затем используйте 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")

Answer 2

predict.glm имеет необязательный аргумент se.fit, который обычно устанавливается на FALSE. Установите значение TRUE, а затем вы можете использовать прогноз +/- 1,96 * std.error для расчета доверительных интервалов Вальда.

Построение графика - зависит от того, хотите ли вы линии или заштрихованные области, но lines или polygon должны вас охватить.