Full test reproducibility appears to be failing on internal tests

[gowerc]

https://statistical-engineering.pages.roche.com/-/systems-integration-tests/-/jobs/77276548/artifacts/public/logs/check-rbmi.html

  ── Failure ('test-reproducibility.R:75:9'): Results are Reproducible ───────────
  `drawobj` (`actual`) not equal to `drawobj2` (`expected`).
  
      actual$samples[[1]]$beta | expected$samples[[1]]$beta    
  [1] 4.844                    - 5.504                      [1]
  [2] 8.212                    - 6.047                      [2]
  [3] 0.724                    - -0.615                     [3]
  [4] 0.448                    - -0.599                     [4]
  [5] 3.090                    - 3.010                      [5]
  [6] 2.247                    - 2.715                      [6]
  [7] -1.141                   - 1.856                      [7]
  [8] -0.486                   - 2.180                      [8]
  
  actual$samples[[1]]$sigma$A vs expected$samples[[1]]$sigma$A
                                             [,1]        [,2]        [,3]
  - actual$samples[[1]]$sigma$A[1, ]    3.8996258  1.81540872  0.18346780
  + expected$samples[[1]]$sigma$A[1, ]  4.9890982  1.04694767 -0.13357694
  - actual$samples[[1]]$sigma$A[2, ]    1.8154087  1.88882220 -0.04483718
  + expected$samples[[1]]$sigma$A[2, ]  1.0469477  1.02407796 -0.03493389
  - actual$samples[[1]]$sigma$A[3, ]    0.1834678 -0.04483718  0.37916568
  + expected$samples[[1]]$sigma$A[3, ] -0.1335769 -0.03493389  0.35201823

Haven't checked but I think this is from the FULL test suite that doesn't get checked via the GitHub actions pipelines

[gowerc]

@gravesti / @danielinteractive / @luwidmer

I can't make up my mind if this is actually an issue or not, essentially there is a slight edge case with how we default the random seed which can cause unexpected reproducibility mismatches. The test here that is failing basically works as follows:

run_test <- function(method) {

    <preamble setup code>

    set.seed(101)
    drawobj1 <- draws(
        data = dat,
        data_ice = dat_ice,
        vars = vars,
        method = method,
        quiet = TRUE
    )

    set.seed(101)
    drawobj2 <- draws(
        data = dat,
        data_ice = dat_ice,
        vars = vars,
        method = method,
        quiet = TRUE
    )

    expect_equal(drawobj1, drawobj2)
}

run_test(method_bayes(n_samples = 5))

The "issue" is that due to R's lazy evaluation the call to method_bayes() doesn't get evaluated until its first usage within the first call to draws() which is after the set.seed() command, the default arguments to method_bayes() then does a random sample to set the internal Stan seed which advances the .Random.seed. Thus when we come to the second call to draw() the method argument has already been evaluated and doesn't need to be evaluated again thus the seed isn't advanced thus the two invocations of draws() end up running on different seeds producing different results.

My gut is saying this is just a quirk of R's lazy evaluation and that there is nothing for us to do here (other than re-write the test code to not step into this trap) but just wanted to double check if you guys agree or not.

[gowerc]

The following is a repex that demonstrates the issue in terms of pure base R stuff:

run_test <- function(mu) {
    set.seed(101)
    x <- rnorm(10, mu)
    
    set.seed(101)
    y <- rnorm(10, mu)

    all.equal(x, y)
}

run_test(runif(1))
[1] "Mean relative difference: 0.9508834"

[gravesti]

@gowerc That is a bit tricky! Given your base R example I agree it's not an rbmi problem but rather an improper test.

[luwidmer]

@gowerc & @gravesti in the reprex you could avoid this behavior by calling force(mu) in the beginning (or doing sanity checking on the argument in a more realistic use-case). Is idea portable to the rbmi case? Would be nice to avoid such headaches for the users as far as possible...

[gowerc]

@ludwidmer - I'm not sure there is much we can do to guard against this from within rbmi though. The only thing I can think of is if we make the Stan seed an explicit argument with no default instead of defaulting to a random number.

Realistically the only scenario that this would occur is when method is being provided as an expression to a function that makes multiple calls to draws() as shown in the above test. For example the following (more likely) scenario is completely fine:

run_rbmi <- function(method) {
   set.seed(101)
   draws(method)
}
run_rbmi()
run_rbmi()

I guess in theory you could do something like capture the random seed, force() the arguments and then restore the seed but I fear this may have other unforeseen consequences as I don't think its a design pattern that is used elsewhere. My gut feeling is still that this is just an unfortunate quirk of R's lazy evaluation rules and that we shouldn't try to code around it.

[gowerc]

For reference I merged a fix to the broken test but am leaving the issue open encase we do decide that we want additional changes to rbmi to address this potential issue

[gowerc]

This came up again today with a user testing the reproducibility of rbmi. For example these all produce different results

1

set.seed(101)
method <- method_bayes()
draws(..., method = method)

1.5

set.seed(101)
draws(..., method = method)  # E.g. reuse method from (1) I think this is equivalent to (2)

2

set.seed(101)
method <- method_bayes()
set.seed(101)
draws(..., method = method)

3

set.seed(101)
draws(..., method =  method_bayes())

4

set.seed(101)
draws(..., method =  method_bayes(seed = 101))

None of this is technically an issue, that is it is expected that these would give different results, but it is very confusing to end users. Again this is more just a quirk of how RNG and lazy evaluation work within R as method_bayes() is doing a random number sample.

At a minimium I think we should update our documentation just to try and make this a bit clearer and provide recommendations for how to guarentee reproducibility.

That being said I am wondering if we should issue a warning if users haven't provided an explicit seed. E.g.

warning - Stan seed not provided, generating a random one, please set an explicit seed if you need to guarantee reproducibility

Any thoughts @danielinteractive / @gravesti / @luwidmer ?

[danielinteractive]

I think the warning idea is a good one! Otherwise, I would not try to change things or workaround the R evaluation behavior, that won't be sustainable I think.

[luwidmer]

@danielinteractive in v1.6.0 we (@gowerc and I) explicitly set the Stan seed in control_bayes in the quickstart vignette - that should reduce the likelihood of new users having to worry about lazy evaluation and where the RNG is used exactly in the code. Do you both think this is good enough (and then close this)?

[danielinteractive]

Thanks @luwidmer for calling this out. Yes I think documentation is a good solution path for this. Maybe the quick start vignette could have an additional paragraph which explicitly highlights the control_bayes usage with the seed? Currently there is just a code comment, but in order for new users to really remember this or see this, it might not be enough.