ETC5521 Worksheet Week 4

Assessing significance of patterns

Author

Prof. Di Cook

Packages and preparation

This worksheet requires R version 4.5.0 or later (so that penguins is available without any extra package) plus three packages. Install them once from CRAN if you have not already done so:

Code 
install.packages(c("tidyverse", "nullabor", "MASS"))

Then load them at the start of your session:

Code 
library(tidyverse)   # ggplot2, dplyr, and friends
library(nullabor)    # lineup(), rorschach(), pvisual(), null_permute(), null_lm()
library(MASS)        # mvrnorm() — multivariate normal simulation (Exercise 3 only)
Package Role in this worksheet Version used
tidyverse [@Wickham2019] Data wrangling (dplyr, tidyr) and plotting (ggplot2) ≥ 2.0
nullabor [@Buja2009] Lineup and Rorschach protocols, null-generating functions, pvisual() ≥ 0.3.10
MASS [@Venables2002] mvrnorm() for bivariate normal simulation in Exercise 3 ≥ 7.3

The penguins dataset is part of the datasets package that ships with base R (≥ 4.5.0) [@Kaye2025], so no additional install is needed for the data.

“The most important maxim for data analysis … is that it is better to look at the data than to stare at the output of a statistical test.” — John W. Tukey

A lineup is a way of making this idea rigorous. You hide the real data among a set of plots generated under the assumption that nothing interesting is happening (the null hypothesis), then ask an observer — you, your classmates, or anyone — to identify which plot is real. If the real data is easy to find, that is evidence something interesting is happening. If it blends in with the rest, the data are consistent with the null hypothesis.

Before we put real data in a lineup, we need to know what “nothing interesting” looks like. That is what the Rorschach protocol is for: it shows a set of plots that are all generated under the null, with no real data at all. This calibrates your eye.

This worksheet introduces both protocols using the penguins dataset, which has been part of base R since version 4.5.0. No extra packages are needed to load the data.

R version note. The base R penguins dataset uses shorter variable names than the palmerpenguins package version. The names used in this worksheet are: bill_len, bill_dep, flipper_len, and body_mass (without unit suffixes).

🎯 Learning Objectives

By the end of this worksheet you should be able to:

  1. Explain the purpose of the Rorschach protocol and use it to calibrate your eye to natural random variation.
  2. Explain the purpose of the lineup protocol and use it to decide whether a pattern in data is surprising.
  3. Choose an appropriate null-generating mechanism (null_permute, null_dist, or null_lm) given a description of the plot and the hypothesis being tested.
  4. Construct a lineup using nullabor::lineup() and display it with facet_wrap().
  5. Compute a visual p-value with pvisual() and state a conclusion.

Setup: meet the data

Code 
# The data loads automatically in R >= 4.5.0
glimpse(penguins)
Rows: 344
Columns: 8
$ species     <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Ad…
$ island      <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgersen, Tor…
$ bill_len    <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, 42.0, …
$ bill_dep    <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, 20.2, …
$ flipper_len <int> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186, 180,…
$ body_mass   <int> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, 4250, …
$ sex         <fct> male, female, female, NA, female, male, female, male, NA, …
$ year        <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007…

There are three species of penguin — Adélie, Chinstrap, and Gentoo — measured on four physical variables (bill_len, bill_dep, flipper_len, body_mass) along with sex, island, and year.

Exercise 1 — The Rorschach: what does nothing look like?

Background

Before you can judge whether a pattern in your data is real, you need to know what patterns arise by chance alone. The Rorschach protocol generates several plots that are all pure noise — no real data anywhere. Looking at them trains your eye before you see a lineup.

The question for this exercise is: do the three penguin species have different flipper lengths? The plot we will use is a boxplot of flipper_len by species.

The null hypothesis is: species labels are exchangeable — knowing which species a penguin belongs to tells you nothing about its flipper length. We generate null data by randomly shuffling (permuting) the species column, so each penguin keeps its actual flipper length but gets a randomly assigned species label.

The null-generating function

null_permute("species") generates one dataset where the species labels have been randomly shuffled. We call lineup() in Rorschach mode (no real data) to get 20 such shuffled datasets.

a. In every panel, the species labels are randomly shuffled, so any difference between the three boxes is pure chance. Do the three boxes look roughly the same height in most panels, or do they vary quite a bit?

b. Find a panel where one species looks noticeably higher or lower than the others. Write down the panel number. Would you have been suspicious of that panel if you had seen it alone, without the other 19?

c. Now think about what the real data would have to look like to stand out from these panels. Describe in one sentence what you would need to see in the real data plot for it to be clearly distinguishable from the null panels.

Exercise 2 — The lineup: do species differ in flipper length?

Background

Now we embed the real penguin data among 19 null panels generated by the same permutation mechanism. This is the lineup protocol.

a. Look carefully at all 20 panels. Write down which panel you think contains the real data and the specific feature that identified it.

