26  Testing Data Analysis Workflows

Program testing can be used to show the presence of bugs, but never to show their absence.

Edsger W. Dijkstra, Notes on Structured Programming (1970)

NoteSources

Adapted from author’s lecture notes and supporting materials for a graduate practicum in biostatistics.

26.1 Prerequisites

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 26.17.

  1. Why would you test a one-off data analysis script that will only ever run on a single dataset?
  2. What is the difference between a unit test, an integration test, and an end-to-end test in the context of a data analysis pipeline?
  3. What does tinytest::expect_equal_to_reference() capture, and when is a reference test preferable to a direct value-comparison test?

26.2 Learning objectives

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

  • Write unit tests for analytic helper functions with tinytest.
  • Add an integration test that exercises the full data pipeline on a small synthetic dataset.
  • Use expect_equal_to_reference() to capture complex outputs (printed objects, summary tables).
  • Run the test suite locally via tinytest::test_all() and in CI via GitHub Actions.
  • Identify test smells (tests that pass for the wrong reason, tests that hide bugs, flaky tests).
  • Apply visual regression testing to ggplot output via vdiffr.

26.3 Orientation

Analysis code rarely gets tested, and the reason is a respectable one. Software is written to run on inputs its author has never seen; an analysis script is written to run on one dataset, once, and then to be thrown away. Testing something that runs once against known data can look like ceremony.

The reasoning has two holes. The first is that the script does not run once. It runs every time a data extract is refreshed, an exclusion is revised, or a reviewer asks for a sensitivity analysis. By then its author has forgotten which of its branches were ever exercised. The second is more serious: the claim ‘it runs on one known dataset’ is doing no work, because the failures that matter in analysis code do not stop it running. A filter that removes the wrong rows, a join that duplicates a patient, a unit conversion applied twice, all produce output. The question a test answers is not ‘did it run’ but ‘would I have noticed if it were wrong’, and for most analysis scripts the honest answer is no.

We bring in this chapter the habits of software engineering to bear on that question.

The framework used in this chapter is tinytest, because it is what projects scaffolded with zzcollab ship with and what the rest of this book assumes. It has no dependencies of its own, which matters when the point of the surrounding chapters is to pin an environment; each test file runs on its own; and the assertion vocabulary is small enough to hold in your head.

NoteA note on testthat

testthat is the framework you are most likely to meet elsewhere, and the one the companion volume uses in the package-development context. The discipline this chapter teaches, testing statistical content rather than plumbing, is identical under either, and most of the translation is mechanical. Three differences are worth knowing before you read testthat code:

  • There is no wrapper. testthat groups assertions inside test_that("description", { ... }). In tinytest each expect_*() call is itself a test, and the description goes in an info= argument. Files stay flatter, and a failure reports the file and line rather than a block name.
  • The vocabulary is smaller. tinytest has no expect_named(), expect_setequal(), expect_type() or expect_s3_class(). Use expect_equal(names(x), ...), expect_equal(sort(a), sort(b)), expect_equal(typeof(x), ...) and expect_inherits(). This is a feature rather than a gap: each is a plain comparison whose failure output shows the two values.
  • Snapshots work differently. expect_equal_to_reference() stores an .rds of the object rather than printed text.

Where a project already depends on testthat, stay with it. Mixing frameworks in one suite buys nothing.

26.4 The statistician’s contribution

Writing a test is five minutes of work. Writing a test that would actually fail if the code were wrong is the part that requires thought:

Test the statistical content, not just the plumbing. A test that ‘returns a tibble of the right shape’ is a structural check, not a test of correctness. A test that the regression coefficient on a known synthetic dataset matches the closed-form answer is a test of correctness.

Edge cases are where bugs live. Empty input, single observation, all-NA, perfect collinearity, boundary value (zero variance, all-zero counts, class with one observation). Test these deliberately.

Test the pipeline, not just the helpers. Unit tests catch helper-function bugs. Integration tests on synthetic data catch pipeline bugs, including ones where helpers work individually but the composition fails.

