---
title: "Introducción a la regresión lineal múltiple (II)"
subtitle: "Ejercicio integrador"
output: 
  html_document:
    toc: TRUE
    toc_float: TRUE
author: "Germán Rosati"

date: "`r format(Sys.time(), '%d %B, %Y')`"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=FALSE)
library(tidyverse)
library(broom)
library(leaps)
library(car)
```

```{r message=FALSE, warning=FALSE, include=FALSE, paged.print=FALSE}
knitr::opts_chunk$set(echo=TRUE, message=FALSE, warning=FALSE, collapse=TRUE, highlight=TRUE, paged.print=FALSE, prompt=TRUE, strip.white=FALSE, tidy = TRUE)
```

## Introducción
Vamos a trabajar con dos set de datos. El primero, es un dataset (`.wb_bank_data_2019.csv`) con información construida por el Banco Mundial acerca de la distribución de la población ocupada según grandes sectores de actividad y categoría ocupacional:

- `SL.AGR.EMPL.ZS`: % total el empleo en agricultura: consistente en las actividades de agricultura, silvicultura, caza y pesca  
- `SL.IND.EMPL.ZS`: % total de empleo en industria:  en actividades de minas y canteras, manufacturas, construcción y energía, gas y agua.
- `SL.SERV.EMPL.ZS` % total de empleo el servicios comprende comercio al por mayor y menor y restaurantes y hoteles, transporte, almacenamiento y comunicaciones, finanzas, seguros, servicios inmobiliarios y a las empresas y servicios personales, sociales y comunales.
- `SL.FAM.WORK.ZS`: % total de trabajadores familiares
- `SL.EMP.MPYR.ZS`: % total de empleadores  
- `SL.EMP.SELF.ZS`: % total de independientes (cuenta propias)
- `SL.EMP.WORK.ZS`: % total de trabajadores asalariados


El segundo (`HDI_HDR2020_040722.csv`) son las componentes del Human Development Index (HDI) elaborado por UNDP para el período 1992-2019. El HDI es un promedio a nivel país de cuatro variables fundamentales:

![](./imgs/hdiRoadMap.png)

En el dataset va a encontrarlas llamadas de la siguiente forma:
- `le_2019`: expectativa de vida al nacer de la población
- `eys_2019`: años "esperados" de escolaridad de la población
- `mys_2019`: años promedios de escolaridad
- `gnipc_2019`: producto bruto interno per cápita
- `hdi_2019`: HDI para el país

En ambos datasets van a encontrar los códigos de país (`iso3c`) y las etiquetas de país (`country`)

## Preprocesamiento
### Consigna 1.
Cargue los dos datasets. ¿Cuál es la estructura de los datos? ¿Qué formato tienen? ¿Cuál es la unidad de análisis y cuáles las variables? ¿Cuál es el alcance temporal y geográfico?

```{r}
df <- read_csv('./data/wb_bank_data_2019.csv')
hdi <- read_csv('./data/HDI_HDR2020_040722.csv')

df <- df %>%
     pivot_wider(id_cols = c(iso3c, iso2c, country), names_from = indicatorID, values_from = value)
```

### Consigna 2.
Seleccione del dataset del HDI solamente las columnas que corresponden al año 2019.
```{r}
hdi <- hdi %>%
        select(iso3c, country, contains("_2019"))

```


### Consigna 3. 
Con el dataset recortado, realice un join para tener todas las variables en una sola tabla. ¿Qué variable debería usar para vincular ambas tablas?

```{r}
final <- hdi %>%
            left_join(
              df %>% select(-country)
              )

final <- final %>% drop_na()

```


## Análisis descriptivo
### Consigna 4. 
'¿Qué relacion existe (si es que existe) entre el HDI y la proporción de población asalariada? ¿Y entre el HDI y la proporción de población en el sector agropecuario?

```{r}
final %>%
  ggplot() + 
    geom_point(aes(x=SL.EMP.WORK.ZS, y=hdi_2019)) + 
    ylim(0,1) + 
    xlim(0,100) +
    theme_minimal()
```

```{r}
final %>%
  ggplot() + 
    geom_point(aes(x=SL.AGR.EMPL.ZS, y=hdi_2019)) + 
    ylim(0,1) + 
    xlim(0,100) +
    theme_minimal()
```

Correlación

```{r}
final %>%
  select(SL.AGR.EMPL.ZS,hdi_2019) %>%
  cor(., method='pearson')
```


## Modelo de regresión lineal múltiple
### Consigna 5
Queremos construir información sobre la relación entre el PBI per cápita (`gnipc_2019`) y variables relativas a la educación, a la salud y al grado de desarrollo capitalista del país. Entrenar una regresión lineal múltiple con las variables que considere relevantes. 

```{r}
lm_1 <- final %>%
          lm(gnipc_2019 ~ le_2019 + mys_2019 + SL.AGR.EMPL.ZS + SL.SRV.EMPL.ZS + SL.EMP.WORK.ZS, data = .)

lm_1 <- lm(gnipc_2019 ~ le_2019 + mys_2019 + SL.AGR.EMPL.ZS + SL.SRV.EMPL.ZS + SL.EMP.WORK.ZS, data = final)


summary(lm_1)
```


## Análisis de supuestos del modelo
### Consigna 6.
Chequeen los residuos del modelo. ¿Qué herramienta gráfica pueden usar para eso? ¿Son aleatorios?
```{r}
resids <- tibble(
  iso3c = final$iso3c,
  gnipc = final$gnipc_2019,
  gnipc_pred = predict(lm_1, final),
  resid =  gnipc - gnipc_pred
) 

resids %>%
  ggplot(aes(x=gnipc_pred, y=resid, color=resid)) + 
    geom_point() +
    scale_color_viridis_c() +
    geom_hline(yintercept = 0, linetype='dashed') + 
    theme_minimal()

```

### Consigna 7.
Realicen un histograma de los residuos. ¿Qué forma tienen?
```{r}
resids %>%
  ggplot() + 
    geom_histogram(aes(x=resid)) +
    theme_minimal()
```

```{r}
plot(lm_1)
```


### Consigna 8.
¿Parece que la condición de variabilidad constante del modelo es respetada? ¿Por qué?
```{r}
###
```

## Análisis de resultados
### Consigna 9. 
¿Cuál es el ajuste del modelo?
```{r}
summary(lm_1)$adj.r.squared
```

### Consigna 10. 
¿Cómo interpreta los resultados del modelo? ¿Cuál es, en su opinion, la variable más importante? 
```{r}
lm_scaled <- final %>%
        select(gnipc_2019, le_2019, mys_2019, SL.AGR.EMPL.ZS, SL.SRV.EMPL.ZS, SL.EMP.WORK.ZS) %>%
        mutate(across(everything(), scale)) %>%
        lm(gnipc_2019 ~ ., data=.)
```

```{r}
tidy(lm_scaled)
```

