Figures for publication

Honest, accessible and ready for print

Dr Massimiliano Canzi

Figures for publication

Did you know: the arctic fox has the most insulating fur measured in any mammal. Its lower critical temperature is below −40 °C: above that, it keeps warm without raising its metabolic rate.

Source: Prestrud (1991), Arctic

Where we are

You already know how to build a plot in ggplot2, and the details are one search away. This lecture is about the decisions that no function makes for you: what to show, how to show it honestly, and how to get it onto the page of a thesis or a journal.

A figure is part of the evidence. Readers should be able to judge it for themselves, including those who print in greyscale or see colours differently.

Today

  1. Show the data
  2. Clutter and colour
  3. Models and uncertainty
  4. Publication-ready

The data

The lexical decision data from lecture 1: 30 participants (15 L1, 15 L2 speakers of English), 40 words each.

set.seed(2027)
participants <- tibble(participant = sprintf("P%02d", 1:30),
                       group = rep(c("L1", "L2"), each = 15),
                       p0 = rnorm(30, 0, 60),
                       p1 = rnorm(30, 0, 12))
words <- tibble(word = sprintf("w%02d", 1:40),
                log_freq = runif(40, 1, 6),
                length = pmin(pmax(round(11 - 1.1 * log_freq + rnorm(40, 0, 1.6)), 3), 13),
                w0 = rnorm(40, 0, 30))
lex <- crossing(participants, words) |>
  mutate(slope = -25 + p1 + if_else(group == "L2", -25, 0),
         rt = round(700 + p0 + w0 + if_else(group == "L2", 120, 0) +
                    slope * (log_freq - 3.5) + 12 * (length - 7) + rnorm(n(), 0, 70))) |>
  select(participant, group, word, log_freq, length, rt)

Show the data

Did you know: an arctic fox’s winter fur is nearly twice as long as its summer fur.

Source: Underwood & Reynolds (1980), report for the Office of Naval Research

What a reader needs

To judge a difference between groups, a reader needs the size of the difference, the variability in the data and the sample size.

A bar chart of group means shows only the size of the difference.

The dynamite plot

lex |>
  summarise(mean = mean(rt), se = sd(rt) / sqrt(n()), .by = group) |>
  ggplot(aes(x = group, y = mean)) +
  geom_col(fill = pandan, width = 0.6) +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.15, colour = mist)

A dynamite plot shows two numbers for 1,200 responses from 30 people. The bars suggest values from zero up to the mean, and the SE treats all 600 responses per group as independent.

What bars hide

Three sentence types with the same mean acceptability: one bell shape, two clusters of speakers, and only four ratings.

Which points?

p_trials <- ggplot(lex, aes(group, rt)) +
  geom_jitter(width = 0.15, alpha = 0.1, colour = kingfisher) + labs(title = "Trials")
p_people <- lex |>
  summarise(rt = mean(rt), .by = c(participant, group)) |>
  ggplot(aes(group, rt)) +
  geom_jitter(width = 0.15, colour = kingfisher, size = 2.5) + labs(title = "Participants")
p_trials + p_people

Plot the unit of analysis your claim is about. A difference between groups of speakers is a claim about participants: show participant means.

Show the data

participant_means <- summarise(lex, rt = mean(rt), .by = c(participant, group))
ggplot(participant_means, aes(x = group, y = rt)) +
  geom_jitter(width = 0.1, height = 0, colour = muted, size = 2.5) +
  stat_summary(fun.data = mean_se, fun.args = list(mult = 1.96),
               colour = kingfisher, size = 0.9, linewidth = 1.3,
               position = position_nudge(x = 0.25)) +
  labs(x = NULL, y = "Mean RT per participant (ms)")

Every participant, plus the group mean and its 95% confidence interval, moved to the side so it hides no points.

Which error bar?

  • The SD describes the spread of the data; the SE and the CI describe the precision of the mean
  • The same 15 participants give three very different bars: always say in the caption which one you show

