Data Visualization

1. Why Visualize Before You Model?

Anscombe’s Quartet (1973) is the classic demonstration that four datasets with identical summary statistics can look completely different. The same principle applies constantly in environmental data: means and standard deviations are not enough.

[]

# Anscombe's Quartet is built into base R
data(anscombe)

# Reshape to long format for ggplot
ansc_long <- data.frame(
  set = rep(1:4, each = 11),
  x   = c(anscombe$x1, anscombe$x2, anscombe$x3, anscombe$x4),
  y   = c(anscombe$y1, anscombe$y2, anscombe$y3, anscombe$y4)
)

# All four datasets have (nearly) identical summary statistics
ansc_long |>
  group_by(set) |>
  summarise(
    mean_x = round(mean(x), 2),
    mean_y = round(mean(y), 2),
    sd_x   = round(sd(x),   2),
    sd_y   = round(sd(y),   2),
    r      = round(cor(x, y), 2)
  ) |>
  knitr::kable(
    col.names = c("Dataset", "Mean X", "Mean Y", "SD X", "SD Y", "Correlation r"),
    caption   = "Anscombe's Quartet: four structurally different datasets with identical summary statistics."
  )

Anscombe’s Quartet: four structurally different datasets with identical summary statistics.
Dataset	Mean X	Mean Y	SD X	SD Y	Correlation r
1	9	7.5	3.32	2.03	0.82
2	9	7.5	3.32	2.03	0.82
3	9	7.5	3.32	2.03	0.82
4	9	7.5	3.32	2.03	0.82

Anscombe’s Quartet. Identical means, SDs, and correlations — but completely different structures. Dataset 1 alone suits linear regression; the others reveal a curve, an outlier, and a leverage point. Without plotting, you would never know.

📌 Visualize → Summarize → Model

This is the correct order of operations. Plotting first lets you: detect outliers and data entry errors before they corrupt an analysis; check whether a linear model is even appropriate; identify subgroups or bimodal distributions that a single mean obscures; and communicate results to an audience that cannot read a regression table.

2. The Grammar of Graphics

ggplot2 is built on a systematic grammar — a set of independent components that can be composed to describe any statistical graphic. Understanding the grammar means you are never guessing syntax; you are assembling known parts.

The seven layers

Every ggplot2 plot is built from up to seven layers stacked on top of each other:

Layer	Purpose	Example
Data	The data frame being plotted	`data = river_df`
Aesthetics `aes()`	Map variables to visual properties: x, y, colour, fill, size, shape	`aes(x = temp, y = DO, colour = site)`
Geometry `geom_*()`	The geometric object: point, line, bar, boxplot, ribbon…	`geom_point()`, `geom_line()`, `geom_boxplot()`
Statistics `stat_*()`	Statistical transformation before drawing: bin, smooth, summarise	`stat_smooth()`, `stat_bin()`
Scales `scale_*()`	Control axis limits, colour palettes, axis labels	`scale_colour_manual()`, `scale_x_log10()`
Facets `facet_*()`	Split the plot into small multiples by a grouping variable	`facet_wrap(~ year)`, `facet_grid(site ~ season)`
Theme `theme()`	All non-data ink: fonts, grid lines, background, margins	`theme_minimal()`, `theme(legend.position = 'bottom')`

Minimal working example

set.seed(1)
river <- data.frame(
  temp_C = runif(50, 12, 24),
  DO_mgl = 14 - 0.35 * runif(50, 12, 24) + rnorm(50, 0, 0.5),
  site   = sample(c("Upstream", "Downstream"), 50, replace = TRUE)
)

# Layer 1: data.  Layer 2: aesthetics.  Layer 3: geometry.
ggplot(data = river,                        # Layer 1 — data
       aes(x = temp_C, y = DO_mgl,          # Layer 2 — aesthetics
           colour = site)) +
  geom_point(size = 2.5, alpha = 0.8)       # Layer 3 — geometry

A minimal ggplot2 scatterplot using only three required layers: data, aesthetics, and geometry. Every additional layer is optional refinement.

Now add refinement layers one at a time:

ggplot(river, aes(x = temp_C, y = DO_mgl, colour = site)) +
  geom_point(size = 2.2, alpha = 0.75) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 0.8) +   # stat layer
  scale_colour_manual(                                          # scale layer
    values = c("Upstream" = "#4a7c5e", "Downstream" = "#c8832a"),
    name   = "Site"
  ) +
  labs(                                                         # labels
    x     = "Water temperature (°C)",
    y     = "Dissolved oxygen (mg/L)",
    title = "DO vs. Temperature — Two River Sites"
  ) +
  theme_minimal(base_size = 12) +                              # theme layer
  theme(
    plot.title  = element_text(family = "serif", face = "bold",
                               colour = "#1e3a2f"),
    legend.position = "top"
  )

The same plot with scales, a smoother, and a theme added. Each layer is independent — remove any one without breaking the others.

⚙️ The ggplot2 template — memorise this
ggplot(data, aes(…)) + geom_() + scale_() + labs() + theme_()

Only ggplot() and one geom_() are required. Everything else is optional. Add layers with +, never |>.

3. One Continuous Variable

When you have a single continuous measurement, your first goal is to understand its distribution — shape, centre, spread, and any unusual features.

Histogram

The histogram bins the data and plots counts (or density) per bin. The choice of bin width is consequential: too wide hides structure, too narrow produces noise.

set.seed(5)
# Daily maximum ozone (ppb) at a monitoring station — one year
ozone <- data.frame(
  ppb  = c(rnorm(200, mean = 38, sd = 9),    # typical summer days
           rnorm(165, mean = 22, sd = 6)),    # winter / clean days
  season = c(rep("Summer", 200), rep("Winter", 165))
)

p1 <- ggplot(ozone, aes(ppb)) +
  geom_histogram(bins = 10, fill = "#4a7c5e", colour = "white", alpha = 0.9) +
  labs(title = "bins = 10 (too coarse)", x = "Ozone (ppb)", y = "Count") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(family = "serif", colour = "#1e3a2f"))

p2 <- ggplot(ozone, aes(ppb)) +
  geom_histogram(binwidth = 2, fill = "#4a7c5e", colour = "white", alpha = 0.9) +
  labs(title = "binwidth = 2 (reveals bimodality)", x = "Ozone (ppb)", y = "Count") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(family = "serif", colour = "#1e3a2f"))

gridExtra::grid.arrange(p1, p2, ncol = 2)

Distribution of daily maximum ozone. Using bins = 30 (the default) vs. a carefully chosen binwidth — the narrower bins reveal the bimodal structure that the coarser histogram obscures.

Density plot

A kernel density estimate (KDE) is a smoothed histogram — useful when you want to compare distributions without the bin-width problem. Overlay multiple groups using transparency (alpha).

ggplot(ozone, aes(ppb, fill = season, colour = season)) +
  geom_density(alpha = 0.35, linewidth = 0.9) +
  scale_fill_manual(values  = c("Summer" = "#c8832a", "Winter" = "#4a7c5e")) +
  scale_colour_manual(values = c("Summer" = "#c8832a", "Winter" = "#4a7c5e")) +
  labs(x = "Ozone (ppb)", y = "Density",
       title  = "Ozone Distribution by Season",
       fill   = NULL, colour = NULL) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title      = element_text(family = "serif", face = "bold", colour = "#1e3a2f"),
    legend.position = "top"
  )

Kernel density estimates of ozone by season. Overlapping semi-transparent fills allow direct comparison without stacking. The bimodal overall distribution is clearly driven by the two seasons.

💡 Histogram vs. density plot

Use a histogram when showing actual counts matters — for instance, to communicate sample size visually alongside distribution shape. Use a density plot when overlaying two or more groups for comparison; raw count histograms of different-sized groups are visually misleading. To get the best of both: geom_histogram(aes(y = after_stat(density))) puts histograms on the density scale so an overlaid geom_density() aligns correctly.

4. One Categorical Variable

When you need to show counts or proportions across categories, the bar chart is the natural choice — provided you follow a few conventions that prevent distortion.

set.seed(8)
# Land-cover classification from remote sensing (n = 300 parcels)
landcover <- data.frame(
  class = sample(
    c("Forest", "Shrubland", "Grassland", "Wetland", "Developed"),
    size = 300, replace = TRUE,
    prob = c(0.34, 0.22, 0.18, 0.12, 0.14)
  )
)

