4  Functional Programming

4.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 4.20.

  1. What is the primary purpose of the purrr package? In what situations would you reach for it rather than base R’s lapply/sapply?
  2. Among map(), map_dbl(), map_chr(), and map_lgl(), how do they differ in return type, and what happens if an element of the result does not match the declared type?
  3. In map(1:5, ~ .x^2), what does the tilde (~) represent, and how would you rewrite this using the \() anonymous function syntax?

4.2 Learning objectives

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

  • Explain why functions are first-class objects in R, and what ‘passing a function as an argument’ buys you.
  • Write anonymous functions with \(x) ... and with the ~ .x shortcut, and know when each is appropriate.
  • Replace a for loop with the appropriate purrr::map_*() variant, choosing the typed version when you know the expected return type.
  • Use map2() and pmap() for iteration over two or more parallel vectors (as in simulation studies).
  • Use purrr::reduce() and purrr::accumulate() to collapse a list into a single value.
  • Apply purrr::safely() and purrr::possibly() to iterate robustly over messy input that includes failures.

4.3 Orientation

A recurring question in applied work is: how do I apply the same operation across many elements of a data structure? Write the computation once as a function, then pass it to an iteration tool that takes care of the mechanics. Once you can write functions, you can write functions that take functions as arguments, and that is where the leverage comes from.

R is fundamentally a functional language. Almost everything you do in R can be expressed as a composition of functions applied to data. The two common iteration toolkits are base R’s apply family (lapply, sapply, vapply, tapply, apply) and the tidyverse’s purrr. Both do the same core job. purrr does it with more consistent naming and type guarantees; apply is ubiquitous and has no dependencies. This chapter covers both and recommends purrr for most work.

4.4 The statistician’s contribution

map, lapply, and friends are boring mechanics. Anyone can look them up. The interesting questions are about what you iterate over, what the function returns, and what you do when things go wrong. These shape whether an analysis scales, whether failures are loud or silent, and whether the code is honest about what it computed.

What does the iteration unit represent? A bootstrap resample, a patient stratum, a simulation parameter combination, a file of raw data. The iteration unit is a statistical object, not just a loop index. Naming it well in code, resample, stratum, params, file_path, is worth the minute it takes, because the same code will be read dozens of times over a project’s life.

What type should the function return? A single number, a logical, a character, a tibble with one row per group, a fitted model object. The statistician’s judgement is to pick the type that makes downstream work easy, not the type that was convenient to produce. If every iteration returns a tibble, map() followed by list_rbind() gives a tidy result in one step. If iterations return mixed types, the downstream code has to re-sort the zoo.

What happens when an iteration fails? Raw iteration aborts on the first failure and returns nothing. For a simulation study with 10,000 parameter combinations, one failing combination should not destroy the other 9,999 results. safely() and possibly() are the right tools here, but they force a judgement: do you want to record the error and continue (so the failure is visible later), substitute a default value (so the downstream code does not need to special-case missing results), or abort (because the failure indicates a bug in the computation)? Choosing well separates a defensible analysis from one that either silently drops half the data or aborts mysteriously on a corner case.

When to resist the refactor. A straightforward for loop with a clear loop body is perfectly acceptable R code. The right reason to replace it with map_*() is that the functional version is clearer or safer, not ‘because loops are bad’. Dogmatic vectorisation that compresses a readable loop into an unreadable chain of |>-threaded anonymous functions is a cost without a benefit.

These are judgement calls. They do not have right answers that an LLM can look up. They are what makes functional code legibly statistical rather than showily clever.

4.5 Functions as objects

In R, functions are ordinary objects. You can assign a function to a variable, store it in a list, pass it as an argument, and return it from another function. This is what ‘functional programming’ means in practice.

square <- function(x) x^2
square(5)
#> [1] 25

# store functions in a list
funs <- list(sq = square, cb = function(x) x^3)
funs$sq(5)
#> [1] 25
funs$cb(2)
#> [1] 8

This flexibility underlies every iteration pattern in this chapter. map(x, f) works because f is just another object that can be named and passed.

4.6 Anonymous functions: three ways

For short functions you use only once, defining them inline avoids cluttering your namespace with single-use names. R offers three forms; they are equivalent for most purposes.

# classic anonymous function (verbose)
map_dbl(my_list, function(x) mean(x, na.rm = TRUE))