Distributions: raincloud

ggplot(lex, aes(x = rt, y = group, fill = group)) +
  geom_density_ridges(jittered_points = TRUE, position = "raincloud", colour = mist,
                      alpha = 0.6, scale = 0.7, point_colour = mist, point_alpha = 0.2, point_size = 0.7,
                      quantile_lines = TRUE, quantiles = 2) +
  labs(x = "RT (ms)", y = NULL) +
  theme(legend.position = "none")

A raincloud plot from {ggridges} shows the shape of the distribution above, the observations below, and a line at the median.

Axes

  • A bar shows a value by its length, so its axis must start at zero
  • Points and lines show values by position, so their axis can zoom in on the range that matters

Free scales in facets and log scales are fine as long as the axis title or caption says so.

Say what you did

In the caption or axis titles, tell the reader:

  • What the error bars or ribbons are: SD, SE or a confidence interval, and for which unit
  • Whether points are raw data, participant means or model predictions
  • How many observations each group has
  • Whether an axis is on a log scale or does not start at zero
  • Whether points are jittered or lines smoothed
  • Whether any data were excluded, how many, and why

Clutter and colour

Did you know: arctic foxes come in two genetically fixed colour phases. The white phase is white only in winter; the blue phase is dark all year.

Sources: Prestrud (1991), Arctic; Norwegian Polar Institute

Before and after

What changed

  • The legend repeated the x axis, so it went
  • Red and green separated nothing the axis did not already separate, and fail for many readers
  • Axis titles became words with units; x needs no title when the labels say it all
  • The grey background, heavy grid and outlines added ink but no information
  • We kept the points, because they show the data

Colour for every reader

  • Choose a qualitative, sequential or diverging palette from the type of variable: categories, low to high, or distance from a meaningful midpoint
  • Red–green colour-vision deficiencies are common, especially among men, and figures are printed in greyscale and shown on dim projectors
  • The Okabe–Ito palette keeps eight colours distinct for the common deficiencies, and is built into R 4.0 or later
okabe <- unname(palette.colors(palette = "Okabe-Ito"))
group_colours <- c(L1 = okabe[2], L2 = okabe[3])

Redundant coding

word_group_means <- summarise(lex, rt = mean(rt), .by = c(word, log_freq, group))
ggplot(word_group_means, aes(x = log_freq, y = rt, fill = group, shape = group)) +
  geom_point(size = 3, colour = floor_col, stroke = 0.6) +
  scale_shape_manual(values = c(L1 = 21, L2 = 24)) +
  scale_fill_manual(values = group_colours) +
  labs(x = "log10 frequency", y = "Mean RT (ms)", fill = NULL, shape = NULL) +
  theme(legend.position = "right")

Colour and shape show the same variable (redundant coding), so the groups survive greyscale. The named vector keeps each group’s colour across all figures.

Order categories

word_means <- summarise(lex, rt = mean(rt), .by = c(word, log_freq))
twelve <- filter(word_means, word %in% sprintf("w%02d", 1:12))
p_alpha <- ggplot(twelve, aes(x = rt, y = word)) + geom_point(colour = kingfisher, size = 3)
p_ordered <- ggplot(twelve, aes(x = rt, y = fct_reorder(word, rt))) +
  geom_point(colour = kingfisher, size = 3) + labs(y = "word")
p_alpha + p_ordered & labs(x = "Mean RT (ms)")

Alphabetical order is rarely the order that matters. fct_reorder() sorts the categories by a value, here the mean RT.

Spaghetti

ggplot(lex, aes(x = log_freq, y = rt, group = participant, colour = group)) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 0.6) +
  labs(x = "log10 frequency", y = "Fitted RT (ms)", colour = NULL)

In this spaghetti plot, thirty overlapping lines make it impossible to follow any one participant.

Small multiples with context

