R Regresion Logistica

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. 1. Cargar datos y contamos las frecuencias de casos sin y con defectos:

datos <- read.table("http://www.uam.es/joser.berrendero/datos/challenger.txt", header = TRUE)
datos

Observe que los datos estan etiquetados en exito=1 y fracaso=0 variable categórica “Y” donde si se encontraron defectos=1 o no=0 en los propulsores ## Tabla de Frecuencias de los datos

tabla<- table(datos$defecto)
tabla

 0  1 
16  7 

Según la tabla anterior se presentarón 7 fallas de 23 ensayos o inspeciones y 16 de 23 no presentaron fallas en las turbinas ## Una representación gráfica de los datos, puede obtenerse mediante:

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

Parece razonable, a la vista de los datos, pensar que la temperatura puede influir en la probabilidad de que los propulsores tengan defectos o No. Luego la idea es ajustar por medio de la regresión logistica,un modelo de regresión logística para estudiar la posible relación entre Y=1 dado diferentes valores de temperatura: P(Y=1|X). Para ajustar el modelo se usa el comando glm (para modelos lineales generalizados) indicando que la respuesta es binomial mediante el argumento family:

Modelo Logistico en R

regresion <- glm(defecto ~ temp, data = datos, family = binomial)
summary(regresion)

Call:
glm(formula = defecto ~ temp, family = binomial, data = datos)

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

Analisis del Modelo

En el modelo de regresión logística: \(P(Y=1|X)=\dfrac{e^{15.0429-0.2322x}}{1+e^{15.0429-0.2322x}}\)

En la tabla coeficientes, la temperatura es significativa con pvalor= 0.0320 el cual es menor que el valor de significancia de 0.05. Lo anterior nos permite concluir que la esis nula se rechaza

\(H_0:\beta_{1}=0\) \(H_1:\beta_{1} \neq 0\)

Con respecto al intercepto: También es sinificativo dado que p valor= 0.0415< 0.05. Lo que significa que se rechaza la hipotesis nula, que afirma el intercepto es cero. Más aun, significa que bajísimos, lo que demuestra que ambas variables son importantes para explicar la variable dependiente (probabilidad de éxito).

\(H_0:\alpha_{0}=0\) \(H_1:\alpha_{0} \neq 0\)

¿La deviance residual con un valor de Residual deviance: 20.315 on 21 degrees of freedom de libertad, puede tomarse como indicio de buen ajuste del modelo?.

Los Odds Rations

Los “Odds” indica cuanto más probable es el éxito que el fracaso dado un vector de valores de variables explicativas (X’s). Por ejemplo, los Odds ratio representa el cociente de proporciones entre las fallas por cada exito entre el grupo de turbinas con un factor de riesgo y el grupo sin dicho factor de riesgo.

\(e^{\alpha_{0}}\)la probabilidad de que la característica está presente en una observación \(i\) cuando \(X_{i} 0\),es decir,al inicio del estudio. Es decir, cuanto más probable es el éxito que el fracaso. En particular, para este ejemplo es la probabilidad de que se presente el defecto cuando la temperatura es cero.

\(e^{\beta_{1}}\) es la razón de Odds al aumentar una unidad la temperatura, (manteniendo estables el resto de las variables en el modelo, si las hay)

exp(-0.2322 )

[1] 0.7927875

Es decir, al aumentar una unidad la temperatura, la Odds o ventaja de falla es casi lo mismo (aproximando a un);es 0.7927875 podríamos decir que el riesgo de falla frente al de no falla se multiplica por 0.8

Los Residuos

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 - p_i ). En la salida anterior estas cantidades se denominan deviance residuals. Para calcular estos pseudo-residuos, podemos ejecutar res = resid(reg).

res<- resid(regresion)
res

         1          2          3          4          5 
-1.0611168  1.7145343 -0.7996042 -0.8817559 -0.9690847 
         6          7          8          9         10 
