Setting seeds and pinning package versions in R

Calling set.seed() makes the random steps of an R script repeatable, but only under the same random number generator settings, the same R version and the same package versions. Change any of the three and the same seed can give different numbers. This post shows that on one bootstrap interval, then lists what to record next to the seed and how to pin package versions with renv.

The worked example is a deterministic Celeus engine run, and the R script it generated is below the article. Everything it prints is read from that run. A few small extra checks were run on the same R and package versions as the engine and are marked as such.

What a seed does and does not fix

R draws random numbers from a generator, and set.seed(n) puts that generator in a known starting state. Every later draw is then a fixed sequence. Run the script again from the same state and you get the same draws.

A seed is a statement about the generator only. It says nothing about which generator is in use, which function turns its output into a sample, or which version of the package wraps the whole thing. The seed is one of at least four inputs:

1. The seed, the integer you pass to set.seed().
2. The generator settings, which RNGkind() reports: the uniform generator, the normal generator and the method sample() uses.
3. The R version and the package versions, because a release can change an algorithm or a default.
4. The position of the call, because any random draw before the step you care about moves the stream along.

A reviewer who is given only the seed has one of the four. The guide to reproducible statistical analysis lists the full package a reviewer needs, and the reproducible research checklist shows how a reviewer checks it item by item; this post goes deeper on the random part and on versions.

The example: one bootstrap interval

The example uses chickwts, a data frame that ships with R. It records the weight of 71 chicks at six weeks, each fed one of six supplements, so the six groups are independent. The question is whether chick weight differs between feeds.

A Kruskal-Wallis test answers it (Kruskal and Wallis, 1952). The choice here is partly practical: the effect size reported with it, epsilon squared for ranks, gets its confidence interval from a percentile bootstrap, so it depends on random draws and shows the problem. The engine’s Shapiro-Wilk check found no evidence against normality in any of the six groups, so for these data a one-way ANOVA or Welch test is a reasonable alternative, covered in the guide to one-way ANOVA, Welch and Kruskal-Wallis. The step-by-step reading of the test output is in the Kruskal-Wallis test in R post.

The engine ran the test with seed 20261111. The test statistic is 37.34 on 5 degrees of freedom, with p below 0.001. The effect size, from the engine’s results file, is:

