Analissidel accidente del transbordador espacial challenger
- Antecedentes
poner imagen del challenger ## Regresión logistica
challenger <- read.table("http://verso.mat.uam.es/~joser.berrendero/datos/challenger.txt", header = TRUE)
table(challenger$defecto)##
## 0 1
## 16 7
- Representar los defectos vs no defestos de forma gráfica
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"))- 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ó 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 devíance residuals. Para calcular estos pseudo-residuos, podemos ejecutar res = resid(reg).
res <- resid(reg)
res## 1 2 3 4 5 6 7
## -1.0611168 1.7145343 -0.7996042 -0.8817559 -0.9690847 -0.5865724 -0.5267644
## 8 9 10 11 12 13 14
## -0.7229433 0.5506685 1.0063470 1.7145343 -0.3018707 -0.9690847 0.3540506
## 15 16 17 18 19 20 21
## -0.9690847 -0.4229076 -0.7229433 -0.2143127 -0.3782680 -0.2694144 2.2175345
## 22 23
## -0.3782680 0.6127353
- Pasaremos a calcular 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")
#Gráfia
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=3)