# R 4.1+ lambda shorthand (recommended in base R)
map_dbl(my_list, \(x) mean(x, na.rm = TRUE))

# purrr tilde shorthand (widely used in tidyverse code)
map_dbl(my_list, ~ mean(.x, na.rm = TRUE))

The \(x) ... form is the modern recommendation and is easier to read when the argument name matters (\(patient), \(fold), \(params)). The ~ .x form is more compact for single-line expressions and is ubiquitous in older tidyverse code. Neither is wrong; consistency within a file matters more than which one you pick.

.x is only meaningful inside ~ formulas. Inside \(x), the argument is whatever you named it. This is a small gotcha that catches people moving code between the two forms.

4.7 The base R apply family

The base R iteration tools predate purrr and remain widely used. A quick tour:

  • lapply(X, FUN) applies FUN to each element of X and always returns a list. Predictable but often inconvenient when you wanted a vector.
  • sapply(X, FUN) is lapply with an attempt to ‘simplify’ the result into a vector or matrix when possible. Convenient interactively, dangerous in scripts: its return type depends on the input (empty list returns list(), not a numeric vector). Avoid in production code.
  • vapply(X, FUN, FUN.VALUE) is the safe alternative to sapply. You declare the expected return type (numeric(1), character(1)) and vapply errors if the function returns something else. This is the recommended base R idiom for scripts.
  • tapply(X, INDEX, FUN) applies FUN to subgroups defined by the factor INDEX. Conceptually similar to group_by |> summarise in dplyr.
  • apply(X, MARGIN, FUN) works on matrices and arrays. MARGIN = 1 iterates over rows, MARGIN = 2 iterates over columns. Note that apply on a data frame coerces to a matrix first, which is usually not what you want.
my_list <- list(a = 1:10, b = 11:20, c = 21:30)

lapply(my_list, mean)              # list of three means
sapply(my_list, mean)              # named numeric vector
vapply(my_list, mean, numeric(1))  # same, type-checked

mat <- matrix(1:12, nrow = 3)
apply(mat, 2, sum)                 # column sums
apply(mat, 1, sum)                 # row sums

The biggest practical wart in the apply family is sapply. Its return type depends on the run-time shape of the result. If every element happens to be a length-1 numeric, you get a vector. If they vary in length, you get a matrix or a list, sometimes silently. In a script that runs on new data, this difference between runs can produce errors far from the sapply call.

4.8 The purrr toolkit

purrr reorganises the same operations under a consistent naming scheme. The core principle: the function name tells you the output type.

library(purrr)

# untyped map: always returns a list
map(my_list, mean)

# typed variants
map_dbl(my_list, mean)    # numeric vector
map_int(my_list, length)  # integer vector
map_chr(my_list, ~ paste("mean:", round(mean(.x), 2)))
map_lgl(my_list, ~ all(.x > 5))

The typed variants check the return type element-by-element and error immediately if something is wrong. This is the main reason to prefer them over sapply:

# errors immediately: "result must be length 1"
map_dbl(list(1:3, "a"), mean)

The same call with sapply would silently coerce or return a list, producing a downstream bug that is hard to trace.

Rule of thumb: if you know what type you expect, use the typed variant. Use map() only when you genuinely need a list back (for example, when each iteration produces a fitted model object).

Question. Why is map_dbl() preferred over sapply() in scripts and packages?

Answer.

sapply()’s return type depends on the run-time shape of the output. For a normal case it may return a numeric vector; for an edge case (empty input, one element returning a different length, mixed types) it silently returns a list, a matrix, or something else. Downstream code that expected a vector then fails or silently coerces. map_dbl() declares the expected output type up front: each element must be a length-1 numeric, or the call errors loudly at the point of failure. Loud errors at the source are cheaper to debug than silent type changes that surface later.

4.9 map2() and pmap(): parallel iteration

map() iterates over one input. Real analyses often iterate over parallel inputs: a vector of sample sizes and a matching vector of means, a tibble of simulation parameters with one row per combination.

# map2: two parallel inputs
means <- c(0, 5, 10)
sds   <- c(1, 2, 3)

samples <- map2(means, sds, \(m, s) rnorm(100, m, s))
length(samples)
#> [1] 3

