A function another officer can run — R and the Tidyverse for Programme Data — Data Analysis

The handover test

An analysis is finished when someone else can run it on next quarter’s export without editing the code. Not “without much editing” — without editing.

That standard rules out the three habits every R script accumulates: a path written for one machine, a threshold changed by hand between runs, and a step that only works if you already ran the lines above it in the right order.

From script to function

The cleaning from the earlier lessons, as one function:

# R/clean_register.R

clean_register <- function(muac, gam_mm = 125, sam_mm = 115) {
  stopifnot(is.data.frame(muac))
  required <- c("child_id", "commune", "muac_mm", "oedema")
  missing <- setdiff(required, names(muac))
  if (length(missing) > 0) {
    stop("Input is missing columns: ", paste(missing, collapse = ", "))
  }

  muac |>
    dplyr::mutate(
      # A MUAC below 40 is centimetres in a millimetre field. Applied only
      # below a threshold no plausible measurement reaches.
      muac_mm = dplyr::if_else(muac_mm < 40, muac_mm * 10L, muac_mm),
      muac_mm = dplyr::if_else(dplyr::between(muac_mm, 80L, 220L), muac_mm, NA_integer_),
      oedema  = dplyr::recode(
        tolower(trimws(oedema)),
        "true" = TRUE, "y" = TRUE, "yes" = TRUE,
        "false" = FALSE, "n" = FALSE, "no" = FALSE,
        .default = NA
      ),
      assessed = !is.na(muac_mm) | !is.na(oedema),
      gam = assessed & (muac_mm < gam_mm | oedema),
      sam = assessed & (muac_mm < sam_mm | oedema)
    )
}

Four things about that shape are deliberate.

It takes a data frame and returns one. No file reading, no writing, no <<-. That is what makes it testable and what lets a notebook and a script share it.

The thresholds are arguments with defaults. A threshold that is only ever changed by editing a line will eventually be changed and not changed back. One that is an argument is visible in the call that produced the output.

It validates its input. stop() with a message naming the missing columns beats a NULL propagating into a mutate and surfacing as an unrelated error two functions later.

Nothing is grouped on the way out. A grouped frame escaping a function is the bug from the grouping lesson, arriving from a new direction.

`stopifnot()` for what cannot happen

stopifnot(
  "every row needs an identifier" = !any(is.na(muac$child_id)),
  "MUAC must be plausible after cleaning" =
    all(is.na(muac$muac_mm) | dplyr::between(muac$muac_mm, 80, 220))
)

Named expressions became the error message, which is the difference between a useful failure and Error: ... is not TRUE.

Use stop() for “this input is not what I was promised” and stopifnot() for “this cannot happen unless the code is wrong”. The distinction matters to whoever reads the failure at seven in the morning.

Passing a column name: tidy evaluation

This is the one piece of R that surprises people coming from Python, and it has exactly one thing to learn.

# Does not work.
rate_by <- function(df, group_col) {
  df |> dplyr::summarise(n = dplyr::n(), .by = group_col)
}

rate_by(muac, commune)
#> Error: object 'commune' not found

dplyr verbs evaluate their arguments inside the data frame. Your function’s group_col is a variable in the function, not a column, so dplyr looks for a column called group_col and fails.

The fix is one operator:

rate_by <- function(df, group_col) {
  df |> dplyr::summarise(n = dplyr::n(), .by = {{ group_col }})
}

rate_by(muac, commune)

{{ }} — “embrace” — means take what the caller wrote and evaluate it in the data. It is the answer for almost every function you will write around dplyr.

For a column name arriving as a string, which is what a command-line argument gives you:

rate_by_name <- function(df, group_col) {
  df |> dplyr::summarise(n = dplyr::n(), .by = dplyr::all_of(group_col))
}

rate_by_name(muac, "commune")

all_of() errors if the column does not exist; any_of() silently skips it. For a script whose output is reported, all_of() is the one you want.

To name an output column from an argument:

count_as <- function(df, group_col, out_name) {
  df |> dplyr::summarise("{out_name}" := dplyr::n(), .by = {{ group_col }})
}

The "{...}" := form is glue syntax inside a tidyverse verb. It is worth recognising; you will need it rarely.

The indicator function

