---
title: "Visualizing Quantitative Distributions"
subtitle: "Lecture 4"
author: "Shannon Gallagher"
date: "September 2, 2026"
date-format: long
format:
  revealjs:
    theme: ../mads-36613-assignments/lecture_drafts/theme.scss
    slide-number: c/t
    chalkboard: true
    code-fold: show
    code-line-numbers: true
    smaller: true
    linestretch: 1.25
    fig-width: 8
    fig-height: 4.5
    footer: "MaDS Data Visualization · Fall 2026"
execute:
  echo: true
  warning: false
  message: false
---

```{r setup}
#| include: false
library(readr)
library(dplyr)
library(forcats)
library(ggplot2)
library(scales)

theme_set(theme_minimal(base_size = 14))

gold <- "#d99a00"
deep_gold <- "#9b6b00"
pale_gold <- "#fff3c4"
charcoal <- "#2f2a1f"
blue_gray <- "#52646d"

local_data <- "../data/noaa_storm_events_2024_clean.csv.gz"
public_data <- paste0(
  "https://stat.cmu.edu/~sgallagh/courses/mads-fall-2026/",
  "datavis-36613/data/noaa_storm_events_2024_clean.csv.gz"
)
data_source <- if (file.exists(local_data)) local_data else public_data

events <- read_csv(data_source, show_col_types = FALSE)
positive_damage <- events |> filter(prop_usd > 0)
```

## Reminders, previously, and today...

> **Graded checkpoint:** Due **Wednesday, September 2 at 11:59 p.m. ET** on Canvas.

- see how histogram choices affect a story,
- read an empirical cumulative distribution function,
- understand what a density estimate smooths,
- compare conditional distributions without hiding their shape.

## Revisit histograms

The visible shape depends on:

- bin width,
- bin origin or boundary,
- scale or transformation,
- inclusion rules,
- sample size.

Defaults are suggestions, not findings.

## What happens as we change the bin width?

```{r histogram-bins, fig.alt="Two histograms of log property damage using different bin widths."}
#| layout-ncol: 2
#| output-location: slide
p1 <- positive_damage |>
  ggplot(aes(log10(prop_usd))) +
  geom_histogram(binwidth = 0.1, fill = gold, color = "white") +
  labs(title = "Bin width 0.1", x = NULL, y = NULL)

p2 <- positive_damage |>
  ggplot(aes(log10(prop_usd))) +
  geom_histogram(binwidth = 0.6, fill = blue_gray, color = "white") +
  labs(title = "Bin width 0.6", x = NULL, y = NULL)

p1
p2
```

## A feature should survive reasonable bin widths

```{r histogram-grid, fig.alt="Four histograms of the same NOAA damage data using different bin widths."}
#| layout-ncol: 2
#| output-location: slide
for (width in c(0.05, 0.15, 0.35, 0.75)) {
  print(
    positive_damage |>
      ggplot(aes(log10(prop_usd))) +
      geom_histogram(binwidth = width, fill = gold, color = "white") +
      labs(title = paste("Bin width", width), x = NULL, y = NULL)
  )
}
```

Stable conclusions should not depend on one convenient bin choice.

## Data or expertise can guide bin width

There are two useful starting points:

- **Numerical selection:** use a data-driven rule or algorithm, much like bandwidth selection for a kernel density estimate.
- **Subject-matter selection:** choose widths or boundaries that make meaningful differences easy to see.

On a $\log_{10}$ scale:

- a width of **1** compares orders of magnitude (each step is $\times 10$),
- a width of **0.1** compares roughly 26% multiplicative steps ($10^{0.1} \approx 1.26$).

Treat any choice as a starting point, then check whether the story survives nearby widths.

## What happens with a different sample?

```{r histogram-sample-variation, fig.alt="Histograms from four random samples of positive NOAA damage reports."}
#| layout-ncol: 2
#| output-location: slide
set.seed(36613)
for (sample_id in 1:4) {
  print(
    positive_damage |>
      slice_sample(n = min(500, nrow(positive_damage))) |>
      ggplot(aes(log10(prop_usd))) +
      geom_histogram(binwidth = 0.25, fill = blue_gray, color = "white") +
      labs(title = paste("Sample", sample_id), x = NULL, y = NULL)
  )
}
```