Don’t test what you don’t trust. Mocking out a function call inside a test of that function makes the test tautological. The test passes because you made it pass, not because the code is correct.

These judgments are what make tests useful rather than performative.

26.5 Why test analyses?

The Orientation argued the general case. Three specific occasions make it concrete, and each is a moment that arrives in ordinary projects rather than unusual ones.

The first is refactoring. You restructure a cleaning script to use dplyr where it previously used base::merge, and the output should be identical. Note what is peculiar about this situation: you have a precise expectation, you have no way to check it by looking, and the consequence of being wrong is a silently different cohort. A test that compares the output against a saved baseline settles the question in a second, and without one the honest position after a refactor is that you do not know whether it changed anything.

The second is upstream drift, and it is the case tests catch that nothing else does. Your script reads an extract from an institutional warehouse, and one quarter the warehouse begins reporting a duration in seconds where it previously reported days. Nothing errors. The column is still numeric, the join still succeeds, the model still fits, and every downstream number is wrong by a factor of 86,400. A test asserting that the column falls in a plausible range catches this on the first run after the change; no amount of reading your own code ever will, because your code is not what changed.

The third is documentation. A test of the form ‘given input X, this function returns Y’ is a statement about intent that cannot rot, because it is executed. Comments drift away from the code they describe; assertions cannot, since a false one fails the build.

Against these the cost is a few minutes per test, weighed against a silent miscalculation in something you have published under your name.

Tests come in layers, and the layers catch different failures. Figure 26.1 is the shape, with the important qualification that the pyramid familiar from software engineering is somewhat flatter for analysis code: the integration layer, a small synthetic dataset run through the whole pipeline, does more work here than it does in most applications.

flowchart TD
  U["<b>Unit</b><br/><i>one function, known input</i><br/>catches: helper bugs<br/>cost: seconds"]
  I["<b>Integration</b><br/><i>whole pipeline, 20-line<br/>synthetic dataset</i><br/>catches: composition bugs,<br/>silent filter and join errors<br/>cost: seconds"]
  E["<b>End to end</b><br/><i>raw data to final figures</i><br/>catches: everything left<br/>cost: minutes to hours"]
  U --> I --> E
Figure 26.1: The three test layers and what each one catches. Unit tests are fast and narrow; end-to-end tests are slow and broad; the integration layer, where a synthetic dataset is run through the whole cleaning pipeline, has the best cost-benefit ratio for analysis code and is the layer most often missing.

Question. A cleaning script filters to adults with filter(age >= 18). Six months later someone discovers the export encodes unknown age as -1, so those rows were silently dropped as minors. Which layer of the diagram would have caught it, and what would the test have asserted?

Answer.

The integration layer, and only if the synthetic dataset contained a -1.

That qualification is the whole lesson. A unit test on the filter would pass: given ages 17, 25 and 40, it correctly keeps two. The function is not broken. The data violated an assumption the function never stated, and no test of the function alone can see that.

An integration test catches it if the synthetic fixture exercises the case, which is why the guidance says the fixture should exercise every branch, including the branches you do not think exist. A fixture built from plausible values only would miss it too.

There are two assertions worth writing, and they fail at different times. One on the input, asserting that age falls in a plausible range before the filter runs, which fails the day the export changes. One on the output, asserting that the number of rows dropped is what the codebook predicts, which fails the day the meaning changes.

The first is a claim about the data; the second is a claim about your understanding of it. Analysis pipelines mostly break for the second reason, and mostly get tested only for the first.

26.6 Unit tests with tinytest

For a function in R/clean.R:

# R/clean.R
#' Compute age groups from numeric age
#' @param age numeric vector
age_group <- function(age) {
  cut(age,
      breaks = c(0, 18, 40, 65, Inf),
      right = FALSE,
      labels = c("under-18", "18-39", "40-64", "65+"))
}

Tests in inst/tinytest/test-clean.R:

library(tinytest)

expect_equal(as.character(age_group(c(20, 50, 75))),
             c("18-39", "40-64", "65+"),
             info = "typical adult ages bin correctly")
