---
title: "t-SNE, UMAP, and Visualizing Trends + Time Series Data"
subtitle: "Lecture 8"
author: "Shannon Gallagher"
date: "September 21, 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")
purple <- "#6A3D9A"
embedding_metrics <- state_metrics |>
  select(
    state = state_upper,
    event_records,
    injuries,
    deaths,
    damage = property_damage,
    duration = mean_duration
  ) |>
  mutate(across(c(injuries, deaths, damage), log1p))
embedding_x <- embedding_metrics |> select(-state) |> scale()
```

## Today

:::: {.columns}
::: {.column width="38%"}
### Logistics

- **Homework 4:** due Tuesday, September 22 at 11:59 p.m.
- **Final project:** keep working on it and bring questions or blockers to class.
:::

::: {.column width="62%"}
### Topics

- distinguish global and local structure in an embedding,
- compare the objectives and settings of t-SNE and UMAP,
- explain why repeated runs and settings matter,
- begin visualizing ordered observations as a time series.
:::
::::

## Linear and nonlinear objectives

PCA asks for linear directions that preserve variance.

t-SNE asks for a low-dimensional arrangement that preserves local neighborhoods, using
probabilities rather than raw distance.

These are different objectives. Neither is universally “better.”

## A spiral exposes the linear limitation

```{r spiral-data, fig.alt="Two-dimensional spiral colored by position along the curve."}
set.seed(36613)
spiral <- tibble(
  t = seq(0, 4 * pi, length.out = 500),
  x = t * cos(t) + rnorm(500, sd = 0.12),
  y = t * sin(t) + rnorm(500, sd = 0.12)
)

spiral |>
  ggplot(aes(x, y, color = t)) +
  geom_point(size = 1.6) +
  scale_color_gradient(low = purple, high = deep_gold) +
  coord_equal() +
  labs(title = "Nearby points follow a curved manifold", color = "Position")
```

## The spiral needs both principal components

```{r spiral-pca-code}
spiral_pca <- prcomp(spiral |> select(x, y), scale. = TRUE)
spiral_variance <- spiral_pca$sdev^2 / sum(spiral_pca$sdev^2)

set.seed(36613)
pca_scores <- as_tibble(spiral_pca$x[, 1, drop = FALSE]) |>
  mutate(t = spiral$t, y = 0) |>
  slice_sample(prop = 1)
```

The curve has one intrinsic coordinate, `t`, but it bends through two linear directions.
PC1 explains **`r percent(spiral_variance[1], accuracy = 1)`** of the variance and PC2
explains **`r percent(spiral_variance[2], accuracy = 1)`**.

## Keeping PC1 alone loses the spiral

```{r spiral-pca, echo=FALSE, fig.alt="The first principal-component score for each point on the spiral, displayed along one horizontal dimension. Purple and gold points from different parts of the spiral overlap."}
pca_scores |>
  ggplot(aes(PC1, y, color = t)) +
  geom_point(size = 1.8, alpha = 0.65) +
  scale_color_gradient(low = purple, high = deep_gold) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(
    title = "One linear axis collapses distant parts of the curve",
    subtitle = paste0(
      "PC1 keeps ", percent(spiral_variance[1], accuracy = 1),
      " and drops ", percent(spiral_variance[2], accuracy = 1)
    ),
    color = "Position"
  )
```

```{r spiral-distance-setup}
#| include: false
spiral_sample_rows <- unique(round(seq(1, nrow(spiral), length.out = 50)))
spiral_sample <- spiral[spiral_sample_rows, ]
spiral_sample_x <- scale(
  spiral_sample |> select(x, y),
  center = spiral_pca$center,
  scale = spiral_pca$scale
)
spiral_sample_pc1 <- predict(
  spiral_pca,
  newdata = spiral_sample |> select(x, y)
)[, 1, drop = FALSE]

