Selección de variables
2025-05-11
library(kableExtra)
library(dplyr)
library(PerformanceAnalytics)
library(corrplot)
library(knitr)
library(tidyr)
library(MASS) #para el stepAIC
library(olsrr) # Para ols_mallows_cp
library(MPV)
library(car)CARGA DE DATOS
Primero, procedemos a cargar los datos que utilizaremos para el problema planteado:
data <- data.frame(
Calidad = c(9.8, 12.6, 11.9, 11.1, 13.3, 12.8, 12.8, 12.0, 13.6, 13.9, 14.4, 12.3, 16.1, 16.1, 15.5, 15.5, 13.8, 13.8, 11.3, 7.9, 15.1, 13.5, 10.8, 9.5, 12.7, 11.6, 11.7, 11.9, 10.8, 8.5, 10.7, 9.1, 12.1, 14.9, 13.5, 12.2, 10.3, 13.2),
Claridad = c(1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.5, 0.8, 0.7, 1.0, 0.9, 1.0, 1.0, 1.0, 0.9, 0.9, 1.0, 0.7, 0.7, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.8, 1.0, 1.0, 0.8, 0.8, 0.8, 0.8),
Aroma = c(3.3, 4.4, 3.9, 3.9, 5.6, 4.6, 4.8, 5.3, 4.3, 4.3, 5.1, 3.3, 5.9, 7.7, 7.1, 5.5, 6.3, 5.0, 4.6, 3.4, 6.4, 5.5, 4.7, 4.1, 6.0, 4.3, 3.9, 5.1, 3.9, 4.5, 5.2, 4.2, 3.3, 6.8, 5.0, 3.5, 4.3, 5.2),
Cuerpo = c(2.8, 4.9, 5.3, 2.6, 5.1, 4.7, 4.8, 4.5, 4.3, 3.9, 4.3, 5.4, 5.7, 6.6, 4.4, 5.6, 5.4, 5.5, 4.1, 5.0, 5.4, 5.3, 4.1, 4.0, 5.4, 4.6, 4.0, 4.9, 4.4, 3.7, 4.3, 3.8, 3.5, 5.0, 5.7, 4.7, 5.5, 4.8),
Sabor = c(3.1, 3.5, 4.8, 3.1, 5.5, 5.0, 4.8, 4.3, 3.9, 4.7, 4.5, 4.3, 7.0, 6.7, 5.8, 5.6, 4.8, 5.5, 4.3, 3.4, 6.6, 5.3, 5.0, 4.1, 5.7, 4.7, 5.1, 5.0, 5.0, 2.9, 5.0, 3.0, 4.3, 6.0, 5.5, 4.2, 3.5, 5.7),
Roble = c(4.1, 3.9, 4.7, 3.6, 5.1, 4.1, 3.3, 5.2, 2.9, 3.9, 3.6, 3.6, 4.1, 3.7, 4.1, 4.4, 4.6, 4.1, 3.1, 3.4, 4.8, 3.8, 3.7, 4.0, 4.7, 4.9, 5.1, 5.1, 4.4, 3.9, 6.0, 4.7, 4.5, 5.2, 4.8, 3.3, 5.8, 3.5)
)
datos<-data.frame(
Claridad = c(1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.5, 0.8, 0.7, 1.0, 0.9, 1.0, 1.0, 1.0, 0.9, 0.9, 1.0, 0.7, 0.7, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.8, 1.0, 1.0, 0.8, 0.8, 0.8, 0.8),
Aroma = c(3.3, 4.4, 3.9, 3.9, 5.6, 4.6, 4.8, 5.3, 4.3, 4.3, 5.1, 3.3, 5.9, 7.7, 7.1, 5.5, 6.3, 5.0, 4.6, 3.4, 6.4, 5.5, 4.7, 4.1, 6.0, 4.3, 3.9, 5.1, 3.9, 4.5, 5.2, 4.2, 3.3, 6.8, 5.0, 3.5, 4.3, 5.2),
Cuerpo = c(2.8, 4.9, 5.3, 2.6, 5.1, 4.7, 4.8, 4.5, 4.3, 3.9, 4.3, 5.4, 5.7, 6.6, 4.4, 5.6, 5.4, 5.5, 4.1, 5.0, 5.4, 5.3, 4.1, 4.0, 5.4, 4.6, 4.0, 4.9, 4.4, 3.7, 4.3, 3.8, 3.5, 5.0, 5.7, 4.7, 5.5, 4.8),
Sabor = c(3.1, 3.5, 4.8, 3.1, 5.5, 5.0, 4.8, 4.3, 3.9, 4.7, 4.5, 4.3, 7.0, 6.7, 5.8, 5.6, 4.8, 5.5, 4.3, 3.4, 6.6, 5.3, 5.0, 4.1, 5.7, 4.7, 5.1, 5.0, 5.0, 2.9, 5.0, 3.0, 4.3, 6.0, 5.5, 4.2, 3.5, 5.7),
Roble = c(4.1, 3.9, 4.7, 3.6, 5.1, 4.1, 3.3, 5.2, 2.9, 3.9, 3.6, 3.6, 4.1, 3.7, 4.1, 4.4, 4.6, 4.1, 3.1, 3.4, 4.8, 3.8, 3.7, 4.0, 4.7, 4.9, 5.1, 5.1, 4.4, 3.9, 6.0, 4.7, 4.5, 5.2, 4.8, 3.3, 5.8, 3.5)
)Ahora, veamos la clasificación de cada una de las variables involucradas en el ejercicio:
# Crear el dataframe CORREGIDO (eliminé la coma sobrante después del último paréntesis)
operacionalizacion <- data.frame(
Variable = c("Calidad", "Claridad", "Aroma", "Cuerpo", "Sabor", "Roble"),
Definición = c(
"Puntuación global de calidad asignada por expertos",
"Grado de transparencia y ausencia de partículas en el vino",
"Intensidad y complejidad de los olores característicos",
"Sensación de densidad y textura en boca",
"Balance e intensidad de los sabores principales",
"Notas derivadas del envejecimiento en barricas de roble"
),
Naturaleza = c("Cuantitativa continua", "Cuantitativa continua",
"Cuantitativa continua", "Cuantitativa continua",
"Cuantitativa continua", "Cuantitativa continua"),
Nivel_de_Medicion = c("De razón", "De intervalo", "De intervalo",
"De intervalo", "De intervalo", "De intervalo"),
Criterio_Clasificacion = c(
"Escala numérica de 7.9 a 16.1 puntos",
"Escala de 0.5 (baja) a 1.0 (alta)",
"Escala de 3.3 (suave) a 7.7 (intenso)",
"Escala de 2.6 (ligero) a 6.6 (robusto)",
"Escala de 2.9 (simple) a 7.0 (complejo)",
"Escala de 2.9 (sutil) a 6.0 (marcado)" # Corregí "subtíl" a "sutil"
)
)
operacionalizacion %>%
kbl(
caption = "Tabla de Operacionalización de Variables - Estudio de Calidad del Vino",
) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed", "responsive"),
font_size = 14, position = "center") %>%
row_spec(0, bold = TRUE, color = "white", background = "#8B0000") %>%
row_spec(1:nrow(operacionalizacion), background = "#FFF5EE") %>%
footnote(
general = "Datos adaptados del estudio de Montgomery, Peck y Vining (2003)",
general_title = "Fuente:",
footnote_as_chunk = TRUE
)| Variable | Definición | Naturaleza | Nivel_de_Medicion | Criterio_Clasificacion |
|---|---|---|---|---|
| Calidad | Puntuación global de calidad asignada por expertos | Cuantitativa continua | De razón | Escala numérica de 7.9 a 16.1 puntos |
| Claridad | Grado de transparencia y ausencia de partículas en el vino | Cuantitativa continua | De intervalo | Escala de 0.5 (baja) a 1.0 (alta) |
| Aroma | Intensidad y complejidad de los olores característicos | Cuantitativa continua | De intervalo | Escala de 3.3 (suave) a 7.7 (intenso) |
| Cuerpo | Sensación de densidad y textura en boca | Cuantitativa continua | De intervalo | Escala de 2.6 (ligero) a 6.6 (robusto) |
| Sabor | Balance e intensidad de los sabores principales | Cuantitativa continua | De intervalo | Escala de 2.9 (simple) a 7.0 (complejo) |
| Roble | Notas derivadas del envejecimiento en barricas de roble | Cuantitativa continua | De intervalo | Escala de 2.9 (sutil) a 6.0 (marcado) |
| Fuente: Datos adaptados del estudio de Montgomery, Peck y Vining (2003) |
Para evitar confusiones, se define cada variable explicativa con su respectiva representación:
Y: Calidad (variable respuesta)
X1: Claridad
X2: Aroma
X3: Cuerpo
X4: Sabor
X5: Roble
y<-data$Calidad
x1<-data$Claridad
x2<-data$Aroma
x3<-data$Cuerpo
x4<-data$Sabor
x5<-data$RoblePRUEBA NORMALIDAD
Primero, realizaremos una prueba de normalidad para ver que métodos utilizar para en análisis descriptivo:
shapiro.test(y)##
## Shapiro-Wilk normality test
##
## data: y
## W = 0.98414, p-value = 0.8559shapiro.test(x1)##
## Shapiro-Wilk normality test
##
## data: x1
## W = 0.67094, p-value = 5.805e-08shapiro.test(x2)##
## Shapiro-Wilk normality test
##
## data: x2
## W = 0.95374, p-value = 0.1183shapiro.test(x3)##
## Shapiro-Wilk normality test
##
## data: x3
## W = 0.97357, p-value = 0.4965shapiro.test(x4)##
## Shapiro-Wilk normality test
##
## data: x4
## W = 0.97471, p-value = 0.5333shapiro.test(x5)##
## Shapiro-Wilk normality test
##
## data: x5
## W = 0.97495, p-value = 0.5411Como p-valor de la variable claridad es menor a 0.05, rechazamos la hipótesis nula, es decir, los datos no siguen una distribución normal. Por lo cual, procederemos a utilizar el método de “Spearman”.
DESCRIPTIVO
A continuación procederemos a hacer un análisis descriptivo de los datos, analizando correlación.
tabla_estadisticos <- operacionalizacion %>%
summarise(
Calidad_Min = min(data$Calidad, na.rm = TRUE),
Calidad_Q1 = quantile(data$Calidad, 0.25, na.rm = TRUE),
Calidad_Mediana = median(data$Calidad, na.rm = TRUE),
Calidad_Media = mean(data$Calidad, na.rm = TRUE),
Calidad_Q3 = quantile(data$Calidad, 0.75, na.rm = TRUE),
Calidad_Max = max(data$Calidad, na.rm = TRUE),
Calidad_Sd = sd(data$Calidad, na.rm = TRUE),
Calidad_IQR = IQR(data$Calidad, na.rm = TRUE),
Claridad_Min = min(data$Claridad, na.rm = TRUE),
Claridad_Q1 = quantile(data$Claridad, 0.25, na.rm = TRUE),
Claridad_Mediana = median(data$Claridad, na.rm = TRUE),
Claridad_Media = mean(data$Claridad, na.rm = TRUE),
Claridad_Q3 = quantile(data$Claridad, 0.75, na.rm = TRUE),
Claridad_Max = max(data$Claridad, na.rm = TRUE),
Claridad_Sd = sd(data$Claridad, na.rm = TRUE),
Claridad_IQR = IQR(data$Claridad, na.rm = TRUE),
Aroma_Min = min(data$Aroma, na.rm = TRUE),
Aroma_Q1 = quantile(data$Aroma, 0.25, na.rm = TRUE),
Aroma_Mediana = median(data$Aroma, na.rm = TRUE),
Aroma_Media = mean(data$Aroma, na.rm = TRUE),
Aroma_Q3 = quantile(data$Aroma, 0.75, na.rm = TRUE),
Aroma_Max = max(data$Aroma, na.rm = TRUE),
Aroma_Sd = sd(data$Aroma, na.rm = TRUE),
Aroma_IQR = IQR(data$Aroma, na.rm = TRUE),
Cuerpo_Min = min(data$Cuerpo, na.rm = TRUE),
Cuerpo_Q1 = quantile(data$Cuerpo, 0.25, na.rm = TRUE),
Cuerpo_Mediana = median(data$Cuerpo, na.rm = TRUE),
Cuerpo_Media = mean(data$Cuerpo, na.rm = TRUE),
Cuerpo_Q3 = quantile(data$Cuerpo, 0.75, na.rm = TRUE),
Cuerpo_Max = max(data$Cuerpo, na.rm = TRUE),
Cuerpo_Sd = sd(data$Cuerpo, na.rm = TRUE),
Cuerpo_IQR = IQR(data$Cuerpo, na.rm = TRUE),
Sabor_Min = min(data$Sabor, na.rm = TRUE),
Sabor_Q1 = quantile(data$Sabor, 0.25, na.rm = TRUE),
Sabor_Mediana = median(data$Sabor, na.rm = TRUE),
Sabor_Media = mean(data$Sabor, na.rm = TRUE),
Sabor_Q3 = quantile(data$Sabor, 0.75, na.rm = TRUE),
Sabor_Max = max(data$Sabor, na.rm = TRUE),
Sabor_Sd = sd(data$Sabor, na.rm = TRUE),
Sabor_IQR = IQR(data$Sabor, na.rm = TRUE),
Roble_Min = min(data$Roble, na.rm = TRUE),
Roble_Q1 = quantile(data$Roble, 0.25, na.rm = TRUE),
Roble_Mediana = median(data$Roble, na.rm = TRUE),
Roble_Media = mean(data$Roble, na.rm = TRUE),
Roble_Q3 = quantile(data$Roble, 0.75, na.rm = TRUE),
Roble_Max = max(data$Roble, na.rm = TRUE),
Roble_Sd = sd(data$Roble, na.rm = TRUE),
Roble_IQR = IQR(data$Roble, na.rm = TRUE)
) %>%
pivot_longer(everything(), names_to = c("Variable", "Estadístico"), names_sep = "_") %>%
pivot_wider(names_from = Estadístico, values_from = value)
# Crear tabla con kableExtra
tabla_estadisticos %>%
kbl(caption = "Resumen Estadístico de las Variables de Calidad del Vino", digits = 2, align = "c") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE) %>%
column_spec(1, bold = TRUE, color = "white", background = "#8B0000") %>% # Resaltar Categorías
column_spec(2:9, background = "#FFF5EE")| Variable | Min | Q1 | Mediana | Media | Q3 | Max | Sd | IQR |
|---|---|---|---|---|---|---|---|---|
| Calidad | 7.9 | 11.15 | 12.45 | 12.44 | 13.75 | 16.1 | 2.05 | 2.60 |
| Claridad | 0.5 | 0.83 | 1.00 | 0.92 | 1.00 | 1.0 | 0.12 | 0.17 |
| Aroma | 3.3 | 4.12 | 4.65 | 4.85 | 5.45 | 7.7 | 1.08 | 1.33 |
| Cuerpo | 2.6 | 4.15 | 4.75 | 4.68 | 5.38 | 6.6 | 0.82 | 1.23 |
| Sabor | 2.9 | 4.23 | 4.80 | 4.77 | 5.50 | 7.0 | 1.03 | 1.27 |
| Roble | 2.9 | 3.70 | 4.10 | 4.26 | 4.77 | 6.0 | 0.74 | 1.07 |
Análisis de Calidad del Vino
El puntaje de calidad presentó una media de 12.44 ± 2.05, donde el 75% de los vinos obtuvo un valor máximo de 13.75. En cuanto a la claridad, la mediana fue de 1.0, con un rango intercuartílico de 0.17, indicando homogeneidad, y el tercer cuartil confirmó que el 75% de las muestras alcanzaron el valor máximo (1.0) con una alta claridad. El aroma registró una media de 4.85 ± 1.08, con un 75% de los vinos ubicándose hasta en un puntaje de 5.45, es decir un aroma casi intenso. Para el cuerpo, la media fue 4.68 ± 0.82, y el 75% de los casos no superó 5.38, es decir tienen una densidad un poco robusta. El sabor mostró una media de 4.77 ± 1.03, con el tercer cuartil en 5.50 con sabores complejos, mientras que la influencia de roble promedió 4.26 ± 0.74, con el 75% de los vinos presentando notas no tan marcadas de a lo más 4.77 de puntaje. Estos resultados destacan la consistencia en claridad y la mayor variabilidad en calidad.
chart.Correlation(data, histogram = TRUE, method = "spearman")
mat_cor <- cor(data, method = "spearman") # Calcula matriz de correlación
mat_cor1<-cor(datos, method = "spearman")
significancia1<- cor.mtest(datos,
method = "spearman",
conf.level = .95)
corrplot(mat_cor1,
p.mat = significancia1$p, #llamado del p-valor para cada coeficiente r
sig.level = 0.05) #definición del nivel de significancia
En cuanto a la distribución, todas parecen seguir una distribución normal, excepto por claridad en la cual se observa una ligera asimetría a la izquierda, lo que explica que no haya pasado la prueba de normalidad. A continuación, se realiza un análisis exhaustivo de la relación de las variables explicativas en torno a la variable respuesta.
Calidad vs Claridad: En el diagrama de dispersión no es clara una correlación, pues parecen estar acumulados cuando la claridad es de 1.0, lo que es consistente con lo dicho previamente de claridad, pues la mayoría de los datos tienen un puntaje de 1.0. Por otro lado, tienen una correlación débil, pues el coeficiente de correlación tiene un valor de 0.02.
Calidad vs Aroma: En el diagrama de dispersión se observa una correlación positiva, a medida que aumenta la intensidad del aroma, la calidad del vino también incrementa, lo cual se confirma con un coeficiente de correlación de 0.68, lo que muestra una correlación positiva moderada, asimismo es evidente una significancia alta.
Calidad vs Cuerpo: En el diagrama de dispersión se evidencia una correlación positiva, a medida que incrementa la sensación de densidad del vino se vuelve más robusta, la calidad de este también aumenta, esto se confirma con un coeficiente de correlación moderado con un valor de 0.54 y una significancia alta.
Calidad vs Sabor: En el diagrama de dispersión se aprecia una correlación positiva, a medida que aumenta la complejidad del sabor, la calidad del vino también incrementa, lo cual se evidencia con un coeficiente de correlación alto de 0.71, el mayor de todos los coeficientes, y una significancia alta.
Calidad vs Roble: En el diagrama de dispersión no se aprecia una correlación o patrón, los datos de la calidad parecen no oscilar tanto a medida que aumenta la influencia del roble, lo que se constata con un coeficiente de correlación débil de 0.27 con una significancia muy baja.
Siguiendo otro orden de ideas, es importante resaltar que existe una correlación positiva moderada entre cuerpo y aroma con coeficiente de correlación de 0.51 y una significancia un poco alta. Asimismo, se percibe una correlación positiva moderada entre cuerpo y sabor con un coeficiente de correlación de 0.59 y una significancia alta. Finalmente, la más preocupante es la correlación positiva fuerte entre aroma y sabor con un coeficiente de correlación de 0.74 con alta significancia. Esto sugiere que, a la hora de hacer el modelo estas variables estrictamente relacionadas entre sí es posible que no contribuyan de manera independiente a explicar la variable respuesta, ya que al estar tan correlacionadas pueden resultar redundantes y generar problemas de multicolinealidad.
SELECCIÓN DE VARIABLES
Para la selección de variables, se empleó el método both basado en el Criterio de Información de Akaike (AIC), que optimiza el balance entre la bondad de ajuste, pues valúa qué tan bien el modelo explica la variabilidad de la calidad del vino, y la parsimonia del modelo, ya que castiga la inclusión de variables innecesarias, evitando el sobreajuste. Hay que tener en cuenta que cuando la diferencia entre AICs es muy pequeña (ej. 0.5), en ocasiones se prioriza la parsimonia.
Backward
A continuación, aplicamos el método backward. con p-value————————————————————-
modelo <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = data)
summary(modelo)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.85552 -0.57448 -0.07092 0.67275 1.68093
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.9969 2.2318 1.791 0.082775 .
## x1 2.3395 1.7348 1.349 0.186958
## x2 0.4826 0.2724 1.771 0.086058 .
## x3 0.2732 0.3326 0.821 0.417503
## x4 1.1683 0.3045 3.837 0.000552 ***
## x5 -0.6840 0.2712 -2.522 0.016833 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.163 on 32 degrees of freedom
## Multiple R-squared: 0.7206, Adjusted R-squared: 0.6769
## F-statistic: 16.51 on 5 and 32 DF, p-value: 4.703e-08La variable con el valor-p más alto es x3 = 0.418, así que es la menos significativa y la que debemos eliminar primero.
modelo2 <- lm(y ~ x1 + x2 + x4 + x5, data = data)
summary(modelo2)##
## Call:
## lm(formula = y ~ x1 + x2 + x4 + x5, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.56069 -0.51239 0.00782 0.69037 1.80377
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.9855 1.8701 2.666 0.0118 *
## x1 1.7942 1.5949 1.125 0.2687
## x2 0.5300 0.2649 2.000 0.0537 .
## x4 1.2643 0.2798 4.519 7.55e-05 ***
## x5 -0.6589 0.2681 -2.457 0.0194 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.157 on 33 degrees of freedom
## Multiple R-squared: 0.7147, Adjusted R-squared: 0.6801
## F-statistic: 20.67 on 4 and 33 DF, p-value: 1.316e-08La variable con el valor-p más alto ahora es x1 = 0.2687, así que es la siguiente candidata a ser eliminada del modelo.
modelo3 <- lm(y ~ x2 + x4 + x5, data = data)
summary(modelo3)##
## Call:
## lm(formula = y ~ x2 + x4 + x5, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5707 -0.6256 0.1521 0.6467 1.7741
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.4672 1.3328 4.852 2.67e-05 ***
## x2 0.5801 0.2622 2.213 0.033740 *
## x4 1.1997 0.2749 4.364 0.000113 ***
## x5 -0.6023 0.2644 -2.278 0.029127 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.161 on 34 degrees of freedom
## Multiple R-squared: 0.7038, Adjusted R-squared: 0.6776
## F-statistic: 26.92 on 3 and 34 DF, p-value: 4.203e-09Todos los predictores tienen valores-p menores a 0.05, lo que significa que todos son estadísticamente significativos. El modelo tiene: R² ajustado: 0.6776 (ligeramente menor que el anterior: 0.6801) Error estándar residual: 1.161 (prácticamente igual al anterior: 1.157)
modelo_final <- lm(y ~ x2 + x4 + x5, data = data)La eliminación de x1 no mejoró el R² ajustado, pero lo mantuvo prácticamente igual, y como x1 no era significativo, su exclusión está justificada.
Con AIC—————————————————————————
modelo_completo <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = data)
AIC(modelo_completo)## [1] 126.7552Ahora, elimina una variable a la vez y calcula el AIC de cada modelo reducido.
m1<-lm(y ~ x2 + x3 + x4 + x5, data = data) # sin x1
AIC(m1)## [1] 126.8556m2<-lm(y ~ x1 + x3 + x4 + x5, data = data) # sin x2
AIC(m2)## [1] 128.309m3<-lm(y ~ x1 + x2 + x4 + x5, data = data) # sin x3
AIC(m3)## [1] 125.548m4<-lm(y ~ x1 + x2 + x3 + x5, data = data) # sin x4
AIC(m4)## [1] 139.1385m5<-lm(y ~ x1 + x2 + x3 + x4, data = data) # sin x5
AIC(m5)## [1] 131.6454El único modelo con AIC menor al del modelo completo (126.7552) es el m3, con AIC = 125.548. Entonces eliminamos x3.
modelo_base <- lm(y ~ x1 + x2 + x4 + x5, data = data)
AIC(modelo_base) ## [1] 125.548m11<-lm(y ~ x2 + x4 + x5, data = data)
AIC(m11)## [1] 124.9781m22<-lm(y ~ x1 + x4 + x5, data = data)
AIC(m22)## [1] 127.897m44<-lm(y ~ x1 + x2 + x5, data = data)
AIC(m44)## [1] 141.8531m55<-lm(y ~ x1 + x2 + x4, data = data)
AIC(m55)## [1] 129.9337Eliminamos x1, porque es la que al eliminarla reduce más el AIC (de 125.548 → 124.9781).
modelo_base2 <- lm(y ~ x2 + x4 + x5, data = data)
AIC(modelo_base2)## [1] 124.9781m222<-lm(y ~ x4 + x5, data = data)
AIC(m222)## [1] 128.0901m444<-lm(y ~ x2 + x5, data = data)
AIC(m444)## [1] 139.882m555<-lm(y ~ x2 + x4, data = data)
AIC(m555)## [1] 128.3761Ninguna eliminación mejora el AIC actual (124.9781). De hecho, todos los AICs aumentan si eliminamos cualquier variable.
El mejor modelo según el criterio de AIC por eliminación backward es:
modelo_final <- lm(y ~ x2 + x4 + x5, data = data)
summary(modelo_final)##
## Call:
## lm(formula = y ~ x2 + x4 + x5, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5707 -0.6256 0.1521 0.6467 1.7741
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.4672 1.3328 4.852 2.67e-05 ***
## x2 0.5801 0.2622 2.213 0.033740 *
## x4 1.1997 0.2749 4.364 0.000113 ***
## x5 -0.6023 0.2644 -2.278 0.029127 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.161 on 34 degrees of freedom
## Multiple R-squared: 0.7038, Adjusted R-squared: 0.6776
## F-statistic: 26.92 on 3 and 34 DF, p-value: 4.203e-09A partir de los coeficientes estimados, la ecuación del modelo de regresión lineal es:
\[ y = 6.467 + 0.5801*x2 + 1.1997*x4 -
0.6023*x5\] Interpretación de los
coeficientes
- Intercepto (6.4672): Valor esperado de y cuando x2, x4 y x5 son cero.
Su valor es significativo (p < 0.001).
- x2 (0.5801): Por cada unidad que aumenta el aroma, la calidad del vino
aumenta en promedio 0.58 unidades, manteniendo constantes las otras
variables. Es significativo (p ≈ 0.034).
- x4 (1.1997): El más influyente positivamente. Por cada unidad que
aumenta el sabor, la calidad del vino sube en promedio 1.20 unidades.
Muy significativo (p < 0.001).
- x5 (−0.6023): Tiene un efecto negativo. Por cada unidad que aumenta la
influencia del roble en el vino, su calidad disminuye en promedio 0.60
unidades, con las otras constantes. Es significativo (p ≈ 0.029).
Calidad del ajuste - R² = 0.7038: El modelo explica aproximadamente el 70.38% de la variabilidad de y, lo cual es un buen ajuste. - R² ajustado = 0.6776: Considerando el número de predictores, sigue siendo sólido. - Error estándar residual = 1.161: En promedio, los valores predichos se desvían en ±1.161 unidades de los valores reales de y.
Significancia global del modelo F-statistic = 26.92, con p < 0.00001: El modelo como un todo es altamente significativo, lo que indica que al menos una variable predictora explica y.
modelo_completo <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = data)
modelo_backward <- stepAIC(modelo_completo, direction = "backward", trace = TRUE)## Start: AIC=16.92
## y ~ x1 + x2 + x3 + x4 + x5
##
## Df Sum of Sq RSS AIC
## - x3 1 0.9118 44.160 15.709
## <none> 43.248 16.916
## - x1 1 2.4577 45.706 17.016
## - x2 1 4.2397 47.488 18.470
## - x5 1 8.5978 51.846 21.806
## - x4 1 19.8986 63.147 29.299
##
## Step: AIC=15.71
## y ~ x1 + x2 + x4 + x5
##
## Df Sum of Sq RSS AIC
## - x1 1 1.6936 45.853 15.139
## <none> 44.160 15.709
## - x2 1 5.3545 49.514 18.058
## - x5 1 8.0807 52.241 20.094
## - x4 1 27.3280 71.488 32.014
##
## Step: AIC=15.14
## y ~ x2 + x4 + x5
##
## Df Sum of Sq RSS AIC
## <none> 45.853 15.139
## - x2 1 6.6026 52.456 18.251
## - x5 1 6.9989 52.852 18.537
## - x4 1 25.6888 71.542 30.043summary(modelo_backward)##
## Call:
## lm(formula = y ~ x2 + x4 + x5, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5707 -0.6256 0.1521 0.6467 1.7741
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.4672 1.3328 4.852 2.67e-05 ***
## x2 0.5801 0.2622 2.213 0.033740 *
## x4 1.1997 0.2749 4.364 0.000113 ***
## x5 -0.6023 0.2644 -2.278 0.029127 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.161 on 34 degrees of freedom
## Multiple R-squared: 0.7038, Adjusted R-squared: 0.6776
## F-statistic: 26.92 on 3 and 34 DF, p-value: 4.203e-09modelo_backward$anova## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## y ~ x1 + x2 + x3 + x4 + x5
##
## Final Model:
## y ~ x2 + x4 + x5
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 32 43.24801 16.91587
## 2 - x3 1 0.911827 33 44.15983 15.70872
## 3 - x1 1 1.693582 34 45.85341 15.13881Nota: El AIC que devuelve la función AIC() directamente sobre un modelo lineal de lm() está ajustado usando la deviance residual (específicamente: n * log(RSS / n) + 2k, donde k es el número de parámetros). - En cambio, step() trabaja con una versión escalada del log-likelihood, lo que puede hacer que el valor numérico del AIC sea diferente, aunque la comparación entre modelos sea consistente.
Desarrollar el proceso de forma manual permitió comprender con más profundidad cómo influye cada variable en el modelo. Por otro lado, la función automática (step) agilizó el proceso y lo sistematizó, lo que es muy útil cuando se trabaja con más variables o modelos complejos. En este caso, ambos caminos convergieron, confirmando que el modelo más adecuado según el AIC es aquel que incluye únicamente x2, x4 y x5.
Fordward
A continuación se procederá a aplicar el método forward con cada una de las variables, utilizando como métrica de selección el AIC. Primero, lo realizaremos de forma “manual”.
Para empezar, aplicamos el primer paso, es decir, creamos el modelo nulo:
M0 <- lm(y ~ 1)
summary(M0)##
## Call:
## lm(formula = y ~ 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5368 -1.2868 0.0132 1.3132 3.6632
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.4368 0.3318 37.48 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.045 on 37 degrees of freedomEl modelo nulo se compone únicamente de su intercepto y su error.
Procederemos con el paso 2:
#Creamos diversos modelos con el M0
C1 <- lm(y ~ x1)
C2 <- lm(y ~ x2)
C3 <- lm(y ~ x3)
C4 <- lm(y ~ x4)
C5 <- lm(y ~ x5)
#Hacemos summary con cada uno y calculamos su AIC
summary(C1)##
## Call:
## lm(formula = y ~ x1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5257 -1.3227 0.0947 1.2773 3.7681
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.0034 2.5610 4.687 3.89e-05 ***
## x1 0.4692 2.7486 0.171 0.865
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.073 on 36 degrees of freedom
## Multiple R-squared: 0.0008089, Adjusted R-squared: -0.02695
## F-statistic: 0.02914 on 1 and 36 DF, p-value: 0.8654AIC(C1)## [1] 167.1786summary(C2)##
## Call:
## lm(formula = y ~ x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4726 -0.8574 -0.0091 0.8346 2.2563
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.9583 1.1050 5.392 4.51e-06 ***
## x2 1.3365 0.2226 6.004 6.87e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.466 on 36 degrees of freedom
## Multiple R-squared: 0.5003, Adjusted R-squared: 0.4864
## F-statistic: 36.04 on 1 and 36 DF, p-value: 6.871e-07AIC(C2)## [1] 140.8463summary(C3)##
## Call:
## lm(formula = y ~ x3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.9669 -0.8386 0.0620 1.2204 3.4502
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.0580 1.6441 3.685 0.000748 ***
## x3 1.3618 0.3458 3.938 0.000361 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.734 on 36 degrees of freedom
## Multiple R-squared: 0.3011, Adjusted R-squared: 0.2817
## F-statistic: 15.51 on 1 and 36 DF, p-value: 0.0003612AIC(C3)## [1] 153.5973summary(C4)##
## Call:
## lm(formula = y ~ x4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.38583 -0.72226 -0.00756 0.62006 2.52822
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.9414 0.9911 4.986 1.57e-05 ***
## x4 1.5719 0.2033 7.732 3.68e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.271 on 36 degrees of freedom
## Multiple R-squared: 0.6242, Adjusted R-squared: 0.6137
## F-statistic: 59.79 on 1 and 36 DF, p-value: 3.683e-09AIC(C4)## [1] 130.0214summary(C5)##
## Call:
## lm(formula = y ~ x5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6483 -1.3886 -0.0527 1.2907 3.6429
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.9916 1.9918 6.522 1.4e-07 ***
## x5 -0.1304 0.4614 -0.283 0.779
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.071 on 36 degrees of freedom
## Multiple R-squared: 0.002213, Adjusted R-squared: -0.0255
## F-statistic: 0.07984 on 1 and 36 DF, p-value: 0.7791AIC(C5)## [1] 167.1251Luego en el paso 3, para determinar cuál modelo es más significativo, nos basamos en el valor del AIC calculado para cada uno. Seleccionamos el modelo cuyo AIC sea el más pequeño y lo comparamos con el modelo anterior. En este primer paso, como el modelo anterior es el modelo nulo (que solo incluye el intercepto), no es necesario realizar una comparación, simplemente seleccionamos el modelo con el AIC más bajo y continuamos con el proceso.
El modelo con el AIC más pequeño fue el de la variable explicativa Sabor (\(x_4\)), por lo tanto, nuestro nuevo modelo queda de la siguiente manera:
M1 <- lm(y ~ x4)
summary(M1)##
## Call:
## lm(formula = y ~ x4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.38583 -0.72226 -0.00756 0.62006 2.52822
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.9414 0.9911 4.986 1.57e-05 ***
## x4 1.5719 0.2033 7.732 3.68e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.271 on 36 degrees of freedom
## Multiple R-squared: 0.6242, Adjusted R-squared: 0.6137
## F-statistic: 59.79 on 1 and 36 DF, p-value: 3.683e-09AIC(M1)## [1] 130.0214Repetimos el paso 2, pero esta vez utilizando como base el modelo M1. Es decir:
#Creamos diversos modelos con el M1
C1 <- lm(y ~ x4 + x1)
C2 <- lm(y ~ x4 + x2)
C3 <- lm(y ~ x4 + x3)
C4 <- lm(y ~ x4 + x5)
#Hacemos summary y calculamos su AIC
summary(C1)##
## Call:
## lm(formula = y ~ x4 + x1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.32580 -0.81428 0.07715 0.92907 2.42101
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.3941 1.9241 1.764 0.0865 .
## x4 1.5882 0.2044 7.771 3.97e-09 ***
## x1 1.5908 1.6946 0.939 0.3543
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.273 on 35 degrees of freedom
## Multiple R-squared: 0.6334, Adjusted R-squared: 0.6125
## F-statistic: 30.24 on 2 and 35 DF, p-value: 2.362e-08AIC(C1)## [1] 131.0765summary(C2)##
## Call:
## lm(formula = y ~ x4 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.19048 -0.60300 -0.03203 0.66039 2.46287
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.3462 1.0091 4.307 0.000127 ***
## x4 1.1702 0.2905 4.027 0.000288 ***
## x2 0.5180 0.2759 1.877 0.068849 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.229 on 35 degrees of freedom
## Multiple R-squared: 0.6586, Adjusted R-squared: 0.639
## F-statistic: 33.75 on 2 and 35 DF, p-value: 6.811e-09AIC(C2)## [1] 128.3761summary(C3)##
## Call:
## lm(formula = y ~ x4 + x3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.55118 -0.71186 0.07937 0.57229 2.51758
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.5846 1.2475 3.675 0.00079 ***
## x4 1.4883 0.2694 5.524 3.27e-06 ***
## x3 0.1613 0.3361 0.480 0.63426
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.285 on 35 degrees of freedom
## Multiple R-squared: 0.6266, Adjusted R-squared: 0.6053
## F-statistic: 29.37 on 2 and 35 DF, p-value: 3.254e-08AIC(C3)## [1] 131.7721summary(C4)##
## Call:
## lm(formula = y ~ x4 + x5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.75329 -0.73280 -0.01534 0.82829 2.05325
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.9122 1.3889 4.977 1.72e-05 ***
## x4 1.6418 0.1990 8.249 1.01e-09 ***
## x5 -0.5414 0.2772 -1.953 0.0588 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.224 on 35 degrees of freedom
## Multiple R-squared: 0.6611, Adjusted R-squared: 0.6417
## F-statistic: 34.14 on 2 and 35 DF, p-value: 5.97e-09AIC(C4)## [1] 128.0901Básicamente, lo que hicimos en este paso fue evaluar cada una de las variables explicativas restantes en combinación con el modelo M1, y observar cuál de ellas producía el menor valor de AIC. Luego, comparamos ese AIC con el del modelo anterior. Lo que nos dio como resultado:
\[128.09 (C_4) < 130.02(M_1)\] Como el AIC del modelo \(C_4\) es menor que el AIC del modelo \(M_1\), concluimos que hemos encontrado un modelo más significativo.
En este caso, la variable más significativa fue Roble (\(x_5\)), por lo cual se agrega al modelo:
M2 <- lm(y ~ x4 + x5)
summary(M2)##
## Call:
## lm(formula = y ~ x4 + x5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.75329 -0.73280 -0.01534 0.82829 2.05325
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.9122 1.3889 4.977 1.72e-05 ***
## x4 1.6418 0.1990 8.249 1.01e-09 ***
## x5 -0.5414 0.2772 -1.953 0.0588 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.224 on 35 degrees of freedom
## Multiple R-squared: 0.6611, Adjusted R-squared: 0.6417
## F-statistic: 34.14 on 2 and 35 DF, p-value: 5.97e-09AIC(M2)## [1] 128.0901Volvemos a repetir el procedimiento, ahora utilizando como base el modelo M2:
C1 <- lm(y ~ x4 + x5 + x1)
C2 <- lm(y ~ x4 + x5 + x2)
C3 <- lm(y ~ x4 + x5 + x3)
summary(C1)##
## Call:
## lm(formula = y ~ x4 + x5 + x1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.71984 -0.85335 -0.01423 0.85779 1.89034
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.9374 1.9507 2.531 0.0162 *
## x4 1.6761 0.1977 8.480 6.66e-10 ***
## x5 -0.6218 0.2790 -2.228 0.0326 *
## x1 2.3310 1.6401 1.421 0.1644
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.207 on 34 degrees of freedom
## Multiple R-squared: 0.6801, Adjusted R-squared: 0.6519
## F-statistic: 24.1 on 3 and 34 DF, p-value: 1.527e-08AIC(C1)## [1] 127.897summary(C2)##
## Call:
## lm(formula = y ~ x4 + x5 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5707 -0.6256 0.1521 0.6467 1.7741
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.4672 1.3328 4.852 2.67e-05 ***
## x4 1.1997 0.2749 4.364 0.000113 ***
## x5 -0.6023 0.2644 -2.278 0.029127 *
## x2 0.5801 0.2622 2.213 0.033740 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.161 on 34 degrees of freedom
## Multiple R-squared: 0.7038, Adjusted R-squared: 0.6776
## F-statistic: 26.92 on 3 and 34 DF, p-value: 4.203e-09AIC(C2)## [1] 124.9781summary(C3)##
## Call:
## lm(formula = y ~ x4 + x5 + x3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.95511 -0.60797 -0.00465 0.81722 2.09109
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.5172 1.5526 4.198 0.000183 ***
## x4 1.5435 0.2606 5.923 1.09e-06 ***
## x5 -0.5494 0.2801 -1.961 0.058106 .
## x3 0.1916 0.3235 0.592 0.557620
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.236 on 34 degrees of freedom
## Multiple R-squared: 0.6646, Adjusted R-squared: 0.635
## F-statistic: 22.45 on 3 and 34 DF, p-value: 3.385e-08AIC(C3)## [1] 129.7001En este caso, el menor AIC se obtiene al agregar la variable Aroma (\(x_2\)), y comparando el AIC del modelo anterior con este nuevo obtenemos que:
\[124.97 < 128.09\]
Por lo cual se incluye en el nuevo modelo:
M3 <- lm(y ~ x4 + x5 + x2)
summary(M3)##
## Call:
## lm(formula = y ~ x4 + x5 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5707 -0.6256 0.1521 0.6467 1.7741
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.4672 1.3328 4.852 2.67e-05 ***
## x4 1.1997 0.2749 4.364 0.000113 ***
## x5 -0.6023 0.2644 -2.278 0.029127 *
## x2 0.5801 0.2622 2.213 0.033740 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.161 on 34 degrees of freedom
## Multiple R-squared: 0.7038, Adjusted R-squared: 0.6776
## F-statistic: 26.92 on 3 and 34 DF, p-value: 4.203e-09AIC(M3)## [1] 124.9781Por último, volvemos a aplicar el mismo proceso sobre el modelo M3:
C1 <- lm(y ~ x4 + x5 + x2 + x1)
C2 <- lm(y ~ x4 + x5 + x2 + x3)
summary(C1)##
## Call:
## lm(formula = y ~ x4 + x5 + x2 + x1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.56069 -0.51239 0.00782 0.69037 1.80377
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.9855 1.8701 2.666 0.0118 *
## x4 1.2643 0.2798 4.519 7.55e-05 ***
## x5 -0.6589 0.2681 -2.457 0.0194 *
## x2 0.5300 0.2649 2.000 0.0537 .
## x1 1.7942 1.5949 1.125 0.2687
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.157 on 33 degrees of freedom
## Multiple R-squared: 0.7147, Adjusted R-squared: 0.6801
## F-statistic: 20.67 on 4 and 33 DF, p-value: 1.316e-08AIC(C1)## [1] 125.548summary(C2)##
## Call:
## lm(formula = y ~ x4 + x5 + x2 + x3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6814 -0.6830 0.1559 0.5860 1.7725
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.2670 1.4833 4.225 0.000177 ***
## x4 1.1567 0.3081 3.754 0.000672 ***
## x5 -0.6053 0.2681 -2.258 0.030706 *
## x2 0.5682 0.2682 2.118 0.041769 *
## x3 0.1016 0.3110 0.327 0.746074
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.177 on 33 degrees of freedom
## Multiple R-squared: 0.7047, Adjusted R-squared: 0.6689
## F-statistic: 19.69 on 4 and 33 DF, p-value: 2.292e-08AIC(C2)## [1] 126.8556Comparando los valores de AIC obtenidos, observamos que el menor corresponde al modelo C1. Al compararlo con el modelo anterior M3, obtenemos:
\[125.54(C1) > 124.97 (M3)\]
Como en esta iteración ningún modelo presenta un AIC menor que el del modelo anterior, se concluye el proceso de selección. Nuestro modelo final queda definido como:
\[y = 6.46 + 0.58x_2 + 1.19x_4 - 0.60x_5\]
Todo este proceso se realizó de forma manual, sin embargo, existe una función que permite simplificar significativamente la selección del modelo:
Modelo_Nulo <- lm(Calidad ~ 1, data = data) #Utilizamos ~ 1, para decir que solo sea la variable respuesta con el intercepto.
Modelo_Completo <- lm(Calidad ~ ., data = data) #La expresión ~ ., representa el modelo con todas las variables explicativas.
modelo_forward <- stepAIC(Modelo_Nulo, scope = list(lower = formula(Modelo_Nulo), upper = formula(Modelo_Completo)), direction = "forward") #Utilizamos stepAIC para que su criterio de selección sea por el AIC.## Start: AIC=55.37
## Calidad ~ 1
##
## Df Sum of Sq RSS AIC
## + Sabor 1 96.615 58.173 20.182
## + Aroma 1 77.442 77.347 31.007
## + Cuerpo 1 46.603 108.186 43.758
## <none> 154.788 55.370
## + Roble 1 0.343 154.446 57.286
## + Claridad 1 0.125 154.663 57.339
##
## Step: AIC=20.18
## Calidad ~ Sabor
##
## Df Sum of Sq RSS AIC
## + Roble 1 5.7174 52.456 18.251
## + Aroma 1 5.3212 52.852 18.537
## <none> 58.173 20.182
## + Claridad 1 1.4286 56.745 21.237
## + Cuerpo 1 0.3803 57.793 21.933
##
## Step: AIC=18.25
## Calidad ~ Sabor + Roble
##
## Df Sum of Sq RSS AIC
## + Aroma 1 6.6026 45.853 15.139
## + Claridad 1 2.9416 49.514 18.058
## <none> 52.456 18.251
## + Cuerpo 1 0.5356 51.920 19.861
##
## Step: AIC=15.14
## Calidad ~ Sabor + Roble + Aroma
##
## Df Sum of Sq RSS AIC
## <none> 45.853 15.139
## + Claridad 1 1.69358 44.160 15.709
## + Cuerpo 1 0.14769 45.706 17.016modelo_forward #Lo llamamos.##
## Call:
## lm(formula = Calidad ~ Sabor + Roble + Aroma, data = data)
##
## Coefficients:
## (Intercept) Sabor Roble Aroma
## 6.4672 1.1997 -0.6023 0.5801Como vemos el resultado coincide con el proceso hecho manualmente. Procedamos a desglosar este modelo:
modelo_forward$anova## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## Calidad ~ 1
##
## Final Model:
## Calidad ~ Sabor + Roble + Aroma
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 37 154.78842 55.36997
## 2 + Sabor 1 96.614981 36 58.17344 20.18202
## 3 + Roble 1 5.717446 35 52.45599 18.25076
## 4 + Aroma 1 6.602579 34 45.85341 15.13881Se pueden destacar varias observaciones importantes:
Primera Fila: Representa el modelo nulo, es decir, donde solo se incluye el intercepto. Aparece el número de grados de libertad (\(Df\)), que en este caso es 37, recordando que se calcula n - p (donde n es el número de observaciones y p el número de variables explicativas). También observamos la Deviancia Residual de 154.78, que mide qué tan mal se ajusta el modelo a los datos: mientras menor sea, mejor el ajuste. Finalmente, tenemos el valor del AIC, que también buscamos minimizar al añadir variables al modelo.
Segunda Fila: Se observa una disminución significativa en la deviancia residual,lo cual indica que al agregar la variable explicativa Sabor, el modelo mejora considerablemente en su capacidad de predecir los datos observados. Además, el AIC también disminuye de forma notable, lo cual refuerza esta mejora.
Tercera fila: Continúa la tendencia observada: al agregar otra variable (por ejemplo, Roble), tanto la deviancia residual como el AIC siguen disminuyendo, indicando que el nuevo modelo se ajusta mejor.
Cuarta Fila: Se repite el mismo patrón antes mencionado, al agregar la variable Aroma, disminuye la deviancia residual y el AIC.
Cabe destacar que, en general, cuanto más cercana a cero sea la Deviancia Residual, mejor será el ajuste del modelo a los datos observados. Asimismo, una disminución progresiva del AIC sugiere que cada nueva variable está aportando valor predictivo sin sobreajustar el modelo.
Observemos más de cerca el modelo, para esto hacemos un summary:
summary(modelo_forward)##
## Call:
## lm(formula = Calidad ~ Sabor + Roble + Aroma, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5707 -0.6256 0.1521 0.6467 1.7741
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.4672 1.3328 4.852 2.67e-05 ***
## Sabor 1.1997 0.2749 4.364 0.000113 ***
## Roble -0.6023 0.2644 -2.278 0.029127 *
## Aroma 0.5801 0.2622 2.213 0.033740 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.161 on 34 degrees of freedom
## Multiple R-squared: 0.7038, Adjusted R-squared: 0.6776
## F-statistic: 26.92 on 3 and 34 DF, p-value: 4.203e-09Both
Se procede a aplicar el método forward con cada una de las variables.
m1<-lm(y~x1)
AIC(m1)## [1] 167.1786m2<-lm(y~x2)
AIC(m2)## [1] 140.8463m3<-lm(y~x3)
AIC(m3)## [1] 153.5973m4<-lm(y~x4) #se escoge este
AIC(m4)## [1] 130.0214m5<-lm(y~x5)
AIC(m5)## [1] 167.1251El menor AIC es el de 130.0214 del modelo 4 con la variable sabor, por tanto lo elegimos y aplicamos backward.
m4<-lm(y~x4)
AIC(m4)## [1] 130.0214m5<-lm(y~1)
AIC(m5)## [1] 165.2093Al eliminar a sabor no mejora el AIC, por tanto se mantiene en el modelo. Se aplica el método forward.
m6<-lm(y~x4+x1)
AIC(m6)## [1] 131.0765m7<-lm(y~x4+x2)
AIC(m7)## [1] 128.3761m8<-lm(y~x4+x3)
AIC(m8)## [1] 131.7721m9<-lm(y~x4+x5) #se escoge este
AIC(m9)## [1] 128.0901El mejor AIC es el de 128.0901 que corresponde al modelo con las variables sabor y roble. Se verifica con backward si eliminar alguna mejora el AIC.
m9<-lm(y~x4+x5)
AIC(m9)## [1] 128.0901m10<-lm(y~x4)
AIC(m10)## [1] 130.0214m11<-lm(y~x5)
AIC(m11)## [1] 167.1251El AIC empeora en cualquiera de los dos casos, se mantienen las dos variables y se repite el proceso con el método forward.
m12<-lm(y~x4+x5+x1)
AIC(m12)## [1] 127.897m13<-lm(y~x4+x5+x2) #se escoge este
AIC(m13)## [1] 124.9781m14<-lm(y~x4+x5+x3)
AIC(m14)## [1] 129.7001El menor AIC es el del modelo con las variables sabor, roble y aroma. Se realiza el método backward:
m13<-lm(y~x4+x5+x2)
AIC(m13)## [1] 124.9781m14<-lm(y~x5+x2) #sin x4
AIC(m14)## [1] 139.882m15<-lm(y~x4+x2) #sin x5
AIC(m15)## [1] 128.3761m16<-lm(y~x4+x5) #sin x2
AIC(m16)## [1] 128.0901En ninguno de los casos el AIC disminuye, por tanto se mantienen las variables y se realiza el método forward.
m17<-lm(y~x4+x5+x2+x1)
AIC(m17)## [1] 125.548m18<-lm(y~x4+x5+x2+x3)
AIC(m18)## [1] 126.8556El AIC no disminuye, por tanto el proceso termina, pues no se puede añadir ninguna variable que reduzca el AIC ni se puede eliminar ninguna variable para mejorar el AIC. Por tanto el modelo final incluye a las variables sabor, roble y aroma.
No obstante, existe una función que condensa todo este proceso, se puede realizar de la siguiente forma:
modelo_inicial<-lm(y~1)
scope = list(lower = ~1, upper = y ~ x1 + x2 + x3 + x4 + x5)
modelo_both <- stepAIC(modelo_inicial, trace=TRUE, direction="both", scope=scope, k=2)## Start: AIC=55.37
## y ~ 1
##
## Df Sum of Sq RSS AIC
## + x4 1 96.615 58.173 20.182
## + x2 1 77.442 77.347 31.007
## + x3 1 46.603 108.186 43.758
## <none> 154.788 55.370
## + x5 1 0.343 154.446 57.286
## + x1 1 0.125 154.663 57.339
##
## Step: AIC=20.18
## y ~ x4
##
## Df Sum of Sq RSS AIC
## + x5 1 5.717 52.456 18.251
## + x2 1 5.321 52.852 18.537
## <none> 58.173 20.182
## + x1 1 1.429 56.745 21.237
## + x3 1 0.380 57.793 21.933
## - x4 1 96.615 154.788 55.370
##
## Step: AIC=18.25
## y ~ x4 + x5
##
## Df Sum of Sq RSS AIC
## + x2 1 6.603 45.853 15.139
## + x1 1 2.942 49.514 18.058
## <none> 52.456 18.251
## + x3 1 0.536 51.920 19.861
## - x5 1 5.717 58.173 20.182
## - x4 1 101.990 154.446 57.286
##
## Step: AIC=15.14
## y ~ x4 + x5 + x2
##
## Df Sum of Sq RSS AIC
## <none> 45.853 15.139
## + x1 1 1.6936 44.160 15.709
## + x3 1 0.1477 45.706 17.016
## - x2 1 6.6026 52.456 18.251
## - x5 1 6.9989 52.852 18.537
## - x4 1 25.6888 71.542 30.043modelo_both$anova ## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## y ~ 1
##
## Final Model:
## y ~ x4 + x5 + x2
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 37 154.78842 55.36997
## 2 + x4 1 96.614981 36 58.17344 20.18202
## 3 + x5 1 5.717446 35 52.45599 18.25076
## 4 + x2 1 6.602579 34 45.85341 15.13881m13<-lm(y~x4+x5+x2)
AIC(m13)## [1] 124.9781summary(m13)##
## Call:
## lm(formula = y ~ x4 + x5 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5707 -0.6256 0.1521 0.6467 1.7741
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.4672 1.3328 4.852 2.67e-05 ***
## x4 1.1997 0.2749 4.364 0.000113 ***
## x5 -0.6023 0.2644 -2.278 0.029127 *
## x2 0.5801 0.2622 2.213 0.033740 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.161 on 34 degrees of freedom
## Multiple R-squared: 0.7038, Adjusted R-squared: 0.6776
## F-statistic: 26.92 on 3 and 34 DF, p-value: 4.203e-09Mediante el summary se observa que todos los coeficientes son estadísticamente significativos, pues tienen un p-valor < 0.05.
\[ y = 6.467 + 0.5801x2 + 1.1997x4 - 0.6023x5\]
Por otro lado, aunque el modelo fue seleccionado mediante AIC, se procede a evaluar su calidad predictiva y parsimonia mediante el Cp de Mallows y el PRESS.
modelo_completo<-lm(y~x1+x2+x3+x4+x5)
cp <- ols_mallows_cp(m13, modelo_completo)
cp## [1] 3.92779El modelo final mostró un Cp cercano al número de parámetros (Cp=3.9 vs p=4), indicando que no omite variables relevantes ni sobreajusta.
PRESS<-PRESS(m13)
SSE <- sum(residuals(m13)^2)
resultado<-PRESS/SSE
resultado## [1] 1.222426El modelo demuestra buena capacidad predictiva, evidenciada por un ratio PRESS/SSE de 1.22, lo que indica que su desempeño en datos nuevos es comparable al ajuste en entrenamiento.
Finalmente, se valida que no haya multicolinealidad:
vifs <- vif(m13)
print(vifs)## x4 x5 x2
## 2.190770 1.044707 2.209831Todos los VIF están por debajo de 5, lo que indica que no hay problemas de multicolinealidad, todas las variables son independientes entre sí y aportan información única.
vifs <- vif(modelo_completo)
print(vifs)## x1 x2 x3 x4 x5
## 1.266390 2.381143 2.056492 2.682277 1.096731Todos los VIF están por debajo de 5, lo que indica que no hay problemas de multicolinealidad, todas las variables son independientes entre sí y aportan información única.
VERIFICACIÓN DE SUPUESTOS
Primero verifiquemos los residuos:
Residuos <- m13$residualsVeamos si son normales:
shapiro.test(Residuos)##
## Shapiro-Wilk normality test
##
## data: Residuos
## W = 0.96487, p-value = 0.2722Como el p-valor es mayor a 0.05, los residuos siguen una distribución aproximadamente normal.
plot(m13)
Podemos destacar varias cosas de los gráficos generados:
EL gráfico de “Residuals Vs Fitted” no se observa una forma curva o patrón, lo que sugiere que cumple con la suposición de linealidad e independencia de los errores. Por otro lado, se observa que la mayoría de sus datos están relativamente cerca del 0, sin embargo, las observaciones 20, 23 y 24 aparecen marcados, por lo que podrían ser valores atípicos.
El gráfico de “Q- Q Residuals” se destaca que la mayoría de los datos se encuentran sobre la línea teórica de la normalidad, lo que sugiere que los residuos son normales (comprobado por la prueba shapiro). Sin embargo, otra vez las observaciones 20, 23 y 24 se marcan y se encuentran en los extremos, lo cual puede estar desviando a la normalidad.
El gráfico de “Scale-Location” no se observa ningún patrón a simple vista, lo cual sugiere que cumple con la suposición de linealidad. No obstante, otra vez se destacan las observaciones 20, 23 y 24.
El gráfico de “Leverage” se observa que los puntos 9, 20 y 30 son de alto leverage.
COMPARACIÓN DE MODELOS
De acuerdo al modelo obtenido por los métodos Stepwise y entrada forzada vamos a compararlos utilizando la métrica del \(R^2\) ajustado y AIC:
De acuerdo con los resultados obtenidos al comparar el coeficiente de determinación ajustado (\(R^2\) ajustado) entre el modelo seleccionado por stepwise y el modelo completo, observamos que:
summary(m13)$adj.r.squared## [1] 0.677629summary(Modelo_Completo)$adj.r.squared## [1] 0.6769428Es decir:
\[R^2_{ajustado, stepwise} = 0.6776 > 0.6769 = R^2_{ajustado, completo}\]
Esto indica que el modelo seleccionado mediante el procedimiento stepwise explica ligeramente mejor la variabilidad de la variable respuesta que el modelo completo, utilizando menos variables.
Ahora, comparando por el AIC, observamos que:
AIC(m13)## [1] 124.9781AIC(Modelo_Completo)## [1] 126.7552Es decir:
\[AIC_{stepwise} = 124.97 < 126.75 = AIC_{completo}\]
En este caso, el modelo por stepwise presenta un AIC ligeramente menor, lo cual refleja un ajuste más eficiente y preciso al modelo, en comparación con el modelo completo.
Por último, realizaremos la comparación con la desviación estándar residual:
summary(m13)$sigma## [1] 1.161305summary(Modelo_Completo)$sigma## [1] 1.16254Es decir:
\[σ_{stepwise} = 1.1613 < 1.1625 = σ_{completo}\]
Esto implica que el modelo seleccionado por Stepwise presenta menores errores de predicción que el modelo completo, ya que tiene una desviación estándar residual ligeramente menor.
Como se puede observar, en todos las métricas el modelo hecho por el método Stepwise fue superior al modelo hecho por el método de entrada forzada.