epsilon-squared = 0.5335 (95% CI 0.4256 to 0.6973)
(two-sided percentile-bootstrap CI at the run's confidence level)

The interval comes from 1000 bootstrap resamples. The script line that produces it starts with the seed:

set.seed(20261111)
effectsize::rank_epsilon_squared(data$weight, data$feed, ci = 0.95,
                                 alternative = "two.sided", iterations = 1000L)

The test statistic and p-value involve no random draw at all. Only the interval does. That is the usual shape of the problem: a seed protects the one step that needs it, and nothing else.

Same seed, same interval

Rerunning the analysis with seed 20261111 gives the same interval. The worked example for this post was checked that way: it was run twice and the two results files were identical, and the generated script was rerun and gave the same statistic and p-value. This is the case a seed is good at. Same data, same script, same environment, same numbers.

A different seed, a different interval

Now change only the seed. This check and the ones after it were run on the same R 4.6.2 and package versions as the engine, outside the engine itself, so you can reproduce them by rerunning the snippet above with a different seed or setting.

seed      generator setting    epsilon squared   95% interval
20261111  default              0.5335            0.4256 to 0.6973
20261111  default, run again   0.5335            0.4256 to 0.6973
20261112  default              0.5335            0.4390 to 0.6977

The point estimate is identical, because it is computed from the data and involves no random draw. The interval bounds move: the lower bound moves in the second decimal place and the upper bound barely moves between the first and third rows. Nothing is wrong with either interval. Each is a valid draw from the same procedure, and the gap between them is the Monte Carlo noise of 1000 resamples.

A bootstrap interval reported without its seed cannot be checked. A reviewer who reruns your script with a different seed, or with none, gets a different second decimal and cannot tell whether that is noise or an error. Report the seed, and the number of resamples, as the reproducibility guide recommends.

The generator settings travel with the seed

Changing the sample.kind setting is the check that surprises people. Keep the seed at 20261111 and change only the method sample() uses to turn random numbers into indices. In R this is the sample.kind argument of RNGkind().

R 3.6.0 changed the default method sample() uses, and the release notes say the previous method can be requested with RNGkind() or RNGversion() to reproduce old results (R Core Team, 2019). The RNGkind() documentation names the old method "Rounding" and the new default "Rejection". Under the old setting the same seed gives a different interval:

seed      sample.kind              epsilon squared   95% interval
20261111  "Rejection" (default)    0.5335            0.4256 to 0.6973
20261111  "Rounding" (pre-3.6.0)   0.5335            0.4220 to 0.6932

The seed is the same, the data are the same and the point estimate is the same. The interval differs because the resampling step draws indices through sample().

The effect is easiest to see without a bootstrap. In the same session, set.seed(123); sample(1:10, 3) returned 3, 10 and 2 under the current default, and 3, 8 and 4 under "Rounding". A script written before R 3.6.0 and rerun today on the same seed can therefore resample differently. The documentation of RNGkind() describes all three settings and RNGversion() for requesting an older R version’s behaviour.

Two practical rules follow:

• Write the settings down. RNGkind() prints the three of them, and the R version belongs next to it.
• If you need to reproduce an old result, request the old behaviour explicitly rather than assuming the default has not moved.

You can also pass kind, normal.kind and sample.kind to set.seed() itself, so one line states the generator as well as the seed:

set.seed(20261111, kind = "Mersenne-Twister",
         normal.kind = "Inversion", sample.kind = "Rejection")

Naming the settings in the call means a future change to R’s defaults cannot move your numbers unnoticed. The cost is a longer line.

Where you place the call matters

A seed is a position in a stream, not a label on a function. If a script draws random numbers before the bootstrap, then adding one more draw earlier in the script changes every number after it, even though the seed has not changed.

Two habits keep this under control:

• Call set.seed() immediately before the random step you want to repeat, so that step does not depend on how much randomness came earlier in the script. The engine does this for its bootstrap: the generated script re-seeds on the line that computes the interval, so reordering the blocks around it cannot change the interval.
• Do not seed once and assume parallel code is covered. Depending on the backend, workers either start from generator states the parent seed does not control, or all inherit the same state and draw identical numbers. R’s parallel package provides independent L’Ecuyer-CMRG streams for workers, set with RNGkind("L'Ecuyer-CMRG") and clusterSetRNGStream() (R Core Team, parallel documentation). Check the documentation of the parallel tool you use before relying on a single seed.

Package versions are the other half

Everything so far assumed that effectsize was the same version. The interval above came from a specific implementation of the bootstrap in effectsize (Ben-Shachar et al., 2020), and a different version could use a different default number of resamples, a different interval method or a different way of drawing the indices. The same seed does not protect against that.

The numbers in this post came from effectsize 1.0.2 on R 4.6.2. If you rerun them on other versions and get a slightly different interval, the versions are the first suspect, and the settings in the previous section are the second.

Version drift is not limited to random steps. A test statistic, a default adjustment or a rounding rule can change with a package update. This is the reason the reproducibility guide says to keep a lockfile rather than rely on memory (Sandve et al., 2013).

Pinning versions with renv

renv gives a project its own library and a lockfile named renv.lock. The workflow is three calls:

renv::init()       # create a project library and renv.lock
renv::snapshot()   # record the package versions you are using now
renv::restore()    # install the versions recorded in renv.lock

You run init() once in the project and snapshot() whenever the set of packages changes, and you commit renv.lock with the script. A reviewer runs restore() and gets the same package versions you used (Ushey and Wickham, renv documentation).

Know what the lockfile does not do:

• It records the R version but does not install R. If the reviewer has a different R version, renv::status() reports the mismatch, but matching R is still manual. This is why the R version belongs in your notes and in the generator settings above.
• It does not record your operating system or system libraries. A package that links to a system library can still behave differently on another machine.
• It records what you used, not whether the analysis is right. A pinned wrong model is still the wrong model.

When a lockfile is not practical, save sessionInfo() output with the results. It cannot reinstall anything, but it tells the reviewer what to try to match.

Sort order and locale

Locale-dependent sort order is a source of silent change that has nothing to do with random numbers. Sorting character strings depends on the locale’s collation rules. In the same R session, sort(c("b", "A", "a", "B"), method = "shell") returned “a”, “A”, “b”, “B” under the en_US.UTF-8 locale, and “A”, “B”, “a”, “b” under the C locale. method = "radix" sorted in the C order regardless of locale (R Core Team, sort documentation).

Ordering matters whenever it decides something: the reference level of a factor created with factor() from text, the order of rows before a split, or the order in which a script assigns labels. The same script on a different machine, or in a container with a different locale, can order the groups differently, and then a model’s reference group changes, and with it the meaning and sign of the group coefficients.

Sort order is fixed the same way as the generator, by stating it. Set factor levels explicitly with factor(x, levels = c(...)), and note the locale in the session record.

Which steps need a seed at all

Not every number in an analysis is random, and seeding everything adds noise to the record. A t-test, a Kruskal-Wallis test, a correlation or a linear model fitted by least squares is a calculation on the data: the same data and the same software give the same statistic and p-value, with no seed involved. In this example the statistic, degrees of freedom and p-value were identical across every row of the tables above.

A seed matters for anything that draws random numbers inside the analysis:

• bootstrap intervals, as here;
• permutation tests and Monte Carlo p-values;
• multiple imputation and any random-start algorithm such as k-means;
• random train and test splits, and cross-validation folds;
• simulated data in a power or sensitivity analysis.

Read each function’s documentation to find out whether it draws. When it does, seed immediately before it, and record the number of draws next to the seed, since 1000 resamples and 5000 resamples are different procedures.

Check it in a fresh session

The only convincing test of a seed is to run the script from the top in a new R process, with nothing left over in the workspace. Objects created interactively, a seed set in the console an hour ago and a package attached by hand can all make a script look reproducible when it is not.

From a terminal, Rscript analysis.R starts a clean process and runs the file. Run it twice, compare the printed numbers, and compare them again after restarting your computer or moving the project folder. If the second run differs, something outside the script was supplying state.

Pair this with renv::status(), which reports whether the project library and the lockfile agree. A mismatch means the numbers you just reproduced came from versions the lockfile does not describe, and the lockfile needs a new snapshot() before you share it.

A checklist to keep next to the seed

Put these in one file, reproducibility-notes.txt or the header of the script, and update them when you rerun:

1. The seed, and the line it is called on.
2. The output of RNGkind(), or the three arguments passed to set.seed().
3. The R version and the operating system.
4. renv.lock committed with the script, or sessionInfo() saved with the results.
5. The number of resamples, permutations or imputations.
6. Explicit factor levels wherever the order matters.

For the interval above, the note would read: seed 20261111, Mersenne-Twister, Inversion, Rejection, R 4.6.2, effectsize 1.0.2, 1000 resamples. With the same data and script, that is enough for a reviewer to reproduce the interval exactly rather than only approximately.

Reporting it

In a methods section, name the interval’s method, the resamples and the seed, and point to where the versions are recorded:

The 95% confidence interval for epsilon squared was computed by a percentile bootstrap with 1000 resamples, using seed 20261111, the default Mersenne-Twister generator with rejection sampling, in R 4.6.2. Package versions are recorded in the lockfile in the supplementary analysis package.

Run it as a Celeus analysis

A Celeus run fixes a seed for the analysis, records it in the run specification, and gives you the generated R script, the results and the engine version in one package, so the seed and the version travel together. The method reference for the Kruskal-Wallis test lists what the engine computes, and the guide to reproducible statistical analysis covers the rest of what a reviewer needs.

# Results

- Test: Kruskal-Wallis test
- n: 71
- statistic: 37.34
- p-value: < 0.001
- effect size: epsilon-squared = 0.5335 (95% CI 0.4256 to 0.6973)
  (two-sided percentile-bootstrap CI at the run's confidence level)
- seed: 20261111

## Assumption checks
- normality by group (Shapiro-Wilk): context only, passed=TRUE (no evidence against normality - normality is not assumed by this rank-based test)
- sample size: passed=TRUE (adequate n)

## Post-hoc (Dunn's test; correction: Holm)
- casein vs horsebean: adj p = 2.081e-05 | rank-biserial r 0.95 [0.8708, 0.9811]
- casein vs linseed: adj p = 0.01033 | rank-biserial r 0.7847 [0.527, 0.9103]
- horsebean vs linseed: adj p = 0.583 | rank-biserial r -0.6667 [-0.8618, -0.2998]
- casein vs meatmeal: adj p = 0.6274 | rank-biserial r 0.4242 [-0.02965, 0.7331]
- horsebean vs meatmeal: adj p = 0.009217 | rank-biserial r -0.8727 [-0.9518, -0.6849]
- linseed vs meatmeal: adj p = 0.4815 | rank-biserial r -0.5227 [-0.7867, -0.09727]
- casein vs soybean: adj p = 0.09938 | rank-biserial r 0.6429 [0.2999, 0.8387]
- horsebean vs soybean: adj p = 0.08339 | rank-biserial r -0.8286 [-0.9305, -0.6078]
- linseed vs soybean: adj p = 0.99 | rank-biserial r -0.2798 [-0.6298, 0.1647]
- meatmeal vs soybean: adj p = 0.99 | rank-biserial r 0.2727 [-0.1831, 0.632]
- casein vs sunflower: adj p = 0.99 | rank-biserial r 0.006944 [-0.4338, 0.445]
- horsebean vs sunflower: adj p = 9.174e-06 | rank-biserial r -0.9833 [-0.9938, -0.9557]
- linseed vs sunflower: adj p = 0.00625 | rank-biserial r -0.8889 [-0.9552, -0.7376]
- meatmeal vs sunflower: adj p = 0.583 | rank-biserial r -0.4545 [-0.75, -0.007905]
- soybean vs sunflower: adj p = 0.0715 | rank-biserial r -0.7381 [-0.8853, -0.4563]

groups: casein (a), horsebean (b), linseed (bc), meatmeal (ac), soybean (abc), sunflower (a) - groups sharing a letter are not significantly different at the adjusted threshold.

## Power
- Power is not calculated for this test. The matching parametric test gives a rough guide.

# Celeus - reproduce.R  (run from this package folder; bundled data in data/)
# Test: Kruskal-Wallis test  |  engine 0.53.2  |  catalog c40-dc08
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261111)
.data_path <- "data/set-seed-pin-versions-r.csv"
if (!file.exists(.data_path)) stop("Run reproduce.R from the package folder (the bundled data is in data/).")
data <- readr::read_csv(.data_path, show_col_types = FALSE)  # needs the readr package
clean_cells <- function (df) 
{
    ctrl_collapse <- "[\r\n\t\v\f]+"
    zero_width_chars <- intToUtf8(c(8203L, 65279L), multiple = TRUE)
    uni_space_chars <- intToUtf8(c(160L, 8239L, 8199L), multiple = TRUE)
    clean_text <- function(x) {
        out <- x
        for (ch in zero_width_chars) out <- gsub(ch, "", out, 
            fixed = TRUE)
        out <- gsub(ctrl_collapse, " ", out)
        for (ch in uni_space_chars) out <- gsub(ch, " ", out, 
            fixed = TRUE)
        trimws(gsub(" +", " ", out))
    }
    changed_cols <- character(0)
    cells_cleaned <- 0L
    retyped <- character(0)
    for (j in seq_along(df)) {
        col <- df[[j]]
        if (!is.character(col)) 
            next
        cleaned <- clean_text(col)
        n_changed <- sum(!is.na(col) & !is.na(cleaned) & col != 
            cleaned, na.rm = TRUE)
        if (n_changed > 0L) {
            changed_cols <- c(changed_cols, names(df)[j])
            cells_cleaned <- cells_cleaned + n_changed
        }
        df[[j]] <- cleaned
        non_blank <- cleaned[!is.na(cleaned) & nzchar(cleaned)]
        if (length(non_blank) && all(grepl("^[+-]?(0|[1-9][0-9]*)([.][0-9]+)?$", 
            non_blank))) {
            parsed <- suppressWarnings(as.numeric(cleaned))
            was_blank <- is.na(cleaned) | !nzchar(cleaned)
            lost <- (is.na(parsed) | is.infinite(parsed)) & !was_blank
            if (!any(lost)) {
                finite <- parsed[!is.na(parsed)]
                all_whole <- all(finite == floor(finite))
                in_range <- all(abs(finite) <= .Machine$integer.max)
                df[[j]] <- if (all_whole && in_range) 
                  as.integer(parsed)
                else parsed
                retyped <- c(retyped, names(df)[j])
            }
        }
    }
    names(df) <- clean_text(names(df))
    if (length(changed_cols)) {
        attr(df, "ingest_clean_columns") <- changed_cols
        attr(df, "ingest_clean_cells_cleaned") <- cells_cleaned
    }
    if (length(retyped)) 
        attr(df, "ingest_clean_retyped") <- retyped
    df
}
data <- clean_cells(data)
# Complete-case analysis: drop rows missing any analysis value.
.acols <- c("weight", "feed")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- kruskal.test(`weight` ~ `feed`, data = data)
print(result)

# Effect size (epsilon_squared) and its confidence interval.
effect <- { set.seed(20261111); .es_d <- data; .es_d[["feed"]] <- as.factor(.es_d[["feed"]]); effectsize::rank_epsilon_squared(.es_d[["weight"]], .es_d[["feed"]], ci = 0.95, alternative = "two.sided", iterations = 1000L) }
print(effect)

# Post-hoc pairwise comparisons (Dunn's test; correction: Holm).
data[["feed"]] <- as.factor(data[["feed"]])
posthoc <- { .ph_d <- data; .ph_d[["feed"]] <- as.factor(.ph_d[["feed"]]); PMCMRplus::kwAllPairsDunnTest(`weight` ~ `feed`, data = .ph_d, p.adjust.method = "holm") }
print(posthoc)