Regresión Lineal

  1. Con la intención de comparar el desempeño de dos clases de discos duros (0 : SDD, 1: HDD). Este desempeño es medido a través de la variable Y: tiempo de respuesta del disco (segundos), la cual se relaciona, posiblemente bajo una dependencia no lineal, de X: la carga del sistema (Número de consultas por minuto).

    Se han realizado múltiples ensayos bajo ambas configuraciones y bajo variación de la carga del sistema. Los resultados se presentan en la siguiente tabla:

    Conf Carga Tiempo
    1 1.0 0.9
    0 2.0 0.3
    1 2.4 2.0
    0 3.1 0.8
    1 4.0 2.7
    1 4.3 2.6
    0 5.8 2.5
    0 6.6 3.2
    0 7.5 3.7
    1 8.0 3.9
    0 9.0 5.3
    1 9.2 4.2
    1 10.2 3.9

    Conf Carga Tiempo
    1 1.8 1.1
    1 2.0 1.5
    0 2.5 0.5
    0 3.9 1.5
    0 4.2 1.6
    1 5.5 3.3
    0 6.4 3.3
    1 7.0 3.5
    0 8.0 4.3
    1 8.2 4.0
    1 9.1 4.3
    0 9.5 5.8
    1. Represente gráficamente la relación observada entre el tiempo de respuesta y la carga de trabajo, para los dos tipos de disco duro. ¿Se evidencia una relación lineal? Mida la fuerza de esta relación para ambos tipos de disco a través de los coeficientes de correlación.

      Se realiza gráfica de comparación para ambos tipos de dicos:

      Ambos tipos de disco muestran una relación lineal entre la carga y el tiempo de respuesta de los dispositivos. Los valores de correlación, por encima de 0.95 indican una fuerza de relación lineal fuerte. Los p-value, que son muy pequeños, muestran una alta significancia de la relación lineal, descartando así la no correlación de las variables.

    2. Ajuste un primer modelo de regresión simple (Modelo 1) que reproduzca la relación entre la carga y el tiempo de respuesta, sin incluir la configuración del disco duro. Evalúe la bondad de ajuste de este modelo e interprete los resultados obtenidos.

      ## 
      ## Call:
      ## lm(formula = tiempo ~ carga, data = data)
      ## 
      ## Residuals:
      ##      Min       1Q   Median       3Q      Max 
      ## -1.16824 -0.40281 -0.03945  0.43541  1.07627 
      ## 
      ## Coefficients:
      ##             Estimate Std. Error t value Pr(>|t|)    
      ## (Intercept)  0.04838    0.26321   0.184    0.856    
      ## carga        0.49214    0.04177  11.783 3.18e-11 ***
      ## ---
      ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      ## 
      ## Residual standard error: 0.5837 on 23 degrees of freedom
      ## Multiple R-squared:  0.8579, Adjusted R-squared:  0.8517 
      ## F-statistic: 138.8 on 1 and 23 DF,  p-value: 3.177e-11

      Se puede observar un p-value = 3.177e-11, un valor muy pequeño, lo cual indica que el modelo explica gran parte de la variabilidad del tiempo de respuesta del disco. R-squared es el porcentaje de la variabilidad del tiempo de respuesta, explicada por el modelo, 85%. Tambien podemos observar que mientras la carga es importante para explicar la variabilidad del tiempo de respuesta, el intercepto no lo es tanto, lo que significa que cuando la carga del disco tiende a cero, el modelo no explica muy bien el tiempo de respuesta del disco.

    3. Obtenga un nuevo modelo (Modelo 2) en el que incluya el tipo de disco (Variable Dummy) y su interacción con la carga del equipo. Evalué la bondad de ajuste del nuevo modelo, e interprete los coeficientes del Modelo 2. Recom. Note que la pendiente y el intercepto no son los mismos para los dos tipos de discos

      Dado que la pendiente y el intercepto lucen diferente para los dos tipos de discos, propondría el siguiente modelo.

      \[ Tiempo = (1-conf)(c1 + c2*carga) + conf*(c3 + c4*carga) + \epsilon \]

      Donde \(conf=1\) cuando el disco es de tipo HDD y \(conf=0\) cuando el disco es de tipo SSD.

      \(c1\) y \(c2\) son el intercepto y la razon de cambio de la carga, cuando el disco es de tipo SSD.

      \(c3\) y \(c4\) son el intercepto y la razon de cambio de la carga, cuando el disco es de tipo HDD.

      Sin embargo hay que hacer alguna transformación a este modelo para quede de la forma

      \[ y=\beta_0 + \beta_1*x_1 + \beta_2*x_2 + \beta_3*x_1*x_2 + \epsilon \]

      que es adecuada para el algoritmo de regresión lineal. En el siguiente bloque se muestra esta transformación.

      \[ Tiempo = c1*(1-conf) + c2*(1-conf)*carga + c3*conf + c4*conf*carga + \epsilon \] \[ Tiempo = c1 - c1 *conf + c2*carga - c2*conf*carga + c3*conf + c4*conf*carga + \epsilon \] \[ Tiempo = c1 + (c3-c1)*conf + c2*carga + (c4-c2)*conf*carga + \epsilon \]

      Si hacemos los siguientes reemplazos \(\beta_0=c1\), \(\beta_1=c3-c1\), \(\beta_2=c2\), y \(\beta_3=c4-c2\) , tendriamos el siguiente modelo:

      \[ Tiempo = \beta_0 + \beta_1*conf + \beta_2*carga + \beta_3*conf*carga + \epsilon \]

      El cual puede ser utilizado para que el algoritmo de regresión lineal calcule el valor de las contantes \(\beta\).

      Se aplica el modelo sobre los datos y se obtiene los siguientes resultados:

      ## 
      ## Call:
      ## lm(formula = tiempo ~ conf + carga + conf * carga, data = data)
      ## 
      ## Residuals:
      ##      Min       1Q   Median       3Q      Max 
      ## -0.68547 -0.11333  0.06881  0.15302  0.41807 
      ## 
      ## Coefficients:
      ##               Estimate Std. Error t value Pr(>|t|)    
      ## (Intercept)   -1.37549    0.20902  -6.581 1.62e-06 ***
      ## confHDD        2.26391    0.26520   8.536 2.86e-08 ***
      ## carga          0.71979    0.03367  21.376 9.88e-16 ***
      ## confHDD:carga -0.35734    0.04227  -8.454 3.36e-08 ***
      ## ---
      ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      ## 
      ## Residual standard error: 0.2844 on 21 degrees of freedom
      ## Multiple R-squared:  0.9692,    Adjusted R-squared:  0.9648 
      ## F-statistic: 220.2 on 3 and 21 DF,  p-value: 5.042e-16

      De aquí podemos observar que el p-value es muy pequeño, lo que argumenta la relevancia del modelo, respecto de la hipotesis nula. F-statistic es mayor en el modelo 2 que en el anterior modelo, lo cual muestra una mejora en la explicabilidad de la varibilidad del tiempo de respuesta. R-cuadrado también muestra un mejor rendimiento de este modelo, con 96% de la explicabilidad del tiempo de carga, en comparación con el anterior modelo que solo explicaba el 85%. Tambien se puede observar que todos los componentes del modelo son estimados como relevantes.

      Teniendo en cuenta las definiciones de las constantes \(\beta\) y las constantes \(c\), podemos interpretar que:

      • \(c1=\beta_0=-1.37\) es el valor promedio del tiempo de respuesta cuando no hay carga y es un disco tipo SSD.

      • \(c2=\beta_2=0.72\) es la razón de cambio del tiempo de respuesta respecto a la carga, cuando es un disco tipo SSD

      • \(c3=\beta_1+c1=2.26-1.37=0.89\) es el valor promedio del tiempo de respuesta cuando no hay carga y es un disco tipo HDD

      • \(c4=\beta_3+c2=-0.36+0.72=0.36\) es la razón de cambio del tiempo de respuesta respecto de la carga, cuando el disco es de tipo HDD.

    4. Mediante el test ANOVA correspondiente, pruebe que la inclusión de la variable cualitativa configuración del disco y su interacción con la carga mejora significativamente el ajuste del modelo.

      Se realiza un test ANOVA y se obtienen los siguientes resultados:

      ## Analysis of Variance Table
      ## 
      ## Model 1: tiempo ~ carga
      ## Model 2: tiempo ~ conf + carga + conf * carga
      ##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
      ## 1     23 7.8375                                  
      ## 2     21 1.6990  2    6.1386 37.938 1.067e-07 ***
      ## ---
      ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      • Se observa que el modelo 2 a pesar de que tiene menos grados de libertad, puesto que usa más parametros que el modelo 1, en realidad tiene mucho menos RSS, que mide la variabilidad de que no es explicada por el modelo. También podemos ver un p-value muy pequeño, lo que indica que la mejora de pasar del modelo 1 al modelo 2 es estadisticamente significativa.
    5. Represente gráficamente el ajuste del Modelo 2 y evalúe el cumplimiento de los supuestos sobre el termino error.

      El siguiente gráfico muestra los valores predichos, comprobando así que el modelo tiene una tendencia promedio similar a los datos de origen.

      Sin embargo parece que las tendencias no fuesen lineales del todo. Mientras que para los discos tipo SSD, el tiempo de respuesta pare crecer exponencialmente respecto a la carga, por otro lado, para los discos HDD el tiempo de espera parece crecer logaritmicamente respecto de la carga.

      Los siguientes gráficos muestran los residuos estandarizados (rstudent) comparados con los valores ajustados del tiempo esperado, la carga, y el tipo de disco.

      Esto nos muestra que varios de los supuestos acerca del error no se están cumpliendo para este modelo:

      • se puede observar que los residuos no tienen una media igual 0 en todos los valores predichos del tiempo de respuesta, lo cual indica que no hay homocedásticidad en los residuos.

      • se puede observar una tendencia curva negativa entre los residuos y la carga del disco. Lo que significa que hay un componente no lineal aún no descubierto en la variabilidad del tiempo de espera relacionado con la carga del disco. Esto tambien nos indica que la varianza de los residuos no es constante y por lo tanto los residuos no son homocedásticos.

      • Respecto al tipo de disco, los residuos parecen distribuirse normalmente, con media=0, para los discos tipo SSD, pero para los discos HDD no hay asimetría y la media es mayor que 0. Esto refuerza la idea de que el error es heterocedástico (incumpliendo nuevamente uno de los supuestos del error). También significa que el modelo se comporta bien para predecir el tiempo de espera de los discos SSD.

      • En las siguientes gráficas se puede observar efectivamente que no hay homocedásticidad interna ni para los discos SSD ni para los discos HDD. Para los discos SSD aquí se puede encontrar de nuevo evidencias de que existe una aparente relación exponencial y para los discos HDD hay una aparente relación logaritmica.

      • Se realiza test de media = 0 para el error:

        ## 
        ##  One Sample t-test
        ## 
        ## data:  data$residuos_rstudent
        ## t = -0.14652, df = 24, p-value = 0.8847
        ## alternative hypothesis: true mean is not equal to 0
        ## 95 percent confidence interval:
        ##  -0.4994321  0.4332192
        ## sample estimates:
        ##   mean of x 
        ## -0.03310641

        El p-value de 0.88, nos indica en este caso que no se puede rechazar la hipotesis nula, la cual es que la media de los residuos es 0. Este test indica que este supuesto si se cumple. El intervalo de confianza incluye el 0 y la media estimada también es cercana a 0.

      • Se realiza test de Homogeneidad de varianza:

        ## 
        ##  studentized Breusch-Pagan test
        ## 
        ## data:  modelo_2
        ## BP = 2.6825, df = 3, p-value = 0.4432

        Un p-value alto en esta prueba indica que no se puede rechar la hipotesis nula, que en este caso siginifica que no hay evidencia clara para rechazar la hipotesis de que la varianza de los residuos es constante y no está relacionada con ninguna de las variables. Debemos tener en cuenta que deacuerdo a los gráficos mostrados anteriormente, hay patrones heterocedasticidad que este test no detectó.

      • Se realiza test de independencia de los residuos

        ## 
        ##  Durbin-Watson test
        ## 
        ## data:  modelo_2
        ## DW = 1.3285, p-value = 0.03421
        ## alternative hypothesis: true autocorrelation is greater than 0

        Este test con un p-value muy pequeño, rechaza la hipotesis nula, es decir, rechaza la hipetesis de independencia de los residuos. El valor DW < 2, indica que hay una autocorrelación positiva entre los residuos del modelo.

      • Se realiza test de normalidad de los residuos

        ## 
        ##  Shapiro-Wilk normality test
        ## 
        ## data:  data$residuos_rstudent
        ## W = 0.92407, p-value = 0.06348

        En este test se puede observar un p-value que aunque es mayor a 0.05, en realidad está muy cercano a este humbral. Junto con el valor de W=0.92, podríamos concluir en terminos prácticos los residuos del modelo son razonablemente normales. En todo caso, como ya se pudo observar en gráficas anteriores, el modelo aún se podría mejorar para tener suficiente evidencia de la distribución normal de los residuos.

    6. Concluya de forma general.

      1. El modelo 2 nos permite acercarnos mucho mejor al ajuste del tiempo de espera, dado que considera interceptos y razones de cambios diferentes para cada tipo de disco.

      2. El modelo 2 aun puede ser mejorado para ajustarse mejor a la variabilidad del tiempo de espera. En particular parece ser que la relación del tiempo de espera para los discos SSD es exponencial, mientras que para los discos HDD es logarimica

        • En realidad me sorprende esta hipótesis, los discos tipo SSD al ser de una tecnología más reciente, parecido a la usada en las memorias usb, son discos más rápidos. Se esperaría un tiempo de espera lineal, si no, logaritmo para este tipo de discos. Por su contra parte, los discos tipo HDD son de una tecnología antigua que usaba discos mecanicos para la lectura y escritura de datos. Se esperaría un tiempo de espera que crece de forma exponencial este segundo tipo de discos.
        • Lo que podría estar pasando es que hubo un error en la recolección de datos y se invirtieron las etiquetas de los tipo de discos. La recomendación sería verificar (y corregir si es el caso) los metodos de recolección de datos y recolectar un volumen mayor de datos para asegurarnos de capturar mejor la tendencia.
  2. Una compañía de seguros de automóvil desea caracterizar la siniestralidad de sus asegurados durante el último año. Para ello dispone información de una muestra aleatoria de 35 asegurados con la siguiente información (accidentes.xlsx):

    Acc: haber tenido algún accidente en el último año (0:no; 1:sí).
    Exp: años de experiencia.
    Edad: edad del conductor.
    Pot: potencia del motor.
    Sexo: 1 (mujer), 2 (hombre).
    1. Con herramientas del análisis exploratorio, estudie la asociación entre la siniestralidad y el conjunto de variables predictoras (Edad, Experiencia, Potencia del motor y Sexo).

      Importamos el conjunto de datos y miramos su estructura.

      ## tibble [35 × 5] (S3: tbl_df/tbl/data.frame)
      ##  $ Acc : num [1:35] 0 0 0 1 0 1 0 1 0 1 ...
      ##  $ Exp : num [1:35] 10 15 7 1 10 2 8 20 18 4 ...
      ##  $ Edad: num [1:35] 30 40 25 21 29 20 40 25 43 23 ...
      ##  $ Pot : num [1:35] 90 85 95 145 70 120 95 135 85 110 ...
      ##  $ Sexo: num [1:35] 1 1 1 2 1 2 1 2 1 2 ...

      Luego aseguramos que la variabla Acc y la variable Sexo sean categoricas y vemos un resumen de los datos.

      ##  Acc          Exp              Edad         Pot            Sexo   
      ##  no:20   Min.   : 1.000   Min.   :20   Min.   : 70.0   mujer :21  
      ##  si:15   1st Qu.: 6.500   1st Qu.:25   1st Qu.: 90.0   hombre:14  
      ##          Median : 9.000   Median :29   Median : 95.0              
      ##          Mean   : 9.543   Mean   :31   Mean   :101.6              
      ##          3rd Qu.:12.000   3rd Qu.:36   3rd Qu.:110.0              
      ##          Max.   :20.000   Max.   :56   Max.   :150.0

      Los siguientes gráficos nos proveen un resumen estadistico univarido y bivariado del conjunto de datos de accidentes, identificando por color los datos de hombres versus los de mujeres.

      En la diagonal tenemos el análisis univariado de cada una de las varibles. Podemos observar que se tienen registrados casos ligeramente desbalanciados hacia los casos de no accidente. Así como tambien un ligero desbalance a tener más datos casos de mujeres que de hombres. Esto puede estar relacionado con el algoritmo de recolección utilizado. Se podría realizar varios experimentos de recolección de nuevos datos aleatorios, para garantizar que efectivamente se está modelando la tendencia.

      Respecto a la experiencia, se puede observar que las mujeres, en los datos recolectados, tienden a tener 9 años de experiencia, mientras que en los hombres parece tener una distribución de años de experiencia un poco más uniforme, aunque con una leve tendencia hacia casos de hombres con 4 años de experiencia.

      Respecto de la edad, la tendencia de hombres y mujeres esta en el rango entre los 25 y 35 años de edad. La tendencia en los casos de hombres es de 25 años de edad (con otro pico de tedencia menos pronunciado a los 45 años). La tendencia es mujeres es de casos de mujeres con 29 años de edad.

      Respecto a la potencia del motor, en hombres se observa un distribución uniforme en diferentes tipos de potencia del motor (mediana de 110 unidades de potencia), mientras que en mujeres se observa una tendencia pronunciada a motores de 90 unidades de fuerza.

      Comparando la varibale Acc con las otras variables, podemos observar que:

      1. Son más los hombres, que mujeres, los que se accidentan.

      2. La mediana de la experiencia de los hombres que se accidentaron es de 5 años, mientras de que los casos de hombres que no se accidentaron, la medianda de la experiencia es de 15 años. La tendencia de la experiencia de las mujeres accidentadas y no accidentadas es entre 9 y 10 años de experiencia. Aunque este indicador en este caso puede estar influenciado por la distribución de años de experiencia en mujeres del conjunto, se puede observar que las mujeres no accidentadas tienen ligeramente más edad que las mujeres accidentadas.

      3. La edad, dada su distribución univaridad, influencia los promedios de edad tanto para los casos de accidente, como para los casos de no accidente, Aunque se puede evidenciar que en promedio la edad de los no accidentados (tanto para hombres como para mujeres) es mayor que la edad de las personas que sufrieron accidentes.

      4. Respecto a la potencia, se puede observar que para los casos de accidente, los vehiculos tienden a tener más potencia, comparado con la potencia promedio de los vehiculos de los casos en que no hubo accidentes.

      Respecto al analisis bivariado del resto de variables.

      1. Se puede una relación positiva, no muy fuerte entre la edad y la experiencia. Esto puede debe ser efecto de la variabilidad de la edad en que las personas deciden aprender a conducir un vehiculo. Por ejemplo en mi caso, con 37 años, aún no he conducido un auto. Sin embargo, una vez se inicia a conducir, a medidad que aumente la edad, también aumentan los años de experiencia.

      2. Se puede observar una relación negativa, pero debil, entre la experiencia y la potencia del auto. Es decir, se observa que a medida que los conductores que tienen más experiencia, tienden a usar vehiculos con menos potencia (al menos en este conjunto de datos).

      3. Respecto a la edad y la potencia del motor, no se logra observar una relación relevante entre estas dos variables.

    2. Utilice la función glm, del software R, para ajustar los siguientes modelos de regresión logística:

      Modelo 1: Acc ~ Exp Modelo 2: Acc ~ Exp + genero

      Represente gráficamente el ajuste de los 2 modelos (observados vs predichos).

      Se realiza el ajuste del modelo clasificación para una sola variable predictoria, Exp, y se obtiene los siguientes resultados:

      Se puede observar que a medidad que aumenta la experiencia, disminuye la probabilidad de tener un accidente. Para el modelo_2, se usa adicionalmente la variable \(EsMujer = (Sexo == mujer)\), para modelar el genero. El gráfico a continuación muestra las observaciones y las probabilidades del modelo 2.

      Aqui se puede observar un comportamiento similar, tanto para hombres como para mujeres, donde a medida de que el conductor tiene más experiencia, hay menos probabilidad de tener un accidente. De la gráfica también podemos observar que la distribución de probabilidad de accidentes para hombres se mantiene uniformemente mayor para hombres que para mujeres, en todos los rangos de experiencia observados.

    3. Escriba las ecuaciones de pronóstico asociadas a los 2 modelos.

      De acuerdo a los coeficientes \(\beta\) del modelo 1, mostrado a continuación:

      ## (Intercept)         Exp 
      ##    1.941925   -0.245607

      La ecuación de pronóstico es:

      \[ P(Acc=1 | X = x) = {1 \over {1 + e^{-(1.94 - 0.25*Exp)}}} \]

      Mientras que dados los coeficientes del modelo 2:

      ## (Intercept)         Exp   EsMujersi 
      ##   3.8756040  -0.2399985  -2.9865699

      La ecuación pronostico del modelo 2 es:

      \[ P(Acc=1 | X = x) = {1 \over {1 + e^{-(3.88 - 0.24*Exp + 2.99 * Sexo)}}} \]

    4. A través de indicadores de bondad de ajuste (incluyendo Deviance, AIC, la curva ROC, el AUC y los test de razón de verosimilitud correspondientes), evalúe y compare el ajuste de los 2 modelos anteriores.

      Se calcula los indicadores de bondad para tanto para el modelo 1 como para el modelo 2. Se obtienen los siguientes resultados:

      ##               AIC deviance        R2       AUC
      ## Modelo 1 44.00583 40.00583 0.1631205 0.7983333
      ## Modelo 2 35.24875 29.24875 0.3881471 0.8683333

      El indicador AIC (log de verosimilitud penalizada por parametros) es menor para el modelo 2, indicando que el modelo 2 hace mejores predicciones que el modelo 1. Por otra parte la devianza, que es un indicador de distancia desde el modelo al modelo perfecto, aporta más evidencias de que el modelo 2 se comporta mejor. El \(R^2\), que es la reducción porcentual de la logverosimilitud, es más grande en el modelo 2, indicando que este tiene más poder de predicción. Finalmente el AUC, que es el area entre la curva ROC y el modelo aleatorio, es también mayor en el modelo 2, indicando que usando este modelo se podría obtener un mejor clasificador.

      Se realiza el test de verosimilitud del modelo 2 respecto del modelo 1:

      ## Analysis of Deviance Table
      ## 
      ## Model 1: Acc ~ Exp
      ## Model 2: Acc ~ Exp + EsMujer
      ##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
      ## 1        33     40.006                        
      ## 2        32     29.249  1   10.757 0.001039 **
      ## ---
      ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

      Se puede observar que el modelo 2 mejora la devianza en 10.76 unidades y que este indicador es estadisticamente muy relevante. Se realiza el test de verosimilitud de los parametros del modelo 2:

      ## Analysis of Deviance Table
      ## 
      ## Model: binomial, link: logit
      ## 
      ## Response: Acc
      ## 
      ## Terms added sequentially (first to last)
      ## 
      ## 
      ##         Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
      ## NULL                       34     47.804            
      ## Exp      1   7.7977        33     40.006 0.005231 **
      ## EsMujer  1  10.7571        32     29.249 0.001039 **
      ## ---
      ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

      Se puede observar que tanto la varible Exp como la varible Sexo ayudan a acercar el modelo 2 al modelo perfecto y que esta estimación es relevante estadisticamente.

    5. Seleccione el mejor de los modelos anteriores, interprete los coeficientes estimados y valide su significancia.

      Se selecciona entonces el modelo 2 como el mejor. Verificamos el resumen de los datos del modelo 2:

      ## 
      ## Call:
      ## glm(formula = Acc ~ Exp + EsMujer, family = "binomial", data = accidentes)
      ## 
      ## Coefficients:
      ##             Estimate Std. Error z value Pr(>|z|)   
      ## (Intercept)   3.8756     1.6936   2.288  0.02211 * 
      ## Exp          -0.2400     0.1176  -2.040  0.04131 * 
      ## EsMujersi    -2.9866     1.0683  -2.796  0.00518 **
      ## ---
      ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      ## 
      ## (Dispersion parameter for binomial family taken to be 1)
      ## 
      ##     Null deviance: 47.804  on 34  degrees of freedom
      ## Residual deviance: 29.249  on 32  degrees of freedom
      ## AIC: 35.249
      ## 
      ## Number of Fisher Scoring iterations: 5

      Se puede observar que todos los coficientes \(\beta\) tienen relevancia estadistica. Obervamos de nuevo los coeficientes del modelo 2 y los intervalos de confianza:

      ##                   Coef      2.5 %      97.5 %
      ## (Intercept)  3.8756040  1.2646915  8.36401717
      ## Exp         -0.2399985 -0.5285605 -0.04017264
      ## EsMujersi   -2.9865699 -5.5884081 -1.12921263

      Pero para la interpretación debemos calcular los odds de este modelo (el cociente de probabilidades a favor y en contra de que ocurra un accidente). A continuación el resultado.

      ##                  e-beta       2.5 %       97.5 %
      ## (Intercept) 48.21180847 3.541999825 4289.8934217
      ## Exp          0.78662902 0.589452891    0.9606236
      ## EsMujersi    0.05046023 0.003740979    0.3232877
      • \(e^{\beta_0} = 48.21 > 1\) indica que cuando la experiencia y EsMujer tienden a 0 (En el caso de es mujer, el valor de referencia es ser hombre), entonces hay mayor probabilidad de tener un accidente.

      • \(e^{\beta_1}=0.79 < 1\) indica que entre más experiencia, menos probabilidad de accidente

      • y \(e^{\beta_2}=0.05 < 1\), indica que cuándo la obsevación es de una mujer, entonces hay menos probabilidad de que haya habido un accidente.

      • Para los 3 indicadores odds, se puede observar que los intervalos de confianza no contienen el 1, lo cual facilita su interpretación.

    6. Para el modelo seleccionado en el punto v. evalúe los indicadores de bondad de clasificación (luego de identificar el mejor punto de corte).

      Verificamos los gráficos de densidad del ajuste de probabilidad de accidente, identificando su comportamiento para los casos en que si hubo accidente y en los que no hubo accidente.

      Se puede observar que un punto de corte entre 0.4 y 0.5 optimizaría el acierto en la mayoría de los casos. Se verifica en la curva ROC cuál sería el mejor punto de corte:

      En este caso, podemos observar que el algoritmo detecto un punto de corte diferente para la probabilidad de tener un accidente en \(0.845\). Esto puede ocurrir por el desbalance que hay de casos de no accidente, por lo cual el modelo favorece optimizar la especificidad, máximixando así la detección de casos de no accidentes. Dependiendo de los objetivos del negocio se podría dar más prioridad a la detección de casos de accidente, por ejemplo, podría no ser tan impactante para mi negocio si obtengo un falso positivo, es decir, predigo un accidente y entonces tomo medidas para evitarlo (aunque en la realidad no hubiese ocurrido), pero por el contrario si sea muy impactante para el negocio obtener falsos negativos, es decir, dejar pasar casos clasificados como no accidente y en la realidad comprobamos que si lo hubo.

      En todo caso continuamos con el punto de corte en \(0.845\) para el ejercicio y a continuación calculamos entonces los indicadores de bondad de clasificación para el modelo 2.

      Matriz de confusión:

      ##           Reference
      ## Prediction no si
      ##         no 20  6
      ##         si  0  9
      • De aquí podemos observar que el modelo predijo correctamente 20 casos de no accidente (verdaderos negativos)

      • También podemos observar que el modelo predijo correctamente 9 casos de accidente (verdaderos positivos)

      • El modelo no cometió errores prediciendo accidentes donde no los había, es decir, 0 falsos positivos.

      • y finalmente el modelo predijo incorrectamente que 6 casos eran de no accidente, pero realmente si eran casos de accidente (falsos negativos).

      Como era de esperarse, el modelo favorecio la disminución de falsos positivos. A continuación calculamos los indicadores generales de bondad de clasificación:

      ##                        RL2
      ## Accuracy       0.828571429
      ## Kappa          0.631578947
      ## AccuracyLower  0.663501700
      ## AccuracyUpper  0.934378199
      ## AccuracyNull   0.571428571
      ## AccuracyPValue 0.001201771
      ## McnemarPValue  0.041226833

      Se observan los siguiente indicadores:

      • La exactitud del modelo (accuracy) es del 82.86%, es decir, el modelo predijo correctamente el 82.86% de los casos, lo cual indica un buen nivel de rendimiento general para clasificar tanto casos de accidente como de no accidente.

      • El indice Kappa es del 0.63, que deacuerdo con las convenciones establecidas (vease cohens-kappa-statistic), un indice Kappa entre el 0.61 y el 0.80 denota una concordancia sustancial, es decir, hay una concordancia sustancial entre las clasificaciones del modelo y los valores verdaderos.

      Los siguientes indicadores de bondad verifican el rendimiento de la clasificación, por cada clase.

      ##                            RL2
      ## Sensitivity          0.6000000
      ## Specificity          1.0000000
      ## Pos Pred Value       1.0000000
      ## Neg Pred Value       0.7692308
      ## Precision            1.0000000
      ## Recall               0.6000000
      ## F1                   0.7500000
      ## Prevalence           0.4285714
      ## Detection Rate       0.2571429
      ## Detection Prevalence 0.2571429
      ## Balanced Accuracy    0.8000000

      De aqui podemos observar:

      • La sensibilidad o tasa de verdaderos positivos es del 60%, es decir, el 60% de los casos de accidente real fueron detectados correctamente.

      • La especificidad o tasa de verdaderos negativos es del 100%, es decir, el 100% de los caso no accidente real fueron identificados correctamente.

      • El \(F_\beta\) Score es un indicador que combina la sensibilidad y la especificidad para medir las bondades del modelo. En este caso \(F_1\) Score, se refiere a que se busca un balance entre la presición y la sensibilidad, al calcular el parametro, que se estimo en 75%. Este parametro podría ser util, si utilizamos un \(\beta =0\) para dar mayor prioridad a la reducción de falsos negativos.

    7. Determine si existe una mejora significativa en el modelo seleccionado, cuando se adicionan las variables edad y potencia del motor.

      Se entrenan otros 3 modelos adicionales, uno que adicione la edad, otro que adicione solo la potencia del motor, y un ultimo modelo que adicione ambas varibales. También se calcula los indicadores de bondad de ajuste para cada uno de estos modelos.

      ##               AIC deviance        R2       AUC
      ## Modelo 1 44.00583 40.00583 0.1631205 0.7983333
      ## Modelo 2 35.24875 29.24875 0.3881471 0.8683333
      ## Modelo 3 37.24052 29.24052 0.3883194 0.8600000
      ## Modelo 4 20.78321 12.78321 0.7325888 0.9650000
      ## Modelo 5 22.69953 12.69953 0.7343392 0.9666667

      Podemos observar que hay una clara mejora con los modelo 4 y 5, respecto del modelo 2. El modelo 3, por el contrario, empeora la bondad de ajuste, por lo cual se descarta este modelo. Se realizan entonces los test de verosimilitud de pasar del modelo 2 al modelo 4:

      ## Analysis of Deviance Table
      ## 
      ## Model 1: Acc ~ Exp + EsMujer
      ## Model 2: Acc ~ Exp + EsMujer + Pot
      ##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
      ## 1        32     29.249                          
      ## 2        31     12.783  1   16.465 4.954e-05 ***
      ## ---
      ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

      De aquí podemos observar que el modelo 4 mejora en 16.46 unidades la devianza del modelo, indicando que el modelo se acerca más al modelo perfecto. Al observar el p-value, podemos concluir que este indicador tiene una alta relevancia estadistica.

      Ahora realizamos test de verosimilitud de pasar del modelo 4 al modelo 5.

      ## Analysis of Deviance Table
      ## 
      ## Model 1: Acc ~ Exp + EsMujer + Pot
      ## Model 2: Acc ~ Exp + EsMujer + Pot + Edad
      ##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
      ## 1        31     12.783                     
      ## 2        30     12.700  1 0.083677   0.7724

      Hay una mejora en la devianza del modelo, pero muy pequeña. Además vemos que p-value al estar muy lejos del 0.05, indica que esta hipotesis de mojora en la devianza no es estadisticamente relevante.

      Vemos el test de verosimilitud del modelo 5 respecto del modelo nulo:

      ## Analysis of Deviance Table
      ## 
      ## Model: binomial, link: logit
      ## 
      ## Response: Acc
      ## 
      ## Terms added sequentially (first to last)
      ## 
      ## 
      ##         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
      ## NULL                       34     47.804              
      ## Exp      1   7.7977        33     40.006  0.005231 ** 
      ## EsMujer  1  10.7571        32     29.249  0.001039 ** 
      ## Pot      1  16.4655        31     12.783 4.954e-05 ***
      ## Edad     1   0.0837        30     12.700  0.772374    
      ## ---
      ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

      Y podemos comprobar que efectivamente la variable edad no agrega nada de valor al modelo 5. Es posible que esto este relacionado con la correlación debil que encontramos entre la edad y la experiencia.

      Pero si verificamos el resumen de los indicadores del modelo 4:

      ## 
      ## Call:
      ## glm(formula = Acc ~ Exp + EsMujer + Pot, family = "binomial", 
      ##     data = accidentes)
      ## 
      ## Coefficients:
      ##             Estimate Std. Error z value Pr(>|z|)  
      ## (Intercept) -16.8520     9.1690  -1.838   0.0661 .
      ## Exp          -0.4800     0.3227  -1.487   0.1369  
      ## EsMujersi    -3.1699     2.2643  -1.400   0.1615  
      ## Pot           0.2390     0.1026   2.330   0.0198 *
      ## ---
      ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      ## 
      ## (Dispersion parameter for binomial family taken to be 1)
      ## 
      ##     Null deviance: 47.804  on 34  degrees of freedom
      ## Residual deviance: 12.783  on 31  degrees of freedom
      ## AIC: 20.783
      ## 
      ## Number of Fisher Scoring iterations: 8

      Podemos observar que en realidad el modelo tiene unos p-value muy elevados para cada variable predictora, reflejando entonces una baja relevancia estadistica de los componentes del modelo.

      Se verifica en todo caso los odds del modelo 4:

      ##                   e-beta        2.5 %     97.5 %
      ## (Intercept) 4.800287e-08 1.837084e-18 0.04136065
      ## Exp         6.188137e-01 2.034211e-01 0.93702368
      ## EsMujersi   4.200812e-02 4.837072e-05 1.33970799
      ## Pot         1.269987e+00 1.092340e+00 1.66855510
      • \(e^{\beta_0} = 4.8*10^{-8} < 1\) No tiene sentido, dice que cuando todas las variables tienen a 0, entonces hay menor probabilidad de accidente.

      • \(e^{\beta_2}=0.6 < 1\) indica que entre más experiencia, menos probabilidad de accidente hay, lo cual tienen sentido.

      • \(e^{\beta_2}=0.042 < 1\), indica que cuándo la obsevación es de una mujer, entonces hay menos probabilidad de que haya habido un accidente. El intervalo de confianza contiene el 1.

      • y \(e^{\beta_3}=1.26 > 1\), indica que a mayor potencia del moto, mayor probabilidad de accidente.

    8. Haciendo uso de sus habilidades de modelación, genera un breve reporte de sus hallazgos en el cual oriente a la compañía sobre los factores que afectan la siniestralidad.

      De acuerdo al estudio realizado, se encontro que el mejor modelo de regresión logistica para clasificar la accidentalidad es aquel que uso como variables predictoras: la experiencia del conductor y saber si el conductor es mujer. En particular, los factores que aumentan el riesgo de accidentalidad es tener pocos años de experiencia y ser hombre.

      Por otro lado, los resultados del analisis no dieron una evidencia concluyente respecto a la influencia de la potencia del motor en la probabilidad de accidente. Pero si se decidiera incluir esta variable, su influencia sería, a medida que el motor tiene más potencia, entonces hay más probabilidad de accidente.