Divide-aplica-combina (split-apply-combine)

Muchos problemas de análisis de datos involucran la aplicación de la estrategia divide-aplica-combina, (Hadley Whickam, 2011) esta consiste en romper un problema en pedazos (de acuerdo a una variable de interés), operar sobre cada subcomjunto de manera independiente (ej. calcular la media de cada grupo, ordenar observaciones por grupo, estandarizar por grupo) y después unir los pedazos nuevamente. El siguiente diagrama ejemplifiaca el paradigma de divide-aplica-combina:

En esta sección trabajaremos con las siguientes bases de datos para ejemplifcar las funciones de divide-aplica-combina:

flights <- tbl_df(read.csv("data/flights.csv", stringsAsFactors = FALSE))
flights
## Source: local data frame [227,496 x 14]
## 
##                   date  hour minute   dep   arr dep_delay arr_delay
##                  (chr) (int)  (int) (int) (int)     (int)     (int)
## 1  2011-01-01 12:00:00    14      0  1400  1500         0       -10
## 2  2011-01-02 12:00:00    14      1  1401  1501         1        -9
## 3  2011-01-03 12:00:00    13     52  1352  1502        -8        -8
## 4  2011-01-04 12:00:00    14      3  1403  1513         3         3
## 5  2011-01-05 12:00:00    14      5  1405  1507         5        -3
## 6  2011-01-06 12:00:00    13     59  1359  1503        -1        -7
## 7  2011-01-07 12:00:00    13     59  1359  1509        -1        -1
## 8  2011-01-08 12:00:00    13     55  1355  1454        -5       -16
## 9  2011-01-09 12:00:00    14     43  1443  1554        43        44
## 10 2011-01-10 12:00:00    14     43  1443  1553        43        43
## ..                 ...   ...    ...   ...   ...       ...       ...
## Variables not shown: carrier (chr), flight (int), dest (chr), plane (chr),
##   cancelled (int), time (int), dist (int)
class(flights$date)
## [1] "character"
flights$date <- as.Date(flights$date)

