Proyecto final - Exploratorio

Regresión Múltiple - IIMAS Logotipo

Author

Angel Escalante

library(tidyverse)
library(rworldmap)
library(corrplot)
library(classInt)
library(RColorBrewer)
library(janitor)
library(GGally)
library(here)

setwd(here("Proyecto"))

# Definiendo colores para consistencia
BLUE <- "#076fa2"
RED <- "#e63946"

datos <- read_delim(file = "BD03P.txt", delim = "|") |> clean_names()

Un vistazo a los datos

glimpse(datos)

Rows: 156
Columns: 9
$ puesto_ranking              <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,…
$ region                      <chr> "Finland", "Norway", "Denmark", "Iceland",…
$ puntaje_felicidad           <dbl> 7.632, 7.594, 7.555, 7.495, 7.487, 7.441, …
$ pib_percapita               <dbl> 1.305, 1.456, 1.351, 1.343, 1.420, 1.361, …
$ apoyo_social                <dbl> 1.592, 1.582, 1.590, 1.644, 1.549, 1.488, …
$ esperanza_de_vida_saludable <dbl> 0.874, 0.861, 0.868, 0.914, 0.927, 0.878, …
$ libertad                    <dbl> 0.681, 0.686, 0.683, 0.677, 0.660, 0.638, …
$ generosidad                 <dbl> 0.202, 0.286, 0.284, 0.353, 0.256, 0.333, …
$ percepcion_corrupcion       <dbl> 0.393, 0.340, 0.408, 0.138, 0.357, 0.295, …

datos <- datos |> select(-puesto_ranking)

Tenemos la siguiente información sobre los datos provista por nuestro querido profesor

Información de un estudio sobre el indice de felicidad en distintas regiones. Las variables se encuentran “estandarizadas” según el criterio del estudio, análogo a puntajes para cada variable.

Mapa sobe el puntaje de la felicidad 📍

mapa_datos <- joinCountryData2Map(
  datos, 
  joinCode = "NAME",
  nameJoinColumn = "region",
  nameCountryColumn = "region"
)

153 codes from your data successfully matched countries in the map
3 codes from your data failed to match with a country code in the map
90 codes from the map weren't represented in your data

par(mai = c(0, 0, 0.2, 0), xaxs = "i", yaxs = "i")

# Dividiendo los puntajes de felicidad en 7 categorías
class_int <- classIntervals(mapa_datos[["puntaje_felicidad"]], n = 7, style = "jenks")
cat_method <- class_int[["brks"]]

# Paleta de colores verde
color_palette <- brewer.pal(7, "Blues")

parametros_mapa <- mapCountryData(
  mapa_datos, 
  nameColumnToPlot = "puntaje_felicidad",
  addLegend = FALSE,
  catMethod = cat_method,
  colourPalette = color_palette
)

do.call(
  addMapLegend,
  c(
    parametros_mapa,
    legendLabels = "all",
    legendWidth = 0.5,
    legendIntervals = "data",
    legendMar = 2
  )
)

¿Cúales son los países más felices?

Vemos que los países occidentales tienden a “ser más felices”. Puede deberse a distintos factores, es difícil decir porque desconocemos la construcción del puntaje, por ahora, tenemos las siguiente dudas sobre la métrica y el estudio en general:

¿Fue una entrevista?
Si lo anterior es afirmativo, ¿cómo estuvo constituida la entrevista?
¿Se realizó a en diferentes momentos en el tiempo?
¿En qué año se realizó el estudio?

Sabemos que las diferencias culturales pueden impactar en los resultados pero sin conocer las condiciones del estudio es difícil hacer afirmaciones seguras. Por ahora, podemos limitarnos a analizar la información asumiendo el mejor esfuerzo y la mayor calidad de los datos.

Comportamiento de la métrica de felicidad 🙌

Densidad de la felicidad agregada y desagregada por continente

mapa_continente <- 
  tibble(
    continente = as.character(mapa_datos$REGION), 
    region = as.character(mapa_datos$region)
  ) |> drop_na(region)

datos <- datos |> left_join(mapa_continente, by = "region")

Hay tres países a los que no pudimos asignar un continente de forma programática, veámos cuáles son y asignemos uno manualmente

datos |> filter(is.na(continente)) |> pull(region)

