Capítulo 6 Loops (purrr) y bibliografía (rticles)

6.1 Paquetes necesarios para este capítulo

Para este capítulo necesitas tener instalado el paquete tidyverse.

Probablemente uno de los puntos que marca la diferencia entre ser un usuario de un lenguaje de programación y un alguién que realmente programa. Es el momento en que una persona aprende a hacer loops. Los loops son una acción repetitiva en la cual una misma acción es realizada por el computador ahorrandonos mucho tiempo de escribir código y muchas veces tiempo de computación tambien.

Existen varias formas de como realizar loops en R, los for loops, la familia de los apply y más recientemente el uso del paquete purrr (Henry and Wickham 2018) presente en el tidyverse. En este capítulo nos enfocaremos principalmente en el uso de este paquete, pero también explicaremos levemente el caso de los for loops.

Dado que este libro es un apoyo para el curso BIO4022, esta clase puede también ser seguida en este link. El video de la clase se encontrará disponible en este link.

6.2 Generando una receta

Como hacer un loop, es una repetición de un código multiples veces, generalmente lo que más nos combiene es generar la receta tomando en cuenta el primer elemento y luego repetirlo en un loop.

6.2.1 Dioxido de nitrógeno en Madrid

Supongamos que queremos estudiar la concentración de dióxido de Nitrógeno en madrid en distintas estaciones, la base de datos puede ser encontrada en el siguiente link. Dentro de esta base de datos tenemos una carpeta con la calidad de aire de estaciones en Madrid, con un archivo para cada año. Supongamos que se quiere hacer lo siguiente, limitandose a las estaciones de Cuatro Caminos, El Pardo, Escuelas Aguirre, Moratalaz y Tres Olivos, calcular los promedios de \(NO_2\) para cada mes y cada año en estas estaciones.

6.2.1.1 Generando la receta

Esto lo podemos hacer con un loop, pero antes generemos la receta tomando en cuenta solo el 2017.

Para esto hacemos lo siguiente:

  • Tomemos la base de datos de calidad de aire de Madrid
  • Leeamos el año 2017
  • Limitemonos a las estaciones de Cuatro Caminos, El Pardo, Escuelas Aguirre, Moratalaz y Tres Olivos
  • Agreguemos una columna con el año y una con el mes
  • Calculemos los promedios de \(NO_2\) para cada mes
  • Eliminemos las columnas innecesarias para estudiar el efecto del \(NO_2\) en Madrid

Vamos paso a paso

6.2.1.1.1 leyendo la base de datos

El primer paso es leer la base de datos, para esto usamos el tidyverse y cargamos además lubridate por si tenemos que trabajar con las fechas. En la tabla 6.1 vemos los resultados del código a continuación.