map2 passes element i of each input to the function, iterating in lockstep. Typed variants exist: map2_dbl, map2_chr, map2_lgl, map2_int.

For three or more parallel inputs, use pmap():

params <- list(
  n    = c(100, 200, 300),
  mean = c(0,   5,  10),
  sd   = c(1,   2,   3)
)

samples <- pmap(params, \(n, mean, sd) rnorm(n, mean, sd))

pmap matches arguments by name when the input list is named (which makes the call site readable), or by position otherwise. The input list is conceptually a tibble: each element is a column, each row is one parameter combination.

For simulation studies, pmap on a tibble of parameters is the canonical pattern. Each row is one simulation scenario; each column is a tuning parameter.

4.10 imap(): iteration with an index or name

imap(.x, .f) passes both the value and its name (or index, if unnamed) to the function:

imap_chr(my_list, \(x, nm) paste(nm, ":", mean(x)))
#>          a          b          c
#>   "a : 5.5" "b : 15.5" "c : 25.5"

Useful when the iteration needs to know which element it is processing, for example, labelling plots or files by group name.

4.11 walk() for side effects

map() and its typed variants are for iteration that produces a value. Iteration that exists only for its side effect, printing, saving a file, drawing a plot, is cleaner with walk():

walk2(
  plots,
  paste0("plot_", names(plots), ".png"),
  \(p, path) ggsave(path, p, width = 6, height = 4)
)

walk() returns its input invisibly, so the call does not clutter the console when used interactively, and so it composes nicely inside a pipe.

4.12 reduce() and accumulate()

reduce() collapses a list into a single value by repeatedly applying a binary function:

reduce(1:5, `+`)
#> [1] 15        (equivalent to 1 + 2 + 3 + 4 + 5)

reduce(list(c(1,2), c(2,3), c(3,4)), union)
#> [1] 1 2 3 4

reduce2() takes two parallel lists. accumulate() is like reduce() but returns all the intermediate values, not just the final one, useful for running totals, cumulative products, and anything else whose history you want to keep.

For most statistical work, reduce() is used for combining results: merging a list of data frames into one, unioning a list of sets, collapsing a list of summary statistics into an aggregate.

4.13 Error handling: safely() and possibly()

Iteration that aborts on the first failure is fragile. In a simulation with 10,000 runs, one failing run should not destroy the other 9,999. purrr offers two wrappers that convert a function into one that handles its own errors.

safely(f) returns a function that always returns a list with two components: result (the output, or NULL on failure) and error (NULL on success, or the error object on failure).

safe_log <- safely(log)

x <- list(1, 2, "a", 4)
results <- map(x, safe_log)

# separate successes from failures
map_lgl(results, \(r) is.null(r$error))
#> [1]  TRUE  TRUE FALSE  TRUE

# just the values for successful cases
map(results, "result") |> compact()

possibly(f, otherwise = NA) returns a function that returns otherwise on error instead of aborting. Simpler than safely when you do not need the error object:

maybe_log <- possibly(log, otherwise = NA_real_)
map_dbl(x, maybe_log)
#> [1] 0.0000000 0.6931472        NA 1.3862944

The judgement call: safely() preserves information about why each iteration failed, at the cost of a more complex return structure. possibly() substitutes a default value, losing diagnostic information but producing a clean numeric vector. For one-off exploratory work, possibly() often suffices. For a production simulation where failures carry information, use safely() and inspect the errors.

Question. When would you use safely() rather than possibly()?

Answer.

Use safely() when the error messages matter, when you need to inspect later which iterations failed and why, or when the pattern of failures is itself informative about the data (e.g., ‘every run with n > 10,000 ran out of memory’). Use possibly() when you know the fallback behaviour is appropriate and you do not need the diagnostic information. possibly() returns a clean vector of the same type as the successful cases; safely() returns a list of success/error pairs that must be post-processed. For polished, reproducible simulation code, safely() is usually the right choice: the small extra complexity is worth the diagnostic richness.

4.14 Worked example: fitting multiple models

The ‘fit the same model to each stratum’ pattern is a daily occurrence in biostatistics. purrr makes it a one-liner.

library(dplyr)
library(purrr)

models <- mtcars |>
  split(mtcars$cyl) |>
  map(\(df) lm(mpg ~ wt + hp, data = df))

