---
title: "High-Dimensional Data"
subtitle: "Lecture 7"
author: "Shannon Gallagher"
date: "September 16, 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}
source("_setup.R")
metrics <- state_metrics |>
  transmute(
    state = state_upper,
    event_records,
    injuries = log1p(injuries),
    deaths = log1p(deaths),
    damage = log1p(property_damage),
    duration = mean_duration
  )
x <- metrics |> select(-state) |> scale()
```

## Reminders, previously, and today...

- explain why distance depends on scale and variables,
- read a multidimensional-scaling display,
- understand PCA as a rotated coordinate system,
- use loadings and explained variance to interpret components.

## Thinking about distance...

Euclidean distance between observations $i$ and $j$ is:

$$
d(i,j)=\sqrt{\sum_{k=1}^{p}(x_{ik}-x_{jk})^2}
$$

The formula is easy. Choosing variables and scales is the hard part.

## Distances in general

- **Euclidean:** straight-line separation; sensitive to scale.
- **Manhattan:** total coordinate-wise separation; less dominated by one large step.
- **Correlation distance:** similarity of profiles after ignoring overall level.
- **Domain-specific distance:** weights or transforms chosen for the actual decision.

There is no universal distance for “similar states.”

## Inspect the distance matrix

```{r distance-matrix}
distance_matrix <- as.matrix(dist(x))
round(distance_matrix[1:6, 1:6], 2)
```

The matrix is symmetric, has zero on the diagonal, and contains one distance for every
pair of observations.

## Dollars dominate without scaling

If we combine raw event counts, deaths, hours, and billions of dollars, damage controls
almost every distance.

Standardization makes one standard deviation count equally in each variable—but that is
also a substantive choice.

## Checkpoint 1 · Define “similar”

::: {.checkpoint}
Two states have similar total damage but very different deaths and event counts. Are they
similar? Choose the variables and scaling that match one plausible client question.
:::

## Multi-dimensional scaling (MDS)

MDS starts from pairwise distances and finds a low-dimensional arrangement that preserves
those distances as closely as possible.

Nearby points are similar under the chosen distance. Axis directions usually have no
standalone meaning.

## MDS example with NOAA state profiles

```{r mds, fig.alt="Multidimensional scaling plot of state-level NOAA profiles."}
mds <- cmdscale(dist(x), k = 2)
mds_df <- tibble(state = metrics$state, dim1 = mds[, 1], dim2 = mds[, 2])

ggplot(mds_df, aes(dim1, dim2, label = state)) +
  geom_hline(yintercept = 0, color = "gray85") +
  geom_vline(xintercept = 0, color = "gray85") +
  geom_text(color = deep_gold, check_overlap = TRUE) +
  labs(title = "Nearby states have similar standardized NOAA profiles", x = "MDS dimension 1", y = "MDS dimension 2")
```

## View structure with additional variables

```{r mds-scaling, fig.alt="Two MDS plots comparing raw and standardized state-level NOAA metrics."}
#| layout-ncol: 2
raw_x <- metrics |> select(-state)
raw_mds <- cmdscale(dist(raw_x), k = 2)

ggplot(tibble(state = metrics$state, x = raw_mds[, 1], y = raw_mds[, 2]), aes(x, y, label = state)) +
  geom_text(color = blue_gray, check_overlap = TRUE) +
  labs(title = "Raw units", x = NULL, y = NULL)

ggplot(mds_df, aes(dim1, dim2, label = state)) +
  geom_text(color = deep_gold, check_overlap = TRUE) +
  labs(title = "Standardized variables", x = NULL, y = NULL)
```

## Evaluate how well two dimensions preserve distance

```{r mds-distance-fit, fig.alt="Scatterplot of original distances and distances in the two-dimensional MDS solution."}
embedded_distance <- as.vector(dist(mds))
original_distance <- as.vector(dist(x))

tibble(original_distance, embedded_distance) |>
  ggplot(aes(original_distance, embedded_distance)) +
  geom_point(alpha = 0.5, color = deep_gold) +
  geom_abline(slope = 1, intercept = 0, color = blue_gray) +
  coord_equal() +
  labs(title = "Departures from the diagonal show distortion", x = "Original standardized distance", y = "Distance in the 2D map")
