1 Load packages

library(tidyverse)  # data wrangling
library(digest)
library(rstanarm)

2 Motivation

Reproducibility of results is a major concern in science and industry alike. However, there are numerous pitfalls which may threaten reproducibility. This post explores one possible issue - the fixation of start values of random numbers drawn by R.

Reasons that seeds in R across different systems may (possibly) vary (I don’t know if they do, but I heard so) may include:

different OSes
different version of OSes
different versions of R

3 User error

Maybe the most straight forward explanation is a user error. Don’t forget to call set.seed() right before the calculation. Be sure to use a pre defined seed value.

4 Your help needed

Please let me know whether (or not) you get the same random numbers, once the seed is fixed. In order to check the possible influence source, please note your OS, and the R version (or paste the session info).

Thanks!

5 Same random numbers

What set.seed() does is the fix the value of .Random.seed, which is in turn used to initiate the random number generator of R.

.Random.seed %>% head()
#> [1]      10403        624 -709689667  906139226 -133236717 1588108856

5.1 Without seed

Without seed, the random numbers drawn by R will be transient, ie., subject of change:

rnorm(1)
#> [1] -1.051114

replicate(n = 10, expr = rnorm(1))
#>  [1]  0.01176695 -1.21472604  1.86946816  0.72725355 -2.53246585 -0.17330070
#>  [7] -1.61645738  0.44063028  0.73718828 -1.89312377

5.2 With seed

set.seed(42)

rnorm(1)
#> [1] 1.370958

set.seed(42)
replicate(n = 10, expr = rnorm(1))
#>  [1]  1.37095845 -0.56469817  0.36312841  0.63286260  0.40426832 -0.10612452
#>  [7]  1.51152200 -0.09465904  2.01842371 -0.06271410

5.3 Using a hash

Let’s use a hash number to control the output more succinctly:

set.seed(42)
x <- replicate(n = 10, expr = rnorm(1))

digest(x)
#> [1] "8cbf5bafbba3d3c63d39c586f29b3a3d"

Let’s recompute and compare the hash values:

set.seed(42)
x2 <- replicate(n = 10, expr = rnorm(1))

digest(x2)
#> [1] "8cbf5bafbba3d3c63d39c586f29b3a3d"

identical(x, x2)
#> [1] TRUE

Identical, as expected.

6 Seeds in regression models

6.1 lm

set.seed(42)
lm1 <- lm(mpg ~ hp, data = mtcars)
lm1_coef <- coef(lm1)
lm1_coef
#> (Intercept)          hp 
#> 30.09886054 -0.06822828

Rerun to double check:

set.seed(42)
lm1 <- lm(mpg ~ hp, data = mtcars)
lm1_coef <- coef(lm1)
lm1_coef
#> (Intercept)          hp 
#> 30.09886054 -0.06822828

6.2 Stan mtcars

lm2 <- stan_glm(mpg ~ hp, data = mtcars, seed = 42, refresh = 0)
lm2_coef <- coef(lm2)
lm2_coef
#> (Intercept)          hp 
#> 30.11668130 -0.06820988

Rerun:

lm2 <- stan_glm(mpg ~ hp, data = mtcars, seed = 42, refresh = 0)
lm2_coef <- coef(lm2)
lm2_coef
#> (Intercept)          hp 
#> 30.11668130 -0.06820988

6.3 Stan penguins

Import data and run the model:

penguins <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/palmerpenguins/penguins.csv")

library(rstanarm)
lm3 <- stan_glm(body_mass_g ~ bill_length_mm, 
                seed = 42,
                refresh = 0,
                data = penguins)
coef(lm3)
#>    (Intercept) bill_length_mm 
#>      356.57465       87.45281

By the way:

Copy-paste the chunk above (which is self-containing in terms of data and R packages) and then type reprex::reprex() in the console (you need to have the package reprex installed). You’ll get a reproducible example on your clipboard.