# extract R-squared for each
map_dbl(models, \(m) summary(m)$r.squared)
#>         4         6         8
#> 0.6867517 0.7591355 0.7299243

# extract coefficients as a tibble
map_dfr(models, \(m) as_tibble(t(coef(m))), .id = "cyl")

The .id = "cyl" argument of map_dfr (or, in newer purrr, map(...) |> list_rbind(names_to = "cyl")) inserts the list names as a new column, so the grouping variable travels with the results.

A related pattern iterates over bootstrap resamples, cross- validation folds, or Monte Carlo replicates. In each case, the iteration unit is a statistical object, and the function returns a summary (coefficient, R², AUC, …) per iteration. Collecting them with map_dfr or map(...) |> list_rbind() gives a tidy tibble ready for downstream plotting or summarisation.

Question. When you split a data frame by a grouping variable and fit the same model to each subset, why is the result naturally a list rather than a data frame?

Answer.

Each fitted model is a complex object with dozens of components (coefficients, residuals, fitted values, model frame, call, terms, …). It does not fit neatly into a single column of a data frame. A list preserves each model as a first-class object and lets downstream code extract whatever components matter for the analysis, R², coefficients, predicted values, diagnostic plots. Once you have extracted a summary per model, that summary (often a tibble or a scalar) does fit into a data frame, and map_dfr or list_rbind produces the tidy result. The list-of-models → tibble-of-summaries pattern is the canonical functional shape for this kind of analysis.

4.15 Collaborating with an LLM on iteration

LLMs are good at producing first drafts of purrr code; they are less good at diagnosing why a pipeline produced unexpected output. Three patterns are worth learning explicitly.

Prompt 1: rewriting a loop. Paste a for loop and ask: ‘rewrite this using purrr so the output is a tibble, one row per iteration. Preserve the iteration index as a column.’

What to watch for. The rewrite will probably use map_dfr() or map() |> list_rbind(). Verify the column types are what you expect, LLMs sometimes emit code that implicitly coerces a factor to character, or a Date to numeric, because each iteration returns a slightly different structure. Check with dplyr::glimpse() on both the original loop’s output and the rewritten version.

Verification. Run both versions on at least three inputs (normal case, edge case, one empty iteration) and compare with waldo::compare(original, rewritten). Treat any difference, no matter how trivial-looking, as a bug to investigate.

Prompt 2: choosing between map variants. Describe the shape of your input (a list of vectors, a named list of tibbles, two parallel vectors of parameters) and the shape you want the output to take, and ask: ‘which purrr function is appropriate, and why?’

What to watch for. LLMs sometimes default to map() and map_dfr() even when a typed variant would catch bugs earlier. If you know the output is numeric, ask explicitly for map_dbl(); if the iteration is over a tibble of parameters, ask for pmap(). A second common pattern: LLMs sometimes suggest map2() when a plain map() over one input would do, or vice versa.

Verification. Trace through the call mentally. map2(a, b, f) requires a and b to have the same length; pmap(list(a, b), f) is equivalent. map_dbl requires every element to be a length-1 numeric. Check these constraints against your real data.

Prompt 3: error handling. Describe a function that sometimes fails and ask: ‘wrap this with safely() so I can iterate over 10,000 inputs without losing the successful results when one fails. Return a tibble with columns for the result and the error message.’

What to watch for. The main risks are shape mismatches: safely(f) returns a list with result and error components, so the downstream tibble assembly needs to flatten that structure carefully. LLMs sometimes produce code that silently discards the error component, making failures invisible.

Verification. Run the pipeline on inputs that are known to include failures (for example, divide by zero, pass a negative number to sqrt, pass a missing file path to read_csv). Confirm that the output tibble has the expected number of rows, that the error column is populated for the known-failing inputs, and that the result column is non-NULL only for inputs you expect to succeed.

The meta-lesson: purrr code is short, which makes it tempting to accept LLM output at a glance. Short does not mean obviously-correct. Read every map_*() choice and every anonymous function with the same care you would give a named function.

4.16 Principle in use

Three habits separate effective functional code from clever functional code:

  1. Pick the iteration unit deliberately. Name it. map(resamples, fit_model) reads better than map(x, f).
  2. Declare return types. Typed map_*() variants catch bugs at the source. Resort to untyped map() only when iterations genuinely produce heterogeneous outputs.
  3. Plan for failure. For any non-trivial iteration, wrap the function with safely() or possibly(). Silent failures in iteration are the silent failures most likely to reach production.