some <- filter(lex, participant %in% c("P01", "P02", "P16", "P17"))
others <- rename(lex, id = participant)
ggplot(some, aes(x = log_freq, y = rt)) +
  geom_smooth(data = others, aes(group = id), method = "lm", se = FALSE,
              colour = vine, linewidth = 0.4) +
  geom_smooth(aes(colour = group), method = "lm", se = FALSE, linewidth = 1.4) +
  facet_wrap(~ participant, nrow = 1) +
  labs(x = "log10 frequency", y = "Fitted RT (ms)", colour = NULL)

In these small multiples, the grey layer has no participant column, so it appears in every panel.

Things to avoid

  • Pie charts ask readers to compare angles and areas, which they do badly; bars or points work better
  • With two different y axes on one plot, the choice of scaling decides what looks related
  • 3D effects, shadows and textures add ink without information
  • Rainbow palettes for continuous values
  • Bars that do not start at zero, and bars for means without the data

Models and uncertainty

Did you know: arctic foxes prepare for winter food shortages by storing fat and caching food during summer and autumn.

Source: Prestrud (1991), Arctic

Raw data or model?

  • Raw data show what was observed, including what the model leaves out
  • A model shows the estimated effects, adjusted for the other predictors, with their uncertainty
  • A good figure often shows both: the model as lines and ribbons, the data behind it

The ribbon of geom_smooth()

ggplot(lex, aes(x = log_freq, y = rt, colour = group, fill = group)) +
  geom_point(alpha = 0.05) +
  geom_smooth(method = "lm") +
  labs(x = "log10 frequency", y = "RT (ms)", colour = NULL, fill = NULL)

This ribbon comes from a separate lm() for each group that treats all 600 responses as independent. It ignores the repeated measures, so it is too narrow.

The model

lex <- lex |>
  mutate(freq_c = log_freq - mean(log_freq),
         length_c = length - mean(length),
         group = factor(group))
contrasts(lex$group) <- contr.sum(2)
m4 <- lmer(rt ~ freq_c * group + length_c + (1 + freq_c | participant) + (1 | word),
           data = lex)

The same model as in lecture 1, now to be shown in a figure.

Predictions with {ggeffects}

pred <- ggpredict(m4, terms = c("freq_c", "group")) |>
  as.data.frame() |>
  as_tibble() |>
  mutate(log_freq = x + mean(lex$log_freq))
pred
# A tibble: 22 × 7
      x predicted std.error conf.low conf.high group log_freq
  <dbl>     <dbl>     <dbl>    <dbl>     <dbl> <fct>    <dbl>
1  -2.5      735.      22.9     690.      779. L1       0.910
2  -2.5      975.      22.9     930.     1020. L2       0.910
3  -2        726.      21.1     685.      767. L1       1.41 
4  -2        950.      21.1     908.      991. L2       1.41 
# ℹ 18 more rows

Predicted RTs for each group across the frequency range, with length_c held at its mean (0), and 95% confidence intervals from the fixed effects.

Your own prediction plot

ggplot(pred, aes(x = log_freq, y = predicted, colour = group, fill = group)) +
  geom_point(data = word_group_means, aes(y = rt), alpha = 0.4, size = 1.5) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.25, colour = NA) +
  geom_line(linewidth = 1.2) +
  labs(x = "log10 frequency", y = "RT (ms)", colour = NULL, fill = NULL,
       caption = "Lines and ribbons: mixed-model predictions with 95% CIs. Points: word means.") +
  theme(legend.position = "right")

ggpredict() returns a data frame, so the plot is an ordinary ggplot() that you can style like any other.

A coefficient plot

coef(summary(m4)) |>
  as_tibble(rownames = "term") |>
  filter(term != "(Intercept)") |>
  mutate(low = Estimate - 1.96 * `Std. Error`, high = Estimate + 1.96 * `Std. Error`) |>
  ggplot(aes(x = Estimate, y = fct_rev(fct_inorder(term)))) +
  geom_vline(xintercept = 0, colour = muted, linetype = "dashed") +
  geom_pointrange(aes(xmin = low, xmax = high), colour = kingfisher, size = 0.8, linewidth = 1.2) +
  labs(x = "Estimate (ms) with 95% Wald interval", y = NULL)

