---
title: "2D Quantitative Data"
subtitle: "Lecture 5"
author: "Shannon Gallagher"
date: "September 9, 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
source("_setup.R")
library(hexbin)

hail <- events |>
  filter(event_type == "Hail", magnitude > 0, prop_usd > 0) |>
  mutate(diameter = magnitude)

hail_model <- lm(log10(prop_usd) ~ diameter, data = hail)
slope <- coef(summary(hail_model))["diameter", ]
slope_ci <- confint(hail_model)["diameter", ]
slope_estimate <- unname(slope["Estimate"])
slope_se <- unname(slope["Std. Error"])
dollar_multiplier <- 10^slope_estimate
dollar_multiplier_se <- log(10) * dollar_multiplier * slope_se

duration_damage <- events |>
  filter(duration_hours > 0, duration_hours <= 168, prop_usd > 0)
```

## Logistics and final-project reminders

:::: {.columns}
::: {.column width="49%"}
**Due dates and feedback**

- **Checkpoint:** tonight, Wednesday, September 9, at **11:59 p.m. ET**
- **HW 3:** Tuesday, September 15, at **11:59 p.m. ET**
- **HW 1 grades:** expected in Gradescope by the end of Friday, September 11
:::

::: {.column width="51%"}
**Final project**

- Your teammate assignment should now be available.
- Keep the **two rubrics** separate: Data Visualization assesses the written report;
  Professional Skills assesses the companion presentation.
- Optional inspiration: [NOAA Shiny demo](../../noaa-dashboard.html). Shiny can help
  communicate an idea, but it is **not required** and we will not cover it explicitly.
:::
::::

## Reminders, previously, and today...

- make scatterplots that reveal rather than conceal,
- use linear and smooth trend lines carefully,
- read residual plots as diagnostic graphics,
- replace overplotted points with bins or contours.

## 2D quantitative data

For hail reports with positive damage:

> How does the distribution of reported property damage change with reported maximum
> hailstone diameter?

- A **marginal distribution** describes one variable by itself: all hailstone diameters,
  or all positive damage reports.
- A **conditional distribution** describes damage among reports at a given diameter:
  $Y\mid X=x$.
- The scatterplot shows the **joint distribution**; a regression will summarize one
  feature of the conditional distribution, $E(Y\mid X=x)$.

## Making scatterplots with `geom_point()`

:::: {.columns}
::: {.column width="46%"}
```{r hail-scatter-code, eval=FALSE}
hail |>
  ggplot(aes(diameter, prop_usd)) +
  geom_point(
    alpha = 0.2,
    color = deep_gold
  ) +
  geom_smooth(
    method = "lm",
    formula = y ~ sqrt(x),
    se = FALSE,
    color = blue_gray
  ) +
  scale_y_log10(labels = label_dollar())
```
:::

::: {.column width="54%"}
```{r hail-scatter}
#| echo: false
#| fig-width: 5.4
#| fig-height: 4.5
#| fig-alt: "Scatterplot of maximum hailstone diameter and positive property damage, with a square-root trend."
hail |>
  ggplot(aes(diameter, prop_usd)) +
  geom_point(alpha = 0.2, color = deep_gold) +
  geom_smooth(
    method = "lm", formula = y ~ sqrt(x), se = FALSE,
    color = blue_gray, linewidth = 1
  ) +
  scale_y_log10(labels = label_dollar()) +
  labs(
    title = "Damage tends to rise with diameter",
    subtitle = "Poisson-like spread; square-root trend in blue",
    x = "Reported maximum hailstone diameter (inches)",
    y = "Reported property damage"
  )
```
:::
::::

## Rugs connect the joint and marginal distributions

```{r hail-rug}
#| echo: false
#| fig-alt: "Scatterplot of hailstone diameter and positive property damage with rugs along the top and right margins."
hail |>
  ggplot(aes(diameter, prop_usd)) +
  geom_point(alpha = 0.18, color = deep_gold) +
  geom_rug(sides = "tr", alpha = 0.25, color = blue_gray) +
  scale_y_log10(labels = label_dollar()) +
  labs(
    title = "The rugs show each variable's marginal distribution",
    subtitle = "A vertical slice of points approximates damage conditional on diameter",
    x = "Reported maximum hailstone diameter (inches)",
    y = "Reported property damage"
  )
