library(nycflights13) library(tidyverse) #> ── Attaching core tidyverse packages ──────────────── tidyverse 1.3.2.9000 ── #> ✔ dplyr 1.0.99.9000 ✔ readr 2.1.3 #> ✔ forcats 0.5.2.9000 ✔ stringr 1.5.0.9000 #> ✔ ggplot2 3.4.0.9000 ✔ tibble 3.1.8 #> ✔ lubridate 1.9.0 ✔ tidyr 1.2.1.9001 #> ✔ purrr 1.0.1 #> ── Conflicts ─────────────────────────────────────── tidyverse_conflicts() ── #> ✖ dplyr::filter() masks stats::filter() #> ✖ dplyr::lag() masks stats::lag() #> ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
Case study: aggregates and sample size
Whenever you do any aggregation, it’s always a good idea to include a count (n()). That way, you can ensure that you’re not drawing conclusions based on very small amounts of data. For example, let’s look at the planes (identified by their tail number) that have the highest average delays:
delays <- flights |>
filter(!is.na(arr_delay), !is.na(tailnum)) |>
group_by(tailnum) |>
summarize(
delay = mean(arr_delay, na.rm = TRUE),
n = n()
)
ggplot(delays, aes(x = delay)) +
geom_freqpoly(binwidth = 10)
Wow, there are some planes that have an average delay of 5 hours (300 minutes)! That seems pretty surprising, so lets draw a scatterplot of number of flights vs. average delay:
ggplot(delays, aes(x = n, y = delay)) + geom_point(alpha = 1/10)
Not surprisingly, there is much greater variation in the average delay when there are few flights for a given plane. The shape of this plot is very characteristic: whenever you plot a mean (or other summary) vs. group size, you’ll see that the variation decreases as the sample size increases*cough* the central limit theorem *cough*..
When looking at this sort of plot, it’s often useful to filter out the groups with the smallest numbers of observations, so you can see more of the pattern and less of the extreme variation in the smallest groups:
delays |> filter(n > 25) |> ggplot(aes(x = n, y = delay)) + geom_point(alpha = 1/10) + geom_smooth(se = FALSE)
Note the handy pattern for combining ggplot2 and dplyr. It’s a bit annoying that you have to switch from |> to +, but it’s not too much of a hassle once you get the hang of it.
There’s another common variation on this pattern that we can see in some data about baseball players. The following code uses data from the Lahman package to compare what proportion of times a player hits the ball vs. the number of attempts they take:
batters <- Lahman::Batting |>
group_by(playerID) |>
summarize(
perf = sum(H, na.rm = TRUE) / sum(AB, na.rm = TRUE),
n = sum(AB, na.rm = TRUE)
)
batters
#> # A tibble: 20,166 × 3
#> playerID perf n
#> <chr> <dbl> <int>
#> 1 aardsda01 0 4
#> 2 aaronha01 0.305 12364
#> 3 aaronto01 0.229 944
#> 4 aasedo01 0 5
#> 5 abadan01 0.0952 21
#> 6 abadfe01 0.111 9
#> # … with 20,160 more rows
When we plot the skill of the batter (measured by the batting average, ba) against the number of opportunities to hit the ball (measured by at bat, ab), you see two patterns:
As above, the variation in our aggregate decreases as we get more data points.
There’s a positive correlation between skill (
perf) and opportunities to hit the ball (n) because obviously teams want to give their best batters the most opportunities to hit the ball.
batters |>
filter(n > 100) |>
ggplot(aes(x = n, y = perf)) +
geom_point(alpha = 1 / 10) +
geom_smooth(se = FALSE)
This also has important implications for ranking. If you naively sort on desc(ba), the people with the best batting averages are clearly lucky, not skilled:
batters |> arrange(desc(perf)) #> # A tibble: 20,166 × 3 #> playerID perf n #> <chr> <dbl> <int> #> 1 abramge01 1 1 #> 2 alberan01 1 1 #> 3 banisje01 1 1 #> 4 bartocl01 1 1 #> 5 bassdo01 1 1 #> 6 birasst01 1 2 #> # … with 20,160 more rows
You can find a good explanation of this problem and how to overcome it at http://varianceexplained.org/r/empirical_bayes_baseball/ and https://www.evanmiller.org/how-not-to-sort-by-average-rating.html.