landcover |>
  count(class) |>
  mutate(class = forcats::fct_reorder(class, n, .desc = TRUE)) |>
  ggplot(aes(class, n, fill = class)) +
  geom_col(width = 0.6, colour = "#1e3a2f", alpha = 0.88) +
  geom_text(aes(label = n), vjust = -0.4,
            family = "mono", size = 3.5, colour = "#1e3a2f") +
  scale_fill_manual(values = c(
    "Forest"    = "#1e3a2f",
    "Shrubland" = "#4a7c5e",
    "Grassland" = "#8aad96",
    "Wetland"   = "#2e7d94",
    "Developed" = "#c8832a"
  )) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.12))) +
  labs(x = NULL, y = "Number of parcels",
       title = "Land-Cover Classification — 300 Parcels") +
  theme_minimal(base_size = 12) +
  theme(
    plot.title      = element_text(family = "serif", face = "bold", colour = "#1e3a2f"),
    legend.position = "none",
    panel.grid.major.x = element_blank()
  )

Land-cover class frequencies. Categories are sorted by count (most to least common) using fct_infreq(). A sorted bar chart lets the reader rank categories at a glance; alphabetical ordering forces the eye to search.

⚠️ The “bar chart with error bars” problem

A bar chart of means with error bars — sometimes called a “dynamite plot” — is one of the most overused and misleading plot types in ecology. The bar implies that data begin at zero, which focuses attention on the wrong feature (the mean height, not the spread). It also hides whether the distribution is symmetric, skewed, bimodal, or contains outliers. See Section 5 for better alternatives when the goal is to compare group means.

5. Continuous Variable by Group

This is the most common plotting task in environmental science: comparing a measurement across two or more groups. A family of plot types exists, each with different strengths.

set.seed(14)
# Macroinvertebrate EPT richness (sensitive taxa) at four stream sites
ept <- data.frame(
  site     = rep(c("Reference", "Urban", "Agricultural", "Mining"), each = 30),
  richness = c(
    round(rnorm(30, 12, 2.1)),   # Reference  — highest quality
    round(rnorm(30,  7, 2.8)),   # Urban
    round(rnorm(30,  9, 2.4)),   # Agricultural
    round(rnorm(30,  4, 1.9))    # Mining — most impaired
  )
) |>
  mutate(site = factor(site,
    levels = c("Reference", "Agricultural", "Urban", "Mining")))

Boxplot — the workhorse

Shows median, IQR, whiskers (1.5×IQR), and outliers. Always add the raw data as jittered points so the reader sees both distribution and sample size.

ggplot(ept, aes(site, richness, fill = site)) +
  geom_boxplot(width = 0.45, outlier.shape = NA,     # hide outlier points —
               colour = "#1e3a2f", alpha = 0.75,     # we show all raw data below
               linewidth = 0.6) +
  geom_jitter(aes(colour = site), width = 0.18,
              size = 1.6, alpha = 0.45) +
  scale_fill_manual(values = c(
    "Reference"    = "#1e3a2f",
    "Agricultural" = "#8aad96",
    "Urban"        = "#c8832a",
    "Mining"       = "#7a3a2f"
  )) +
  scale_colour_manual(values = c(
    "Reference"    = "#1e3a2f",
    "Agricultural" = "#4a7c5e",
    "Urban"        = "#c8832a",
    "Mining"       = "#7a3a2f"
  )) +
  labs(x = NULL, y = "EPT richness (no. of taxa)",
       title = "Macroinvertebrate EPT Richness by Stream Site") +
  theme_minimal(base_size = 12) +
  theme(
    plot.title         = element_text(family = "serif", face = "bold", colour = "#1e3a2f"),
    legend.position    = "none",
    panel.grid.major.x = element_blank()
  )

EPT richness across four stream sites. Boxes show the IQR; horizontal bars the median; jittered points show all individual observations. Colour is redundant with the x-axis position — it aids quick group identification.

Violin plot — shows the full distribution shape

When sample size is large enough (n ≥ 25 per group), a violin plot is superior to a boxplot because it reveals bimodality and skewness that the box hides.

