Quantile dotplots and HOPs of snowfall

Setup

The following libraries are needed:

library(tidyverse)
library(tidybayes)
library(gganimate)

Introduction

This is a silly example based on this tweet:

To generate a quantile dotplot from this, I am going to start by approximating an inverse cumulative distribution function, F_inv, for these predictions (aka a quantile function). I am doing this rather than just turning the above histogram directly into dots because the bin widths above are not equal, and because it will make a smoother-looking plot. This is probably overkill.

F_inv = splinefun(     # or use approxfun for linear interpolation
  c(0, .12, .49, .91, .99,  1),
  c(0,   3,   6,  10,  15, 20),
  method = "monoH.FC"
  )

The quantile function I’ve derived looks like this:

curve(F_inv(x), xlim = c(0,1))

Quantile dotplots

We can now use the inverse CDF to generate quantiles from the predictive distribution and construct quantile dotplots. E.g. a 50-dot dotplot:

tibble(
  p = ppoints(50),
  snowfall = F_inv(p)
) %>%
  ggplot(aes(x = snowfall)) +
  geom_dotplot(method = "histodot", binwidth = 1, origin = 0, dotsize = .96, stackratio = 1.03) +
  scale_y_continuous(breaks = NULL) +
  scale_x_continuous(breaks = seq(0, 18, by = 2)) +
  coord_cartesian(expand = FALSE, xlim = c(0,18)) +
  xlab("Predicted snowfall in inches. Each dot represents a 1/50\nchance of that level of snowfall.") +
  ylab(NULL) +
  theme_tidybayes() +
  theme(axis.title.x.bottom = element_text(hjust = 0), axis.line.x.bottom = element_line(color = "gray75"))

Or a 100-dot dotplot:

tibble(
  p = ppoints(100),
  snowfall = F_inv(p)
) %>%
  ggplot(aes(x = snowfall)) +
  geom_dotplot(method = "histodot", binwidth = 1, dotsize = .9, stackratio = 1.05, origin = 0) +
  scale_y_continuous(breaks = NULL) +
  scale_x_continuous(breaks = seq(0, 18, by = 2)) +
  coord_cartesian(expand = FALSE, xlim = c(0,18)) +
  xlab("Predicted snowfall in inches. Each dot represents a 1/100\nchance of that level of snowfall.") +
  ylab(NULL) +
  theme_tidybayes() +
  theme(axis.title.x.bottom = element_text(hjust = 0), axis.line.x.bottom = element_line(color = "gray75"))

Vince Baumel pointed out that this might look like a pile of snow, so why not push the metaphor with color?

tibble(
  p = ppoints(100),
  snowfall = F_inv(p)
) %>%
  ggplot(aes(x = snowfall)) +
  geom_dotplot(method = "histodot", binwidth = 1, dotsize = .9, stackratio = 1.05, origin = 0, fill = "white", color = NA) +
  scale_y_continuous(breaks = NULL) +
  scale_x_continuous(breaks = seq(0, 18, by = 2)) +
  coord_cartesian(expand = FALSE, xlim = c(0,18)) +
  xlab("Predicted snowfall in inches. Each dot represents a 1/100\nchance of that level of snowfall.") +
  ylab(NULL) +
  theme_tidybayes() +
  theme(
    axis.title.x.bottom = element_text(hjust = 0), 
    axis.line.x.bottom = element_line(color = "gray75"),
    panel.background = element_rect(fill = "skyblue")
  )

HOPs

We could also randomly re-order the quantiles to generate a hypothetical outcome plot:

anim = tibble(
  p = ppoints(100),
  snowfall = F_inv(p)
) %>%
  sample_n(100) %>%
  mutate(id = 1:n()) %>%
  ggplot() +
  geom_vline(aes(xintercept = snowfall)) +
  scale_y_continuous(breaks = NULL) +
  scale_x_continuous(breaks = seq(0, 18, by = 2)) +
  coord_cartesian(expand = FALSE, xlim = c(0,18)) +
  xlab("Predicted snowfall in inches. Each line is an approximately\nequally likely predicted outcome.") +
  ylab(NULL) +
  theme_tidybayes() +
  theme(axis.title.x.bottom = element_text(hjust = 0), axis.line.x.bottom = element_line(color = "gray75")) +
  transition_manual(id)

animate(anim, fps = 2.5, res = 100, width = 400, height = 200)

Or if we’re getting excessive, we can make a HOPs version with falling snowballs. We’ll also use the experimental ggfx package to blur the snowballs to make them look a bit more authentic:

k = 50
nframes = k * 10

anim = tibble(
    p = ppoints(k),
    snowfall = F_inv(p)
  ) %>%
  sample_n(k) %>%
  mutate(
    id = 1:n(),
    time = 1:n(),
    y = list(c(0,1))
  ) %>%
  unnest(y) %>%
  mutate(
    time = time + y
  ) %>%
  ggplot(aes(x = snowfall, y = 1 - y, group = id)) +
  ggfx::with_kernel(
    # sigma = 3,
    geom_point(size = 3, color = "white"),
    "Comet:0x10-90"
  ) +
  scale_y_continuous(breaks = NULL) +
  scale_x_continuous(breaks = seq(0, 18, by = 2)) +
  coord_cartesian(expand = FALSE, xlim = c(0,18)) +
  xlab("Predicted snowfall in inches. Each snowball is an\napproximately equally likely predicted outcome.") +
  ylab(NULL) +
  theme_tidybayes() +
  theme(
    axis.title.x.bottom = element_text(hjust = 0), 
    axis.line.x.bottom = element_line(color = "gray75"),
    panel.background = element_rect(fill = "skyblue")
  ) +
  transition_components(time = time, exit_length = .5) +
  ease_aes("quadratic-in") +
  exit_fade() +
  shadow_wake(wake_length = 0.25/k)

animate(anim, nframes = nframes, fps = nframes / k * 2.5, res = 150, width = 600, height = 300)