set.seed(36613)
spiral_sample_tsne <- Rtsne::Rtsne(
  spiral_sample_x,
  dims = 1,
  perplexity = 10,
  check_duplicates = FALSE
)$Y

original_distance_matrix <- as.matrix(dist(spiral_sample_x))
pc1_distance_matrix <- as.matrix(dist(spiral_sample_pc1))
tsne_distance_matrix <- as.matrix(dist(spiral_sample_tsne))
pair_index <- which(upper.tri(original_distance_matrix), arr.ind = TRUE)

nearest_neighbor <- t(apply(
  original_distance_matrix,
  1,
  function(distance_row) rank(distance_row, ties.method = "first") <= 6
))

spiral_pair_distances <- tibble(
  original = original_distance_matrix[pair_index],
  pc1 = pc1_distance_matrix[pair_index],
  tsne = tsne_distance_matrix[pair_index],
  local_pair = nearest_neighbor[pair_index] |
    nearest_neighbor[cbind(pair_index[, 2], pair_index[, 1])]
)
```

## PC1 changes distances on the spiral

```{r spiral-pca-distance, echo=FALSE, fig.alt="Scatterplot comparing original standardized distances among 50 sampled spiral points with their distances after projection onto PC1. Many points fall well below the perfect-preservation diagonal."}
spiral_pair_distances |>
  ggplot(aes(original, pc1)) +
  geom_abline(slope = 1, intercept = 0, color = purple, linewidth = 0.9) +
  geom_point(alpha = 0.35, color = deep_gold) +
  coord_equal() +
  labs(
    title = "Projection onto PC1 shortens many pairwise distances",
    x = "Original standardized distance",
    y = "Distance after keeping PC1"
  )
```

Perfect preservation would place every pair on the purple line. The points far below it
represent relationships that collapse when we discard PC2.

## t-SNE prioritizes local proximity

- t-SNE gives the closest pairs the greatest weight.
- **Its central goal is to preserve local proximity:** neighbors in the original space
  should remain neighbors in two dimensions.
- The algorithm represents closeness through probabilities, not exact metric distances.
- Global separation and cluster area can change substantially.

## t-SNE distances on the spiral

```{r spiral-tsne-distance, echo=FALSE, fig.alt="Scatterplot comparing original standardized distances among 50 sampled spiral points with their distances in a one-dimensional t-SNE map. Pairs that were nearest neighbors in the original space are highlighted in purple and tend to remain close in the map."}
spiral_pair_distances |>
  ggplot(aes(original, tsne)) +
  geom_point(color = "gray75", alpha = 0.25) +
  geom_point(
    data = \(data) filter(data, local_pair),
    color = purple,
    alpha = 0.8
  ) +
  labs(
    title = "Original nearest neighbors usually remain close",
    subtitle = "Purple pairs include one of the five nearest neighbors of either point",
    x = "Original standardized distance",
    y = "Distance in the 1D t-SNE map"
  )
```

t-SNE does not target the purple diagonal from the PCA slide. Its vertical scale has no
original unit, and distant pairs can move substantially.


## t-distributed stochastic neighbor embedding

- exploring possible local groupings,
- separating nearby neighborhoods,
- creating hypotheses for later investigation.

It is not a reliable measurement of cluster size, inter-cluster distance, or global
geometry.

## Perplexity sets the neighborhood scale

Think of perplexity as the **effective number of neighbors** each point consults.

:::: {.columns}
::: {.column width="50%"}
### Lower perplexity

- narrower neighborhood
- emphasizes very local detail
- can produce small, fragmented groups
:::

::: {.column width="50%"}
### Higher perplexity

- broader neighborhood
- emphasizes larger-scale patterns
- can smooth over small local groups
:::
::::

For our 50 states, we will compare 5, 10, and 15. The useful scale depends on the number
of observations and the structure in the data.

## Two places to turn on caching

:::: {.columns}
::: {.column width="48%"}
### Every chunk in the document

```yaml
---
execute:
  cache: true