Sampling variability affects the picture even when the plotting rule stays fixed.

## Checkpoint 1 · What survives?

::: {.checkpoint}
Compare the two histograms. List one feature that appears in both and one feature that
depends on the binning. Which statement would you be willing to put in a report?
:::

## Display the full distribution with an ECDF plot

The empirical cumulative distribution function at value $x$ is:

$$
\widehat{F}(x) = \frac{\#\{X_i \le x\}}{n}
$$

Read it as: the share of observations at or below $x$.

## Build the fitted-normal reference

```{r ecdf-normal-reference}
log_damage <- log10(positive_damage$prop_usd)
normal_mean <- mean(log_damage)
normal_sd <- sd(log_damage)

normal_reference <- tibble(
  log_damage = seq(min(log_damage), max(log_damage), length.out = 400)
) |>
  mutate(
    prop_usd = 10^log_damage,
    cdf = pnorm(
      log_damage,
      mean = normal_mean,
      sd = normal_sd
    )
  )
```

## Check both sides of each ECDF jump

```{r ecdf-ks-candidates}
sorted_damage <- sort(log_damage)
n_damage <- length(sorted_damage)
fitted_normal_cdf <- pnorm(
  sorted_damage,
  mean = normal_mean,
  sd = normal_sd
)

ks_candidates <- bind_rows(
  tibble(
    log_damage = sorted_damage,
    empirical_cdf = seq_len(n_damage) / n_damage,
    fitted_cdf = fitted_normal_cdf
  ),
  tibble(
    log_damage = sorted_damage,
    empirical_cdf = (seq_len(n_damage) - 1) / n_damage,
    fitted_cdf = fitted_normal_cdf
  )
)
```

## Keep the largest gap for annotation

```{r ecdf-ks-mark}
ks_mark <- ks_candidates |>
  mutate(gap = abs(empirical_cdf - fitted_cdf)) |>
  slice_max(gap, n = 1, with_ties = FALSE) |>
  mutate(
    prop_usd = 10^log_damage,
    label_x = prop_usd * 1.35,
    label_y = (empirical_cdf + fitted_cdf) / 2,
    label = paste0("KS gap: D = ", number(gap, accuracy = 0.001))
  ) |>
  select(prop_usd, empirical_cdf, fitted_cdf, label_x, label_y, label)
```

## An ECDF answers threshold questions

```{r damage-ecdf, fig.alt="Empirical CDF of positive property damage and a fitted normal CDF on a logarithmic dollar axis. A red vertical segment at 5,000 dollars marks the maximum gap, the Kolmogorov-Smirnov statistic D equals 0.128."}
#| output-location: slide
ggplot() +
  stat_ecdf(
    data = positive_damage,
    aes(prop_usd, color = "Fn: empirical"),
    geom = "step",
    linewidth = 1
  ) +
  geom_line(
    data = normal_reference,
    aes(prop_usd, cdf, color = "F0: fitted normal"),
    linewidth = 1,
    linetype = "dashed"
  ) +
  geom_segment(
    data = ks_mark,
    aes(
      x = prop_usd,
      xend = prop_usd,
      y = empirical_cdf,
      yend = fitted_cdf
    ),
    inherit.aes = FALSE,
    color = "#b31b34",
    linewidth = 1.6
  ) +
  geom_point(
    data = ks_mark,
    aes(x = prop_usd, y = empirical_cdf),
    inherit.aes = FALSE,
    color = "#b31b34",
    size = 2.7
  ) +
  geom_point(
    data = ks_mark,
    aes(x = prop_usd, y = fitted_cdf),
    inherit.aes = FALSE,
    color = "#b31b34",
    size = 2.7
  ) +
  geom_text(
    data = ks_mark,
    aes(x = label_x, y = label_y, label = label),
    inherit.aes = FALSE,
    color = "#b31b34",
    fontface = "bold",
    hjust = 0,
    size = 4.2
  ) +
  scale_x_log10(
    breaks = c(1e3, 1e6, 1e9),
    labels = c("$1K", "$1M", "$1B")
  ) +
  scale_y_continuous(labels = percent) +
  scale_color_manual(
    values = c("Fn: empirical" = deep_gold, "F0: fitted normal" = blue_gray),
    breaks = c("Fn: empirical", "F0: fitted normal")
  ) +
  labs(
    title = "The largest fitted-normal gap occurs at $5,000",
    subtitle = "Normal fit to log10(property damage) using the sample mean and SD",
    x = "Reported property damage",
    y = "Share at or below amount",
    color = NULL
  )
```

## One-Sample Kolmogorov–Smirnov Test

The one-sample KS statistic is the largest vertical gap between the empirical and
reference cumulative distributions:

$$
D_n = \sup_x |F_n(x)-F_0(x)|.
$$

- $F_n$: the empirical CDF from the observed sample.
- $F_0$: a fixed reference CDF, or one informed by independent data.

Ideally, $F_0$ is specified in advance or estimated independently of the data behind
$F_n$. If the same data inform both curves—as they do here—treat the KS comparison as
an approximation.

## Boxplots compress aggressively

A boxplot shows median, quartiles, and a conventional outlier rule.

It does not show:

- multimodality,
- gaps,
- detailed tail behavior,
- every observation.

Use it for compact comparison, not as a complete distribution.

## Boxplots enable compact comparisons

```{r damage-boxplot, fig.alt="Boxplots of positive damage for five event types on a log scale."}
#| output-location: slide
compare_types <- positive_damage |>
  filter(event_type %in% c("Tornado", "Flash Flood", "Hail", "High Wind", "Flood"))

compare_types |>
  ggplot(aes(prop_usd, fct_reorder(event_type, prop_usd, median))) +
  geom_boxplot(fill = pale_gold, outlier.alpha = 0.15) +
  scale_x_log10(labels = label_dollar()) +
  labs(
    title = "Positive damage distributions differ by event type",
    x = "Reported property damage",
    y = NULL
  )
```

## Faceted histograms preserve each group's shape

```{r conditional-histograms, fig.alt="Faceted histograms of log damage for five NOAA event types."}
#| output-location: slide
compare_types |>
  ggplot(aes(log10(prop_usd))) +
  geom_histogram(binwidth = 0.3, fill = gold, color = "white") +
  facet_wrap(~ event_type, ncol = 2, scales = "free_y") +
  labs(title = "Small multiples separate shape from overlap", x = "log10(property damage)", y = "Reports")
```

## Density curves estimate the PDF

```{r conditional-density-lines, fig.alt="Density curves of log damage for selected NOAA event types."}
#| output-location: slide
compare_types |>
  ggplot(aes(log10(prop_usd), color = event_type)) +
  geom_density(linewidth = 1) +
  labs(
    title = "Each KDE estimates a PDF with total area one",
    x = "log10(property damage)",
    y = "Estimated density",
    color = "Event type"
  )
```

The curves compare estimated shapes. Because each integrates to one, height does not
encode the number of observations in a group.

## Checkpoint 2 · What vanished?

::: {.checkpoint}
Choose one event type in the boxplot. Name two facts the plot supports and two features
you would need a histogram, ECDF, or raw points to inspect.
:::

## Bandwidth controls the bias–variance tradeoff

A kernel density estimate places a small smooth bump at every observation and adds them.

Bandwidth determines how much those bumps overlap:

- **Smaller bandwidth:** lower smoothing bias, higher sampling variance; noise can look meaningful.
- **Larger bandwidth:** higher smoothing bias, lower sampling variance; real structure can disappear.

No single bandwidth reveals the truth. Look for conclusions that survive a reasonable range.

## Let R choose a starting bandwidth

```{r default-density, fig.alt="Kernel density estimate of log property damage using the default bandwidth selected by R's density function."}
#| output-location: slide
#| fig-height: 3.8
#| fig-cap: '[How `bw.nrd0()` calculates the default bandwidth](https://stat.ethz.ch/R-manual/R-patched/library/stats/html/bandwidth.html)'
default_kde <- density(log_damage)

default_kde_data <- tibble(
  log_damage = default_kde$x,
  estimated_pdf = default_kde$y
)

default_kde_data |>
  ggplot(aes(log_damage, estimated_pdf)) +
  geom_line(color = deep_gold, linewidth = 1.2) +
  labs(
    title = paste0(
      "density() chose bandwidth ",
      number(default_kde$bw, accuracy = 0.001)
    ),
    subtitle = 'R uses bw = "nrd0" by default',
    x = "log10(property damage)",
    y = "Estimated density"
  )
```

::: {.notes}
[Sources]

- [R: Bandwidth Selectors for Kernel Density Estimation](https://stat.ethz.ch/R-manual/R-patched/library/stats/html/bandwidth.html)
:::

## Bandwidth changes the apparent structure

```{r density-bandwidth, fig.alt="Density estimates of log property damage using two bandwidth adjustments."}
#| output-location: slide
positive_damage |>
  sample_n(min(10000, nrow(positive_damage))) |>
  ggplot(aes(log10(prop_usd))) +
  geom_density(aes(color = "Less smoothing"), adjust = 0.35, linewidth = 1) +
  geom_density(aes(color = "More smoothing"), adjust = 1.8, linewidth = 1) +
  scale_color_manual(values = c("Less smoothing" = deep_gold, "More smoothing" = blue_gray)) +
  labs(title = "Bandwidth changes which features look real", x = "log10(property damage)", y = "Density", color = NULL)
```

## We should NOT fill overlapping density curves

Multiple filled densities create an occlusion problem. Prefer:

- lines for two or three groups,
- small multiples,
- violins plus raw points,
- ECDFs when thresholds matter.

The comparison task should determine the display.

## Violin plots show smoothed conditional shape

```{r conditional-violins, fig.alt="Violin plots of log damage for selected NOAA event types."}
#| output-location: slide
compare_types |>
  ggplot(aes(fct_reorder(event_type, prop_usd, median), log10(prop_usd))) +
  geom_violin(fill = pale_gold, color = deep_gold, trim = FALSE) +
  geom_boxplot(width = 0.12, outlier.shape = NA, fill = "white") +
  coord_flip() +
  labs(title = "Violins show smoothed shape; the box marks robust summaries", x = NULL, y = "log10(property damage)")
```

## Raw points reveal within-group variation

```{r conditional-jitter, fig.alt="Jittered points and boxplots of sampled log damage by NOAA event type."}
#| output-location: slide
compare_types |>
  group_by(event_type) |>
  slice_sample(prop = 1) |>
  slice_head(n = 250) |>
  ungroup() |>
  ggplot(aes(fct_reorder(event_type, prop_usd, median), log10(prop_usd))) +
  geom_jitter(width = 0.15, alpha = 0.22, color = blue_gray) +
  geom_boxplot(width = 0.18, outlier.shape = NA, fill = NA, color = charcoal) +
  coord_flip() +
  labs(title = "A sample of raw points shows variation that summaries compress", x = NULL, y = "log10(property damage)")
```

## Checkpoint 3 · Choose the display

::: {.checkpoint}
Choose one question:

1. What share of flood records reported less than $50,000?
2. Which event types have the widest positive-damage distribution?
3. Is there evidence of multiple damage regimes?

Pick ECDF, boxplot, histogram, or density and defend the choice.
:::

## Recap and next steps

- Histogram features must survive reasonable bin choices.
- ECDFs are excellent for thresholds and comparisons.
- KDEs estimate PDFs; bandwidth trades bias against variance.
- Conditional displays need an explicit comparison task.

**After class:** Complete the graded checkpoint on Canvas by **Wednesday, September 2
at 11:59 p.m. Eastern.**

Next: relationships between two quantitative variables.