Objetivo 3: Regresión
Elena Castán
2023-04-26
Librerías que usaremos:
library(DataExplorer)
library(dplyr)
library(highcharter)
library(treemap)
library(viridis)
library(ggplot2)
library(treemapify)
library(viridisLite)
library(knitr)
library(gridExtra)
library(corrplot)Observar base de datos
- Cargamos los datos que nos ha proporcionado el objetivo 1. En este archivo tenemos 4 dataframes, usaremos la base de datos original con todas las variables imputadas.
load("base_objetivo1.RData")Incluimos las siguientes variables en el modelo de regresión:
Salario Medio
Índice felicidad
Gasto en educación
Gasto en salud
Consumo eléctrico
PIB
Son las mismas variables que hemos introducido en el ANOVA multifactorial. Hemos elegido meterlas todas independientemente de si su influencia ha resultado ser significativa o no. De esta manera compararemos las variables de nuestro modelo de regresión final y veremos si ambos estudios son coherentes.
Creamos nuestro nuevo dataframe con las variables elegidas. También hemos decidido escalar las variables para que tengan media 0 y desviación típica 1. El resultado de la predicción del coste de vida será el mismo que si no se escalaran previamente, únicamente cambia el valor de los regresores, especialmente de \(b_{0}\), que será nulo.
Como el objetivo principal de esta regresión no es la predicción sino la interpretación de resultados para cuantificar la influencias de las variables regresoras a la variable respuesta (coste de vida), nos interesa mucho escalar los datos. De esta manera el valor numérico de los regresores no se verá afectado por la magnitud de su variable asociada. Si no escalaramos los regresores de las variables de magnitudes grandes como el PIB serían más pequeños y nos daría una idea errónea de su efecto en el coste de vida.
data <- select(datos_imputados, Indice_de_Costo_de_Vida, Salario_Medio, Indice_Felicidad, Gasto_Publico, Gasto_Educacion, Gasto_Salud, Consumo_GWh, PIB)
data <- data.frame(scale(data))
head(data)## Indice_de_Costo_de_Vida Salario_Medio Indice_Felicidad Gasto_Publico
## 1 3.866864 4.0155214 -0.1093252 -0.01367741
## 2 2.825705 2.3748759 -0.1077708 -0.27771714
## 3 2.778772 1.8266465 -0.1011088 -0.04958975
## 4 2.315375 -0.8898684 3.1683778 -0.28705829
## 5 1.974975 -0.8880796 -0.1526282 -0.25541038
## 6 1.826084 1.1812483 -0.2863123 1.90799784
## Gasto_Educacion Gasto_Salud Consumo_GWh PIB
## 1 -0.2905838 -0.2878116 0.29909789 2.9247839
## 2 -0.2887913 -0.2838930 7.70613120 2.2168813
## 3 -0.2904777 -0.2821163 3.06427294 2.5783849
## 4 -0.2951020 -0.2876207 0.09314938 0.6993215
## 5 -0.2978490 -0.2879943 0.96283872 4.1940744
## 6 -0.2981672 -0.2773342 0.48334230 0.9645584Regresión
Modelo inicial
Metemos todas las variables en un modelo de regresión lineal. Usaremos el método mixto de regresión para generar nuestro primer modelo.
modelo <- lm(Indice_de_Costo_de_Vida ~ Salario_Medio + Indice_Felicidad + Gasto_Educacion + Gasto_Salud + Consumo_GWh + PIB, data = data )
summary(modelo)##
## Call:
## lm(formula = Indice_de_Costo_de_Vida ~ Salario_Medio + Indice_Felicidad +
## Gasto_Educacion + Gasto_Salud + Consumo_GWh + PIB, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.02650 -0.35668 -0.05786 0.33922 1.84744
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.998e-16 4.641e-02 0.000 1.00000
## Salario_Medio 2.566e-01 7.862e-02 3.264 0.00145 **
## Indice_Felicidad 9.204e-02 5.195e-02 1.772 0.07914 .
## Gasto_Educacion -2.858e-02 5.672e-02 -0.504 0.61534
## Gasto_Salud 1.186e-01 6.159e-02 1.925 0.05675 .
## Consumo_GWh 4.895e-02 6.380e-02 0.767 0.44451
## PIB 6.089e-01 8.473e-02 7.186 7.65e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5084 on 113 degrees of freedom
## Multiple R-squared: 0.7546, Adjusted R-squared: 0.7415
## F-statistic: 57.9 on 6 and 113 DF, p-value: < 2.2e-16Tenemos una R cuadrada alta, de 0.7546. Esto significa que nuestro modelo de regresión es capaz de explicar el 75% de la variabilidad del coste de vida. El p-valor del modelo es significativo lo que nos asegura que el modelo no es por azar, pero hay algunos p-valores no significativos lo que puede significar que algunas variables no están contribuyendo al modelo y haciéndolo más complejo sin necesidad.
Cálculo de posibles modelos
A continuación haremos varios modelos de regresión usando el método del stepwise mixto para elegir nuestras variables.
step(object = modelo, direction = "both", trace = 1)## Start: AIC=-155.57
## Indice_de_Costo_de_Vida ~ Salario_Medio + Indice_Felicidad +
## Gasto_Educacion + Gasto_Salud + Consumo_GWh + PIB
##
## Df Sum of Sq RSS AIC
## - Gasto_Educacion 1 0.0656 29.273 -157.30
## - Consumo_GWh 1 0.1522 29.359 -156.95
## <none> 29.207 -155.57
## - Indice_Felicidad 1 0.8113 30.018 -154.28
## - Gasto_Salud 1 0.9577 30.165 -153.70
## - Salario_Medio 1 2.7542 31.961 -146.76
## - PIB 1 13.3461 42.553 -112.41
##
## Step: AIC=-157.3
## Indice_de_Costo_de_Vida ~ Salario_Medio + Indice_Felicidad +
## Gasto_Salud + Consumo_GWh + PIB
##
## Df Sum of Sq RSS AIC
## - Consumo_GWh 1 0.1459 29.419 -158.70
## <none> 29.273 -157.30
## - Indice_Felicidad 1 0.8671 30.140 -155.80
## + Gasto_Educacion 1 0.0656 29.207 -155.57
## - Gasto_Salud 1 0.9841 30.257 -155.33
## - Salario_Medio 1 2.8305 32.103 -148.22
## - PIB 1 13.4506 42.723 -113.93
##
## Step: AIC=-158.7
## Indice_de_Costo_de_Vida ~ Salario_Medio + Indice_Felicidad +
## Gasto_Salud + PIB
##
## Df Sum of Sq RSS AIC
## <none> 29.419 -158.70
## + Consumo_GWh 1 0.1459 29.273 -157.30
## - Gasto_Salud 1 0.8828 30.302 -157.16
## + Gasto_Educacion 1 0.0594 29.359 -156.95
## - Indice_Felicidad 1 1.0262 30.445 -156.59
## - Salario_Medio 1 3.1268 32.545 -148.58
## - PIB 1 16.5195 45.938 -107.22##
## Call:
## lm(formula = Indice_de_Costo_de_Vida ~ Salario_Medio + Indice_Felicidad +
## Gasto_Salud + PIB, data = data)
##
## Coefficients:
## (Intercept) Salario_Medio Indice_Felicidad Gasto_Salud
## 2.140e-16 2.691e-01 1.014e-01 9.535e-02
## PIB
## 6.332e-01Modelo final
El mejor modelo de regresión ha resultado ser el siguiente:
Indice_de_Costo_de_Vida ~ Salario_Medio + Indice_Felicidad + Gasto_Salud + PIBHa eliminado algunas variables, seguramente porque su contribución al modelo no era suficiente. Las variables que el algoritmo ha incluido en el modelo final coinciden completamente con las variables cuya influencia ha resultado significativa en el ANOVA. Ambos estudios normalmente van de la manos pero al ser algoritmos diferentes suele haber alguna que otra discrepancia entre ellos pero en nuestro caso son completamente coherentes.
Creamos el modelo final:
modelo <- lm(Indice_de_Costo_de_Vida ~ Salario_Medio + Indice_Felicidad + Gasto_Salud + PIB, data = data )
summary(modelo)##
## Call:
## lm(formula = Indice_de_Costo_de_Vida ~ Salario_Medio + Indice_Felicidad +
## Gasto_Salud + PIB, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.07973 -0.33100 -0.05676 0.32866 1.81812
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.140e-16 4.617e-02 0.000 1.000000
## Salario_Medio 2.691e-01 7.697e-02 3.496 0.000672 ***
## Indice_Felicidad 1.014e-01 5.064e-02 2.003 0.047541 *
## Gasto_Salud 9.535e-02 5.133e-02 1.858 0.065768 .
## PIB 6.332e-01 7.879e-02 8.036 9.06e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5058 on 115 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7442
## F-statistic: 87.55 on 4 and 115 DF, p-value: < 2.2e-16Volvemos a tener una R cuadrada alta, de 0.7528. Esto significa que nuestro modelo de regresión es capaz de explicar el 75% de la variabilidad del coste de vida. Porque la R cuadrada sea menor que el modelo inicial no significa que sea peor, para comparar el rendimiento de modelos con cantidad de regresores diferentes se utiliza la R cuadrada corregida que en este modelo es mayor que en la anterior (0.7442 opuesto a 0.7415). Además, en este caso concreto hemos decidido sacrificar unas centésimas para deshacernos de la complejidad innecesaria de modelo inicial.
Para una mayor información calculamos los intervalos de confianza de los coeficientes regresores pero nos basaremos sobre todo en su valor numérico a la hora de interpretarlos.
confint(modelo)## 2.5 % 97.5 %
## (Intercept) -0.091456488 0.09145649
## Salario_Medio 0.116633457 0.42155992
## Indice_Felicidad 0.001118490 0.20173158
## Gasto_Salud -0.006318248 0.19702463
## PIB 0.477092838 0.78923593Interpretación
Una de las finalidades de este objetivo es cuantificar cómo condicionan las variables elegidas el coste de vida de un país. Para ello usaremos los coeficientes regresores de nuestro modelo final. Creamos un dataframe con los valores que nos interesan: los coeficientes del modelo.
coeff <- data.frame(coefficients = modelo$coefficients)
coeff$variables <- rownames(coeff)
kable(coeff[1])| coefficients | |
|---|---|
| (Intercept) | 0.0000000 |
| Salario_Medio | 0.2690967 |
| Indice_Felicidad | 0.1014250 |
| Gasto_Salud | 0.0953532 |
| PIB | 0.6331644 |
Cuanto mayor sea el coeficiente de regresión, mayor será la contribución de su variable asociada al coste de vida (nuestra variable respuesta) y mayor por lo tanto la dependencia del coste de vida en esa variable. Para mayor comprensión se podría hablar de “importancia” de la variable (p.e. el gasto en salud es más importante que el salario medio) pero sería erróneo pues una variable no es más importante que la otra, simplemente tiene mayor influencia.
A pesar de la interpretación anterior de los coeficientes numéricos (más exacta) vamos a observar unos gráficos que nos permitan hacernos una idea más a simple vista de las influencias.
hchart(coeff, "column", hcaes(x = variables, y = coefficients))Para comparar las variables entre sí quizás es más útil el siguiente gráfico de áreas en el que se mostrarán más grandes las variables con mayor influencia.
ggplot(coeff, aes(area = coefficients, fill = variables, label = variables)) +
geom_treemap() +
geom_treemap_text(color = "white",
place = "centre",
size = 15)
Hemos tomado la decisión de usar gráficos interactivos en nuestro vídeo de presentación del proyecto. Repetimos el gráfico anterior pero con la característica de la interactividad.
hchart(coeff, "treemap", hcaes(x = variables, value = coefficients, color = coefficients)) %>%
hc_colorAxis(stops = color_stops(colors = viridis::inferno(4)))Gráficos de residuos
Ahora vemos la forma que adoptan los residuos de nuestro modelo final.
Indice_de_Costo_de_Vida ~ Salario_Medio + Indice_Felicidad + Gasto_Salud + PIB
plot1 <- ggplot(data = data, aes(Salario_Medio , modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot2 <- ggplot(data = data, aes(Indice_Felicidad , modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot3 <- ggplot(data = data, aes(Gasto_Salud, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot4 <- ggplot(data = data, aes(Consumo_GWh, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot1, plot2, plot3, plot4)## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Complementaremos esto con el gráfico de normalidad Q-Q :
qqnorm(modelo$residuals)
qqline(modelo$residuals)
Todas las observaciones se acoplan bastante bien a la línea, tenemos algún dato que se desvía arriba a la derecha pero no parece lo suficientemente extremo como para ser preocupante.
shapiro.test(modelo$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo$residuals
## W = 0.97373, p-value = 0.01886El test de Shapiro-Wilk confirma la normalidad de los residuos usando un nivel de confianza de 99% (p-valor mayor que 0.01).
Por último, en el siguiente gráfico controlamos que no hayan correlaciones altas entre las variables regresoras que hayan podido afectar en el cálculo del modelo.
corrplot(cor(dplyr::select(data, Salario_Medio, Indice_Felicidad, Gasto_Salud, Consumo_GWh)), method = "number", tl.col = "black")
Observamos una pequeña correlación entre el salario medio y el consumo eléctrico. Aun así no es lo suficientemente alta como para comprometer los resultados de la regresión, más si tenemos en cuenta que el consumo eléctrico no se encuentra como variable regresora en el modelo final.