#> ----- PASSED      : <-->
#>  call| expect_equal(as.character(age_group(c(20, 50, 75))), c("18-39", 
#>  call| "40-64", "65+"), info = "typical adult ages bin correctly")
#>  info| typical adult ages bin correctly

# the interval starts are inclusive, because right = FALSE
expect_equal(as.character(age_group(c(18, 40, 65))),
             c("18-39", "40-64", "65+"),
             info = "boundary values fall in the upper bin")
#> ----- PASSED      : <-->
#>  call| expect_equal(as.character(age_group(c(18, 40, 65))), c("18-39", 
#>  call| "40-64", "65+"), info = "boundary values fall in the upper bin")
#>  info| boundary values fall in the upper bin
expect_equal(as.character(age_group(c(17, 39, 64))),
             c("under-18", "18-39", "40-64"),
             info = "values just below a break stay below")
#> ----- PASSED      : <-->
#>  call| expect_equal(as.character(age_group(c(17, 39, 64))), c("under-18", 
#>  call| "18-39", "40-64"), info = "values just below a break stay below")
#>  info| values just below a break stay below

# NA_real_, not NA: bare NA is logical, and cut() rejects it. The
# distinction is invisible until the test is actually run.
expect_true(is.na(age_group(NA_real_)),
            info = "numeric NA propagates rather than erroring")
#> ----- PASSED      : <-->
#>  call| expect_true(is.na(age_group(NA)), info = "numeric NA propagates rather than erroring")
#>  info| numeric NA propagates rather than erroring

expect_error(age_group("twenty"),
             info = "non-numeric input is rejected")
#> ----- PASSED      : <-->
#>  call| expect_error(age_group("twenty"), info = "non-numeric input is rejected")
#>  info| non-numeric input is rejected

Each expect_*() call is one test. There is no enclosing block, which is the most visible difference from testthat: the info= string does the work test_that()’s description does, and a failure reports the file and line number alongside it.

Run all tests:

tinytest::test_all()          # from the project root
# or, for a package:
tinytest::test_package("practicum")

The layout tinytest expects is inst/tinytest/ for the test files and a one-line tests/tinytest.R that runs them, which is what zzcollab scaffolds. Shipping tests under inst/ has a consequence worth knowing: they are installed with the package, so a user who has your compendium can run its test suite without your source tree. For an analysis whose credibility rests on being checkable by someone else, that is the right default.

26.7 Integration tests on synthetic data

Unit tests catch helper-function bugs. They do not catch bugs where the helpers work individually but the pipeline composition fails.

Integration tests run the full pipeline on a small synthetic dataset:

# inst/tinytest/test-pipeline.R
# In a package the test file would load run_pipeline from R/; it is
# defined here so the example executes as written.
run_pipeline <- function(raw) {
  raw |>
    janitor::clean_names() |>
    dplyr::filter(age >= 18) |>
    dplyr::mutate(age_group = age_group(age)) |>
    tidyr::pivot_longer(dplyr::starts_with("bp_"),
                        names_to = "visit", names_prefix = "bp_v",
                        values_to = "bp") |>
    dplyr::filter(!is.na(bp)) |>
    dplyr::mutate(visit = as.integer(visit)) |>
    dplyr::select(patient_id, age_group, sex, treatment, visit, bp)
}

# synthetic data that exercises every branch
raw <- tibble::tibble(
  patient_id = c(1, 2, 3, 4),
  age        = c(25, 50, 75, 17),     # one excluded
  bp_v1      = c(120, 140, 160, NA),
  bp_v2      = c(118, 138, NA,  NA),
  sex        = c("M", "F", "F", "M"),
  treatment  = c("placebo", "active", "placebo", "active")
)

clean <- run_pipeline(raw)

