LABORATORIO 4 – REGRESIÓN LINEAL

Punto 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:

Analisis

El ejercicio busca evaluar la relación entre la carga del sistema y el tiempo de respuesta de dos tipos de discos duros:

  • disk SDD (Conf = 0)
  • disk HDD (Conf = 1)

Las variables son:

  • Load: Número de consultas por minuto, representa la carga del sistema.
  • Time: Tiempo de respuesta del disco (en segundos), la variable objetivo.
  • Config: Variable categórica (Dummy) que distingue entre SDD (0) y HDD (1).

Cargando los datos

#load tiene el valor de la variable carga
load <- c(1, 2, 2.4, 3.1, 4, 4.3, 5.8, 6.6, 7.5, 8, 9, 9.2, 10.2, 
          1.8, 2, 2.5, 3.9, 4.2, 5.5, 6.4, 7, 8, 8.2, 9.1, 9.5)

#Tiempo de carga
time <- c(0.9, 0.3, 2, 0.8, 2.7, 2.6, 2.5, 3.2, 3.7, 3.9, 5.3, 4.2, 3.9,
          1.1, 1.5, 0.5, 1.5, 1.6, 3.3, 3.3, 3.5, 4.3, 4, 4.3, 5.8)

#Creando vector de 1 y 0 que asemejan el tipo de disco duro
config <- c(1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 
          1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0) # Configuración del disco

#Dataframe para generar la tabla
disks <- data_frame(Config = config, Load = load, Time = time)

Representación gráfica inicial

Se busca representar gráficamente cómo varía el tiempo de respuesta con la carga del sistema para los dos tipos de disco.

Propósito:

  • Visualizar si existe una relación lineal entre las variables Carga y Tiempo.
  • Observar si esta relación varía según el tipo de disco.
ggplot(disks, aes(x = Load, y = Time, color = factor (Config)))+
  geom_point()+
  geom_smooth(method = "lm", se = FALSE)#Aplico el método

  labs(
    title = "Load vs Time", color = "Configuration",
    x = "System Load (Query/min)", y = "Response Time (s)" # Variables
  )+
  theme_minimal()
## NULL

Calculando las fuerzas de correlación

#Fuerza de la relación Load vs Time en cada tipo de disco calculando coeficiente de correlación
cor(disks$Load[disks$Config == 0], disks$Time[disks$Config == 0]) #Para los discos SdD
## [1] 0.9938293
cor(disks$Load[disks$Config == 1], disks$Time[disks$Config == 1]) #Para los discos HSD
## [1] 0.9640003

Evaluando el Resultado

  • 0.9938293 para SDD indica una relación positiva fuerta ya que es un valor cercano a 1.
  • 0.9640003 para HDD indica una relación positiva fuerta ya que es un valor cercano a 1.

Ajustando un primer módelo de regresión simple (Modelo1)

Se ajusta un modelo de regresión simple que relaciona el tiempo de respuesta (Time) con la carga (Load), ignorando el tipo de disco.

Modelo 1

# Modelo 1: Regresión simple
modelo1 <- lm(Time ~ load, data = disks)
summary(modelo1)
## 
## Call:
## lm(formula = Time ~ load, data = disks)
## 
## 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    
## load         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

Ecuación del modelo1:

Time = 𝛽0 + 𝛽1 ⋅ Load

  • β0: Es el intercepto que indica el tiempo de respuesta cuando la carga es cero.
  • β1: Es la pendiente que nos indica cómo cambia el tiempo por cada unidad adicional de carga.

Ahora evaluemos el R^2 para la bondad de ajuste. Representa la variabilidad del tiempo por la carga

par(mfrow = c(2, 2))
plot(modelo1)

# Pendiente la interpretación de resultados

Modelo 2

Ahora el modelo incluirá el tipo de disco. Se ajusta un modelo que incluye la interacción entre Load (Carga) y Config (tipo de disco) lopermite que la relación entre carga y tiempo sea distinta para SDD y HDD.

# Modelo 2: Incluir variable Dummy e Interacción
modelo2 <- lm(Time ~ Load * config, data = disks)
summary(modelo2) # Mostrar resumen
## 
## Call:
## lm(formula = Time ~ Load * config, data = disks)
## 
## 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 ***
## Load         0.71979    0.03367  21.376 9.88e-16 ***
## config       2.26391    0.26520   8.536 2.86e-08 ***
## Load:config -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