[1] "Trinidad & Tobago"       "Northern Cyprus"        
[3] "Palestinian Territories"

Tenemos:

Trinidad y Tobago
Norte de Chipre
Territorios Palestinos

Por simplicidad, asignameros Trinidad y Tobago al continente Norteaméricano y al resto de países (Norte de Chipre y Territoruos Palestinos) al continente Asiático.

# Asignando continentes a los países faltantes
datos <- datos |> 
  mutate(
    continente = case_when(
      region == "Trinidad & Tobago" ~ "North America",
      region %in% c("Northern Cyprus", "Palestinian Territories") ~ "Asia",
      TRUE ~ continente
    )
  )

datos |>
  ggplot(aes(x = puntaje_felicidad)) + 
  geom_histogram(aes(y = ..density..), fill = BLUE, alpha = 0.5) +
  geom_density() + 
  theme_minimal() +
  geom_vline(xintercept = mean(datos$puntaje_felicidad), linetype = "dashed", color = RED) +
  labs(y = "", x = "Puntaje de felicidad")

datos |> 
  ggplot(aes(x = puntaje_felicidad)) + 
  geom_histogram(aes(y = ..density.., fill = as.factor(continente)), alpha = 0.4) +
  geom_density() +
  theme_minimal() +
  facet_wrap(~ continente, scales = "free_y") +
  labs(y = "", x = "Puntaje de felicidad") +
  guides(fill = "none")

Parece haber una tendencia de puntajes altos en los continentes que presentan países con mayor desarrollo capitalista o con tradición occidental, esto es, en los continentes Europeo, Norteamericano, Oceanía (Australia) y Sudamericano. Lo interesante es, Sudamerica no tiene un gran desarrollo capital por lo que valdría la pena explorar que otros componentes ayudan a explicar su métrica de felicidad.

Relación entre variables numéricas 🔎

Correlación

col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))

cor_datos <- datos |> select_if(is.numeric) |> cor()

test_significancia_cor <- function(mat, ...) {
    mat <- as.matrix(mat)
    n <- ncol(mat)
    p_mat <- matrix(NA, n, n)
    diag(p_mat) <- 0
    for (i in 1:(n - 1)) {
        for (j in (i + 1):n) {
            tmp <- cor.test(mat[, i], mat[, j], ...)
            p_mat[i, j] <- p_mat[j, i] <- tmp$p.value
        }
    }
  colnames(p_mat) <- rownames(p_mat) <- colnames(mat)
  p_mat
}

# Matriz del p-value de la correlation
p_mat <- datos |> select_if(is.numeric) |> test_significancia_cor()

corrplot(
  cor_datos, 
  method = "color",
  col = col(200),  
  type = "upper",
  order = "hclust", 
  addCoef.col = "black",
  tl.col = "black",
  tl.srt = 45,
  p.mat = p_mat,
  sig.level = 0.01,
  insig = "blank", 
  diag = FALSE 
)

Notamos en esta matrix de correlación que hay muchas relaciones fuertes y positivas entre las variables. Sin embargo, el coeficiente de correlación de la variable de generosidad no es significativa, y por tanto se excluyen sus valores de la gráfica.

Se puede apreciar que puntaje_felicidad tiene mayor correlación con pib_percapita y esperanza_de_ _vida_saludable (con \(\%80\) y \(\%78\) respectivamente), esto sugiere la relación entre el puntaje_felicidad y los países desarrollados que aparecen en el Top 20.

ajuste_fun <- function(data, mapping, ...) {
  p <- ggplot(data = data, mapping = mapping) + 
    geom_point(alpha = 0.3) + 
    geom_smooth(method = loess, fill = RED, color = RED, se = FALSE,...) +
    geom_smooth(method = lm, fill = BLUE, color = BLUE, se = FALSE, ...)
  p
}

datos |> 
  select_if(is.numeric) |> 
  ggpairs(lower = list(continuous = ajuste_fun)) +
  theme_minimal()

Si analizamos la columna de puntaje_felicidad, nuestra variable dependiende, vemos que hay una relación lineal positiva con pib_percapita, apoyo_social, esperanza_de_vida_saludable y libertad, mientras que no hay mucha claridad con generosidad y percepcion_corrupcion. En esta última, sobre todo, parece existir una relación no lineal, observamos una clara curva en el ajuste de la línea roja, el cuál es un ajuste LOESS (Local Polynomial Regression). Ver esta referencia para más detalles.