library(tidyverse)
library(lubridate)
Madrid2017 <- read_csv("csvs_per_year/madrid_2017.csv")
Tabla 6.1: Los primeros 20 datos de calidad de aire del 2017 en Madrid para todas las estaciones.
date BEN CH4 CO EBE NMHC NO NO_2 NOx O_3 PM10 PM25 SO_2 TCH TOL station
2017-06-01 01:00:00 NA NA 0.3 NA NA 4 38 NA NA NA NA 5 NA NA 28079004
2017-06-01 01:00:00 0.6 NA 0.3 0.4 0.08 3 39 NA 71 22 9 7 1.40 2.9 28079008
2017-06-01 01:00:00 0.2 NA NA 0.1 NA 1 14 NA NA NA NA NA NA 0.9 28079011
2017-06-01 01:00:00 NA NA 0.2 NA NA 1 9 NA 91 NA NA NA NA NA 28079016
2017-06-01 01:00:00 NA NA NA NA NA 1 19 NA 69 NA NA 2 NA NA 28079017
2017-06-01 01:00:00 0.1 NA 0.3 0.2 NA 1 26 NA 70 26 NA 1 NA 0.3 28079018
2017-06-01 01:00:00 0.3 NA 0.2 0.1 0.17 1 19 NA 79 23 9 3 0.86 1.8 28079024
2017-06-01 01:00:00 NA NA NA NA NA 1 9 NA 87 NA NA NA NA NA 28079027
2017-06-01 01:00:00 NA NA 0.3 NA NA 3 30 NA 70 NA NA NA NA NA 28079035
2017-06-01 01:00:00 NA NA 0.1 NA NA 1 15 NA NA 22 NA 10 NA NA 28079036
2017-06-01 01:00:00 0.7 NA NA 1.0 NA 1 25 NA NA 21 10 2 NA 3.5 28079038
2017-06-01 01:00:00 NA NA 0.2 NA NA 1 21 NA 75 NA NA NA NA NA 28079039
2017-06-01 01:00:00 NA NA NA NA NA 2 17 NA NA 24 NA 9 NA NA 28079040
2017-06-01 01:00:00 NA NA NA NA NA 1 22 NA NA 23 15 NA NA NA 28079047
2017-06-01 01:00:00 NA NA NA NA NA 2 30 NA NA 17 9 NA NA NA 28079048
2017-06-01 01:00:00 NA NA NA NA NA 1 12 NA 74 NA NA NA NA NA 28079049
2017-06-01 01:00:00 NA NA NA NA NA 2 15 NA NA 16 12 NA NA NA 28079050
2017-06-01 01:00:00 NA NA NA NA NA 3 12 NA 84 NA NA NA NA NA 28079054
2017-06-01 01:00:00 0.2 NA NA 0.6 0.08 1 12 NA NA 15 NA NA 1.16 1.5 28079055
2017-06-01 01:00:00 NA NA 0.1 NA NA 9 47 NA 59 NA NA NA NA NA 28079056
6.2.1.1.2 Limitemonos a las estaciones seleccionadas

Revisando el archivo stations.csv, podemos ver que el código de estaciones que estudiaremos son 28079036, 28079008,28079058, 28079060 y 28079038, por lo que lo ponemos en un filter. El resultado de esto lo podemos ver en la tabla 6.2

Madrid2017 <- read_csv("csvs_per_year/madrid_2017.csv") %>% filter(station %in% 
    c(28079036, 28079008, 28079058, 28079060, 28079038))
Tabla 6.2: Los primeros 20 datos de calidad de aire del 2017 en Madrid después de filtrar según estación.
date BEN CH4 CO EBE NMHC NO NO_2 NOx O_3 PM10 PM25 SO_2 TCH TOL station
2017-06-01 01:00:00 0.6 NA 0.3 0.4 0.08 3 39 NA 71 22 9 7 1.40 2.9 28079008
2017-06-01 01:00:00 NA NA 0.1 NA NA 1 15 NA NA 22 NA 10 NA NA 28079036
2017-06-01 01:00:00 0.7 NA NA 1.0 NA 1 25 NA NA 21 10 2 NA 3.5 28079038
2017-06-01 01:00:00 NA NA NA NA NA 1 10 NA 66 NA NA NA NA NA 28079058
2017-06-01 01:00:00 NA NA NA NA NA 1 26 NA 79 86 NA NA NA NA 28079060
2017-06-01 02:00:00 0.3 NA 0.3 0.2 0.07 2 27 NA 72 16 7 7 1.40 2.3 28079008
2017-06-01 02:00:00 NA NA 0.1 NA NA 1 13 NA NA 17 NA 10 NA NA 28079036
2017-06-01 02:00:00 0.2 NA NA 0.5 NA 9 20 NA NA 13 4 2 NA 1.3 28079038
2017-06-01 02:00:00 NA NA NA NA NA 1 11 NA 68 NA NA NA NA NA 28079058
2017-06-01 02:00:00 NA NA NA NA NA 1 9 NA 90 62 NA NA NA NA 28079060
2017-06-01 03:00:00 0.3 NA 0.3 0.2 0.08 2 25 NA 73 14 7 7 1.40 2.0 28079008
2017-06-01 03:00:00 NA NA 0.1 NA NA 1 11 NA NA 18 NA 10 NA NA 28079036
2017-06-01 03:00:00 0.1 NA NA 0.4 NA 6 20 NA NA 11 6 2 NA 1.8 28079038
2017-06-01 03:00:00 NA NA NA NA NA 1 6 NA 68 NA NA NA NA NA 28079058
2017-06-01 03:00:00 NA NA NA NA NA 1 8 NA 86 19 NA NA NA NA 28079060
2017-06-01 04:00:00 0.3 NA 0.2 0.2 0.08 2 22 NA 75 15 10 6 1.41 1.4 28079008
2017-06-01 04:00:00 NA NA 0.1 NA NA 1 14 NA NA 12 NA 10 NA NA 28079036
2017-06-01 04:00:00 0.2 NA NA 0.5 NA 1 12 NA NA 10 6 2 NA 1.7 28079038
2017-06-01 04:00:00 NA NA NA NA NA 1 6 NA 60 NA NA NA NA NA 28079058
2017-06-01 04:00:00 NA NA NA NA NA 1 11 NA 75 8 NA NA NA NA 28079060
6.2.1.1.3 Agreguemos aparte el mes, el año y el nombre de la estación