#Interpretación

Ecuación del modelo2:

Time = 𝛽0 + 𝛽1 ⋅ Load +𝛽2 ⋅ Config + β3⋅(Load ⋅ Config)

Interpretación de los coeficientes

  • β0: Representa el intercepto para el SDD
  • β1: Es la pendiente de la relación entre Load(Carga) y Time(Tiempo).
  • β2: Es la diferencia de la pendiente entre HDD y SDD.

Test ANOVA

Ahora se usa el test ANOVA para evaluar si la inclusión de la variable cualitativa de configuración del disco y la interacción de la carga mejora significativamente el ajuste del modelo.

# Test ANOVA para comparar modelos
anova(modelo1, modelo2) # Leer mas para entenderlo mejor
## Analysis of Variance Table
## 
## Model 1: Time ~ load
## Model 2: Time ~ Load * config
##   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

Resultado:

Si el valor ppp es pequeño (por ejemplo, p<0.05p, el modelo con interacción (Modelo 2) se ajusta mejor.

Gráfico del modelo2 separado por tipo de disco.

par(mfrow = c(1, 1))
ggplot(disks, aes(x = Load, y = Time, color = factor(Config))) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE) +
  labs(title = "Modelo 2: Regression with Interaction",
       x = "System Load (Query/min)", y = "Response Time (s)", #Variables
       color = "Disk type") + #Color por tipo de disco
  theme_minimal()

Evaluación de supuestos

# Evaluación de supuestos
par(mfrow = c(2, 2))
plot(modelo2)

Conclusiones

Según la gráfica se puede notar que en los discos HDD y SDD, el tiempo de respuesta a la acrga del sistema varía segun el tipo de disco. La respuesta de los discos SDD es muy menor en comparación a los discos SDD.

A medida que aumenta la carga del sistema , el tiempo de carga tambien aumenta pero el la pendiente de la relación entre los dos tiempo de carga es mas pornuciada para los disco HDD, lo que indica que estos disco tienden a ser mas sensibles a estos cambios.

La inclusión de la interacción en el modelo permite mejorar el ajuste, ya que al dejar que las pendientes y los interceptos varíen entre los dos tipos de discos. Esto es consistente de acuerdo a las diferencias observadas en los datos y en la representación gráfica.

Punto 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).
#Cargando el archivo con los datos
accidents <<- readxl::read_excel("accidentes.xlsx")
## Resumen
df <- accidents
summary(df)
##       Acc              Exp              Edad         Pot             Sexo    
##  Min.   :0.0000   Min.   : 1.000   Min.   :20   Min.   : 70.0   Min.   :1.0  
##  1st Qu.:0.0000   1st Qu.: 6.500   1st Qu.:25   1st Qu.: 90.0   1st Qu.:1.0  
##  Median :0.0000   Median : 9.000   Median :29   Median : 95.0   Median :1.0  
##  Mean   :0.4286   Mean   : 9.543   Mean   :31   Mean   :101.6   Mean   :1.4  
##  3rd Qu.:1.0000   3rd Qu.:12.000   3rd Qu.:36   3rd Qu.:110.0   3rd Qu.:2.0  
##  Max.   :1.0000   Max.   :20.000   Max.   :56   Max.   :150.0   Max.   :2.0
head(accidents)
## # A tibble: 6 × 5
##     Acc   Exp  Edad   Pot  Sexo
##   <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     0    10    30    90     1
## 2     0    15    40    85     1
## 3     0     7    25    95     1
## 4     1     1    21   145     2
## 5     0    10    29    70     1
## 6     1     2    20   120     2

Ahora, usando herramientas de análisis exploratorio, estudie la asociación entre la siniestralidad y el conjunto de variables predictoras (Edad, Experiencia, Potencia del motor y Sexo).

boxplot(Exp ~ Acc, data = df, main = "Experiencia y Accidentes")

boxplot(Edad ~ Acc, data = df, main = "Edad y Accidentes")

boxplot(Pot ~ Acc, data = df, main = "Potencia del Motor y Accidentes")

Usando la función gml del software R

Función glm para los siguinetes modelos de regresión logística

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

Modelo 1