A coefficient plot has one row per fixed effect. It works best when there are many predictors on the same scale.

Which plot for a model?

Predictions show effects in the units of the outcome, which makes them the better choice for interactions and non-linear effects. Coefficients are compact, but each estimate depends on coding and centering, so readers also need the model specification.

In both cases, say which predictors are held constant, and at which values.

Publication-ready

Did you know: among arctic mammals up to the size of an arctic fox, bigger species have longer winter fur. Larger species, such as wolves and caribou, all have about the same fur length.

Source: Scholander et al. (1950), as summarised by Underwood & Reynolds (1980)

Titles that say something

labs(title = "Frequency speeds up L2 speakers more than L1 speakers",
     subtitle = "Mixed-model predictions with 95% confidence intervals",
     caption = "30 participants, 40 words. Points: word means.")
  • On slides and posters, a title that states the finding tells the audience what to look for
  • In papers, the journal’s caption usually takes that role: put the finding in the first sentence of the caption

Reference lines and notes

ggplot(participant_means, aes(x = group, y = rt)) +
  geom_hline(yintercept = mean(participant_means$rt), colour = muted, linetype = "dashed") +
  geom_jitter(width = 0.1, height = 0, colour = kingfisher, size = 2.5) +
  annotate("text", x = 2.45, y = mean(participant_means$rt), label = "grand mean",
           colour = muted, vjust = -0.5, size = 5) +
  scale_x_discrete(expand = expansion(add = c(0.6, 1))) +
  labs(x = NULL, y = "Mean RT per participant (ms)")

annotate() adds one element at fixed coordinates, without a data frame.

Labels that move aside

slowest <- slice_max(word_means, rt, n = 5)
ggplot(word_means, aes(x = log_freq, y = rt)) +
  geom_point(colour = kingfisher, size = 2.5) +
  geom_text_repel(data = slowest, aes(label = word), colour = mist, size = 5,
                  min.segment.length = 0, box.padding = 0.6, seed = 1) +
  labs(x = "log10 frequency", y = "Mean RT per word (ms)")

{ggrepel} places labels so they do not overlap each other or the points. Label only the points the text talks about.

Combining figures

p_pred <- ggplot(pred, aes(x = log_freq, y = predicted, colour = group, fill = group)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.25, colour = NA) +
  geom_line(linewidth = 1.2) + labs(x = "log10 frequency", y = "Predicted RT (ms)") +
  theme(legend.position = "right")
p_raw <- ggplot(participant_means, aes(x = group, y = rt, colour = group)) +
  geom_jitter(width = 0.1, height = 0, size = 2.5) + labs(x = NULL, y = "Mean RT (ms)") +
  theme(legend.position = "none")
(p_raw | p_pred) + plot_layout(widths = c(1, 2)) +
  plot_annotation(tag_levels = "A") & labs(colour = NULL, fill = NULL)

{patchwork}: | side by side, / on top, & for every panel. plot_layout(guides = "collect") merges identical legends.

A theme for the paper

theme_paper <- function(base_size = 11, base_family = "") {
  theme_minimal(base_size = base_size, base_family = base_family) +
    theme(plot.background = element_rect(fill = "white", colour = NA),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "grey90"),
          axis.title = element_text(colour = "grey30"),
          strip.text = element_text(face = "bold", hjust = 0),
          plot.title.position = "plot",
          legend.position = "bottom")
}

Write the theme function once, keep it in a file, and source() it at the top of every script for the same paper. theme_set(theme_paper()) makes it the default.

The theme in use

p_pred +
  scale_colour_manual(values = group_colours) +
  scale_fill_manual(values = group_colours) +
  labs(colour = NULL, fill = NULL) +
  theme_paper(base_size = 16) +
  theme(legend.position = "right")