Usando mutate y las funciones monthy year de lubridate podemos agregar el més y el año para cada observación, además usando left_joint, podemos agreagar el nombre de las estaciones usando la base de datos stations.csv. El resultado de esto lo podemos ver en la tabla 6.3

stations <- read_csv("stations.csv") %>% rename(station = id)
Madrid2017 <- read_csv("csvs_per_year/madrid_2017.csv") %>% filter(station %in% 
    c(28079036, 28079008, 28079058, 28079060, 28079038)) %>% 
    mutate(month = month(date), year = year(date)) %>% left_join(stations)
Tabla 6.3: Los primeros 20 datos de calidad de aire del 2017 en Madrid después de filtrar según estación, con mes, año y nombre.
date BEN CH4 CO EBE NMHC NO NO_2 NOx O_3 PM10 PM25 SO_2 TCH TOL station month year name address lon lat elevation
2017-06-01 01:00:00 0.6 NA 0.3 0.4 0.08 3 39 NA 71 22 9 7 1.40 2.9 28079008 6 2017 Escuelas Aguirre Entre C/ Alcalá y C/ O’ Donell -3.682319 40.42156 670
2017-06-01 01:00:00 NA NA 0.1 NA NA 1 15 NA NA 22 NA 10 NA NA 28079036 6 2017 Moratalaz Avd. Moratalaz esq. Camino de los Vinateros -3.645306 40.40795 685
2017-06-01 01:00:00 0.7 NA NA 1.0 NA 1 25 NA NA 21 10 2 NA 3.5 28079038 6 2017 Cuatro Caminos Avda. Pablo Iglesias esq. C/ Marqués de Lema -3.707128 40.44554 698
2017-06-01 01:00:00 NA NA NA NA NA 1 10 NA 66 NA NA NA NA NA 28079058 6 2017 El Pardo Avda. La Guardia -3.774611 40.51806 615
2017-06-01 01:00:00 NA NA NA NA NA 1 26 NA 79 86 NA NA NA NA 28079060 6 2017 Tres Olivos Plaza Tres Olivos -3.689761 40.50059 715
2017-06-01 02:00:00 0.3 NA 0.3 0.2 0.07 2 27 NA 72 16 7 7 1.40 2.3 28079008 6 2017 Escuelas Aguirre Entre C/ Alcalá y C/ O’ Donell -3.682319 40.42156 670
2017-06-01 02:00:00 NA NA 0.1 NA NA 1 13 NA NA 17 NA 10 NA NA 28079036 6 2017 Moratalaz Avd. Moratalaz esq. Camino de los Vinateros -3.645306 40.40795 685
2017-06-01 02:00:00 0.2 NA NA 0.5 NA 9 20 NA NA 13 4 2 NA 1.3 28079038 6 2017 Cuatro Caminos Avda. Pablo Iglesias esq. C/ Marqués de Lema -3.707128 40.44554 698
2017-06-01 02:00:00 NA NA NA NA NA 1 11 NA 68 NA NA NA NA NA 28079058 6 2017 El Pardo Avda. La Guardia -3.774611 40.51806 615
2017-06-01 02:00:00 NA NA NA NA NA 1 9 NA 90 62 NA NA NA NA 28079060 6 2017 Tres Olivos Plaza Tres Olivos -3.689761 40.50059 715
2017-06-01 03:00:00 0.3 NA 0.3 0.2 0.08 2 25 NA 73 14 7 7 1.40 2.0 28079008 6 2017 Escuelas Aguirre Entre C/ Alcalá y C/ O’ Donell -3.682319 40.42156 670
2017-06-01 03:00:00 NA NA 0.1 NA NA 1 11 NA NA 18 NA 10 NA NA 28079036 6 2017 Moratalaz Avd. Moratalaz esq. Camino de los Vinateros -3.645306 40.40795 685
2017-06-01 03:00:00 0.1 NA NA 0.4 NA 6 20 NA NA 11 6 2 NA 1.8 28079038 6 2017 Cuatro Caminos Avda. Pablo Iglesias esq. C/ Marqués de Lema -3.707128 40.44554 698
2017-06-01 03:00:00 NA NA NA NA NA 1 6 NA 68 NA NA NA NA NA 28079058 6 2017 El Pardo Avda. La Guardia -3.774611 40.51806 615
2017-06-01 03:00:00 NA NA NA NA NA 1 8 NA 86 19 NA NA NA NA 28079060 6 2017 Tres Olivos Plaza Tres Olivos -3.689761 40.50059 715
2017-06-01 04:00:00 0.3 NA 0.2 0.2 0.08 2 22 NA 75 15 10 6 1.41 1.4 28079008 6 2017 Escuelas Aguirre Entre C/ Alcalá y C/ O’ Donell -3.682319 40.42156 670
2017-06-01 04:00:00 NA NA 0.1 NA NA 1 14 NA NA 12 NA 10 NA NA 28079036 6 2017 Moratalaz Avd. Moratalaz esq. Camino de los Vinateros -3.645306 40.40795 685
2017-06-01 04:00:00 0.2 NA NA 0.5 NA 1 12 NA NA 10 6 2 NA 1.7 28079038 6 2017 Cuatro Caminos Avda. Pablo Iglesias esq. C/ Marqués de Lema -3.707128 40.44554 698
2017-06-01 04:00:00 NA NA NA NA NA 1 6 NA 60 NA NA NA NA NA 28079058 6 2017 El Pardo Avda. La Guardia -3.774611 40.51806 615
2017-06-01 04:00:00 NA NA NA NA NA 1 11 NA 75 8 NA NA NA NA 28079060 6 2017 Tres Olivos Plaza Tres Olivos -3.689761 40.50059 715

