Nonlinear regression and more quantitative data visualization


36-315: Statistical Graphics and Visualization, Summer 2026

Beyond linear regression

  • The truth is (almost) never linear

  • But often the linearity assumption is good enough

  • Many nonparametric/nonlinear regression methods

    • splines

    • local regression

    • generalized additive models

    • kernel regression

    • wavelets

Nonlinear regression

  • Assume normality, but not linearity where \(f(x)\) is some unknown function \[Y_i \stackrel{\textsf{iid}}{\sim} N(f(X_i), \sigma^2)\]

  • Do not make explicit assumptions about the functional form of \(f\) (i.e., \(f\) is estimated from the data)

  • Intuition: Any nonlinear function is locally linear

  • Local linear regressions fits a bunch of… local linear regressions, and then glues them together

  • Local linear regression is basically weighted linear regression, where only “local units” get weight

Weighted linear regression

  • In linear regression, we solve

\[\arg \min_{\beta_0,\beta_1} \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i)^2\]

  • In weighted linear regression, we solve

\[\arg \min_{\beta_0,\beta_1} \sum_{i=1}^n w_i \cdot (Y_i - \beta_0 - \beta_1 X_i)^2\]

  • Local linear regression is exactly the same, except the weights depend on which \(x\) we want to estimate \(f(x)\)

    • Fit a different linear regression everywhere, weighting the data points by how close they are to the point of interest

Local linear regression algorithm

Local regression at \(X=x_0\)

  1. Compute the fraction \(k/n\) of observations who \(x_i\) are the closest to \(x_0\)

  2. Assign a weight \(w_i = w(x_i, x_0)\) to each point in this local neighborhood, so that the closest/furthest has the highest/lowest weights. All but these \(k\) nearest neighbors get weight zero

  3. Fit a weighted linear (least squares) regression

\[\text{arg }\underset{\beta_0, \beta_1}{\text{min}} \sum_i^n w_i(x_0) \cdot \big(Y_i - \beta_0 - \beta_1 X_i \big)^2\]

Local linear regression via LOESS

In local linear regression, we estimate \(f(X_i)\) and solve

\[\text{arg }\underset{\beta_0, \beta_1}{\text{min}} \sum_i^n w_i(x) \cdot \big(Y_i - \beta_0 - \beta_1 X_i \big)^2\]

  • Notice the weights depend on \(x\), so observations close to \(x\) gets more weight in estimating \(f(x)\)

geom_smooth() uses tri-cubic weighting:

\[w_i(d_i) = \begin{cases} (1 - |d_i|^3)^3, \text{ if } i \in \text{neighborhood of } x\\ 0 \text{ if } i \notin \text{neighborhood of } x \end{cases}\]

  • \(d_i\): distance between \(x\) and \(X_i\) scaled to be between 0 and 1

  • span: controls fraction of observations in neighborhood around the observation of interest

    • Default: span = 0.75 (use 75% of the observations to fit local linear regression with weights)

Local linear regression, animated

Displaying trend lines: LOESS

library(tidyverse)
theme_set(theme_light())
penguins |> 
  ggplot(aes(x = flipper_len, y = body_mass)) +
  geom_point(alpha = 0.5) + 
  geom_smooth()

  • LOESS is default behavior when \(n \leq 1000\)

  • Default span = 0.75

  • For \(n > 1000\), mgcv::gam() is used with
    formula = y ~ s(x, bs = "cs") and
    method = "REML"

Adjusting span

Pairs plots

  • Create a pairs plot to see all pairwise relationships in one plot

  • Pairs plot can include the various kinds of pairwise plots

    • Two quantitative variables: scatterplot

    • One categorical, one quantitative: side-by-side violins, stacked histograms, overlaid densities

    • Two categorical: stacked bars, side-by-side bars, mosaic plots

Pairs plots

library(GGally)
penguins |> 
  select(bill_len, bill_dep, flipper_len, body_mass) |> 
  ggpairs()

Pair plots of quantitative variables

  • Diagonal: marginal distributions

  • Off-diagonal: joint (pairwise) distributions or statistical summaries (e.g., correlation)

  • Matrix of plots is symmetric

penguins |> 
  select(bill_len, bill_dep, flipper_len, body_mass) |> 
  ggpairs(aes(alpha = 0.2))

Syntax note: change alpha within aes (annoying…)

Pair plots with categorical variables

  • Can either use select() before plotting or specify columns directly
