Basic syntax

sites <- list(
  kern   = c(12.1, 14.3, 11.8, 15.6, 13.2),
  mojave = c(8.4,   7.9,  9.1,  8.6),
  coast  = c(18.0, 19.4, 17.5, 20.1, 18.8, 19.0)
)

# Pre-allocate output — never grow a vector inside a loop
site_means <- vector("double", length(sites))
names(site_means) <- names(sites)

for (i in seq_along(sites)) {
  site_means[i] <- mean(sites[[i]])
}
site_means

##   kern mojave  coast 
##   13.4    8.5   18.8

The three-part boilerplate — pre-allocate, index, assign — is required every single time. That is the price for full visibility.

Iterating over names instead of indices

# Often cleaner to iterate over names directly
for (nm in names(sites)) {
  cat(nm, "→ n =", length(sites[[nm]]),
      "  mean =", round(mean(sites[[nm]]), 2), "\n")
}

## kern → n = 5   mean = 13.4 
## mojave → n = 4   mean = 8.5 
## coast → n = 6   mean = 18.8

next and break — flow control inside loops

# Skip sites with fewer than 5 observations; stop at first site with mean > 17
for (nm in names(sites)) {
  if (length(sites[[nm]]) < 5) { cat("Skipping", nm, "(too few obs)\n"); next }
  if (mean(sites[[nm]]) > 17)  { cat("Found high-mean site:", nm, "\n");  break }
  cat(nm, "passed filter, mean =", round(mean(sites[[nm]]), 2), "\n")
}

## kern passed filter, mean = 13.4 
## Skipping mojave (too few obs)
## Found high-mean site: coast

When a for loop is genuinely the right choice

A loop is preferable to map() when each iteration depends on the result of the previous one — map() has no mechanism for carry-forward state:

# Simulate a random walk — each step depends on the last
set.seed(42)
n_steps  <- 200
position <- numeric(n_steps)
position[1] <- 0

for (i in 2:n_steps) {
  position[i] <- position[i - 1] + rnorm(1, mean = 0, sd = 1)
}

plot(position, type = "l", col = "#4a7c5e", lwd = 1.5,
     xlab = "Step", ylab = "Position",
     main = "Random walk — carry-forward state requires a loop",
     frame.plot = FALSE)
abline(h = 0, col = "#c8832a", lty = 2)

while and repeat

Use while when the stopping criterion is a condition, not a count:

# Newton–Raphson: find square root of 50 by iteration
x   <- 25          # initial guess
tol <- 1e-10
iter <- 0

while (abs(x^2 - 50) > tol) {
  x    <- x - (x^2 - 50) / (2 * x)   # Newton step
  iter <- iter + 1
}
cat(sprintf("sqrt(50) ≈ %.10f  (converged in %d iterations)\n", x, iter))

## sqrt(50) ≈ 7.0710678119  (converged in 6 iterations)

sqrt(50)           # confirm with built-in

## [1] 7.071068

apply() — rows and columns of a matrix

# apply() is for matrices; MARGIN = 1 → rows, MARGIN = 2 → columns
m <- matrix(c(3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8), nrow = 4)
m

##      [,1] [,2] [,3]
## [1,]    3    5    5
## [2,]    1    9    3
## [3,]    4    2    5
## [4,]    1    6    8

apply(m, MARGIN = 1, FUN = mean)   # row means

## [1] 4.333333 4.333333 3.666667 5.000000

apply(m, MARGIN = 2, FUN = sum)    # column sums

## [1]  9 22 21

apply() does not work well on data frames with mixed types — use purrr::map_dbl() or dplyr::summarise() instead.

lapply() — list in, list out (always)

# lapply() always returns a list of the same length
lapply(sites, mean)

## $kern
## [1] 13.4
## 
## $mojave
## [1] 8.5
## 
## $coast
## [1] 18.8

# To get a plain vector you must unlist() or use sapply()
unlist(lapply(sites, mean))

##   kern mojave  coast 
##   13.4    8.5   18.8

sapply() — convenient but treacherous

# Works fine when every result is a scalar
sapply(sites, mean)

##   kern mojave  coast 
##   13.4    8.5   18.8

# ⚠ Silent failure: returns a list when simplification is impossible
# (e.g. results have different lengths — hard to catch downstream)
sapply(sites, range)   # 2-row matrix this time — not always obvious

##      kern mojave coast
## [1,] 11.8    7.9  17.5
## [2,] 15.6    9.1  20.1

vapply() — the safe alternative to sapply()

# Declare the expected return type via FUN.VALUE
vapply(sites, mean, FUN.VALUE = numeric(1))

##   kern mojave  coast 
##   13.4    8.5   18.8

# vapply() errors if the result does not match — safety at zero extra cost
# vapply(sites, range, FUN.VALUE = numeric(1))  # ← would error: length 2 ≠ 1

mapply() — multiple inputs in parallel

# mapply() is like map2() / pmap() but function-first, data-second
site_means_vec <- c(kern = 12, mojave = 8, coast = 18)
site_sds_vec   <- c(kern =  2, mojave = 1, coast =  3)

set.seed(42)
mapply(function(mu, s) round(mean(rnorm(30, mu, s)), 2),
       mu = site_means_vec,
       s  = site_sds_vec)

##   kern mojave  coast 
##  12.14   7.88  18.57

