openpharma / mmrm

Mixed Models for Repeated Measures (MMRM) in R.
https://openpharma.github.io/mmrm/
Other
122 stars 22 forks source link

mmrm 0.3.13's first fit is different than subsequent ones #472

Open luwidmer opened 4 hours ago

luwidmer commented 4 hours ago

Summary mmrm 0.3.13 gives a different result on the first call than subsequent ones. Note that with TMB 1.9.15 and mmrm 0.3.12 this does not happen (the result is the same every time).

I've been able to narrow it down to the TMB::config call in the warning about optimize.instantly (https://github.com/openpharma/mmrm/blob/main/R/utils.R#L536 - also, if I manually set optimize.instantly to 1, the warning does not trigger?), so this seems to be a side effect of calling TMB::config to get the configuration even without changing the config? Reproducible example that allows for switching mmrm versions easily below:

mmrm_version <- "0.3.13" # "0.3.12" for comparison
mmrm_libpath <- paste0(".libs-", mmrm_version)

if(!dir.exists(mmrm_libpath)){
  dir.create(mmrm_libpath)
}
.libPaths(c(normalizePath(mmrm_libpath), .libPaths()))

if (!dir.exists(paste0(mmrm_libpath, "/mmrm"))) {
  remotes::install_version("mmrm", version = mmrm_version, type="source", repos="https://stat.ethz.ch/CRAN", dependencies="never")
}

library(mmrm)
library(digest)
options(digits = 22)

res <- character()
cat(as.character(packageVersion("mmrm")), "\n")
for (i in 1:5) {
  set.seed(123)
  mmrmFit <- mmrm(
    formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
    data = fev_data
  )
  res[i] <- digest(coef(mmrmFit))
  cat("Res ", i, ": ", res[i], "\n")
}

Output for mmrm 0.3.13:

> source("~/.active-rstudio-document")
0.3.13 
Res  1 :  a005171052e35c971e920321dea782f6 
Res  2 :  903104ec4747c253d773e48e6ab96b32 
Res  3 :  903104ec4747c253d773e48e6ab96b32 
Res  4 :  903104ec4747c253d773e48e6ab96b32 
Res  5 :  903104ec4747c253d773e48e6ab96b32

Output for mmrm 0.3.12:

> source("~/.active-rstudio-document")
0.3.12 
Res  1 :  903104ec4747c253d773e48e6ab96b32 
Res  2 :  903104ec4747c253d773e48e6ab96b32 
Res  3 :  903104ec4747c253d773e48e6ab96b32 
Res  4 :  903104ec4747c253d773e48e6ab96b32 
Res  5 :  903104ec4747c253d773e48e6ab96b32 

OS / Environment Multiple (R 4.4 on ARM64 / MacOS, R 4.3.1 on RHEL 7)

@danielinteractive @kaskr any ideas?

danielinteractive commented 4 hours ago

Oh no 😮...

kaskr commented 1 hour ago

@luwidmer Surprisingly, this doesn't seem to be a problem with TMB::config? If I out-comment this line:

rlang::warn(msg, .frequency = "once", .frequency_id = "tmb_warn_optimization")

the results become reproducible. Seems like rlang::warn has side effects?

So, why is this warning triggered at all? That's because your are testing for the wrong flag. Should have been:

tmb_config <- TMB::config(DLL = "mmrm")
if ( ! tmb_config$tmbad_deterministic_hash) { ... ## NOTE: NULL on old versions of TMB }
danielinteractive commented 1 hour ago

Thanks @kaskr , wow that is really weird that rlang::warn seems to be the problem. I mean then I can also just go with a usual warning() - it will spam the console of the user but I guess that is less problematic than side effects. Thanks for the pointer to the right flag, I guess I was a bit too naive there when assuming the option has the same name as the list element.

luwidmer commented 1 hour ago

@kaskr super interesting, and sorry for blaming TMB, my bad 😢! I did re-run 0.3.13 through RDvalgrind with level 2 valgrind instrumentation and no uninitialized memory is used (also not by rlang) there as far as this instrumentation can detect. Commenting the warning and re-installing also did the trick for me, but I'm a bit at a loss of why rlang::warn would no longer make things reproducible... that's worrying at least 🤔

luwidmer commented 1 hour ago

@danielinteractive is it possible a different optimizer is used when this warning is fired? It doesn't show up for me in the above code, even though it should?! I very much doubt the issue is in rlang::warn

luwidmer commented 59 minutes ago

@danielinteractive yup I think that might be it, changing to warning() gives the following result:

> source("~/.active-rstudio-document")
0.3.13.9000 
Error in refit_multiple_optimizers(fit = fit, control = control) : 
  No optimizer led to a successful model fit. Please try to use a different covariance structure or other covariates.
danielinteractive commented 58 minutes ago

@luwidmer ah yeah true, the warning will be captured internally and then another optimizer will be chosen, causing the other result. Alright, then I just did a mistake there putting the warning in the TMB fitting function, I guess we will need to move it higher up into mmrm().