Ejemplo de RLM
Ángel R.
MnPyR 2020-2
El Problema del Vino
Se cree que la calidad del vino Pinot Noir se relaciona con algunas propiedades. El dataset “vino.xlsx” contiene información de dicho vino.
Ajuste un modelo de regresión lineal múltiple que relacione la calidad del vino con esas variables regresoras.
calidad aroma sabor cuerpo
1 9.8 3.3 3.1 2.8
2 12.6 4.4 3.5 4.9
3 11.9 3.9 4.8 5.3
4 11.1 3.9 3.1 2.6
5 13.3 5.6 5.5 5.1
6 12.8 4.6 5.0 4.7
Solución
Ajustaremos un modelo de RLM considerando lo siguiente:
Utilizaremos el método de todos los posibles modelos.
El criterio de selección del mejor modelo estará basado en \( R^2 \) ajustada.
Primer Paso: Relación entre Variables (Indicios)
Podemos hacer gráficas de dispersión para ver el comportamiento de los datos.
La gráfica anterior muestra que:
Las variables regresoras: aroma y sabor, tienen una correlación alta con la variable respuesta: calidad.
Posiblemente tengamos problemas de colinealidad, ya que, por ejemplo, cuerpo y sabor presentan correlación medianamente alta.
Posiblemente tengamos una o dos variables regresoras en el modelo para poder explicar la variable respuesta.
Segundo Paso: Generación de Modelos
El método de todos los posibles modelos consiste en generar \( 2^K \) modelos posibles y establecer, bajo algún criterio(o criterios), el mejor que pueda explicar la variable respuesta \( y \). Estamos asumiendo de antemano que todos los modelos que podamos generar incluyen el intercepto (\( \beta_0 \))
En este caso \( K=3 \) (3 regresores), entonces debemos generar \( 2^3=8 \) modelos, es decir:
\( y=\beta_0 \)
\( y=\beta_0 + \beta_1 x_1 \)
\( y=\beta_0 + \beta_2 x_2 \)
\( y=\beta_0 + \beta_3 x_3 \)
\( y=\beta_0 + \beta_1 x_1 +\beta_2 x_2 \)
\( y=\beta_0 + \beta_1 x_1 +\beta_3 x_3 \)
\( y=\beta_0 + \beta_2 x_2 +\beta_3 x_3 \)
\( y=\beta_0 + \beta_1 x_1 + \beta_2 x_2 +\beta_3 x_3 \)
Teóricamente tenemos 8 modelos, no obstante, haremos el ajuste de 7 modelos pues es evidente que \( y=\beta_0 \) no tiene mucho sentido.
Ajuste de los modelos
Hacemos el ajuste de los 7 modelos:
modelo1=lm(calidad~aroma)
modelo2=lm(calidad~sabor)
modelo3=lm(calidad~cuerpo)
modelo4=lm(calidad~aroma+sabor)
modelo5=lm(calidad~aroma+cuerpo)
modelo6=lm(calidad~sabor+cuerpo)
modelo7=lm(calidad~aroma+sabor+cuerpo)
mod1=summary(modelo1)
mod2=summary(modelo2)
mod3=summary(modelo3)
mod4=summary(modelo4)
mod5=summary(modelo5)
mod6=summary(modelo6)
mod7=summary(modelo7)
Algunos resultados que podemos observar son:
Call:
lm(formula = calidad ~ cuerpo)
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 ***
cuerpo 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.0003612
Call:
lm(formula = calidad ~ sabor + cuerpo)
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 ***
sabor 1.4883 0.2694 5.524 3.27e-06 ***
cuerpo 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-08
Call:
lm(formula = calidad ~ aroma + sabor + cuerpo)
Residuals:
Min 1Q Median 3Q Max
-2.26940 -0.55519 -0.03745 0.70825 2.45892
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.18465 1.22966 3.403 0.00172 **
aroma 0.50855 0.28254 1.800 0.08075 .
sabor 1.13707 0.32602 3.488 0.00137 **
cuerpo 0.07793 0.32906 0.237 0.81422
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.246 on 34 degrees of freedom
Multiple R-squared: 0.6591, Adjusted R-squared: 0.629
F-statistic: 21.91 on 3 and 34 DF, p-value: 4.436e-08
Tercer Paso: Elección del Mejor Modelo
La elección del mejor modelo será bajo el criterio de \( R^2 \) ajustada. Diremos que el mejor modelo es aquel que tiene \( R^2 \) ajustada más grande.
El mejor modelo es el modelo 4:
Call:
lm(formula = calidad ~ aroma + sabor)
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 ***
aroma 0.5180 0.2759 1.877 0.068849 .
sabor 1.1702 0.2905 4.027 0.000288 ***
---
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-09
Consideraciones
1.- El mejor modelo (modelo 4) es el que considera las variables regresoras: aroma ysabor, no obstante, la prueba parcial \( t \) muestra que la variable aroma podría ser no significativa (específicamente si \( \alpha=5 \)%)
2.- El \( R^2 \) ajustada del modelo es del 63.9%
3.- Haremos la validación de supuestos y veremos si es necesario modificar el modelo.
Cuarto Paso: Validación de Supuestos
Relación lineal de los regresores y la variable respuesta
La relación lineal se puede observar haciendo gráficas de los residuales:
Todos los supuestos sobre los Residuales
Pruebas Formales
Anderson-Darling normality test
data: rstandard(modelo4)
A = 0.23232, p-value = 0.7852
studentized Breusch-Pagan test
data: modelo4
BP = 4.1476, df = 2, p-value = 0.1257
Durbin-Watson test
data: modelo4
DW = 0.86864, p-value = 5.466e-05
alternative hypothesis: true autocorrelation is greater than 0
Multicolinealidad
Para ver si hay multicolinealidad en las variables regresoras que estamos considerando, podemos calcular el \( VIF \):
library("car")
vif(modelo4)
aroma sabor
2.185899 2.185899
De acuerdo a la literatura, si los \( VIF \) de las variables regresoras es mayor a 5 o 10, hay indicio de que los coeficientes asociados a la regresión están mal estimados debido a la multicolinealidad.
También podemos calcular el número de condición (coeficiente \( \kappa \)):
X1=scale(datos[,c(-1)])
A=t(X1)%*%X1
kappa=max(eigen(A)$values)/min(eigen(A)$values)
kappa
[1] 9.312234
La literatura indica que, en general, si el número de condición es menor que 100 entonces no hay problemas graves de multicolinealidad; si el número de condición oscila entre 100 a 1000 implican multicolinealidad de moderada a fuerte; si \( \kappa>1000 \) entonces hay indicio de una fuerte multicolinealidad.
Conclusiones Parciales
¿Cuáles son las conclusiones parciales?
Quinto Paso: Modificaciones del Modelo(Opcional)
Habíamos dicho que la variable regresora aroma podría ser no significativa a cierta \( \alpha \). Lo que haremos es quitar dicha variable y ver qué pasa.
El modelo sin la variable aroma
Si quitamos la variable aroma entonces tenemos el modelo \( calidad=\beta_0+ \beta_2 *sabor \) , el cual corresponde al modelo 2.
El modelo 2 presenta un \( R^2 \) ajustada del 61.37%
Pruebas Formales
Anderson-Darling normality test
data: rstandard(modelo2)
A = 0.24184, p-value = 0.7545
studentized Breusch-Pagan test
data: modelo2
BP = 2.2252, df = 1, p-value = 0.1358
Durbin-Watson test
data: modelo2
DW = 0.93935, p-value = 0.0001516
alternative hypothesis: true autocorrelation is greater than 0
Conclusiones Respecto al Problema
¿Cuáles son las conclusiones respecto al problema?
Conclusiones Respecto a la Metodología Empleada
El método de todos los posibles modelos es computacionalmente costoso. Si \( K \) es grande, el número de modelos posibles crece muy rápido.
La \( R^2 \) ajustada puede ser una buena estadística para comparar modelos en términos de variabilidad explicada. Lo recomendable es usarlo junto con otros criterios.
Existen otros criterios de comparación y métodos de generación de modelos.