# structural assertions. tinytest has no expect_s3_class or
# expect_named, and does not need them: both are ordinary
# comparisons, and a failure prints the two values.
expect_inherits(clean, "tbl_df", info = "returns a tibble")
#> ----- PASSED      : <-->
#>  call| expect_inherits(clean, "tbl_df", info = "returns a tibble")
#>  info| returns a tibble
expect_equal(names(clean),
             c("patient_id", "age_group", "sex",
               "treatment", "visit", "bp"),
             info = "column names and order are as specified")
#> ----- PASSED      : <-->
#>  call| expect_equal(names(clean), c("patient_id", "age_group", "sex", 
#>  call| "treatment", "visit", "bp"), info = "column names and order are as specified")
#>  info| column names and order are as specified

# row count: 3 adults x 2 visits, minus NAs
expect_equal(nrow(clean), 5,
             info = "patient 4 excluded; patient 3 has one NA bp")
#> ----- PASSED      : <-->
#>  call| expect_equal(nrow(clean), 5, info = "patient 4 excluded; patient 3 has one NA bp")
#>  info| patient 4 excluded; patient 3 has one NA bp

# specific values
expect_equal(clean$bp[clean$patient_id == 1 & clean$visit == 1],
             120,
             info = "a known cell survives the pipeline unchanged")
#> ----- PASSED      : <-->
#>  call| expect_equal(clean$bp[clean$patient_id == 1 & clean$visit == 
#>  call| 1], 120, info = "a known cell survives the pipeline unchanged")
#>  info| a known cell survives the pipeline unchanged

The synthetic dataset is small (4 patients, 2 visits) but exercises each cleaning rule: age exclusion (patient 4), missing values (patient 3 visit 2). A test that runs in a fraction of a second covers the pipeline’s logic.

For each cleaning rule (filter, derive, pivot, recode), the synthetic dataset should have at least one row that exercises it.

Question. Your test asserts expect_equal(nrow(result), 100). The test passes. Does this mean the result is correct?

Answer.

No. The test asserts only that the result has 100 rows. It says nothing about which rows, what values, or whether the right rows were filtered. A bug that produces 100 rows of the wrong patients passes the test silently. To test correctness:

result        <- data.frame(patient_id = 1:100, age = rep(40:49, 10))
expected_ids  <- 1:100
expected_ages <- rep(40:49, 10)

expect_equal(nrow(result), 100)            # structural
#> ----- PASSED      : <-->
#>  call| expect_equal(nrow(result), 100)
expect_equal(sort(result$patient_id),      # correctness
             sort(expected_ids))
#> ----- PASSED      : <-->
#>  call| expect_equal(sort(result$patient_id), sort(expected_ids))
expect_equal(result$age, expected_ages)
#> ----- PASSED      : <-->
#>  call| expect_equal(result$age, expected_ages)

The combination tests structure and content. The structural test alone is a false positive generator. This pattern, ‘looks right at the surface, wrong underneath’, is the most common testing failure mode in data analysis.

26.8 Reference tests

For complex outputs that are hard to assert literally, compare against a stored reference:

fit <- lm(mpg ~ wt + hp, data = mtcars)
expect_equal_to_reference(coef(fit), "ref/mtcars-coefs.rds",
                          info = "coefficients unchanged")

expect_equal_to_reference(summary(fit)$r.squared,
                          "ref/mtcars-r2.rds")

On first run the object is written to the named file and the test passes trivially. On subsequent runs the value is compared against it; if it differs the test fails, and you decide whether the change is intended (delete the file and regenerate) or a regression.

One difference from testthat is worth understanding rather than glossing. expect_snapshot() captures printed output as text, so it also detects changes in formatting. expect_equal_to_reference() stores the object, so it compares values and ignores how they print. For analysis code the object comparison is usually what you want: a coefficient that changed matters, and a print method that gained a blank line does not.

Reference tests are useful when:

  • Output is complex and tedious to assert literally.
  • Output is human-readable formatting (tables, printed summaries, plot text).
  • You care more about ‘is this what was reviewed’ than ‘is this exactly the right value’.