---
```
:::

::: {.column width="52%"}
### One expensive R chunk

```r
#| cache: true
#| cache.extra: !expr file.mtime("data.csv")

fit <- Rtsne::Rtsne(x)
```
:::
::::

For this deck, chunk-level caching targets the slow t-SNE fits. `cache.extra` forces a
new result when the input file changes.

## Fit t-SNE in R

```{r tsne-code, eval=FALSE}
set.seed(36613)

fit <- Rtsne::Rtsne(
  embedding_x,
  dims = 2,
  perplexity = 10,
  check_duplicates = FALSE
)
```

`fit$Y` is a 50 × 2 matrix: one row per state and one column per plotted dimension.

## Compare several settings on the NOAA state profiles

```{r tsne-settings, fig.alt="t-SNE embeddings of state-level NOAA profiles using several perplexity values."}
#| echo: false
#| layout-ncol: 3
#| cache: true
#| cache.extra: !expr file.mtime("../data/noaa_storm_events_2024_clean.csv.gz")
for (perplexity_value in c(5, 10, 15)) {
  set.seed(36613)
  fit <- Rtsne::Rtsne(
    embedding_x,
    dims = 2,
    perplexity = perplexity_value,
    check_duplicates = FALSE
  )
  print(
    tibble(state = embedding_metrics$state, x = fit$Y[, 1], y = fit$Y[, 2]) |>
      ggplot(aes(x, y, label = state)) +
      geom_text(color = deep_gold, check_overlap = TRUE) +
      labs(title = paste("Perplexity", perplexity_value), x = NULL, y = NULL)
  )
}
```

## Checkpoint 1 · Read with restraint

::: {.checkpoint}
Someone says, “Cluster A is twice as far from B as from C on the t-SNE plot.” Explain why
that claim is unsafe and propose a defensible statement instead.
:::

## Random starts change the embedding

t-SNE commonly begins from a random arrangement. Repeated runs can differ.

A pattern that disappears across reasonable seeds or settings is not strong evidence.
Stability checks belong beside the visualization.

## Compare repeated runs

```{r tsne-seeds, fig.alt="Three t-SNE embeddings using different random seeds."}
#| echo: false
#| layout-ncol: 3
#| cache: true
#| cache.extra: !expr file.mtime("../data/noaa_storm_events_2024_clean.csv.gz")
for (seed_value in c(1, 2, 3)) {
  set.seed(seed_value)
  fit <- Rtsne::Rtsne(
    embedding_x,
    dims = 2,
    perplexity = 10,
    check_duplicates = FALSE
  )
  print(
    tibble(state = embedding_metrics$state, x = fit$Y[, 1], y = fit$Y[, 2]) |>
      ggplot(aes(x, y, label = state)) +
      geom_text(color = blue_gray, check_overlap = TRUE) +
      labs(title = paste("Seed", seed_value), x = NULL, y = NULL)
  )
}
```

## UMAP also starts with local neighborhoods

- UMAP builds a weighted graph that connects each observation to nearby observations.
- It searches for a low-dimensional layout with similar local connections.
- UMAP often runs faster than t-SNE as the number of observations grows.
- The axes, orientation, and absolute distances still have no direct interpretation.

UMAP may retain more broad structure than t-SNE, but neither method guarantees a faithful
global map.

## UMAP settings

:::: {.columns}
::: {.column width="50%"}
### Neighborhood size

**`n_neighbors`**

- lower values emphasize smaller neighborhoods
- higher values consider broader structure
- similar role to perplexity, but not the same quantity
:::

::: {.column width="50%"}
### Visual spacing

**`min_dist`**

- lower values allow tighter groups
- higher values keep nearby points more spread out
- changes the visual packing, not the input data
:::
::::

## Fit UMAP in R

```{r umap-code, eval=FALSE}
#| cache: true
set.seed(36613)

