Numerical reproducibility issues

[tobiasmuetze]

Describe the bug
We are experiencing reproducibility issues when method_bayes() is used even when all the conditions for reproducibility outlined in the stan manual are fulfilled.
Our understanding is that the cause for this is the starting value for the MCMC are the REML parameter estimates that are obtained using the R package mmrm. We have observed that the model files from the R package mmrm are not numerically reproducible; see corresponding mmrm issue.
We have reproduced this issue on rbmi versions 1.2.3 and 1.2.6. We haven't observed this issue in v1.1.3, i.e., prior to the switch to the R package mmrm.
This is joint work with @luwidmer and @bailliem.

To Reproduce

We have used the code from the vignette to generate the issue.

library(rbmi)
library(tidyverse)

data("antidepressant_data")
dat <- antidepressant_data

# Use expand_locf to add rows corresponding to visits with missing outcomes to the dataset
dat <- expand_locf(
  dat,
  PATIENT = levels(dat$PATIENT), # expand by PATIENT and VISIT 
  VISIT = levels(dat$VISIT),
  vars = c("BASVAL", "THERAPY"), # fill with LOCF BASVAL and THERAPY
  group = c("PATIENT"),
  order = c("PATIENT", "VISIT")
)

# create data_ice and set the imputation strategy to JR for
# each patient with at least one missing observation
dat_ice <- dat %>% 
  arrange(PATIENT, VISIT) %>% 
  filter(is.na(CHANGE)) %>% 
  group_by(PATIENT) %>% 
  slice(1) %>%
  ungroup() %>% 
  select(PATIENT, VISIT) %>% 
  mutate(strategy = "JR")

vars <- set_vars(
  outcome = "CHANGE",
  visit = "VISIT",
  subjid = "PATIENT",
  group = "THERAPY",
  covariates = c("BASVAL*VISIT", "THERAPY*VISIT")
)

# Define which imputation method to use (here: Bayesian multiple imputation with 150 imputed datsets) 
method <- method_bayes(
  burn_in = 200,
  burn_between = 5,
  n_samples = 150,
  seed = 675442751
)

# Create samples for the imputation parameters by running the draws() function
seed <- 987
set.seed(seed)
# seed <- NULL
drawObj <- draws(
  data = dat,
  data_ice = dat_ice,
  vars = vars,
  method = method,
  quiet = TRUE
)

imputeObj <- impute(
  drawObj,
  references = c("DRUG" = "PLACEBO", 
                 "PLACEBO" = "PLACEBO")
)
imputeObj



anaObj <- analyse(
  imputations = imputeObj,
  fun = ancova,
  vars = set_vars(
    subjid = "PATIENT",
    outcome = "CHANGE",
    visit = "VISIT",
    group = "THERAPY",
    covariates = c("BASVAL")
  )
)

poolObj <- pool(
  anaObj, 
  conf.level = 0.95, 
  alternative = "two.sided"
)

Environment:

• Platform: x86_64-pc-linux-gnu (64-bit)
• R version 4.3.1
• rbmi version 1.2.3

[gowerc]

Just to say from our prior conversation that we are pretty confident this is due to the mmrm upstream issue. Once openpharma/mmrm#472 and openpharma/mmrm#462 are resolved this should be fixed.

[gowerc]

@tobiasmuetze - Given the up stream fix thats been applied to mmrm I was planning on closing this. Please let me know if you have any outstanding concerns / issues

[gowerc]

Please feel free to re-open if this comes up again 😄