Challenge # 2
En esta práctica, vamos a validar un modelo de Regresión Lineal Múltiple (RLM) mediante una función. Utilizaremos un modelo ajustado que fue desarrollado como ejemplo durante la sesión #3 del módulo de Análitica Avanzada de Datos.
Es importante destacar que el contenido de este documento, hasta el apartado Modelo Ajustado, fue extraído del documento session-03-rlm.qmd con autoría de Jorge I. Vélez, PhD. Esto se hizo con el propósito de implementar la función de validación de supuestos.
The Tire Company (TTC) fabrica compuestos para llantas de F1. En la actualidad, TTC desarrolla un nuevo compuesto para la temporada 2023, con miras a superar nuevamente a su más desafiante competidor.
En el proceso de fabricación existen 3 compuestos \(x_1, x_2\) y \(x_3\) que conforman cada llanta. De acuerdo con las regulaciones, una llanta de F1 debe pesar no más de 13 kilogramos sin incluir el rin y los sensores. Además, se sabe que el costo por cada gramo de material \(j\), en dólares, es \(c_j = \sqrt{j}\), \(j=1,2,3\). TTC está convencida de que modificando la composición de la llanta es posible alcanzar un mejor rendimiento (medido en kilómetros recorridos), de tal forma que no es necesario realizar tantas paradas en pits durante una competencia regular. En la fase de experimentación, el equipo de analítica, donde usted trabaja, encontró que la formulación actual permite recorrer \(26 \pm 7\) kilómetros con un set de llantas. Se sabe que, en libras, \(x_1\in(5, 12)\), \(x_2\in(3, 20)\) y \(x_3\in(1, 6)\). Por lo tanto, el resto del peso de la llanta debe completarse con caucho sin procesar.
El día de ayer, usted fue asignado a este proyecto con miras a optimizar la composición de las llantas en TTC, de tal manera que el número de kilómetros recorridos por cada set de llantas sea considerablemente mayor que en la actualidad.
Las primeras 6 filas de la base de datos d son
## primeras 6 líneas
head(d) y x1 x2 x3
1 23.72364 9.574030 10 2.451304
2 30.99403 9.685377 7 2.333709
3 17.03163 6.430698 10 1.181157
4 26.08888 9.152238 7 1.982518
5 30.49450 8.208728 9 3.635287
6 23.71911 7.595480 10 3.791815Analicemos inicialmente la correlación entre las variables disponibles:
## matriz de correlación
par(mfrow = c(1,1), mar = c(.1, .1, .1, .1))
qgraph(cor(d), graph = "cor", layout = "spring",
sampleSize = nrow(d),
legend.cex = 1, alpha = 0.05)
Numéricamente, la matriz de correlación es
## matriz de correlación
cor(d) y x1 x2 x3
y 1.0000000 0.29127317 -0.762485290 0.288541500
x1 0.2912732 1.00000000 0.040792570 0.049891558
x2 -0.7624853 0.04079257 1.000000000 -0.003470425
x3 0.2885415 0.04989156 -0.003470425 1.000000000Estos resultados indican que la correlación entre \(y\) y \(x_1\) es 0.2912, entre \(y\) y \(x_2\) es -0.762, y entre \(y\) y \(x_3\) es 0.288. En cuanto a la correlaciones entre las variables independientes, podemos concluir que estas son pequeñas, lo cual sugiere que, efectivamente, \(x_1, x_2\) y \(x_3\) son independientes.
El gráfico 3D entre \(x_1\), \(x_2\) y \(y\) sería:
attach(d)
p <- plot_ly(x = ~x1, y = ~x2, z = ~y, type = "scatter3d",
mode = "markers",
marker = list(size = 8, color = toRGB("blue"),
opacity = 0.8,
line = list(width = 1,
color = toRGB("black")))) %>%
layout(title = 'Distancia recorrida (km) vs (cantidad de compuesto (g) de x1,x2)',
scene = list(xaxis = list(title = 'Compuesto x1 (g)'),
yaxis = list(title = 'Compuesto x2 (g)'),
zaxis = list(title = 'Distancia Recorrida (km)')))
pEl modelo ajustado es:
fit <- lm(y~., data = d)
summary(fit)
La ecuación del modelo ajustado es
\[ \hat{y} = 25.24 + 1.193 x_1 - 1.651x_2 + 1.735x_3 \]
summary(fit)
Call:
lm(formula = y ~ ., data = d)
Residuals:
Min 1Q Median 3Q Max
-7.7744 -2.2090 -0.0784 1.8554 7.5166
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 25.24199 1.48841 16.959 < 2e-16 ***
x1 1.19280 0.15717 7.589 3.49e-12 ***
x2 -1.65137 0.08684 -19.016 < 2e-16 ***
x3 1.73492 0.26130 6.640 5.78e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.815 on 146 degrees of freedom
Multiple R-squared: 0.7584, Adjusted R-squared: 0.7535
F-statistic: 152.8 on 3 and 146 DF, p-value: < 2.2e-16summary(fit)$fstatistic[4]<NA>
NA # Llamar a los grados de libertad
F_statistic <- summary(fit)$fstatistic[1]
DF_numerador <- summary(fit)$fstatistic[2]
DF_denominador <- summary(fit)$fstatistic[3]
1 - pf(F_statistic, DF_numerador, DF_denominador)value
0 Considerese el modelo de regresión lineal múltiple dado anteriormente, cuya ecuación general se encuentra dada por:
\[ y_i = \beta_0 + \beta_1X_{1i} + \beta_2X_{2i}+, \dots , \beta_kx{ki} + \varepsilon_i, \space \space i = 1,2, \dots, n \]
Donde: \(\varepsilon_i \sim \mathcal{N}(0,\sigma^2), \space \space i = 1,2,\dots,n\).
Para la correcta implementación de un modelo de RLM, es fundamental asegurar que se cumplan los supuestos relacionados con los residuos del modelo (\(\varepsilon_i\)).
Normalidad
Independencia
Varianza Constante
A continuación se define una función para automatizar la validación de los supuestos antes mencionados y la validación del modelo:
## Defición de función:
val_model <- function(model) {
#Definición de datos importantes
ei <- residuals(model) # Extraer los residuales ordinarios
ri <- rstudent(model) # Extraer residuos estudentizados
r_i <- which(ri > 3 | ri < -3) # Almacenar residuos mayores a 3 o menores a -3
r2 <- summary(model)$adj.r.squared #Coeficiente de determinación
F_statistic <- summary(model)$fstatistic[1]
DF_num <- summary(model)$fstatistic[2]
DF_den <- summary(model)$fstatistic[3]
f_sta <- 1 - pf(F_statistic, DF_num, DF_den)# FStadistico
#Coeficiente de determinación R^2
coef_d <- ifelse(r2 >= 0.70,
paste('Las variables predictoras explican en un', round(r2 * 100, 2), '% la variabilidad de la respuesta'),
'El modelo tiene un desempeño deficiente')
# Prueba de significancia Global
s_global <- ifelse(f_sta<= 0.05, 'El modelo tiene un efecto globalmente significativo','El modelo no tiene un efecto globalmente significativo' )
#Prueba de significancia Marginal
var_predictoras <- names(coef(model))[-1] # Obtener los nombres de las variables predictoras, excluyendo el intercepto
# Verificar individualmente los p-valores para cada variable predictora
resultado <- sapply(var_predictoras, function(var_nombre) {
p_val <- summary(model)$coefficients[var_nombre, "Pr(>|t|)"]
if (!is.na(p_val) && p_val <= 0.05) {
return(TRUE)
} else {
return(FALSE)
}
})
# Verificar si todas las variables cumplen con el criterio
if (all(resultado)) {
s_marginal <- "Controlando las variables"
s_marginal <- paste(s_marginal, paste(var_predictoras[resultado], collapse = ","))
s_marginal <- paste(s_marginal, " en el proceso, permitiría modificar la respuesta de la variable dependiente")
} else {
s_marginal <- "Al menos una de las variables predictoras no cumple con el criterio"
}
#Supuesto de Normalidad
normalidad <- ifelse(shapiro.test(ri)$p.value >= 0.05,
'Los residuales siguen una distribución normal',
'Los residuales no siguen una distribución normal')
#Supuesto de Independencia
independencia <- ifelse(durbinWatsonTest(model)$p >= 0.05,
'Los residuales son independientes',
'Los residuales no son independientes')
#Supuesto de Homocedasticidad
homocedasticidad <- ifelse(ncvTest(model)$p >= 0.05,
'Los residuales tienen varianza constante',
'Los residuales no tienen varianza constante')
#Detección de Outliers
n_outliers <- ifelse(length(r_i) == 0, 0, length(r_i))
OutLier <- ifelse(n_outliers >= 0, 'No existen Outliers',
paste('Existen', n_outliers, 'Outliers'))
return(list(coef_d = coef_d,
s_global = s_global,
s_marginal = s_marginal,
normalidad = normalidad,
independencia = independencia,
homocedasticidad = homocedasticidad,
OutLier= OutLier))
}df_val_model <- data.frame(val_model(fit))
names(df_val_model) <- c('Coeficiente de Determinación', 'Significancia Global', 'Significancia Marginal', 'Test de Normalidad' , 'Test de Independencia', 'Test de Homocedasticidad', 'Verificación de Outliers')
# Generar la tabla con kableExtra
kbl(df_val_model,
align = "ccccc",
position = "top") %>%
kable_classic(full_width = FALSE,
html_font = "Verdana",
bootstrap_options = "hover")| Coeficiente de Determinación | Significancia Global | Significancia Marginal | Test de Normalidad | Test de Independencia | Test de Homocedasticidad | Verificación de Outliers | |
|---|---|---|---|---|---|---|---|
| value | Las variables predictoras explican en un 75.35 % la variabilidad de la respuesta | El modelo tiene un efecto globalmente significativo | Controlando las variables x1,x2,x3 en el proceso, permitiría modificar la respuesta de la variable dependiente | Los residuales siguen una distribución normal | Los residuales son independientes | Los residuales tienen varianza constante | No existen Outliers |
Conclusión : Dando seguimiento a los resultados previos podemos confirmar la validez de nuestro modelo para realizar predicciones, debido a que se cumplen con los supuestos de normalidad, independencia y homocedasticidad de los residuales del modelo.
Para más detalles sobre el modelo de RLM se sugiere consultar el Capítulo 3 del texto Modelos de Regresión: Una aproximación práctica con R.
Montgomery, D. C., Peck, E. A.„ Vining, G. G. (2012). Introduction to Linear Regression Analysis (5th ed.). Wiley & Sons