-0.5865724 -0.5267644 -0.7229433  0.5506685  1.0063470 
        11         12         13         14         15 
 1.7145343 -0.3018707 -0.9690847  0.3540506 -0.9690847 
        16         17         18         19         20 
-0.4229076 -0.7229433 -0.2143127 -0.3782680 -0.2694144 
        21         22         23 
 2.2175345 -0.3782680  0.6127353 

Cuál observación presenta mayoy residual? También como en el modelo de regresión lineal se tienen los valores ajustados

 fit<- fitted.values(regresion)
fit

         1          2          3          4          5 
0.43049313 0.22996826 0.27362105 0.32209405 0.37472428 
         6          7          8          9         10 
0.15804910 0.12954602 0.22996826 0.85931657 0.60268105 
        11         12         13         14         15 
0.22996826 0.04454055 0.37472428 0.93924781 0.37472428 
        16         17         18         19         20 
0.08554356 0.22996826 0.02270329 0.06904407 0.03564141 
        21         22         23 
0.08554356 0.06904407 0.82884484 

Comparación

Teniendo en cuenta lo anterior cuál es la temperatura que tiene el mayor probabilidad de falla y corresponde a que temperatura? Las prediciones: predice los valores de los logits (valores de los logits ajustados)

prediciones<-predict(regresion)
prediciones

         1          2          3          4          5 
-0.2798395 -1.2084904 -0.9763277 -0.7441650 -0.5120022 
         6          7          8          9         10 
-1.6728159 -1.9049787 -1.2084904  1.8096252  0.4166488 
        11         12         13         14         15 
-1.2084904 -3.0657924 -0.5120022  2.7382762 -0.5120022 
        16         17         18         19         20 
-2.3693042 -1.2084904 -3.7622806 -2.6014669 -3.2979551 
        21         22         23 
-2.3693042 -2.6014669  1.5774625 

probabilidad= predict(regresion,type="response") 
probabilidad

         1          2          3          4          5 
0.43049313 0.22996826 0.27362105 0.32209405 0.37472428 
         6          7          8          9         10 
0.15804910 0.12954602 0.22996826 0.85931657 0.60268105 
        11         12         13         14         15 
0.22996826 0.04454055 0.37472428 0.93924781 0.37472428 
        16         17         18         19         20 
0.08554356 0.22996826 0.02270329 0.06904407 0.03564141 
        21         22         23 
0.08554356 0.06904407 0.82884484 

Generar una data frame que permita ver las variables originales con la probabilidad

comparacion<-data.frame( probabilidad,datos) 
comparacion

Predicción para valores nuevos con el modelo ajustado

Notas importantes: Tenga en cuenta que si ha usado en la fórmula nombres de variables que implican el data.frame que las contiene, el nuevo data.frame debe llamarse igual y con los mismos nombres de variables. Es decir, si la fórmula es a\(y~a\)x, para predecir debe usar el mismo nombre, a, para el data.frame, y por supuesto los mismos nombres de variables, aunque los valores sean distintos. Si usa sólo el nombre de las variables en la fórmula debe construir el data.frame nuevo con las columnas cuyos nombres representen a dichas variables. ##Ejemplo

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): 1. Primero generamos una nueva secuencia de datos, el cula llamaremos datos nuevos

datos_nuevos <- data.frame(temp = seq(50, 85, 0.1))
datos_nuevos

Calculo de las nuevas probabilidades

probabilidades_nuevas <- predict(regresion, datos_nuevos, type = "response")
probabilidades_nuevas

          1           2           3           4           5 
0.968773521 0.968063501 0.967337881 0.966596344 0.965838565 
          6           7           8           9          10 
0.965064218 0.964272967 0.963464474 0.962638394 0.961794377 
         11          12          13          14          15 
0.960932066 0.960051102 0.959151116 0.958231738 0.957292588 
         16          17          18          19          20 
0.956333284 0.955353437 0.954352653 0.953330531 0.952286667 
         21          22          23          24          25 