Snapshots are not useful when:

  • Output changes for trivial reasons (random seed, system locale, an unstable sort). These produce false positives.
  • The stored object is too large to inspect when it changes, which makes ‘accept or investigate’ a coin toss.

26.9 Visual regression testing with vdiffr

For ggplot output, vdiffr saves an SVG of the rendered plot and compares on subsequent runs:

library(vdiffr)

fit <- lm(mpg ~ wt, data = mtcars)
p <- ggplot(broom::augment(fit), aes(.fitted, .resid)) +
      geom_point() +
      geom_hline(yintercept = 0)
expect_doppelganger("regression-residual-plot", p)

vdiffr handles cross-platform rendering issues (font differences, anti-aliasing) better than raw image comparison. Useful for catching unintended changes to plot code (an axis label change, a geom swap).

Note that vdiffr is built on testthat and expects to run inside it, so this is the one place in the chapter where the two frameworks do not substitute for each other. If visual regression matters to a project, that is a reason to use testthat there; tinysnapshot is the tinytest-native alternative and is younger.

26.10 Continuous integration

usethis::use_github_action("check-standard") sets up GitHub Actions to run R CMD check on every push and pull request, on multiple OS / R-version combinations. (Older material calls this use_github_action_check_standard(), a function usethis has since removed in favor of naming the workflow as an argument.) For a non-package analysis project, the equivalent is a workflow that runs tests:

# .github/workflows/test.yaml
name: tests

on: [push, pull_request]

jobs:
  test:
    runs-on: ubuntu-latest
    steps:
      - uses: actions/checkout@v4
      - uses: r-lib/actions/setup-r@v2
      - uses: r-lib/actions/setup-r-dependencies@v2
        with:
          packages: |
            local::.
            any::tinytest
      - name: Run tests
        run: Rscript -e 'tinytest::test_all()'

Push the workflow; on every commit, GitHub runs the tests and reports pass/fail. Regressions caught at push time are far cheaper to fix than regressions caught in production.

26.11 Test smells

Common patterns that look like tests but produce false confidence:

Tautological tests.

# foo ignores its argument and always returns 42. It is wrong for
# every input, and the test below still passes.
foo <- function(x) 42

expect_equal(foo(1), foo(1),      # always passes
             info = "foo works")
#> ----- PASSED      : <-->
#>  call| expect_equal(foo(1), foo(1), info = "foo works")
#>  info| foo works

The test cannot fail because it tests foo() against itself. Useless.

Tests that mock the tested function.

fit_model <- function(d) list()          # stub, for the test
expect_equal(typeof(fit_model(d)), "list")

You replaced the function under test with a stub; the test passes because the stub returns a list, not because fit_model works. tinytest has no mocking facility, which makes this smell harder to write by accident than it is under testthat’s with_mocked_bindings().

Brittle snapshots.

expect_equal_to_reference(do_thing(today()),
                          "ref/do_thing.rds")   # date in output

The stored object includes the current date, so the test fails every day after the first. Either take the date out of the returned object or compare the parts that should be stable, which usually means asserting on specific fields rather than on the whole object.

Overly permissive tolerances.

result <- 0.9                                # nowhere near 0.5

expect_equal(result, 0.5, tolerance = 0.5)   # passes for any value
#> ----- PASSED      : <-->
#>  call| expect_equal(result, 0.5, tolerance = 0.5)

A 100% tolerance is no test.

Tests that test the test framework, not the code.

expect_true(TRUE)
#> ----- PASSED      : <-->
#>  call| expect_true(TRUE)

These are placeholders that should be filled in. Watch for them in code reviews.

26.12 Worked example: testing a cleaning pipeline

# R/cleaning.R
clean_visits <- function(raw) {
  raw |>
    janitor::clean_names() |>
    dplyr::filter(age >= 18) |>
    tidyr::pivot_longer(
      cols = dplyr::starts_with("bp_"),
      names_to = "visit",
      names_prefix = "bp_v",
      values_to = "bp"
    ) |>
    dplyr::filter(!is.na(bp))
}