ggplot(ept, aes(site, richness, fill = site)) +
  geom_violin(trim = TRUE, colour = "#1e3a2f",
              alpha = 0.7, linewidth = 0.5) +
  geom_boxplot(width = 0.08, fill = "white",
               colour = "#1e3a2f", outlier.shape = NA, linewidth = 0.7) +
  scale_fill_manual(values = c(
    "Reference"    = "#1e3a2f",
    "Agricultural" = "#8aad96",
    "Urban"        = "#c8832a",
    "Mining"       = "#7a3a2f"
  )) +
  labs(x = NULL, y = "EPT richness (no. of taxa)",
       title = "Violin + Boxplot: Full Distribution Shape") +
  theme_minimal(base_size = 12) +
  theme(
    plot.title         = element_text(family = "serif", face = "bold", colour = "#1e3a2f"),
    legend.position    = "none",
    panel.grid.major.x = element_blank()
  )

$Violin plots of EPT richness. The width of each violin at any given y-value is proportional to the density of observations at that value. The embedded boxplot preserves the median and IQR summary.$

Violin plots of EPT richness. The width of each violin at any given y-value is proportional to the density of observations at that value. The embedded boxplot preserves the median and IQR summary.

Point + CI plot — for inference

When your goal is comparison of means (i.e., you are about to run a t-test or ANOVA), show the mean and 95% confidence interval explicitly. This is more honest than a “dynamite plot” and more focused than a boxplot.

ept_summary <- ept |>
  group_by(site) |>
  summarise(
    n    = n(),
    mean = mean(richness),
    se   = sd(richness) / sqrt(n),
    t_crit = qt(0.975, df = n - 1),
    lower  = mean - t_crit * se,
    upper  = mean + t_crit * se
  )

ggplot(ept_summary, aes(site, mean, colour = site)) +
  geom_hline(yintercept = ept_summary$mean[ept_summary$site == "Reference"],
             linetype = "dashed", colour = "#aaaaaa", linewidth = 0.6) +
  geom_errorbar(aes(ymin = lower, ymax = upper),
                width = 0.15, linewidth = 1.1) +
  geom_point(size = 4) +
  scale_colour_manual(values = c(
    "Reference"    = "#1e3a2f",
    "Agricultural" = "#4a7c5e",
    "Urban"        = "#c8832a",
    "Mining"       = "#7a3a2f"
  )) +
  labs(x = NULL, y = "Mean EPT richness ± 95% CI",
       title = "Mean EPT Richness with 95% Confidence Intervals") +
  theme_minimal(base_size = 12) +
  theme(
    plot.title      = element_text(family = "serif", face = "bold", colour = "#1e3a2f"),
    legend.position = "none"
  )

Mean EPT richness ± 95% CI by site. Non-overlapping intervals between Reference and Mining provide strong visual evidence of a meaningful difference. This plot directly mirrors what a t-test or ANOVA tests.

📌 When to use each group-comparison plot

Boxplot: default choice for most ecological data; works well for n ≥ 10 per group; always layer raw data on top.
Violin: upgrade when n ≥ 25 and distribution shape matters — detecting bimodality in species abundance, for instance.
Point + CI: use when your primary message is inferential — “these groups are or are not significantly different”. Pairs naturally with ANOVA results.

6. Two Continuous Variables

Scatterplot with smoother

Use geom_smooth() to overlay a trend line. The default method "loess" is non-parametric and reveals non-linear patterns. Use method = "lm" when you are specifically asserting a linear relationship.

set.seed(21)
n <- 80
soil <- data.frame(
  organic_pct = runif(n, 1, 8),
  moisture    = NA,
  bulk_density = NA
)
soil$moisture     <- 12 + 4.5 * soil$organic_pct + rnorm(n, 0, 3)
soil$bulk_density <- 1.85 - 0.08 * soil$organic_pct + rnorm(n, 0, 0.08)
soil$land_use     <- sample(c("Forest", "Cropland"), n, replace = TRUE,
                             prob = c(0.5, 0.5))