0.951220648 0.950132060 0.949020480 0.947885482 0.946726633 
         26          27          28          29          30 
0.945543495 0.944335626 0.943102578 0.941843899 0.940559130 
         31          32          33          34          35 
0.939247809 0.937909469 0.936543638 0.935149839 0.933727593 
         36          37          38          39          40 
0.932276413 0.930795812 0.929285297 0.927744372 0.926172536 
         41          42          43          44          45 
0.924569286 0.922934118 0.921266521 0.919565984 0.917831994 
         46          47          48          49          50 
0.916064035 0.914261590 0.912424138 0.910551162 0.908642138 
         51          52          53          54          55 
0.906696548 0.904713868 0.902693578 0.900635158 0.898538088 
         56          57          58          59          60 
0.896401851 0.894225931 0.892009815 0.889752993 0.887454957 
         61          62          63          64          65 
0.885115205 0.882733239 0.880308564 0.877840693 0.875329145 
         66          67          68          69          70 
0.872773443 0.870173122 0.867527720 0.864836788 0.862099883 
         71          72          73          74          75 
0.859316573 0.856486439 0.853609070 0.850684068 0.847711050 
         76          77          78          79          80 
0.844689643 0.841619491 0.838500253 0.835331602 0.832113229 
         81          82          83          84          85 
0.828844843 0.825526170 0.822156956 0.818736965 0.815265985 
         86          87          88          89          90 
0.811743821 0.808170305 0.804545289 0.800868649 0.797140287 
         91          92          93          94          95 
0.793360131 0.789528132 0.785644271 0.781708557 0.777721025 
         96          97          98          99         100 
0.773681741 0.769590801 0.765448331 0.761254490 0.757009466 
        101         102         103         104         105 
0.752713482 0.748366792 0.743969688 0.739522490 0.735025558 
        106         107         108         109         110 
0.730479285 0.725884100 0.721240466 0.716548885 0.711809895 
        111         112         113         114         115 
0.707024069 0.702192019 0.697314391 0.692391873 0.687425185 
        116         117         118         119         120 
0.682415087 0.677362374 0.672267881 0.667132476 0.661957064 
        121         122         123         124         125 
0.656742588 0.651490025 0.646200385 0.640874717 0.635514100 
        126         127         128         129         130 
0.630119647 0.624692505 0.619233851 0.613744894 0.608226871 
        131         132         133         134         135 
0.602681051 0.597108727 0.591511222 0.585889883 0.580246082 
        136         137         138         139         140 
0.574581214 0.568896695 0.563193963 0.557474474 0.551739701 
        141         142         143         144         145 
0.545991136 0.540230281 0.534458656 0.528677790 0.522889220 
        146         147         148         149         150 
0.517094496 0.511295171 0.505492803 0.499688956 0.493885192 
        151         152         153         154         155 
0.488083076 0.482284169 0.476490030 0.470702213 0.464922263 
        156         157         158         159         160 
0.459151719 0.453392110 0.447644951 0.441911745 0.436193982 
        161         162         163         164         165 
0.430493132 0.424810651 0.419147974 0.413506514 0.407887664 
        166         167         168         169         170 
0.402292795 0.396723249 0.391180346 0.385665379 0.380179611 
        171         172         173         174         175 
0.374724277 0.369300583 0.363909702 0.358552778 0.353230920 
        176         177         178         179         180 
0.347945206 0.342696676 0.337486341 0.332315173 0.327184110 
        181         182         183         184         185 
0.322094054 0.317045870 0.312040387 0.307078398 0.302160657 
        186         187         188         189         190 
0.297287883 0.292460756 0.287679922 0.282945986 0.278259519 
        191         192         193         194         195 
0.273621055 0.269031090 0.264490086 0.259998466 0.255556621 
        196         197         198         199         200 
0.251164905 0.246823638 0.242533104 0.238293555 0.234105211 
        201         202         203         204         205 
0.229968258 0.225882849 0.221849108 0.217867127 0.213936969 
        206         207         208         209         210 