weather <- tbl_df(read.csv("data/weather.csv", stringsAsFactors = FALSE))
weather$date <- as.Date(weather$date)
weather 
## Source: local data frame [8,723 x 14]
## 
##          date  hour  temp dew_point humidity pressure visibility wind_dir
##        (date) (int) (dbl)     (dbl)    (int)    (dbl)      (dbl)    (chr)
## 1  2011-01-01     0  59.0      28.9       32    29.86         10      NNE
## 2  2011-01-01     1  57.2      28.4       33    29.88         10      NNE
## 3  2011-01-01     2  55.4      28.4       36    29.93         10      NNW
## 4  2011-01-01     3  53.6      28.4       38    29.94         10    North
## 5  2011-01-01     4    NA        NA       NA    29.99         10      NNW
## 6  2011-01-01     5    NA        NA       NA    30.02         10    North
## 7  2011-01-01     6  53.1      17.1       24    30.05         10    North
## 8  2011-01-01     7  53.1      16.0       23    30.07         10    North
## 9  2011-01-01     8  54.0      18.0       24    30.09         10    North
## 10 2011-01-01     9  55.4      17.6       23    30.09         10      NNE
## ..        ...   ...   ...       ...      ...      ...        ...      ...
## Variables not shown: wind_dir2 (int), wind_speed (dbl), gust_speed (dbl),
##   precip (dbl), conditions (chr), events (chr)
planes <- tbl_df(read.csv("data/planes.csv", stringsAsFactors = FALSE))
planes
## Source: local data frame [2,853 x 9]
## 
##     plane  year               mfr          model no.eng no.seats speed
##     (chr) (int)             (chr)          (chr)  (int)    (int) (int)
## 1  N576AA  1991 MCDONNELL DOUGLAS DC-9-82(MD-82)      2      172    NA
## 2  N557AA  1993        MARZ BARRY      KITFOX IV      1        2    NA
## 3  N403AA  1974             RAVEN           S55A     NA        1    60
## 4  N492AA  1989 MCDONNELL DOUGLAS DC-9-82(MD-82)      2      172    NA
## 5  N262AA  1985 MCDONNELL DOUGLAS DC-9-82(MD-82)      2      172    NA
## 6  N493AA  1989 MCDONNELL DOUGLAS DC-9-82(MD-82)      2      172    NA
## 7  N477AA  1988 MCDONNELL DOUGLAS DC-9-82(MD-82)      2      172    NA
## 8  N476AA  1988 MCDONNELL DOUGLAS DC-9-82(MD-82)      2      172    NA
## 9  N504AA    NA AUTHIER ANTHONY P      TIERRA II      1        2    NA
## 10 N565AA  1987 MCDONNELL DOUGLAS DC-9-83(MD-83)      2      172    NA
## ..    ...   ...               ...            ...    ...      ...   ...
## Variables not shown: engine (chr), type (chr)
airports <- tbl_df(read.csv("data/airports.csv", stringsAsFactors = FALSE))
airports
## Source: local data frame [3,376 x 7]
## 
##     iata              airport             city state country      lat
##    (chr)                (chr)            (chr) (chr)   (chr)    (dbl)
## 1    00M             Thigpen       Bay Springs    MS     USA 31.95376
## 2    00R Livingston Municipal       Livingston    TX     USA 30.68586
## 3    00V          Meadow Lake Colorado Springs    CO     USA 38.94575
## 4    01G         Perry-Warsaw            Perry    NY     USA 42.74135
## 5    01J     Hilliard Airpark         Hilliard    FL     USA 30.68801
## 6    01M    Tishomingo County          Belmont    MS     USA 34.49167
## 7    02A          Gragg-Wade           Clanton    AL     USA 32.85049
## 8    02C              Capitol       Brookfield    WI     USA 43.08751
## 9    02G    Columbiana County   East Liverpool    OH     USA 40.67331
## 10   03D     Memphis Memorial          Memphis    MO     USA 40.44726
## ..   ...                  ...              ...   ...     ...      ...
## Variables not shown: long (dbl)

Cuando pensamos como implementar la estrategia divide-aplica-combina es natural pensar en iteraciones, por ejemplo utilizar un ciclo for para recorrer cada grupo de interés y aplicar las funciones, sin embargo la aplicación de ciclos for desemboca en código difícil de entender. Adicionalmente, dplyr es mucho más veloz.

Estudiaremos las siguientes funciones:

Estas funciones trabajan de manera similar, el primer argumento que reciben es un data frame (usualmente en formato limpio), los argumentos que siguen indican que operación se va a efectuar y el resultado es un nuevo data frame.

Veamos con ejemplos.

Filtrar

Creamos una base de datos de juguete para mostrar el funcionamiento de cada instrucción:

df_ej <- data.frame(genero = c("mujer", "hombre", "mujer", "mujer", "hombre"), 
                    estatura = c(1.65, 1.80, 1.70, 1.60, 1.67))
df_ej
##   genero estatura
## 1  mujer     1.65
## 2 hombre     1.80
## 3  mujer     1.70
## 4  mujer     1.60
## 5 hombre     1.67
filter(df_ej, genero == "mujer")
##   genero estatura
## 1  mujer     1.65
## 2  mujer     1.70
## 3  mujer     1.60
filter(df_ej, estatura > 1.65 & estatura < 1.75)
##   genero estatura
## 1  mujer     1.70
## 2 hombre     1.67

Algunos operadores importantes para filtrar son:

x > 1
x >= 1
x < 1
x <= 1
x != 1
x == 1
x %in% ("a", "b")

# Conjuntos
a | b
a & b
a & !b
xor(a, b)

Ejercicio

Seleccionar

Elegir columnas de un conjunto de datos.

df_ej
##   genero estatura
## 1  mujer     1.65
## 2 hombre     1.80
## 3  mujer     1.70
## 4  mujer     1.60
## 5 hombre     1.67
select(df_ej, genero)
##   genero
## 1  mujer
## 2 hombre
## 3  mujer
## 4  mujer
## 5 hombre
select(df_ej, -genero)
##   estatura
## 1     1.65
## 2     1.80
## 3     1.70
## 4     1.60
## 5     1.67
select(df_ej, starts_with("g"))
select(df_ej, contains("g"))

