Repeat & package: loops, matrices, functions

The R Statistical Language · Session 3

A. Stamm

Department of Mathematics Jean Leray, UMR CNRS 6629, Nantes University

1 Welcome back

1.1 Where we left off

You can already explore a dataset: shape it with the pipe and dplyr/tidyr, and read it back with ggplot2’s grammar.

  • the pipe |> and the dplyr/tidyr verbs (filter, select, mutate, arrange, group_by/summarise);
  • base-R plots and ggplot2: histograms, scatters, facets, boxplots;
  • the transform → plot → interpret EDA loop.

The move today. Your EDA loop is full of moves you repeat: grouped summaries, reshapes, simulations. This session looks underneath them: we leave the tidyverse and work in base R, replacing copy-paste with loops, reaching into two dimensions with matrices, and turning the moves you repeat into reusable functions you keep.

1.2 Today’s map

1Foundations
2Explore your data
3Repeat & package

Session 3 is the “repeat & package” move of the course: stop doing by hand what R can do by itself:

  • repeat & branch: loops (for, while) and conditionals (if/else, ifelse);
  • reach 2D: matrices, the vector’s 2D cousin, plus linear algebra and apply;
  • package the logic: writing your own functions, scope-aware, vectorised, tested.

The thread. The Monte-Carlo experiment ties the session together: we loop over thousands of samples, and we package that loop into one reusable, tested function. Matrices give us the two-dimensional structures the algebra often lives in.

2 Loops and conditionals

2.1 The problem loops solve

Often you want to repeat the same step on many things:

  • your function models ten species;
  • you average over 1000 simulations;
  • you transform every column of a table.

You could copy-paste… or you let R repeat for you.

2.2 for loops

Repeat a block a fixed number of times: the body runs once per value:

for (i in 1:5) {
  # i takes each value of the sequence in turn
  print(i^2)
}
[1] 1
[1] 4
[1] 9
[1] 16
[1] 25
total <- 0
for (x in c(2, 4, 6, 8)) {
  # the body runs once per value
  total <- total + x
}
total
[1] 20

Afterwards i and x are normal variables: they still hold the loop’s last value.

2.3 Pre-allocate for speed

Growing a vector inside a loop is slow. Pre-allocate the result first:

n <- 1000L
# a numeric vector of n zeros
out <- numeric(length = n) 
# character(length = n): character vector with n ""
# logical(length = n): logical vector with n FALSE
# integer(length = n): integer vector with n 0L 
# vector(mode = "list", length = n): list with n NULL
for (i in 1:n) {
  out[i] <- rnorm(1)
}
par(bg = "white")
hist(out)

Rule of thumb. Decide the output’s length before the loop, create the container, then fill slots by index (out[i] <- ...). This is much faster than out <- c(out, ...) inside the loop, a habit that matters once your data grows.

2.4 while loops

Repeat while a condition holds: useful when you don’t know the number of steps in advance:

x <- 1
while (x < 1000) {
  x <- x * 2
}
x
[1] 1024

Guard against infinite loops

If the condition never becomes FALSE, the loop runs forever. Always add a safeguard, e.g.:

x <- 1
max_iter <- 10
counter <- 0
while (x < 1000 && counter < max_iter) {
  x <- x * 2
  counter <- counter + 1
}

2.5 Conditionals: if / else if / else

Branch on a condition: only one branch runs:

score <- 82
grade <- if (score >= 90) {
  "A"
} else if (score >= 80) {
  "B"
} else {
  "C"
}
grade
[1] "B"

Long vectors break if

if needs a single TRUE/FALSE: use ifelse() for vectors.

2.6 Element-wise choice: ifelse

For vectorised decisions, use ifelse(condition, when_true, when_false):

temps <- c(-2, 15, 28, 33)
ifelse(temps > 25, "hot", "not hot")
[1] "not hot" "not hot" "hot"     "hot"

It applies the branch element-wise to a whole vector: no loop needed. This is the idiomatic replacement for “for each value, choose one of two things”.

if vs ifelse.

  • Use if/else to choose which block of code runs once;
  • Use ifelse to map a vector to one of two values.

Mixing them up is a classic beginner bug.

2.7 Which tool when?

for: repeat a block a known number of times (loop over columns, pre-allocated slots).

while: repeat until a condition flips (unknown count of steps).

if/else: choose one block of code based on a condition.