Panel B again, with the paper’s theme and the Okabe–Ito colours.

Size on the page

A figure made 17 cm wide with 11 pt text, then shrunk or stretched to fit: the text shrinks or grows with it. Make the figure at its final size.

Exporting with ggsave()

fig1 <- (p_raw | p_pred) + plot_annotation(tag_levels = "A")
ggsave("figures/figure1.pdf", fig1, width = 17, height = 8, units = "cm")
ggsave("figures/figure1.png", fig1, width = 17, height = 8, units = "cm", dpi = 300)
  • Set width and height to the size on the page: often about 8.5 cm for one column and 17 cm for the full width; check the journal’s guidelines
  • Use the same base_size in your theme for every figure, so all text matches the body text

Vector or raster?

Vector (.pdf, .svg) Raster (.png, .tiff)
Made of shapes and text pixels
Zoom sharp at any size blurry when enlarged
Best for most plots plots with very many points; some journals
Size setting width and height width, height and dpi (300 for print)

Never paste a screenshot of the RStudio plot pane into a thesis or a paper.

Before you submit

  • Does each figure answer one clear question?
  • Can the reader see the data, the variability and the sample size?
  • Do the caption and axes say what the points, error bars and ribbons are?
  • Are axes labelled in words, with units?
  • Does it work in greyscale and for colour-blind readers?
  • Is the text readable at the final size?
  • Is the code saved, so you can redraw the figure when a reviewer asks?

Summary

  • Show the data, at the level of the unit your claim is about, with variability and sample size
  • Say what error bars and ribbons are; SD, SE and CI answer different questions
  • Remove what carries no information; order categories by something meaningful; use small multiples instead of spaghetti
  • Choose palettes that work for colour-blind readers and in greyscale, and code key groups twice
  • Plot model predictions with their uncertainty over the raw data; geom_smooth() ignores repeated measures
  • Make figures at their final size, in a vector format, with one theme function for the whole paper

References

  • Allen, M., Poggiali, D., Whitaker, K., Marshall, T. R., van Langen, J., & Kievit, R. A. (2019). Raincloud plots: A multi-platform tool for robust data visualization. Wellcome Open Research, 4, 63.
  • Lüdecke, D. (2018). ggeffects: Tidy data frames of marginal effects from regression models. Journal of Open Source Software, 3(26), 772.
  • Okabe, M., & Ito, K. (2008). Color universal design (CUD): How to make figures and presentations that are friendly to colorblind people. https://jfly.uni-koeln.de/color/
  • Tufte, E. R. (1983). The visual display of quantitative information. Graphics Press.
  • Weissgerber, T. L., Milic, N. M., Winham, S. J., & Garovic, V. D. (2015). Beyond bar and line graphs: Time for a new data presentation paradigm. PLoS Biology, 13(4), e1002128.
  • Wilke, C. O. (2019). Fundamentals of data visualization. O’Reilly. https://clauswilke.com/dataviz/

Sources of the animal facts

  • Prestrud, P. (1991). Adaptations by the arctic fox (Alopex lagopus) to the polar winter. Arctic, 44(2), 132–138.
  • Scholander, P. F., Walters, V., Hock, R., & Irving, L. (1950). Body insulation of some arctic and tropical mammals and birds. Biological Bulletin, 99(2), 225–236.
  • Norwegian Polar Institute. Arctic fox (Vulpes lagopus). https://npolar.no/?p=5053
  • Underwood, L. S., & Reynolds, P. (1980). Photoperiod and fur lengths in the arctic fox (Alopex lagopus L.). Report to the Office of Naval Research, University of Alaska Anchorage. https://apps.dtic.mil/sti/tr/pdf/ADA085976.pdf

Slides and code

QR code linking to these slides

These slides

mcanzi.de/ss27_advanced_statistics/02_figures.html

Code and all lectures

github.com/mcanzi/ss27_advanced_statistics