# inst/tinytest/test-cleaning.R
# In the package this line loads the function above; here it is
# already in scope, so the tests run as written.
# source("../../R/cleaning.R")

# the age filter
raw <- tibble::tibble(
  patient_id = 1:3,
  age        = c(25, 17, 30),
  bp_v1      = c(120, 110, 130),
  bp_v2      = c(122, 112, 128)
)
result <- clean_visits(raw)
expect_equal(sort(unique(result$patient_id)), c(1, 3),
             info = "cleaning excludes minors")
#> ----- PASSED      : <-->
#>  call| expect_equal(sort(unique(result$patient_id)), c(1, 3), info = "cleaning excludes minors")
#>  info| cleaning excludes minors

# the NA filter
raw <- tibble::tibble(
  patient_id = c(1, 1, 2),
  age        = c(25, 25, 30),
  bp_v1      = c(120, NA, 130),
  bp_v2      = c(NA,  118, 128)
)
result <- clean_visits(raw)
expect_equal(nrow(result), 4,
             info = "6 cells, 2 NA, so 4 rows survive")
#> ----- PASSED      : <-->
#>  call| expect_equal(nrow(result), 4, info = "6 cells, 2 NA, so 4 rows survive")
#>  info| 6 cells, 2 NA, so 4 rows survive
expect_true(all(!is.na(result$bp)),
            info = "no NA blood pressures remain")
#> ----- PASSED      : <-->
#>  call| expect_true(all(!is.na(result$bp)), info = "no NA blood pressures remain")
#>  info| no NA blood pressures remain

# the reshape
raw <- tibble::tibble(
  patient_id = 1,
  age        = 25,
  bp_v1      = 120,
  bp_v2      = 118
)
result <- clean_visits(raw)
expect_equal(sort(names(result)),
             sort(c("patient_id", "age", "visit", "bp")),
             info = "long format has the expected columns")
#> ----- PASSED      : <-->
#>  call| expect_equal(sort(names(result)), sort(c("patient_id", "age", 
#>  call| "visit", "bp")), info = "long format has the expected columns")
#>  info| long format has the expected columns
expect_equal(nrow(result), 2,
             info = "one patient with two visits gives two rows")
#> ----- PASSED      : <-->
#>  call| expect_equal(nrow(result), 2, info = "one patient with two visits gives two rows")
#>  info| one patient with two visits gives two rows

The tests cover: the age filter, the NA filter, the pivot. Each test is a single named scenario; each asserts specific values, not just structure. Running the test suite catches any regression to the cleaning logic.

26.13 Collaborating with an LLM on tests

LLMs draft tests well; the judgment about what to test needs human input.

Prompt 1: drafting tests. Paste the function and ask: ‘write tinytest tests covering happy path, edge cases, and error handling. Make the assertions specific to values, not just structure.’

What to watch for. Default LLM tests tend to be shape-checks. Push for value-checks. Edge cases: empty, NA, single observation, type mismatch.

Verification. Introduce a bug into the function and re-run the tests. If the bug is caught, the tests are useful. If not, add a more specific test.

Prompt 2: integration test on synthetic data. Describe the pipeline; ask the LLM to generate a small synthetic dataset that exercises every branch.

What to watch for. The LLM may generate a dataset that exercises the happy path but not the edges. Verify each branch (filter rules, missing-data paths, factor levels) is covered.

Verification. Add a deliberate bug to one branch; the test should fail. If it does not, the synthetic data does not exercise that branch.

Prompt 3: diagnosing a flaky test. Paste the test and the failure; ask: ‘what’s the source of non-determinism?’

What to watch for. Common causes: random number generation without seed, parallel processing without seed-aware RNG, system-time-dependent output, locale issues. The LLM should know these.

Verification. Apply the fix and run the test multiple times. If it stops being flaky, fixed.

26.14 Principle in use

Three habits separate a suite that catches bugs from one that reports green:

  1. Test the math, not just the shape. Structural checks alone produce false confidence.
  2. Integration tests on synthetic data. A 20-line synthetic dataset that exercises every branch catches pipeline bugs that unit tests miss.
  3. CI as a tripwire. GitHub Actions running tests on every push catches regressions at the commit, not at the publication.

