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.
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)
Encuentra todos los vuelos hacia SFO ó OAK.
Los vuelos con un retraso mayor a una hora.
En los que el retraso de llegada es más del doble que el retraso de salida.
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"))
Ve la ayuda de select (?select
) y escribe tres maneras de seleccionar las variables de retraso (delay).
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
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 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
Calcula la velocidad en millas por hora a partir de la variable tiempo y la distancia (en millas). ¿Quá vuelo fue el más rápido?
Crea una nueva variable que muestre cuánto tiempo se ganó o perdió durante el vuelo.
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
Calcula el retraso promedio por fecha.
¿Qué otros resúmenes puedes hacer para explorar el retraso por fecha?
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.
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”.
¿Qué destinos tienen el promedio de retrasos más alto?
¿Qué vuelos (compañía + vuelo) ocurren diario?
En promedio, ¿cómo varían a lo largo del día los retrasos de vuelos no cancelados? (pista: hour + minute / 60)
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)
¿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)
¿Qué condiciones climáticas están asociadad con retrasos en las salidas de Houston?
Explora si los aviones más viejos están asociados a mayores retrasos, responde con una gráfica.
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
Tidy Data, Hadley Wickham.
The Slit-Apply-Combine Strategy for Data Analysis, Hadley Wickham 2011.
Data Wrangling Cheat Sheet, RStudio.