```

- **Top rug:** the marginal distribution of reported diameter, $X$.
- **Right rug:** the marginal distribution of positive damage, $Y$.
- **Vertical point cloud at $x$:** the conditional distribution, $Y\mid X=x$.

## Checkpoint 1 · Audit the scatterplot

:::: {.columns}
::: {.column width="54%"}
```{r hail-checkpoint-plot}
#| echo: false
#| fig-width: 5.2
#| fig-height: 4
#| fig-alt: "Scatterplot of hailstone diameter and positive property damage with one-inch and two-inch vertical slices highlighted."
hail |>
  ggplot(aes(diameter, prop_usd)) +
  geom_point(alpha = 0.18, color = deep_gold) +
  geom_vline(xintercept = c(1, 2), color = blue_gray, linetype = "dashed") +
  geom_rug(sides = "tr", alpha = 0.2, color = blue_gray) +
  scale_y_log10(labels = label_dollar()) +
  labs(
    x = "Maximum hailstone diameter (inches)",
    y = "Reported property damage"
  )
```
:::

::: {.column width="46%"}
::: {.checkpoint}
Work in three steps:

1. Use the **rugs** to describe each marginal distribution.
2. Compare the conditional damage distributions near **1 inch** and **2 inches**.
3. Identify one design choice that reveals structure, one that could mislead, and one
   omitted variable that could explain the vertical spread.
:::
:::
::::

## Displaying trend lines: linear regression

**Measurement note:** NOAA's hail magnitude is the reported maximum hailstone
**diameter** in inches—not radius or accumulation.
([NOAA recording guidance](https://www.ncei.noaa.gov/stormevents/pd01016005curr.pdf))

:::: {.columns}
::: {.column width="46%"}
```{r hail-lm-code, eval=FALSE}
hail |>
  ggplot(aes(diameter, log10(prop_usd))) +
  geom_point(
    alpha = 0.15,
    color = blue_gray
  ) +
  geom_smooth(
    method = "lm",
    se = TRUE,
    color = deep_gold
  )
```
:::

::: {.column width="54%"}
```{r hail-lm}
#| echo: false
#| fig-width: 5.4
#| fig-height: 4.3
#| fig-alt: "Scatterplot of maximum hailstone diameter and log property damage with a fitted regression line."
hail |>
  ggplot(aes(diameter, log10(prop_usd))) +
  geom_point(alpha = 0.15, color = blue_gray) +
  geom_smooth(method = "lm", se = TRUE, color = deep_gold) +
  labs(
    title = "Expected log damage rises with diameter",
    subtitle = "Line = fitted conditional mean; band = 95% CI",
    x = "Reported maximum hailstone diameter (inches)",
    y = "log10(reported property damage)"
  )
