| SOCR ≫ | DSPA ≫ | DSPA3 Topics ≫ |
DSPA3: Data Science and Predictive Analytics (UMich HS650)
Chapter 2: Data Quality, Missingness, and Exploratory Visual
Analytics
Chapter 2: Data Quality, Missingness, and Exploratory Visual Analytics
SOCR/MCAIM (Ivo Dinov)
SOCR/MCAIM (Ivo Dinov)
August 25, 2026
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork) # panel composition, replaces gridExtra / par(mfrow)
library(plotly) # interactive figures and ALL 3-D graphics
library(DT)How this chapter uses graphics
Every two-dimensional figure is drawn with
ggplot2and rendered statically. Immediately after each one you will find the equivalentplot_ly()code in a chunk markedeval=FALSE, echo=TRUE, the code is printed in these notes but not executed, so you can copy it into a live session for an interactive demonstration without inflating the size of this document.Every three-dimensional figure, surfaces, 3-D scatter, image stacks, volumetric neuroimaging, interactive dendrograms and scatterplot matrices — is drawn with
plot_ly()and evaluated, because rotation and zoom carry information no static image can. Those figures are live in the rendered HTML.The rationale is perceptual, not aesthetic: static images are reproducible, printable, and diffable; interaction earns its cost only when the reader must navigate rather than read (§2.10).
1 Learning objectives
After completing this chapter you will be able to:
- Ingest tabular, compressed, foreign-format (SPSS), spreadsheet, and HTML-table data, and verify what you loaded via checksum and structure.
- Choose an appropriate measure of centrality and dispersion by reasoning about the breakdown point and the AM–GM–HM inequality, rather than by habit.
- Derive the bias–variance tradeoff that determines optimal histogram bin width and kernel bandwidth, and apply Sturges, Scott, Freedman–Diaconis, and Silverman’s rules.
- Analyze a contingency table correctly: test statistic, cell contributions, standardized residuals, effect size, and validity conditions.
- Classify missingness as MCAR / MAR / MNAR using formal notation, and prove why single mean imputation deflates variance.
- Implement EM imputation, state its convergence property, and analyze its computational complexity.
- Pool multiple imputations using Rubin’s rules and interpret the between- vs. within-imputation variance decomposition.
- Rebalance an imbalanced cohort without leaking information from the test set, and quantify the optimism that leakage induces.
- Select a visualization for a given analytic question using the Cleveland–McGill accuracy ordering.
- Fit univariate and mixture distribution models by maximum likelihood and EM, and assess fit with methods valid under estimated parameters.
Estimated time: 8–10 hours including exercises.
Prerequisites: Chapter 1, in particular the
reproducibility conventions (§1.6.6) and the dspa_read()
loader (§1.6.7), both of which this chapter reuses.
2 PART I, DATA HANDLING, QUALITY, AND MISSINGNESS
3 Ingesting data
Every analysis begins by moving bytes into memory and losing information in the process. The question is only whether you know which information you lost.
3.1 Native R serialization
#> [1] "data.frame"# Chapter 1 convention: never write into the working directory.
rds_path <- file.path(tempdir(), "iris.rds")
saveRDS(iris, rds_path)
iris2 <- readRDS(rds_path)
identical(iris, iris2)#> [1] TRUE
The iris
data records four morphological measurements, sepal and petal length
and width, in centimetres, on 50 flowers from each of three species
(setosa, versicolor, virginica). Fisher used
them in his 1936
discriminant-analysis paper; they remain the standard small
multivariate test case.
3.2 Delimited text
The World Drinking Water study (Case Study 07) reports, by country and year, the share of the population using improved drinking-water sources and improved sanitation facilities.
water <- dspa_read(
url = "https://umich.instructure.com/files/399172/download?download_frd=1",
name = "CaseStudy07_WorldDrinkingWater_Data.csv",
header = TRUE,
encoding = "UTF-8", # fileEncoding = "UTF-8",
refresh = TRUE # <--- THIS FORCES RE-DOWNLOADING
)
names(water) <- c("year", "region", "country", "residence_area",
"improved_water", "sanitation_facilities")
dim(water)#> [1] 3331 6#> MD5: 014a58df26d129872337a1eebe3858efNotes.
read.csv()’s default separator is
sep = ",", a comma, with no space. Passing
sep = ", " makes the space part of the delimiter, so every
field after the first acquires a leading space and every numeric column
returns as character. The same mistake with
write.csv()/write.table() may produce files
that will not round-trip.
Encoding is not optional metadata. This file
CaseStudy07_WorldDrinkingWater_Data.csv intentionally
contains non-ASCII country names; declaring fileEncoding
explicitly is the difference between Côte d'Ivoire and
C te d'Ivoire. When the encoding, is unknowne, using
dreadr::guess_encoding() reports it.
out <- file.path(tempdir(), "water_subset.csv")
write.csv(head(water, 20), out, row.names = FALSE) # sep="," and no row names
readLines(out, n = 2)#> [1] "\"year\",\"region\",\"country\",\"residence_area\",\"improved_water\",\"sanitation_facilities\""
#> [2] "1990,\"Africa\",\"Algeria\",\"Rural\",88,77"3.3 Compressed and foreign formats
Health data frequently arrives as an SPSS .sav inside a
.zip. The National
Ambulatory Medical Care Survey is a good stress test: 28,332
encounters × 1,096 variables.
library(foreign)
zip_path <- file.path(tempdir(), "namcs2015.zip")
download.file("https://umich.instructure.com/files/8111611/download?download_frd=1",
zip_path, mode = "wb")
sav <- unzip(zip_path, files = "namcs2015-spss.sav",
exdir = tempdir(), overwrite = TRUE)
namcs <- read.spss(sav, to.data.frame = TRUE)
dim(namcs) # 28332 x 1096foreign::read.spss() is the base option;
haven::read_sav() is generally preferable now, it preserves
SPSS value labels as labelled vectors rather than silently
converting them to factors, which matters when a code of 9
means “refused” rather than the number nine.
Other formats: readxl::read_excel() /
openxlsx::read.xlsx() for spreadsheets,
arrow::read_parquet() for columnar data at scale,
jsonlite::fromJSON() for APIs, DBI + a backend
driver for databases (Appendix A), and
rvest::html_table() for HTML (§2.17).
3.4 Look before you compute
#> 'data.frame': 3331 obs. of 6 variables:
#> $ year : int 1990 1990 1990 1990 1990 1990 1990 1990 1990 1990 ...
#> $ region : chr "Africa" "Africa" "Africa" "Africa" ...
#> $ country : chr "Algeria" "Angola" "Benin" "Botswana" ...
#> $ residence_area : chr "Rural" "Rural" "Rural" "Rural" ...
#> $ improved_water : num 88 42 49 86 39 67 34 46 37 83 ...
#> $ sanitation_facilities: num 77 7 0 22 2 42 27 12 4 11 ...#> water[, c("year", "improved_water", "sanitation_facilities")]
#>
#> 3 Variables 3331 Observations
#> --------------------------------------------------------------------------------
#> year
#> n missing distinct Info Mean Gmd
#> 3331 0 6 0.972 2002 8.788
#>
#> Value 1990 1995 2000 2005 2010 2012
#> Frequency 520 561 570 570 556 554
#> Proportion 0.156 0.168 0.171 0.171 0.167 0.166
#>
#> For the frequency table, variable is rounded to the nearest 0
#> --------------------------------------------------------------------------------
#> improved_water
#> n missing distinct Info Mean Gmd .05 .10
#> 3299 32 124 0.988 84.9 19.29 42 54
#> .25 .50 .75 .90 .95
#> 77 93 99 100 100
#>
#> lowest : 3 5 9 11 13
#> highest: 99.8836 99.8903 99.897 99.8983 100
#> --------------------------------------------------------------------------------
#> sanitation_facilities
#> n missing distinct Info Mean Gmd .05 .10
#> 3196 135 106 0.993 68.87 34.55 10 19
#> .25 .50 .75 .90 .95
#> 42 81 97 100 100
#>
#> lowest : 0 1 2 3 4, highest: 96 97 98 99 100
#> --------------------------------------------------------------------------------str() gives types and shape;
Hmisc::describe() adds distinct-value counts, missingness,
quantiles, and the highest and lowest observed values, the last of which
is the fastest way to spot a sentinel code (-99,
999) masquerading as data.
4 Centrality: which average, and why
4.1 Three means and an inequality
For a sample \(x_1,\dots,x_n\) of strictly positive values, three classical means are defined:
\[ \text{AM} = \underbrace{\frac{1}{n}\sum_{i=1}^{n} x_i}_{\text{Arithmetic Mean}}, \qquad \text{GM} = \underbrace{\left(\prod_{i=1}^{n} x_i\right)^{1/n}}_{\text{Geometric Mean}}, \qquad \text{HM} = \underbrace{\frac{n}{\sum_{i=1}^{n} 1/x_i}}_{\text{Harmonic Mean}}. \]
They are ordered, always:
\[\text{HM} \;\le\; \text{GM} \;\le\; \text{AM},\]
with equality iff all \(x_i\) are equal. (The AM–GM step follows from Jensen’s inequality applied to the concave function \(\log\); the GM–HM step is AM–GM applied to \(1/x_i\).)
Note the domain: positive, not “non-negative”. A single zero drives the geometric mean to zero regardless of every other value, and \(\log 0\) is undefined, so the usual computational route \(\mathrm{GM}=\exp\!\big(\tfrac1n\sum\log x_i\big)\) fails outright.
x <- c(40, 56, 99)
am <- mean(x)
gm <- exp(mean(log(x))) # numerically safer than prod(x)^(1/n)
hm <- length(x) / sum(1 / x)
c(HM = hm, GM = gm, AM = am)#> HM GM AM
#> 56.64850 60.52866 65.00000#> [1] TRUEWhich one is right depends on what “combining” means for your quantity.
| Use | When the quantity is | Example |
|---|---|---|
| Arithmetic | Additive | Total dose across visits |
| Geometric | Multiplicative / ratio-scaled | Antibody titers, fold-changes, growth rates |
| Harmonic | A rate whose denominator is fixed | Average speed over equal distances; \(F_1\) score |
Averaging log-normally distributed lab values arithmetically overstates the typical value. The geometric mean is the natural center because the logs are symmetric.
4.2 Robustness and the breakdown point
The breakdown point of an estimator is the smallest fraction of observations that, if replaced by arbitrary values, can drive the estimate to infinity. It formalizes “how much contamination can this statistic survive.”
| Estimator | Breakdown point |
|---|---|
| Mean | \(1/n \to 0\) |
| \(\alpha\)-trimmed mean | \(\alpha\) |
| Median | \(0.5\) (the maximum attainable) |
| Standard deviation | \(1/n \to 0\) |
| IQR | \(0.25\) |
| MAD | \(0.5\) |
One corrupted value is enough to move the mean anywhere. The median needs half the sample.
set.seed(11)
clean <- rnorm(100, mean = 50, sd = 5)
dirty <- clean; dirty[1] <- 1e6 # one data-entry error
rbind(
clean = c(mean = mean(clean), trim10 = mean(clean, trim = 0.10),
median = median(clean), sd = sd(clean), IQR = IQR(clean),
MAD = mad(clean)),
dirty = c(mean = mean(dirty), trim10 = mean(dirty, trim = 0.10),
median = median(dirty), sd = sd(dirty), IQR = IQR(dirty),
MAD = mad(dirty))
) |> round(2)#> mean trim10 median sd IQR MAD
#> clean 49.38 49.22 49.07 4.57 6.61 4.92
#> dirty 10048.91 49.34 49.14 99995.06 6.65 5.00The mean and standard deviation are destroyed; the trimmed mean,
median, IQR, and MAD are essentially unchanged. This is why
mean(x, trim = 0.08) appears in exploratory code, it is not
a stylistic choice but an \(8\%\)
breakdown guarantee.
4.3 The mode
The mode is the most frequent value. It is the only centrality measure defined for nominal data, and the only one that can be multimodal, itself an important signal, usually that the sample mixes two populations (§2.14).
sapply(water[c("year", "region", "residence_area")],
\(v) { tb <- table(v); paste(names(tb)[tb == max(tb)], collapse = ", ") })#> year region residence_area
#> "2000, 2005" "Europe" "Urban"region and residence_area are unimodal;
year is bimodal (2000 and 2005 tie), which tells you
something about survey administration rather than about water quality.
See the SOCR
EDA centrality module for further discussion.
5 Dispersion and distributional shape
5.1 The five-number summary
Minimum, \(Q_1\), median, \(Q_3\), maximum. The interquartile range \(\mathrm{IQR}=Q_3-Q_1\) spans the middle half of the data and has a breakdown point of \(0.25\).
iw <- water$improved_water
c(min = min(iw, na.rm = TRUE),
Q1 = unname(quantile(iw, 0.25, na.rm = TRUE)),
med = median(iw, na.rm = TRUE),
Q3 = unname(quantile(iw, 0.75, na.rm = TRUE)),
max = max(iw, na.rm = TRUE),
IQR = IQR(iw, na.rm = TRUE),
range_width = diff(range(iw, na.rm = TRUE))) |> round(2)#> min Q1 med Q3 max IQR
#> 3 77 93 99 100 22
#> range_width
#> 97#> 0% 20% 40% 60% 80% 100%
#> 3 71 89 97 100 100Read the asymmetry: the gap from the minimum to \(Q_1\) is enormous while the gap from \(Q_3\) to the maximum is about one percentage point. The lower quarter is stretched, the upper quarter compressed against a ceiling of 100%. That is a left-skewed, ceiling-bounded variable, and it is why \(\text{mean} < \text{median}\) here. Percentages bounded at 100 behave this way routinely — which is an argument for a logit or arcsine transform before any normality-assuming procedure.
\[ \text{Var}(X)=\sigma^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar x)^2, \qquad \text{SD}(X)=\sigma=\sqrt{\text{Var}(X)}. \]
5.2 Histograms: how wide should the bins be?
A histogram with bin width \(h\) estimates the density as
\[\hat f_h(x)=\frac{\#\{i:\; x_i \in B_j\}}{n\,h},\qquad x\in B_j .\]
This is a genuine estimator with a genuine bias–variance tradeoff:
- Bias. Within a bin the estimate is constant, so it cannot track the slope of \(f\). A Taylor expansion gives \(\text{Bias}\,\hat f_h(x)=O(h)\), wide bins oversmooth.
- Variance. The count in a bin is roughly \(\text{Binomial}(n, hf(x))\), so \(\text{Var}\,\hat f_h(x)= O\!\big(\tfrac{1}{nh}\big)\), narrow bins are noisy.
Integrating the squared bias plus variance gives the asymptotic mean integrated squared error
\[\mathrm{AMISE}(h)=\underbrace{\frac{h^{2}}{12}R(f')}_{\text{bias}^2}+\underbrace{\frac{1}{nh}}_{\text{variance}},\qquad R(g)=\int g(u)^2\,du .\]
Minimizing over \(h\):
\[\boxed{\;h^{\star}=\left(\frac{6}{n\,R(f')}\right)^{1/3}\propto n^{-1/3}\;}\]
so the optimal bin width shrinks like \(n^{-1/3}\) and the number of bins grows like \(n^{1/3}\). The three rules you will meet in practice are all attempts to estimate \(R(f')\):
| Rule | Formula | Basis |
|---|---|---|
| Sturges | \(k=\lceil \log_2 n\rceil+1\) bins | Binomial approximation; assumes normality, undersmooths large \(n\) |
| Scott | \(h=3.49\,\hat\sigma\,n^{-1/3}\) | AMISE-optimal under normality |
| Freedman–Diaconis | \(h=2\,\mathrm{IQR}\,n^{-1/3}\) | Same rate, robust scale, preferred for skewed or heavy-tailed data |
Note the \(\sqrt{n}\) rule sometimes
quoted is not on this list. It grows the bin count like \(n^{1/2}\) rather than \(n^{1/3}\), so it oversmooths small samples
and undersmooths large ones. R’s hist() defaults to
Sturges; nclass.scott() and nclass.FD() give
the others.
iw_c <- iw[!is.na(iw)]
n <- length(iw_c)
data.frame(
rule = c("sqrt (avoid)", "Sturges", "Scott", "Freedman-Diaconis"),
bins = c(ceiling(sqrt(n)), nclass.Sturges(iw_c),
nclass.scott(iw_c), nclass.FD(iw_c))
)mk <- function(bins, lab)
ggplot(data.frame(x = iw_c), aes(x)) +
geom_histogram(bins = bins, fill = "steelblue", colour = "white") +
labs(title = sprintf("%s (%d bins)", lab, bins), x = NULL, y = "count") +
theme_dspa(10)
(mk(ceiling(sqrt(n)), "sqrt") | mk(nclass.Sturges(iw_c), "Sturges")) /
(mk(nclass.scott(iw_c), "Scott") | mk(nclass.FD(iw_c), "Freedman-Diaconis")) +
plot_annotation(title = "Bin-width rules on the same data (improved water, %)",
theme = theme_dspa())
# --- Interactive equivalent (run in a live session) ------------------------
plot_ly(x = ~iw_c, type = "histogram", nbinsx = nclass.FD(iw_c),
name = "Freedman-Diaconis") |>
add_histogram(x = ~iw_c, nbinsx = nclass.Sturges(iw_c), name = "Sturges",
opacity = 0.6) |>
layout(barmode = "overlay", bargap = 0.05,
title = "Improved water (%), bin-rule comparison",
xaxis = list(title = "Percent of population"),
yaxis = list(title = "Frequency"),
legend = list(orientation = "h"))5.3 Kernel density estimation: how wide should the kernel be?
A histogram is discontinuous and depends on bin origin as well as width. The kernel density estimator removes both problems:
\[\hat f_h(x)=\frac{1}{nh}\sum_{i=1}^{n}K\!\left(\frac{x-x_i}{h}\right), \qquad \int K(u)\,du=1,\;\int uK(u)\,du=0,\;\mu_2(K)=\int u^2K(u)\,du<\infty .\]
Because \(K\) is smooth, the bias expansion reaches second order:
\[\mathrm{AMISE}(h)=\frac{h^{4}}{4}\mu_2(K)^{2}R(f'')+\frac{R(K)}{nh} \;\;\Longrightarrow\;\; \boxed{\;h^{\star}=\left(\frac{R(K)}{\mu_2(K)^2R(f'')\,n}\right)^{1/5}\propto n^{-1/5}\;}\]
The optimal bandwidth shrinks more slowly than the optimal histogram bin width (\(n^{-1/5}\) vs \(n^{-1/3}\)), and the attainable error rate is better — \(O(n^{-4/5})\) versus \(O(n^{-2/3})\). That is the precise sense in which kernel smoothing beats binning.
Plugging a Gaussian reference for \(f\) gives Silverman’s rule of thumb:
\[h_{\text{Silverman}} = 0.9\,\min\!\left(\hat\sigma,\;\frac{\mathrm{IQR}}{1.34}\right)n^{-1/5},\]
which is bw.nrd0(), R’s default. The min
with a scaled IQR is a robustness guard against the outliers that would
otherwise inflate \(\hat\sigma\) and
oversmooth.
The kernel shape barely matters. Relative to the Epanechnikov kernel (which is AMISE-optimal), the Gaussian kernel has efficiency \(\approx 0.951\) and the uniform \(\approx 0.930\). The bandwidth matters enormously; the kernel choice is nearly irrelevant. Beginners usually get this backwards.
bw_sil <- bw.nrd0(iw_c)
ggplot(data.frame(x = iw_c), aes(x)) +
geom_histogram(aes(y = after_stat(density)), bins = nclass.FD(iw_c),
fill = "grey88", colour = "white") +
geom_density(aes(colour = "h/4 (undersmoothed)"), bw = bw_sil / 4, linewidth = 0.7) +
geom_density(aes(colour = "h (Silverman)"), bw = bw_sil, linewidth = 1.1) +
geom_density(aes(colour = "4h (oversmoothed)"), bw = bw_sil * 4, linewidth = 0.7) +
scale_colour_manual(values = c("h/4 (undersmoothed)" = "firebrick",
"h (Silverman)" = "black",
"4h (oversmoothed)" = "steelblue")) +
labs(title = "Kernel density estimate: bandwidth is the whole game",
subtitle = sprintf("Silverman bandwidth h = %.2f", bw_sil),
x = "Improved water (%)", y = "Density", colour = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
d_lo <- density(iw_c, bw = bw_sil / 4)
d_md <- density(iw_c, bw = bw_sil)
d_hi <- density(iw_c, bw = bw_sil * 4)
plot_ly() |>
add_histogram(x = ~iw_c, histnorm = "probability density",
name = "Histogram", opacity = 0.35,
marker = list(color = "grey")) |>
add_lines(x = d_lo$x, y = d_lo$y, name = "h/4 (undersmoothed)") |>
add_lines(x = d_md$x, y = d_md$y, name = "h (Silverman)",
line = list(width = 4)) |>
add_lines(x = d_hi$x, y = d_hi$y, name = "4h (oversmoothed)") |>
layout(bargap = 0.05,
title = "Density of (%) improved water quality",
xaxis = list(title = "Percent"), yaxis = list(title = "Density"),
legend = list(orientation = "h"))5.4 Boxplots and what “outlier” means
The box spans \(Q_1\) to \(Q_3\) with the median inside. The whiskers
extend to the most extreme observation within \(1.5\times\mathrm{IQR}\) of the nearer
hinge; anything beyond is drawn as a point. There is no \(3\times\) anywhere in what R draws, a \(3\times\mathrm{IQR}\) fence is a separate
convention for extreme outliers, and boxplot()’s
range = argument would have to be changed to produce
it.
The fence is worth calibrating. Under exact normality, \(\mathrm{IQR}=2\times 0.6745\,\sigma \approx 1.349\sigma\), so the fences sit at
\[Q_1 - 1.5\,\mathrm{IQR} \approx \mu - 2.698\sigma, \qquad Q_3 + 1.5\,\mathrm{IQR} \approx \mu + 2.698\sigma,\]
giving a per-observation false-flag rate of
\[2\,\Phi(-2.698)\approx 0.70\% .\]
#> [1] 0.006976603In a clean normal sample of \(n=1{,}000\) you should therefore expect about seven “outliers”. A flagged point is a hypothesis, not a verdict, never delete on the strength of a whisker alone.
water_long <- water |>
dplyr::select(improved_water, sanitation_facilities) |>
pivot_longer(everything(), names_to = "indicator", values_to = "percent") |>
filter(!is.na(percent))
ggplot(water_long, aes(indicator, percent, fill = indicator)) +
geom_boxplot(outlier.alpha = 0.35, width = 0.55) +
scale_fill_manual(values = c("improved_water" = "#4C9BE8",
"sanitation_facilities" = "#E8A24C")) +
labs(title = "Improved water quality and sanitation facilities",
subtitle = "Whiskers at 1.5 x IQR; points beyond are flagged, not condemned",
x = NULL, y = "Percent of population") +
theme_dspa() + theme(legend.position = "none")
# --- Interactive equivalent ------------------------------------------------
plot_ly(y = ~water$improved_water, type = "box",
name = "improved water quality") |>
add_trace(y = ~water$sanitation_facilities, name = "sanitation") |>
layout(title = "Boxplots of Improved Water Quality and Sanitation Facilities",
yaxis = list(title = "Percent"))5.5 Skewness, and reading it off the plot
For a left-skewed variable the long tail is on the left, and typically \(\text{mean} < \text{median}\). Right skew reverses both.
set.seed(2026)
N <- 10000
sk <- bind_rows(
data.frame(x = rnbinom(N, size = 5, prob = 0.1), shape = "Right skewed NB(5, 0.1)"),
data.frame(x = rnorm(N, mean = 15, sd = 3.7), shape = "Symmetric N(15, 3.7)"),
data.frame(x = 100 - rgamma(N, shape = 2, rate = 0.4), shape = "Left skewed 100 - Gamma(2, 0.4)")
)
sk_stats <- sk |> summarise(mean = mean(x), median = median(x), .by = shape)
ggplot(sk, aes(x)) +
geom_histogram(aes(y = after_stat(density)), bins = 60,
fill = "grey85", colour = NA) +
geom_density(colour = "grey20", linewidth = 0.8) +
geom_vline(data = sk_stats, aes(xintercept = mean, colour = "mean"),
linewidth = 0.9) +
geom_vline(data = sk_stats, aes(xintercept = median, colour = "median"),
linewidth = 0.9, linetype = "dashed") +
scale_colour_manual(values = c(mean = "firebrick", median = "steelblue")) +
facet_wrap(~ shape, scales = "free", ncol = 3) +
labs(title = "Skew moves the mean away from the median, in the direction of the tail",
x = NULL, y = "Density", colour = NULL) +
theme_dspa(10)
# --- Interactive equivalent ------------------------------------------------
xr <- rnbinom(10000, 5, 0.1); fr <- density(xr)
plot_ly(x = xr, type = "histogram", name = "Data histogram") |>
add_trace(x = fr$x, y = fr$y, type = "scatter", mode = "lines",
fill = "tozeroy", opacity = 0.3, yaxis = "y2",
name = "Density (rnbinom(N, 5, 0.1))") |>
layout(title = "Right Skewed Process",
yaxis2 = list(overlaying = "y", side = "right"),
legend = list(orientation = "h"))6 Reference distributions
6.1 Uniform and normal models
Under a uniform model every value in the support is equally likely, so the histogram is flat up to sampling noise. Under a normal model the density is unimodal, symmetric, and light-tailed.
set.seed(7)
u <- runif(1000, 1, 50)
g <- rnorm(1000, 0, 1)
p_u <- ggplot(data.frame(x = u), aes(x)) +
geom_histogram(aes(y = after_stat(density)), bins = 25,
fill = "lightblue", colour = "white") +
geom_hline(yintercept = 1 / 49, colour = "firebrick", linewidth = 1) +
labs(title = "Uniform(1, 50)", subtitle = "Red line: theoretical density 1/49",
x = NULL, y = "Density") + theme_dspa(10)
p_g <- ggplot(data.frame(x = g), aes(x)) +
geom_histogram(aes(y = after_stat(density)), bins = 30,
fill = "lightblue", colour = "white") +
stat_function(fun = dnorm, args = list(mean = 0, sd = 1),
colour = "firebrick", linewidth = 1) +
labs(title = "Normal(0, 1)", subtitle = "Red curve: theoretical density",
x = NULL, y = "Density") + theme_dspa(10)
p_u | p_g
# --- Interactive equivalents -----------------------------------------------
plot_ly(x = ~u, type = "histogram", histnorm = "probability",
name = "proportion", showlegend = FALSE) |>
add_lines(x = ~u, y = ~0.038, mode = "lines") |>
layout(bargap = 0.1, title = "Uniform(1, 50) Histogram",
yaxis = list(title = "probability"))
normDensity <- density(g, bw = 0.5)
dens <- data.frame(x = normDensity$x, y = normDensity$y)
plot_ly(dens) |>
add_histogram(x = g, name = "Normal Histogram") |>
add_lines(data = dens, x = ~x, y = ~y, yaxis = "y2",
line = list(width = 3), name = "N(0,1)") |>
layout(bargap = 0.1, title = "Normal(0,1)",
yaxis2 = list(overlaying = "y", side = "right",
range = c(0, max(dens$y)), showgrid = FALSE,
zeroline = FALSE),
legend = list(orientation = "h"))6.2 The
d/p/q/r family,
visually
For every distribution R supplies four functions. For the standard normal:
z <- seq(-4, 4, by = 0.01)
q <- seq(0.001, 0.999, by = 0.001)
std <- data.frame(z = z, density = dnorm(z), cdf = pnorm(z))
qtl <- data.frame(p = q, quantile = qnorm(q))
p1 <- ggplot(std, aes(z, density)) +
geom_area(fill = "steelblue", alpha = 0.35) + geom_line(linewidth = 0.9) +
labs(title = "dnorm, density", x = "z", y = "f(z)") + theme_dspa(10)
p2 <- ggplot(std, aes(z, cdf)) + geom_line(linewidth = 0.9, colour = "darkgreen") +
labs(title = "pnorm, CDF", x = "z", y = "F(z)") + theme_dspa(10)
p3 <- ggplot(qtl, aes(p, quantile)) + geom_line(linewidth = 0.9, colour = "firebrick") +
labs(title = "qnorm, inverse CDF", x = "p", y = expression(F^-1*(p))) + theme_dspa(10)
p1 | p2 | p3
# --- Interactive equivalents -----------------------------------------------
plot_ly(x = z, y = dnorm(z), mode = "lines", name = "Normal Density Curve") |>
layout(title = "Normal Density Curve",
xaxis = list(title = "critical values"),
yaxis = list(title = "Density"), legend = list(orientation = "h"))
plot_ly(x = z, y = pnorm(z), mode = "lines", name = "Normal CDF") |>
layout(title = "Normal Distribution",
xaxis = list(title = "critical values"),
yaxis = list(title = "Cumulative Distribution"),
legend = list(orientation = "h"))
plot_ly(x = q, y = qnorm(q), mode = "lines", name = "Normal Quantile Function") |>
layout(title = "Normal Quantile (Inverse CDF)",
xaxis = list(title = "probability values"),
yaxis = list(title = "Critical Values"),
legend = list(orientation = "h"))Explore further with the SOCR distribution calculators, the Distributome navigator, and the SOCR EBook probability chapter.
6.3 The 68–95–99.7 rule, and checking whether it applies
Under \(X\sim N(\mu,\sigma^2)\),
\[P(|X-\mu|\le k\sigma)=2\Phi(k)-1 =\begin{cases} 0.6827 & k=1\\ 0.9545 & k=2\\ 0.9973 & k=3. \end{cases}\]
zz <- seq(-4, 4, by = 0.001)
band <- function(k) data.frame(x = zz[abs(zz) <= k], k = paste0(k, "sd"))
areas <- bind_rows(band(3), band(2), band(1)) |>
mutate(y = dnorm(x),
k = factor(k, levels = c("3sd", "2sd", "1sd")))
ggplot() +
geom_area(data = areas, aes(x, y, fill = k), position = "identity") +
geom_line(data = data.frame(x = zz, y = dnorm(zz)), aes(x, y), linewidth = 0.9) +
geom_vline(xintercept = c(-3:-1, 1:3), linetype = "dashed",
colour = "grey45", linewidth = 0.4) +
scale_fill_manual(values = c("1sd" = "#F0A868", "2sd" = "#7FB3D5", "3sd" = "#A9DFBF"),
labels = c("1sd" = "68.3%", "2sd" = "95.4%", "3sd" = "99.7%")) +
scale_x_continuous(breaks = -3:3,
labels = c(expression(mu-3*sigma), expression(mu-2*sigma),
expression(mu-sigma), expression(mu),
expression(mu+sigma), expression(mu+2*sigma),
expression(mu+3*sigma))) +
labs(title = "The 68-95-99.7 rule", x = NULL, y = "Density", fill = NULL) +
theme_dspa()
# --- Interactive equivalent (annotated, with segment guides) ---------------
set.seed(7); norm <- rnorm(1000, 0, 1)
nd <- density(norm, bw = 0.5); dens <- data.frame(x = nd$x, y = nd$y)
xLabels <- c("μ-3σ","μ-2σ","μ-σ","μ",
"μ+σ","μ+2σ","μ+3σ")
labelColors <- c("green","red","orange","black","orange","red","green")
plot_ly(dens) |>
add_histogram(x = norm, name = "Normal Histogram") |>
add_lines(data = dens, x = ~x, y = ~y + 0.05, yaxis = "y2",
line = list(width = 3), name = "N(0,1)") |>
add_annotations(x = -3:3, y = 0.2, mode = "text", text = xLabels,
showarrow = FALSE, textposition = "middle right",
textfont = list(color = labelColors, size = 16)) |>
add_segments(x = -3, xend = 3, y = 100, yend = 100, name = "99.7%",
line = list(dash = "dash", color = "green")) |>
add_segments(x = -2, xend = 2, y = 90, yend = 90, name = "95%",
line = list(dash = "dash", color = "red")) |>
add_segments(x = -1, xend = 1, y = 80, yend = 80, name = "68%",
line = list(dash = "dash", color = "orange")) |>
layout(bargap = 0.1, title = "Normal 68-95-99.7% Rule",
yaxis = list(title = "density/frequency"),
yaxis2 = list(overlaying = "y", side = "right",
range = c(0, max(dens$y) + 0.1),
showgrid = FALSE, zeroline = FALSE),
legend = list(orientation = "h"))The rule is a consequence of normality, not a description of data. Before quoting \(\mu\pm k\sigma\) intervals, check whether normality holds, and report the empirical coverage alongside the theoretical one.
baseball <- dspa_read(
url = "https://umich.instructure.com/files/330381/download?download_frd=1",
name = "01a_data.txt",
reader = read.table, header = TRUE
)
str(baseball)#> 'data.frame': 1034 obs. of 6 variables:
#> $ Name : chr "Adam_Donachie" "Paul_Bako" "Ramon_Hernandez" "Kevin_Millar" ...
#> $ Team : chr "BAL" "BAL" "BAL" "BAL" ...
#> $ Position: chr "Catcher" "Catcher" "Catcher" "First_Baseman" ...
#> $ Height : int 74 74 72 72 73 69 69 71 76 71 ...
#> $ Weight : int 180 215 210 210 188 176 209 200 231 180 ...
#> $ Age : num 23 34.7 30.8 35.4 35.7 ...empirical_rule <- function(v, label) {
v <- v[!is.na(v)]; m <- mean(v); s <- sd(v)
data.frame(
variable = label,
mean = round(m, 1), sd = round(s, 1),
k = 1:3,
lower = round(m - 1:3 * s, 1),
upper = round(m + 1:3 * s, 1),
theoretical = round(2 * pnorm(1:3) - 1, 4),
empirical = round(sapply(1:3, \(k) mean(abs(v - m) <= k * s)), 4)
)
}
rbind(empirical_rule(baseball$Weight, "Weight (lb)"),
empirical_rule(baseball$Height, "Height (in)"))Rounding to one decimal is deliberate: quoting an interval endpoint
as 180.7168 lb, implies a precision that neither the scale
nor the normal approximation supports.
qq_panel <- function(v, label) {
v <- v[!is.na(v)]
ggplot(data.frame(s = v), aes(sample = s)) +
stat_qq(size = 0.7, alpha = 0.5) +
stat_qq_line(colour = "firebrick", linewidth = 0.8) +
labs(title = label, x = "Theoretical N(0,1) quantiles",
y = "Sample quantiles") + theme_dspa(10)
}
qq_panel(baseball$Weight, "MLB weight") | qq_panel(baseball$Height, "MLB height")
Weight shows mild right-tail departure; height is close to normal. That is what justifies, or qualifies, the intervals in the table above.
bb <- baseball |>
dplyr::select(Weight, Height) |>
pivot_longer(everything(), names_to = "measure", values_to = "value")
bb_par <- bb |> summarise(m = mean(value), s = sd(value), .by = measure)
ggplot(bb, aes(value)) +
geom_histogram(aes(y = after_stat(density)), bins = 30,
fill = "grey85", colour = "white") +
geom_density(colour = "steelblue", linewidth = 0.9) +
purrr::pmap(list(bb_par$m, bb_par$s, bb_par$measure), \(m, s, g)
stat_function(data = subset(bb, measure == g), fun = dnorm,
args = list(mean = m, sd = s),
colour = "firebrick", linetype = "dashed", linewidth = 0.8)) +
facet_wrap(~ measure, scales = "free") +
labs(title = "MLB player measurements vs. fitted normal models",
subtitle = "Blue: kernel density estimate. Red dashed: N(mean, sd) model",
x = NULL, y = "Density") + theme_dspa(11)
# --- Interactive equivalent ------------------------------------------------
xw <- rnorm(10000, mean(baseball$Weight, na.rm = TRUE), sd(baseball$Weight, na.rm = TRUE))
fitw <- density(xw, bw = 10)
plot_ly(x = ~baseball$Weight, type = "histogram", name = "Weight Histogram",
histnorm = "probability") |>
add_trace(x = ~fitw$x, y = ~5 * fitw$y, type = "scatter", mode = "lines",
opacity = 0.3, fill = "tozeroy", name = "Normal Density") |>
layout(title = "Baseball Weight Histogram & Model Normal Distribution",
xaxis = list(title = "Weight"),
yaxis = list(title = "relative frequency/density"),
legend = list(orientation = "h"))
xh <- rnorm(10000, mean(baseball$Height, na.rm = TRUE), sd(baseball$Height, na.rm = TRUE))
fith <- density(xh, bw = 1)
plot_ly(x = ~baseball$Height, type = "histogram", name = "Height Histogram",
histnorm = "probability") |>
add_trace(x = ~fith$x, y = ~fith$y, type = "scatter", mode = "lines",
opacity = 0.3, fill = "tozeroy", name = "Normal Density") |>
layout(title = "Baseball Height Histogram & Model Normal Distribution",
xaxis = list(title = "Height"),
yaxis = list(title = "relative frequency/density"),
legend = list(orientation = "h"))7 Categorical data and contingency tables
7.1 One-way tables
#>
#> 1990 1995 2000 2005 2010 2012
#> 520 561 570 570 556 554#>
#> Africa Americas Eastern Mediterranean
#> 797 613 373
#> Europe South-East Asia Western Pacific
#> 910 191 447#>
#> Rural Total Urban
#> 1095 1109 1127#>
#> 1990 1995 2000 2005 2010 2012
#> 15.6 16.8 17.1 17.1 16.7 16.6Note the argument: prop.table() takes a
table, not a data vector. Passing the raw vector
divides each observation by the sum of all observations and returns
something that is not a frequency distribution (Chapter 1, §1.7.12).
7.2 Two-way tables done correctly
Is residence-area composition different in the WHO African region than elsewhere?
water$africa <- factor(ifelse(water$region == "Africa", "Africa", "Rest of world"),
levels = c("Rest of world", "Africa"))
tab <- table(residence_area = water$residence_area, group = water$africa)
tab#> group
#> residence_area Rest of world Africa
#> Rural 828 267
#> Total 845 264
#> Urban 861 266Four quantities are routinely confused. They are different objects:
\[ \underbrace{E_{ij}=\frac{n_{i\cdot}n_{\cdot j}}{n}}_{\text{expected count}}, \qquad \underbrace{c_{ij}=\frac{(O_{ij}-E_{ij})^2}{E_{ij}}}_{\text{cell contribution}}, \qquad \underbrace{X^2=\sum_{i,j}c_{ij}}_{\text{test statistic}}, \qquad \underbrace{p=P(\chi^2_{(r-1)(c-1)}\ge X^2)}_{\text{p-value}} . \]
The cell contribution \(c_{ij}\) is not a probability. It is a squared standardized residual on \([0,\infty)\), unbounded above. This is not “the probability that the difference is due to chance”, which is a common misreading of contingency-table output.
For signed, approximately-\(N(0,1)\) diagnostics use the standardized Pearson residual
\[r_{ij}=\frac{O_{ij}-E_{ij}}{\sqrt{E_{ij}\,(1-p_{i\cdot})(1-p_{\cdot j})}},\]
which chisq.test() returns as $stdres.
Values beyond \(\pm 2\) mark cells
driving the association.
#>
#> Pearson's Chi-squared test
#>
#> data: tab
#> X-squared = 0.19964, df = 2, p-value = 0.905#>
#> Expected counts:#> group
#> residence_area Rest of world Africa
#> Rural 833.0 262.0
#> Total 843.7 265.3
#> Urban 857.3 269.7cat("\nCell contributions (O-E)^2/E:\n"); print(round((cs$observed - cs$expected)^2 / cs$expected, 3))#>
#> Cell contributions (O-E)^2/E:#> group
#> residence_area Rest of world Africa
#> Rural 0.030 0.096
#> Total 0.002 0.007
#> Urban 0.016 0.050#>
#> Standardized residuals:#> group
#> residence_area Rest of world Africa
#> Rural -0.432 0.432
#> Total 0.116 -0.116
#> Urban 0.314 -0.314Validity conditions. The \(\chi^2\) approximation requires all
expected counts \(\ge 1\) and at least
80% of them \(\ge 5\). Check, don’t
assume, and when it fails, use fisher.test() or
chisq.test(..., simulate.p.value = TRUE).
#> min_expected pct_expected_ge5
#> 261.9979 100.00007.3 Effect size, and why the p-value is not enough
With \(n \approx 3{,}300\), almost any departure from independence is “significant”. Report a magnitude. Cramér’s \(V\) rescales \(X^2\) to \([0,1]\):
\[V=\sqrt{\frac{X^2}{n\,\min(r-1,\,c-1)}} .\]
cramers_v <- function(tb) {
cs <- suppressWarnings(chisq.test(tb))
sqrt(unname(cs$statistic) / (sum(tb) * (min(dim(tb)) - 1)))
}
c(n = sum(tab),
X2 = round(unname(cs$statistic), 3),
p = signif(cs$p.value, 3),
Cramers_V = round(cramers_v(tab), 4))#> n X2 p Cramers_V
#> 3331.0000 0.2000 0.9050 0.0077tab_df <- as.data.frame(tab) |>
mutate(prop = Freq / sum(Freq), .by = group)
ggplot(tab_df, aes(group, prop, fill = residence_area)) +
geom_col(position = position_dodge(width = 0.75), width = 0.7) +
geom_text(aes(label = scales::percent(prop, accuracy = 0.1)),
position = position_dodge(width = 0.75), vjust = -0.35, size = 3) +
scale_y_continuous(labels = scales::percent, expand = expansion(c(0, 0.12))) +
labs(title = "Residence-area composition, Africa vs. rest of world",
subtitle = sprintf("Cramer's V = %.4f, statistically detectable, substantively small",
cramers_v(tab)),
x = NULL, y = "Share within group", fill = "Residence area") +
theme_dspa()
# --- Interactive equivalent, plus the classic gmodels cross-tabulation -----
gmodels::CrossTable(x = water$residence_area, y = water$africa)
plot_ly(tab_df, x = ~group, y = ~prop, color = ~residence_area, type = "bar") |>
layout(barmode = "group",
title = "Residence area by region group",
yaxis = list(title = "Share within group", tickformat = ".0%"),
xaxis = list(title = ""), legend = list(orientation = "h"))How to state the conclusion. Cramér’s \(V\) near \(0.02\) means the two groups’ residence-area compositions are, for practical purposes, indistinguishable — but that is an argument from effect size, not from a large p-value. “We failed to reject independence” never licenses “the groups are the same”: that is accepting the null, and it is the same error as reading a large Kolmogorov–Smirnov p-value as proof of distributional equality (Chapter 1, §1.12.5). If you want to claim equivalence, run an equivalence test with a pre-specified margin.
8 Bivariate relationships
w2 <- water |> filter(!is.na(improved_water), !is.na(sanitation_facilities))
ggplot(w2, aes(sanitation_facilities, improved_water)) +
geom_point(alpha = 0.25, size = 1.1, colour = "steelblue") +
geom_smooth(method = "loess", formula = y ~ x, se = TRUE,
colour = "firebrick", linewidth = 0.9) +
labs(title = "Improved water quality vs. sanitation facilities",
subtitle = sprintf("Pearson r = %.3f | Spearman rho = %.3f | n = %d",
cor(w2$sanitation_facilities, w2$improved_water),
cor(w2$sanitation_facilities, w2$improved_water,
method = "spearman"),
nrow(w2)),
x = "Sanitation facilities (%)", y = "Improved water (%)") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = ~water$sanitation_facilities, y = ~water$improved_water,
type = "scatter", mode = "markers",
text = ~paste0(water$country, " (", water$year, ")"),
marker = list(opacity = 0.4)) |>
layout(title = "Scatterplot: Improved Water Quality vs. Sanitation Facilities",
xaxis = list(title = "Sanitation facilities (%)"),
yaxis = list(title = "Improved water (%)"))Reporting both Pearson and Spearman is cheap insurance. Pearson measures linear association and is sensitive to outliers and to the ceiling effect visible here; Spearman measures monotone association and is invariant to any increasing transformation. A large gap between them is a signal, usually nonlinearity, ties, or a few influential points.
9 Missing data
Up to here we have written na.rm = TRUE and moved on.
That is a modeling decision disguised as a keyword argument, and it is
time to make it explicit.
9.1 A formal taxonomy
Let \(X=(X_{obs},X_{mis})\) denote the complete data and let \(R\) be the response indicator matrix, \(R_{ij}=1\) if \(X_{ij}\) is observed. The missingness mechanism is the conditional law \(P(R\mid X,\psi)\) (Rubin, 1976).
| Mechanism | Condition | Reading |
|---|---|---|
| MCAR | \(P(R\mid X,\psi)=P(R\mid \psi)\) | Missingness independent of all data, observed or not |
| MAR | \(P(R\mid X,\psi)=P(R\mid X_{obs},\psi)\) | Missingness depends only on what you can see |
| MNAR | \(P(R\mid X,\psi)\) depends on \(X_{mis}\) | Missingness depends on the unseen values themselves |
Nesting: \(\text{MCAR}\Rightarrow\text{MAR}\), never the converse.
The distinction is not academic, it determines what is estimable:
- Under MCAR, complete-case analysis is unbiased (merely inefficient).
- Under MAR, complete-case analysis is generally biased, but likelihood and multiple-imputation methods that condition on \(X_{obs}\) recover valid inference. The mechanism is then ignorable.
- Under MNAR, no method recovers the truth from the observed data alone. You must model the mechanism, and the model is unverifiable. All you can responsibly do is sensitivity analysis across plausible mechanisms.
Concretely: a lab machine that fails at random is MCAR; low-SES participants declining to report income when SES is recorded is MAR; heavier respondents declining to report weight is MNAR. You cannot test MAR against MNAR using the observed data, that is a theorem, not a limitation of current software.
9.2 Why single mean imputation is worse than it looks
The simplest fix replaces each missing entry with its column mean.
Proposition (variance deflation). Let \(x_1,\dots,x_n\) be a sample with \(n_{mis}\) values missing completely at random and \(n_{obs}=n-n_{mis}\) observed. Fill every missing entry with \(\bar x_{obs}\). Then the completed data have sample mean exactly \(\bar x_{obs}\) and sample variance \[\hat\sigma^2_{\text{imp}}=\frac{n_{obs}-1}{\,n-1\,}\,\hat\sigma^2_{obs}\;\approx\;(1-p)\,\sigma^{2},\qquad p=\frac{n_{mis}}{n}.\]
Proof. The completed mean is \(\frac{1}{n}\!\left(\sum_{obs}x_i+n_{mis}\bar x_{obs}\right)=\bar x_{obs}\). Each imputed cell therefore contributes \((\bar x_{obs}-\bar x_{obs})^2=0\) to the sum of squares, so \(\sum_{i=1}^n (x_i-\bar x)^2=\sum_{obs}(x_i-\bar x_{obs})^2=(n_{obs}-1)\hat\sigma^2_{obs}\). Dividing by \(n-1\) gives the result. \(\blacksquare\)
The consequences are mechanical, not probabilistic:
- Variance understated by a factor \(\approx (1-p)\).
- Standard errors understated by \(\approx\sqrt{1-p}\); confidence intervals too narrow; type-I error inflated.
- Covariances with every other variable attenuated by \(\approx (1-p)\), so correlations shrink toward zero and regression slopes are biased downward.
At the 30% missingness used in the simulation below, that is a 30% loss of variance and roughly 16% understatement of every standard error.
demo_deflation <- function(v, p, seed = 99) {
set.seed(seed)
v <- v[!is.na(v)]
n <- length(v)
idx <- sample(n, floor(p * n))
vm <- v; vm[idx] <- NA
vi <- vm; vi[is.na(vi)] <- mean(vm, na.rm = TRUE)
n_obs <- sum(!is.na(vm))
c(p = p,
var_complete = var(v),
var_imputed = var(vi),
ratio_observed = var(vi) / var(v),
ratio_predicted = (n_obs - 1) / (n - 1))
}
do.call(rbind, lapply(c(0.1, 0.2, 0.3, 0.5),
\(p) demo_deflation(water$improved_water, p))) |> round(4)#> p var_complete var_imputed ratio_observed ratio_predicted
#> [1,] 0.1 369.178 335.9256 0.9099 0.9002
#> [2,] 0.2 369.178 297.5748 0.8060 0.8002
#> [3,] 0.3 369.178 257.9415 0.6987 0.7001
#> [4,] 0.5 369.178 184.5310 0.4998 0.5000Observed ratio and predicted ratio agree to three decimals. The bias is deterministic and therefore correctable in principle, which is exactly what multiple imputation does (§2.7.6).
# Vectorized mean imputation. Note: no hard-coded row count.
impute_mean <- function(df, cols) {
for (cl in cols) {
m <- mean(df[[cl]], na.rm = TRUE)
df[[cl]][is.na(df[[cl]])] <- m
}
df
}
water_imp <- impute_mean(water, c("improved_water", "sanitation_facilities"))
rbind(original = sapply(water[c("improved_water", "sanitation_facilities")],
\(v) c(mean = mean(v, na.rm = TRUE), sd = sd(v, na.rm = TRUE))),
imputed = sapply(water_imp[c("improved_water", "sanitation_facilities")],
\(v) c(mean = mean(v), sd = sd(v)))) |> round(3)#> improved_water sanitation_facilities
#> mean 84.905 68.872
#> sd 19.214 31.353
#> mean 84.905 68.872
#> sd 19.121 30.7109.3 The modern missing-data toolbox
| Package | Approach | Best for |
|---|---|---|
naniar |
Tidy exploration and ggplot2-native missingness visuals | First look; gg_miss_var(),
gg_miss_upset() |
VIM |
Graphical diagnostics + kNN/hot-deck | Aggregation plots, matrix plots |
mice |
Multivariate Imputation by Chained Equations | The default choice under MAR |
Amelia |
Bootstrapped EM under multivariate normality | Fast MI for continuous, roughly-normal data |
missForest |
Random-forest, non-parametric | Mixed types, nonlinear structure |
missRanger |
ranger + predictive mean matching |
A much faster missForest |
missMDA |
Imputation inside PCA / MCA | When the goal is dimension reduction |
jomo |
Multilevel joint modeling | Clustered / hierarchical data |
rMCAIM |
Denoising autoencoders | Large, high-dimensional data |
mi |
Bayesian MI with rich metadata | Detailed per-variable model control |
Rule of thumb: mice under MAR with moderate
dimensionality; missForest / missRanger when
relationships are nonlinear or types are mixed; Amelia when
multivariate normality is plausible and speed matters.
9.4 Simulating a dataset with known missingness
set.seed(123)
# MCAR generator: blank a fixed proportion of each column, independently.
create_missing <- function(data, pct_mis = 10) {
stopifnot(pct_mis >= 0, pct_mis <= 100)
n <- nrow(data)
n_mis <- floor(n * pct_mis / 100)
for (j in seq_len(ncol(data))) {
if (n_mis > 0) data[sample(seq_len(n), n_mis), j] <- NA
}
as.data.frame(data)
}
n <- 1000
u1 <- rbinom(n, 1, 0.5); v1 <- log(rnorm(n, 5, 1)); x1 <- u1 * exp(v1)
u2 <- rbinom(n, 1, 0.5); v2 <- log(rnorm(n, 5, 1)); x2 <- u2 * exp(v2)
x3 <- rbinom(n, 1, 0.45)
x4 <- ordered(sample(1:5, n, replace = TRUE))
x5 <- sample(letters[1:10], n, replace = TRUE)
x6 <- trunc(runif(n, 1, 10))
x7 <- rnorm(n)
x8 <- factor(sample(1:10, n, replace = TRUE))
x9 <- runif(n, 0.1, 0.99)
x10 <- rpois(n, 4)
y <- x1 + x2 + x7 + x9 + rnorm(n)
sim_data <- data.frame(y, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10)
sim_30 <- create_missing(sim_data, pct_mis = 30)
colSums(is.na(sim_30))#> y x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
#> 300 300 300 300 300 300 300 300 300 300 300missForest::prodNA() does the same thing in one call;
writing it out makes the mechanism visible, which matters because
which mechanism you simulate determines which methods should
work.
library(naniar)
gg_miss_var(sim_30, show_pct = TRUE) +
labs(title = "Missingness by variable (MCAR, 30% per column)") +
theme_dspa(11)
vis_miss(sim_30, warn_large_data = FALSE) +
labs(title = "Missingness map: rows = observations, columns = variables") +
theme_dspa(10) + theme(axis.text.x = element_text(angle = 60, hjust = 0))
The map is the diagnostic. Under MCAR it looks like
static. Under MAR it shows horizontal banding, blocks
of rows missing together, and under MNAR the missing
cells concentrate at one end of another variable’s range. Compare
against mi::image(missing_data.frame(sim_30)) if you prefer
the classic rendering.
9.5 EM imputation
9.5.1 The algorithm
Expectation–Maximization maximizes a likelihood that cannot be maximized directly by alternating two steps. With observed data \(X\), latent/missing data \(Y\), and parameters \(\theta\), the marginal likelihood is
\[L(\theta\mid X)=p(X\mid\theta)=\int p(X,Y\mid\theta)\,dY,\]
usually intractable. EM instead iterates:
- E-step. \(\displaystyle Q(\theta\mid\theta^{(t)})=E_{Y\mid X,\theta^{(t)}}\big[\log L(\theta\mid X,Y)\big]\)
- M-step. \(\displaystyle \theta^{(t+1)}=\arg\max_\theta Q(\theta\mid\theta^{(t)})\)
Monotone ascent. \(\ell(\theta^{(t+1)}\mid X)\ge \ell(\theta^{(t)}\mid X)\) at every iteration.
The proof is a one-line consequence of Jensen’s inequality: the difference in observed-data log-likelihoods equals the increase in \(Q\) plus a Kullback–Leibler divergence, which is non-negative. So the likelihood never decreases, but it may converge to a local maximum or a saddle point, and the rate is only linear, governed by the fraction of missing information. More missingness, slower convergence.
9.5.2 The multivariate normal case
Assume \(Z\sim N_d(\mu,\Sigma)\). Partition each observation into observed (\(o\)) and missing (\(m\)) parts:
\[\mu=\begin{pmatrix}\mu_{o}\\ \mu_{m}\end{pmatrix},\qquad \Sigma=\begin{pmatrix}\Sigma_{oo} & \Sigma_{om}\\ \Sigma_{mo} & \Sigma_{mm}\end{pmatrix}.\]
The conditional distribution of the missing block given the observed block is Gaussian, which makes the E-step closed form:
\[\boxed{\;E(Y\mid X)=\mu_{m}+\Sigma_{mo}\Sigma_{oo}^{-1}\left(X-\mu_{o}\right)\;}\] \[\operatorname{Cov}(Y\mid X)=\Sigma_{mm}-\Sigma_{mo}\Sigma_{oo}^{-1}\Sigma_{om}\;\equiv\;C .\]
The M-step re-estimates from the completed sufficient statistics:
\[\mu^{(t+1)}=\frac1n\sum_{i=1}^n E(Z_i\mid X_i),\qquad \Sigma^{(t+1)}=\frac1n\sum_{i=1}^n\Big[\big(\hat z_i-\mu^{(t+1)}\big)\big(\hat z_i-\mu^{(t+1)}\big)^{\!\top}+C_i\Big].\]
The \(C_i\) term is not
optional. Dropping it, computing
Sigma <- var(imputed) on the filled-in matrix turns EM
into iterated conditional-mean imputation. That algorithm still
converges to a fixed point, but the fixed point is not
the MLE. Rather, it systematically under-estimates \(\Sigma\), for exactly the reason
established in §2.7.2. Conditional
means carry no conditional variance. We implement both so the size of
the error is visible.
9.5.3 Computational complexity
A naive implementation loops over rows, inverting \(\Sigma_{oo}\) for each:
\[\text{cost per iteration}=O\!\big(n\,d_{obs}^{3}\big).\]
But \(\Sigma_{oo}\) depends only on which entries are missing, not on the row’s values. Rows sharing a missingness pattern share a linear system. Group by pattern and the cost becomes
\[O\!\big(P\,d^{3}+n\,d^{2}\big),\]
with \(P\) the number of distinct patterns. Since \(P\ll n\) in essentially all real data, patterns arise from instrument failures and skip logic, not from individuals, this is a large constant-factor win and changes the practical ceiling on \(d\). It is the single most valuable optimization in the routine.
Also: use rcond() before inverting. \(\Sigma_{oo}\) becomes ill-conditioned as
\(d\) approaches \(n_{obs}\), and forcing a solve through with
tol = 1e-40 does not fix the conditioning, it only removes
the warning that would have told you the result is numerical noise.
9.5.4 Implementation
em_impute <- function(X, tol = 1e-4, max_iter = 500,
ridge = 1e-8, correct_cov = TRUE, verbose = FALSE) {
X <- as.matrix(X); storage.mode(X) <- "double"
n <- nrow(X); d <- ncol(X)
M <- is.na(X)
if (!any(M)) return(list(imputed = as.data.frame(X), iterations = 0L))
# OLD
# cc <- stats::complete.cases(X)
# if (sum(cc) <= d)
# stop("Need more complete cases than columns to initialise (have ",
# sum(cc), ", need > ", d, ").", call. = FALSE)
#
# mu <- colMeans(X[cc, , drop = FALSE])
# Sigma <- stats::var(X[cc, , drop = FALSE])
cc <- stats::complete.cases(X)
if (sum(cc) <= d) {
# Fallback: Initialize mu using available values per column,
# and use a regularized covariance or identity matrix
warning(
"Need more complete cases than columns to initialise (have ",
sum(cc), ", need > ", d, "). Fallback to initializing mu using available values per column, and using a regularized covariance or identity matrix ",
call. = FALSE
)
mu <- colMeans(X, na.rm = TRUE)
# Fill remaining NAs with column means for initial covariance calculation
X_temp <- X
for (j in seq_len(d)) X_temp[is.na(X_temp[, j]), j] <- mu[j]
Sigma <- stats::var(X_temp) + diag(ridge, d)
} else {
mu <- colMeans(X[cc, , drop = FALSE])
Sigma <- stats::var(X[cc, , drop = FALSE])
}
# --- group rows by missingness pattern: one linear solve per pattern -----
key <- apply(M, 1L, \(r) paste0(as.integer(r), collapse = ""))
patterns <- split(seq_len(n), key)
Z <- X
for (j in seq_len(d)) Z[M[, j], j] <- mu[j] # warm start at column means
delta_trace <- numeric(0); min_rcond <- Inf
for (it in seq_len(max_iter)) {
Z_old <- Z
C_accum <- matrix(0, d, d)
## ---------------- E step ----------------
for (rows in patterns) {
miss <- M[rows[1L], ]
if (!any(miss)) next
obs <- !miss
if (!any(obs)) { # entirely missing row
Z[rows, miss] <- rep(mu[miss], each = length(rows)); next
}
Soo <- Sigma[obs, obs, drop = FALSE]
Soo <- Soo + diag(ridge, nrow(Soo)) # tiny ridge for stability
rc <- rcond(Soo); min_rcond <- min(min_rcond, rc)
B <- Sigma[miss, obs, drop = FALSE] %*% solve(Soo) # Sigma_mo Sigma_oo^-1
Xo <- Z[rows, obs, drop = FALSE]
Z[rows, miss] <- sweep(sweep(Xo, 2L, mu[obs], "-") %*% t(B),
2L, mu[miss], "+")
if (correct_cov) { # conditional covariance C
Cmm <- Sigma[miss, miss, drop = FALSE] - B %*% Sigma[obs, miss, drop = FALSE]
C_accum[miss, miss] <- C_accum[miss, miss] + length(rows) * Cmm
}
}
## ---------------- M step ----------------
mu <- colMeans(Z)
Zc <- sweep(Z, 2L, mu, "-")
Sigma <- (crossprod(Zc) + C_accum) / n # MLE scaling (1/n)
delta <- max(abs(Z - Z_old))
delta_trace <- c(delta_trace, delta)
if (verbose) cat(sprintf("iter %3d max|dZ| = %.3e\n", it, delta))
if (delta < tol) break
}
list(imputed = as.data.frame(Z), mu = mu, Sigma = Sigma,
iterations = it, delta_trace = delta_trace,
min_rcond = min_rcond, correct_cov = correct_cov)
}set.seed(202227)
# Multivariate normal with heteroscedastic, uncorrelated features
p_dim <- 20; n_obs <- 200
mu_true <- rep(2, p_dim)
sig_true <- diag(1:p_dim)
sim_mvn <- MASS::mvrnorm(n_obs, mu_true, sig_true) +
MASS::mvrnorm(n_obs, rep(0, p_dim), diag(p_dim))
sim_complete <- sim_mvn
# Introduce ~500 MCAR cells out of 4,000 (sampling without replacement)
n_cells <- length(sim_mvn)
miss_idx <- sample(n_cells, 500)
sim_mvn[miss_idx] <- NA
sim_df <- as.data.frame(sim_mvn)
c(cells = n_cells, missing = sum(is.na(sim_df)),
pct = round(100 * mean(is.na(sim_df)), 2))#> cells missing pct
#> 4000.0 500.0 12.5Note that sampling without replacement gives exactly
500 missing cells. Sampling with replacement
(e1071::rdiscrete()) produces a random number
below 500 through repeated indices, which makes the missingness
rate itself a random quantity for no benefit.
fit_corr <- em_impute(sim_df, tol = 1e-4, correct_cov = TRUE)
fit_naive <- em_impute(sim_df, tol = 1e-4, correct_cov = FALSE)
c(iterations_corrected = fit_corr$iterations,
iterations_naive = fit_naive$iterations,
min_rcond = signif(fit_corr$min_rcond, 3))#> iterations_corrected iterations_naive min_rcond
#> 18.0000 27.0000 0.0218conv <- data.frame(
iteration = c(seq_along(fit_corr$delta_trace), seq_along(fit_naive$delta_trace)),
delta = c(fit_corr$delta_trace, fit_naive$delta_trace),
variant = rep(c("with conditional-covariance correction", "conditional-mean only"),
c(length(fit_corr$delta_trace), length(fit_naive$delta_trace))))
ggplot(conv, aes(iteration, delta, colour = variant)) +
geom_line(linewidth = 0.9) + geom_point(size = 1.2) +
scale_y_log10() +
labs(title = "EM convergence: monotone, linear rate",
subtitle = "Vertical axis is log-scaled; a straight line indicates linear convergence",
x = "Iteration", y = expression(max~group("|",Delta*Z,"|")), colour = NULL) +
theme_dspa()
Does the covariance correction matter? Compare each variant’s estimated variances against the truth we simulated:
var_true <- diag(sig_true) + 1 # signal variance + noise variance
var_corr <- diag(fit_corr$Sigma)
var_naive <- diag(fit_naive$Sigma)
data.frame(
feature = 1:p_dim,
true = var_true,
corrected = round(var_corr, 2),
naive = round(var_naive, 2),
bias_corr = round(var_corr - var_true, 2),
bias_naive = round(var_naive - var_true, 2)
) |> head(10)c(mean_relative_bias_corrected = round(mean((var_corr - var_true) / var_true), 4),
mean_relative_bias_naive = round(mean((var_naive - var_true) / var_true), 4))#> mean_relative_bias_corrected mean_relative_bias_naive
#> 0.0340 -0.0704The conditional-mean-only version is biased downward across the
board, the same deflation proved in §2.7.2, now appearing inside an
algorithm that looks far more sophisticated. Sophistication is
not protection. From here on we use fit_corr.
sim_imputed <- fit_corr$imputed
plot_pair <- function(i, j) {
was_missing <- is.na(sim_df[[i]]) | is.na(sim_df[[j]])
d <- data.frame(x = sim_imputed[[i]], y = sim_imputed[[j]],
status = ifelse(was_missing, "imputed", "observed"))
ggplot(d, aes(x, y)) +
stat_ellipse(type = "norm", colour = "#000099", alpha = 0.6, linewidth = 0.5) +
geom_point(aes(colour = status, size = status, alpha = status)) +
scale_colour_manual(values = c(observed = "grey25", imputed = "magenta")) +
scale_size_manual(values = c(observed = 0.7, imputed = 1.8)) +
scale_alpha_manual(values = c(observed = 0.45, imputed = 0.95)) +
labs(x = paste0("V", i), y = paste0("V", j), colour = NULL) +
guides(size = "none", alpha = "none") + theme_dspa(10)
}
((plot_pair(1, 2) | plot_pair(5, 6)) / (plot_pair(13, 20) | plot_pair(18, 19))) +
plot_layout(guides = "collect") +
plot_annotation(title = "EM-imputed values sit inside the joint distribution",
subtitle = "Grey: observed. Magenta: at least one coordinate imputed. Ellipse: 95% normal contour",
theme = theme_dspa())
# --- Interactive equivalent ------------------------------------------------
was_missing <- is.na(sim_df$V1) | is.na(sim_df$V5)
plot_ly(sim_imputed, x = ~V1, y = ~V5, type = "scatter", mode = "markers",
color = ~ifelse(was_missing, "imputed", "observed"),
colors = c(observed = "grey30", imputed = "magenta"),
marker = list(size = 7)) |>
layout(title = "EM imputation: V1 vs V5",
xaxis = list(title = "V1"), yaxis = list(title = "V5"))Notice what imputation cannot do: every magenta point lies close to the conditional mean surface, so the imputed cloud is less dispersed than the observed cloud. A single completed dataset always understates uncertainty — which is the entire motivation for the next section.
9.6 Multiple imputation and Rubin’s rules
Single imputation treats an estimate as if it were data. Multiple imputation generates \(m\) completed datasets that differ by exactly the amount the missing data are uncertain, analyzes each, and pools, propagating that uncertainty into the standard errors.
Let \(\hat Q^{(k)}\) and \(U^{(k)}\) be the estimate and its variance from imputation \(k\). Rubin’s rules give
\[\bar Q=\frac1m\sum_{k=1}^m \hat Q^{(k)},\qquad \bar U=\frac1m\sum_{k=1}^m U^{(k)},\qquad B=\frac{1}{m-1}\sum_{k=1}^m\big(\hat Q^{(k)}-\bar Q\big)^2,\]
\[\boxed{\;T=\underbrace{\bar U}_{\text{within}}+\underbrace{\Big(1+\tfrac1m\Big)B}_{\text{between}}\;}\]
with approximate degrees of freedom \(\nu=(m-1)\big(1+\frac{\bar U}{(1+1/m)B}\big)^2\). The fraction of missing information is \(\lambda=\frac{(1+1/m)B}{T}\), the share of total uncertainty attributable to missingness, and the number you should report.
The \((1+1/m)\) factor is the finite-\(m\) correction; the rule of thumb is to take \(m\) at least as large as the percentage of incomplete cases (White et al., 2011).
library(mice)
set.seed(4242)
sim_small <- sim_30[, c("y", "x1", "x2", "x7", "x9", "x10")]
imp <- mice(sim_small, m = 5, maxit = 10, printFlag = FALSE, seed = 4242)
fit <- with(imp, lm(y ~ x1 + x2 + x7 + x9 + x10))
pooled <- pool(fit)
summary(pooled) |> as.data.frame() |>
transform(estimate = round(estimate, 4), std.error = round(std.error, 4),
statistic = round(statistic, 3), p.value = signif(p.value, 3))# Rubin's decomposition, made explicit
pool_detail <- pooled$pooled[, c("term", "m", "estimate", "ubar", "b", "t", "dfcom", "riv", "fmi")]
pool_detail |>
transform(estimate = round(estimate, 4), ubar = signif(ubar, 3),
b = signif(b, 3), t = signif(t, 3),
riv = round(riv, 3), fmi = round(fmi, 3))Read the last column. fmi is the fraction of missing
information per coefficient; riv is the relative increase
in variance due to nonresponse. If fmi is far above your
missingness rate, the imputation model is weak for that variable.
library(mice)
# This will correctly dispatch to mice's internal plot.mids method during a knit
plot(imp, layout = c(2, 2))
Chain traces should mix without trend. Systematic drift means
maxit is too small. Convergence is assessed with \(\hat R\); note that the modern threshold is
\(\hat R<1.01\) (Vehtari et al., 2021), not
the \(1.1\) quoted in older
references.
# --- The classic `mi` workflow, retained for reference ---------------------
library(betareg); library(mi)
mdf <- missing_data.frame(sim_30)
show(mdf)
image(mdf) # missingness pattern
mdf <- change(mdf, y = "x1", what = "imputation_method", to = "pmm")
imputations <- mi(mdf, n.iter = 30, n.chains = 5, verbose = FALSE)
round(mipply(imputations, mean, to.matrix = TRUE), 3)
Rhats(imputations, statistic = "moments") # target < 1.01
plot(imputations); hist(imputations); image(imputations)
model_results <- pool(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10,
data = imputations, m = 5)
display(model_results); summary(model_results)9.7 Comparing imputation methods
Amelia
draws multiple imputations from a bootstrapped EM under multivariate
normality, a natural benchmark for our own EM.
library(Amelia)
set.seed(505)
am <- amelia(sim_df, m = 5, p2s = 0)
am_imputed <- am$imputations[[5]]
dim(am_imputed)#> [1] 200 20# NOTE: earlier editions hard-coded aes(X1, X2) for the Amelia layer, so the
# two methods were plotted on DIFFERENT variables. Both layers now use i, j.
compare_pair <- function(i, j) {
miss <- is.na(sim_df[[i]]) | is.na(sim_df[[j]])
base <- data.frame(x = sim_imputed[[i]], y = sim_imputed[[j]])
em <- data.frame(x = sim_imputed[[i]][miss], y = sim_imputed[[j]][miss])
amp <- data.frame(x = am_imputed[[i]][miss], y = am_imputed[[j]][miss])
ggplot(base, aes(x, y)) +
geom_point(alpha = 0.35, size = 0.7, colour = "grey30") +
stat_ellipse(type = "norm", colour = "#000099", alpha = 0.6, linewidth = 0.5) +
geom_point(data = em, aes(x, y, colour = "EM (manual)"), size = 2.2, shape = 16) +
geom_point(data = amp, aes(x, y, colour = "Amelia"), size = 2.2, shape = 18) +
scale_colour_manual(values = c("EM (manual)" = "magenta", "Amelia" = "#FF9933")) +
labs(x = paste0("V", i), y = paste0("V", j), colour = NULL) +
theme_dspa(10)
}
(compare_pair(2, 4) | compare_pair(17, 18)) +
plot_layout(guides = "collect") +
plot_annotation(title = "EM vs. Amelia imputations, same variables in both layers",
theme = theme_dspa())
dens_compare <- function(idx) {
miss <- is.na(sim_df[[idx]])
bind_rows(
data.frame(value = sim_df[[idx]][!miss], source = "Observed"),
data.frame(value = sim_imputed[[idx]][miss], source = "EM (manual)"),
data.frame(value = am_imputed[[idx]][miss], source = "Amelia")
) |> mutate(feature = paste0("V", idx))
}
bind_rows(dens_compare(5), dens_compare(9), dens_compare(10)) |>
ggplot(aes(value, colour = source, fill = source)) +
geom_density(alpha = 0.18, linewidth = 0.9) +
facet_wrap(~ feature, scales = "free", ncol = 3) +
labs(title = "Observed vs. imputed distributions",
subtitle = "Under MCAR these should agree in location; imputed spread is narrower by construction",
x = NULL, y = "Density", colour = NULL, fill = NULL) +
theme_dspa(11)
# --- Interactive equivalent ------------------------------------------------
my_plotly <- function(index) {
miss <- is.na(sim_df[[index]])
observed <- sim_df[[index]][!miss]
em <- sim_imputed[[index]][miss]
amelia_i <- am_imputed[[index]][miss]
plot_ly() |>
add_lines(x = ~density(em)$x, y = ~density(em)$y, name = "EM") |>
add_lines(x = ~density(amelia_i)$x, y = ~density(amelia_i)$y, name = "Amelia") |>
add_lines(x = ~density(observed)$x, y = ~density(observed)$y, name = "Observed") |>
layout(title = sprintf("Distributions: Feature V%d", index),
xaxis = list(title = "Measurements"),
yaxis = list(title = "Densities"),
legend = list(title = list(text = "Distributions"), orientation = "h"))
}
my_plotly(5); my_plotly(9); my_plotly(10)The imputed densities are centerd correctly and narrower than the observed ones, the visual signature of the deflation theorem. It is why you pool \(m\) imputations rather than trusting one.
9.8 Case study: traumatic brain injury
tbi <- dspa_read(
url = "https://umich.instructure.com/files/720782/download?download_frd=1",
name = "08_EpiBioSData_Incomplete.csv",
na.strings = c("", ".", "NA")
)
dim(tbi)#> [1] 46 19#> id age sex mechanism field.gcs er.gcs
#> 0 0 0 0 2 2
#> icu.gcs worst.gcs X6m.gose X2013.gose skull.fx temp.injury
#> 1 1 5 0 0 0
#> surgery spikes.hr min.hr max.hr acute.sz late.sz
#> 0 18 18 18 0 0
#> ever.sz
#> 0Clinical assessment scores (GCS, GOSE) are documented in DOI 10.1080/02699050701727460; the raw table is on the SOCR wiki.
The UpSet plot is strictly more informative than a marginal bar chart: it shows which variables go missing together. Co-occurring blocks point to a shared cause, an instrument, a protocol step, a skip pattern, and that is evidence about the mechanism, which is what determines whether MAR is defensible.
tbi_vars <- c("ever.sz", "surgery", "worst.gcs", "sex", "age", "spikes.hr")
tbi_sub <- tbi[, intersect(tbi_vars, names(tbi))]
tbi_sub$sex <- factor(tbi_sub$sex)
set.seed(808)
tbi_imp <- mice(tbi_sub, m = 5, maxit = 10, printFlag = FALSE, seed = 808)
tbi_fit <- with(tbi_imp, glm(ever.sz ~ surgery + worst.gcs + sex + age,
family = binomial()))
summary(pool(tbi_fit)) |> as.data.frame() |>
transform(estimate = round(estimate, 3), std.error = round(std.error, 3),
p.value = signif(p.value, 3))The pooled logistic regression is fit five times and combined by Rubin’s rules — never fit once on a single completed dataset, which would treat imputed values as observations.
10 Class imbalance and cohort rebalancing
10.1 Why imbalance is a problem, and when it isn’t
With cohorts of 434 patients and 166 controls, a classifier that predicts “patient” for everyone attains 72% accuracy while learning nothing. Imbalance causes trouble because:
- Accuracy stops being informative. Use balanced accuracy, \(F_1\), MCC, or AUC-PR instead (§2.16).
- The decision boundary migrates toward the minority class under most loss functions, since minority errors contribute less total loss.
- Variance estimates degrade for minority-class parameters, simply because \(n\) is small there.
Rule of thumb: worry when one class is more than an order of magnitude larger than another. Do not resample reflexively, if your metric is threshold-free (AUC) and your model outputs calibrated probabilities, rebalancing can hurt by distorting those probabilities.
10.2 SMOTE
SMOTE (Synthetic Minority Over-sampling TEchnique) creates new minority observations by interpolating between existing ones. For a minority point \(x_i\), pick one of its \(k\) nearest minority neighbours \(x_{nn}\) and generate
\[x_{\text{new}} = x_i + \delta\,(x_{nn}-x_i),\qquad \delta\sim\text{Uniform}(0,1).\]
The synthetic point lies on the segment joining two real minority points, so it inherits their local structure. Complexity is \(O(n_{min}\,\log n_{min}\cdot d)\) with a k-d tree, \(O(n_{min}^2 d)\) brute force.
Two properties follow directly from the construction and matter enormously:
- Every synthetic point is a convex combination of two real points.
- Therefore, if those two real points end up in different folds of a train/test split, a near-duplicate of a test observation is sitting in training.
10.3 The leakage trap
Rule. Resample inside the training fold. Never before the split.
Running ubBalance() on the full dataset and then
cross-validating has a measurable consequence, as we show below.
ppmi <- dspa_read(
url = "https://umich.instructure.com/files/330400/download?download_frd=1",
name = "06_PPMI_ClassificationValidationData_Short.csv",
header = TRUE
)
table(ppmi$ResearchGroup)#>
#> Control PD SWEDD
#> 166 434 61ppmi_base <- ppmi |>
filter(VisitID == 1) |>
mutate(PD = factor(ifelse(ResearchGroup == "Control", "Control", "Patient"),
levels = c("Control", "Patient"))) |>
dplyr::select(-any_of(c("ResearchGroup", "X", "FID_IID", "VisitID")))
ppmi_base <- ppmi_base[, c(which(sapply(ppmi_base, is.numeric)), which(names(ppmi_base) == "PD"))]
ppmi_base <- ppmi_base[stats::complete.cases(ppmi_base), ]
c(n = nrow(ppmi_base), p = ncol(ppmi_base) - 1)#> n p
#> 422 100#>
#> Control Patient
#> 122 300library(themis) # SMOTE that composes with recipes / tidymodels
library(recipes)
library(rsample)
set.seed(2718)
fit_eval <- function(train, test) {
m <- suppressWarnings(
glm(PD ~ ., data = train, family = binomial(),
control = list(maxit = 60)))
p <- predict(m, newdata = test, type = "response")
pred <- factor(ifelse(p > 0.5, "Patient", "Control"),
levels = c("Control", "Patient"))
tp <- sum(pred == "Patient" & test$PD == "Patient")
tn <- sum(pred == "Control" & test$PD == "Control")
sens <- tp / sum(test$PD == "Patient")
spec <- tn / sum(test$PD == "Control")
c(accuracy = mean(pred == test$PD),
balanced_accuracy = (sens + spec) / 2)
}
smote_it <- function(df) {
recipe(PD ~ ., data = df) |>
step_smote(PD, over_ratio = 1) |>
prep() |> juice() |> as.data.frame()
}
folds <- vfold_cv(ppmi_base, v = 5, strata = PD)
## --- WRONG: resample the whole dataset, then split ------------------------
ppmi_leaked <- smote_it(ppmi_base)
folds_leaked <- vfold_cv(ppmi_leaked, v = 5, strata = PD)
res_leaked <- vapply(folds_leaked$splits,
\(s) fit_eval(analysis(s), assessment(s)), numeric(2))
## --- RIGHT: split first, resample only the training fold ------------------
res_clean <- vapply(folds$splits, \(s) {
tr <- smote_it(analysis(s)) # SMOTE applied to TRAINING data only
fit_eval(tr, assessment(s))
}, numeric(2))
rbind(
`SMOTE before split (LEAKY)` = round(rowMeans(res_leaked), 4),
`SMOTE inside fold (correct)` = round(rowMeans(res_clean), 4)
)#> accuracy balanced_accuracy
#> SMOTE before split (LEAKY) 0.9850 0.9850
#> SMOTE inside fold (correct) 0.9716 0.9703The leaky protocol reports the better number. It is not the better model, it is the same model, evaluated on data it has partially memorized. Because Chapters 5, 6, 9, and 11 reuse this pattern, the correction belongs here.
Diagnostic habit. If resampled performance is dramatically better than performance on a held-out set you never touched, suspect leakage before celebrating.
cmp <- bind_rows(
data.frame(value = ppmi_base[[1]], source = "Original"),
data.frame(value = ppmi_leaked[[1]], source = "After SMOTE")
)
p_dens <- ggplot(cmp, aes(value, colour = source, fill = source)) +
geom_density(alpha = 0.2, linewidth = 0.9) +
labs(title = "Marginal distribution", x = names(ppmi_base)[1], y = "Density",
colour = NULL, fill = NULL) + theme_dspa(10)
qq <- as.data.frame(qqplot(ppmi_base[[1]], ppmi_leaked[[1]], plot.it = FALSE))
p_qq <- ggplot(qq, aes(x, y)) +
geom_point(size = 0.8, alpha = 0.6) +
geom_abline(slope = 1, intercept = 0, colour = "firebrick", linetype = "dashed") +
labs(title = "Q-Q: original vs. rebalanced", x = "Original", y = "Rebalanced") +
theme_dspa(10)
p_dens | p_qq
# --- Interactive equivalent ------------------------------------------------
QQ <- qqplot(ppmi_base[[1]], ppmi_leaked[[1]], plot.it = FALSE)
plot_ly(x = ~QQ$x, y = ~QQ$y, type = "scatter", mode = "markers",
showlegend = FALSE) |>
add_lines(x = range(QQ$x), y = range(QQ$x), showlegend = FALSE) |>
layout(title = "QQ-Plot Original vs. Rebalanced Data",
xaxis = list(title = "original data"),
yaxis = list(title = "rebalanced data"))Whne contrasting the original against rebalanced values, recall that the two samples are not independent, as the rebalanced set contains the original majority observations verbatim plus synthetic points interpolated from the original minority observations. So, we need to compare distributions descriptively (as above), or test on a genuinely held-out sample.
Alternatives. Class weights
(glm(weights=), ranger(class.weights=))
reweight the loss without fabricating data. Threshold tuning on a
validation set often achieves what resampling attempts. ADASYN,
Borderline-SMOTE, and Tomek-link cleaning are refinements available in
themis and smotefamily. And note that SMOTE
interpolates in Euclidean space, so it is poorly defined for categorical
features, SMOTE-NC exists for that case.
11 PART II, EXPLORATORY VISUAL ANALYTICS
Part I asked whether the data are trustworthy. Part II asks what they say.
Exploratory data analysis is not a decorative stage before “real” modeling. It is the stage where you discover that a variable is bimodal, that two instruments disagree, that a ceiling effect makes your linear model untenable, all things no summary statistic reports and no model will warn you about.
Reminder on graphics policy (§2.0): every 2-D figure is
ggplot2, with theplot_ly()equivalent printed but not evaluated; every 3-D figure is liveplot_ly().
12 The grammar and the perception of graphics
12.1 What questions drive a visualization?
- What exploratory techniques let me interrogate this data structure?
- How do I examine paired associations in a multivariate dataset?
- Which display makes the comparison I actually care about easy to read accurately?
The third question has an empirical answer, and it is the one usually skipped.
12.2 Classifying visualization methods
Visualizations can be organized along several axes:
- Data type, structured/unstructured, small/large, complete/incomplete, temporal/spatial, Euclidean/non-Euclidean.
- Task type, how the analyst, the data, and the display interact.
- Scalability, how many observations the technique tolerates before it saturates.
- Dimensionality, how many attributes are encoded simultaneously.
- Positioning, whether relative placement itself carries meaning (as in correlation displays).
- Investigative need, composition, distribution, comparison, or relationship. This is the axis we organize by.
12.3 The grammar of graphics
ggplot2 implements Wilkinson’s grammar of
graphics: a plot is a mapping from data through an
aesthetic mapping to geometric
objects, under scales, optionally split by
facets, positioned by a coordinate
system, and styled by a theme.
\[\text{plot}=\text{data}+\underbrace{\text{aes}(x,y,\text{colour},\dots)}_{\text{mapping}}+\underbrace{\text{geom}}_{\text{marks}}+\text{scales}+\text{facets}+\text{coord}+\text{theme}\]
The payoff is compositional: once you know the grammar, an unfamiliar plot type is a recombination rather than a new API.
12.4 Cleveland–McGill: not all encodings are equal
Cleveland and McGill (1984) measured how accurately people extract quantitative information from different visual encodings. The resulting ordering, most accurate first:
| Rank | Elementary perceptual task | Encodes value as |
|---|---|---|
| 1 | Position along a common scale | Dot plot, scatterplot, bar height on a shared axis |
| 2 | Position along non-aligned scales | Small multiples / facets |
| 3 | Length, direction, angle | Unaligned bars, slopes |
| 4 | Area | Bubble size, treemap |
| 5 | Volume, curvature | 3-D bars |
| 6 | Shading, colour saturation | Heatmap, choropleth |
Two operational consequences:
- Prefer position. A dot plot or a bar on a shared axis is read more accurately than a pie chart (angle/area) or a heatmap (shading).
- Colour is for categories and for coarse magnitude, not for values a reader must compare precisely. When a heatmap is the right structure, add the numbers, or provide the underlying values in a table.
This ranking is why the pie charts in §2.11.2 come with a caveat and a replacement, and why the association heatmaps in §2.11.3 are rescaled and given a visible colour bar.
13 Composition
13.1 Histograms and density
A histogram estimates the distribution of a sample; the theory is in §2.4.2. The graphical form long predates the name: William Playfair used bar charts in 1786 and André-Michel Guerry in 1833. Karl Pearson coined the term “histogram” in his 1891 lectures (in print, 1895) — naming, rather than inventing, the device.
The instructive display superimposes three things: the sample histogram, a kernel density estimate, and the theoretical density of the model you believe generated the data.
set.seed(1234)
N <- 10000; mu <- 15; sd <- 3.7
x <- rnorm(N, mean = mu, sd = sd)
ggplot(data.frame(x = x), aes(x)) +
geom_histogram(aes(y = after_stat(density)), bins = nclass.FD(x),
fill = "grey85", colour = "white") +
geom_density(aes(colour = "Kernel density estimate"), linewidth = 1) +
stat_function(aes(colour = "Theoretical N(15, 3.7)"), fun = dnorm,
args = list(mean = mu, sd = sd), linewidth = 1, linetype = "dashed") +
scale_colour_manual(values = c("Kernel density estimate" = "steelblue",
"Theoretical N(15, 3.7)" = "firebrick")) +
labs(title = "Sample histogram, kernel density estimate, and model density",
subtitle = sprintf("n = %s; Freedman-Diaconis bins = %d",
format(N, big.mark = ","), nclass.FD(x)),
x = NULL, y = "Density", colour = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
fit <- density(x)
z <- seq(mu - 4 * sd, mu + 4 * sd, 0.1)
plot_ly(x = x, type = "histogram", name = "Data Histogram",
histnorm = "probability") |>
add_trace(x = fit$x, y = fit$y, type = "scatter", mode = "lines",
opacity = 0.3, fill = "tozeroy", yaxis = "y2",
name = "Kernel density estimate") |>
add_trace(x = z, y = dnorm(z, mu, sd), type = "scatter", mode = "lines",
opacity = 0.3, fill = "tozeroy", yaxis = "y2",
name = "Normal(15, 3.7)") |>
layout(title = "Data Histogram, Density Estimate & Theoretical Model",
yaxis2 = list(overlaying = "y", side = "right"),
legend = list(orientation = "h"))13.2 Pie charts, and what to use instead
Pie charts encode value as angle and area, ranks 3–4 on the Cleveland–McGill scale. They are readable for two or three categories of very different size and become unreadable quickly after that. The SOCR Latin letter frequency data has 26 categories, which is the textbook case against.
We show both, because seeing the failure is more convincing than being told about it.
library(rvest)
letter <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_LetterFrequencyData") |>
html_nodes("table") |> _[[1]] |> html_table()
str(letter[, 1:6])#> tibble [27 × 6] (S3: tbl_df/tbl/data.frame)
#> $ Letter : chr [1:27] "a" "b" "c" "d" ...
#> $ English : num [1:27] 0.08 0.01 0.03 0.04 0.13 0.02 0.02 0.06 0.07 0 ...
#> $ French : num [1:27] 0.08 0.01 0.03 0.04 0.15 0.01 0.01 0.01 0.08 0.01 ...
#> $ German : num [1:27] 0.07 0.02 0.03 0.05 0.17 0.02 0.03 0.05 0.08 0 ...
#> $ Spanish : num [1:27] 0.13 0.01 0.05 0.06 0.14 0.01 0.01 0.01 0.06 0 ...
#> $ Portuguese: num [1:27] 0.15 0.01 0.04 0.05 0.13 0.01 0.01 0.01 0.06 0 ...langs <- c("English", "Spanish", "Swedish", "Polish")
lf <- letter |>
dplyr::select(Letter, all_of(langs)) |>
pivot_longer(-Letter, names_to = "language", values_to = "freq") |>
mutate(freq = as.numeric(freq)) |>
filter(!is.na(freq))
p_pie <- ggplot(lf, aes(x = "", y = freq, fill = Letter)) +
geom_col(width = 1, colour = "white", linewidth = 0.1) +
coord_polar(theta = "y") +
facet_wrap(~ language, nrow = 2) +
labs(title = "26-slice pie charts: which letter is 4th most common in Polish?") +
theme_void(base_size = 10) + theme(legend.position = "none",
plot.title = element_text(face = "bold"))
p_bar <- lf |>
mutate(Letter = forcats::fct_reorder(Letter, freq, .fun = sum)) |>
ggplot(aes(freq, Letter, fill = language)) +
geom_col(position = position_dodge(width = 0.8), width = 0.75) +
labs(title = "Same data as position on a common scale",
x = "Relative frequency", y = NULL, fill = NULL) +
theme_dspa(9)
p_pie / p_bar + plot_layout(heights = c(1, 1.35))
The question in the pie title is essentially unanswerable; on the bar chart it takes about a second. If a reader must rank or compare magnitudes, use position. Reserve pies for part-to-whole with very few, very unequal parts.
# --- Interactive equivalent: the original 2x2 pie and donut panels ---------
plot_ly(letter, labels = ~Letter, values = ~English, type = "pie",
name = "English", textposition = "inside", textinfo = "label+percent",
showlegend = FALSE, domain = list(row = 0, column = 0)) |>
add_pie(labels = ~Letter, values = ~Spanish, name = "Spanish",
textposition = "inside", textinfo = "label+percent",
showlegend = FALSE, domain = list(row = 0, column = 1)) |>
add_pie(labels = ~Letter, values = ~Swedish, name = "Swedish",
textposition = "inside", textinfo = "label+percent",
showlegend = FALSE, domain = list(row = 1, column = 0)) |>
add_pie(labels = ~Letter, values = ~Polish, name = "Polish",
textposition = "inside", textinfo = "label+percent",
showlegend = FALSE, domain = list(row = 1, column = 1)) |>
add_annotations(x = 0.01, y = 0.99, text = "English", showarrow = FALSE) |>
add_annotations(x = 0.58, y = 0.99, text = "Spanish", showarrow = FALSE) |>
add_annotations(x = 0.01, y = 0.01, text = "Swedish", showarrow = FALSE) |>
add_annotations(x = 0.58, y = 0.01, text = "Polish", showarrow = FALSE) |>
layout(title = "Pie Charts of English, Spanish, Swedish & Polish Letters",
grid = list(rows = 2, columns = 2),
xaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
yaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE))
# Donut variant (hole = 0.5) for German / Italian / Dutch / Esperanto
plot_ly(letter, labels = ~Letter, values = ~German, type = "pie", name = "German",
textposition = "inside", textinfo = "label+percent", showlegend = FALSE,
domain = list(row = 0, column = 0), hole = 0.5) |>
add_pie(labels = ~Letter, values = ~Italian, name = "Italian",
domain = list(row = 0, column = 1), hole = 0.5, showlegend = FALSE) |>
add_pie(labels = ~Letter, values = ~Dutch, name = "Dutch",
domain = list(row = 1, column = 0), hole = 0.5, showlegend = FALSE) |>
add_pie(labels = ~Letter, values = ~Esperanto, name = "Esperanto",
domain = list(row = 1, column = 1), hole = 0.5, showlegend = FALSE) |>
layout(title = "Pie Charts of German, Italian, Dutch & Esperanto Letters",
grid = list(rows = 2, columns = 2))13.3 Heat maps done responsibly
Heat maps encode magnitude as shading, the least accurate channel. They earn their place when the structure of a matrix matters more than any individual cell.
Our example is a neuroimaging-genetics study: rows are SNPs, columns are brain regions of interest, and each cell is the p-value of that SNP–ROI association, for three cohorts, Alzheimer’s disease (AD), mild cognitive impairment (MCI), and normal controls (NC).
read_assoc <- function(id, nm)
dspa_read(sprintf("https://umich.instructure.com/files/%s/download?download_frd=1", id),
nm, reader = read.table, header = TRUE, row.names = 1,
sep = ",", dec = ".")
AD <- read_assoc("330387", "AD_SNP_ROI_Assoc_P_values.csv")
MCI <- read_assoc("330390", "MCI_SNP_ROI_Assoc_P_values.csv")
NC <- read_assoc("330391", "NC_SNP_ROI_Assoc_P_values.csv")
dim(AD) # rows = SNPs, columns = ROIs#> [1] 28 45Three corrections over the naive rendering.
First, scale. p-values live on \((0,1]\) and everything interesting is crushed near zero. Plot \(-\log_{10} p\), on which \(p=0.05\mapsto 1.3\) and \(p=10^{-6}\mapsto 6\).
Second, multiplicity. This is a \(28 \times 45\) grid of tests. At \(\alpha=0.05\) you expect 5% of cells to look “significant” under a complete null. Apply Benjamini–Hochberg and mark the threshold.
Third, show the colour bar. Hiding it,
hide_colorbar(), removes the only way to read a value off a
shading encoding.
tidy_assoc <- function(mat, label) {
m <- as.matrix(mat); storage.mode(m) <- "double"
data.frame(
SNP = rep(rownames(m), times = ncol(m)),
ROI = rep(colnames(m), each = nrow(m)),
p = as.vector(m),
cohort = label
)
}
assoc <- bind_rows(tidy_assoc(AD, "AD"), tidy_assoc(MCI, "MCI"), tidy_assoc(NC, "NC")) |>
filter(!is.na(p)) |>
mutate(q = p.adjust(p, method = "BH"), .by = cohort) |>
mutate(neglog10p = -log10(pmax(p, .Machine$double.xmin)))
# BH threshold on the -log10 scale, per cohort
bh_line <- assoc |>
summarise(thr = { s <- sort(p); k <- sum(s <= 0.05 * seq_along(s) / length(s))
if (k > 0) -log10(s[k]) else NA_real_ }, .by = cohort)
bh_lineggplot(assoc, aes(ROI, SNP, fill = neglog10p)) +
geom_raster() +
scale_fill_viridis_c(option = "magma", direction = -1,
name = expression(-log[10](p))) +
facet_wrap(~ cohort, nrow = 1) +
labs(title = "SNP x ROI association strength by cohort",
subtitle = "Darker = stronger association. Colour bar shown deliberately; see BH thresholds above.",
x = "ROI imaging biomarker", y = "SNP") +
theme_dspa(9) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 5),
axis.text.y = element_text(size = 4),
legend.position = "right")
assoc |>
summarise(tests = dplyr::n(),
raw_p_lt_05 = sum(p < 0.05),
expected_by_chance = round(0.05 * dplyr::n()),
BH_q_lt_05 = sum(q < 0.05),
.by = cohort)Compare raw_p_lt_05 against
expected_by_chance. If they are close, the apparent signal
is what a complete null would produce. BH_q_lt_05 is the
count that survives multiplicity control, the number worth talking
about.
# --- Interactive equivalents (one per cohort) ------------------------------
AD_mat <- as.matrix(AD); class(AD_mat) <- "numeric"
MCI_mat <- as.matrix(MCI); class(MCI_mat) <- "numeric"
NC_mat <- as.matrix(NC); class(NC_mat) <- "numeric"
heat <- function(mat, ttl)
plot_ly(x = ~colnames(mat), y = ~rownames(mat),
z = ~-log10(pmax(mat, .Machine$double.xmin)),
type = "heatmap", colorbar = list(title = "-log10(p)")) |>
layout(title = ttl,
xaxis = list(title = "ROI Imaging Biomarkers"),
yaxis = list(title = "SNPs"))
heat(AD_mat, "AD Neuroimaging-Genomic Associations")
heat(MCI_mat, "MCI Neuroimaging-Genomic Associations")
heat(NC_mat, "(Normal) HC Neuroimaging-Genomic Associations")Differences in the block structure across cohorts suggest different genetic traits, or different brain regions, are implicated in the three clinical groups. That is a hypothesis-generating observation, the heatmap does not test it.
14 Comparison
14.1 Scatterplots
A scatterplot places each observation at \((x_i,y_i)\), position on a common scale, the most accurate encoding available. When neither variable is clearly the response, either assignment is legitimate and the display shows association, not causation.
set.seed(21)
N <- 50
sc <- data.frame(ind = 1:N, x = runif(N), y = runif(N), z = runif(N))
ggplot(sc[1:20, ], aes(x, y, colour = z)) +
geom_point(size = 3) +
scale_colour_viridis_c(name = "z") +
labs(title = "Random scatterplot", x = "X", y = "Y") +
theme_dspa()
# --- Interactive equivalent (hover text carries the point identity) --------
hoverText <- paste0("Point ", sc$ind, ": (", round(sc$x, 3), ", ", round(sc$y, 3), ")")
plot_ly(x = ~sc$x[1:20], y = ~sc$y[1:20], type = "scatter", mode = "markers",
color = ~sc$z[1:20], size = 2, text = hoverText[1:20]) |>
layout(title = "Random Scatterplot",
xaxis = list(title = "X"), yaxis = list(title = "Y")) |>
hide_colorbar()14.2 Scatterplot matrices
With \(d\) variables there are \(\binom{d}{2}\) pairwise plots. A SPLOM shows them all at once; the cost is \(O(d^2)\) panels, so it saturates around \(d\approx 8\).
set.seed(31)
Nn <- 1000
dfm <- data.frame(
x = runif(Nn), y = runif(Nn), z = runif(Nn),
w = rnorm(Nn), u = rpois(Nn, 1.7),
class = sample(LETTERS[1:3], Nn, replace = TRUE, prob = c(0.2, 0.5, 0.3))
)
GGally::ggpairs(dfm, columns = 1:5, aes(colour = class, alpha = 0.4),
upper = list(continuous = GGally::wrap("cor", size = 2.6)),
lower = list(continuous = GGally::wrap("points", size = 0.4))) +
theme_dspa(8)
Because a SPLOM is fundamentally about brushing, selecting points in one panel and seeing where they land in the others, this is a case where the interactive version carries genuinely more information, so it is evaluated:
pl_colorscale <- list(c(0.0, "#19d3f3"), c(0.333, "#19d3f3"),
c(0.333, "#e763fa"), c(0.666, "#e763fa"),
c(0.666, "#636efa"), c(1, "#636efa"))
axis_style <- list(showline = FALSE, zeroline = FALSE,
gridcolor = "#ffff", ticklen = 4)
plot_ly(dfm) |>
add_trace(
type = "splom",
dimensions = list(list(label = "X", values = ~x), list(label = "Y", values = ~y),
list(label = "Z", values = ~z), list(label = "w", values = ~w),
list(label = "U", values = ~u)),
text = ~class,
marker = list(color = as.integer(factor(dfm$class)),
colorscale = pl_colorscale, size = 5,
line = list(width = 1, color = "rgb(230,230,230)"))) |>
layout(title = "Random Data Pairs Plot (drag to brush and link)",
hovermode = "closest", dragmode = "select",
plot_bgcolor = "rgba(240,240,240,0.95)",
xaxis = axis_style, yaxis = axis_style,
xaxis2 = axis_style, xaxis3 = axis_style, xaxis4 = axis_style,
yaxis2 = axis_style, yaxis3 = axis_style, yaxis4 = axis_style)14.2.1 A real example: mental-health services survey
The 2010
National Mental Health Services Survey covers 10,374 facilities. Two
variables of interest: supp (count of specialty and support
services offered) and qual (count of quality indicators
present).
mh <- dspa_read(
"https://umich.instructure.com/files/399128/download?download_frd=1",
"Case_MentalHealthServices.txt", reader = read.table, header = TRUE)
dim(mh)#> [1] 10374 10#> STFIPS majorfundtype FacilityType Ownership Focus
#> 0 4522 0 0 630
#> PostTraum GLBT num qual supp
#> 1275 1304 1594 2190 2099Note we compute missingness rather than discovering it through a
warning from pairs(). And note there is no
attach() here. Attaching a data frame and then
reassigning it leaves the search path pointing at a stale copy, so bare
column names silently resolve to old values.
mh_cc <- mh |> filter(!is.na(qual), !is.na(supp))
ggplot(mh_cc, aes(qual, supp)) +
geom_point(alpha = 0.12, size = 0.9, colour = "steelblue") +
geom_smooth(method = "loess", formula = y ~ x, span = 0.7,
colour = "firebrick", fill = "grey70", linewidth = 1) +
labs(title = "Support services vs. quality indicators",
subtitle = sprintf("LOESS (span = 0.7) with 95%% confidence band; n = %s complete cases",
format(nrow(mh_cc), big.mark = ",")),
x = "Quality indicators", y = "Support services") +
theme_dspa()
# --- Interactive equivalents: scatter, SPLOM, and LOESS with ribbon --------
plot_ly(mh, x = ~qual, y = ~supp, type = "scatter", mode = "markers",
size = 2, color = ~num, text = ~STFIPS) |>
layout(title = "2010 National Mental Health Services Survey",
xaxis = list(title = "Quality Indicators"),
yaxis = list(title = "Support Services")) |>
hide_colorbar()
plot_ly(mh) |>
add_trace(type = "splom",
dimensions = list(list(label = "FacilityType", values = ~FacilityType),
list(label = "Ownership", values = ~Ownership),
list(label = "Focus", values = ~Focus),
list(label = "PostTraum", values = ~PostTraum),
list(label = "num", values = ~num)),
text = ~STFIPS,
marker = list(color = as.integer(mh$qual), colorscale = pl_colorscale,
size = 7, line = list(width = 1, color = mh$qual))) |>
layout(title = "Mental Health Services Survey Pairs Plot (color = qual)",
hovermode = "closest", dragmode = "select",
plot_bgcolor = "rgba(240,240,240,0.95)")
ll <- loess(supp ~ qual, data = mh_cc, span = 0.7)
pred <- predict(ll, se = TRUE)
ord <- order(mh_cc$qual)
plot_ly(x = mh_cc$qual, y = mh_cc$supp, type = "scatter", mode = "markers",
name = "Data", marker = list(opacity = 0.2)) |>
add_lines(x = mh_cc$qual[ord], y = pred$fit[ord], name = "LOESS mean",
line = list(color = "gray", width = 4)) |>
add_ribbons(x = mh_cc$qual[ord],
ymin = (pred$fit - 1.96 * pred$se.fit)[ord],
ymax = (pred$fit + 1.96 * pred$se.fit)[ord],
name = "95% CI",
line = list(opacity = 0.4, width = 1, color = "lightgray")) |>
layout(title = "LOESS Model (supp ~ qual) with Confidence Band",
xaxis = list(title = "Quality Indicator"),
yaxis = list(title = "Support Services"))14.3 Jitter: seeing through overplotting
When many observations share coordinates, integer or rounded data especially — points stack invisibly. Two fixes: jitter (add small random displacement) and alpha (make points translucent so density reads as darkness).
quake <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_Dinov_021708_Earthquakes") |>
html_nodes("table") |> _[[2]] |> html_table()
str(quake)#> tibble [168 × 13] (S3: tbl_df/tbl/data.frame)
#> $ Date_(YYYY/MM/DD): chr [1:168] "1969/10/02" "1969/10/02" "1972/02/24" "1974/11/28" ...
#> $ Time : chr [1:168] "04:56:45.30" "06:19:56.39" "15:56:50.99" "23:01:24.59" ...
#> $ Latitude : num [1:168] 38.5 38.5 36.6 36.9 40.5 ...
#> $ Longitude : num [1:168] -123 -123 -121 -121 -124 ...
#> $ Depth : num [1:168] 0.22 5.14 4.18 5.48 23.48 ...
#> $ Mag : num [1:168] 5.6 5.7 5.1 5.2 5.3 5.7 5.1 5.2 5 6.3 ...
#> $ Magt : chr [1:168] "ML" "ML" "ML" "ML" ...
#> $ Nst : int [1:168] 38 53 10 51 15 24 7 7 13 13 ...
#> $ Gap : int [1:168] 104 139 128 61 176 116 102 146 319 318 ...
#> $ Clo : int [1:168] 52 58 6 4 5 14 2 2 91 73 ...
#> $ RMS : num [1:168] 0.22 0.22 0.06 0.13 0.04 0.07 0.06 0.03 0.15 0.09 ...
#> $ SRC : chr [1:168] "NCSN" "NCSN" "NCSN" "NCSN" ...
#> $ EventID : int [1:168] -1003132 -1003135 -1009260 -1021953 -1024134 -1025002 -1025150 -1025155 -1027524 -1032452 ...p_raw <- ggplot(quake, aes(Depth, Latitude, colour = Magt)) +
geom_point(size = 1.4) +
labs(title = "No jitter, no transparency", x = "Depth (km)", y = "Latitude") +
theme_dspa(10)
p_jit <- ggplot(quake, aes(Depth, Latitude, colour = Magt)) +
geom_point(position = position_jitter(width = 0.3, height = 0.3),
alpha = 0.45, size = 1.4) +
labs(title = "Jitter (0.3) + alpha (0.45)", x = "Depth (km)", y = "Latitude") +
theme_dspa(10)
p_raw | p_jit
# --- Interactive equivalent: symbol-coded magnitude types ------------------
glyph_for <- function(name)
dplyr::case_when(name == "Md" ~ "diamond-open", name == "ML" ~ "circle-open",
name == "Mw" ~ "square-open", name == "Mx" ~ "x-open",
TRUE ~ "triangle-up")
quake$glyph <- glyph_for(quake$Magt)
plot_ly(quake) |>
add_markers(x = ~Longitude, y = ~Latitude, color = ~Magt,
marker = list(size = ~Depth, symbol = ~glyph,
line = list(color = "black", width = 2)),
text = ~paste0("Mag ", Mag, " / depth ", Depth, " km")) |>
layout(title = "California Earthquakes (1969-2007)",
xaxis = list(title = "Longitude"), yaxis = list(title = "Latitude"))14.4 Bar charts and error bars
ggplot2::diamonds holds 53,940 records with price
($326–$18,823), carat, cut, colour (D best → J worst), clarity (I1 worst
→ IF best), dimensions, depth, and table.
dm <- ggplot2::diamonds
p_box <- ggplot(dm, aes(cut, log(price), fill = cut)) +
geom_boxplot(outlier.alpha = 0.05, width = 0.6) +
labs(title = "log(price) by cut", x = "Cut", y = "log(price)") +
theme_dspa(10) + theme(legend.position = "none")
p_grp <- ggplot(dm, aes(clarity, log(price), fill = color)) +
geom_boxplot(outlier.shape = NA, width = 0.75) +
labs(title = "log(price) by clarity and colour", x = "Clarity", y = "log(price)",
fill = "Colour") +
theme_dspa(10)
p_box / p_grp
# --- Interactive equivalents ----------------------------------------------
plot_ly(dm, x = ~cut, y = ~price, type = "bar", color = ~clarity, text = ~clarity)
plot_ly(dm, y = ~log(price), color = ~cut, type = "box") |>
layout(title = "Boxplot of Diamond (log) Price by Cut",
xaxis = list(title = "Diamond Cut"))
plot_ly(dm, x = ~clarity, y = ~log(price), color = ~color, type = "box") |>
layout(boxmode = "group",
title = "Grouped Boxplot of Diamond (log) Price by Clarity and Color",
legend = list(title = list(text = "<b> Diamond Color </b>")),
xaxis = list(title = "Diamond Clarity"))
# For a jittered box: boxpoints = "all", jitter = 0.3, pointpos = -1.814.4.1 Child trauma: group means with dispersion
Case 04 examines post-traumatic psychopathology and service utilization among trauma-exposed children.
trauma <- dspa_read(
"https://umich.instructure.com/files/399129/download?download_frd=1",
"Case_04_ChildTrauma.txt", reader = read.table, header = TRUE)
str(trauma)#> 'data.frame': 1000 obs. of 9 variables:
#> $ id : int 1 2 3 4 5 6 7 8 9 10 ...
#> $ sex : int 1 1 0 0 1 0 0 1 0 1 ...
#> $ age : int 6 14 6 11 7 9 12 9 9 13 ...
#> $ ses : int 0 0 0 0 0 0 0 0 1 0 ...
#> $ race : chr "black" "black" "black" "black" ...
#> $ traumatype: chr "sexabuse" "sexabuse" "sexabuse" "sexabuse" ...
#> $ ptsd : int 1 0 0 0 1 1 1 0 1 1 ...
#> $ dissoc : int 1 0 1 1 1 0 1 1 1 0 ...
#> $ service : int 17 12 9 11 15 6 9 10 11 13 ...trauma_stats <- trauma |>
summarise(n = dplyr::n(),
mean_service = mean(service, na.rm = TRUE),
sd_service = sd(service, na.rm = TRUE),
se_service = sd_service / sqrt(n),
.by = c(traumatype, race))
head(trauma_stats, 8)ggplot(trauma_stats, aes(traumatype, mean_service, fill = race)) +
geom_col(position = position_dodge(width = 0.85), width = 0.78) +
geom_errorbar(aes(ymin = mean_service - se_service,
ymax = mean_service + se_service),
position = position_dodge(width = 0.85), width = 0.2,
linewidth = 0.4) +
labs(title = "Mean service utilization by trauma type and race",
subtitle = "Error bars: +/- 1 standard error of the mean (not SD)",
x = "Trauma type", y = "Mean services used", fill = "Race") +
theme_dspa()
Say what the bars mean. \(\pm 1\) SD describes the spread of individuals; \(\pm 1\) SE describes the precision of the mean and is narrower by \(\sqrt{n}\). An unlabelled error bar is uninterpretable, and the two are routinely confused in published figures.
# --- Interactive equivalent ------------------------------------------------
plot_ly(data = subset(trauma_stats, race == "black"),
x = ~traumatype, y = ~mean_service, type = "bar", name = "Black",
error_y = ~list(array = se_service, color = "#000000")) |>
add_trace(data = subset(trauma_stats, race == "hispanic"), name = "Hispanic") |>
add_trace(data = subset(trauma_stats, race == "other"), name = "Other") |>
add_trace(data = subset(trauma_stats, race == "white"), name = "White") |>
layout(title = "Statistical Barplots (Child Trauma Dataset)",
legend = list(title = list(text = "<b> Race </b>")))ggplot(trauma, aes(race, fill = race)) +
geom_bar() +
facet_grid(. ~ traumatype) +
labs(title = "Counts of trauma type by race", x = NULL, y = "Count") +
theme_dspa(10) +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")
14.5 Trees and dendrograms
A graph is an ordered pair \(G=(V,E)\) of vertices and edges; a tree is a connected acyclic graph. Hierarchical clustering builds a tree by repeatedly merging the two closest clusters.
Agglomerative clustering with \(n\)
objects requires the full \(\binom{n}{2}\) distance matrix, \(O(n^2)\) memory, and naive implementations
cost \(O(n^3)\) time, reducible to
\(O(n^2\log n)\) with a priority queue
and to \(O(n^2)\) for single/complete
linkage via SLINK/CLINK. The quadratic memory requirement is
what limits hclust() in practice, not the
time.
The linkage function determines what “closest” means:
\[ d_{\text{single}}(A,B)=\min_{a\in A,b\in B}d(a,b),\quad d_{\text{complete}}(A,B)=\max_{a\in A,b\in B}d(a,b),\quad d_{\text{average}}(A,B)=\frac{1}{|A||B|}\sum_{a,b}d(a,b). \]
Single linkage chains; complete linkage produces compact spheres; average (UPGMA) is a compromise. Ward’s method minimizes within-cluster variance and is usually the best default for numeric data.
nof1 <- dspa_read("https://umich.instructure.com/files/330385/download?download_frd=1",
"02_Nof1_Data.csv", sep = ",", header = TRUE)
dim(nof1); head(nof1, 3)#> [1] 900 10hc <- hclust(dist(scale(nof1)), method = "average")
mem <- cutree(hc, k = 10)
cent <- do.call(rbind, lapply(1:10, \(k) colMeans(nof1[mem == k, , drop = FALSE])))
hc1 <- hclust(dist(cent), method = "average", members = table(mem))
hc1$labels <- paste("Cluster", 1:10)
dd <- ggdendro::dendro_data(hc1, type = "rectangle")
ggplot() +
geom_segment(data = ggdendro::segment(dd),
aes(x, y, xend = xend, yend = yend), linewidth = 0.5) +
geom_text(data = ggdendro::label(dd), aes(x, y, label = label),
hjust = 1, angle = 90, size = 3, nudge_y = -0.15) +
expand_limits(y = -1.2) +
labs(title = "Hierarchical clustering restarted from 10 cluster centroids",
subtitle = "Average linkage on standardized features",
x = NULL, y = "Merge height") +
theme_dspa(10) +
theme(axis.text.x = element_blank(), panel.grid.major.x = element_blank())
Because dendrograms invite interrogation, which subject sits in which cluster, at what height, the interactive version is evaluated.
plot_ly() |>
add_segments(data = ggdendro::segment(dd),
x = ~x, y = ~y, xend = ~xend, yend = ~yend,
line = list(color = "black"), hoverinfo = "none",
showlegend = FALSE) |>
layout(title = "Re-start from 10 clusters (hover and zoom)",
xaxis = list(title = "", tickvals = seq_along(dd$labels$label),
ticktext = dd$labels$label),
yaxis = list(title = "Merge height"))14.6 Correlation displays
# NOTE the margin. The file is read with row.names = 1, so ROWS are SNPs and
# COLUMNS are ROIs. cor() operates on COLUMNS, so this is an ROI x ROI matrix.
NC_mat <- as.matrix(NC); storage.mode(NC_mat) <- "double"
M_roi <- cor(NC_mat, use = "pairwise.complete.obs") # ROI x ROI
M_snp <- cor(t(NC_mat), use = "pairwise.complete.obs") # SNP x SNP
dim(M_roi); dim(M_snp)#> [1] 47 47#> [1] 36 36M_long <- as.data.frame(as.table(M_roi)) |>
setNames(c("ROI_1", "ROI_2", "r"))
ggplot(M_long, aes(ROI_1, ROI_2, fill = r)) +
geom_raster() +
scale_fill_gradient2(low = "#B2182B", mid = "white", high = "#2166AC",
midpoint = 0, limits = c(-1, 1), name = "r") +
labs(title = "Correlation between ROIs of their SNP-association profiles",
subtitle = "cor() works on columns: this is ROI x ROI, not SNP x SNP",
x = NULL, y = NULL) +
theme_dspa(8) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 5),
axis.text.y = element_text(size = 5))
corrplot::corrplot(M_roi, method = "ellipse", type = "upper", order = "hclust",
diag = FALSE, tl.cex = 0.45, tl.col = "black",
mar = c(0, 0, 1, 0), title = "ROI correlation (hclust-ordered)")
corrplot offers "circle",
"square", "ellipse", "number",
"shade", "color", and "pie"
glyphs, plus corrplot.mixed(). Two choices matter more than
the glyph: order = "hclust", which
permutes rows and columns so correlated blocks become visible, and a
diverging palette centerd at zero, so sign is readable
at a glance.
15 Relationships
15.1 Line plots
Line charts connect successive points, which asserts that the intervening values are meaningful. Only use them when the \(x\) axis has a natural order, time, dose, position. Connecting unordered categories draws a trend that does not exist.
ggplot(quake, aes(Longitude, Latitude, colour = Magt, shape = Magt)) +
geom_point(size = 2.2, alpha = 0.75) +
labs(title = "California earthquakes, 1969-2007",
subtitle = "Magnitude type (Magt) as colour and glyph",
x = "Longitude", y = "Latitude", colour = "Magt", shape = "Magt") +
coord_fixed(ratio = 1) +
theme_dspa()
# --- Interactive equivalent -----------------------------------------------
plot_ly(quake) |>
add_markers(x = ~Longitude, y = ~Latitude, color = ~Magt,
marker = list(size = ~Depth, symbol = ~glyph_for(Magt),
line = list(color = "black", width = 2))) |>
layout(title = "California Earthquakes (1969 - 2007)")15.2 Density comparison across groups
ggplot(quake, aes(Latitude, colour = Magt, fill = Magt)) +
geom_density(alpha = 0.12, linewidth = 1) +
geom_vline(xintercept = 37.5, linetype = "dashed", colour = "grey35") +
annotate("text", x = 37.6, y = Inf, label = "37.5 deg N", hjust = 0, vjust = 1.6,
size = 3, colour = "grey35") +
labs(title = "Latitude distribution by earthquake magnitude type",
x = "Latitude (degrees N)", y = "Density", colour = NULL, fill = NULL) +
theme_dspa()
The local (ML) magnitude type peaks near 37–38 degrees
north, a concentration that reflects the seismograph network’s
coverage as much as the underlying seismicity. A density is a
picture of the sampling process as much as of the
phenomenon.
# --- Interactive equivalent -----------------------------------------------
dens_by <- lapply(split(quake$Latitude, quake$Magt), density)
p <- plot_ly()
for (nm in names(dens_by))
p <- add_lines(p, x = dens_by[[nm]]$x, y = dens_by[[nm]]$y, name = nm)
p |> layout(title = "Latitude density by magnitude type",
xaxis = list(title = "Latitude"), yaxis = list(title = "Density"),
legend = list(orientation = "h"))15.3 LOESS: local
regression, and what span controls
LOESS (locally estimated scatterplot smoothing) fits a low-degree polynomial in a moving neighbourhood, weighting nearby points more heavily.
At a target \(x_0\), let \(q=\lceil \alpha n\rceil\) be the
neighbourhood size (\(\alpha\) =
span) and \(d(x_0)\) the
distance to the \(q\)-th nearest
neighbour. Each observation receives the tricube
weight
\[w_i(x_0)=\begin{cases}\left(1-\left|\dfrac{x_i-x_0}{d(x_0)}\right|^{3}\right)^{3}, & |x_i-x_0|<d(x_0)\\[2mm] 0,&\text{otherwise,}\end{cases}\]
and the local fit solves the weighted least-squares problem
\[\hat\beta(x_0)=\arg\min_{\beta}\sum_{i=1}^{n}w_i(x_0)\Big(y_i-\beta_0-\beta_1(x_i-x_0)-\cdots-\beta_p(x_i-x_0)^p\Big)^2,\]
with \(\hat
g(x_0)=\hat\beta_0(x_0)\). Degree \(p=1\) or \(2\); loess() defaults to
2.
span is the bias–variance dial. Small
\(\alpha\) → few points per fit → low
bias, high variance (wiggly). Large \(\alpha\) → many points → high bias, low
variance (over-smoothed, approaching a global polynomial). Because LOESS
is a linear smoother, \(\hat
y = L y\), its effective degrees of freedom are \(\operatorname{tr}(L)\), and it can be tuned
by generalized cross-validation.
Cost: one weighted least-squares solve per evaluation point, so \(O(n\,q\,p^2)\) in general;
loess() uses a k-d tree with vertex interpolation
(surface = "interpolate") to keep this practical.
spans <- c(0.15, 0.4, 0.75)
ggplot(mh_cc, aes(qual, supp)) +
geom_point(alpha = 0.10, size = 0.8, colour = "grey40") +
purrr::map(spans, \(s)
geom_smooth(method = "loess", formula = y ~ x, se = FALSE,
span = s, aes(colour = factor(s)), linewidth = 1)) +
scale_colour_brewer(palette = "Set1", name = "span") +
labs(title = "LOESS: span is the bias-variance dial",
subtitle = "Small span tracks noise; large span flattens real structure",
x = "Quality indicators", y = "Support services") +
theme_dspa()
sapply(spans, \(s) {
m <- loess(supp ~ qual, data = mh_cc, span = s)
c(span = s, enp = round(m$enp, 2), residual_sd = round(m$s, 4))
}) |> t()#> span enp residual_sd
#> [1,] 0.15 14.15 4.9596
#> [2,] 0.40 9.87 4.9598
#> [3,] 0.75 5.31 4.9624enp is the equivalent number of parameters, the
effective degrees of freedom spent. Watch it fall as span
rises: that is the bias–variance tradeoff with a number attached.
16 Distributions and distributional modeling
16.1 Empirical and theoretical, discrete and continuous
There is a duality between the theoretical mass/density/distribution functions and their empirical counterparts. For a discrete process the theoretical object is a probability mass function and must be evaluated only at its support points.
set.seed(1234)
pois_sample <- rpois(1000, lambda = 1)
emp <- as.data.frame(table(pois_sample)) |>
setNames(c("k", "count")) |>
mutate(k = as.integer(as.character(k)), prop = count / sum(count))
# dpois() must be evaluated at INTEGERS. Feeding it seq(0, 6, by = 0.01)
# returns 0 with a warning at every non-integer point.
theo <- data.frame(k = 0:6, pmf = dpois(0:6, lambda = 1))
kde <- density(pois_sample, kernel = "gaussian")
ggplot() +
geom_col(data = emp, aes(k, prop, fill = "Sample relative frequency"),
width = 0.65, alpha = 0.75) +
geom_point(data = theo, aes(k, pmf, colour = "Theoretical Poisson(1) mass"),
size = 3) +
geom_segment(data = theo, aes(x = k, xend = k, y = 0, yend = pmf,
colour = "Theoretical Poisson(1) mass"),
linewidth = 0.9) +
geom_line(data = data.frame(x = kde$x, y = kde$y),
aes(x, y, colour = "Gaussian KDE (inappropriate here)"),
linewidth = 0.9, linetype = "dashed") +
scale_fill_manual(values = c("Sample relative frequency" = "#9EC5E8")) +
scale_colour_manual(values = c("Theoretical Poisson(1) mass" = "magenta",
"Gaussian KDE (inappropriate here)" = "steelblue")) +
labs(title = "Poisson(1): sample, theoretical mass, and a kernel density estimate",
subtitle = "The KDE smooths a DISCRETE distribution into a continuous one - a category error",
x = "k", y = "Probability", fill = NULL, colour = NULL) +
theme_dspa()
The dashed curve is included as a warning. A Gaussian kernel density estimate applied to count data spreads mass onto non-integers where the distribution places none, and it places positive density below zero. For discrete data use the empirical mass function; if you need smoothing, use a discrete kernel.
# --- Interactive equivalent -----------------------------------------------
h <- hist(pois_sample, breaks = seq(-0.5, max(pois_sample) + 0.5, by = 1), plot = FALSE)
plot_ly(x = h$mids, y = h$density, type = "bar", name = "Sample Histogram") |>
add_lines(x = 0:6, y = dpois(0:6, 1), name = "Theoretical Poisson mass") |>
add_lines(x = density(pois_sample)$x, y = density(pois_sample)$y,
name = "Gaussian KDE (for contrast)") |>
layout(bargap = 0.1, title = "Histogram (Simulated Poisson Data)",
legend = list(orientation = "h"))Note the equal-width breaks, instead of using hard-coded
break-points, e.g., c(-0.5, 0.5, 1.5, 2.5, 3.5, 6.5), which
may yield a final bin three units wide. A bar of that (wider) bin’s
density drawn at the same width as the others misrepresents its
mass.
16.2 The Data Modeler
A recurring task is choosing a parametric family for an observed feature. Our example is crystallographic spectral data from the Ivanova Lab, nine samples, each a length spectrum with its own sampling rate.
crystal <- dspa_read(
"https://umich.instructure.com/files/11653615/download?download_frd=1",
"crystallography_length.csv", header = TRUE)
col_names <- colnames(crystal); col_names#> [1] "AC1338" "AC1432" "AC1593" "AC1679" "AC1860" "AC1874" "AC1881" "AC1903"
#> [9] "Rec"#> AC1338 AC1432 AC1593 AC1679 AC1860 AC1874 AC1881 AC1903 Rec
#> 453 358 477 399 337 489 591 319 653crystal_long <- crystal |>
pivot_longer(everything(), names_to = "sample", values_to = "value") |>
filter(!is.na(value))
ggplot(crystal_long, aes(value, colour = sample)) +
geom_density(bw = 15, linewidth = 0.8) +
labs(title = "Crystallography sample densities",
subtitle = "Common bandwidth (15) so the curves are comparable",
x = "Spectral length", y = "Density", colour = "Sample") +
theme_dspa()
# --- Interactive equivalent -----------------------------------------------
cc <- crystal[complete.cases(crystal), ]
dl <- apply(cc, 2, density, kernel = "gaussian", bw = 15)
df <- data.frame(x = unlist(lapply(dl, `[[`, "x")),
y = unlist(lapply(dl, `[[`, "y")),
sample = rep(names(dl), each = length(dl[[1]]$x)))
plot_ly(df, x = ~x, y = ~y, color = ~sample, type = "scatter", mode = "lines") |>
layout(title = "Crystallography Sample Densities",
legend = list(title = list(text = "<b> Samples </b>")),
xaxis = list(title = "X"), yaxis = list(title = "Density"))16.3 Fitting univariate models by maximum likelihood
We fit three positive-support families to each sample. Their densities, means, and modes:
| Family | Density \(f(x)\), \(x>0\) | Mean | Mode |
|---|---|---|---|
| Weibull\((k,\lambda)\) | \(\dfrac{k}{\lambda}\left(\dfrac{x}{\lambda}\right)^{k-1}e^{-(x/\lambda)^k}\) | \(\lambda\,\Gamma\!\left(1+\tfrac1k\right)\) | \(\lambda\left(\tfrac{k-1}{k}\right)^{1/k}\) if \(k>1\), else \(0\) |
| Gamma\((\alpha,\beta)\) | \(\dfrac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x}\) | \(\alpha/\beta\) | \((\alpha-1)/\beta\) if \(\alpha\ge 1\), else no interior mode |
| Log-normal\((\mu,\sigma)\) | \(\dfrac{1}{x\sigma\sqrt{2\pi}}\exp\!\left(-\dfrac{(\ln x-\mu)^2}{2\sigma^2}\right)\) | \(e^{\mu+\sigma^2/2}\) | \(e^{\mu-\sigma^2}\) |
In this example, the domain restrictions are critically important.
library(fitdistrplus)
fit_W <- fit_G <- fit_LN <- vector("list", col_num)
for (i in seq_len(col_num)) {
v <- crystal[[i]]; v <- v[!is.na(v) & v > 0]
fit_W[[i]] <- fitdist(v, "weibull")
fit_G[[i]] <- fitdist(v, "gamma")
fit_LN[[i]] <- fitdist(v, "lnorm")
}
weibull_mode <- function(k, lambda) if (k > 1) lambda * ((k - 1) / k)^(1 / k) else 0
gamma_mode <- function(a, b) if (a >= 1) (a - 1) / b else NA_real_
model_stats <- do.call(rbind, lapply(seq_len(col_num), \(i) {
kW <- fit_W[[i]]$estimate[["shape"]]; lW <- fit_W[[i]]$estimate[["scale"]]
aG <- fit_G[[i]]$estimate[["shape"]]; bG <- fit_G[[i]]$estimate[["rate"]]
mL <- fit_LN[[i]]$estimate[["meanlog"]]; sL <- fit_LN[[i]]$estimate[["sdlog"]]
data.frame(
sample = col_names[i],
W_shape = kW, W_scale = lW,
W_mean = lW * gamma(1 + 1 / kW),
W_mode = weibull_mode(kW, lW),
W_sd = lW * sqrt(gamma(1 + 2/kW) - gamma(1 + 1/kW)^2),
G_mean = aG / bG, G_mode = gamma_mode(aG, bG), G_sd = sqrt(aG) / bG,
LN_mean = exp(mL + sL^2 / 2), LN_mode = exp(mL - sL^2),
LN_sd = sqrt((exp(sL^2) - 1) * exp(2 * mL + sL^2)),
AIC_W = fit_W[[i]]$aic, AIC_G = fit_G[[i]]$aic, AIC_LN = fit_LN[[i]]$aic)
}))
model_stats |>
dplyr::select(sample, W_mean, W_mode, G_mean, G_mode, LN_mean, LN_mode) |>
mutate(across(where(is.numeric), \(z) round(z, 2)))Mean and mode differ substantially for every sample, which is the
point. On a right-skewed density they are different questions,
and a legend that says m= without saying which is
uninformative.
grid_for <- function(i, n = 400) {
v <- crystal[[i]]; v <- v[!is.na(v) & v > 0]
z <- seq(min(v), max(v), length.out = n)
kW <- fit_W[[i]]$estimate; aG <- fit_G[[i]]$estimate; mL <- fit_LN[[i]]$estimate
bind_rows(
data.frame(x = z, dens = dweibull(z, kW[["shape"]], kW[["scale"]]), model = "Weibull"),
data.frame(x = z, dens = dgamma(z, aG[["shape"]], aG[["rate"]]), model = "Gamma"),
data.frame(x = z, dens = dlnorm(z, mL[["meanlog"]], mL[["sdlog"]]), model = "Log-normal")
) |> mutate(sample = col_names[i])
}
dens_all <- bind_rows(lapply(seq_len(col_num), grid_for))
ggplot() +
geom_histogram(data = crystal_long, aes(value, after_stat(density)),
bins = 25, fill = "grey87", colour = "white") +
geom_line(data = dens_all, aes(x, dens, colour = model), linewidth = 0.8) +
facet_wrap(~ sample, scales = "free", ncol = 3) +
scale_colour_manual(values = c(Weibull = "firebrick", Gamma = "forestgreen",
`Log-normal` = "steelblue")) +
labs(title = "Maximum-likelihood fits: Weibull, Gamma, Log-normal",
x = "Spectral length", y = "Density", colour = NULL) +
theme_dspa(10)
# --- Interactive equivalent: one panel per sample, subplotted --------------
pl_list <- lapply(seq_len(col_num), function(i) {
v <- crystal[[i]]; v <- v[!is.na(v) & v > 0]
z <- seq(min(v), max(v), length.out = 400)
kW <- fit_W[[i]]$estimate; aG <- fit_G[[i]]$estimate; mL <- fit_LN[[i]]$estimate
lg <- c(sprintf("Weibull(k=%.2f, lambda=%.1f)", kW[["shape"]], kW[["scale"]]),
sprintf("Gamma(alpha=%.2f, beta=%.3f)", aG[["shape"]], aG[["rate"]]),
sprintf("LogNormal(mu=%.2f, sigma=%.2f)", mL[["meanlog"]], mL[["sdlog"]]))
plot_ly(x = ~v, type = "histogram", histnorm = "probability density",
name = col_names[i], nbinsx = 25, opacity = 0.6, showlegend = FALSE,
marker = list(color = "lightgray",
line = list(color = "darkgray", width = 1))) |>
add_trace(x = z, y = dweibull(z, kW[["shape"]], kW[["scale"]]),
type = "scatter", mode = "lines", name = lg[1],
line = list(color = "red", width = 2)) |>
add_trace(x = z, y = dgamma(z, aG[["shape"]], aG[["rate"]]),
type = "scatter", mode = "lines", name = lg[2],
line = list(color = "green", width = 2)) |>
add_trace(x = z, y = dlnorm(z, mL[["meanlog"]], mL[["sdlog"]]),
type = "scatter", mode = "lines", name = lg[3],
line = list(color = "blue", width = 2)) |>
layout(bargap = 0.1,
xaxis = list(title = col_names[i]), yaxis = list(title = "Density"))
})
plotly::subplot(pl_list, nrows = 3, titleX = TRUE, titleY = TRUE) |>
layout(title = "Distribution models of crystallography data (interactive)")
# fitdistrplus also supplies four diagnostic panels:
# denscomp(); cdfcomp(); qqcomp(); ppcomp()
# NOTE: do NOT call windows() to size these - that device is Windows-only and
# is invisible to knitr on every platform. Use chunk fig.width / fig.height.ecdf_df <- crystal_long |>
arrange(sample, value) |>
mutate(F_emp = seq_along(value) / dplyr::n(), .by = sample)
cdf_all <- bind_rows(lapply(seq_len(col_num), \(i) {
v <- crystal[[i]]; v <- v[!is.na(v) & v > 0]
z <- seq(min(v), max(v), length.out = 400)
kW <- fit_W[[i]]$estimate; aG <- fit_G[[i]]$estimate; mL <- fit_LN[[i]]$estimate
bind_rows(
data.frame(x = z, F = pweibull(z, kW[["shape"]], kW[["scale"]]), model = "Weibull"),
data.frame(x = z, F = pgamma(z, aG[["shape"]], aG[["rate"]]), model = "Gamma"),
data.frame(x = z, F = plnorm(z, mL[["meanlog"]], mL[["sdlog"]]), model = "Log-normal")
) |> mutate(sample = col_names[i])
}))
ggplot() +
geom_step(data = filter(ecdf_df, sample %in% col_names[1:3]),
aes(value, F_emp), colour = "grey25", linewidth = 0.5) +
geom_line(data = filter(cdf_all, sample %in% col_names[1:3]),
aes(x, F, colour = model), linewidth = 0.8) +
facet_wrap(~ sample, ncol = 3) +
scale_colour_manual(values = c(Weibull = "firebrick", Gamma = "forestgreen",
`Log-normal` = "steelblue")) +
labs(title = "Empirical CDF (step) against fitted model CDFs",
subtitle = "The vertical gap is the Kolmogorov-Smirnov statistic D_n",
x = "Spectral length", y = "F(x)", colour = NULL) +
theme_dspa(10)
16.4 Goodness of fit when parameters are estimated
The Kolmogorov–Smirnov statistic is the largest vertical gap above:
\[D_n=\sup_x\left|F_n(x)-F(x;\theta)\right| .\]
Its classical null distribution is distribution-free only when \(\theta\) is specified in advance. Here \(\hat\theta\) was estimated from the same sample, so the fitted CDF is pulled toward the empirical CDF, \(D_n\) is too small, and the nominal p-value is anti-conservative, it will fail to reject models it should reject.
The fix is a parametric bootstrap (the Lilliefors construction):
- Fit \(\hat\theta\) to the data; compute \(D_n\).
- For \(b=1,\dots,B\): simulate \(n\) values from \(F(\cdot;\hat\theta)\), refit to get \(\hat\theta^{*}_b\), compute \(D^{*}_b\).
- \(\displaystyle \hat p=\frac{1+\#\{b:D^{*}_b\ge D_n\}}{B+1}\).
The \(+1\) in numerator and denominator keeps \(\hat p>0\) and makes the test valid at finite \(B\).
ks_bootstrap <- function(v, dist = c("weibull", "gamma", "lnorm"), B = 200, seed = 1) {
dist <- match.arg(dist)
set.seed(seed)
v <- v[!is.na(v) & v > 0]; n <- length(v)
pfun <- switch(dist, weibull = pweibull, gamma = pgamma, lnorm = plnorm)
rfun <- switch(dist, weibull = rweibull, gamma = rgamma, lnorm = rlnorm)
fit0 <- fitdist(v, dist)
D0 <- suppressWarnings(
ks.test(v, pfun, fit0$estimate[[1]], fit0$estimate[[2]])$statistic)
Dstar <- replicate(B, {
vb <- rfun(n, fit0$estimate[[1]], fit0$estimate[[2]])
fb <- try(fitdist(vb, dist), silent = TRUE)
if (inherits(fb, "try-error")) return(NA_real_)
suppressWarnings(
ks.test(vb, pfun, fb$estimate[[1]], fb$estimate[[2]])$statistic)
})
Dstar <- Dstar[!is.na(Dstar)]
naive <- suppressWarnings(
ks.test(v, pfun, fit0$estimate[[1]], fit0$estimate[[2]])$p.value)
c(D = unname(D0),
p_naive = unname(naive),
p_bootstrap = (1 + sum(Dstar >= D0)) / (length(Dstar) + 1))
}
gof <- do.call(rbind, lapply(seq_len(3), \(i) {
v <- crystal[[i]]
data.frame(sample = col_names[i],
model = c("weibull", "gamma", "lnorm"),
t(sapply(c("weibull", "gamma", "lnorm"), \(d) ks_bootstrap(v, d, B = 200))))
}))
gof |> mutate(across(where(is.numeric), \(z) signif(z, 4)))Compare the two p-value columns. The naive column is systematically larger — that is the anti-conservatism, quantified. Whenever a naive KS p-value is comfortably above 0.05, ask whether it survives the bootstrap.
For comparing models, use information criteria, not p-values. A p-value asks whether one model is refutable; AIC and BIC ask which fits best per parameter spent:
\[\mathrm{AIC}=-2\hat\ell+2k,\qquad \mathrm{BIC}=-2\hat\ell+k\log n .\]
model_stats |>
dplyr::select(sample, AIC_W, AIC_G, AIC_LN) |>
mutate(best = c("Weibull", "Gamma", "Log-normal")[
max.col(-as.matrix(across(c(AIC_W, AIC_G, AIC_LN))))],
across(where(is.numeric), \(z) round(z, 1)))A difference of \(\Delta\mathrm{AIC}<2\) is not evidence
of a preference; \(>10\) is strong.
fitdistrplus::gofstat() also reports Cramér–von Mises and
Anderson–Darling statistics, the latter being more sensitive in the
tails — usually where distributional choice actually matters.
16.5 Mixture distributions
A single family often fails because the sample blends several subpopulations. A finite mixture with \(K\) components is
\[f(x\mid\Theta)=\sum_{k=1}^{K}\lambda_k\,f_k(x\mid\theta_k),\qquad \lambda_k\ge 0,\;\sum_k\lambda_k=1 .\]
Direct maximization is hard because the log of a sum does not separate. EM solves it by introducing latent labels \(Z_i\in\{1,\dots,K\}\).
E-step, responsibilities. The posterior probability that observation \(i\) came from component \(k\):
\[\gamma_{ik}^{(t)}=\frac{\lambda_k^{(t)}f_k\!\left(x_i\mid\theta_k^{(t)}\right)}{\sum_{j=1}^{K}\lambda_j^{(t)}f_j\!\left(x_i\mid\theta_j^{(t)}\right)} .\]
M-step. Maximize the expected complete-data log-likelihood \(\sum_i\sum_k\gamma_{ik}\big[\log\lambda_k+\log f_k(x_i\mid\theta_k)\big]\). The weights always update as
\[\lambda_k^{(t+1)}=\frac1n\sum_{i=1}^{n}\gamma_{ik}^{(t)},\]
and for Gaussian components the parameters have closed forms:
\[\mu_k^{(t+1)}=\frac{\sum_i\gamma_{ik}x_i}{\sum_i\gamma_{ik}},\qquad \big(\sigma_k^{2}\big)^{(t+1)}=\frac{\sum_i\gamma_{ik}\big(x_i-\mu_k^{(t+1)}\big)^2}{\sum_i\gamma_{ik}} .\]
For Weibull and Gamma components there is no closed form, so the M-step is itself iterative, a generalized EM.
Three properties you must know before trusting a mixture fit.
Complexity. \(O(nK)\) density evaluations per iteration for the E-step, so \(O(TnK)\) overall for Gaussian components; more when the M-step is iterative.
Local maxima. The mixture likelihood is multimodal. Different starting values give different answers. Always restart from several initializations and keep the best log-likelihood.
Unboundedness. For Gaussian mixtures with free variances the likelihood is unbounded: place \(\mu_k\) exactly on a data point and let \(\sigma_k\to 0\), and the likelihood diverges. The MLE in the strict sense does not exist; what EM finds is a well-behaved local maximum. Variance floors or priors are the standard guards.
Label switching. Component labels are not identifiable, permuting them gives the same mixture. Never interpret “component 1” across restarts without imposing an ordering constraint.
library(mixtools)
K <- 3
mix_N <- mix_G <- mix_W <- vector("list", col_num)
invisible(capture.output(
for (i in seq_len(col_num)) {
v <- crystal[[i]]; v <- v[!is.na(v) & v > 0]
# Multiple restarts: keep the best log-likelihood (guards against local maxima)
best <- NULL
for (s in 1:5) {
set.seed(100 + s)
cand <- try(normalmixEM(v, k = K, maxit = 2000, epsilon = 1e-6),
silent = TRUE)
if (!inherits(cand, "try-error") &&
(is.null(best) || cand$loglik > best$loglik)) best <- cand
}
mix_N[[i]] <- best
set.seed(202); mix_G[[i]] <- try(gammamixEM(v, k = K), silent = TRUE)
set.seed(203); mix_W[[i]] <- try(weibullRMM_SEM(v, k = K, verb = FALSE),
silent = TRUE)
}))
c(components = K,
loglik_sample1 = round(mix_N[[1]]$loglik, 2),
weights_sample1 = paste(round(sort(mix_N[[1]]$lambda, decreasing = TRUE), 3),
collapse = " / "))#> components loglik_sample1 weights_sample1
#> "3" "-2295.33" "0.579 / 0.404 / 0.018"mix_curves <- function(i) {
m <- mix_N[[i]]; v <- crystal[[i]]; v <- v[!is.na(v) & v > 0]
z <- seq(min(v), max(v), length.out = 500)
comp <- bind_rows(lapply(seq_len(K), \(k)
data.frame(x = z, dens = m$lambda[k] * dnorm(z, m$mu[k], m$sigma[k]),
curve = paste0("Component ", k))))
total <- data.frame(x = z, dens = rowSums(sapply(seq_len(K), \(k)
m$lambda[k] * dnorm(z, m$mu[k], m$sigma[k]))), curve = "Mixture")
bind_rows(comp, total) |> mutate(sample = col_names[i])
}
mx <- bind_rows(lapply(1:2, mix_curves))
ggplot() +
geom_histogram(data = filter(crystal_long, sample %in% col_names[1:2]),
aes(value, after_stat(density)), bins = 30,
fill = "pink", colour = "white") +
geom_line(data = mx, aes(x, dens, colour = curve,
linewidth = curve == "Mixture")) +
scale_linewidth_manual(values = c(`FALSE` = 0.7, `TRUE` = 1.3), guide = "none") +
scale_colour_manual(values = c("Component 1" = "#3B4CC0", "Component 2" = "#8C8CC0",
"Component 3" = "#C03B3B", "Mixture" = "black")) +
facet_wrap(~ sample, scales = "free", ncol = 2) +
labs(title = sprintf("Mixture of %d normal components, fitted by EM", K),
subtitle = "Coloured: weighted components. Black: their sum",
x = "Spectral length", y = "Density", colour = NULL) +
theme_dspa(11)
# --- Interactive equivalent -----------------------------------------------
i <- 1
m <- mix_N[[i]]; v <- crystal[[i]]; v <- v[!is.na(v) & v > 0]
z <- seq(min(v), max(v), length.out = 500)
p <- plot_ly(x = ~v, type = "histogram", histnorm = "probability density",
name = "Data", opacity = 0.55, nbinsx = 30)
for (k in seq_len(K))
p <- add_lines(p, x = z, y = m$lambda[k] * dnorm(z, m$mu[k], m$sigma[k]),
name = sprintf("Component %d (lambda=%.2f)", k, m$lambda[k]))
p |> add_lines(x = z,
y = rowSums(sapply(seq_len(K), \(k)
m$lambda[k] * dnorm(z, m$mu[k], m$sigma[k]))),
name = "Mixture", line = list(color = "black", width = 4)) |>
layout(bargap = 0.1,
title = sprintf("Mixture of %d Normal Models of %s", K, col_names[i]),
xaxis = list(title = "Intensities"), yaxis = list(title = "Density"))mixture_summary <- function(i) {
m <- mix_N[[i]]
ord <- order(m$mu) # order by mean: guards label switching
data.frame(sample = col_names[i], component = seq_len(K),
weight = round(m$lambda[ord], 3),
mean = round(m$mu[ord], 2),
sd = round(m$sigma[ord], 2),
loglik = round(m$loglik, 1))
}
bind_rows(lapply(1:2, mixture_summary))Components are sorted by mean before reporting. Without an ordering
constraint, “component 1” means nothing across samples or restarts, the
label-switching problem, made harmless by one order()
call.
17 Two- and three-dimensional density surfaces
Everything in this section is interactive and evaluated. A surface must be rotated to be read; a static projection of a 3-D object discards precisely the information the third dimension was added to convey.
17.1 2-D kernel density estimation
The bivariate kernel density estimator generalizes §2.4.3 with a bandwidth matrix \(H\):
\[\hat f_H(\mathbf{x})=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{\sqrt{\det H}}\, K\!\left(H^{-1/2}(\mathbf{x}-\mathbf{x}_i)\right),\qquad \mathbf{x}\in\mathbb R^2 .\]
MASS::kde2d() uses a diagonal \(H=\mathrm{diag}(h_1^2,h_2^2)\) and
evaluates on an \(n\times n\) grid,
returning x (grid abscissae), y (grid
ordinates), and a matrix z with \(z_{jk}=\hat f_H(x_j,y_k)\).
Interpreting
z. Mathematically,zis the estimated density evaluated on the grid, a nonlinear functional of the data, not of the grid vectors.
Computational cost is \(O(n\,g^2)\)
for \(n\) observations on a \(g\times g\) grid, which is why
kde2d() defaults to \(g=25\); the FFT-based
KernSmooth::bkde2D() reduces this to \(O(n+g^2\log g)\) by binning first.
#> List of 3
#> $ x: num [1:50] 0.833 0.928 1.022 1.116 1.21 ...
#> $ y: num [1:50] 43 44.3 45.7 47 48.3 ...
#> $ z: num [1:50, 1:50] 9.07e-13 1.81e-12 3.43e-12 6.11e-12 1.03e-11 ...#> [,1] [,2] [,3] [,4]
#> [1,] 9.068691e-13 4.238943e-12 1.839285e-11 7.415672e-11
#> [2,] 1.814923e-12 8.473636e-12 3.671290e-11 1.477410e-10
#> [3,] 3.428664e-12 1.599235e-11 6.920273e-11 2.780463e-10
#> [4,] 6.114498e-12 2.849475e-11 1.231748e-10 4.942437e-10with(kd, plot_ly(x = x, y = y, z = z, type = "surface")) |>
layout(title = "Old Faithful: 2-D kernel density of eruption duration vs. waiting time",
scene = list(xaxis = list(title = "Duration (min)"),
yaxis = list(title = "Waiting (min)"),
zaxis = list(title = "Density")))
Rotate the surface: the two peaks are the point. Old Faithful’s eruptions form two clusters, short duration / short wait and long duration / long wait — which a marginal histogram of either variable alone would blur into one broad mode. This is a mixture (§2.14.4) made visible.
17.2 Deterministic 3-D surfaces
R’s built-in volcano records elevation of Maunga Whau on
an \(87\times 61\) grid.
#> [1] 87 61plot_ly(z = volcano, type = "surface") |>
layout(title = "Maunga Whau volcano, topographic surface",
scene = list(xaxis = list(title = "Grid (east)"),
yaxis = list(title = "Grid (north)"),
zaxis = list(title = "Elevation (m)")))
Here z really is all you need: plot_ly()
supplies integer grid coordinates for x and y
when they are omitted.
17.3 Confidence surfaces on an image
A grayscale image is a 2-D array of intensities, that is, a surface. Stacking a smoothed version with offset copies gives a visual analogue of a confidence band in two dimensions.
library(jpeg)
img_url <- "https://umich.instructure.com/files/1627149/download?download_frd=1"
img_file <- file.path(tempdir(), "MRI_ImageHematoma.jpg")
ok <- tryCatch({ download.file(img_url, img_file, mode = "wb", quiet = TRUE); TRUE },
error = function(e) FALSE)
if (ok) {
img <- readJPEG(img_file)[, , 1] # first channel of RGB -> univariate 2-D array
dim(img)
}#> [1] 256 256img_s <- as.matrix(spatstat.explore::blur(spatstat.geom::as.im(img), sigma = 10))
plot_ly(z = img, type = "surface", showscale = FALSE, name = "Observed") |>
add_trace(z = img_s + 1, type = "surface", showscale = FALSE, opacity = 0.95,
name = "Upper band") |>
add_trace(z = img_s - 1, type = "surface", showscale = FALSE, opacity = 0.95,
name = "Lower band") |>
layout(title = "MRI intensity surface with smoothed upper/lower envelopes",
scene = list(zaxis = list(title = "Intensity")))spatstat was split into spatstat.geom,
spatstat.explore, and siblings; blur() now
lives in spatstat.explore and as.im() in
spatstat.geom.
17.4 3-D and 4-D neuroimaging
Many biomedical data are intrinsically multidimensional. The body is a 3-D solid (sMRI) that changes over time (fMRI), so a functional acquisition is a 4-D hypervolume \(I(x,y,z,t)\).
The SOCR BrainViewer demonstrates browser-based navigation of 2-D cross-sections, volume rendering, and embedded 1-D and 2-D models in a shared 3-D scene.
# --- Volumetric rendering with brainR + rgl ---------------------------------
# Not evaluated during knitting: rgl's WebGL output fails to load on some
# browser/OS combinations. Runs well in an interactive session.
library(brainR)
brainURL <- "https://socr.umich.edu/HTML5/BrainViewer/data/TestBrain.nii.gz"
brainFile <- file.path(tempdir(), "TestBrain.nii.gz")
download.file(brainURL, brainFile, quiet = TRUE)
brainVolume <- oro.nifti::readNIfTI(brainFile, reorient = FALSE)
dims <- dim(brainVolume); dims
# Lower isosurface levels give smoother surfaces; see ?contour3d
contour3d(brainVolume, level = 20, alpha = 0.1, draw = TRUE)
# Multiple levels reveal nested shells (e.g. hyper-intense white matter)
contour3d(brainVolume, level = c(10, 120), alpha = c(0.3, 0.5),
add = TRUE, color = c("yellow", "red"))
text3d(x = dims[1] / 2, y = dims[2] / 2, z = dims[3] * 0.98, text = "Top")
text3d(x = dims[1] * 0.98, y = dims[2] / 2, z = dims[3] / 2, text = "Right")Additional volumes in NIfTI format: sMRI (3-D structural), fMRI (4-D functional), PET (3-D perfusion).
library(oro.nifti)
fmri_url <- "https://socr.umich.edu/HTML5/BrainViewer/data/fMRI_FilteredData_4D.nii.gz"
fmri_file <- file.path(tempdir(), "fMRI_FilteredData_4D.nii.gz")
download.file(fmri_url, fmri_file, quiet = TRUE)
fmri <- readNIfTI(fmri_file, reorient = FALSE)
dim(fmri) # 64 x 64 x 21 x 180 voxels; 4mm x 4mm x 6mm x 3s#> [1] 64 64 21 180# Orthographic view: axial slice through the thalamus centers the crosshair
library(oro.nifti)
orthographic(fmri, xyz = c(34, 29, 10), zlim = range(fmri) * 0.9)
vox <- as.vector(fmri[, , , 1])
ggplot(data.frame(v = vox), aes(v)) +
geom_histogram(bins = 80, fill = "steelblue", colour = NA) +
scale_y_log10() +
labs(title = "fMRI voxel intensity distribution (volume 1)",
subtitle = "Log-scaled count axis; the huge spike near zero is background outside the head",
x = "Intensity", y = "Count (log scale)") +
theme_dspa()
# --- Interactive equivalent -----------------------------------------------
h <- hist(fmri, plot = FALSE)
plot_ly(x = h$mids, y = h$density, type = "bar") |>
layout(bargap = 0.1, title = "fMRI Histogram")
stat_fmri <- ifelse(fmri > 15000, fmri, NA)
h2 <- hist(stat_fmri, plot = FALSE)
plot_ly(x = h2$mids, y = h2$density, type = "bar") |>
layout(bargap = 0.1, title = "fMRI Histogram (high intensities)")The time course of a single voxel is a 1-D series embedded in the 4-D volume. Because the reader will want to zoom into individual acquisitions, this one is interactive:
tt <- seq_len(dim(fmri)[4])
ts_ <- fmri[30, 30, 10, ]
lo <- loess(ts_ ~ tt, span = 0.25) # genuine LOESS
med <- stats::smooth(ts_) # Tukey running median
ks <- ksmooth(tt, ts_, kernel = "normal", bandwidth = 5)
plot_ly(x = tt, y = ts_, type = "scatter", mode = "lines",
name = "Raw fMRI", line = list(width = 1)) |>
add_lines(x = tt, y = predict(lo), name = "LOESS (span 0.25)",
line = list(width = 3)) |>
add_lines(x = tt, y = as.numeric(med), name = "Tukey running median") |>
add_lines(x = ks$x, y = ks$y, name = "Gaussian kernel smoother (bw 5)") |>
layout(title = "Time series of voxel (x=30, y=30, z=10)",
xaxis = list(title = "Acquisition (TR = 3 s)"),
yaxis = list(title = "BOLD intensity"),
legend = list(orientation = "h"))This examples shows three different smoothers. Chapter 12 develops time-series and longitudinal modeling properly. DSPA Appendix 3 covers parametric and implicit manifold visualization.
18 Parsing and visualizing web data
Much public health and geoscience data lives in HTML tables.
rvest::html_table() turns one into a data frame.
quake_tbl <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_Dinov_021708_Earthquakes") |>
html_nodes("table") |> _[[2]] |> html_table()
c(rows = nrow(quake_tbl), cols = ncol(quake_tbl))#> rows cols
#> 168 13#> [1] 5.00 7.39Scraping is fragile by construction.
[[2]] is a positional dependency on page layout: an editor
inserting a table above shifts the index and your code silently loads
the wrong data. Prefer a CSS selector or an id, validate
the column names after loading, and cache the result (§1.6.7).
18.1 From points to a gridded surface
To render magnitudes as a surface we need a matrix on a regular grid. Rounding coordinates to integer degrees and aggregating is transparent and correct:
# Earlier editions built this with sparseMatrix() using `180 - x` in one place
# and `180 + x` in another, two different matrices for the same data, and
# passed dimnames = list(levels(x), levels(y)) on numeric vectors, which is
# always list(NULL, NULL). xtabs() is clearer and keeps real dimnames.
quake_grid <- quake |>
transmute(lon = as.integer(round(Longitude)),
lat = as.integer(round(Latitude)),
mag = as.numeric(Mag))
mag_sum <- xtabs(mag ~ lon + lat, data = quake_grid) # summed magnitude
mag_max <- tapply(quake_grid$mag,
list(quake_grid$lon, quake_grid$lat), max, default = 0)
dim(mag_sum)#> [1] 12 11#> lat
#> lon 33 34 35 36 37
#> -128 0 0 0 0 0
#> -126 0 0 0 0 0
#> -125 0 0 0 0 0
#> -124 0 0 0 5 0as.data.frame(mag_sum) |>
mutate(lon = as.numeric(as.character(lon)),
lat = as.numeric(as.character(lat))) |>
ggplot(aes(lon, lat, fill = Freq)) +
geom_raster() +
scale_fill_viridis_c(option = "inferno", name = "Summed\nmagnitude") +
coord_fixed() +
labs(title = "California seismicity, 1969-2007",
subtitle = "Total magnitude per 1-degree cell",
x = "Longitude", y = "Latitude") +
theme_dspa()
z_mat <- t(as.matrix(mag_sum)) # rows = latitude, cols = longitude
plot_ly(x = as.numeric(colnames(z_mat)),
y = as.numeric(rownames(z_mat)),
z = z_mat, type = "surface") |>
layout(title = "Aggregate earthquake magnitude surface",
scene = list(xaxis = list(title = "Longitude"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "Summed magnitude")))plot_ly(quake, x = ~Longitude, y = ~Latitude, z = ~Mag,
type = "scatter3d", mode = "markers",
color = ~Depth, size = ~Mag,
text = ~paste0("Mag ", Mag, " | depth ", Depth, " km | type ", Magt),
marker = list(sizemode = "diameter", opacity = 0.75)) |>
layout(title = "Earthquake epicenters in 3-D: location and magnitude",
scene = list(xaxis = list(title = "Longitude"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "Magnitude")))kq <- with(quake, MASS::kde2d(Longitude, Latitude, n = 50))
with(kq, plot_ly(x = x, y = y, z = z, type = "surface")) |>
layout(title = "Kernel density of epicenter locations",
scene = list(xaxis = list(title = "Longitude"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "Density")))19 Computational complexity summary
\(n\) = observations, \(d\) = features, \(T\) = iterations, \(m\) = imputations, \(K\) = mixture components, \(P\) = distinct missingness patterns, \(g\) = grid size.
| Method | Time | Memory | Practical bottleneck |
|---|---|---|---|
| Histogram | \(O(n)\) | \(O(k)\) | Bin-rule choice, not compute |
| Quantiles, five-number summary | \(O(n\log n)\) | \(O(n)\) | Sorting |
| KDE, direct | \(O(n\,g)\) | \(O(g)\) | Bandwidth selection |
| KDE, binned/FFT | \(O(n+g\log g)\) | \(O(g)\) | , |
2-D KDE (kde2d) |
\(O(n\,g^{2})\) | \(O(g^{2})\) | Grid resolution |
| \(\chi^2\) contingency test | \(O(n+rc)\) | \(O(rc)\) | Sparse expected counts |
| Correlation matrix | \(O(n d^{2})\) | \(O(d^{2})\) | \(d^2\) memory |
| Mean imputation | \(O(nd)\) | \(O(1)\) | Bias, not cost |
| EM imputation, per-row | \(O(T\,n\,d_{obs}^{3})\) | \(O(d^{2})\) | Repeated identical solves |
| EM imputation, pattern-grouped | \(O(T(P d^{3}+n d^{2}))\) | \(O(d^{2})\) | \(P\ll n\) makes this the win |
| MICE (chained equations) | \(O(T\,m\,d\cdot c_{\text{reg}})\) | \(O(mnd)\) | \(m\) full passes |
| Amelia (bootstrap EM) | \(O(m\,T\,n\,d^{2})\) | \(O(nd)\) | Normality assumption |
| SMOTE (k-d tree) | \(O(n_{min}\log n_{min}\cdot d)\) | \(O(n_{min}d)\) | Degrades in high \(d\) |
| Hierarchical clustering | \(O(n^{2}\log n)\) | \(\mathbf{O(n^{2})}\) | Memory, at \(n\gtrsim 10^4\) |
| LOESS (interpolated) | \(O(n\log n)\) | \(O(n)\) | Span selection |
MLE fit (fitdist) |
\(O(T n)\) | \(O(n)\) | Optimizer starting values |
| Mixture EM | \(O(T\,n\,K)\) | \(O(nK)\) | Local maxima; restart |
| Parametric-bootstrap GOF | \(O(B(n\log n+c_{\text{fit}}))\) | \(O(n)\) | \(B\) refits |
Two entries deserve emphasis. Pattern-grouped EM
turns \(O(n d^3)\) into \(O(P d^3 + n d^2)\), the same output, a
different asymptotic class. Hierarchical clustering’s \(O(n^2)\) memory is what stops
hclust() at around \(10^4\) objects, long before time becomes
the issue; beyond that, use fastcluster, mini-batch
k-means, or HDBSCAN.
20 Common pitfalls
| # | Pitfall | Consequence | Fix |
|---|---|---|---|
| 1 | sep = ", " in
read.csv()/write.csv() |
Every field gains a leading space; numerics become character | sep = "," |
| 2 | prop.table(v) on a raw vector |
Returns non-proportions without complaint | prop.table(table(v)) |
| 3 | Reading the \(\chi^2\) cell contribution as a probability | Nonsense inference | Contribution \(\ne\) residual \(\ne\) p-value |
| 4 | Concluding “no difference” from a large p-value | Accepting the null | Report effect size; use an equivalence test |
| 5 | Mean imputation, then computing SEs | Variance deflated by \((1-p)\); CIs too narrow | Multiple imputation |
| 6 | Sigma <- var(imputed) in an EM M-step |
Omits the conditional covariance; not the MLE | Add \(C_i\) |
| 7 | solve(S, tol = 1e-40) |
Silences the warning, not the ill-conditioning | Check rcond(); ridge |
| 8 | SMOTE before the train/test split | Optimistically biased performance | Resample inside the fold |
| 9 | Wilcoxon test comparing original vs. resampled data | Samples are not independent | Compare descriptively |
| 10 | Naive KS p-value with estimated parameters | Anti-conservative; fails to reject | Parametric bootstrap |
| 11 | \(\sqrt n\) bins | Wrong growth rate | Scott or Freedman–Diaconis |
| 12 | dpois() on a continuous grid |
Zero almost everywhere, plus warnings | Evaluate at integers |
| 13 | Unequal histogram breaks compared to a density | Misrepresents the wide bin’s mass | Equal-width breaks |
| 14 | Gaussian KDE on count data | Density where no mass exists; negative support | Empirical mass function |
| 15 | attach() then reassigning the object |
Bare names silently use stale data | with(), or $ |
| 16 | windows() for figure sizing |
Errors off Windows; invisible to knitr | Chunk fig.width/fig.height |
| 17 | 26-slice pie chart | Unreadable (angle/area encoding) | Bar or dot plot |
| 18 | Raw p-value heatmap, no colour bar | Signal invisible, values unreadable | \(-\log_{10}p\), BH-FDR, show the bar |
| 19 | cor() without stating the margin |
Correlates the wrong entities | Name rows vs. columns |
| 20 | Unlabelled error bars | SD and SE differ by \(\sqrt n\) | Say which, in the caption |
| 21 | Hard-coded 1:3331 |
Silent truncation when data change | seq_len(nrow(x)) |
| 22 | Interpreting mixture “component 1” across restarts | Label switching | Order components by a parameter |
21 Practice problems
21.1 Problem 1, Bin-width rules from scratch
Implement Scott’s and Freedman–Diaconis’ rules without
nclass.*(), apply them to a strongly right-skewed sample,
and explain which is preferable and why.
Solution
set.seed(4)
v <- rlnorm(2000, meanlog = 1, sdlog = 1.2)
h_scott <- 3.49 * sd(v) * length(v)^(-1/3)
h_fd <- 2 * IQR(v) * length(v)^(-1/3)
c(h_scott = round(h_scott, 3), h_fd = round(h_fd, 3),
bins_scott = ceiling(diff(range(v)) / h_scott),
bins_fd = ceiling(diff(range(v)) / h_fd),
check_scott = nclass.scott(v), check_fd = nclass.FD(v))#> h_scott h_fd bins_scott bins_fd check_scott check_fd
#> 2.337 0.773 53.000 159.000 53.000 159.000(ggplot(data.frame(v), aes(v)) +
geom_histogram(binwidth = h_scott, fill = "steelblue", colour = "white") +
labs(title = "Scott") + theme_dspa(10)) |
(ggplot(data.frame(v), aes(v)) +
geom_histogram(binwidth = h_fd, fill = "firebrick", colour = "white") +
labs(title = "Freedman-Diaconis") + theme_dspa(10))
21.2 Problem 2, Does the deflation theorem survive MAR?
§2.7.2 proves the deflation factor under MCAR. Simulate a MAR mechanism, missingness in \(x_2\) depending on observed \(x_1\), and check whether \(\hat\sigma^2_{\text{imp}}/\sigma^2 \approx (1-p)\) still holds.
Solution
set.seed(77)
n <- 5000
x1 <- rnorm(n)
x2 <- 0.8 * x1 + rnorm(n, sd = 0.6) # correlated with x1
# MAR: probability of missing x2 rises with x1
p_miss <- plogis(1.5 * x1 - 0.5)
mis <- runif(n) < p_miss
x2_mar <- x2; x2_mar[mis] <- NA
x2_imp <- x2_mar; x2_imp[is.na(x2_imp)] <- mean(x2_mar, na.rm = TRUE)
# MCAR comparison at the same overall rate
mis_mcar <- sample(n, sum(mis))
x2_mcar <- x2; x2_mcar[mis_mcar] <- NA
x2_mcar_imp <- x2_mcar; x2_mcar_imp[is.na(x2_mcar_imp)] <- mean(x2_mcar, na.rm = TRUE)
p <- mean(mis)
rbind(
MCAR = c(rate = p, predicted = 1 - p,
observed = var(x2_mcar_imp) / var(x2),
mean_bias = mean(x2_mcar_imp) - mean(x2)),
MAR = c(rate = p, predicted = 1 - p,
observed = var(x2_imp) / var(x2),
mean_bias = mean(x2_imp) - mean(x2))
) |> round(4)#> rate predicted observed mean_bias
#> MCAR 0.4124 0.5876 0.5849 -0.0108
#> MAR 0.4124 0.5876 0.4932 -0.326821.3 Problem 3, Contingency-table quantities by hand
For the residence-area × Africa table, compute expected counts, cell
contributions, and standardized Pearson residuals from first principles,
and verify against chisq.test().
Solution
O <- as.matrix(tab)
n <- sum(O)
ri <- rowSums(O); cj <- colSums(O)
E <- outer(ri, cj) / n # expected counts
Cc <- (O - E)^2 / E # contributions
X2 <- sum(Cc)
R <- (O - E) / sqrt(E * outer(1 - ri/n, 1 - cj/n)) # standardized residuals
cat("X2 =", round(X2, 4), " df =", (nrow(O)-1)*(ncol(O)-1),
" p =", signif(pchisq(X2, (nrow(O)-1)*(ncol(O)-1), lower.tail = FALSE), 4), "\n\n")#> X2 = 0.1996 df = 2 p = 0.905#> group
#> residence_area Rest of world Africa
#> Rural -0.432 0.432
#> Total 0.116 -0.116
#> Urban 0.314 -0.314#> [1] TRUE#> [1] TRUECc is non-negative and
unbounded, R is signed and roughly \(N(0,1)\), and the p-value is a tail
probability of \(\chi^2_{(r-1)(c-1)}\).
Only the last is a probability.
21.4 Problem 4, Verify EM’s monotone ascent
Extend em_impute() to record the observed-data
log-likelihood at each iteration and confirm it never decreases.
Solution
obs_loglik <- function(X, mu, Sigma) {
# Sum over rows of the observed-margin multivariate normal log-density
M <- is.na(X); ll <- 0
key <- apply(M, 1L, \(r) paste0(as.integer(r), collapse = ""))
for (rows in split(seq_len(nrow(X)), key)) {
obs <- !M[rows[1L], ]
if (!any(obs)) next
S <- Sigma[obs, obs, drop = FALSE]
Xo <- X[rows, obs, drop = FALSE]
ctr <- sweep(Xo, 2L, mu[obs], "-")
Si <- solve(S)
quad <- rowSums((ctr %*% Si) * ctr)
ll <- ll + sum(-0.5 * (sum(obs) * log(2 * pi) +
determinant(S, logarithm = TRUE)$modulus + quad))
}
as.numeric(ll)
}
# Re-run EM, capturing the likelihood every iteration
Xm <- as.matrix(sim_df)
ll_trace <- numeric(0)
for (k in 1:25) {
f <- em_impute(Xm, tol = 0, max_iter = k, correct_cov = TRUE)
ll_trace[k] <- obs_loglik(Xm, f$mu, f$Sigma)
}
c(monotone = all(diff(ll_trace) >= -1e-6),
first = round(ll_trace[1], 1), last = round(ll_trace[25], 1))#> monotone first last
#> 1.0 -8880.6 -8879.7ggplot(data.frame(it = seq_along(ll_trace), ll = ll_trace), aes(it, ll)) +
geom_line(linewidth = 0.9, colour = "steelblue") + geom_point(size = 1.3) +
labs(title = "Observed-data log-likelihood is non-decreasing",
x = "Iteration", y = "log L(theta | X)") + theme_dspa(10)
21.5 Problem 5, Quantify leakage optimism
Using synthetic data where you control the imbalance ratio, measure how the gap between leaky and correct cross-validated accuracy varies with imbalance.
Solution
library(rsample); library(recipes); library(themis)
leak_gap <- function(n_major, n_minor, d = 8, seed = 1) {
set.seed(seed)
X <- matrix(rnorm((n_major + n_minor) * d), ncol = d)
# Weak true signal, so leakage has room to show
X[seq_len(n_minor), ] <- X[seq_len(n_minor), ] + 0.35
df <- data.frame(X, y = factor(rep(c("min", "maj"), c(n_minor, n_major))))
smote_it2 <- function(z) recipe(y ~ ., data = z) |>
step_smote(y, over_ratio = 1) |> prep() |> juice() |> as.data.frame()
ev <- function(tr, te) {
m <- suppressWarnings(glm(y ~ ., data = tr, family = binomial()))
p <- predict(m, te, type = "response")
mean(factor(ifelse(p > 0.5, levels(tr$y)[2], levels(tr$y)[1]),
levels = levels(te$y)) == te$y)
}
leaked <- vfold_cv(smote_it2(df), v = 5)
a_leak <- mean(vapply(leaked$splits, \(s) ev(analysis(s), assessment(s)), 0))
clean <- vfold_cv(df, v = 5, strata = y)
a_ok <- mean(vapply(clean$splits, \(s) ev(smote_it2(analysis(s)), assessment(s)), 0))
c(ratio = n_major / n_minor, leaky = a_leak, correct = a_ok, gap = a_leak - a_ok)
}
res <- do.call(rbind, lapply(c(2, 5, 10, 20), \(r) leak_gap(600, round(600 / r))))
round(res, 4)#> ratio leaky correct gap
#> [1,] 2 0.7000 0.6878 0.0122
#> [2,] 5 0.7158 0.6972 0.0186
#> [3,] 10 0.7475 0.7045 0.0430
#> [4,] 20 0.7050 0.6905 0.0145ggplot(as.data.frame(res), aes(ratio, gap)) +
geom_line(linewidth = 1, colour = "firebrick") + geom_point(size = 2.4) +
labs(title = "Optimism from resampling before the split",
subtitle = "Gap grows with imbalance: more synthetic points, more memorization",
x = "Majority : minority ratio", y = "Leaky accuracy - correct accuracy") +
theme_dspa()
21.6 Problem 6, Cross-validated bandwidth
Implement least-squares cross-validation for KDE bandwidth and compare against Silverman’s rule on a bimodal sample.
Solution
set.seed(9)
bim <- c(rnorm(400, -2, 0.6), rnorm(400, 2, 0.9))
# ucv() minimizes the unbiased CV criterion; bcv() the biased CV criterion
h_sil <- bw.nrd0(bim); h_ucv <- bw.ucv(bim); h_sj <- bw.SJ(bim)
c(silverman = round(h_sil, 3), ucv = round(h_ucv, 3), sheather_jones = round(h_sj, 3))#> silverman ucv sheather_jones
#> 0.509 0.224 0.224ggplot(data.frame(x = bim), aes(x)) +
geom_histogram(aes(y = after_stat(density)), bins = 50,
fill = "grey88", colour = "white") +
geom_density(aes(colour = "Silverman"), bw = h_sil, linewidth = 1) +
geom_density(aes(colour = "Cross-validated"), bw = h_ucv, linewidth = 1) +
geom_density(aes(colour = "Sheather-Jones"), bw = h_sj, linewidth = 1) +
stat_function(fun = \(z) 0.5 * dnorm(z, -2, 0.6) + 0.5 * dnorm(z, 2, 0.9),
colour = "black", linetype = "dashed", linewidth = 0.8) +
labs(title = "Bandwidth selection on a bimodal sample",
subtitle = "Dashed black: the true mixture density",
x = NULL, y = "Density", colour = NULL) +
theme_dspa()
21.7 Problem 7, Choosing \(K\) by BIC
Fit Gaussian mixtures with \(K=1,\dots,6\) to a crystallography sample and select \(K\) by BIC. Explain why the log-likelihood alone cannot be used.
Solution
v <- crystal[[1]]; v <- v[!is.na(v) & v > 0]
pick_K <- function(v, Kmax = 6) {
do.call(rbind, lapply(1:Kmax, \(k) {
if (k == 1) {
ll <- sum(dnorm(v, mean(v), sd(v), log = TRUE)); npar <- 2
} else {
best <- NULL
for (s in 1:5) {
set.seed(300 + s)
cand <- try(mixtools::normalmixEM(v, k = k, maxit = 2000,
epsilon = 1e-6), silent = TRUE)
if (!inherits(cand, "try-error") &&
(is.null(best) || cand$loglik > best$loglik)) best <- cand
}
ll <- best$loglik; npar <- 3 * k - 1 # k means, k sds, k-1 free weights
}
data.frame(K = k, loglik = ll, npar = npar,
AIC = -2 * ll + 2 * npar,
BIC = -2 * ll + npar * log(length(v)))
}))
}
sel <- pick_K(v); sel |> mutate(across(where(is.numeric), \(z) round(z, 1)))#> number of iterations= 159
#> number of iterations= 186
#> number of iterations= 127
#> number of iterations= 185
#> number of iterations= 158
#> number of iterations= 139
#> number of iterations= 493
#> number of iterations= 723
#> number of iterations= 198
#> number of iterations= 213
#> number of iterations= 247
#> number of iterations= 190
#> number of iterations= 689
#> number of iterations= 1266
#> number of iterations= 614
#> number of iterations= 429
#> WARNING! NOT CONVERGENT!
#> number of iterations= 2000
#> number of iterations= 606
#> WARNING! NOT CONVERGENT!
#> number of iterations= 2000
#> WARNING! NOT CONVERGENT!
#> number of iterations= 2000
#> number of iterations= 139
#> number of iterations= 1799
#> number of iterations= 424
#> WARNING! NOT CONVERGENT!
#> number of iterations= 2000
#> number of iterations= 1007#> BIC-optimal K = 3sel |> pivot_longer(c(AIC, BIC), names_to = "criterion", values_to = "value") |>
ggplot(aes(K, value, colour = criterion)) +
geom_line(linewidth = 1) + geom_point(size = 2.4) +
scale_x_continuous(breaks = 1:6) +
labs(title = "Selecting the number of mixture components",
x = "K", y = "Criterion (lower is better)", colour = NULL) +
theme_dspa()
22 Checkpoint
- A colleague fills 25% missing values with the column mean and reports a 95% CI. In which direction is the interval wrong, and by roughly what factor?
- Why does the optimal histogram bin width shrink like \(n^{-1/3}\) while the optimal kernel bandwidth shrinks like \(n^{-1/5}\)?
- A cell in a contingency table has contribution \(8.4\). Is that a p-value, a residual, or neither, and what would you report instead?
- You are told SMOTE improved cross-validated accuracy from 0.71 to 0.94. What is the first thing to check?
- A KS test gives \(p=0.41\) for a Weibull fit whose parameters were estimated from the same data. Why is that number not usable, and what replaces it?
- Name two things a 26-category pie chart makes harder than a bar chart, using the Cleveland–McGill ordering.
Answers
- Too narrow. Variance is deflated by \(\approx(1-p)=0.75\), so the SE is understated by \(\sqrt{0.75}\approx 0.87\), roughly a 13% narrowing, and type-I error is inflated. Under MAR the point estimate is biased as well.
- The histogram is piecewise constant, so its bias is first order in \(h\); balancing \(h^2\) against \(1/(nh)\) gives \(n^{-1/3}\). A smooth symmetric kernel cancels the first-order bias term, leaving \(O(h^2)\) bias and \(O(h^4)\) squared bias; balancing \(h^4\) against \(1/(nh)\) gives \(n^{-1/5}\). The smoother estimator tolerates a wider window and achieves a better rate, \(O(n^{-4/5})\) vs \(O(n^{-2/3})\).
- Neither. It is \((O-E)^2/E\), one term of the \(\chi^2\) statistic — non-negative and unbounded. Report the standardized Pearson residual (signed, approximately \(N(0,1)\)) to identify which cells drive the association, the overall statistic with its df and p-value, and an effect size such as Cramér’s \(V\).
- Whether SMOTE ran before the train/test split. Synthetic minority points are convex combinations of real minority points; if the parents straddle the split, near-duplicates of test points appear in training. Re-run with resampling inside the training fold only.
- The KS null distribution is distribution-free only for a fully specified reference. With estimated parameters the fitted CDF tracks the empirical CDF too closely, \(D_n\) is too small, and \(p\) is anti-conservative. Replace it with a parametric bootstrap (Lilliefors), and compare candidate models by AIC/BIC rather than by p-value.
- A pie encodes value as angle and area (ranks 3–4), so (i) ranking similar-sized slices is unreliable, and (ii) comparing across pies means comparing angles at different orientations with no shared baseline. A bar chart uses position on a common scale (rank 1) for both tasks.
23 Summary
Part I, data handling and quality
- Ingestion loses information; checksums,
str(), and encoding declarations are how you find out what was lost. - Choose a centrality measure from the scale of the quantity (AM/GM/HM) and the contamination you expect (breakdown point).
- Histogram bin width and kernel bandwidth are estimation problems with derivable optima, \(n^{-1/3}\) and \(n^{-1/5}\). Bandwidth matters; kernel shape barely does.
- A contingency table produces four distinct quantities. Only one of them is a probability.
- Missingness is MCAR, MAR, or MNAR, and the distinction is not testable from the observed data. Mean imputation deflates variance by exactly \((1-p)\).
- EM’s likelihood never decreases, converges linearly, and costs \(O(T(Pd^3+nd^2))\) when you group rows by missingness pattern. Dropping the conditional-covariance term makes it not-EM and biases \(\Sigma\) downward.
- Multiple imputation propagates uncertainty through Rubin’s rules; \(T=\bar U+(1+1/m)B\) is the entire idea.
- Resample inside the training fold. Always.
Part II, exploratory visual analytics
- The grammar of graphics makes plots compositional; Cleveland–McGill makes encoding choice empirical rather than aesthetic.
- Position beats angle beats area beats shading. Choose accordingly.
- p-value displays need \(-\log_{10}\) scaling, multiplicity control, and a visible colour bar.
- LOESS’s
spanis a bias–variance dial with a measurable cost in effective degrees of freedom. - Goodness-of-fit with estimated parameters requires a bootstrap; model choice requires an information criterion.
- Mixture likelihoods are multimodal and unbounded. Restart, order components, and select \(K\) by BIC.
- Reach for 3-D only when rotation carries information a projection would destroy.
Where these threads continue
| Thread | Continues in |
|---|---|
| Missing data in modeling pipelines | Feature selection |
| What accuracy, kappa, ROC actually mean | Model assessment |
| Matrix computing behind covariance and PCA | Linear algebra & regression |
| Mixtures as clustering | Unsupervised clustering |
| Kernel methods and the bandwidth analogy | Black-box methods |
| Time series and 4-D imaging | Longitudinal analysis |
| EM as constrained optimization | Function optimization |
| Streaming, out-of-core, deployment | Specialized topics |
24 Chapter roadmap
Forward references in this chapter resolve here.
- Chapter 1,
Foundations. R toolchain, reproducibility conventions,
dspa_read(), simulation. - Linear algebra, matrix computing, regression. Covariance structure, least squares, regression and model trees.
- Model assessment, validation, improvement. Accuracy, sensitivity/specificity, \(\kappa\), ROC/AUC, cross-validation, tuning.
- Variable importance and feature selection. Filters, wrappers, embedded methods, regularization, FDR control.
- Unsupervised clustering. k-means, hierarchical, spectral, and Gaussian mixture models.
- Black-box methods. Neural networks, SVM (kernels and bandwidth), ensembles.
- Longitudinal and time-series analysis. ARIMA/ARIMAX, mixed models, GEE, recurrent architectures.
- Function optimization. Gradient descent, Newton methods, EM as an optimization scheme, Bayesian optimization.
- Specialized topics. Databases, streaming, parallel and out-of-core computing, deployment.
25 Appendix A: databases
See DSPA
Appendix 5 for detail. The modern interface is DBI plus
a backend driver; RODBC is legacy.
library(DBI)
con <- dbConnect(RSQLite::SQLite(), file.path(tempdir(), "demo.sqlite"))
dbWriteTable(con, "iris", iris)
dbListTables(con)
dbGetQuery(con, "SELECT Species, AVG(`Sepal.Length`) AS mean_sl
FROM iris GROUP BY Species")
dbDisconnect(con)
# Other backends: RPostgres::Postgres(), RMariaDB::MariaDB(), odbc::odbc(),
# duckdb::duckdb() for larger-than-memory analytics.
# dbplyr lets you write dplyr verbs that are translated to SQL and executed
# in the database, so the data never enters R's memory.26 Appendix B: case study, traumatic brain injury
tmp <- tempfile(fileext = ".xlsx")
download.file("https://umich.instructure.com/files/416270/download?download_frd=1",
tmp, mode = "wb", quiet = TRUE)
df_TBI <- openxlsx::read.xlsx(tmp, sheet = "Sheet1", skipEmptyRows = TRUE)
dim(df_TBI)#> [1] 46 19df_clean <- df_TBI |>
naniar::replace_with_na_all(condition = ~ .x %in% c("NA", ".")) |>
as.data.frame()
df_clean <- df_clean[, -c(3:4)] |> tidyr::drop_na()
rownames(df_clean) <- as.character(df_clean[[1]])
ids <- rownames(df_clean)
df_clean <- df_clean[, -1]
df_clean <- as.data.frame(lapply(df_clean, as.numeric))
# When as.data.frame() converts the data, R automatically sanitized any column names that started with digits (like 6m.gose and 2013.gose) by prepending an "X" to make them valid variable names. (Notice how 6m.gose also became X6m.gose).
colnames(df_clean) <- sub("^X", "", colnames(df_clean))
rownames(df_clean) <- ids
df_clean <- df_clean[, c("age", "2013.gose", "skull.fx", "temp.injury",
"surgery", "acute.sz")]
df_scaled <- as.data.frame(scale(df_clean))
dim(df_scaled)#> [1] 23 6hc_tbi <- hclust(dist(df_scaled), method = "average")
dd_tbi <- ggdendro::dendro_data(as.dendrogram(hc_tbi), type = "rectangle")
ggplot() +
geom_segment(data = ggdendro::segment(dd_tbi),
aes(x, y, xend = xend, yend = yend), linewidth = 0.5) +
labs(title = "TBI subjects: hierarchical clustering (average linkage, scaled)",
x = NULL, y = "Merge height") +
theme_dspa(10) + theme(axis.text.x = element_blank())
# Interactive: brushing and persistent highlighting of subtrees
plotly::plot_dendro(as.dendrogram(hc_tbi), height = 550) |>
layout(title = "TBI dendrogram (click a node to highlight its subtree)") |>
hide_legend() |>
highlight(persistent = TRUE, dynamic = TRUE)#> $`k=2`
#>
#> 1 2
#> 19 4
#>
#> $`k=3`
#>
#> 1 2 3
#> 6 13 4
#>
#> $`k=4`
#>
#> 1 2 3 4
#> 6 11 4 2
#>
#> $`k=5`
#>
#> 1 2 3 4 5
#> 6 11 3 1 2groups3 <- cutree(hc_tbi, k = 3)
aggregate(df_scaled, list(cluster = groups3), median) |>
mutate(across(where(is.numeric), \(z) round(z, 2)))fit_tbi <- lm(`2013.gose` ~ age, data = df_scaled)
ggplot(df_scaled, aes(age, `2013.gose`)) +
geom_point(size = 2, alpha = 0.7, colour = "steelblue") +
geom_smooth(method = "lm", formula = y ~ x, colour = "firebrick", se = TRUE) +
labs(title = sprintf("Standardized GOSE vs. age (r = %.3f)",
cor(df_scaled$`2013.gose`, df_scaled$age)),
subtitle = "Variables are z-scored, so units are standard deviations",
x = "Age (z)", y = "2013 GOSE (z)") +
theme_dspa()
# --- Interactive equivalent -----------------------------------------------
plot_ly(df_scaled, x = ~age, y = ~`2013.gose`, type = "scatter",
mode = "markers", name = "Data") |>
add_lines(x = ~age, y = fitted(fit_tbi), name = "Linear model") |>
layout(title = paste0("Correlation(2013.gose, age) = ",
round(cor(df_scaled$`2013.gose`, df_scaled$age), 3)))27 Appendix C: additional
ggplot2 examples
27.1 Housing price index
hpi <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_Dinov_091609_SnP_HomePriceIndex") |>
html_nodes("table") |> _[[1]] |> html_table()
hpi <- hpi |>
mutate(Date = lubridate::parse_date_time(paste(Year, Month), "ym")) |>
dplyr::select(-Year, -Month, -1)
hpi_long <- hpi |>
pivot_longer(-Date, names_to = "city", values_to = "index") |>
mutate(index = as.numeric(index)) |>
filter(!is.na(index))
ggplot(hpi_long, aes(Date, index, colour = city)) +
geom_line(linewidth = 0.7) +
labs(title = "US housing price index, 1991-2009",
x = NULL, y = "HPI", colour = NULL) +
theme_dspa(10) + theme(legend.position = "right",
legend.text = element_text(size = 7))
# --- Interactive equivalent -----------------------------------------------
plot_ly(hpi_long, x = ~Date, y = ~index, color = ~city,
type = "scatter", mode = "lines+markers") |>
layout(title = "US Housing Price Index (1991-2009)",
yaxis = list(title = "HPI"), legend = list(orientation = "h"))hpi_sf <- hpi |> dplyr::select(la = `CA-LosAngeles`, sf = `CA-SanFrancisco`) |>
mutate(across(everything(), as.numeric)) |> filter(!is.na(la), !is.na(sf))
ggplot(hpi_sf, aes(la, sf)) +
geom_point(alpha = 0.6, colour = "steelblue") +
geom_smooth(method = "lm", formula = y ~ x, colour = "magenta", linewidth = 1.2) +
labs(title = "San Francisco vs. Los Angeles home price index",
subtitle = sprintf("OLS fit; R-squared = %.3f",
summary(lm(sf ~ la, hpi_sf))$r.squared),
x = "Los Angeles HPI", y = "San Francisco HPI") +
theme_dspa()
pairs_df <- hpi[, 10:15] |> mutate(across(everything(), as.numeric))
names(pairs_df) <- c("Atlanta", "Chicago", "Boston", "Detroit",
"Minneapolis", "Charlotte")
GGally::ggpairs(pairs_df) + theme_dspa(8)
27.2 Los Angeles neighborhoods
la <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_LA_Neighborhoods_Data") |>
html_nodes("table") |> _[[2]] |> html_table()
la <- la |> mutate(across(c(Longitude, Latitude, Population, Income), as.numeric))
# Note: bare column names inside aes(), never aes(x = df$col), the latter
# breaks faceting and is warned about in current ggplot2.
ggplot(la, aes(Longitude, Latitude)) +
geom_point(aes(size = Population, fill = Income), shape = 21,
stroke = 0.2, alpha = 0.75, colour = "grey20") +
geom_text(aes(label = LA_Nbhd), size = 1.6, hjust = 0.5, vjust = 2,
check_overlap = TRUE) +
scale_size_area(max_size = 9) +
scale_fill_distiller(palette = "RdBu", na.value = "white", name = "Income") +
coord_fixed(ratio = 1) +
labs(title = "LA neighborhoods: location, population, income") +
theme_dspa(10)
27.3 Latin letter frequency
letter_long <- letter |>
pivot_longer(-Letter, names_to = "language", values_to = "freq") |>
mutate(freq = as.numeric(freq)) |> filter(!is.na(freq))
ggplot(letter_long, aes(Letter, freq, fill = language)) +
geom_col() +
labs(title = "Latin letter frequency across languages",
x = NULL, y = "Count", fill = NULL) +
theme_dspa(10) + theme(legend.text = element_text(size = 7))
# --- Interactive equivalent -----------------------------------------------
plot_ly(letter_long, x = ~Letter, y = ~freq, type = "bar",
name = ~language, color = ~language) |>
layout(yaxis = list(title = "Count"), barmode = "stack")More at the SOCR letter-frequency page.
28 Session information
#> R version 4.3.3 (2024-02-29 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 26200)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=English_United States.utf8
#> [2] LC_CTYPE=English_United States.utf8
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.utf8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] oro.nifti_0.11.4 rvest_1.0.4 mice_3.17.0 naniar_1.1.0
#> [5] DT_0.33 plotly_4.12.0 patchwork_1.3.0 tidyr_1.3.1
#> [9] dplyr_1.1.4 ggplot2_4.0.1
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 ggdendro_0.2.0 rstudioapi_0.18.0
#> [4] jsonlite_1.8.9 shape_1.4.6.1 magrittr_2.0.3
#> [7] spatstat.utils_3.0-5 jomo_2.7-6 farver_2.1.2
#> [10] corrplot_0.92 nloptr_2.1.1 rmarkdown_2.31
#> [13] vctrs_0.6.5 minqa_1.2.7 spatstat.explore_3.2-7
#> [16] base64enc_0.1-3 htmltools_0.5.8.1 forcats_1.0.0
#> [19] curl_6.2.0 broom_1.0.6 Formula_1.2-5
#> [22] mitml_0.4-5 sass_0.4.9 bslib_0.9.0
#> [25] htmlwidgets_1.6.4 plyr_1.8.9 lubridate_1.9.3
#> [28] cachem_1.1.0 lifecycle_1.0.5 iterators_1.0.14
#> [31] pkgconfig_2.0.3 Matrix_1.6-5 R6_2.6.1
#> [34] fastmap_1.2.0 digest_0.6.37 selectr_0.4-2
#> [37] colorspace_2.1-1 GGally_2.2.1 ps_1.9.0
#> [40] tensor_1.5 crosstalk_1.2.1 Hmisc_5.1-3
#> [43] labeling_0.4.3 spatstat.sparse_3.0-3 timechange_0.3.0
#> [46] httr_1.4.7 polyclip_1.10-6 abind_1.4-5
#> [49] mgcv_1.9-1 compiler_4.3.3 withr_3.0.2
#> [52] htmlTable_2.4.2 S7_0.2.1 backports_1.5.0
#> [55] viridis_0.6.5 UpSetR_1.4.0 ggstats_0.6.0
#> [58] dendextend_1.17.1 pan_1.9 MASS_7.3-60.0.1
#> [61] tools_4.3.3 chromote_0.4.0 foreign_0.8-87
#> [64] otel_0.2.0 zip_2.3.1 visdat_0.6.0
#> [67] nnet_7.3-19 goftest_1.2-3 glue_1.8.0
#> [70] nlme_3.1-165 promises_1.3.2 grid_4.3.3
#> [73] checkmate_2.3.1 cluster_2.1.6 generics_0.1.3
#> [76] gtable_0.3.6 spatstat.data_3.1-2 websocket_1.4.1
#> [79] data.table_1.16.4 xml2_1.3.6 spatstat.geom_3.2-9
#> [82] foreach_1.5.2 pillar_1.10.1 stringr_1.5.1
#> [85] later_1.4.1 splines_4.3.3 lattice_0.22-6
#> [88] survival_3.7-0 deldir_2.0-4 tidyselect_1.2.1
#> [91] knitr_1.51 gridExtra_2.3 xfun_0.52
#> [94] stringi_1.8.4 fftwtools_0.9-11 lazyeval_0.2.2
#> [97] yaml_2.3.10 boot_1.3-30 evaluate_1.0.3
#> [100] codetools_0.2-20 RNifti_1.6.1 tibble_3.2.1
#> [103] cli_3.6.3 rpart_4.1.23 processx_3.8.6
#> [106] jquerylib_0.1.4 Rcpp_1.0.14 spatstat.random_3.2-3
#> [109] bitops_1.0-7 lme4_1.1-35.5 glmnet_4.1-8
#> [112] viridisLite_0.4.2 scales_1.4.0 openxlsx_4.2.5.2
#> [115] purrr_1.0.2 rlang_1.1.5