Ejercicio

Ve la ayuda de select (?select) y escribe tres maneras de seleccionar las variables de retraso (delay).

Arreglar

Arreglar u ordenar de acuerdo al valor de una o más variables:

arrange(df_ej, genero)
##   genero estatura
## 1 hombre     1.80
## 2 hombre     1.67
## 3  mujer     1.65
## 4  mujer     1.70
## 5  mujer     1.60
arrange(df_ej, desc(estatura))
##   genero estatura
## 1 hombre     1.80
## 2  mujer     1.70
## 3 hombre     1.67
## 4  mujer     1.65
## 5  mujer     1.60

Ejercicio

Ordena los vuelos por fecha de salida y hora.

¿Cuáles son los vuelos con mayor retraso?

¿Qué vuelos ganaron más tiempo en el aire?

Mutar

Mutar consiste en crear nuevas variables:

mutate(df_ej, estatura_cm = estatura * 100) 
##   genero estatura estatura_cm
## 1  mujer     1.65         165
## 2 hombre     1.80         180
## 3  mujer     1.70         170
## 4  mujer     1.60         160
## 5 hombre     1.67         167
mutate(df_ej, estatura_cm = estatura * 100, estatura_in = estatura_cm * 0.3937) 
##   genero estatura estatura_cm estatura_in
## 1  mujer     1.65         165     64.9605
## 2 hombre     1.80         180     70.8660
## 3  mujer     1.70         170     66.9290
## 4  mujer     1.60         160     62.9920
## 5 hombre     1.67         167     65.7479

Ejercicio

Summarise y resúmenes por grupo

Summarise sirve para crear nuevas bases de datos con resúmenes o agregaciones de los datos originales.

summarise(df_ej, promedio = mean(estatura))
##   promedio
## 1    1.684

Podemos hacer resúmenes por grupo, primero creamos una base de datos agrupada:

by_genero <- group_by(df_ej, genero)
by_genero
## Source: local data frame [5 x 2]
## Groups: genero [2]
## 
##   genero estatura
##   (fctr)    (dbl)
## 1  mujer     1.65
## 2 hombre     1.80
## 3  mujer     1.70
## 4  mujer     1.60
## 5 hombre     1.67

y después operamos sobre cada grupo, creando un resumen a nivel grupo y uniendo los subconjuntos en una base nueva:

summarise(by_genero, promedio = mean(estatura))
##   promedio
## 1    1.684

Ejercicio

Algunas funciones útiles con summarise son min(x), median(x), max(x), quantile(x, p), n(), sum(x), sum(x > 1), mean(x > 1), sd(x).

by_date <- flights %>%
  group_by(date)
no_miss <- by_date %>%
  filter(!is.na(dep))
delays <- no_miss %>%
  dplyr::summarise(
    mean_delay = mean(dep_delay),
    n = n())

Nota importante: existen conflictos con la función summarise. Por lo tanto, puede ser necesario indicar que se debe utilizar la función summarise del paquete dplyr. Esto se hace anteponiendo a la llamada de la función dos símbolos de ‘dos puntos’ :: y antecediendo el nombre del paquete que se desea utilizar, en este caso, dplyr. A veces se necesita hacer lo mismo con la función mutate. La razón de esto es que hay un conflicto entre los paquetes plyr y dplyr. En principio, plyr ya no se usa porque dplyr es la versión nueva (y mejor) de dplyr. Pero aquí usamos este paquete para hacer uso de las funciones laply, dlply, y ldply.

Operador pipeline

Cuando uno hace varias operaciones es difícil leer y entender el código:

hourly_delay <- filter(dplyr::summarise(group_by(filter(flights, !is.na(dep_delay)), date, hour), delay = mean(dep_delay), n = n()), n > 10)