mod_log_1 <- glm(Acc ~ Exp, family = binomial, data = df)
summary(mod_log_1)
## 
## Call:
## glm(formula = Acc ~ Exp, family = binomial, data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   1.9419     0.9816   1.978   0.0479 *
## Exp          -0.2456     0.1044  -2.354   0.0186 *
## ---
## 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: 40.006  on 33  degrees of freedom
## AIC: 44.006
## 
## Number of Fisher Scoring iterations: 4
pred1 <- ifelse(fitted(mod_log_1) > 0.5, 1, 0)
ggplot(data = df, aes(x = Acc, y = fitted(mod_log_1))) +
  geom_point(alpha = 0.5) +
  geom_jitter(width = 0.2, height = 0.01, color = "green") +
  labs(title = "Modelo 1: Observados vs Predichos",
       x = "Valores Observados (Acc)",
       y = "Valores Predichos (Probabilidad)") +
  theme_minimal()

Ecuación de pronóstico modelo 1

La ecuación general para el modelo 1 es:

\[ \log\left(\frac{P(\text{Acc} = 1)}{1 - P(\text{Acc} = 1)}\right) = \beta_0 + \beta_1 \cdot \text{Exp} \]

Modelo 2

mod_log_2 <- glm(Acc ~ Exp + Sexo, family = binomial, data = df)
summary(mod_log_2)
## 
## Call:
## glm(formula = Acc ~ Exp + Sexo, family = binomial, data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -2.0975     1.5546  -1.349  0.17725   
## Exp          -0.2400     0.1176  -2.040  0.04131 * 
## Sexo          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
pred2 <- ifelse(fitted(mod_log_2) > 0.5, 1, 0)
ggplot(data = df, aes(x = Acc, y = fitted(mod_log_2))) +
  geom_point(alpha = 0.5) +
  geom_jitter(width = 0.2, height = 0.01, color = "blue") +
  labs(title = "Modelo 2: Observados vs Predichos",
       x = "Valores Observados (Acc)",
       y = "Valores Predichos (Probabilidad)") +
  theme_minimal()

Ecuación de pronóstico modelo 2

La ecuación general para el modelo 2 es:

\[ \log\left(\frac{P(\text{Acc} = 1)}{1 - P(\text{Acc} = 1)}\right) = \beta_0 + \beta_1 \cdot \text{Exp} + \beta_2 \cdot \text{Sexo} \]

Bonda de ajuste

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), se evalúa y se compara el ajuste de los 2 modelos anteriores.

Modelo AIC

AIC(mod_log_1,mod_log_2)
##           df      AIC
## mod_log_1  2 44.00583
## mod_log_2  3 35.24875

Para el indicador AIC, el modelo 2 tiene un menor AIC que el modelo 1

Curva ROC

roc_mod1 <- roc(df$Acc, fitted(mod_log_1))
roc_mod2 <- roc(df$Acc, fitted(mod_log_2))

ANOVA

anova(mod_log_1, mod_log_2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: Acc ~ Exp
## Model 2: Acc ~ Exp + Sexo
##   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

Luego de realizar el test de verosimilitud con el test “Chisq”, se detecta que incluyendo la variable sexo mejora el p-valor por lo que el modelo 2 es significativamente mejor.

AUC

auc (roc_mod1)
## Area under the curve: 0.7983
auc (roc_mod2)
## Area under the curve: 0.8683

El AUC del modelo 2 es mayor que el del modelo, lo que me indica que tiene una mejor capacidad predictiva

Visualización de la curva ROC

plot(roc_mod1, col = "blue", main = "Curvas ROC de los Modelos")
lines(roc_mod2, col = "red")
legend("bottomright", legend = c("Modelo 1", "Modelo 2"),
       col = c("blue", "red"), lty = 1)

Seleccionando el mejor módelo

De acuerdo al ejercicio anterior, con los indicadores ANOVA y AIC, además del area bajo la curva, selecciono le modelo 2.

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).

Identifiocando mejor punto de corte

# Identificar el mejor punto de corte
mod_selected <- roc(df$Acc, fitted(mod_log_2))
optimal_cutoff <- coords(mod_selected, "best", ret = "threshold")
print(optimal_cutoff)
##   threshold
## 1 0.8449838
# Matriz de confusión: No pude continuar. Revisar documentación por error
# ==  Error en table(Predicted = predicted, Actual = df_used$Acc): 
  #== all arguments must have the same length

#predict_mod <- ifelse(fitted(mod_log_2) > optimal_cutoff, 1, 0)
#conf_matrix <- table(Predict_mod = predict_mod, Actual = df$Acc)
#print(conf_matrix)