umap_fit <- uwot::umap(
  embedding_x,
  n_neighbors = 10,
  min_dist = 0.1,
  n_components = 2,
  n_threads = 1
)
```

`umap_fit` is a 50 × 2 matrix. Install the package once with
`install.packages("uwot")`.

## t-SNE and UMAP on the same state profiles

```{r tsne-umap-comparison, fig.alt="Side-by-side t-SNE and UMAP embeddings of the same 50 standardized NOAA state profiles. States are labeled by two-letter abbreviation."}
#| echo: false
#| layout-ncol: 2
#| cache: true
#| cache.extra: !expr file.mtime("../data/noaa_storm_events_2024_clean.csv.gz")
set.seed(36613)
comparison_tsne <- Rtsne::Rtsne(
  embedding_x,
  dims = 2,
  perplexity = 10,
  check_duplicates = FALSE
)$Y

set.seed(36613)
comparison_umap <- uwot::umap(
  embedding_x,
  n_neighbors = 10,
  min_dist = 0.1,
  n_components = 2,
  n_threads = 1,
  verbose = FALSE
)

tibble(
  state = embedding_metrics$state,
  x = comparison_tsne[, 1],
  y = comparison_tsne[, 2]
) |>
  ggplot(aes(x, y, label = state)) +
  geom_text(color = deep_gold, check_overlap = TRUE) +
  labs(title = "t-SNE", x = NULL, y = NULL)

tibble(
  state = embedding_metrics$state,
  x = comparison_umap[, 1],
  y = comparison_umap[, 2]
) |>
  ggplot(aes(x, y, label = state)) +
  geom_text(color = purple, check_overlap = TRUE) +
  labs(title = "UMAP", x = NULL, y = NULL)
```

Look for neighborhoods that appear in both maps. Do not compare coordinates, axis scales,
or the size of empty spaces across the panels.

## A responsible nonlinear embedding workflow

1. Clean and scale variables deliberately.
2. Try several seeds and neighborhood settings.
3. Compare with PCA or the original variables.
4. Investigate apparent groups in the source data.
5. Label the embedding as exploratory.

## Checkpoint 2 · Stability plan

::: {.checkpoint}
Design a three-run stability check for an embedding of state impact profiles. What would
need to remain similar before you trusted an apparent group?
:::

## From neighborhoods to temporal order

An embedding uses position to show similarity among unordered observations. Time series
already have a meaningful order, and a line makes that order visible.

Before drawing it, identify:

- time unit,
- aggregation rule,
- missing intervals,
- whether the series measures events, observations, or people.

## Example: NOAA event reports by day

```{r daily-counts, fig.alt="Line chart of daily NOAA event record counts during 2024."}
daily <- events |>
  mutate(day = as.Date(begin_dt)) |>
  count(day) |>
  complete(day = seq(min(day), max(day), by = "day"), fill = list(n = 0))

daily |>
  ggplot(aes(day, n)) +
  geom_line(color = deep_gold, linewidth = 0.7) +
  scale_x_date(date_breaks = "2 months", date_labels = "%b") +
  labs(title = "Event reports arrive in bursts", subtitle = "Daily NOAA Storm Events records, 2024", x = NULL, y = "Event records")
```

## Add lines to emphasize order

```{r points-versus-lines, fig.alt="Daily NOAA counts displayed first as points and then as a line."}
#| layout-ncol: 2
daily |>
  ggplot(aes(day, n)) +
  geom_point(alpha = 0.35, color = blue_gray) +
  labs(title = "Points", x = NULL, y = "Event records")

daily |>
  ggplot(aes(day, n)) +
  geom_line(color = deep_gold) +
  labs(title = "A line emphasizes sequence", x = NULL, y = "Event records")
