5  Density plots

TipTheme: theme_wired()

Graphs in this chapter use fmggp2::theme_wired(), built from the WIRED (US) November/December 2024 cover: a black background, a white frame and bold Helvetica Neue title, and red, blue, pink, gray, and steel blue marks.

Requirements:

✅ a numeric (continuous) variable

5.1 Description

A density plot displays data distribution using a smooth curve instead of bars. A histogram divides the x axis into discrete ‘bins’ and counts the values in each one:

Instead of bins, density plots transform the distribution according to a kernel density estimate. geom_density() places a small smooth curve (a bell shaped normal kernel by default) over each value and adds them together into one curve:

The area under the curve represents the total probability (it adds up to one). The curve is high where values are common and low where they are rare. Density plots also help compare multiple sets of data. Legends are used to explain each curve, and different colors are used to differentiate them.

5.2 Set up

PACKAGES:

Install packages.

library(ggplot2)
library(dplyr) # for tibble() and data manipulation

DATA:

WIRED has covered technology and culture since its first issue in 1993. Picture a month of WIRED stories: lots of short news posts, plus long features like the ones in the magazine’s Big Interview issue.

We’ll simulate the reading times of 400 stories in two steps. First, sample() decides whether each story is a news post or a feature. Then rnorm() draws a reading time from that format’s own normal distribution. Setting the seed with set.seed() means you’ll get the same 400 stories every time you run the code below.

set.seed(1993) # same 400 stories every time the code runs
wired_stories <- tibble(
  story_id = 1:400, # 400 different stories, one row each
  format = sample(
    c("news", "feature"), # each story is a news post or a feature
    size = 400, # one format for each story
    replace = TRUE, # formats repeat across stories
    prob = c(0.65, 0.35) # about two thirds are news posts
    ),
  minutes = round(
    if_else(format == "news",
      rnorm(400, mean = 6, sd = 1.5), # news posts center on 6 minutes
      rnorm(400, mean = 30, sd = 7) # features center on 30 minutes
      )
    )
)
glimpse(wired_stories)
#> Rows: 400
#> Columns: 3
#> $ story_id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, …
#> $ format   <chr> "news", "news", "news", "featur…
#> $ minutes  <dbl> 3, 8, 5, 43, 7, 7, 5, 6, 17, 5,…

5.3 Distribution

Every distribution we’ve seen so far has a single peak. Real data doesn’t always cooperate, and the WIRED stories are a good example: they have two peaks, one for news posts and one for features. Density plots are good at showing this kind of shape, but summary statistics can hide it completely. Below we’ll move from the reading times, to a picture, to the formula.

COUNT IT:

Start with the first ten stories.

head(wired_stories, 10)

The reading times fall into two clumps: news posts between 3 and 8 minutes, and features at 17 and 43 minutes. Nothing in this first ten takes between 9 and 16 minutes. Now look at what the mean and median say about all 400 stories.

wired_stories |> 
  summarise(
    mean = mean(minutes), 
    median = median(minutes), 
    stories_12_to_18_min = sum(minutes >= 12 & minutes <= 18))

The mean reading time is about 15 minutes, but only 8 of the 400 stories take between 12 and 18 minutes to read. The mean describes almost none of the stories.

DRAW IT:

A bar model shows the stories as two parts of one whole. News posts make up about two thirds of the bar.

wired_stories |> 
  count(format) |> 
  ggplot(aes(x = n, y = "400 stories", fill = format)) + 
  geom_col(position = position_stack(reverse = TRUE)) + 
  geom_text(aes(label = n), color = "white", 
    position = position_stack(vjust = 0.5, reverse = TRUE)) + 
  labs(x = "Number of stories", y = NULL, fill = "Format")

The density plot below draws each format’s reading times as its own curve. The news curve is tall and narrow (most news posts take close to 6 minutes), and the feature curve is short and wide (feature lengths vary more). The dashed line marks the mean reading time, which falls in the valley between the two peaks.

ggplot(data = wired_stories, 
    aes(x = minutes)) + 
  geom_density(aes(y = after_stat(count), fill = format), 
    alpha = 0.6, color = NA) + 
  geom_density(aes(y = after_stat(count)), fill = NA) + 
  geom_vline(xintercept = mean(wired_stories$minutes), linetype = "dashed") + 
  labs(x = "Reading time (minutes)", y = "Stories per minute", fill = "Format")

The solid line is what geom_density() draws for all 400 stories together. It finds the same two peaks, but it uses one amount of smoothing (the bandwidth) for every story, so the narrow news peak comes out lower and wider than the news curve. You can adjust the smoothing with the adjust argument (e.g., geom_density(adjust = 0.5)).

NAME IT:

Our two step simulation creates a mixture distribution. Each story’s reading time comes from one of two normal distributions, chosen with probabilities 0.65 and 0.35:

\[ f(x) = 0.65 \cdot \text{Normal}(x;\, 6, 1.5) + 0.35 \cdot \text{Normal}(x;\, 30, 7) \]

The weights (0.65 and 0.35) set how much area each peak gets, the means (6 and 30) set where the peaks sit, and the standard deviations (1.5 and 7) set how wide they are. The mean of the mixture is the weighted average of the two means: \(0.65 \times 6 + 0.35 \times 30 = 14.4\) minutes. That’s why the mean lands in the valley.

Try moving the two means closer together in rnorm(). Once the peaks overlap enough, the valley disappears and the curve looks like it has a single peak, even though the data still comes from two groups.

5.4 Grammar

CODE:

Create labels with labs()

Initialize the graph with ggplot() and provide data

Map minutes to the x axis

Add the geom_density() layer

labs_density <- labs(
  title = "WIRED",
  caption = "Simulated data",
  subtitle = "Reading times of 400 WIRED stories",
  x = "Reading time (minutes)")
ggp2_density <- ggplot(data = wired_stories,
        aes(x = minutes)) +
      geom_density()

ggp2_density +
  labs_density

GRAPH:

A downside of density plots is that the y axis is hard to interpret.

The y axis shows probability density, not a count of stories. Its units are 1 per minute, and the total area under the curve is one. That’s why the values are so small (the peak is under 0.07), and why they would change if we measured reading time in hours instead of minutes. We read a density plot by where the curve is high and where it’s low, so it’s common to hide the y axis with scale_y_continuous(NULL, breaks = NULL).

ggp2_density <- ggplot(data = wired_stories,
        aes(x = minutes)) +
      geom_density() +
      scale_y_continuous(NULL, breaks = NULL)

ggp2_density +
  labs_density