penguins |> 
  ggpairs(columns = c("bill_len", "body_mass", "island"),
          aes(alpha = 0.2, color = species))

More customization

  • Can customize which figures are displayed in which panels

  • For more details, check out the package documentation here

penguins |> 
  ggpairs(
    columns = c("bill_len", "body_mass", "island"),
    mapping = aes(alpha = 0.2, color = species), 
    lower = list(continuous = "smooth_lm", 
                 combo = "facetdensitystrip"),
    upper = list(continuous = "cor", 
                 combo = "facethist")
  )

More customization

  • Can also incorporate ggplot2 customization
penguins |>
  ggpairs(columns = c("bill_len", "bill_dep", "flipper_len"),
          aes(color = species, alpha = 0.2)) +
  scale_color_manual(values = c("orange", "blue", "cyan4")) +
  scale_fill_manual(values = c("orange", "blue", "cyan4")) +
  theme_minimal() +
  theme(strip.text = element_text(face = "bold"))

Correlation heatmaps

  • Use correlate() and scretch() functions from the corrr package to compute correlation matrix and pivot

  • Use geom_tile() to create a heatmap (tile map)

library(corrr)
penguins |> 
  select(bill_len, bill_dep, flipper_len, body_mass) |> 
  correlate(diagonal = 1) |> 
  stretch() |> 
  ggplot(aes(x, y, fill = r)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(high = "red", mid = "white", low = "blue", 
                       midpoint = 0, limits = c(-1, 1)) +
  coord_fixed()


Parallel coordinates plots

  • Display each variable side-by-side on standardized axis

  • Connect observations with lines

  • Useful when working with a moderate number of observations and variables
    (can be overwhelming with too many)

  • Adjust the transparency of the lines to handle overlap

library(GGally)
penguins |>
  select(bill_len, bill_dep, flipper_len, body_mass) |> 
  ggparcoord(alphaLines = 0.2)

Parallel coordinates plots

  • Color lines by a categorical variable

  • Can be useful for revealing group structure

  • How do we order the x-axis?

penguins |>
  ggparcoord(columns = 3:6, alphaLines = 0.2, 
             groupColumn = "species", order = c(6, 5, 3, 4)) +
  theme(legend.position = "top")

Parallel coordinates plots

  • By default, the y-axis is scaled by subtracting the mean and dividing by the standard deviation

  • Use scale = "uniminmax" so that minimum of the variable is 0 and the maximum is 1

penguins |>
  ggparcoord(columns = 3:6, alphaLines = 0.2, 
             groupColumn = "species", scale = "uniminmax") +
  theme(legend.position = "top")

Distance between observations

  • Goal: visualize structure among observations using distances matrices
  • What does it means for two or more observations to be similar or different?

  • This requires characterizing the distance between observations

    • Clusters: groups of observations that are “close” together
  • This is easy to do for 2 quantitative variables: just make a scatterplot (possibly with contours or heatmap)

  • Define distance beyond 2D:

    • Let \(\boldsymbol{x}_i = (x_{i1}, \dots, x_{ip})\) be a vector of \(p\) features for observation \(i\)

    • Question of interest: How “far away” is \(\boldsymbol{x}_i\) from \(\boldsymbol{x}_j\)?

    • When looking at a scatterplot, Euclidean distance is used \[d(\boldsymbol{x}_i, \boldsymbol{x}_j) = \sqrt{(x_{i1} - x_{j1})^2 + \dots + (x_{ip} - x_{jp})^2}\]

Distances in general

  • Let \(d(\boldsymbol{x}_i, \boldsymbol{x}_j)\) denote the pairwise distance between two observations \(i\) and \(j\)

    1. Identity: \(\boldsymbol{x}_i = \boldsymbol{x}_j \Leftrightarrow d(\boldsymbol{x}_i, \boldsymbol{x}_j) = 0\)

    2. Non-negativity: \(d(\boldsymbol{x}_i, \boldsymbol{x}_j) \geq 0\)

    3. Symmetry: \(d(\boldsymbol{x}_i, \boldsymbol{x}_j) = d(\boldsymbol{x}_j, \boldsymbol{x}_i)\)

    4. Triangle inequality: \(d(\boldsymbol{x}_i, \boldsymbol{x}_j) \leq d(\boldsymbol{x}_i, \boldsymbol{x}_k) + d(\boldsymbol{x}_k, \boldsymbol{x}_j)\)

  • Distance matrix: matrix \(D\) of all pairwise distances

    • \(D_{ij} = d(\boldsymbol{x}_i, \boldsymbol{x}_j)\)

    • Note: \(D_{ii} = 0\) and \(D_{ij} = D_{ji}\)

\[ D = \begin{pmatrix} 0 & D_{12} & \cdots & D_{1n} \\ D_{21} & 0 & \cdots & D_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ D_{n1} & \cdots & \cdots & 0 \end{pmatrix} \]

Multidimensional scaling (MDS)

  • Find a \(k\)-dimensional projection of the original \(p\)-dimensional dataset (where \(k < p\))

  • Preserve the distances as much as possible

    • i.e., pairs of observations with small distances in the original \(p\)-column dataset still have small distances in the smaller \(k\)-column dataset (approximately equal to the actual distance)

  • Formally, MDS attempts to find \(\mathbf y_1, \ldots, \mathbf y_n \in \mathbb{R}^k\) whose distance matrix is approximately \(\mathbf D=(d_{ij}: \, i,j=1, \ldots , n)\) (the original distance matrix)

    • i.e., for which \(\textsf{distance}(\mathbf y_i, \mathbf y_j) \approx d_{ij}\)

MDS example: Starbucks drinks

starbucks <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2021/2021-12-21/starbucks.csv") |>
  # convert columns to numeric that were saved as character
  mutate(trans_fat_g = as.numeric(trans_fat_g), fiber_g = as.numeric(fiber_g))
names(starbucks)
 [1] "product_name"    "size"            "milk"            "whip"           
 [5] "serv_size_m_l"   "calories"        "total_fat_g"     "saturated_fat_g"
 [9] "trans_fat_g"     "cholesterol_mg"  "sodium_mg"       "total_carbs_g"  
[13] "fiber_g"         "sugar_g"         "caffeine_mg"
  • First, scale the data to ensure the variables are on the same scale (i.e., variances are equal to 1)
starbucks_scaled <- starbucks |> 
  select(serv_size_m_l:caffeine_mg) |> 
  mutate(across(everything(), ~ .x / sd(.x, na.rm = TRUE)))
  • Then, compute distance matrix with dist() (default: method = "euclidean")
starbucks_dist <- dist(starbucks_scaled)
class(starbucks_dist)
[1] "dist"

Implementing MDS

  • Use cmdscale() function and supply a distance matrix d and the number of dimensions k
    (usually k = 2 for visualization purposes)

  • What are the dimensions of the input and output?

  • How do we interpret the 2D projection?

starbucks_mds <- cmdscale(d = starbucks_dist, k = 2)

starbucks <- starbucks |> 
  mutate(mds1 = starbucks_mds[,1],
         mds2 = starbucks_mds[,2])

starbucks |>
  ggplot(aes(x = mds1, y = mds2)) +
  geom_point(alpha = 0.5) +
  labs(x = "MDS coordinate 1", y = "MDS coordinate 2")

Visualizing MDS by categorical/continuous variables

Code 
starbucks |>
  ggplot(aes(x = mds1, y = mds2, color = size)) +
  geom_point(alpha = 0.5) +
  labs(x = "MDS coordinate 1", y = "MDS coordinate 2")

Code 
starbucks |>
  ggplot(aes(x = mds1, y = mds2, color = sugar_g)) +
  geom_point(alpha = 0.5) +
  scale_color_gradient(low = "darkblue", high = "darkorange") +
  labs(x = "MDS coordinate 1", y = "MDS coordinate 2")

Another MDS example: MCU movies

Which of these movies are most similar/different?

Code 
mcu <- read_csv("https://raw.githubusercontent.com/qntkhvn/36-315-summer26/refs/heads/master/data/mcu.csv")
# names(mcu)

mcu_scaled <- mcu |> 
  select(-c(film, category, year)) |> 
  mutate(across(everything(), ~ .x / sd(.x, na.rm = TRUE)))

mcu_mds <- cmdscale(d = dist(mcu_scaled), k = 2)

mcu <- mcu |>
  mutate(mds1 = mcu_mds[, 1], mds2 = mcu_mds[, 2])

library(ggrepel)
mcu |>
  ggplot(aes(x = mds1, y = mds2)) +
  geom_point(color = "darkblue", size = 2, alpha = 0.5) +
  geom_text_repel(aes(label = film), alpha = 0.75, color = "black") +
  labs(x = "MDS coordinate 1", y = "MDS coordinate 2")