```
:::
::::

## Assessing assumptions of linear regression

For outcome $Y$ and predictor $X$:

$$
E(Y\mid X=x)=\beta_0+\beta_1x.
$$

- Here, $Y$ is log property damage and $X$ is reported maximum hailstone diameter.
- $\beta_1$ is the expected change in log damage for a one-inch difference in diameter.
  It describes association, not causation. See your **Linear Models** course for the
  estimation and inferential details.

## Interpreting the slope on two scales

:::: {.columns}
::: {.column width="50%"}
**Modeled scale: log10 dollars**

For reports whose maximum hailstone diameter differs by one inch, we estimate expected
`log10(property damage)` to be **`r sprintf("%.3f", slope_estimate)` higher**
(SE **`r sprintf("%.3f", slope_se)`**; 95% CI
**`r sprintf("%.3f", slope_ci[1])`–`r sprintf("%.3f", slope_ci[2])`**).
:::

::: {.column width="50%"}
**Original scale: dollars**

Under the regression model, expected property damage is
**`r sprintf("%.2f", dollar_multiplier)` times as large** (approx. SE
**`r sprintf("%.2f", dollar_multiplier_se)`**; 95% CI
**`r sprintf("%.2f", 10^slope_ci[1])`–`r sprintf("%.2f", 10^slope_ci[2])`**).
:::
::::

## What the regression model assumes

For $Y_i=\beta_0+\beta_1X_i+\epsilon_i$, we assume:

- $E(\epsilon_i\mid X_i)=0$: errors have conditional mean zero,
- errors are independent across observations,
- **homoskedasticity**: $Var(\epsilon_i\mid X_i)=\sigma^2$,
- approximately normal residuals for the usual small-sample inference,
- a constant error distribution for retransformation to expected dollars.

See your **Linear Models** course for formal diagnostics and the consequences of each
assumption.

## Residual vs fit plots

```{r hail-residuals, fig.alt="Residual versus fitted plot for a linear regression of log damage on maximum hailstone diameter."}
tibble(fitted = fitted(hail_model), residual = resid(hail_model)) |>
  ggplot(aes(fitted, residual)) +
  geom_hline(yintercept = 0, color = "gray60") +
  geom_point(alpha = 0.2, color = deep_gold) +
  geom_smooth(method = "loess", se = FALSE, color = blue_gray) +
  labs(
    title = "Residual structure tests the adequacy of the line",
    x = "Fitted value",
    y = "Residual"
  )
```

## A useful residual plot looks uneventful

Look for:

- curvature, suggesting the mean function is wrong;
- a funnel, suggesting nonconstant spread;
- clusters, suggesting omitted groups;
- isolated points with unusual residuals;
- changing patterns across the fitted range.

Ask: **Are the linear-model assumptions plausible?** An uneventful plot supports
plausibility; strong structure is evidence against the model.

## Locally weighted smoothing with LOESS

`geom_smooth(method = "lm")` answers: what linear trend summarizes the data?

`geom_smooth(method = "loess")` answers: what flexible local pattern appears?

R's `loess()` really does use the unusual **tricubic** distance weight. At each target
value $x_0$, the default fit uses the nearest 75% of observations and weights observation
$i$ by

$$
w_i(x_0)=\left[1-\left(\frac{d_i}{d_{\max}}\right)^3\right]^3.
$$

It then fits a weighted local polynomial—degree 2 by default. A degree-1 fit is local
linear regression. ([R `loess()` documentation](https://stat.ethz.ch/R-manual/R-devel/library/stats/html/loess.html))

## Displaying trend lines: LOESS

:::: {.columns}
::: {.column width="44%"}
```{r hail-line-loess-code, eval=FALSE}
hail |>
  ggplot(aes(
    diameter,
    sqrt(prop_usd)
  )) +
  geom_point(
    alpha = 0.1,
    color = "gray55"
  ) +
  geom_smooth(
    aes(color = "Linear"),
    method = "lm", se = FALSE
  ) +
  geom_smooth(
    aes(color = "LOESS"),
    method = "loess", se = FALSE
  )
```
:::

::: {.column width="56%"}
```{r hail-line-loess}
#| echo: false
#| fig-width: 5.6
#| fig-height: 4.6
#| fig-alt: "Maximum hailstone diameter and square-root property damage with linear and LOESS trend lines."
hail |>
  ggplot(aes(diameter, sqrt(prop_usd))) +
  geom_point(alpha = 0.1, color = "gray55") +
  geom_smooth(aes(color = "Linear"), method = "lm", se = FALSE) +
  geom_smooth(aes(color = "LOESS"), method = "loess", se = FALSE) +
  scale_color_manual(
    values = c("Linear" = deep_gold, "LOESS" = blue_gray)
  ) +
  labs(
    title = "Square-root damage still shows curvature",
    subtitle = "Compare the linear summary with the local LOESS fit",
    x = "Reported maximum hailstone diameter (inches)",
    y = "sqrt(property damage)",
    color = NULL
  )
```
:::
::::

## Comparing linear and curved models

```r
m_linear <- lm(sqrt(prop_usd) ~ diameter, data = hail)
m_quadratic <- lm(sqrt(prop_usd) ~ poly(diameter, 2), data = hail)