ggplot(soil, aes(organic_pct, moisture, colour = land_use)) +
  geom_point(size = 2.2, alpha = 0.7) +
  geom_smooth(aes(group = 1), method  = "loess", span = 0.8,
              colour = "#1e3a2f", fill = "#d4e6da",
              linewidth = 1, se = TRUE) +
  geom_smooth(aes(group = 1), method = "lm", se = FALSE,
              colour = "#c8832a", linetype = "dashed",
              linewidth = 0.8) +
  scale_colour_manual(
    values = c("Forest" = "#4a7c5e", "Cropland" = "#c8832a"),
    name   = "Land use"
  ) +
  labs(
    x     = "Soil organic matter (%)",
    y     = "Volumetric moisture (%)",
    title = "Soil Moisture vs. Organic Matter"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title      = element_text(family = "serif", face = "bold", colour = "#1e3a2f"),
    legend.position = "top"
  )

Soil moisture increases with organic matter content. The loess smoother (blue) reveals a slightly non-linear relationship — the slope flattens at high organic matter. The linear fit (dashed) would miss this. Point colour encodes land use.

Small multiples (facets)

When you have a third grouping variable, facet_wrap() produces small multiples — the same axes and scales repeated across groups, which makes cross-group comparison easier than overlaying everything on one plot.

ggplot(soil, aes(organic_pct, moisture)) +
  geom_point(aes(colour = land_use), size = 1.9, alpha = 0.65) +
  geom_smooth(method = "lm", colour = "#1e3a2f",
              fill = "#d4e6da", linewidth = 1) +
  facet_wrap(~ land_use) +
  scale_colour_manual(values = c("Forest" = "#4a7c5e", "Cropland" = "#c8832a")) +
  labs(x = "Soil organic matter (%)", y = "Volumetric moisture (%)") +
  theme_minimal(base_size = 12) +
  theme(
    strip.text      = element_text(family = "mono", face = "bold",
                                   colour = "#1e3a2f", size = 10),
    legend.position = "none"
  )

The organic matter–moisture relationship faceted by land use. Keeping the axes identical across panels allows direct comparison of both slope and intercept. Cropland soils retain less moisture at a given organic matter level.

Correlation matrix for exploratory analysis

When you have several continuous variables and want a quick overview of their pairwise relationships, a pairs plot is invaluable.

pairs(
  soil[, c("organic_pct", "moisture", "bulk_density")],
  pch    = 20, cex = 0.7,
  col    = ifelse(soil$land_use == "Forest", "#4a7c5e", "#c8832a"),
  labels = c("Organic\nmatter (%)", "Moisture\n(%)", "Bulk\ndensity\n(g/cm³)"),
  panel  = function(x, y, ...) {
    points(x, y, ...)
    lines(lowess(x, y), col = "#1e3a2f", lwd = 1.5)
  },
  main = "Soil Variable Pairs"
)

Lower triangle: scatterplots with loess smoothers. Diagonal: variable name. Upper triangle would show correlation coefficients (see GGally::ggpairs() for a polished version). Bulk density is inversely related to both moisture and organic matter — classic compaction signal.

7. Time Series

Environmental data are often collected over time. Connecting observations with a line communicates temporal continuity; adding points marks the actual measurement occasions.

set.seed(33)
# Monthly mean streamflow (m³/s) — 5 years of a hypothetical CA watershed
months <- seq(as.Date("2019-01-01"), as.Date("2023-12-01"), by = "month")
n_mo   <- length(months)

# Seasonal signal + long-term decline + noise
seasonal   <- 12 * sin((seq_len(n_mo) - 2) * 2 * pi / 12) + 8
trend      <- seq(0, -4, length.out = n_mo)
streamflow <- data.frame(
  date  = months,
  flow  = pmax(0.2, seasonal + trend + rnorm(n_mo, 0, 2)),
  year  = format(months, "%Y")
)

ggplot(streamflow, aes(date, flow)) +
  geom_area(fill = "#d4e6da", colour = NA, alpha = 0.6) +
  geom_line(colour = "#2d5a3d", linewidth = 0.85) +
  geom_point(colour = "#1e3a2f", size = 1.2, alpha = 0.6) +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  labs(
    x     = NULL,
    y     = "Mean monthly streamflow (m³/s)",
    title = "Streamflow — Hypothetical CA Watershed, 2019–2023"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title      = element_text(family = "serif", face = "bold", colour = "#1e3a2f"),
    panel.grid.minor = element_blank()
  )

Monthly mean streamflow over five years. The line reveals the annual wet-season pulse; points mark individual monthly measurements. The declining trend is subtle but visible by comparing peak flows across years.

