¿se llamaría intervalo de predicción? Tal vez estoy usando el término equivocado)
Eso no sería un intervalo de predicción. Un intervalo de predicción incorporaría la incertidumbre en la generación de datos. Es un poco inútil para una regresión logística binaria, ya que sabemos que el resultado será 0 o 1. Un intervalo de predicción puede ser más útil cuando se tienen datos de prueba (por ejemplo, predigo entre 8 y 12 eventos de los 20).
A calcular el intervalo de confianza. En primer lugar, hay que empezar por calcular $\operatorname{Var}(\operatorname{logit}(p))$ . Desde $\operatorname{logit}(p) = \beta_0 + \beta_1x$
$$\operatorname{Var}(\operatorname{logit}(p)) = \operatorname{Var}(\beta_0) + x^2 \operatorname{Var}(\beta_1) + 2x\operatorname{Cov}(\beta_0, \beta_1)$$
Puede encontrar estas cantidades a partir de la matriz de covarianza estimada a partir del modelo. Si se saca la raíz cuadrada de esta cantidad, se obtiene la desviación típica de $\operatorname{logit}(p)$ . Llamemos a esta cantidad $\sigma$ por ahora.
Un intervalo de confianza para $\operatorname{logit}(p)$ es $(\theta_L, \theta_U) = (\operatorname{logit}(p) - 1.96 \sigma,\operatorname{logit}(p) + 1.96 \sigma )$ . Invierta cada punto final para obtener un intervalo de confianza para $p$ .
Así es como se puede realizar esto en R.
library(tidyverse)
#Simulate some data to run a regression on
x = rnorm(100)
eta = -0.8 + 0.2*x
p = plogis(eta)
y = rbinom(100, 1, p)
# Fit a model
model = glm(y~x, family = binomial())
Bigma = vcov(model)
# Formula above for standard deviation of logit p. Vectorize for ease of computation
sig =Vectorize(function(x) sqrt(Bigma[1,1] + x^2*Bigma[2,2] + 2*x*Bigma[1,2]))
# New xs to make prediction on.
new_x = seq(-2, 2, 0.01)
logit_p = predict(model, newdata=list(x=new_x))
theta_L = logit_p - 1.96*sig(new_x)
theta_U = logit_p + 1.96*sig(new_x)
# Invert the estimates
p_L = plogis(theta_L)
p = plogis(logit_p)
p_U = plogis(theta_U)
d = data.frame(x=new_x, p_L, p, p_U)
# Now compare with builtin tools
fits = predict(model, newdata = list(x=new_x), se.fit = T) %>%
bind_cols() %>%
mutate(x=new_x) %>%
mutate(p = plogis(fit),
p_L = plogis(fit - 1.96*se.fit),
p_U = plogis(fit + 1.96*se.fit))
ggplot()+
# R's computation of the confidence interval in Black
geom_line(data=fits, aes(x, p))+
geom_ribbon(data=fits, aes(x=x, ymin=p_L, ymax=p_U), alpha = 0.5)+
# Out computation of the confidence interval in Red
geom_line(data=d, aes(x, p), color = 'red')+
geom_ribbon(data=d, aes(x=x, ymin=p_L, ymax=p_U), alpha = 0.5, fill='red')
Este es un ejemplo de salida
![enter image description here]()
El rojo y el negro están uno encima del otro, lo que significa que el procedimiento manual reproduce la salida de los métodos incorporados de R para calcular el error estándar del ajuste. Puedes ejecutar este código para múltiples semillas aleatorias y verificar que producen las mismas estimaciones (al menos hasta unos pocos decimales, imagino).