La dificultad radica en que usualmente los parámetros se asignan después del nombre de la función usando (). El operador “Forward Pipe” (%>%) cambia este orden, manera que un parámetro que precede a la función es enviado (“piped”) a la función.

Veamos como cambia el código anterior:

hourly_delay <- flights %>%
  filter(!is.na(dep_delay)) %>%
  group_by(date, hour) %>%
  dplyr::summarise(delay = mean(dep_delay), n = n()) %>%
  filter(n > 10)

podemos leer %>% como “después”.

Ejercicio

Variables por grupo

En ocasiones es conveniente crear variables por grupo, por ejemplo estandarizar dentro de cada grupo z = (x - mean(x)) / sd(x).

Veamos un ejemplo:

planes <- flights %>%
  filter(!is.na(arr_delay)) %>%
  group_by(plane) %>%
  filter(n() > 30)

planes %>%
  mutate(z_delay = (arr_delay - mean(arr_delay)) / sd(arr_delay)) %>%
  filter(z_delay > 5)
## Source: local data frame [1,415 x 15]
## Groups: plane [774]
## 
##          date  hour minute   dep   arr dep_delay arr_delay carrier flight
##        (date) (int)  (int) (int) (int)     (int)     (int)   (chr)  (int)
## 1  2011-01-28    15     16  1516  1916       351       326      CO      1
## 2  2011-01-27    18     22  1822  1945       234       210      CO    137
## 3  2011-01-27    21     37  2137  2254       242       219      CO    150
## 4  2011-01-27    22     37  2237   153       227       208      CO    250
## 5  2011-01-27    11      2  1102  1621       127       181      CO    510
## 6  2011-01-27    21     28  2128   136       231       216      CO   1632
## 7  2011-01-26    13     25  1325  1816       160       181      CO    106
## 8  2011-01-26    11     46  1146  1633       171       193      CO    510
## 9  2011-01-26     9     49   949  1436       144       180      CO    776
## 10 2011-01-20     6     35   635   807       780       775      CO     59
## ..        ...   ...    ...   ...   ...       ...       ...     ...    ...
## Variables not shown: dest (chr), plane (chr), cancelled (int), time (int),
##   dist (int), z_delay (dbl)

Verbos de dos tablas

¿Cómo mostramos los retrasos de los vuelos en un mapa?

Para responder esta pregunta necesitamos unir la base de datos de vuelos con la de aeropuertos.

location <- airports %>%
  select(dest = iata, name = airport, lat, long)

flights %>%
  group_by(dest) %>%
  filter(!is.na(arr_delay)) %>%
  dplyr::summarise(
    arr_delay = mean(arr_delay),
    n = n() ) %>%
  arrange(desc(arr_delay)) %>%
  left_join(location) %>%
  tbl_df
## Joining by: "dest"
## Source: local data frame [116 x 6]
## 
##     dest arr_delay     n                                name      lat
##    (chr)     (dbl) (int)                               (chr)    (dbl)
## 1    ANC  26.08065   124 Ted Stevens Anchorage International 61.17432
## 2    CID  17.80049   406                       Eastern Iowa  41.88459
## 3    DSM  15.95110   634            Des Moines International 41.53493
## 4    SFO  14.89036  2800         San Francisco International 37.61900
## 5    BPT  14.33333     3            Southeast Texas Regional 29.95083
## 6    GRR  13.71729   665           Kent County International 42.88082
## 7    DAY  13.67117   444             James M Cox Dayton Intl 39.90238
## 8    VPS  12.45718   864                Eglin Air Force Base 30.48325
## 9    ECP  12.42222   720                                  NA       NA
## 10   SAV  12.33137   851              Savannah International 32.12758
## ..   ...       ...   ...                                 ...      ...
## Variables not shown: long (dbl)

Hay varias maneras de unir dos bases de datos y debemos pensar en el objetivo:

x <- data.frame(name = c("John", "Paul", "George", "Ringo", "Stuart", "Pete"),
                instrument =c("guitar", "bass", "guitar", "drums", "bass", "drums"))