Combine these three habits and purrr code becomes the clearest way to express a large class of statistical computations. Ignore them and it becomes another source of subtle bugs, just at a higher level of abstraction than a for loop would have been.

4.17 Exercises

  1. Without using any apply/map functions, write my_map(): it takes a list x and a function f and returns a list of the same length. Test it on a simple example. Compare your implementation to the source of purrr::map() (it is short).
  2. Use purrr::map_dfr() (or map() |> list_rbind()) to read every CSV in a directory and return a single tibble with a source_file column. Handle the case where one file is malformed using safely().
  3. Use purrr::reduce() to implement sum() for a numeric vector from scratch (without calling sum() or + on more than two numbers at a time).
  4. Write a simulation that varies n, mean, and sd across a grid of values using pmap() over a tibble of parameters. For each combination, fit a one-sample t-test and extract the p-value. Return a tibble of results.
  5. Run the simulation from exercise 4 with a bug that makes some parameter combinations fail (e.g., n < 2). Rewrite the pipeline with safely() so the successful runs are preserved and the failures are reported.

4.18 Further reading

4.19 Practice test

The following multiple-choice questions exercise the chapter’s content. Attempt each question before expanding the answer.

4.19.1 Question 1

What is the primary purpose of the purrr package in R?

    1. Data visualisation
    1. Statistical modelling
    1. Functional programming
    1. Database connectivity

C. purrr is a functional-programming toolkit that provides map_*() and related functions for iterating over vectors and lists with consistent, type-safe interfaces.

4.19.2 Question 2

Which purrr function would you use to apply a function to each element of a list and return the results as a numeric vector?

    1. map()
    1. map_dbl()
    1. map_chr()
    1. map_lgl()

B. map_dbl() returns a double (numeric) vector and errors if any element is not a length-1 numeric.

4.19.3 Question 3

In purrr, what does the tilde (~) operator represent in expressions like map(1:5, ~ .x^2)?

    1. A mathematical operator for exponentiation
    1. A shorthand for creating anonymous functions
    1. A pipe operator similar to %>%
    1. A logical negation operator

B. The tilde constructs a single-argument anonymous function where .x refers to the formal argument.

4.19.4 Question 4

Which of the following is the main risk of using sapply() in scripts and packages?

    1. It is slower than lapply().
    1. Its return type depends on the run-time shape of the result, producing silent type changes between runs on different data.
    1. It does not support anonymous functions.
    1. It is deprecated in recent versions of R.

B. This is the core problem vapply() (in base R) and the typed map_*() variants (in purrr) were designed to fix.

4.19.5 Question 5

You have a function fit_one(stratum) that fits a model and sometimes fails. You want to iterate over 1000 strata, preserve successful fits, and record which strata failed with their error messages. Which tool is appropriate?

    1. map()
    1. purrr::possibly(fit_one, otherwise = NA)
    1. purrr::safely(fit_one)
    1. purrr::walk(strata, fit_one)

C. safely() preserves both the result (for successful strata) and the error (for failing strata), which is what you need when you want to report which strata failed and why. possibly() discards the error; plain map() aborts at the first failure; walk() is for side effects.

4.20 Prerequisites answers

  1. purrr provides a consistent family of functional-programming tools, chiefly the map_*() functions, which apply a function over a vector or list and return a predictable type. You reach for it when you want the behaviour of lapply/sapply but with type safety, consistent naming, and tidyverse-compatible interfaces. For one-off interactive work, sapply is fine; for scripts and packages, map_dbl, map_chr, etc. are safer because they error at the source when the expected type is not produced.
  2. map() always returns a list, regardless of what its function returned. map_dbl(), map_chr(), and map_lgl() return atomic vectors of the named type and require each iteration to return a length-1 value of the matching type; they error immediately if an iteration returns something else. map() accepts heterogeneous output; the typed variants enforce homogeneity.
  3. The tilde is a shortcut for a single-argument anonymous function whose body uses .x as the formal argument. map(1:5, ~ .x^2) is equivalent to map(1:5, \(x) x^2) in R 4.1+. Both compute the squares of 1 through 5 and return the results as a list.