26.15 Exercises

  1. age_group() earlier in this chapter is tested for typical ages, boundaries, NA, and non-numeric input. One gap remains: what the function does with an age outside the range the breaks cover. Find the behavior, decide whether it is the behavior you want, and write the test that pins whichever answer you choose.
  2. Write a summarize_species() function returning mean body mass and count per species from palmerpenguins::penguins. Give it one structural test and one correctness test, and say in a comment which is which and why the structural one is insufficient on its own.
  3. Take the tautological test from this chapter, expect_equal(foo(1), foo(1)), and replace it with an assertion that fails against foo <- function(x) 42 while passing against a correct implementation.

On your own project. No answer key; these need code, a repository, or a colleague you supply.

  1. Add unit tests for two functions from an existing analysis of yours. Aim for at least: a happy-path test, an edge-case test, and an error-path test.
  2. Write an integration test: a 20-line synthetic dataset that exercises every branch of your pipeline. Run the full pipeline in under a second.
  3. Set up GitHub Actions with a workflow that runs tinytest::test_all() on every push. Push and verify the workflow runs green.
  4. Introduce a deliberate bug into your pipeline. Verify the test suite catches it. If it does not, add a test that does.
  5. Apply vdiffr::expect_doppelganger to one of your project’s plots. Modify the plot’s theme slightly; verify the test fails and you can accept or reject the change.
  6. Review a colleague’s test suite rather than their code. For each test, decide whether it checks statistical content or only structure, and find at least one test that would pass even if the function under test were subtly wrong. Propose the value-level assertion that would catch the bug, and deliver the review in writing as a pull-request comment (Chapter 8).

1. What happens outside the range of the breaks.

as.character(age_group(c(-5, 0, 200)))
#> [1] NA         "under-18" "65+"

A negative age returns NA. The breaks begin at 0, and cut() places anything below the first break outside every interval. An age of 200 is absorbed into 65+, because the last break is Inf.

The asymmetry is the finding. An impossible age at the top of the range is silently accepted as plausible; an impossible age at the bottom becomes NA, which is the same value a genuinely missing age produces. Downstream, a data-entry error of -5 is indistinguishable from an age nobody recorded.

Whether that is acceptable is a judgment, and the test should pin whichever judgment you make. Pinning the current behavior:

expect_true(is.na(age_group(-5)),
            info = "age below the first break gives NA")
#> ----- PASSED      : <-->
#>  call| expect_true(is.na(age_group(-5)), info = "age below the first break gives NA")
#>  info| age below the first break gives NA

Preferring an error instead, the test comes first and the function changes to satisfy it:

age_group_strict <- function(age) {
  stopifnot(is.numeric(age))
  if (any(age < 0, na.rm = TRUE))
    stop("negative age")
  cut(age, breaks = c(0, 18, 40, 65, Inf), right = FALSE,
      labels = c("under-18", "18-39", "40-64", "65+"))
}

expect_error(age_group_strict(-5), pattern = "negative age",
             info = "negative age is rejected")
#> ----- PASSED      : <-->
#>  call| expect_error(age_group_strict(-5), pattern = "negative age", 
#>  call| info = "negative age is rejected")
#>  info| negative age is rejected
expect_equal(as.character(age_group_strict(30)), "18-39",
             info = "a valid age still bins correctly")
#> ----- PASSED      : <-->
#>  call| expect_equal(as.character(age_group_strict(30)), "18-39", info = "a valid age still bins correctly")
#>  info| a valid age still bins correctly

Either is defensible. What is not defensible is leaving the behavior untested, because then it is not a decision at all, only whatever cut() happened to do.

2. A structural test and a correctness test.

library(palmerpenguins)

summarize_species <- function(data) {
  data |>
    dplyr::filter(!is.na(body_mass_g)) |>
    dplyr::group_by(species) |>
    dplyr::summarize(mean_mass = mean(body_mass_g),
                     n = dplyr::n(), .groups = "drop")
}