y <- data.frame(name = c("John", "Paul", "George", "Ringo", "Brian"),
                band = c("TRUE", "TRUE", "TRUE",  "TRUE", "FALSE"))
x
##     name instrument
## 1   John     guitar
## 2   Paul       bass
## 3 George     guitar
## 4  Ringo      drums
## 5 Stuart       bass
## 6   Pete      drums
y
##     name  band
## 1   John  TRUE
## 2   Paul  TRUE
## 3 George  TRUE
## 4  Ringo  TRUE
## 5  Brian FALSE
inner_join(x, y)
## Joining by: "name"
## Warning in inner_join_impl(x, y, by$x, by$y): joining factors with
## different levels, coercing to character vector
##     name instrument band
## 1   John     guitar TRUE
## 2   Paul       bass TRUE
## 3 George     guitar TRUE
## 4  Ringo      drums TRUE
left_join(x, y)
## Joining by: "name"
## Warning in left_join_impl(x, y, by$x, by$y): joining factors with different
## levels, coercing to character vector
##     name instrument band
## 1   John     guitar TRUE
## 2   Paul       bass TRUE
## 3 George     guitar TRUE
## 4  Ringo      drums TRUE
## 5 Stuart       bass <NA>
## 6   Pete      drums <NA>
semi_join(x, y)
## Joining by: "name"
## Warning in semi_join_impl(x, y, by$x, by$y): joining factors with different
## levels, coercing to character vector
##     name instrument
## 1   John     guitar
## 2   Paul       bass
## 3 George     guitar
## 4  Ringo      drums
anti_join(x, y)
## Joining by: "name"
## Warning in anti_join_impl(x, y, by$x, by$y): joining factors with different
## levels, coercing to character vector
##     name instrument
## 1   Pete      drums
## 2 Stuart       bass

Resumamos lo que observamos arriba:

Tipo Acción
inner Incluye únicamente las filas que aparecen tanto en x como en y
left Incluye todas las filas en x y las filas de y que coincidan
semi Incluye las filas de x que coincidan con y
anti Incluye las filas de x que no coinciden con y

Ahora combinamos datos a nivel hora con condiciones climáticas, ¿cuál es el tipo de unión adecuado?

hourly_delay <- flights %>%
  group_by(date, hour) %>%
  filter(!is.na(dep_delay)) %>%
  dplyr::summarise(
    delay = mean(dep_delay),
    n = n() ) %>%
  filter(n > 10)

delay_weather <- hourly_delay %>% left_join(weather)
## Joining by: c("date", "hour")
arrange(delay_weather, -delay)
## Source: local data frame [5,759 x 16]
## Groups: date [365]
## 
##          date  hour    delay     n  temp dew_point humidity pressure
##        (date) (int)    (dbl) (int) (dbl)     (dbl)    (int)    (dbl)
## 1  2011-05-12    23 184.2121    33  75.9      73.9       94    29.91
## 2  2011-09-29    22 166.2647    34  84.9      71.1       63    29.85
## 3  2011-09-29    23 164.5455    11  82.9      71.1       67    29.86
## 4  2011-05-12    22 162.7955    44  75.9      73.9       94    29.88
## 5  2011-05-20    23 151.9167    12  79.0      73.9       84    29.78
## 6  2011-05-12    21 145.5250    40  75.9      73.9       94    29.88
## 7  2011-05-12    20 144.9737    38  75.9      73.9       94    29.86
## 8  2011-12-31    20 141.0909    11  53.1      50.0       89    30.05
## 9  2011-05-12    19 130.2143    28  77.0      75.2       94    29.83
## 10 2011-06-22    23 126.8182    11  84.0      75.9       76    29.86
## ..        ...   ...      ...   ...   ...       ...      ...      ...
## Variables not shown: visibility (dbl), wind_dir (chr), wind_dir2 (int),
##   wind_speed (dbl), gust_speed (dbl), precip (dbl), conditions (chr),
##   events (chr)

Ejercicio

Métodos más generales (Do)