Seasonal decomposition plot

Faceting by year allows direct visual comparison of the seasonal cycle across years — much clearer than trying to read inter-annual differences from the single time series above.

streamflow |>
  mutate(month_num = as.integer(format(date, "%m")),
         month_lbl = format(date, "%b")) |>
  mutate(month_lbl = factor(month_lbl, levels = month.abb)) |>
  ggplot(aes(month_lbl, flow, group = year, colour = year)) +
  geom_line(linewidth = 1) +
  geom_point(size = 1.8) +
  scale_colour_manual(
    values = c("2019" = "#1e3a2f", "2020" = "#2d5a3d",
               "2021" = "#4a7c5e", "2022" = "#8aad96", "2023" = "#c8832a"),
    name = "Year"
  ) +
  labs(x = NULL, y = "Streamflow (m³/s)",
       title = "Monthly Streamflow by Year") +
  theme_minimal(base_size = 11) +
  theme(
    plot.title      = element_text(family = "serif", face = "bold", colour = "#1e3a2f"),
    axis.text.x     = element_text(size = 8),
    legend.position = "top"
  )

Streamflow by month, faceted by year. Each panel shares the same x-axis (month of year) and y-axis scale, making it easy to compare peak-flow timing and magnitude across years. The declining late-season flow is especially clear.

8. Design Principles

Technical correctness is necessary but not sufficient. A plot also needs to communicate clearly. These principles have empirical support from perceptual research.

Colour

Colour is the most misused visual channel. Three rules cover most situations:

✅ Do Use a qualitative palette (e.g., ColorBrewer Set2) for nominal categories where no ordering exists. Use a sequential palette (dark→light) for ordered or quantitative data. Use diverging palettes (two hues meeting at a neutral midpoint) for data centred on a meaningful zero, like anomalies from a long-term average. Test your palette with a colour-blindness simulator.

❌ Don’t Use the default rainbow palette — it has no perceptual order and is uninterpretable to red-green colourblind readers (~8% of men). Never use colour alone to encode a data dimension; also vary shape or line type so the plot is legible in black-and-white print and to colourblind readers. Do not use more than ~6–7 distinct colours in one plot.

p_bad <- ggplot(ept, aes(site, richness, fill = site)) +
  geom_boxplot(outlier.shape = NA, linewidth = 0.5) +
  scale_fill_manual(values = c("red", "orange", "green", "blue")) +
  labs(title = "❌ Rainbow colours", x = NULL, y = "EPT richness") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "none",
        plot.title = element_text(colour = "#bf360c", face = "bold"))

p_good <- ggplot(ept, aes(site, richness, fill = site)) +
  geom_boxplot(outlier.shape = NA, linewidth = 0.5) +
  scale_fill_brewer(palette = "Dark2") +
  labs(title = "✅ ColorBrewer 'Dark2'", x = NULL, y = "EPT richness") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "none",
        plot.title = element_text(colour = "#1b5e20", face = "bold"))

gridExtra::grid.arrange(p_bad, p_good, ncol = 2)

Left: rainbow palette — perceptually non-monotone, inaccessible. Right: ColorBrewer ‘Dark2’ qualitative palette — balanced, colourblind-friendly. Both show the same EPT richness data.

Axis scales and the zero baseline

p_trunc <- ggplot(ept_summary, aes(site, mean, colour = site)) +
  geom_errorbar(aes(ymin = lower, ymax = upper),
                width = 0.15, linewidth = 1) +
  geom_point(size = 3.5) +
  coord_cartesian(ylim = c(3, 14)) +
  scale_colour_manual(values = c("#1e3a2f","#4a7c5e","#c8832a","#7a3a2f")) +
  labs(title = "❌ Truncated y-axis", x = NULL, y = "Mean EPT richness") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "none",
        plot.title = element_text(colour = "#bf360c", face = "bold"))

p_zero <- ggplot(ept_summary, aes(site, mean, colour = site)) +
  geom_errorbar(aes(ymin = lower, ymax = upper),
                width = 0.15, linewidth = 1) +
  geom_point(size = 3.5) +
  expand_limits(y = 0) +
  scale_colour_manual(values = c("#1e3a2f","#4a7c5e","#c8832a","#7a3a2f")) +
  labs(title = "✅ Axis from zero", x = NULL, y = "Mean EPT richness") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "none",
        plot.title = element_text(colour = "#1b5e20", face = "bold"))