anova(m_linear, m_quadratic)  # nested-model test
AIC(m_linear, m_quadratic)    # fit with a complexity penalty
BIC(m_linear, m_quadratic)    # stronger complexity penalty
```

- **ANOVA** asks whether the added polynomial term improves a nested model.
- **AIC/BIC** trade goodness of fit against complexity; lower is preferred.
- Residual plots reveal *where* a candidate still misses structure.

These are complementary evidence, not automatic decision rules.

## Checkpoint 2 · Line or smoother?

::: {.checkpoint}
Use the preceding trend and residual plots:

1. Is a straight-line mean relationship plausible?
2. Would the curved model add an important feature or mostly complexity?
3. Defend your choice with one visual diagnostic and one comparison tool: ANOVA, AIC,
   or BIC.
:::

## What about focusing on the joint distribution?

With thousands of observations, a scatterplot can become cluttered, and individual
points can distract from where the data concentrate.

Alternatives:

- two-dimensional bins,
- hexagonal bins,
- density contours,
- small multiples,
- a representative sample plus a summary.

## Visualizing grid heat maps

`bins = 20` requests about **20 intervals in each direction**—up to roughly
$20\times20$ rectangular cells, not 20 cells total.

:::: {.columns}
::: {.column width="45%"}
```{r duration-damage-bin-code, eval=FALSE}
duration_damage |>
  ggplot(aes(
    log10(duration_hours),
    log10(prop_usd)
  )) +
  geom_bin_2d(bins = 20) +
  scale_fill_gradient(
    low = pale_gold,
    high = deep_gold,
    trans = "log10"
  )
```
:::

::: {.column width="55%"}
```{r duration-damage-bin}
#| echo: false
#| fig-width: 5.5
#| fig-height: 4.3
#| fig-alt: "Two-dimensional binned heatmap of event duration and positive property damage."
duration_damage |>
  ggplot(aes(log10(duration_hours), log10(prop_usd))) +
  geom_bin_2d(bins = 20) +
  scale_fill_gradient(low = pale_gold, high = deep_gold, trans = "log10") +
  labs(
    title = "Binning shows where records concentrate",
    x = "log10(duration in hours)",
    y = "log10(property damage)",
    fill = "Records"
  )
```
:::
::::

## Alternative idea: hexagonal binning

Here, `bins = 20` similarly controls the number of hexagons **across each direction**,
not the total number of hexagons drawn.

:::: {.columns}
::: {.column width="45%"}
```{r duration-damage-hex-code, eval=FALSE}
duration_damage |>
  ggplot(aes(
    log10(duration_hours),
    log10(prop_usd)
  )) +
  geom_hex(bins = 20) +
  scale_fill_gradient(
    low = pale_gold,
    high = deep_gold,
    trans = "log10"
  )
```
:::

::: {.column width="55%"}
```{r duration-damage-hex}
#| echo: false
#| fig-width: 5.5
#| fig-height: 4.3
#| fig-alt: "Hexagonal binned plot of duration and positive property damage."
duration_damage |>
  ggplot(aes(log10(duration_hours), log10(prop_usd))) +
  geom_hex(bins = 20) +
  scale_fill_gradient(low = pale_gold, high = deep_gold, trans = "log10") +
  labs(
    title = "Hexagons show local concentration",
    x = "log10(duration in hours)",
    y = "log10(property damage)",
    fill = "Records"
  )
```
:::
::::

## Going from 1D to 2D density estimation

A two-dimensional kernel density estimate places a smooth surface over the plane. A
contour line connects locations with the same estimated density, much as a topographic
line connects locations with the same elevation.

With one shared bandwidth $h$, a simplified estimator is

$$
\widehat f_h(x,y)=\frac{1}{nh^2}\sum_{i=1}^n
K\!\left(\frac{x-x_i}{h},\frac{y-y_i}{h}\right).
$$

The $h^2$ appears because smoothing now occurs in two directions. Smaller $h$ creates
more local detail; larger $h$ produces a smoother surface.

## So how do we display densities for 2D data?

```{r duration-damage-contours, fig.alt="Density contours over duration and positive property damage."}
duration_damage |>
  ggplot(aes(log10(duration_hours), log10(prop_usd))) +
  geom_point(alpha = 0.05, color = "gray55") +
  geom_density_2d(color = deep_gold, linewidth = 0.8) +
  labs(
    title = "Contours summarize where observations concentrate",
    x = "log10(duration in hours)",
    y = "log10(property damage)"
  )