ifelse: map a vector element-wise to one of two values.

*apply family: apply a function to each element (Session 1 lists, matrices below).

Often the cleanest loop is no loop at all: vectorization, ifelse, or lapply/vapply are simpler than an explicit for.

2.8 The *apply family: loop & apply

Say “apply this function to every element” instead of writing the loop by hand:

l <- list(a = 1:3, b = 4:6)
lapply(l, mean)          # a list back — keeps the original structure
$a
[1] 2

$b
[1] 5
l <- list(a = 1:3, b = 4:6)
sapply(l, mean)          # tries to simplify: a named vector
a b 
2 5 
l <- list(a = 1:3, b = 4:6)
vapply(l, mean, 0)       # like sapply, but you promise the output type
a b 
2 5

The family you first met on lists in Session 1 (and will meet again on matrices as apply, MARGIN = ...):

  • lapply always,
  • sapply when simple,
  • vapply when you want the guarantee.

It’s “loop, then apply” collapsed into one verb.

2.9 Your turn: loops & conditionals

  1. for loop printing squares of 1:10.
  2. Loop-sum 1:1000; confirm it equals sum(1:1000).
  3. Given temps <- c(-5, 4, 12, 25, 31), use ifelse to label "freezing"…"hot" by thresholds of your choice.
  4. Monte-Carlo: average 5000 samples of rexp(30) into a pre-allocated vector, then hist() it: bell-shaped despite the exponential source?

3 Matrices

3.1 What is a matrix?

A matrix is a two-dimensional rectangular array of one type (usually numeric): a vector folded into rows and columns. It’s the 2D cousin of the vector.

(m <- matrix(1:12, nrow = 3, ncol = 4))
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
m[2, ] # subset with [row, col], like a data frame
[1]  2  5  8 11
dim(m) # length-2 vector with dimensions
[1] 3 4
nrow(m) # number of rows
[1] 3
ncol(m) # number of columns
[1] 4

Key distinction. A data frame is a list of columns (columns may differ in type); a matrix is one vector in 2D (single type). For tabular data with mixed types, prefer a data frame/tibble; reach for a matrix when you need linear algebra.

3.2 Building & reshaping

A matrix is just a vector with a dim(ension) attribute; you can build or reshape them freely:

# column-major fill
(A <- matrix(1:6, nrow = 2, ncol = 3))
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
# rbind: bind by row
identical(A, rbind(1:3, 4:6))
[1] FALSE
# cbind: bind by column
cbind(A, c(7, 8))
     [,1] [,2] [,3] [,4]
[1,]    1    3    5    7
[2,]    2    4    6    8
# strip the dim -> back to a vector
dim(A) <- NULL 
A
[1] 1 2 3 4 5 6

Column-major by default

matrix() fills column by column unless you pass byrow = TRUE. A classic source of “why is my matrix transposed?” surprises.

3.3 Arithmetic & recycling

Element-wise arithmetic works like it does on vectors, and recycling applies:

(M <- matrix(1:6, nrow = 2, ncol = 3))
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
# scalar multiplies every entry
M * 2
     [,1] [,2] [,3]
[1,]    2    6   10
[2,]    4    8   12
# length-2 vector recycles down the COLUMNS
M + c(10, 20)
     [,1] [,2] [,3]
[1,]   11   13   15
[2,]   22   24   26

Watch the recycling axis. A length-2 vector recycles along columns (down), not across: get this wrong and your standardisation silently miscalculates each column.

3.4 Linear algebra

Matrices exist to do algebra: matrix product, transpose, inverse, determinant:

(A <- matrix(c(2, 4, 1, 3), nrow = 2))
     [,1] [,2]
[1,]    2    1
[2,]    4    3
A %*% t(A) # matrix product
     [,1] [,2]
[1,]    5   11
[2,]   11   25
t(A) # transpose
     [,1] [,2]
[1,]    2    4
[2,]    1    3
solve(A) # inverse
     [,1] [,2]
[1,]  1.5 -0.5
[2,] -2.0  1.0
(A <- matrix(c(2, 4, 1, 3), nrow = 2))
     [,1] [,2]
[1,]    2    1
[2,]    4    3
det(A) # determinant
[1] 2
diag(A) # extract diagonal
[1] 2 3
eigen(A) # eigen decomposition
eigen() decomposition
$values
[1] 4.5615528 0.4384472