indicator_table <- function(muac, gam_mm = 125, sam_mm = 115) {
  clean_register(muac, gam_mm, sam_mm) |>
    dplyr::summarise(
      screened  = dplyr::n(),
      assessed  = sum(assessed, na.rm = TRUE),
      gam_cases = sum(gam, na.rm = TRUE),
      sam_cases = sum(sam, na.rm = TRUE),
      .by = commune
    ) |>
    dplyr::mutate(
      assessment_rate = round(assessed / screened, 3),
      gam_rate = round(gam_cases / assessed, 3),
      sam_rate = round(sam_cases / assessed, 3)
    ) |>
    dplyr::arrange(dplyr::desc(gam_rate))
}

Note what it does not do: it does not drop communes with a low assessment rate. Dropping one would change the district denominator and nothing in the output would say so. Flagging is the caller’s job, and the caller can see the rate.

Testing the part that computes

Once the computation is a function taking a frame, a test is three lines:

# tests/testthat/test-indicator-table.R
test_that("the denominator excludes unmeasured children", {
  muac <- tibble::tibble(
    child_id = c("a", "b", "c"),
    commune  = "X",
    muac_mm  = c(110L, 130L, NA_integer_),
    oedema   = c("false", "false", NA_character_)
  )

  out <- indicator_table(muac)

  expect_equal(out$screened, 3)
  expect_equal(out$assessed, 2)     # the unmeasured child is not a denominator
  expect_equal(out$gam_rate, 0.5)
})

That test encodes the decision from the grouping lesson — which children are in the denominator — in a form that fails if someone changes it. Worth writing for every indicator whose definition you had to argue about.

usethis::use_testthat() sets the directory up; devtools::test() runs it.

The command-line script

#!/usr/bin/env Rscript
# scripts/indicator_table.R

suppressPackageStartupMessages({
  library(optparse)
  library(readr)
})

source(here::here("R", "clean_register.R"))
source(here::here("R", "indicator_table.R"))

option_list <- list(
  make_option("--input", type = "character", help = "Screening register CSV."),
  make_option("--output", type = "character", help = "Directory for the tables."),
  make_option("--gam-threshold", type = "integer", default = 125L,
              help = "MUAC cut-off in mm for GAM [default %default].")
)

opt <- parse_args(OptionParser(option_list = option_list))
if (is.null(opt$input) || is.null(opt$output)) {
  stop("--input and --output are required.")
}

muac <- read_register(opt$input)
message("Read ", nrow(muac), " rows from ", opt$input)

table <- indicator_table(muac, gam_mm = opt$`gam-threshold`)

thin <- table$commune[table$assessment_rate < 0.8]
if (length(thin) > 0) {
  warning("Assessment rate below floor: ", paste(thin, collapse = ", "))
}

dir.create(opt$output, recursive = TRUE, showWarnings = FALSE)
destination <- file.path(opt$output, "gam_by_commune.csv")
write_csv(table, destination)
message("Wrote ", destination, " (", nrow(table), " communes)")

Rscript scripts/indicator_table.R \
  --input data/raw/muac-screening-artibonite-2024.v1.csv \
  --output outputs/tables

message() rather than print(): messages go to standard error, so the script’s narration does not end up mixed into piped output. warning() rather than message() for the thin communes, because it is a condition a caller might want to escalate with options(warn = 2).

Recording what produced the output

run_metadata <- function(opt) {
  list(
    generated_at  = format(Sys.time(), "%Y-%m-%dT%H:%M:%S%z"),
    input         = opt$input,
    gam_threshold = opt$`gam-threshold`,
    r_version     = as.character(getRversion()),
    dplyr         = as.character(packageVersion("dplyr"))
  )
}

jsonlite::write_json(run_metadata(opt), file.path(opt$output, "run.json"),
                     auto_unbox = TRUE, pretty = TRUE)

Six months later, “which threshold produced this table” has an answer that does not depend on anyone’s memory.

Then the real test, from the first lesson:

rm -rf outputs/
Rscript scripts/indicator_table.R --input data/raw/muac.csv --output outputs/tables
git status

If that does not restore the outputs exactly, something is not reproducible.

Where this course ends

You can build a project that reopens a year later with its package versions recorded, read any export without corrupting it, keep the labels a survey file carries, put categories in the order a report needs, aggregate to a numerator and a denominator that came from one call, reshape a flattened repeat group, and hand the result to someone else as a script rather than as a favour.

What this course deliberately did not teach is what to compute — which indicator, which denominator, which threshold, and how to defend it. That is Data Analysis Foundations, and the sector courses under Sector Analysis on the programme roadmap take it further into the classifications your cluster holds you to.

If your team works in Python rather than R, Python for Programme Data covers the same ground in that language — and the table of reversals in lesson 6 is the fastest way to move between them.