0.210058667 0.206232227 0.202457624 0.198734810 0.195063709 
        211         212         213         214         215 
0.191444218 0.187876214 0.184359545 0.180894039 0.177479501 
        216         217         218         219         220 
0.174115715 0.170802445 0.167539433 0.164326403 0.161163064 
        221         222         223         224         225 
0.158049102 0.154984192 0.151967989 0.149000135 0.146080257 
        226         227         228         229         230 
0.143207970 0.140382874 0.137604559 0.134872602 0.132186571 
        231         232         233         234         235 
0.129546023 0.126950506 0.124399560 0.121892715 0.119429497 
        236         237         238         239         240 
0.117009424 0.114632005 0.112296749 0.110003155 0.107750719 
        241         242         243         244         245 
0.105538936 0.103367294 0.101235280 0.099142378 0.097088069 
        246         247         248         249         250 
0.095071835 0.093093156 0.091151509 0.089246374 0.087377230 
        251         252         253         254         255 
0.085543556 0.083744831 0.081980538 0.080250158 0.078553176 
        256         257         258         259         260 
0.076889080 0.075257357 0.073657498 0.072088999 0.070551357 
        261         262         263         264         265 
0.069044072 0.067566648 0.066118593 0.064699419 0.063308641 
        266         267         268         269         270 
0.061945779 0.060610357 0.059301904 0.058019954 0.056764043 
        271         272         273         274         275 
0.055533716 0.054328519 0.053148006 0.051991734 0.050859266 
        276         277         278         279         280 
0.049750171 0.048664022 0.047600398 0.046558883 0.045539068 
        281         282         283         284         285 
0.044540546 0.043562920 0.042605795 0.041668783 0.040751501 
        286         287         288         289         290 
0.039853572 0.038974624 0.038114291 0.037272214 0.036448035 
        291         292         293         294         295 
0.035641406 0.034851984 0.034079428 0.033323406 0.032583590 
        296         297         298         299         300 
0.031859657 0.031151291 0.030458178 0.029780014 0.029116496 
        301         302         303         304         305 
0.028467327 0.027832217 0.027210880 0.026603034 0.026008403 
        306         307         308         309         310 
0.025426717 0.024857707 0.024301114 0.023756680 0.023224153 
        311         312         313         314         315 
0.022703286 0.022193835 0.021695562 0.021208233 0.020731619 
        316         317         318         319         320 
0.020265494 0.019809637 0.019363831 0.018927865 0.018501529 
        321         322         323         324         325 
0.018084618 0.017676933 0.017278277 0.016888457 0.016507284 
        326         327         328         329         330 
0.016134573 0.015770143 0.015413814 0.015065414 0.014724771 
        331         332         333         334         335 
0.014391718 0.014066090 0.013747727 0.013436472 0.013132169 
        336         337         338         339         340 
0.012834669 0.012543823 0.012259486 0.011981516 0.011709773 
        341         342         343         344         345 
0.011444123 0.011184431 0.010930566 0.010682402 0.010439813 
        346         347         348         349         350 
0.010202675 0.009970870 0.009744280 0.009522789 0.009306286 
        351 
0.009094660 

Por defecto calcularía \(log p_{i}/(1-p_{i})\), para calcular \(p_{i}\) usamos el argumento type Interprete el dato 316 Ahora, si queremos ilustrar el ajuste con los datos

plot(datos$temp, datos$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_nuevos$temp, probabilidades_nuevas, col = "blue", lwd = 2)

Preguntas

1. Se puede afirmar a nivel \(\alpha = 0.05\) que la temperatura influye en la probabilidad de que los propulsores tengan defectos? Utilice el test de Wald.

2. Interpreta el valor del coeficiente estimado para la variable temperatura

3. ¿Para qué valores de la temperatura la probabilidad estimada de que se produzcan defectos es menor que 0.1?

4. ¿Para qué valores de la temperatura se predice que se van a producir defectos?