gridExtra::grid.arrange(p_trunc, p_zero, ncol = 2)

Left: y-axis truncated at 3 — the Urban–Mining difference looks enormous. Right: y-axis starting at 0 shows the true proportion of the range occupied by the difference. Truncating axes exaggerates effects and can mislead readers.

⚠️ The zero-baseline exception

The “start at zero” rule applies to bar charts and area charts, where the encoded quantity is the bar’s height or the filled area. For point-and-line plots and scatterplots, forcing a zero baseline often wastes most of the plot area and obscures variation that is scientifically meaningful. Judgment is required; the key is that the axis range should not exaggerate or minimise the effect relative to what is natural for the variable.

The data-ink ratio

Edward Tufte’s principle: every drop of ink should encode data. Remove gridlines, borders, and tick marks that do not help the reader.

base_plot <- ggplot(soil, aes(organic_pct, moisture)) +
  geom_point(colour = "#4a7c5e", alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE, colour = "#1e3a2f") +
  labs(x = "Organic matter (%)", y = "Moisture (%)")

p_default <- base_plot +
  labs(title = "❌ Default grey theme") +
  theme_grey(base_size = 11) +
  theme(plot.title = element_text(colour = "#bf360c", face = "bold"))

p_clean <- base_plot +
  labs(title = "✅ Minimal theme, refined") +
  theme_minimal(base_size = 11) +
  theme(
    plot.title       = element_text(family = "serif", face = "bold",
                                    colour = "#1b5e20"),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_line(colour = "#ede8df", linewidth = 0.5),
    axis.line        = element_line(colour = "#3a3a38", linewidth = 0.4)
  )

gridExtra::grid.arrange(p_default, p_clean, ncol = 2)

Left: default ggplot2 grey theme — heavy grid, bold background, double-encoded axis. Right: theme_minimal() with further refinement — the same data, more signal, less noise.

Labelling and annotation

Axis labels must include units. Titles should state the finding, not just the variable names. Annotations (annotate(), geom_text()) can replace or supplement a legend, reducing eye travel.

soil_means <- soil |>
  group_by(land_use) |>
  summarise(x = mean(organic_pct), y = mean(moisture))

ggplot(soil, aes(organic_pct, moisture, colour = land_use)) +
  geom_point(size = 2, alpha = 0.55) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 1) +
  geom_text(data = soil_means,
            aes(x = x + 0.3, y = y + 2, label = land_use, colour = land_use),
            fontface = "bold", family = "mono", size = 3.8) +
  scale_colour_manual(values = c("Forest" = "#1e3a2f", "Cropland" = "#c8832a")) +
  labs(
    x     = "Soil organic matter (%)",
    y     = "Volumetric moisture (%)",
    title = "Moisture increases with organic matter — more strongly in forest soils"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title      = element_text(family = "serif", face = "bold",
                                   colour = "#1e3a2f", size = 11),
    legend.position = "none"
  )

Direct annotation replaces the legend entirely — the reader’s eye never has to leave the plot. The title states the finding (‘Moisture increases with organic matter’) rather than just naming the axes.

9. Choosing the Right Plot

Plot selection guide. The best choice depends on your data structure and the question you want to answer visually.
Goal	Best choice	Avoid
Understand distribution of one variable	Histogram (`geom_histogram()`)	Pie chart; bars of means
Understand distribution of one variable	Density plot (`geom_density()`)	Histogram with too few or too many bins
Count observations per category	Bar chart, sorted (`geom_col()` + `fct_reorder()`)	Pie chart; 3-D bar chart
Compare distributions across groups	Boxplot + jitter (`geom_boxplot()` + `geom_jitter()`)	Bar chart of means only (‘dynamite plot’)
Compare distributions across groups (large n)	Violin + embedded boxplot (`geom_violin()` + `geom_boxplot()`)	Single mean ± SD bar
Compare means across groups (inference)	Point + CI (`geom_point()` + `geom_errorbar()`)	Bar chart of means without CI
Show relationship between two continuous variables	Scatterplot + smoother (`geom_point()` + `geom_smooth()`)	Line chart (implies continuity between discrete points)
Show three variables simultaneously	Scatterplot, colour/size encodes 3rd variable	3-D scatterplot (almost never interpretable)
Show change over time	Line plot (`geom_line()` + `scale_x_date()`)	Bar chart; connecting non-temporal measurements
Compare many variables pairwise	Pairs plot (`pairs()` or `GGally::ggpairs()`)	Correlation heatmap without scatterplot check