b. The null hypothesis being tested is: species labels are exchangeable — flipper length does not differ between species. What does it mean for the real panel to stand out clearly from the null panels, in terms of this hypothesis?

c. Decrypt the lineup to reveal the true position. Were you correct?

d. Suppose this lineup was shown to 16 students and 14 of them correctly identified the real panel. Compute the visual p-value and state your conclusion. The lineup has \(m = 20\) panels.

e. Describe what the p-value means in plain language, without using the word “probability”. Use the lineup panels to explain your answer.

Exercise 3 — Is bill length and depth jointly normal for Adélie penguins?

Background

When we plot two continuous variables against each other in a scatter plot, we are often implicitly asking whether their joint distribution follows a particular shape. A bivariate normal distribution is the two-dimensional generalisation of the familiar bell curve: it produces elliptically shaped scatter plots with no outliers or curved structure.

The code below will make a plot of bill length (bill_len) against bill depth (bill_dep) for Adélie penguins.

Code 
adelie <- penguins |>
  filter(species == "Adelie",
         !is.na(bill_len), !is.na(bill_dep))

ggplot(adelie, aes(x = bill_len, y = bill_dep)) +
  geom_point(colour = "#FF8C00", alpha = 0.6, size = 2) +
  labs(
    title = "Adélie penguins: bill length vs bill depth",
    x     = "Bill length (mm)",
    y     = "Bill depth (mm)"
  )

Choosing the null-generating mechanism

To test whether this scatter plot is consistent with a bivariate normal distribution, we need to generate null datasets that are bivariate normal with the same means, variances, and correlation as the observed data.

The null_dist() function in nullabor can generate data from a named distribution, but it works on one variable at a time. For a bivariate normal we use null_lm() — not because we are fitting a regression, but because regressing bill_dep on bill_len and rotating the residuals generates null scatter plots that preserve the linear relationship’s strength while removing any non-normal structure.

Alternatively, and more directly, we can write our own null-generating function that simulates from a bivariate normal with parameters estimated from the data.

Code 
library(MASS)  # for mvrnorm()

adelie <- penguins |>
  filter(species == "Adelie",
         !is.na(bill_len), !is.na(bill_dep))

# Estimate bivariate normal parameters from the data
bvn_params <- list(
  mu    = c(mean(adelie$bill_len), mean(adelie$bill_dep)),
  Sigma = cov(cbind(adelie$bill_len, adelie$bill_dep))
)

# Null-generating function: simulate one bivariate normal dataset
null_bvn <- function(data) {
  n    <- nrow(data)
  sims <- mvrnorm(n, mu = bvn_params$mu, Sigma = bvn_params$Sigma)
  data.frame(bill_len = sims[, 1],
             bill_dep = sims[, 2])
}

a. Why is it important that the simulated null data have the same means, variances, and correlation as the observed data? What would go wrong if we simulated from a standard bivariate normal with mean 0 and variance 1?

b. The lineup below uses null_bvn as the null-generating function. Before looking at the panels, predict: if the real data are consistent with a bivariate normal, how hard should it be to find the real panel?

c. Write down which panel you think is the real data and what feature (if any) made it stand out. If you genuinely cannot tell, write that down too — that is also a valid and informative result.

d. Decrypt.

e. Suppose 6 out of 18 students found the real panel. Compute the visual p-value. What does the result tell you about whether the bivariate normal is an adequate model for Adélie bill measurements?

f. The null here is null_bvn — simulation from a bivariate normal. Contrast this with Exercise 2, where the null was null_permute("species"). In your own words, what is the key difference in what these two null mechanisms destroy or preserve?

Exercise 4 — Is there a linear relationship between body mass and flipper length for Gentoo penguins?

Background

A scatter plot of two continuous variables with a smooth line through it is implicitly asking: is there a linear trend, or is the relationship flat? The null hypothesis is that there is no linear relationship — any apparent slope is just random variation.

The standard null-generating mechanism for this question is null_lm(). It fits a linear model to the data, then generates null datasets by rotating the residuals — a technique that preserves the marginal distributions of both variables and the overall spread, while destroying any systematic linear relationship between them.

Code 
gentoo <- penguins |>
  filter(species == "Gentoo",
         !is.na(body_mass), !is.na(flipper_len))

ggplot(gentoo,
       aes(x = flipper_len, y = body_mass)) +
  geom_point(colour = "#009ACD", alpha = 0.6, size = 2) +
  geom_smooth(method = "lm", se = TRUE,
              colour = "#009ACD", fill = "#009ACD", alpha = 0.2) +
  labs(
    title = "Gentoo penguins: flipper length vs body mass",
    x     = "Flipper length (mm)",
    y     = "Body mass (g)"
  )

Choosing the null

For this question the null-generating mechanism is:

null_lm(body_mass ~ flipper_len, method = "rotate")

