Trazado de intervalos de confianza para las probabilidades predichas de una regresión logística

Question

Ok, tengo una regresión logística y he utilizado el predict() para desarrollar una curva de probabilidad basada en mis estimaciones.

## LOGIT MODEL:
library(car)
mod1 = glm(factor(won) ~ as.numeric(bid), data=mydat, family=binomial(link="logit"))

## PROBABILITY CURVE:
all.x <- expand.grid(won=unique(won), bid=unique(bid))
y.hat.new <- predict(mod1, newdata=all.x, type="response")
plot(bid<-000:1000,predict(mod1,newdata=data.frame(bid<-c(000:1000)),type="response"), lwd=5, col="blue", type="l")

Esto está muy bien, pero tengo curiosidad por saber cómo trazar los intervalos de confianza de las probabilidades. He intentado plot.ci() pero no ha tenido suerte. ¿Puede alguien indicarme alguna forma de conseguirlo, preferiblemente con el car paquete o base R.

Answer 1

El código que ha utilizado estima un modelo de regresión logística utilizando el glm función. No has incluido datos, así que me inventaré algunos.

set.seed(1234)
mydat <- data.frame(
    won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
    bid=runif(250, min=0, max=1000)
)
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))

Un modelo de regresión logística modela la relación entre una variable de respuesta binaria y, en este caso, un predictor continuo. El resultado es un logit-transformado probabilidad como una relación lineal con el predictor. En su caso, el resultado es una respuesta binaria que corresponde a ganar o no ganar en el juego y está siendo predicho por el valor de la apuesta. Los coeficientes de mod1 se dan en probabilidades registradas (que son difíciles de interpretar), según:

$$\text{logit}(p)=\log\left(\frac{p}{(1-p)}\right)=\beta_{0}+\beta_{1}x_{1}$$

Para convertir las probabilidades registradas en probabilidades, podemos traducir lo anterior en

$$p=\frac{\exp(\beta_{0}+\beta_{1}x_{1})}{(1+\exp(\beta_{0}+\beta_{1}x_{1}))}$$

Puedes utilizar esta información para configurar la parcela. En primer lugar, necesita un rango de la variable de predicción:

plotdat <- data.frame(bid=(0:1000))

A continuación, utilizando predict puede obtener predicciones basadas en su modelo

preddat <- predict(mod1, newdata=plotdat, se.fit=TRUE)

Obsérvese que los valores ajustados también pueden obtenerse mediante

mod1$fitted

Al especificar se.fit=TRUE También se obtiene el error estándar asociado a cada valor ajustado. El resultado data.frame es una matriz con los siguientes componentes: las predicciones ajustadas ( fit ), los errores estándar estimados ( se.fit ), y un escalar que da la raíz cuadrada de la dispersión utilizada para calcular los errores estándar ( residual.scale ). En el caso de un logit binomial, el valor será 1 (que se puede ver introduciendo preddat$residual.scale en R ). Si quiere ver un ejemplo de lo que ha calculado hasta ahora, puede escribir head(data.frame(preddat)) .

El siguiente paso es preparar el terreno. A mí me gusta configurar primero un área de ploteo en blanco con los parámetros:

with(mydat, plot(bid, won, type="n", 
    ylim=c(0, 1), ylab="Probability of winning", xlab="Bid"))

Ahora puedes ver dónde es importante saber cómo calcular las probabilidades ajustadas. Puedes dibujar la línea correspondiente a las probabilidades ajustadas siguiendo la segunda fórmula anterior. Utilizando la preddat data.frame puede convertir los valores ajustados en probabilidades y utilizarlos para trazar una línea contra los valores de su variable de predicción.

with(preddat, lines(0:1000, exp(fit)/(1+exp(fit)), col="blue"))

Por último, responde a tu pregunta, los intervalos de confianza pueden añadirse al gráfico calculando la probabilidad de los valores ajustados +/- 1.96 veces el error estándar:

with(preddat, lines(0:1000, exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)), lty=2))
with(preddat, lines(0:1000, exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit)), lty=2))

El gráfico resultante (a partir de los datos generados aleatoriamente) debería tener el siguiente aspecto:

Por conveniencia, aquí está todo el código en un trozo:

set.seed(1234)
mydat <- data.frame(
    won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
    bid=runif(250, min=0, max=1000)
)
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))
plotdat <- data.frame(bid=(0:1000))
preddat <- predict(mod1, newdata=plotdat, se.fit=TRUE)
with(mydat, plot(bid, won, type="n", 
    ylim=c(0, 1), ylab="Probability of winning", xlab="Bid"))
with(preddat, lines(0:1000, exp(fit)/(1+exp(fit)), col="blue"))
with(preddat, lines(0:1000, exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)), lty=2))
with(preddat, lines(0:1000, exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit)), lty=2))

(Nota: Esta es una respuesta muy editada en un intento de hacerla más relevante para stats.stackexchange).

Answer 2

Aquí hay una modificación de la solución de @smillig. Aquí utilizo las herramientas de tidyverse, y también uso el linkinv que forma parte del objeto del modelo GLM mod1 . De esta manera, no tiene que invertir manualmente la función logística, y este enfoque funcionará independientemente del GLM específico que ajuste.

library(tidyverse)
library(magrittr)

set.seed(1234)

# create fake data on gambling. Does prob win depend on bid size? 
mydat <- data.frame(
  won=as.factor(sample(c(0, 1), 250, replace=TRUE)), 
  bid=runif(250, min=0, max=1000)
)

# logistic regression model: 
mod1 <- glm(won~bid, data=mydat, family=binomial(link="logit"))

# new predictor values to use for prediction: 
plotdat <- data.frame(bid=(0:1000))

# df with predictions, lower and upper limits of CIs: 
preddat <- predict(mod1,
               type = "link",
               newdata=plotdat,
               se.fit=TRUE) %>% 
  as.data.frame() %>% 
  mutate(bid = (0:1000), 

         # model object mod1 has a component called linkinv that 
         # is a function that inverts the link function of the GLM:
         lower = mod1$family$linkinv(fit - 1.96*se.fit), 
         point.estimate = mod1$family$linkinv(fit), 
         upper = mod1$family$linkinv(fit + 1.96*se.fit)) 

# plotting with ggplot: 
preddat %>% ggplot(aes(x = bid, 
                   y = point.estimate)) + 
  geom_line(colour = "blue") + 
  geom_ribbon(aes(ymin = lower,
                  ymax = upper), 
              alpha = 0.5) + 
  scale_y_continuous(limits = c(0,1))