$vectors
           [,1]       [,2]
[1,] -0.3636591 -0.5392856
[2,] -0.9315321  0.8421229

%*% is not *

  • A * B is element-wise;
  • A %*% B is the matrix product (rows × columns).

Absent the right dimensions, %*% errors, which is usually the clue you wanted a matrix, not a data frame.

3.5 Fold a data frame into a matrix

Any rectangular numeric block can be coerced to a matrix when algebra needs it:

nums <- iris[, 1:4] # 150 x 4, all numeric
head(nums, n = 9L)
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2
3          4.7         3.2          1.3         0.2
4          4.6         3.1          1.5         0.2
5          5.0         3.6          1.4         0.2
6          5.4         3.9          1.7         0.4
7          4.6         3.4          1.4         0.3
8          5.0         3.4          1.5         0.2
9          4.4         2.9          1.4         0.2
M <- as.matrix(nums) # data frame -> matrix
dim(M)
[1] 150   4
class(M)
[1] "matrix" "array"

as.matrix() forces everything to one type: if any column is character, the whole matrix becomes character. Check str() before you coerce, and keep a real data frame for mixed-type data.

3.6 Row/column summaries & apply

Summarise along one axis with the row*/col* helpers, or generalise with apply:

(M <- matrix(1:12, 3, 4))
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
# per row (MARGIN = 1) - see also rowSums(M)
rowMeans(M)
[1] 5.5 6.5 7.5
# per column (MARGIN = 2) - see also colSums(M)
colMeans(M)
[1]  2  5  8 11
# standard deviation of each row
apply(M, 1, sd)
[1] 3.872983 3.872983 3.872983
# range of each column
apply(M, 2, \(x) max(x) - min(x))
[1] 2 2 2 2

apply(M, MARGIN, fun): MARGIN = 1 goes down rows, MARGIN = 2 across columns. It’s lapply’s 2D sibling: a loop over one axis hidden inside a one-liner.

3.7 Standardise a matrix (full example)

Let’s standardize (z-score) every column of a matrix, the machinery behind a real analysis:

M   <- matrix(c(1, 100, 2, 200), 2, 2)
mu  <- colMeans(M)            # mean of each column
sig <- apply(M, 2, sd)        # sd of each column
Z   <- sweep(M, 2, mu,  "-")  # subtract each column's mean
Z   <- sweep(Z, 2, sig, "/")  # divide each column's sd
colMeans(Z)                   # ~ 0
[1] 0 0
apply(Z, 2, sd)               # ~ 1
[1] 1 1

sweep(M, MARGIN, vector, FUN) applies a function across one axis, here subtracting a column mean.

This is the pattern scale(M) automates; you now see what’s underneath.

3.8 Your turn: matrices

  1. Build matrix(1:20, 4, 5); report dim, rowSums, colMeans.
  2. Transpose it with t(); what happens to the dimensions?
  3. A <- matrix(c(4, 2, 1, 3), 2, 2): compute A %*% solve(A); what should you get? det(A)?
  4. Standardise the columns of iris[, 1:4] (as a matrix) so each column has mean 0, sd 1. Check with colMeans / apply(., 2, sd).

4 Writing custom functions

4.1 Why functions?

Copy-pasting code is how bugs multiply. A function lets you:

  • name a block of logic so its intent is readable;
  • reuse it with different inputs;
  • fix a bug once instead of in every copy.

Any analysis you do more than twice deserves a function.

4.2 Anatomy of a function

function_name <- function(arg1, arg2 = default) {
  result <- ...      # body: work on the arguments
  return(result)     # hands a value back
}
  • name: how you call it;
  • arguments: inputs (optionals get defaults);
  • body & return: code that runs & value handed back.

4.3 Scope: what the body sees

Variables created inside a function live only there: they don’t leak out:

g <- function(x) {
  internal <- x * 2    # exists only inside the call
  internal + 1
}
g(5)
[1] 11

Lexical scoping

A function can see variables from the outside, but relying on that is fragile: pass everything it needs as arguments. Explicit inputs in, explicit output out, is what makes functions testable and reusable.

4.4 Functions wrap loops

The Monte-Carlo from the loops section is a loop crying out to be a function:

clt_sim <- function(n = 30, R = 5000) {
  means <- numeric(R)
  for (r in 1:R) means[r] <- mean(rexp(n))
  means
}
par(bg = "white")
hist(clt_sim())