method = "rotate" rotates the residuals from the fitted model, which creates new response values that are close to the fitted values but with no systematic linear trend remaining.

a. Before generating the lineup: look at the scatter plot above. Does the relationship between flipper length and body mass look roughly linear? Are there any obvious outliers or curves?

b. Why is null_permute("flipper_len") not a good null for this question, even though it would also destroy the linear relationship?

Hint: think about what permuting flipper_len does to its marginal distribution and whether values of 180 mm and 230 mm are equally plausible for the same penguin.

Now generate the lineup.

c. Look at the 20 panels. Write down which panel you think is real and the visual feature that led you there.

d. Decrypt.

e. Suppose 15 out of 17 students found the real panel. Compute the visual p-value and state your conclusion.

f. The lineup includes a geom_smooth(method = "lm") line in every panel, including the null panels. Explain why this is useful rather than misleading. What does the slope of the line in the null panels tell you?

g. Review all three lineups in this worksheet and complete the table below:

Exercise Question Null mechanism What the null destroys
2 Do species differ in flipper length? null_permute("species")
3 Are bill dims jointly normal for Adélie? Custom null_bvn()
4 Is there a linear relationship (Gentoo)? null_lm(..., "rotate")

Summary

Protocol What it is When to use it
Rorschach All panels are null — no real data Before a lineup, to calibrate your eye to what chance looks like
Lineup One real panel hidden among null panels To test whether a pattern in your data is surprising
Null mechanism Use when… What it does
null_permute("x") Testing whether group differences or correlations are real Shuffles one column, breaking its relationship with other variables
null_dist("x", dist, params) Testing whether one variable follows a specific distribution Simulates new values from the named distribution
Custom function Testing a bivariate or multivariate distributional assumption You define how to simulate from the null model
null_lm(y ~ x, "rotate") Testing whether a linear trend is real Rotates residuals to remove the linear signal while preserving marginal distributions

The three steps of a lineup test:

  1. Decide what you want to test — state the null hypothesis in plain language.
  2. Choose a null-generating mechanism that matches the null hypothesis.
  3. Show the lineup to observers, record how many find the real panel, and compute pvisual(x, K, m).

References

R packages

Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). “Welcome to the tidyverse.” Journal of Open Source Software, 4(43), 1686. https://doi.org/10.21105/joss.01686

Buja A, Cook D, Hofmann H, Lawrence M, Lee E, Swayne DF, Wickham H (2009). “Statistical inference for exploratory data analysis and model diagnostics.” Philosophical Transactions of the Royal Society A, 367(1906), 4361–4383. https://doi.org/10.1098/rsta.2009.0120

Cook D, Wickham H, Roy Chowdhury N, Hofmann H, Thulin M (2024). nullabor: Tools for Graphical Inference. R package version 0.3.15. https://dicook.github.io/nullabor/

Venables WN, Ripley BD (2002). Modern Applied Statistics with S, 4th ed. Springer, New York. ISBN 0-387-95457-0. https://www.stats.ox.ac.uk/pub/MASS4/

Data

Kaye E, Turner H, Gorman KB, Horst AM, Presmanes Hill A (2025). “Preparing the Palmer Penguins data for the datasets package in R.” Zenodo. https://doi.org/10.5281/zenodo.14902740

Horst AM, Presmanes Hill A, Gorman KB (2022). “Palmer Archipelago penguins data in the palmerpenguins R package — an alternative to Anderson’s irises.” The R Journal, 14(1), 244–254. https://doi.org/10.32614/RJ-2022-020

Gorman KB, Williams TD, Fraser WR (2014). “Ecological sexual dimorphism and environmental variability within a community of Antarctic penguins (genus Pygoscelis).” PLOS ONE, 9(3), e90081. https://doi.org/10.1371/journal.pone.0090081

Visual inference methodology

Wickham H, Cook D, Hofmann H, Buja A (2010). “Graphical inference for infovis.” IEEE Transactions on Visualization and Computer Graphics, 16(6), 973–979. https://doi.org/10.1109/TVCG.2010.161

Majumder M, Hofmann H, Cook D (2013). “Validation of visual statistical inference, applied to linear models.” Journal of the American Statistical Association, 108(503), 942–956. https://doi.org/10.1080/01621459.2013.808157

Hofmann H, Follett L, Majumder M, Cook D (2012). “Graphical tests for power comparison of competing designs.” IEEE Transactions on Visualization and Computer Graphics, 18(12), 2441–2448. https://doi.org/10.1109/TVCG.2012.230

Licence

This worksheet is released under a Creative Commons Attribution-ShareAlike 4.0 International (CC BY-SA 4.0) licence.

CC BY-SA 4.0

You are free to share (copy and redistribute in any medium or format) and adapt (remix, transform, and build upon) this material for any purpose, even commercially, provided you:

  • give appropriate credit and a link to the licence,
  • indicate if changes were made, and
  • distribute your contributions under the same licence as the original.