nlmixr2save: fits you can read, runs you don't repeat

Two problems

Problem 1: Fits you can’t open later. Nonlinear mixed-effects models can take minutes to weeks to fit. Once one is done, you want to keep it. The usual way is saveRDS(), but that ties the result to R’s binary format and to how nlmixr2 stored the fit internally.

That second tie actually broke. nlmixr2 used the qs package to compress fits. When qs was removed from CRAN, nlmixr2est moved to qs2, which broke backward compatibility for saved fits. It later moved to base R serialization, and old qs2 fits can still be read only if qs2 is installed. A fit that depends on how a particular package version stored its internals is not really an archive.

Problem 2: Running the same model again for no reason. During development you source a script, or render a Quarto/R Markdown report, over and over. Each time, every nlmixr2() call runs again, even when nothing about the model, data, or settings changed.

nlmixr2save fixes both.

A model to work with

one.cmt <- function() {
  ini({
    tka <- 0.45; label("Log Ka")
    tcl <- log(c(0, 2.7, 100)); label("Log Cl")
    tv <- 3.45; label("log V")
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.sd <- 0.7
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    linCmt() ~ add(add.sd)
  })
}

Saving a fit as text

saveFit() breaks a fit into its parts and writes each one as text. Models, including the final estimated model, become nlmixr2/rxode2 R source. Tables (estimates, objective function values, the fit table, the data) become CSV. Everything is bundled into one .zip:

fit <- nlmixr2(one.cmt, theo_sd, est = "focei",
               control = list(print = 0))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> done

saveFit(fit)  # writes fit.zip

zip::zip_list("fit.zip")[, c("filename", "uncompressed_size")]
#> # A data frame: 19 × 2
#>    filename            uncompressed_size
#>    <chr>                           <dbl>
#>  1 fit-env.R                       20423
#>  2 fit-etaObf.csv                    971
#>  3 fit-foceiModel.R                 3246
#>  4 fit-iniDf0.csv                    544
#>  5 fit-objDf.csv                     186
#>  6 fit-origData.csv                 3239
#>  7 fit-parFixed.csv                  369
#>  8 fit-parFixedDf.csv                687
#>  9 fit-parHistData.csv             32573
#> 10 fit-phiC.R                       2232
#> 11 fit-phiH.R                       2029
#> 12 fit-ranef.csv                     761
#> 13 fit-scaleInfo.rds                 477
#> 14 fit-sessioninfo.rds               959
#> 15 fit-shrink.csv                    793
#> 16 fit-time.csv                      223
#> 17 fit-ui.R                          679
#> 18 fit.csv                         40955
#> 19 fit.R                            3617

Unzip it and you can read nearly all of the fit in a text editor, a spreadsheet, or another language (a couple of small components, such as the session info, are still .rds). No R is needed to see what was estimated. In R, loadFit() rebuilds the fit:

fit2 <- loadFit("fit")
fit2$parFixedDf
#>        Parameter  Estimate          SE       %RSE Back-transformed   CI Lower
#> tka       Log Ka 0.4720882 0.208170077 44.0955909        1.6033388  1.0661802
#> tcl       Log Cl 1.0126157 0.007282174  0.7191448        2.7527921  2.7137812
#> tv         log V 3.4595932 0.045996816  1.3295441       31.8040373 29.0622807
#> add.sd           0.6954516 0.094032384 13.5210539        0.6954516  0.5111515
#>          CI Upper BSV(CV%) Shrink(SD)%
#> tka     2.4111264 69.90700    1.343056
#> tcl     2.7923639 26.98168    4.412461
#> tv     34.8044532 13.88238   10.430442
#> add.sd  0.8797517       NA          NA

A loaded fit is a working fit, not a snapshot. Diagnostic plots, VPCs, and model piping all work, so you can use it as the starting point for the next model:

plot(fit2)
vpcPlot(fit2)
fit3 <- fit2 |> model(ka <- exp(tka)) |> nlmixr2(est = "focei")

The saved fit also records the nlmixr2est and rxode2 versions that produced it. If you load it later with different versions installed, you get a warning explaining why a rerun might not match exactly.

Caching with :=

The second tool is the one you will use daily. Change <- to :=:

fit <- nlmixr2(one.cmt, theo_sd, est = "saem")   # always runs
fit := nlmixr2(one.cmt, theo_sd, est = "saem")   # runs once, then loads

That is the only change. Here is how := works:

1. It computes a hash of the model, the data, and every estimation option. For random calls, such as simulations, it also records the random seed at the start.
2. If a cache file with that hash (and, for random calls, the same starting seed) exists, := loads it. For random calls, it also restores the ending seed state, so the rest of your script sees the same random numbers it would have seen if the call had run.
3. Otherwise the call runs and the result is saved.

The first call fits the model:

system.time(
  fitS := nlmixr2(one.cmt, theo_sd, est = "saem",
                  control = list(print = 0))
)
#> [====|====|====|====|====|====|====|====|====|====
#> ====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>    user  system elapsed 
#>  23.976   3.322  10.841

The same call again loads from disk:

system.time(
  fitS := nlmixr2(one.cmt, theo_sd, est = "saem",
                  control = list(print = 0))
)
#>    user  system elapsed 
#>   0.384   0.001   0.386

Change anything that matters and it runs again. Here, the number of burn-in iterations changes:

system.time(
  fitS := nlmixr2(one.cmt, theo_sd, est = "saem",
                  control = list(print = 0, nBurn = 250))
)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#>    user  system elapsed 
#>  23.275   3.172   9.875

Hashing the data is smarter than hashing the whole data frame. The hash only uses the columns that affect estimation: standard columns, model covariates, and table$keep columns. If you add a notes column to the data, fitS still loads from the cache. If you change DV, it refits.

set.seed(42)
sim := rxSolve(fit, nStud = 10)

Options for reports and projects

Version 0.2.0 adds options for cached work that lives in a project or a rendered document:

nlmixr2save.check = FALSE is useful when the cache is committed to version control. The document renders the same way on every machine, even with different nlmixr2 versions installed. When you do want to rerun everything, nlmixr2saveInvalidate() deletes every file in nlmixr2save.dir whose name starts with nlmixr2save.prefix. Give the cache its own directory or prefix first; with the defaults, that means every file in the working directory.

Sharing a fit without the data

Often you can share a model but not the subject-level data behind it. nlmixr2saveShare() takes an existing saved fit and writes a smaller copy next to it. The original is not changed:

nlmixr2saveShare("fit")                # fit-noData.zip
nlmixr2saveShare("fit", noFit = TRUE)  # fit-noData-noFit.zip

The first copy removes the data and keeps the predictions and residuals. The second copy also removes the per-observation table, leaving the model, estimates, objective function, eta table, and parameter history. Anything that needs the data won’t work on the copy, including VPCs and refitting. The sharing vignette lists exactly what is lost.

Summary

• saveFit() / loadFit(): fits stored as R source and CSV that stay readable and don’t depend on a binary format.
• fit := nlmixr2(...): rerun only when something changed, with reproducible random seeds.

If you source scripts over and over or render reports with fits in them, try switching <- to :=.

More information:

• nlmixr2save documentation
• Using nlmixr2save
• Version tracking
• GitHub for issues and ideas

References

1. Fidler M, Wilkins JJ, Hooijmaijers R, et al. Nonlinear mixed-effects model development and simulation using nlmixr and related R open-source packages. CPT Pharmacometrics Syst Pharmacol. 2019;8(9):621-633.
2. Wang W, Hallow KM, James DA. A tutorial on RxODE: simulating differential equation pharmacometric models in R. CPT Pharmacometrics Syst Pharmacol. 2016;5(1):3-10.