```

## Area charts emphasize volume

```{r daily-area, fig.alt="Area chart of daily NOAA event counts."}
daily |>
  ggplot(aes(day, n)) +
  geom_area(fill = pale_gold, color = deep_gold, linewidth = 0.35) +
  scale_x_date(date_breaks = "2 months", date_labels = "%b") +
  labs(
    title = "Area emphasizes total volume",
    subtitle = "Fine day-to-day variation becomes harder to see",
    x = NULL,
    y = "Event records"
  )
```

## Lines emphasize trends and sequence

Connecting daily counts helps reveal sequence and bursts. It does not imply that the
underlying hazard changes smoothly between days.

For sparse or irregular observations, points, steps, or explicit gaps may be more honest.

## Several time series

```{r multi-series, fig.alt="Monthly lines for four NOAA event types."}
events |>
  filter(event_type %in% c("Hail", "Tornado", "Heat", "Winter Weather")) |>
  count(month, event_type) |>
  complete(month = 1:12, event_type, fill = list(n = 0)) |>
  ggplot(aes(month, n, color = event_type)) +
  geom_line(linewidth = 1) +
  geom_point(size = 1.5) +
  scale_x_continuous(breaks = 1:12, labels = month.abb) +
  labs(title = "Different hazards have different seasonal signatures", x = NULL, y = "Event records", color = NULL) +
  theme(legend.position = "bottom")
```

## Directly label lines

```{r direct-series-labels, fig.alt="Monthly NOAA event-type lines with labels at their December endpoints."}
series <- events |>
  filter(event_type %in% c("Hail", "Tornado", "Heat", "Winter Weather")) |>
  count(month, event_type) |>
  complete(month = 1:12, event_type, fill = list(n = 0))

ggplot(series, aes(month, n, color = event_type)) +
  geom_line(linewidth = 1) +
  ggrepel::geom_text_repel(
    data = filter(series, month == 12),
    aes(label = event_type),
    direction = "y",
    hjust = 0,
    nudge_x = 0.35,
    segment.color = "gray70",
    min.segment.length = 0,
    max.overlaps = Inf,
    seed = 36613,
    show.legend = FALSE
  ) +
  scale_x_continuous(breaks = 1:12, labels = month.abb, limits = c(1, 14.5)) +
  labs(title = "Label lines where the reader finishes tracing them", x = NULL, y = "Event records", color = NULL) +
  theme(legend.position = "none")
```

## Highlighting one series among many

```{r many-state-lines, fig.alt="Monthly event counts for many states with Texas highlighted."}
state_series <- events |>
  mutate(state = state.abb[match(str_to_title(state_upper), state.name)]) |>
  filter(!is.na(state)) |>
  count(state, month) |>
  complete(state, month = 1:12, fill = list(n = 0))

ggplot(state_series, aes(month, n, group = state)) +
  geom_line(color = "gray82") +
  geom_line(data = filter(state_series, state == "TX"), color = deep_gold, linewidth = 1.4) +
  scale_x_continuous(breaks = 1:12, labels = month.abb) +
  labs(title = "Texas in context", subtitle = "Other states remain visible in gray", x = NULL, y = "Event records")
```

## Time axes and missing intervals

- Show the full relevant time window.
- Make missing intervals explicit.
- Avoid unequal spacing that visually implies equal time.
- Use zero when magnitude comparison requires it; explain justified truncation.
- Separate seasonal patterns from changes in coverage or reporting.

## Checkpoint 3 · What is the series?

::: {.checkpoint}
For the multi-series plot, write what one point represents. Then identify one reporting or
exposure issue that prevents the line from being interpreted as pure weather frequency.
:::

## Recap and next steps

- Nonlinear embeddings are exploratory neighborhood views.
- t-SNE and UMAP use different objectives and settings.
- Settings, seeds, and preprocessing are part of the result.
- Time adds order and makes missing intervals consequential.
- A line chart is only as meaningful as its aggregation rule.

Next: time-series structure and spatial foundations.