Finalmente, agrupamos sacamos el promedio por mes y sacamos las columnas sobrantes al mismo tiempo, como vemos en la tabla 6.4

library(lubridate)
stations <- read_csv("stations.csv") %>% rename(station = id)
## Parsed with column specification:
## cols(
##   id = col_integer(),
##   name = col_character(),
##   address = col_character(),
##   lon = col_double(),
##   lat = col_double(),
##   elevation = col_integer()
## )
Madrid2017 <- read_csv("csvs_per_year/madrid_2017.csv") %>% filter(station %in% 
    c(28079036, 28079008, 28079058, 28079060, 28079038)) %>% 
    mutate(month = month(date), year = year(date)) %>% left_join(stations) %>% 
    group_by(month, name, year) %>% summarise(NO_2 = mean(NO_2, 
    na.rm = TRUE))
## Parsed with column specification:
## cols(
##   date = col_datetime(format = ""),
##   BEN = col_double(),
##   CH4 = col_character(),
##   CO = col_double(),
##   EBE = col_double(),
##   NMHC = col_double(),
##   NO = col_double(),
##   NO_2 = col_double(),
##   NOx = col_character(),
##   O_3 = col_double(),
##   PM10 = col_double(),
##   PM25 = col_double(),
##   SO_2 = col_double(),
##   TCH = col_double(),
##   TOL = col_double(),
##   station = col_integer()
## )
## Joining, by = "station"
## Warning: Detecting old grouped_df format, replacing `vars` attribute by
## `groups`
Tabla 6.4: Los primeros 20 datos de calidad de aire del 2017 en Madrid después de filtrar según estación, con mes, año y nombre.
month name year NO_2
1 Cuatro Caminos 2017 59.88124
1 Cuatro Caminos 2018 5.00000
1 El Pardo 2017 27.89892
1 El Pardo 2018 1.00000
1 Escuelas Aguirre 2017 69.43666
1 Escuelas Aguirre 2018 20.00000
1 Moratalaz 2017 54.17004
1 Moratalaz 2018 10.00000
1 Tres Olivos 2017 55.28706
1 Tres Olivos 2018 5.00000
2 Cuatro Caminos 2017 46.01045
2 El Pardo 2017 17.18955
2 Escuelas Aguirre 2017 59.89686
2 Moratalaz 2017 44.19581
2 Tres Olivos 2017 39.38209
3 Cuatro Caminos 2017 43.71833
3 El Pardo 2017 15.07962
3 Escuelas Aguirre 2017 61.74798
3 Moratalaz 2017 43.01353
3 Tres Olivos 2017 36.68329
6.2.1.1.4 Últimos detalles