The payoff. The loop from “Your turn: loops” is now a named, reusable, parameterisable tool: clt_sim(n = 100) measures the same CLT with bigger samples, no copy-paste.

4.5 Vectorised functions & defaults

Most functions are vectorised when building upon vectorized functions; they work element-wise on vectors:

clamp <- function(x, lo = 0, hi = 1) {
  ifelse(x < lo, lo, ifelse(x > hi, hi, x))
}
clamp(c(-1, 0.5, 2))
[1] 0.0 0.5 1.0
  • ifelse in the body: the element-wise tool from earlier, and the vectorised clamp needs no loop at all;
  • Note the use of default arguments which makes the parameters optional.

4.6 When functions aren’t vectorised: Vectorize() & replicate()

Built-in functions are usually vectorised, but your own (or a third-party) function may only accept one scalar at a time. Two loop-free helpers:

# not vectorised: takes one date string
parse_it <- function(x) strsplit(x, "-")[[1]][1]
# only the first is parsed
parse_it(c("2024-05-01", "2024-06-15"))
[1] "2024"
# wraps it => element-wise
parse_all <- Vectorize(parse_it)
parse_all(c("2024-05-01", "2024-06-15"))
2024-05-01 2024-06-15 
    "2024"     "2024" 
# repeat an expression R times
replicate(4, mean(rexp(10)))
[1] 0.9542195 0.8648132 0.9167770 0.8644224

Loops without loops. Vectorize() turns a scalar-only function into a vector one; replicate(R, expr) re-runs an expression R times and collects results, the function-flavoured cousin of for/sapply. Works with any expression, not just a named function.

4.7 Functions compose into pipelines

The tidyverse pipes data through functions with |>:

penguins |>
  dplyr::filter(!is.na(body_mass_g)) |>
  dplyr::summarize(avg_mass = mean(body_mass_g))

# wrap a computation in your own function and pipe right in
avg_mass_by <- function(df, group) {
  df |>
    dplyr::filter(!is.na(body_mass_g)) |>
    dplyr::group_by({{ group }}) |>
    dplyr::summarize(avg_mass = mean(body_mass_g, na.rm = TRUE))
}

{ var } is the tidyverse “embrace” shortcut that passes a bare column name into a function cleanly. We’ll rely on it, and on writing tidy functions, throughout the Reproducible Workflows course.

4.8 Test your functions

A good habit: check a function against a case you already know the answer to:

stopifnot(all(clamp(c(-1, 0.5, 2)) == c(0, 0.5, 1)))

Why it matters. Writing one-line assertions (stopifnot) turns “trust me” into “checked automatically”. You’ll build on this with a proper testing framework in the Reproducible Workflows course.

4.9 Your turn: functions

  1. zscore(x) = (x - mean(x)) / sd(x), tested on c(1,2,3,4,5).
  2. mad_ratio(x) = mean(abs(x - mean(x))) / sd(x), tested on rnorm(1000).
  3. describe(x) → named list of mean, sd, median, min, max.
  4. (Stretch) Wrap the matrices standardisation into standardize_cols(M) using sweep.

5 Wrap up

5.1 Three moves, one toolbox

Repeat, reach, package.

5.2 Repeat, reach, package

One toolbox for repetition & structure.

Repeat & branch

for/while to loop, if/else and ifelse to branch, *apply and pre-allocation for speed.

Reach 2D with matrices

One vector in two dimensions: build with matrix, do algebra with %*%/solve, summarise rows/columns with apply and sweep.

Package the logic

Turn repeated work into named, vectorised, tested functions: scope-aware, pipe-friendly, reusable.

5.3 Your toolbox: source() your own .R

From one-off script to a toolbox you keep.

Bundle the functions you wrote today

Every function in this session (clt_sim(), zscore(), your loop and matrix helpers) is a tool. Drop them all into one file my_tools.R and source("my_tools.R") when an analysis needs them. That file is your first toolbox: a step toward a real package, and the exact artifact the Reproducible Workflows course picks up: source()d inside your Quarto reports, tracked by git, pinned by renv.

5.4 Next time

Reproducible Workflows: Quarto, reports, git, and publishing.

Between sessions: finish any outstanding labs and collect your best functions into your my_tools.R. The Reproducible Workflows course picks up that toolbox and turns it into reports, source()d in, versioned with git: bring your penguins work and your tools along.