[BUG] Handle random seed in R sessions where `.Random.seed` has not been set

[andrewherren]

Description

The stochtree R package attempts to respect user-provided random seeds without modifying the global RNG state by caching the user's pre-existing .Random.seed and resetting the random number generator based on this seed after stochtree::bart() or stochtree::bcf() has run.

This is a problem if a user's session hasn't set .Random.seed (either by calling set.seed() or running a random number generation function like runif()).

Reproducing

1. Save the following lines to an R script called example.R with

library(stochtree)
data(iris)
y <- iris$Sepal.Length
X <- iris[, c("Sepal.Width", "Petal.Length", "Petal.Width")]
bart_model <- bart(
  X_train = X,
  y_train = y,
  general_params = list(random_seed = 1234)
)

2. Run at the command line

Rscript example.R

Expected behavior

This should not throw an error, bart() and bcf() should be robust to NULL-valued .Random.seed

[SG-TLR]

I had a related issue with reproducibility that I resolved to my satisfaction. I wanted results to be reproducible while preserving desired randomness in repeated calls to the bart function passing in the same data, even if there had been misceallaneous calls to the random number generator in between. I was able to accomplish this by specifying a value for .random seed argument in the general_params list: .random_seed = .Random.seed[2]. I suggest including this tip in the R documentation.

 data(iris)
 y <- iris$Sepal.Length
 X <- iris[, c("Sepal.Width", "Petal.Length", "Petal.Width")]

set.seed(1)
bart_model_1 <- bart(
  X_train = X,
  y_train = y,
  general_params = list(random_seed = .Random.seed[2])
)

#' external call to random number generator
temp <- rnorm(1)

bart_model_2 <- bart(
  X_train = X,
  y_train = y,
  general_params = list(random_seed = .Random.seed[2])
)

#' another external call to random number generator
temp <- rnorm(1)

set.seed(1)
bart_model_1_rep <- bart(
  X_train = X,
  y_train = y,
  general_params = list(random_seed = .Random.seed[2])
)

pred1 <- rowMeans(bart_model_1$y_hat_train)
pred2 <- rowMeans(bart_model_2$y_hat_train)
pred1_rep <- rowMeans(bart_model_1_rep$y_hat_train)

all.equal(pred1, pred2)  # expect these to be different
all.equal(pred1, pred1_rep) # expect TRUE

[andrewherren]

@SG-TLR thank you for writing this up, very useful indeed! I'll be sure to add a reproducibility section in the docs and include this tip