```

## What MDS does not say

- It does not identify causal groups.
- The axes may rotate or flip without changing the solution.
- A two-dimensional view necessarily loses information.
- Clusters can depend strongly on distance and preprocessing.

## Dimension reduction—searching for variance

Principal component analysis finds directions that:

1. are linear combinations of the original variables,
2. capture as much variance as possible,
3. are mutually orthogonal.

PCA rotates the coordinate system; it does not discover truth.

## Principal Component Analysis (PCA)

For centered variables, the first component has the form

$$
PC_1=a_{11}X_1+a_{21}X_2+\cdots+a_{p1}X_p.
$$

The weights are chosen so that the scores on $PC_1$ have maximum variance, subject to
the loading vector having length one. Later components maximize remaining variance while
remaining orthogonal to earlier components.

## What are principal components?

- **Centering** moves every variable to mean zero.
- **Scaling** gives every variable unit standard deviation.
- PCA on a covariance matrix preserves original-unit variance.
- PCA on a correlation matrix treats standardized variables equally.

With counts, deaths, hours, and dollars together, scaling is essential unless the units
are intentionally weighted.

## Computing Principal Components

```{r pca-summary}
pca <- prcomp(x, center = FALSE, scale. = FALSE)
summary(pca)
```

## Making PCs interpretable with loadings

```{r loadings}
round(pca$rotation[, 1:3], 2)
```

A loading describes how strongly an original variable contributes to a component.
Component signs can flip; relative patterns matter.

## Making PCs interpretable with biplots

```{r pca-biplot, fig.alt="PCA biplot of state scores and variable loading directions."}
loading_scale <- 4
loading_df <- as_tibble(pca$rotation[, 1:2], rownames = "variable") |>
  mutate(PC1 = PC1 * loading_scale, PC2 = PC2 * loading_scale)
scores <- as_tibble(pca$x[, 1:2]) |>
  mutate(state = metrics$state)

ggplot(scores, aes(PC1, PC2)) +
  geom_text(aes(label = state), color = deep_gold, check_overlap = TRUE) +
  geom_segment(data = loading_df, aes(x = 0, y = 0, xend = PC1, yend = PC2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.15, "in")), color = blue_gray) +
  geom_text(data = loading_df, aes(PC1, PC2, label = variable), inherit.aes = FALSE,
            color = blue_gray, vjust = -0.5) +
  labs(title = "Scores locate states; arrows show loading directions")
```

Arrows pointing together indicate variables that contribute similarly to these two
components. Their displayed length depends on the biplot scaling convention.

## Checkpoint 2 · Name PC1

::: {.checkpoint}
Read the PC1 loadings. Propose a plain-language name for the component, then list one
reason your name is incomplete or potentially misleading.
:::

## NOAA state profiles: PC1 and PC2

```{r pca-scores, fig.alt="Scatterplot of state PCA scores on the first two components."}
scores <- as_tibble(pca$x[, 1:2]) |>
  mutate(state = metrics$state)

ggplot(scores, aes(PC1, PC2, label = state)) +
  geom_hline(yintercept = 0, color = "gray85") +
  geom_vline(xintercept = 0, color = "gray85") +
  geom_text(color = deep_gold, check_overlap = TRUE) +
  labs(title = "PCA summarizes state profiles on new axes")
```

## How many principal components should we use?

```{r scree, fig.alt="Scree plot showing variance explained by each principal component."}
variance <- pca$sdev^2 / sum(pca$sdev^2)
tibble(component = seq_along(variance), variance) |>
  ggplot(aes(component, variance)) +
  geom_col(fill = gold) +
  geom_line() +
  geom_point() +
  scale_y_continuous(labels = percent) +
  labs(title = "A scree plot makes information loss visible", x = "Principal component", y = "Variance explained")
```

## Create a scree plot (aka "elbow plot") to choose

```{r cumulative-variance, fig.alt="Cumulative variance explained by the principal components."}
tibble(component = seq_along(variance), cumulative = cumsum(variance)) |>
  ggplot(aes(component, cumulative)) +
  geom_line(color = deep_gold, linewidth = 1) +
  geom_point(color = deep_gold, size = 3) +
  geom_hline(yintercept = 0.8, linetype = 2, color = blue_gray) +
  scale_y_continuous(labels = percent, limits = c(0, 1)) +
  labs(title = "Retaining components trades simplicity for information", x = "Number of components", y = "Cumulative variance explained")
```

## PCA: *singular value decomposition (SVD)*

For a centered data matrix $X$,

$$
X=UDV^T.
$$

- columns of $V$ are the loading directions;
- $UD$ contains the component scores;
- squared singular values determine variance explained.

This connection makes PCA computationally stable and links it to many other matrix
methods.

## Checkpoint 3 · Keep how many?

::: {.checkpoint}
Using the scree plot, choose how many components you would retain for exploration. Report
the cumulative variance and one important reason variance explained is not the only
criterion.
:::

## Recap and next steps

- Distance encodes a definition of similarity.
- MDS preserves pairwise distances; its axes are usually not interpreted directly.
- PCA preserves variance along linear directions.
- Loadings, scores, and information loss must be read together.

Next: nonlinear embeddings and the transition to time.