Compare to the purrr equivalent (function-last, which reads better in pipes):

set.seed(42)
map2_dbl(site_means_vec, site_sds_vec,
         ~ round(mean(rnorm(30, .x, .y)), 2))

##   kern mojave  coast 
##  12.14   7.88  18.57

tapply() — grouped summaries from a vector

# tapply(values, groups, FUN) — quick group summaries
ndvi_vals <- c(0.45, 0.48, 0.50, 0.62, 0.65, 0.70, 0.30, 0.28, 0.32)
watershed <- c("W1","W1","W1","W2","W2","W2","W3","W3","W3")

tapply(ndvi_vals, watershed, mean)

##        W1        W2        W3 
## 0.4766667 0.6566667 0.3000000

# Modern equivalent using dplyr:
data.frame(ndvi = ndvi_vals, ws = watershed) |>
  dplyr::summarise(mean_ndvi = mean(ndvi), .by = ws)

##   ws mean_ndvi
## 1 W1 0.4766667
## 2 W2 0.6566667
## 3 W3 0.3000000

Side-by-side: the same task five ways

# Task: compute the coefficient of variation (CV) for each site
cv <- function(x) sd(x) / mean(x) * 100

# 1. for loop
cv_loop <- vector("double", length(sites)); names(cv_loop) <- names(sites)
for (nm in names(sites)) cv_loop[nm] <- cv(sites[[nm]])

# 2. lapply + unlist
cv_lapply <- unlist(lapply(sites, cv))

# 3. sapply  (convenient, but not type-safe)
cv_sapply <- sapply(sites, cv)

# 4. vapply  (type-safe base R)
cv_vapply <- vapply(sites, cv, numeric(1))

# 5. purrr::map_dbl  (type-safe, tidyverse-friendly)
cv_purrr  <- map_dbl(sites, cv)

# All identical?
stopifnot(all.equal(cv_loop, cv_lapply),
          all.equal(cv_loop, cv_sapply),
          all.equal(cv_loop, cv_vapply),
          all.equal(cv_loop, cv_purrr))
cat("All five approaches produce identical results.\n")

## All five approaches produce identical results.

Passing extra arguments

When your function needs more than just the current element, pass extra arguments after .f:

# Trim 10% from each tail when computing mean
map_dbl(sites, mean, trim = 0.1)

##   kern mojave  coast 
##   13.4    8.5   18.8

# Or equivalently with the formula shorthand
map_dbl(sites, ~ mean(.x, trim = 0.1))

##   kern mojave  coast 
##   13.4    8.5   18.8

Fitting a model to each site

A common ecological workflow: split data by group, fit a model to each group, extract coefficients.

# Simulate a rainfall–NDVI dataset across 4 watersheds
set.seed(123)
n_sites <- 4
obs_df <- purrr::map_dfr(seq_len(n_sites), function(s) {
  n   <- 25
  ppt <- runif(n, 200, 600)
  tibble::tibble(
    site  = paste0("WS", s),
    ppt   = ppt,
    ndvi  = 0.3 + s * 0.05 + 0.0008 * ppt + rnorm(n, 0, 0.04)
  )
})

# Split → fit → extract slope, all in one pipeline
slopes <- obs_df |>
  split(~site) |>
  map(~ lm(ndvi ~ ppt, data = .x)) |>
  map_dfr(~ broom::tidy(.x)[2, c("term", "estimate", "p.value")],
          .id = "site")
slopes

## # A tibble: 4 × 4
##   site  term  estimate  p.value
##   <chr> <chr>    <dbl>    <dbl>
## 1 WS1   ppt   0.000834 7.59e-12
## 2 WS2   ppt   0.000660 1.30e- 9
## 3 WS3   ppt   0.000922 1.90e-12
## 4 WS4   ppt   0.000923 5.55e-10

Reading many files at once

# Suppose you have data/site_A.csv, data/site_B.csv, ...
csv_files <- list.files("data", pattern = "\\.csv$", full.names = TRUE)

all_data <- csv_files |>
  set_names(~ tools::file_path_sans_ext(basename(.x))) |>
  map_dfr(readr::read_csv, .id = "site")

The set_names() step labels each element with its filename before mapping, so .id = "site" in map_dfr() can populate a site column automatically. No loop, no rbind() accumulation pattern.

Applying the same plot to many groups

site_plots <- obs_df |>
  split(~site) |>
  imap(function(df, nm) {
    ggplot(df, aes(ppt, ndvi)) +
      geom_point(color = "#4a7c5e", alpha = .7, size = 2) +
      geom_smooth(method = "lm", se = FALSE, color = "#c8832a", linewidth = .8) +
      labs(title = nm, x = "Precipitation (mm)", y = "NDVI") +
      theme_minimal(base_size = 11) +
      theme(plot.title = element_text(family = "serif", face = "bold",
                                      color = "#1e3a2f"))
  })

# Print just the first two
patchwork::wrap_plots(site_plots[1:2], ncol = 2)

imap() is a shortcut for map2(.x, names(.x), .f) — it hands you both the element and its name, which is useful any time you need the label inside the function.

What’s next?

• String Manipulation — stringr and regex for cleaning messy text data
• Data Wrangling — dplyr and tidyr for reshaping data frames
• Visualization — ggplot2 for publication-quality graphics