18 votos

Trazado de los intervalos de confianza para la probabilidad predicha de una regresión logística

Estadísticas y R newbie aquí.

OK, tengo una regresión logística y ha utilizado la función de predecir para desarrollar una curva de probabilidad basada en mis cálculos.

## 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 es genial pero tengo curiosidad sobre el trazado de los intervalos de confianza para las probabilidades. He intentado plot.ci() pero no tuvo mucha suerte. Puede alguien me punto a algunas maneras de conseguir esto, preferiblemente con el paquete de coche o base R.

Gracias.

23voto

Jeff Davis Puntos1999

El código que utilizan las estimaciones de un modelo de regresión logística utilizando el glm función. No se incluyen los datos, así que voy a hacer 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 modelos de la relación entre una variable respuesta binaria y, en este caso, una continua predictor. El resultado es un logit-transormed la probabilidad como una relación lineal para el factor de predicción. En su caso, el resultado es una respuesta binaria correspondiente a ganar o no ganar en el juego de azar y es que es predicho por el valor de la apuesta. Los coeficientes de mod1 se dan en la sesión de probabilidades (que son difíciles de interpet), de acuerdo a:

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

To convert logged odds to probabilities, we can translate the above to

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

Usted puede utilizar esta información para configurar la trama. Primero, usted necesita un rango de la variable de predicción:

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

A continuación, utilizando predict, usted puede obtener las predicciones basadas en el modelo

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

Tenga en cuenta que los valores ajustados también pueden obtenerse a través de

mod1$fitted

Especificando se.fit=TRUE, usted también consigue el error estándar asociado con cada uno equipado valor. El resultado data.frame es una matriz con los siguientes componentes: el conjunto de predicciones (fit), la estimación de los errores estándar (se.fit), y un escalar dando la raíz cuadrada de la dispersión se utiliza para calcular los errores estándar (residual.scale). En el caso de una binomial logit, el valor será 1 (que se puede ver ingresando preddat$residual.scale en R). Si desea ver un ejemplo de lo que hemos calculado hasta ahora, puede escribir head(data.frame(preddat)).

El siguiente paso es configurar la trama. Me gusta para configurar un espacio en blanco área de trazado con los parámetros de:

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

Ahora usted puede ver donde es importante saber cómo calcular el conjunto de probabilidades. Usted puede dibujar la línea correspondiente a la equipada con probabilidades después de la segunda fórmula de arriba. El uso de la preddat data.frame usted puede convertir los valores ajustados a las probabilidades y el uso que trazar una línea en contra de los valores de la variable predictora.

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

Finalmente, la respuesta a su pregunta, los intervalos de confianza pueden ser añadidos a la trama mediante el cálculo de la probabilidad para 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))

La gráfica resultante (de hecho al azar de datos) debe ser algo como esto:

enter image description here

Para la conveniencia del amor, aquí está todo el código en un fragmento:

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 gran editado respuesta en un intento de hacer más pertinentes para las estadísticas.stackexchange.)

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by: