overlapping_points.Rmd

```{r echo = FALSE, message = FALSE}
# run setup script
source("_common.R")

library(lubridate)
```

# Handling overlapping points {#overlapping-points}

When we want to visualize large or very large datasets, we often experience the challenge that simple *x*--*y* scatter plots do not work very well because many points lie on top of each other and partially or fully overlap. And similar problems can arise even in small datasets if data values were recorded with low precision or rounded, such that multiple observations have exactly the same numeric values. The technical term commonly used to describe this situation is "overplotting", i.e., plotting many points on top of each other. Here I describe several strategies you can pursue when encountering this challenge.


## Partial transparency and jittering

We first consider a scenario with only a moderate number of data points but with extensive rounding. Our dataset contains fuel economy during city driving and engine displacement for 234 popular car models released between 1999 and 2008 (Figure \@ref(fig:mpg-cty-displ-solid)). In this dataset, fuel economy is measured in miles per gallon (mpg) and is rounded to the nearest integer value. Engine displacement is measured in liters and is rounded to the nearest deciliter. Due to this rounding, many car models have exactly identical values. For example, there are 21 cars total with 2.0 liter engine displacement, and as a group they have only four different fuel economy values, 19, 20, 21, or 22 mpg. Therefore, in Figure \@ref(fig:mpg-cty-displ-solid) these 21 cars are represented by only four distinct points, so that 2.0 liter engines appear much less popular than they actually are.

(ref:mpg-cty-displ-solid) City fuel economy versus engine displacement, for popular cars released between 1999 and 2008. Each point represents one car. The point color encodes the drive train: front-wheel drive (FWD), rear-wheel drive (RWD), or four-wheel drive (4WD). The figure is labeled "bad" because many points are plotted on top of others and obscure them. 

```{r mpg-cty-displ-solid, fig.width=5.5, fig.asp=.7416, fig.cap='(ref:mpg-cty-displ-solid)'}
p_mpg_solid <- ggplot(mpg, aes(y = cty, x = displ, color = drv, fill = drv)) +
  geom_point(size = 3, shape = 21) + 
  ylab("fuel economy (mpg)") +
  xlab("displacement (l)") +
  scale_color_manual(values=c("#202020", "#E69F00", "#56B4E9"), 
                     name="drive train",
                     breaks=c("f", "r", "4"),
                     labels=c("FWD", "RWD", "4WD")) +
  scale_fill_manual(values=c("#202020", "#E69F00", "#56B4E9"), 
                     name="drive train",
                     breaks=c("f", "r", "4"),
                     labels=c("FWD", "RWD", "4WD")) +
  theme_half_open() +
  theme(legend.position = c(.7, .8))

stamp_bad(p_mpg_solid)
```

One way to ameliorate this problem is to use partial transparency. If we make individual points partially transparent, then overplotted points appear as darker points and thus the shade of the points reflects the density of points in that location of the graph (Figure \@ref(fig:mpg-cty-displ-transp)).

(ref:mpg-cty-displ-transp) City fuel economy versus engine displacement. Because points have been made partially transparent, points that lie on top of other points can now be identified by their darker shade. 

```{r mpg-cty-displ-transp, fig.width=5.5, fig.asp=.7416, fig.cap='(ref:mpg-cty-displ-transp)'}
p_mpg_transp <- ggplot(mpg, aes(y = cty, x = displ, color = drv, fill = drv)) +
  geom_point(size = 3, shape = 21) + 
  ylab("fuel economy (mpg)") +
  xlab("displacement (l)") +
  scale_color_manual(values=c("#202020", "#E69F00", "#56B4E9"), 
                     name="drive train",
                     breaks=c("f", "r", "4"),
                     labels=c("FWD", "RWD", "4WD")) +
  scale_fill_manual(values=c("#20202080", "#E69F0080", "#56B4E980"), 
                     name="drive train",
                     breaks=c("f", "r", "4"),
                     labels=c("FWD", "RWD", "4WD")) +
  theme_half_open() +
  theme(legend.position = c(.7, .8))

stamp_phantom(p_mpg_transp)
```

However, making points partially transparent is not always sufficient to solve the issue of overplotting. For example, even though we can see in Figure \@ref(fig:mpg-cty-displ-transp) that some points have a darker shade than others, it is difficult to estimate how many points were plotted on top of each other in each location. In addition, while the differences in shading are clearly visible, they are not self-explanatory. A reader who sees this figure for the first time will likely wonder why some points are darker than others and will not realize that those points are in fact multiple points stacked on top of each other. A simple trick that helps in this situation is to apply a small amount of jitter to the points, i.e., to displace each point randomly by a small amount in either the *x* or the *y* direction or both. With jitter, it is immediately apparent that the darker areas arise from points that are plotted on top of each other (Figure \@ref(fig:mpg-cty-displ-jitter)).

(ref:mpg-cty-displ-jitter) City fuel economy versus engine displacement. By adding a small amount of jitter to each point, we can make the overplotted points more clearly visible without substantially distorting the message of the plot.

```{r mpg-cty-displ-jitter, fig.width=5.5, fig.asp=.7416, fig.cap='(ref:mpg-cty-displ-jitter)'}
p_mpg_jitter <- ggplot(mpg, aes(y = cty, x = displ, color = drv, fill = drv)) +
  geom_point(size = 3, shape = 21,
             position = position_jitter(width = 0.01 * diff(range(mpg$displ)),
                                        height = 0.01 * diff(range(mpg$cty)),
                                        seed = 7384)) + 
  ylab("fuel economy (mpg)") +
  xlab("displacement (l)") +
  scale_color_manual(values=c("#202020", "#E69F00", "#56B4E9"), 
                     name="drive train",
                     breaks=c("f", "r", "4"),
                     labels=c("FWD", "RWD", "4WD")) +
  scale_fill_manual(values=c("#20202080", "#E69F0080", "#56B4E980"), 
                     name="drive train",
                     breaks=c("f", "r", "4"),
                     labels=c("FWD", "RWD", "4WD")) +
  theme_half_open() +
  theme(legend.position = c(.7, .8))

stamp_phantom(p_mpg_jitter)
```

One downside of jittering is that it does change the data and therefore has to be performed with care. If we jitter too much, we end up placing points in locations that are not representative of the underlying dataset. The result is a misleading visualization of the data. See Figure \@ref(fig:mpg-cty-displ-jitter-extreme) as an example.

(ref:mpg-cty-displ-jitter-extreme) City fuel economy versus engine displacement. By adding too much jitter to the points, we have created a visualization that does not accurately reflect the underlying dataset.

```{r mpg-cty-displ-jitter-extreme, fig.width=5.5, fig.asp=.7416, fig.cap='(ref:mpg-cty-displ-jitter-extreme)'}
p_mpg_jitter_extreme <- ggplot(mpg, aes(y = cty, x = displ, color = drv, fill = drv)) +
  geom_point(size = 3, shape = 21,
             position = position_jitter(width = 0.1 * diff(range(mpg$displ)),
                                        height = 0.1 * diff(range(mpg$cty)))) + 
  scale_x_continuous(breaks = 2:7) +
  ylab("fuel economy (mpg)") +
  xlab("displacement (l)") +
  scale_color_manual(values=c("#202020", "#E69F00", "#56B4E9"), 
                     name="drive train",
                     breaks=c("f", "r", "4"),
                     labels=c("FWD", "RWD", "4WD")) +
  scale_fill_manual(values=c("#20202080", "#E69F0080", "#56B4E980"), 
                     name="drive train",
                     breaks=c("f", "r", "4"),
                     labels=c("FWD", "RWD", "4WD")) +
  theme_half_open() +
  theme(legend.position = c(.7, .8))

stamp_bad(p_mpg_jitter_extreme)
```

## 2d histograms

When the number of individual points gets very large, partial transparency (with or without jittering) will not be sufficient to resolve the overplotting issue. What will typically happen is that areas with high point density will appear as uniform blobs of dark color while in areas with low point density the individual points are barely visible (Figure \@ref(fig:nycflights-points)). And changing the transparency level of individual points will either ameliorate one or the other of these problems while worsening the other; no transparency setting can address both at the same time. 

(ref:nycflights-points) Departure delay in minutes versus the flight departure time, for all flights departing Newark airport (EWR) in 2013. Each dot represents one departure.

```{r nycflights-points, fig.asp = 0.75, fig.cap = '(ref:nycflights-points)'}
# break points along the x axis
breaks_x <- c("0:00", "6:00", "12:00", "18:00", "24:00")

p_flights_base <- ggplot(flight_delays, aes(`departure time`, `departure delay (minutes)`)) + 
  geom_abline(slope = 0, intercept = 0, color="grey80") +
  scale_x_time(breaks = hm(breaks_x),
               labels = breaks_x) +
  theme_half_open()
  
p_flights_scatter <- p_flights_base + geom_point(alpha = 0.2)

stamp_bad(p_flights_scatter)
```

Figure \@ref(fig:nycflights-points) shows departure delays for over 100,000 individual flights, with each dot representing one flight departure. Even though we have made the individual dots fairly transparent, the majority of them just forms a black band between 0 and 300 minutes departure delay. This band obscures whether most flights depart approximately on time or with substantial delay (say 50 minutes or more). At the same time, the most delayed flights (with delays of 400 minutes or more) are barely visible due to the transparency of the dots.

In such cases, instead of plotting individual points, we can make a 2d histogram, where we subdivide the entire *x*--*y* plane into small rectangles, count how many observations fall into each rectangles, and then color the rectangles by that count. Figure \@ref(fig:nycflights-2d-bins) shows the result of this approach for the departure-delay data. This visualization clearly highlights several important features of the flight-departure data. First, the vast majority of departures during the day (6am to about 9pm) actually depart without delay or even early (negative delay). However, a modest number of departures has a substantial delay. Moreover, the later a plane departs in the day the more of a delay it can have. Importantly, the departure time is the actual time of departure, not the scheduled time of departure. So this figure does not necessarily tell us that planes scheduled to depart early never experience delay. What it does tell us, though, is that if a plane departs early it either has little delay or, in very rare cases, a delay of around 900 minutes.

(ref:nycflights-2d-bins) Departure delay in minutes versus the flight departure time. Each colored rectangle represents all flights departing at that time with that departure delay. Coloring represents the number of flights represented by that rectangle.

```{r nycflights-2d-bins, fig.asp = 0.75, fig.cap = '(ref:nycflights-2d-bins)'}
p_flights_2d_bins <- p_flights_base +
  geom_bin2d(bins=50) +
  scale_fill_continuous_sequential(palette = "Blue-Yellow", l2 = 90, c2 = 20) +
  theme(legend.position = c(0.85, .85))
stamp_phantom(p_flights_2d_bins)
```


As an alternative to binning the data into rectangle, we can also bin into hexagons. This approach, first proposed by 
@Carr-et-al-1987, has the advantage that the points in a hexagon are, on average, closer to the hexagon center than the points in an equal-area square are to the center of the square. Therefore, the colored hexagon represents the data slightly more accurately than the colored rectangle does. Figure \@ref(fig:nycflights-hex-bins) shows the flight departure data with hexagon binning rather than rectangular binning.

(ref:nycflights-hex-bins) Departure delay in minutes versus the flight departure time. Each colored hexagon represents all flights departing at that time with that departure delay. Coloring represents the number of flights represented by that hexagon.

```{r nycflights-hex-bins, fig.asp = 0.75, fig.cap = '(ref:nycflights-hex-bins)'}
p_flights_hex_bins <- p_flights_base +
  geom_hex(bins=50) +
  scale_fill_continuous_sequential(palette = "Blue-Yellow", l2 = 90, c2 = 20) +
  theme(legend.position = c(0.85, .85))
stamp_phantom(p_flights_hex_bins)
```


## Contour lines

Instead of binning data points into rectangles or hexagons, we can also estimate the point density across the plot area and indicate regions of different point densities with contour lines. This technique works well when the point density changes slowly across both the *x* and the *y* dimensions.

As an example for this approach, we consider the relationship between population number and area for counties in the Midwest. We have data for 1055 counties, and a scatter plot looks like a cloud of points (Figure \@ref(fig:midwest-scatter)). We can highlight the distribution of points more clearly by making them very small and partially transparent and ploting them on top of contour lines that delineate regions of comparable point density (Figure \@ref(fig:midwest-density-dots)). We can also plot just the contour lines, without the individual points (Figure \@ref(fig:midwest-density-smooth)). In this case, it can be helpful to add a trendline that shows the overall trend in the data. Here, there isn't much of a trend, and the shape of the trendline (approximately flat) reflects this lack of a trend.

(ref:midwest-scatter) Population versus area for counties in midwestern states. Data are taken from the 2010 US census and are shown for 1055 counties covering 12 states. Each dot represents one county. 

```{r midwest-scatter, fig.asp = 0.8, fig.cap = '(ref:midwest-scatter)'}
county_midwest <- left_join(US_census, US_regions) %>%
  filter(region == "Midwest")

p_midwest_base <- ggplot(county_midwest, aes(area, pop2010)) + 
  scale_fill_gradient(low = "grey70", high = "grey30", guide = "none") +
  scale_y_log10(breaks = c(1000, 1e4, 1e5, 1e6, 1e7),
                limits = c(450, 5.2e6),
                labels = label_log10) +
  scale_x_log10(breaks = c(100, 300, 1000, 3000),
                limits = c(50, 7000), expand = c(0, 0)) +
  ylab("population total") +
  xlab("area (square miles)") +
  theme_minimal_grid()

p_midwest_base + geom_point(color = "navy", size = 1, alpha = .7)
```


(ref:midwest-density-dots) Population versus area for counties in midwestern states. Contour lines and shaded areas indicate the density of counties for that combination of population total and area. Individual counties are shown as light blue dots.

```{r midwest-density-dots, fig.asp = 0.8, fig.cap = '(ref:midwest-density-dots)'}
p_midwest_base + 
  stat_density_2d(aes(fill = ..level..), geom = "polygon", color = "black", size = 0.2, alpha = 0.5) +
  geom_point(color = "navy", size = .5, alpha = .4)
```

(ref:midwest-density-smooth) Population versus area for counties in midwestern states. Contour lines and shaded areas indicate the density of counties for that combination of population total and area. Note that some counties lie outside the largest shaded area. The solid blue line highlights the mean relationship between population total and county area. It was obtained via least-square fitting of a general additive model with cubic spline base to the underlying data.

```{r midwest-density-smooth, fig.asp = 0.8, fig.cap = '(ref:midwest-density-smooth)'}
p_midwest_base + 
  stat_density_2d(aes(fill = ..level..), geom = "polygon", color = "black", size = 0.2, alpha = 0.5) +
  geom_smooth(color = "navy", se = FALSE, size = 0.5)
```