A diferencia de summarise, do es una función más lenta pero más general, resulta particularmente útil cuando se usa para ajustar un conjunto de modelos.

Veamos un ejemplo con la base de datos de bateo Lahman. Agruparemos por año y ajustaremos un modelo para explorar la relación entre número de bateos y carreras.

batting <- read.csv("data/batting.csv")

batting_center <- batting %>% 
  group_by(yearID) %>%
  mutate(
    AB_center = AB - mean(AB, na.rm = TRUE) # centramos dentro de cada grupo (año)
  ) 

models <- batting_center %>%
  group_by(yearID) %>%
  do(
    mod = lm(formula = R ~ AB_center, data = .)
  )
models
## Source: local data frame [142 x 2]
## Groups: <by row>
## 
##    yearID     mod
##     (int)   (chr)
## 1    1871 <S3:lm>
## 2    1872 <S3:lm>
## 3    1873 <S3:lm>
## 4    1874 <S3:lm>
## 5    1875 <S3:lm>
## 6    1876 <S3:lm>
## 7    1877 <S3:lm>
## 8    1878 <S3:lm>
## 9    1879 <S3:lm>
## 10   1880 <S3:lm>
## ..    ...     ...
models[[1, 'mod']]
## 
## Call:
## lm(formula = R ~ AB_center, data = .)
## 
## Coefficients:
## (Intercept)    AB_center  
##     39.3348       0.2688

Podemos extraer componentes de los objetos lm guerdados en la variable mod.

models_r <- models %>% 
  dplyr::summarise(r.squared = summary(mod)$r.squared) 

models_2 <- models %>% 
  group_by(yearID) %>% 
  do(data.frame(t(coef(.[[1, 'mod']])))) %>%
  mutate(intercept = X.Intercept.) %>%
  select(-X.Intercept.)
## Warning: Grouping rowwise data frame strips rowwise nature
models_3 <- cbind(models_2, models_r)


ggplot(models_3, aes(x = yearID, y = r.squared)) + geom_point() + geom_smooth()

Veamos las rectas de regresión.

ggplot(batting, aes(x = AB, y = R, color = yearID, group = yearID)) +
  geom_smooth(method = "lm")

Podemos comparar los coeficientes

models_3 <- gather(models_2, coef, value, -yearID)

ggplot(models_3, aes(x = yearID, y = value)) +
  geom_line() + 
  facet_wrap(~ coef, scales = "free_y", nrow = 2)

Podemos ver las diferencias a nivel año en los datos crudos.

ggplot(batting, aes(x = AB, y = R, color = yearID)) + 
  geom_point(alpha = 0.8, size = 1.6) 
## Warning: Removed 6413 rows containing missing values (geom_point).

Es conveniente usar el paquete biglm cuando se ajustan muchos modelos lineales, ya que guarda objetos de menor tamaño.

library(biglm)
## Loading required package: DBI
models_big <- batting %>%
  group_by(yearID) %>%
  do(
    mod = biglm(R ~ AB, data = .)
  )

print(object.size(models), unit = "MB")
## 22.7 Mb
print(object.size(models_big), unit = "MB")
## 0.8 Mb

En la función do puede ser un poco complicado extraer las componentes de los objetos que ajustamos. Una alternativa es usar las funciones del paquete plyr (en particular ddply, ldply) o las funciones que vienen en R estándar (lapply, apply, tapply).

models <- dlply(batting_center, "yearID", function(df){ 
  lm(R ~ AB_center, data = df)     # regresión lineal
})

models[[1]]
## 
## Call:
## lm(formula = R ~ AB_center, data = df)
## 
## Coefficients:
## (Intercept)    AB_center  
##     39.3348       0.2688
coef_list <- lapply(models, coef)
head(do.call(rbind, coef_list))
##      (Intercept) AB_center
## 1871    39.33482 0.2687665
## 1872    35.18026 0.2494250
## 1873    33.12953 0.2409013
## 1874    28.02950 0.2046628
## 1875    25.46118 0.1926905
## 1876    23.37078 0.1728847

Recursos adicionales