Explorando relación de percepción corrupción con la felicidad

datos |> 
  ggplot(aes(x = percepcion_corrupcion, y = puntaje_felicidad)) +
  geom_point(aes(color = continente)) +
  geom_smooth(color = BLUE, se = FALSE) +
  geom_smooth(method = "lm", color = RED, se = FALSE) +
  labs(x = "Percepción corrupción", y = "Puntaje felicidad") +
  theme_minimal()

Desconocemos el sentido de la variable percepcion_corrupcion, ¿a mayor valor, es menor la sensación de corrupción? Nos gustaría asumir que así es, dada la relación positiva que existe entre estas variables.

Un primer modelo

Ajustamos un modelo con la variable de Puntaje de felicidad como variable dependiente, y las siguientes variables como variables independientes (predictoras):

PIB Percapita
Apoyo Social
Esperanza de vida saludable
Percepción de corrupción
Libertad

\[ y_{felicidad}=\beta_0+\beta_{pib\_percapita} + \beta_{apoyo\_social} + \beta_{esperanza\_vida} + \beta_{percepcion\_corrupcion} \]

Sólo generosidad es omitida por no tener una relación clara con el puntaje de felicidad.

modelo_base <- lm(puntaje_felicidad ~ percepcion_corrupcion + pib_percapita + apoyo_social + esperanza_de_vida_saludable + libertad, data = datos)

summary(modelo_base)


Call:
lm(formula = puntaje_felicidad ~ percepcion_corrupcion + pib_percapita + 
    apoyo_social + esperanza_de_vida_saludable + libertad, data = datos)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.85879 -0.29406  0.00114  0.36971  1.07259 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                   1.9581     0.1844  10.621  < 2e-16 ***
percepcion_corrupcion         0.9184     0.5071   1.811   0.0721 .  
pib_percapita                 1.0565     0.2102   5.025 1.41e-06 ***
apoyo_social                  1.0128     0.2024   5.004 1.56e-06 ***
esperanza_de_vida_saludable   0.8228     0.3300   2.493   0.0137 *  
libertad                      1.4355     0.3132   4.584 9.55e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5246 on 150 degrees of freedom
Multiple R-squared:  0.7875,    Adjusted R-squared:  0.7804 
F-statistic: 111.2 on 5 and 150 DF,  p-value: < 2.2e-16

En este modelo base, tenemos un coeficiente \(R^2\) ajustado de 78.04%. Con todos sus coeficientes positivos y significativos a nivel \(\alpha =0.05\), excepto por \(\beta_{percepcion\_corrupcion}\) el cuál, como lo habíamos sospechado antes, puede tener una relación no lineal.

modelo_ajuste_tbl <- tibble(
  puntaje_felicidad = modelo_base$model$puntaje_felicidad,
  fitted = modelo_base$fitted.values,
  error = puntaje_felicidad - fitted
)

ggplot(modelo_ajuste_tbl, aes(x = fitted, y = puntaje_felicidad)) +
  geom_point() + 
  geom_abline(color = RED) +
  labs(x = "Valores predichos", y = "Puntaje felicidad") +
  theme_minimal()

ggplot(modelo_ajuste_tbl, aes(x = error)) +
  geom_histogram(aes(y = ..density..)) + 
  theme_minimal() +
  geom_density(color = RED)

Vemos que la densidad en el gráfico de los errores tienen cierto sesgo hacia la izquierda, lo cuál podría significar que el modelo sobreestima para algunos países. Sin embargo, los errores tienen una forma de campana, con un máximo en 0.

ggplot(modelo_ajuste_tbl, aes(sample = error)) +
  stat_qq(size = 2.5, color = RED, alpha = 0.6) +
  stat_qq_line() +
  theme_minimal() +
  labs(title = "Q-Q Plot")

Este mismo fenómeno se observa en el Q-Q Plot donde para los valores mayores a 2 y menores a -2 distan de los cuantiles teóricos.

Con este modelo tiene un buen ajuste, sin embargo, tenemos opciones para explorar: - Hay multicolinealidad entre variables - El modelo no ajusta muy bien para valores extremos - Podríamos probar relaciones no lineales para variables como percepcion_corrupcion