Gapminder del doctor Hans Rosling

Esta clase nos ayudará a 3 metodoligías diferntes para enfocar diferentes modelos:

Vamos a trabajar con un dataset que nos muestra sobre la esperanza de vida en diferentes paises del mundo. Vamos a intentar modelar a partir de cada visualización, usando sumaries de modelos para identificar outliers,… .

library(tidyverse)
## -- Attaching packages ----------------------------------------------------------------------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 3.1.0     v purrr   0.2.5
## v tibble  1.4.2     v dplyr   0.7.8
## v tidyr   0.8.2     v stringr 1.3.1
## v readr   1.1.1     v forcats 0.3.0
## -- Conflicts -------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(modelr)
library(gapminder)

Vamos a empezar a trabajar con el dataset:

gapminder %>%
  ggplot(aes(year, lifeExp, group = country)) +
    geom_line(alpha =0.2)

Life expectancy respecto los años, parece que ha habido un incremento, pero hay algunos que no. ¿Cómo podemos saberlo delante de un grafico tan complejo como éste? Hacemos filtro de “España”:

Selección de país y modelo

es <- filter(gapminder, country == "Spain") 

es %>%
  ggplot(aes(year, lifeExp)) + geom_line() + ggtitle("Full data") +
    theme(plot.title = element_text(hjust = 0.5))

Vamos a predecir la esperanza de vida en función del año.

es_mod <- lm(lifeExp ~year, data = es)
# Anadimos las predicciones
es %>%
  add_predictions(es_mod) %>%
  ggplot(aes(year, pred)) +
    geom_line() +
    ggtitle("Linear trend") +
    theme(plot.title = element_text(hjust = 0.5))

Hemos analizado la tendencia, pero debemos mirar como han performado los residuos.

es %>%
  add_residuals(es_mod) %>%
  ggplot(aes(year, resid)) +
    geom_hline(yintercept = 0, color = "white", size = 3) +
    geom_line() +
    ggtitle("Residual pattern") +
    theme(plot.title = element_text(hjust = 0.5))

El problema es hacer todos estos modelos para cada uno de los 140 paises. Debemos ver como agrupar y anidar información en un dataset a partir de como está montado desde el inicio.

El data frame anidado

Para crear un dataframe anidado hay que empezar por agrupar el dataframe y después irlo anidando según las variables que nos hacen falta:

by_country <- gapminder %>%
  group_by(country, continent) %>%
  nest() # esta función nos permite crear tibbles dentro del dataframe
by_country

Fijémonos que cada una de las filas contiene un dataframe. Ejemplo de visualización de un dataframe dentro del dataframe:

by_country$data[[1]]

Modelando y operando con data frames anidados

Vamos a crear una función para que cree modelos para cada uno de los datasets del dataset “by_country”.

country_model <- function(df){
  lm(lifeExp ~ year, data = df)
}

Ahora debemos pasarle todos los dataframes, pero como estan en una lista, usaremos la función map():

models <- map(by_country$data, country_model)

Ahora almacenaremos el modelo para cada país en una nueva columna. El hecho que tengamos los valores en listas nos hace el trabajo muy fácil:

by_country <- by_country %>%
  mutate(model = map(data, country_model))
# Visualicemos como ha cambiado el dataset:
head(by_country)

Ahora ha aparecido una nueva columna en el dataset que nos indica que hay un nuevo modelo introducido en él.

Operaciones con el dataframe:

by_country %>%
  filter(continent == "Europe") %>%
  head()
by_country %>%
  arrange(continent, country) %>%
  head()

Añadiendo los residuos

by_country %>%
  mutate(resids = map2(data, model, add_residuals)) -> by_country
head(by_country)

Visualizamos la información

Ahora desagruparemos la información con la función unnest()

resids <- unnest(by_country, resids)
head(resids)

Vamos a ver los resultados de los modelos para cada residuo mirando a cada país:

resids %>%
  ggplot(aes(year, resid)) + 
    geom_line(aes(group = country), alpha = 0.2) +
    geom_smooth(se = FALSE)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Vamos a separalo por continente para verlo mejor:

resids %>%
  ggplot(aes(year, resid)) + 
    geom_line(aes(group = country), alpha = 0.2) +
    geom_smooth(se = FALSE)+
    facet_wrap(~continent)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Podemos ver que hay ciertos patrones que aparecen y que el modelo no explica. Vamos a mejorar la calidad de nuestro modelo:

Midiendo la calidad del modelo con broom

library(broom)
## 
## Attaching package: 'broom'
## The following object is masked from 'package:modelr':
## 
##     bootstrap

Este paquete nos proprciona un conjunto de herramientas para limpiar cualquier dataset. Primero usaremos la función gllance() para medir la calidad del modelo:

glance(es_mod)

Nos muestran una serie de medidas muy representativas. Vamos a introdirle los resultados en el dataset:

glance <- by_country %>%
  mutate(glance = map(model, glance)) %>%
  unnest(glance, .drop = TRUE) # con drop eliminamos el resto de columnas que no nos interessan

head(glance)

Usaremos el factor de r squares para saber los paises que nos performan peor:

glance %>%
  arrange(r.squared) %>%
  head()

Los peores modelos que hemos generado parece que viven en África. Vamos a visualizar los peores resultados.

glance %>%
  ggplot(aes(continent, r.squared)) +
    geom_jitter(width = 0.5)

Podemos ver como así es, la mayoría de modelos malos se encuentran en África. Vamos a intentar mejorar su interpretación:

bad_fit <- filter(glance, r.squared < 0.25) 
bad_fit
gapminder %>%
  semi_join(bad_fit, by = "country") %>% # semi_join nos permite hacer un join de todos los paises originales junto con bad_fit, solo nos quedan los datos con las variables de bad_fit, como si hiceramos una select
  ggplot(aes(year, lifeExp, color = country)) +
    geom_line()

Estas caidas tan pronunciadas son debido a factores sociopolíticos.