```

## How to read contour plots?

A density contour connects locations with equal estimated density.

:::: {.columns}
::: {.column width="50%"}
**It does show:**

- high-density peaks,
- clusters or multiple modes,
- orientation and shape of the joint pattern,
- steep density gradients where lines are close.
:::

::: {.column width="50%"}
**It does not show:**

- a regression line,
- equal probability between contours,
- causality,
- exact counts unless levels are labeled appropriately.
:::
::::

## Checkpoint 3 · Match graph to task

::: {.checkpoint}
Choose scatterplot, fitted line, residual plot, or binned heatmap for each task:

1. find individual outliers;
2. assess whether a line is adequate;
3. show where thousands of points concentrate;
4. summarize average association.
:::

## Checkpoint 3 answers · inspect individual records

:::: {.columns}
::: {.column width="50%"}
```{r checkpoint3-scatter-answer}
#| echo: false
#| fig-width: 5
#| fig-height: 4
#| fig-alt: "Scatterplot for finding individual observations that depart from the main hail-damage pattern."
hail |>
  ggplot(aes(diameter, log10(prop_usd))) +
  geom_point(alpha = 0.25, color = deep_gold) +
  labs(
    title = "1. Find individual outliers",
    subtitle = "Use a scatterplot",
    x = "Maximum hailstone diameter (inches)",
    y = "log10(property damage)"
  )
```
:::

::: {.column width="50%"}
```{r checkpoint3-residual-answer}
#| echo: false
#| fig-width: 5
#| fig-height: 4
#| fig-alt: "Residual-versus-fitted plot for assessing whether the fitted line is adequate."
tibble(fitted = fitted(hail_model), residual = resid(hail_model)) |>
  ggplot(aes(fitted, residual)) +
  geom_hline(yintercept = 0, color = "gray60") +
  geom_point(alpha = 0.25, color = deep_gold) +
  geom_smooth(method = "loess", se = FALSE, color = blue_gray) +
  labs(
    title = "2. Assess line adequacy",
    subtitle = "Use a residual plot",
    x = "Fitted value",
    y = "Residual"
  )
```
:::
::::

## Checkpoint 3 answers · summarize the pattern

:::: {.columns}
::: {.column width="50%"}
```{r checkpoint3-bin-answer}
#| echo: false
#| fig-width: 5
#| fig-height: 4
#| fig-alt: "Binned heatmap for showing where many event-duration and property-damage observations concentrate."
duration_damage |>
  ggplot(aes(log10(duration_hours), log10(prop_usd))) +
  geom_bin_2d(bins = 20) +
  scale_fill_gradient(low = pale_gold, high = deep_gold, trans = "log10") +
  labs(
    title = "3. Show concentration",
    subtitle = "Use a binned heatmap",
    x = "log10(duration in hours)",
    y = "log10(property damage)",
    fill = "Records"
  )
```
:::

::: {.column width="50%"}
```{r checkpoint3-line-answer}
#| echo: false
#| fig-width: 5
#| fig-height: 4
#| fig-alt: "Scatterplot with fitted regression line for summarizing average association between hailstone diameter and log property damage."
hail |>
  ggplot(aes(diameter, log10(prop_usd))) +
  geom_point(alpha = 0.12, color = blue_gray) +
  geom_smooth(method = "lm", se = TRUE, color = deep_gold) +
  labs(
    title = "4. Summarize average association",
    subtitle = "Use a fitted line",
    x = "Maximum hailstone diameter (inches)",
    y = "log10(property damage)"
  )
```
:::
::::

## Recap and next steps

- Transparency and scale choices determine what a scatterplot reveals.
- Trend lines are statistical summaries with assumptions.
- Residuals are a visual model check.
- Bins and contours answer density questions when points overplot.

**Complete today's checkpoint by 11:59 p.m. ET tonight, Wednesday, September 9.**

Next: seeing structure across many variables.