Regresión logística con R

En 1986, el transbordador espacial Challenger tuvo un accidente catastrófico debido a un incendio en una de las piezas de sus propulsores. Era la vez 25 en que se lanzaba un transbordador espacial. En todas las ocasiones anteriores se habían inspeccionado los propulsores de las naves, y en algunas de ellas se habían encontrando defectos. El fichero challenger contiene 23 observaciones de las siguientes variables: defecto, que toma los valores 1 y 0 en función de si se encontraron defectos o no en los propulsores; y temp, la temperatura (en grados Fahrenheit) en el momento del lanzamiento.

Primero leemos los datos y contamos las frecuencias de casos sin y con defectos:

challenger <- read.table('http://verso.mat.uam.es/~joser.berrendero/datos/challenger.txt', header = TRUE)
table(challenger$defecto)

## 
##  0  1 
## 16  7

Una representación gráfica de los datos, puede obtenerse mediante:

colores <- NULL
colores[challenger$defecto==0] <- 'green'
colores[challenger$defecto==1] <- 'red'
plot(challenger$temp, challenger$defecto, pch = 21, bg = colores, xlab = 'Temperatura', ylab = 'Probabilidad de defectos')
legend('bottomleft', c('No defecto', 'Si defecto'), pch = 21, col = c('green', 'red'))

Hemos usado los argumentos pch y bg para mejorar la apariencia del gráfico. También hemos usado el comando legend para incluir una leyenda explicativa.

Parece razonable, a la vista de los datos, pensar que la temperatura puede influir en la probabilidad de que los propulsores tengan defectos. En esta práctica, vamos a ajustar un modelo de regresión logística para estudiar la posible relación. Para ajustar el modelo se usa el comando glm (para modelos lineales generalizados) indicando que la respuesta es binomial mediante el argumento family:

reg <- glm(defecto ~ temp, data = challenger, family=binomial)
summary(reg)

## 
## Call:
## glm(formula = defecto ~ temp, family = binomial, data = challenger)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0611  -0.7613  -0.3783   0.4524   2.2175  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  15.0429     7.3786   2.039   0.0415 *
## temp         -0.2322     0.1082  -2.145   0.0320 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 28.267  on 22  degrees of freedom
## Residual deviance: 20.315  on 21  degrees of freedom
## AIC: 24.315
## 
## Number of Fisher Scoring iterations: 5

En el modelo de regresión logística la raíz de las desviaciones representa el papel de los residuos: $$D_i = \mp \sqrt{-2 [Y_i\log \hat p_i + (1-Y_i)\log(1-\hat p_i)]},$$ donde el signo coincide con el signo de $Y_i - \hat p_i$. En la salida anterior estas cantidades se denominan deviance residuals. Para calcular estos pseudo-residuos, podemos ejecutar res = resid(reg).

Para representar gráficamente la función logística estimada, calculamos las probabilidades de fallo estimadas (usando el comando predict) para un vector adecuado de nuevas temperaturas (entre 50 y 85 grados):

datos <- data.frame(temp = seq(50, 85, 0.1))
probabilidades <- predict(reg, datos, type = 'response')  # por defecto calcularía log p_i/(1-p_i), para calcular p_i usamos el argumento type 
plot(challenger$temp, challenger$defecto, pch = 21, bg = colores, xlab = 'Temperatura', ylab = 'Probabilidad de defectos')
legend('bottomleft', c('No defecto', 'Si defecto'), pch = 21, col = c('green', 'red'))
lines(datos$temp, probabilidades, col = 'blue', lwd = 2)