Vemos que hay algunos valores del 2018, esto parece raro, ya que leimos los archivos del 2017. Al revisar mas con summarize, vemos que en realidad son tan solo unas pocas observaciones las que generan esta anomalía debido a algunas medidas del 1 de enero del 2018.

Para eliminarlas agregamos el siguiente código.

library(lubridate)
stations <- read_csv("stations.csv") %>% rename(station = id)
## Parsed with column specification:
## cols(
##   id = col_integer(),
##   name = col_character(),
##   address = col_character(),
##   lon = col_double(),
##   lat = col_double(),
##   elevation = col_integer()
## )
Madrid2017 <- read_csv("csvs_per_year/madrid_2017.csv") %>% filter(station %in% 
    c(28079036, 28079008, 28079058, 28079060, 28079038)) %>% 
    mutate(month = month(date), year = year(date)) %>% left_join(stations) %>% 
    group_by(month, name, year) %>% summarise(NO_2 = mean(NO_2, 
    na.rm = TRUE), n = n()) %>% filter(n > 500)
## Parsed with column specification:
## cols(
##   date = col_datetime(format = ""),
##   BEN = col_double(),
##   CH4 = col_character(),
##   CO = col_double(),
##   EBE = col_double(),
##   NMHC = col_double(),
##   NO = col_double(),
##   NO_2 = col_double(),
##   NOx = col_character(),
##   O_3 = col_double(),
##   PM10 = col_double(),
##   PM25 = col_double(),
##   SO_2 = col_double(),
##   TCH = col_double(),
##   TOL = col_double(),
##   station = col_integer()
## )
## Joining, by = "station"

Esto nos dá al fin la receta final que usaremos en el loop.

6.3 Empezando el loop

En este capíulo usaremos principalmente la función map del paquete purrr para generar loops, en esta función los dos argumentos generales que necesitamos es un vector o lista (argumento .x) de los elementos que pasarán por una función, y una funcion (argumento .f) que se aplicará a toda esta lista. Es importante establecer que el resultado de map siempre será una lista.

6.3.1 Volvamos a nuestra receta

Veamos el código que usamos para el año 2017

library(lubridate)
stations <- read_csv("stations.csv") %>% rename(station = id)
Madrid2017 <- read_csv("csvs_per_year/madrid_2017.csv") %>% filter(station %in% 
    c(28079036, 28079008, 28079058, 28079060, 28079038)) %>% 
    mutate(month = month(date), year = year(date)) %>% left_join(stations) %>% 
    group_by(month, name, year) %>% summarise(NO_2 = mean(NO_2, 
    na.rm = TRUE), n = n()) %>% filter(n > 500)

La primera parte del código es la lectura del archivo

Madrid2017 <- read_csv("csvs_per_year/madrid_2017.csv")

Para hacer esto por todos los archivos de la base de datos requeriríamos de una lista o vector con los nombres de cada uno de los archivos. ¡Si solo hubiera una función en R que nos permitiera leer los archivos de una carpeta! La función list.files hace eso.

Entonces el código que vemos abajo genera un vector con todos los nombres de los archivos que queremos incorporar:

Archivos <- list.files("csvs_per_year", full.names = TRUE)
Archivos
##  [1] "csvs_per_year/madrid_2001.csv" "csvs_per_year/madrid_2002.csv"
##  [3] "csvs_per_year/madrid_2003.csv" "csvs_per_year/madrid_2004.csv"
##  [5] "csvs_per_year/madrid_2005.csv" "csvs_per_year/madrid_2006.csv"
##  [7] "csvs_per_year/madrid_2007.csv" "csvs_per_year/madrid_2008.csv"
##  [9] "csvs_per_year/madrid_2009.csv" "csvs_per_year/madrid_2010.csv"
## [11] "csvs_per_year/madrid_2011.csv" "csvs_per_year/madrid_2012.csv"
## [13] "csvs_per_year/madrid_2013.csv" "csvs_per_year/madrid_2014.csv"
## [15] "csvs_per_year/madrid_2015.csv" "csvs_per_year/madrid_2016.csv"
## [17] "csvs_per_year/madrid_2017.csv" "csvs_per_year/madrid_2018.csv"

Entonces poner dentro de map, un vector con el nombre de los archivos (Archivos), y una función para leer los archivos (read_csv). Esto es el siguiente código

Madrid <- map(Archivos, read_csv)

Genera una lista donde cada elemento es un data frame de un año de mediciones.

Cuando se agregan otras funciones mas complejas en un loop usando map. Como por ejemplo filter, ponemos un simbolo ~ dentro de map, y un .x dentro de filter para representar a cada dataframe que usaremos.

Madrid <- map(Archivos, read_csv) %>% map(~filter(.x, station %in% 
    c(28079036, 28079008, 28079058, 28079060, 28079038)))

De esta forma podemos seguir la receta creada anteriormente sin ningún problema.

Madrid <- map(Archivos, read_csv) %>% map(~filter(.x, station %in% 
    c(28079036, 28079008, 28079058, 28079060, 28079038))) %>% 
    map(~mutate(.x, month = month(date), year = year(date))) %>% 
    map(~left_join(.x, stations)) %>% map(~group_by(.x, month, 
    name, year)) %>% map(~summarise(.x, NO_2 = mean(NO_2, na.rm = TRUE), 
    n = n())) %>% map(~filter(.x, n > 500))

Pero en este momento tenemos una lista con 17 data frames, en vez de un gran data frame con todos los datos. Para esto debenos unir esta lista usando la función reduce, lo cual nos genera el siguiente código y la tabla 6.5

library(lubridate)
Madrid <- map(Archivos, read_csv) %>% map(~filter(.x, station %in% 
    c(28079036, 28079008, 28079058, 28079060, 28079038))) %>% 
    map(~mutate(.x, month = month(date), year = year(date))) %>% 
    map(~left_join(.x, stations)) %>% map(~group_by(.x, month, 
    name, year)) %>% map(~summarise(.x, NO_2 = mean(NO_2, na.rm = TRUE), 
    n = n())) %>% map(~filter(.x, n > 500)) %>% reduce(bind_rows)
## Warning: Detecting old grouped_df format, replacing `vars` attribute by
## `groups`
Tabla 6.5: Los primeros 20 datos de calidad de aire del todos los años en Madrid después de filtrar según estación, con mes, año y nombre.
month name year NO_2 n
1 Cuatro Caminos 2017 59.88124 743
1 El Pardo 2017 27.89892 743
1 Escuelas Aguirre 2017 69.43666 743
1 Moratalaz 2017 54.17004 743
1 Tres Olivos 2017 55.28706 743
2 Cuatro Caminos 2017 46.01045 672
2 El Pardo 2017 17.18955 672
2 Escuelas Aguirre 2017 59.89686 672
2 Moratalaz 2017 44.19581 672
2 Tres Olivos 2017 39.38209 672
3 Cuatro Caminos 2017 43.71833 744
3 El Pardo 2017 15.07962 744
3 Escuelas Aguirre 2017 61.74798 744
3 Moratalaz 2017 43.01353 744
3 Tres Olivos 2017 36.68329 744
4 Cuatro Caminos 2017 34.35139 720
4 El Pardo 2017 11.48122 720
4 Escuelas Aguirre 2017 51.56485 720
4 Moratalaz 2017 30.07232 720
4 Tres Olivos 2017 25.51253 720

Referencias

Henry, Lionel, and Hadley Wickham. 2018. Purrr: Functional Programming Tools. https://CRAN.R-project.org/package=purrr.