The greatest value of a picture is when it forces us to notice what we never expected to see.
John W. Tukey, Exploratory Data Analysis (1977)
Stat 545 Part VII; blog posts 16-plotsfrompurrr, 38-tableplacementrmarkdown; zzlongplot for longitudinal examples (Chapter 31).
Answer the following questions to see if you can bypass this chapter. You can find the answers at the end of the chapter in Section 22.17.
ggplot objects?patchwork::wrap_plots() provide that facet_wrap() does not?By the end of this chapter you should be able to:
ggplot2 figures for continuous, categorical, and time-varying outcomes.purrr::map() over a list of groups.patchwork and arrange shared legends cleanly.theme_*() function and reuse it across all figures in a project.zzlongplot and zztable1 rather than rolling your own.A figure is the only part of a paper that most readers will examine closely, and it is the part over which the analyst has the most control and exercises the least deliberation. The default output of a plotting call is a decision, made by somebody who had never seen your data, about what the reader’s eye should be drawn to. Accepting it is also a decision.
A figure that survives peer review is one whose data, encoding, and aesthetic choices all serve the single claim the figure exists to make. This chapter is mostly about the mechanics of getting there, since the mechanics are what can be taught in a chapter; the fuller principles of effective visualization are taken up in the companion textbook (chapters 15–16 of Statistical Computing in the Age of AI), and the ways a figure can mislead without printing a false number are taken up later, in Chapter 34.
The combination ggplot2 + purrr + patchwork is the modern stack for going from a single exploratory plot to a polished multi-panel figure: ggplot2 for the plot, purrr for iteration over groups, patchwork for composition.
Software handles the mechanics; ggplot2 will draw whatever you specify. Whether what you specify communicates or misleads is not something it can check:
One plot, one message. A figure trying to convey three relationships succeeds at none. Three figures each conveying one are usually clearer. Resist the ‘pack everything into one figure to save space’ impulse.
Encode the claim, not the data. A scatter plot with 5,000 points obscures any pattern. The same relationship as a hex-binned plot or a 2D density contour shows the structure without the visual noise. Choose encoding for legibility, not faithfulness to the raw data.
Show the distribution, not a summary of it. A bar chart of group means with error bars is the default figure in the clinical literature and one of the least informative. It reports two numbers per group, hides skewness, bimodality, and outliers, and devotes most of its ink to a rectangle whose height is measured from an arbitrary zero. Harrell’s chapter on descriptive statistics and graphics (Harrell, 2025a, ch. 4) makes the case at length: a bar has a high ink-to-information ratio, its width carries no information at all, and graphing the raw data is usually essential rather than optional. With a modest number of observations per group, plot the points: geom_jitter() or ggbeeswarm::geom_quasirandom(), with the summary overlaid rather than substituted for the data. With many observations, geom_boxplot() or geom_violin() still shows a shape where a bar shows none. This is not in tension with the previous judgment: encoding the claim means choosing the display that makes the claim legible, and ‘the groups differ’ is a claim about distributions.
Consistency across the project. All figures in a paper should share fonts, color palette, and theme defaults. A house theme function applied via theme_set() makes consistency the default. Without it, ggplot2 produces nine variants of gray.
Display uncertainty. A point estimate without a CI invites over-trust. A regression line without a confidence band the same. The tools are simple (geom_errorbar, se = TRUE on geom_smooth); the discipline is the analyst’s.
Ask what the reader’s eye is measuring. A bar encodes length, so a truncated baseline makes a small difference fill the panel; a map encodes area, and land does not have a disease; a bubble’s radius is read as area, so a radius proportional to a count exaggerates it. None of these prints a false number, and all of them mislead. Chapter 34 returns to these failure modes and their remedies at length; the habit to carry into every figure is to name the visual channel the reader is decoding before choosing it.
These judgments are what make figures publishable.
The pattern: split data by a grouping variable, apply a plotting function to each group, collect into a named list.
library(tidyverse)
library(palmerpenguins)
penguins_clean <- na.omit(penguins)
# split by species
plots <- penguins_clean |>
split(penguins_clean$species) |>
imap(\(d, name) {
ggplot(d, aes(flipper_length_mm, body_mass_g)) +
geom_point() +
geom_smooth(method = "lm") +
labs(title = name,
x = "Flipper length (mm)", y = "Body mass (g)") +
theme_minimal()
})
# inspect one
plots$Adelie
plots holds all three.
The same iteration writes every panel to disk. This is a side effect, so we reach for walk2 rather than map2:
# save all three panels as publication-quality PDFs
walk2(plots,
paste0("figures/penguin-", names(plots), ".pdf"),
\(p, path) ggsave(path, p, width = 6, height = 4))Three things to internalize:
split() returns a named list. The names propagate through imap and walk2.imap and walk2 use both the value and the name. Useful for titles, file paths, and identifying which plot is which.walk2 for side effects. Returns invisibly; map2 would build a list of ggsave’s NULL returns.For a dplyr-native version:
penguins_clean |>
group_by(species) |>
group_map(\(d, k) {
ggplot(d, aes(flipper_length_mm, body_mass_g)) +
geom_point() +
labs(title = k$species)
})
# returns a list of plots, with the group keys in `k`patchwork composes distinct plots into one figure. Where facet_wrap produces small multiples of the same plot, wrap_plots (or the + and / operators) combines independent plots that may differ in everything.
library(patchwork)
p1 <- ggplot(penguins_clean,
aes(flipper_length_mm, body_mass_g, color = species)) +
geom_point() +
labs(title = "A. Body mass vs flipper length")
p2 <- ggplot(penguins_clean,
aes(species, body_mass_g, fill = species)) +
geom_boxplot() +
labs(title = "B. Body mass by species") +
guides(fill = "none") # hide redundant fill legend
p3 <- ggplot(penguins_clean,
aes(bill_length_mm, bill_depth_mm, color = species)) +
geom_point() +
labs(title = "C. Bill morphology")
# 2x2 grid with shared legend at bottom
(p1 + p2) / (p3 + plot_spacer()) +
plot_layout(guides = "collect") &
theme(legend.position = "bottom")
Operators:
+ juxtaposes side by side./ stacks vertically.| is equivalent to + (horizontal).& applies a theme to every plot in the composition.plot_layout(guides = "collect") collects duplicate legends into one.For more control, plot_layout takes widths, heights, nrow, ncol, byrow, and other arguments.
facet_wrap for small multiples of the same plot; patchwork for distinct plots in one figure.
A function that returns a theme() object, applied once per session:
theme_practicum <- function(base_size = 11) {
theme_minimal(base_size = base_size) +
theme(
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "gray40"),
axis.title = element_text(face = "bold"),
axis.text = element_text(color = "gray20"),
panel.grid.minor = element_blank(),
legend.position = "bottom",
strip.background = element_rect(fill = "gray95", color = NA),
strip.text = element_text(face = "bold")
)
}
# apply globally
theme_set(theme_practicum())
# update geom defaults to match
update_geom_defaults("point", list(size = 1.5, alpha = 0.7))
update_geom_defaults("line", list(linewidth = 0.6))For project consistency, define theme_practicum() in a setup chunk at the top of every analysis file, or in a project-level helper that gets sourced. Every plot then matches without per-plot theme calls.
For palettes, define and reuse:
practicum_palette <- c(
"Adelie" = "#1f4e79",
"Chinstrap" = "#9d2235",
"Gentoo" = "#2e8b57"
)
scale_color_practicum <- function(...)
scale_color_manual(values = practicum_palette, ...)
# in plots
ggplot(...) + ... + scale_color_practicum()For color-blind safety, verify the palette with the Coblis simulator (colorbrewer2.org/learnmore/colorblind-simulator.html) or colorBlindness::cvdPlot().
ggsave() writes to many formats, picked by extension:
# vector for journal submission
ggsave("figure1.pdf", plot = p, width = 6, height = 4,
device = cairo_pdf)
# raster for Word, slides, web
ggsave("figure1.png", plot = p, width = 6, height = 4,
dpi = 300)
# vector for the web (preserves quality on zoom)
ggsave("figure1.svg", plot = p, width = 6, height = 4)DPI considerations:
For LaTeX submission, prefer PDF (vector). For Word, prefer PNG at 300 DPI. For web, SVG (vector).
device = cairo_pdf ensures non-default fonts are embedded so the PDF renders correctly on a system that lacks the font.
For journals, check size requirements (typically column widths of 85–90 mm or 170–180 mm) and produce figures at exactly the target size, not larger.
zzlongplot and zztable1The in-house zzlongplot package provides opinionated longitudinal-study plots: spaghetti plots with confidence ribbons, mean trajectories with group-specific overlays, ICC-aware error bars. The defaults match this practicum’s conventions.
library(zzlongplot)
zzlongplot::lplot(adni_long,
form = ADAS ~ VISIT | DX,
cluster_var = "RID",
xlab = "Visit",
ylab = "ADAS-Cog")When to reach for zzlongplot:
When to write your own ggplot2:
The zztable1 package similarly handles Table 1 generation; covered in the analysis-plan chapter (Chapter 24).
theme_set(theme_practicum())
fit <- lm(body_mass_g ~ flipper_length_mm + species,
data = na.omit(penguins))
diag <- augment(fit)
p_resid <- ggplot(diag, aes(.fitted, .resid)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "loess", se = FALSE,
color = "red", linewidth = 0.6) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "A. Residuals vs. fitted",
x = "Fitted (g)", y = "Residual")
p_qq <- ggplot(diag, aes(sample = .std.resid)) +
geom_qq(alpha = 0.5) +
geom_qq_line(color = "red") +
labs(title = "B. Normal Q-Q",
x = "Theoretical", y = "Standardized residual")
p_cook <- ggplot(diag, aes(seq_len(nrow(diag)), .cooksd)) +
geom_col(width = 0.6) +
geom_hline(yintercept = 4 / nrow(diag),
linetype = "dashed", color = "red") +
labs(title = "C. Cook's distance",
x = "Observation", y = "Cook's distance")
(p_resid + p_qq) / p_cook +
plot_layout(heights = c(1, 0.7))
To write the composed figure to disk for a manuscript, add a ggsave call after building it:
ggsave("figures/diagnostics.pdf", width = 7, height = 6,
device = cairo_pdf)The composition reads naturally: top row two related diagnostics, bottom row a single observation-level diagnostic. House theme applied via theme_set(theme_practicum()). Saved as a publication-quality PDF.
The penguin data are cross-sectional, and cross-sectional scatter plots are the easy case. The figure a biostatistician is asked for most often is harder: the trajectory of a repeated-measures outcome over time, by treatment arm, with the individual patients visible behind the group means. To illustrate, we simulate a small two-arm trial in which a cognitive score (higher is worse, in the manner of ADAS-Cog) declines more slowly on the active arm than on placebo.
set.seed(1)
n_id <- 60
months <- c(0, 6, 12, 18, 24)
sim <- expand_grid(id = seq_len(n_id), month = months) |>
mutate(arm = if_else(id <= n_id / 2, "Placebo", "Active")) |>
group_by(id) |>
mutate(
slope = if_else(first(arm) == "Placebo", 0.45, 0.20) +
rnorm(1, 0, 0.08),
adas = 20 + rnorm(1, 0, 4) + slope * month + rnorm(n(), 0, 1.5)
) |>
ungroup()
arm_means <- sim |>
group_by(arm, month) |>
summarize(mean = mean(adas), se = sd(adas) / sqrt(n()),
.groups = "drop")
pal <- c(Placebo = "#9d2235", Active = "#1f4e79")
ggplot(sim, aes(month, adas, color = arm)) +
geom_line(aes(group = id), alpha = 0.15) +
geom_ribbon(data = arm_means,
aes(month, ymin = mean - 1.96 * se,
ymax = mean + 1.96 * se, fill = arm),
alpha = 0.2, color = NA, inherit.aes = FALSE) +
geom_line(data = arm_means, aes(month, mean, color = arm),
linewidth = 1.1, inherit.aes = FALSE) +
scale_color_manual(values = pal) +
scale_fill_manual(values = pal) +
labs(x = "Month", y = "ADAS-Cog (higher is worse)",
color = "Arm", fill = "Arm",
title = "Simulated cognitive decline by treatment arm")
zzlongplot package produces by default.
The individual trajectories, drawn at low opacity, convey the within-arm heterogeneity that a plot of group means alone would hide. The confidence ribbon on the arm means keeps the reader from over-reading a separation that is within the noise. The real ADNI version of this figure appears in Chapter 31.
Many students arrive with more Python than R, and a good deal of collaborative biostatistics now crosses the two languages. Where a comparison is instructive we show the Python equivalent alongside the R. Quarto executes Python chunks through reticulate, and the r object makes an R data frame available to Python without writing an intermediate file. Here is the body-mass-versus-flipper-length scatter in pandas and matplotlib.
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
df = r.penguins_clean
colors = {"Adelie": "#1f4e79", "Chinstrap": "#9d2235",
"Gentoo": "#2e8b57"}
fig, ax = plt.subplots(figsize=(6, 4))
for species, sub in df.groupby("species", observed=True):
ax.scatter(sub["flipper_length_mm"], sub["body_mass_g"],
label=species, color=colors[species], alpha=0.7, s=18)
ax.set_xlabel("Flipper length (mm)")
ax.set_ylabel("Body mass (g)")
ax.legend(title="Species")
fig.tight_layout()
plt.show()
r object, so no intermediate file is written.
For a reader who wants the grammar of graphics in Python rather than the imperative matplotlib style, the plotnine package reproduces the ggplot2 API almost verbatim (from plotnine import ggplot, aes, geom_point). The choice between them is a matter of taste and of what the rest of the analysis is written in. The judgment the figure requires, one message per plot, uncertainty shown, color chosen for legibility, does not change with the language.
LLMs handle ggplot well; the trap is busy plots that encode too much.
Prompt 1: drafting a plot. Describe the data and the question, ask: ‘write a ggplot that addresses the question. Use a color-blind-safe palette and clear axis labels with units.’
What to watch for. The default LLM plot tends to be busy: too many aesthetics encoded, default ggplot theme. Push for clarity. Multiple iterations of ‘simpler’ tend to improve.
Verification. Render the plot. Ask whether a reader who has never seen the data could state the message in one sentence.
Prompt 2: combining plots. Describe four plots, ask: ‘combine these into a 2x2 grid with shared legend, panel labels A through D.’
What to watch for. patchwork::plot_layout(guides = "collect") for the legend. Panel labels usually go in labs(title = "A. ...") rather than tag_levels. The LLM may use either; both work.
Verification. Render the combined plot. Are legends shared? Are panels labeled in the right order?
Prompt 3: theme function. Ask: ‘write a theme_practicum() function with serif body text, sans-serif axis labels, and a color-blind-safe default palette.’
What to watch for. The output is a starting point. Test with a few plots; iterate.
Verification. Apply to several plots; ensure consistency. Verify palette via Coblis or colorBlindness::cvdPlot().
Three habits carry from the first figure to the last:
theme_set(theme_practicum()) at the top of every analysis script.palmerpenguins data, build a three-panel patchwork: (a) scatter of body mass vs flipper length colored by species; (b) residuals from a linear fit of (a); (c) QQ plot of residuals. Share the species legend across all three.plot_per_site(data, site_col, outcome) that returns a named list of ggplot objects (one per site level) and a helper that saves each to figures/<site>.pdf.theme_practicum() and apply it to three plots from any prior exercise. Verify the plots render consistently in HTML, PDF, and Word output.clauswilke.com/dataviz, effective- visualization principles.patchwork, cowplot, and gganimate package vignettes.ggplot2 graphics chapter written for biomedical reporting, including interactive figures in Quarto output.split(data, data$group) |> map(\(d) ggplot(d, aes(x, y)) + geom_point()) to produce a named list of plots, one per group. Alternatively, in a tidyverse pipeline: data |> group_by(group) |> group_map(\(d, k) ggplot(d, aes(x, y)) + geom_point()) (note: group_map returns a list, with names taken from the group keys if you pass .keep = TRUE and post-process).facet_wrap() splits a single plot into panels by a factor. Every panel shares the same aes mapping, same geom, same scales. patchwork::wrap_plots() composes distinct plots that can differ in every respect (different geoms, different data, different scales). Use facet_wrap() for small-multiples of the same plot; wrap_plots() for multi-panel figures showing different views.