# STRUCTURAL: shape only. Every line below passes even if
# mean_mass is computed as a median, a sum, or the wrong column.
out <- summarize_species(penguins)
expect_inherits(out, "tbl_df", info = "structural: a tibble")
#> ----- PASSED      : <-->
#>  call| expect_inherits(out, "tbl_df", info = "structural: a tibble")
#>  info| structural: a tibble
expect_equal(names(out), c("species", "mean_mass", "n"),
             info = "structural: column names")
#> ----- PASSED      : <-->
#>  call| expect_equal(names(out), c("species", "mean_mass", "n"), info = "structural: column names")
#>  info| structural: column names
expect_equal(nrow(out), 3, info = "structural: one row per species")
#> ----- PASSED      : <-->
#>  call| expect_equal(nrow(out), 3, info = "structural: one row per species")
#>  info| structural: one row per species

# CORRECTNESS: values the test computes independently.
adelie <- penguins$body_mass_g[penguins$species == "Adelie"]
expect_equal(out$mean_mass[out$species == "Adelie"],
             mean(adelie, na.rm = TRUE),
             info = "correctness: the mean is the mean")
#> ----- PASSED      : <-->
#>  call| expect_equal(out$mean_mass[out$species == "Adelie"], mean(adelie, 
#>  call| na.rm = TRUE), info = "correctness: the mean is the mean")
#>  info| correctness: the mean is the mean
expect_equal(out$n[out$species == "Adelie"],
             sum(!is.na(adelie)),
             info = "correctness: n counts non-missing")
#> ----- PASSED      : <-->
#>  call| expect_equal(out$n[out$species == "Adelie"], sum(!is.na(adelie)), 
#>  call| info = "correctness: n counts non-missing")
#>  info| correctness: n counts non-missing

The structural test is worth keeping, because a change in column names breaks callers. It is insufficient because every assertion in it would still pass if mean() were replaced by median(). Only the second test would fail, and only because it recomputes the answer a different way rather than comparing the function to itself.

3. Replacing the tautology.

foo <- function(x) 42          # wrong for every input

expect_equal(foo(1), 2, info = "foo doubles 1")
#> ----- FAILED[data]: <-->
#>  call| expect_equal(foo(1), 2, info = "foo doubles 1")
#>  diff| Expected '2', got '42'
#>  info| foo doubles 1
expect_equal(foo(3), 6, info = "foo doubles 3")
#> ----- FAILED[data]: <-->
#>  call| expect_equal(foo(3), 6, info = "foo doubles 3")
#>  diff| Expected '6', got '42'
#>  info| foo doubles 3

The fix is to assert against a value computed independently of the function, here by hand. Any assertion whose expected value comes from calling the function is circular, and no amount of adding cases repairs that.

26.16 Further reading

  • (Wickham & Bryan, 2023) testing chapters, the canonical testthat reference.
  • vdiffr documentation on CRAN, visual regression testing for ggplot2.
  • The testthat package vignettes.

26.17 Prerequisites answers

  1. Even a one-off analysis benefits from tests: they catch regressions when you refactor mid-analysis, they serve as executable documentation of expected behavior, and they detect silent failures when upstream data changes format. The cost of writing a handful of tests is small; the cost of silently mis-analyzing data is much larger.
  2. A unit test exercises a single function with known inputs and checks the output. An integration test runs a pipeline segment (multiple functions) on synthetic data and checks the pipeline’s output. An end-to-end test runs the full analysis pipeline from raw data to final figures/tables. Unit tests are fast and narrow; end-to-end tests are slow and broad; integration tests sit in between, with good cost-benefit for analytic pipelines.
  3. expect_snapshot() captures the output of an expression (printed text, numerical result, plot) to a file on first run. On subsequent runs it compares the new output against the saved snapshot. Use it when the output is complex (a printed object, a plot, a formatted table) and not easily expressed as a literal value. Review changes manually when they occur; do not auto-accept.