🔑 Key visualization takeaways

Plot first, always. Summary statistics can be identical across structurally different datasets (Anscombe’s Quartet). A histogram or scatterplot takes 2 lines of code and catches errors that hours of modelling will not.
The grammar of graphics is data + aesthetics + geometry + optional refinements. Build plots layer by layer; each component is independent.
Show the raw data whenever n is small enough. Overlay jittered points on boxplots; show individual dots on bar charts. Summary statistics that hide data are dishonest.
Use 95% CI bars, not SD bars, when making inferential comparisons. Label error bars explicitly — readers cannot tell SD from SE from CI by looking.
Colour for communication, not decoration. Use colourblind-safe palettes, never more than ~6 categories, and always a non-colour backup encoding (shape, linetype).
Titles state findings (“Moisture increases with organic matter”), not variable names (“Organic matter vs. moisture”). Your reader should understand the main message without reading the caption.

What’s next?

Describing Data — means, SD, SE, and CIs
Four Frameworks — choosing the right statistical test
t-test in depth — worked example with course data

✏️ Practice Exercises

Exercise 1 — Basic ggplot2

Using airquality, create a scatter plot of Temp (x-axis) vs Ozone (y-axis). Colour points by Month (as a factor). Add a smoothed trend line with geom_smooth(). Give the plot a descriptive title.

Show solution

library(ggplot2)
data(airquality)
ggplot(airquality, aes(x = Temp, y = Ozone, colour = factor(Month))) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "loess", colour = "black", se = TRUE) +
  labs(title = "Temperature vs. Ozone, New York 1973",
       x = "Temperature (°F)", y = "Ozone (ppb)",
       colour = "Month") +
  theme_minimal()

Exercise 2 — Faceting

Extend your scatter plot by faceting on Month using facet_wrap(). Does the relationship between temperature and ozone look the same in each month? Note one month that looks different and speculate why.

Show solution

ggplot(airquality, aes(x = Temp, y = Ozone)) +
  geom_point(alpha = 0.6, colour = "#4a7c5e") +
  geom_smooth(method = "lm", se = FALSE, colour = "#c8832a") +
  facet_wrap(~ Month, labeller = label_both) +
  labs(x = "Temperature (°F)", y = "Ozone (ppb)") +
  theme_minimal()

Month 9 (September) often shows a weaker relationship — lower ozone overall as summer photochemical activity declines, even at high temperatures.

Exercise 3 — Distribution comparison

Create side-by-side boxplots of Ozone by Month, then overlay the raw data points (jittered). Which month has the highest median ozone? Which has the greatest spread?

Show solution

ggplot(airquality, aes(x = factor(Month), y = Ozone, fill = factor(Month))) +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) +
  geom_jitter(width = 0.15, alpha = 0.5, size = 1.5) +
  scale_fill_brewer(palette = "Greens") +
  labs(x = "Month", y = "Ozone (ppb)", fill = "Month") +
  theme_minimal() + theme(legend.position = "none")

Descriptive Statistics ↑ Module Index Four Frameworks

Data Visualization

Marc Los Huertos

EA030 · Compiled 2026-04-14

Data Visualization

1. Why Visualize Before You Model?

2. The Grammar of Graphics

The seven layers

Minimal working example

3. One Continuous Variable

Histogram

Density plot

4. One Categorical Variable

5. Continuous Variable by Group

Boxplot — the workhorse

Violin plot — shows the full distribution shape

Point + CI plot — for inference

6. Two Continuous Variables

Scatterplot with smoother

Small multiples (facets)

Correlation matrix for exploratory analysis

7. Time Series

Seasonal decomposition plot

8. Design Principles

Colour

Axis scales and the zero baseline

The data-ink ratio

Labelling and annotation

9. Choosing the Right Plot

🔑 Key visualization takeaways

What’s next?