Objetivo 3: ANOVA
Elena Castán
2023-04-21
Librerías usadas:
library(dplyr)
library(car)
library(onewaytests)
library(knitr)
library(ggplot2)
library(highcharter)Lectura y transformación de datos
Hemos decidido realizar un ANOVA multifactorial con estas variables. Para ello usaremos estas variables discretizadas para que sirvan como factores y el Índice de coste de vida en su versión numérica puesto que es la variable en la nos queremos focalizar.
load("base_objetivo1.RData")
futuros_factores <- select(datos_imputados, Pais, Salario_Medio, Salario_Minimo, Indice_Felicidad, Gasto_Educacion, Gasto_Salud, Consumo_GWh, PIB)A continuación creamos una función que normaliza las variables y devuelve sus percentiles.
Ahora aplicamos esta función a cada variable de nuestra base de datos que queremos que cuente como factor para nuestro ANOVA.
for (nom_col in colnames(futuros_factores)[2:8]){
futuros_factores[[nom_col]] = percentiles(futuros_factores[[nom_col]])
}
head(futuros_factores)## Pais Salario_Medio Salario_Minimo Indice_Felicidad Gasto_Educacion
## 1 Suiza 0.3500 0.6000 0.6167 0.7917
## 2 Islandia 0.5833 0.6083 0.9917 0.8667
## 3 Noruega 0.6167 0.6750 0.7750 0.2750
## 4 Bahamas 0.4083 0.7000 0.3667 0.2833
## 5 Luxemburgo 0.4917 0.7500 0.9833 0.9250
## 6 Japón 0.9917 0.7750 0.9583 0.2583
## Gasto_Salud Consumo_GWh PIB
## 1 0.9750 0.7583 0.3500
## 2 0.9250 0.6750 0.6750
## 3 0.8500 0.3500 0.7583
## 4 0.9833 0.9667 0.5333
## 5 1.0000 0.5333 0.9667
## 6 0.9167 0.4333 0.4333Finalmente, una vez hechas estas transformaciones discretizamos en factores de 4 niveles. Hemos elegido disminuir el número de niveles en comparación al objetivo 1 para facilitar la comprensión y la interpretación del ANOVA. Al haber hecho estas transformaciones antes de discretizar nos aseguramos de que el número de observaciones de cada factor es parecido si no el mismo. Esto nos garantiza que los resultados del ANOVA más fiables porque esta implementación en R de la herramienta tiene problemas a la hora de aplicarlo a clases desequilibradas.
datos_factores <- futuros_factores
for (nom_var in colnames(futuros_factores)[2:8]){
a = cut(futuros_factores[[nom_var]], breaks = c(0, 0.25, 0.5, 0.75, 1))
levels(a) <- c("Bajo", "Medio-Bajo", "Medio-Alto", "Alto")
datos_factores[[nom_var]] = a
}
head(datos_factores)## Pais Salario_Medio Salario_Minimo Indice_Felicidad Gasto_Educacion
## 1 Suiza Medio-Bajo Medio-Alto Medio-Alto Alto
## 2 Islandia Medio-Alto Medio-Alto Alto Alto
## 3 Noruega Medio-Alto Medio-Alto Alto Medio-Bajo
## 4 Bahamas Medio-Bajo Medio-Alto Medio-Bajo Medio-Bajo
## 5 Luxemburgo Medio-Bajo Medio-Alto Alto Alto
## 6 Japón Alto Alto Alto Medio-Bajo
## Gasto_Salud Consumo_GWh PIB
## 1 Alto Alto Medio-Bajo
## 2 Alto Medio-Alto Medio-Alto
## 3 Alto Medio-Bajo Alto
## 4 Alto Alto Medio-Alto
## 5 Alto Medio-Alto Alto
## 6 Alto Medio-Bajo Medio-BajoAhora concatenamos nuestros factores con las variables de coste de vida. En principion solamente usaremos Índice de coste de vida pero añadiremos también el índice que tiene en cuenta el alquiler por si nos sobra tiempo y queremos desarrollar más el análisis.
variables_respuesta <- select(datos, Pais, Indice_de_Costo_de_Vida, Indice_de_Costo_de_Vida_mas_Alquiler)
datos_ANOVA <- merge(variables_respuesta, datos_factores, by = "Pais")ANOVA
Supuestos básicos
Antes comenzar comprobaremos las tres hipótesis básicas para realizar un ANOVA: la hipótesis de independencia, la hipótesis de normalidad y la hipótesis de homocedasticidad. Usaremos para todas las pruebas una confianza del 99%.
Se supone independencia de las observaciones cuando se aplica muestreo aleatorio.
Comprobamos la hipótesis de normalidad con el test de Shapiro-Wilk.
shapiro.test(datos_ANOVA$Indice_de_Costo_de_Vida)##
## Shapiro-Wilk normality test
##
## data: datos_ANOVA$Indice_de_Costo_de_Vida
## W = 0.92918, p-value = 7.292e-06El p-valor es menor que 0.01, por lo que se rechaza la hipótesis nula de normalidad. Haremos una transformación aplicando el logaritmo y volvemos a aplicar el test.
datos_ANOVA$Indice_de_Costo_de_Vida <- log(datos_ANOVA$Indice_de_Costo_de_Vida)
shapiro.test(datos_ANOVA$Indice_de_Costo_de_Vida)##
## Shapiro-Wilk normality test
##
## data: datos_ANOVA$Indice_de_Costo_de_Vida
## W = 0.98983, p-value = 0.5067Efectivamente, ahora el p-valor es mayor que 0.01 aceptando la hipótesis de normalidad de los datos. Procedemos a comprobar la homocedasticidad de la variable respuesta dentro de los diferentes grupos. Para ello hemos decidido usar el test de Fligner-Killeen porque es menos sensible que el test de Bartlett a posibles asimetrías y nos dará un resultado más robusto.
p.valores = NULL
for (i in 4:10){
factor = colnames(datos_ANOVA)[i]
formula = as.formula(paste("Indice_de_Costo_de_Vida", "~", factor))
test <- fligner.test(formula, data=datos_ANOVA)
p.valores = c(p.valores, test$p.value)
}
names(p.valores) = colnames(datos_ANOVA)[4:10]
kable(p.valores, col.names = "p-valores")| p-valores | |
|---|---|
| Salario_Medio | 0.2075643 |
| Salario_Minimo | 0.0076603 |
| Indice_Felicidad | 0.0422756 |
| Gasto_Educacion | 0.6338567 |
| Gasto_Salud | 0.0746719 |
| Consumo_GWh | 0.1835042 |
| PIB | 0.7085429 |
Como bien se ve en la tabla, todos los p-valores son mayores a 0.01 excepto el salario mínimo, por lo que se acepta la hipótesis nula de homocedasticidad en la mayoría de factores.
Normalmente se meterían todas las variables, pues que no se cumpla la homocedastidad en alguna variable no imposibilita el realizar el ANOVA ni significa que sus resultados serán erróneos, solamente que hay que poner cuidado a la hora de interpretarlos. Sin embargo, previamente se había decidido meter solamente una de las dos variables de salario para evitar inexactitudes en el ANOVA por la alta correlación entre ambas. Lógicamente, eliminamos el salario mínimo que es la variable que no cumple con la hipótesis de homocedastidad.
Tenemos suerte, pues con el tratamiento previo de datos y los resultados de los tests anteriores podemos asumir unos resultados robustos y de la máxima fiabilidad. Esto no siempre ocurre y aunque no ocasiona problemas per se, puede dar lugar a inexactitudes que influyan en las conclusiones finales. Ahora tenemos la ventaja de poder operar en el ANOVA con comodidad y sin la inquietud de que hayan detalles que hayan pasado desapercibidos.
Aceptamos los tres supuestos básicos de un ANOVA: el supuesto de independencia, el supuesto de normalidad y el supuesto de homocedasticidad.
Modelo
Generamos el modelo de ANOVA con nuestros 6 factores y vemos cuáles son significativos y cuáles no.
anova <- aov(Indice_de_Costo_de_Vida ~ Salario_Medio*Indice_Felicidad*Gasto_Educacion*Gasto_Salud*Consumo_GWh*PIB, data = datos_ANOVA)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## Salario_Medio 3 5.632 1.8774 64.775 5.92e-06 ***
## Indice_Felicidad 3 1.646 0.5488 18.934 0.000543 ***
## Gasto_Educacion 3 0.254 0.0848 2.926 0.099835 .
## Gasto_Salud 3 1.447 0.4823 16.639 0.000845 ***
## Consumo_GWh 3 0.263 0.0878 3.029 0.093372 .
## PIB 3 0.628 0.2093 7.221 0.011525 *
## Salario_Medio:Indice_Felicidad 9 0.455 0.0506 1.745 0.222446
## Salario_Medio:Gasto_Educacion 9 0.443 0.0492 1.697 0.234007
## Indice_Felicidad:Gasto_Educacion 9 0.465 0.0516 1.781 0.214103
## Salario_Medio:Gasto_Salud 8 0.249 0.0311 1.074 0.460889
## Indice_Felicidad:Gasto_Salud 9 0.493 0.0547 1.889 0.191518
## Gasto_Educacion:Gasto_Salud 9 0.511 0.0568 1.960 0.177987
## Salario_Medio:Consumo_GWh 9 0.627 0.0696 2.403 0.115692
## Indice_Felicidad:Consumo_GWh 9 0.332 0.0369 1.273 0.372155
## Gasto_Educacion:Consumo_GWh 7 0.430 0.0615 2.121 0.156945
## Gasto_Salud:Consumo_GWh 9 0.881 0.0979 3.377 0.050415 .
## Salario_Medio:PIB 4 0.036 0.0091 0.313 0.861905
## Gasto_Educacion:PIB 1 0.007 0.0066 0.226 0.647047
## Gasto_Salud:PIB 1 0.000 0.0000 0.000 0.999912
## Consumo_GWh:PIB 2 0.170 0.0848 2.925 0.111345
## Residuals 8 0.232 0.0290
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interacción Gasto Salud - Consumo GWh
Investigamos más en profundidad la interacción significativa entre el gasto en salud y el consumo eléctrico.
ggplot(data = datos_ANOVA, aes(x = Consumo_GWh, y = Indice_de_Costo_de_Vida, colour = Gasto_Salud,
group = Gasto_Salud)) +
stat_summary(fun = mean, geom = "point") +
stat_summary(fun = mean, geom = "line") +
labs(y = 'mean (Indice_de_Costo_de_Vida)') +
theme_bw()
Box-plot
Investigamos las medias separadas por factores con gráficos boxplot.
ggplot(data = datos_ANOVA, aes(x = Salario_Medio,
y = Indice_de_Costo_de_Vida,
color = Salario_Medio)) +
geom_boxplot() +
theme_bw()
ggplot(data = datos_ANOVA, aes(x = Indice_Felicidad,
y = Indice_de_Costo_de_Vida,
color = Indice_Felicidad)) +
geom_boxplot() +
theme_bw()
ggplot(data = datos_ANOVA, aes(x = Gasto_Educacion,
y = Indice_de_Costo_de_Vida,
color = Gasto_Educacion)) +
geom_boxplot() +
theme_bw()
ggplot(data = datos_ANOVA, aes(x = Gasto_Salud,
y = Indice_de_Costo_de_Vida,
color = Gasto_Salud)) +
geom_boxplot() +
theme_bw()
ggplot(data = datos_ANOVA, aes(x = Consumo_GWh,
y = Indice_de_Costo_de_Vida,
color = Consumo_GWh)) +
geom_boxplot() +
theme_bw()
ggplot(data = datos_ANOVA, aes(x = PIB,
y = Indice_de_Costo_de_Vida,
color = PIB)) +
geom_boxplot() +
theme_bw()
Hemos tomado la decisión de usar gráficos interactivos en nuestro vídeo de presentación del proyecto. Repetimos los gráficos anteriores pero esta vez con la característica de la interactividad.
hcboxplot(x = datos_ANOVA$Indice_de_Costo_de_Vida, var = datos_ANOVA$Salario_Medio,
name = "Length", color = "#2980b9") %>%
hc_add_theme(hc_theme_bloom()) %>% hc_title(text = "Salario Medio")hcboxplot(x = datos_ANOVA$Indice_de_Costo_de_Vida, var = datos_ANOVA$Indice_Felicidad,
name = "Length", color = "#2980b9") %>%
hc_add_theme(hc_theme_bloom()) %>% hc_title(text = "Índice de Vida")hcboxplot(x = datos_ANOVA$Indice_de_Costo_de_Vida, var = datos_ANOVA$Gasto_Educacion,
name = "Length", color = "#2980b9") %>%
hc_add_theme(hc_theme_bloom()) %>% hc_title(text = "Gasto en Educación")hcboxplot(x = datos_ANOVA$Indice_de_Costo_de_Vida, var = datos_ANOVA$Gasto_Salud,
name = "Length", color = "#2980b9") %>%
hc_add_theme(hc_theme_bloom()) %>% hc_title(text = "Gasto en Salud")hcboxplot(x = datos_ANOVA$Indice_de_Costo_de_Vida, var = datos_ANOVA$Consumo_GWh,
name = "Length", color = "#2980b9") %>%
hc_add_theme(hc_theme_bloom()) %>% hc_title(text = "Consumo eléctrico")hcboxplot(x = datos_ANOVA$Indice_de_Costo_de_Vida, var = datos_ANOVA$PIB,
name = "Length", color = "#2980b9") %>%
hc_add_theme(hc_theme_bloom()) %>% hc_title(text = "PIB")Gráficos de residuos
Observamos los gráficos de residuos en busca de alguna irregularidad que pueda afectar a nuestro análisis.
par(mfrow = c(1,2))
plot(anova, which = 1:4)
Vemos 3 observaciones que destacan. Las observaciones 6, 91 y 66, corresponden a los países Pakistán, Armenia y Jordania. Pakistán destacaba ya en el objetivo 1 por ser el país más bajo respecto al coste de vida y el índice económico.