Getting Started with bbr.bayes and NONMEM Bayes • bbr.bayes

Introduction

This vignette takes you through some basic scenarios for modeling with NONMEM Bayes using bbr.bayes, introducing you to its standard workflow and functionality.

library(bbr)
library(bbr.bayes)
library(dplyr)

bbr.bayes is an extension to bbr which focuses on traceable and reproducible modeling using Bayesian methods. Note that you will also need to load bbr to use most bbr.bayes functionality.

The bbr “Getting Started with bbr and NONMEM” vignette provides an introduction to modeling with NONMEM using bbr. This vignette covers modeling in NONMEM using Bayesian estimation methods, highlighting ways NONMEM Bayes models in bbr.bayes extend or are different from NONMEM models in bbr.

See the “Getting Started with bbr.bayes and Stan” vignette for an introduction to modeling with Stan.

Installing and configuring bbi

bbi is required to work with NONMEM Bayes models. Below is a high-level overview, but please see the “Installing bbi” section of bbr’s “Getting Started with bbr and NONMEM” vignette for detailed instructions.

First, point the bbr.bbi_exe_path option to a path where bbi is or should be installed. This should be done for each R session, so you likely want to include it in your project’s .Rprofile.

options("bbr.bbi_exe_path" = file.path(getwd(), "bin", "bbi"))

Then, if bbi is not already installed at the above location, call use_bbi() to install it.

use_bbi()

Finally, use bbi_init() to initialize a directory in which you plan to store models.

MODEL_DIR <- "model/nonmem"
bbi_init(
  .dir = MODEL_DIR,
  .nonmem_dir = "/opt/NONMEM",
  .nonmem_version = "nm75"
)

This needs to be done only once for a given directory.

Creating a NONMEM Bayes model object

As in bbr, the bbr.bayes workflow is organized around the concept of model objects. These are S3 objects which you will pass to bbr and bbr.bayes functions that submit, summarize, or manipulate models. A model object points to a path on disk; the set of files that define a model is derived from this path.

A NONMEM Bayes model is similar to the “regular” NONMEM model from bbr, but it has a more complex structure underneath that is used to represent each MCMC chain. See the NONMEM Bayes models inbbr reference page for more details.

Creating a model object de novo

A NONMEM bbr.bayes can be created by calling new_model() with the path to your control stream without the file extension, in the same way you would create a standard NONMEM model with bbr, as described in the “Getting Started with bbr and NONMEM” vignette. To tell new_model() to create a NONMEM Bayes model, pass “nmbayes” as its .model_type argument.

mod1000 <- new_model(file.path(MODEL_DIR, 1000), .model_type = "nmbayes")
mod1000$model_type
#> [1] "nmbayes"

Creating a model object from an existing model

While you can create a fresh NONMEM Bayes model with new_model(), more commonly you will create a new NONMEM Bayes model from an existing model, either a regular NONMEM model (bbi_nonmem_model) or a NONMEM Bayes model (bbi_nmbayes_model).

From a regular NONMEM model

To create a NONMEM Bayes model from a regular NONMEM model, use copy_model_as_nmbayes().

mod100 <- read_model(file.path(MODEL_DIR, 100))
mod1100 <- copy_model_as_nmbayes(mod100, 1100, .inherit_tags = TRUE) %>%
  update_model_id() %>%
  replace_tag("FOCE", "BAYES")
#> replacing 100 with 1100 in /tmp/RtmpvzWHhH/nmbayes-vignette-47f074b9e1d1/model/nonmem/bayes/1100.ctl


mod100$model_type
#> [1] "nonmem"

mod1100$model_type
#> [1] "nmbayes"

That will copy over the control stream and YAML , like copy_model_from(), and then switch the model type to “nmbayes”. In addition, it will replace the estimation methods with two estimation records that serve as a template for defining a NONMEM Bayes run.


; TODO: This model was copied by bbr.bayes::copy_model_as_nmbayes().
;       nmbayes models require a METHOD=CHAIN estimation record and a
;       METHOD=BAYES or METHOD=NUTS estimation record. The records
;       below are meant as a starting point.  At the very least, you
;       need to adjust the number of iterations (see NITER option in
;       the second $EST block), but please review all options
;       carefully.
;
;       See ?bbr.bayes::bbr_nmbayes and the NONMEM docs for details.
$EST METHOD=CHAIN FILE=1100.chn NSAMPLE=4 ISAMPLE=0 SEED=1
     CTYPE=0 IACCEPT=0.3 DF=10 DFS=0

$EST METHOD=NUTS SEED=1 NBURN=250 NITER=NNNN
     AUTO=2 CTYPE=0 OLKJDF=2 OVARF=1
     NUTS_DELTA=0.95 PRINT=10 MSFO=1100.msf RANMETHOD=P PARAFPRINT=10000
     BAYES_PHI_STORE=1

From a NONMEM Bayes model

To create a NONMEM Bayes from another one, use copy_model_from(), just as you would in bbr to copy and iterate on a model.

mod1101 <- copy_model_from(mod1100) %>%
  update_model_id()
#> replacing 1100 with 1101 in /tmp/RtmpvzWHhH/nmbayes-vignette-47f074b9e1d1/model/nonmem/bayes/1101.ctl


# Edit control stream...

mod1101 <- add_notes(mod1101, "Switched from foo to bar")

Submitting a model

To execute a model, pass it to submit_model().

submit_model(mod1100)

Underneath the submit_model() method of bbi_nmbayes_model objects creates several sub-models, one for generating the initial values and one for each MCMC chain.

dir(get_output_dir(mod1100))
#>  [1] "1100-1"       "1100-1.ctl"   "1100-1.yaml"  "1100-2"       "1100-2.ctl"  
#>  [6] "1100-2.yaml"  "init"         "init.chn"     "init.ctl"     "init.yaml"   
#> [11] "nmbayes.json"

You shouldn’t often need to directly access these sub-models because bbr.bayes implements custom methods and functions to work with them behind the scenes. For example, to monitor running models, you can use tail_lst() and tail_output() in the same way as you would with regular NONMEM models. They will print the result for each chain.

Summarizing a model

`read_fit_model()`

Once a model finishes running, an important entry point to summarizing the result is the ?posterior::draws_array object returned by read_fit_model().

fit <- read_fit_model(mod1100)
dim(fit)
#> [1]  50   2 985


# Note that real models should use more chains and iterations!
posterior::niterations(fit)
#> [1] 50

posterior::nchains(fit)
#> [1] 2

posterior::nvariables(fit)
#> [1] 985

posterior provides various functions for inspecting and summarizing posterior draws. For example, you can use posterior::summarize_draws() to generate a table of summary statistics and diagnostics:

posterior::subset_draws(fit, "THETA") %>%
  posterior::summarize_draws()
#> # A tibble: 8 × 10
#>   variable    mean  median     sd    mad     q5   q95  rhat ess_bulk ess_tail
#>   <chr>      <dbl>   <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
#> 1 THETA[1]  0.207   0.239  0.155  0.150  -0.128 0.406 1.72      3.72     24.2
#> 2 THETA[2]  3.97    4.00   0.0892 0.0809  3.80  4.07  1.68      3.80     18.8
#> 3 THETA[3]  1.16    1.17   0.0271 0.0281  1.12  1.20  0.998    81.1     116. 
#> 4 THETA[4]  4.31    4.30   0.0642 0.0594  4.24  4.44  1.45      4.57     24.2
#> 5 THETA[5]  1.54    1.53   0.101  0.119   1.40  1.73  2.04      3.24     13.3
#> 6 THETA[6]  0.503   0.506  0.0609 0.0671  0.408 0.598 1.03     30.8      24.8
#> 7 THETA[7] -0.0307 -0.0416 0.0853 0.0727 -0.177 0.108 1.04     55.1      79.9
#> 8 THETA[8]  0.392   0.387  0.0921 0.0908  0.231 0.551 1.05     41.1      78.4

`nm_join_bayes()`

nm_join_bayes() and nm_join_bayes_quick() are functions for joining model output summaries to the input data, in a similar way that bbr::nm_join does for non-Bayesian NONMEM models.

nm_join_bayes() selects a subset of the posterior samples. It uses these as input to a user-defined mrgsolve model to simulate EPRED and IPRED and then summarizes the simulated values in the result. In addition, it feeds the simulated EPRED values to npde to calculate EWRES and NPDE values.

nm_join_bayes_quick() does not simulate quantities from the posterior samples. It instead summarizes the reported values from the table files. As the name suggests, this is useful for taking an initial look at the results, though simulation-based values should be used for more reliable estimates.

nm_join_bayes_quick(mod1100) %>%
  select(NUM, TIME, EPRED, IPRED, NPDE, EWRES)
#> # A tibble: 4,292 × 6
#>      NUM  TIME EPRED IPRED   NPDE  EWRES
#>    <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>
#>  1     1  0      0     0    0      0    
#>  2     2  0.61  59.3  65.7 -0.417 -0.516
#>  3     3  1.15  77.8  90.0  0.189  0.128
#>  4     4  1.73  83.9  98.6  1.42   1.53 
#>  5     5  2.15  84.0  99.0  1.75   1.75 
#>  6     6  3.19  78.3  91.2 -0.130 -0.220
#>  7     7  4.21  69.9  80.4 -1.23  -1.10 
#>  8     8  5.09  62.7  71.5 -1.81  -1.57 
#>  9     9  6.22  54.5  61.6  0.544  0.376
#> 10    10  8.09  43.8  49.2  1.43   1.30 
#> # ℹ 4,282 more rows

Accessing other information

If you need access to output files in the chain submodels, you can use the chain_paths() helper.

lst_files <- chain_paths(mod1100, extension = ".lst")
tab_files <- chain_paths(
  mod1100,
  name = get_model_id(mod1100), extension = ".tab"
)

fs::path_rel(lst_files, get_output_dir(mod1100))
#> 1100-1/1100-1.lst 1100-2/1100-2.lst

file.size(lst_files)
#> [1] 37211 37211

fs::path_rel(tab_files, get_output_dir(mod1100))
#> 1100-1/1100.tab 1100-2/1100.tab

file.size(tab_files)
#> [1] 983103 983103

`model_summary()`

We have not yet implemented model_summary() for bbi_nmbayes_model objects. This is a work in progress, but the ?posterior::draws_array object returned by read_fit_model() has various methods for summarizing the outputs.

Run log and friends

run_log(), config_log(), and summary_log() will all find and list NONMEM Bayes models.

run_log(MODEL_DIR) %>%
  select(run, model_type, tags, notes) %>%
  collapse_to_string(tags, notes)
#> # A tibble: 4 × 4
#>   run   model_type tags                                notes                   
#>   <chr> <chr>      <chr>                               <chr>                   
#> 1 100   nonmem     two-compartment + absorption, FOCE  NA                      
#> 2 1000  nmbayes    NA                                  NA                      
#> 3 1100  nmbayes    two-compartment + absorption, BAYES NA                      
#> 4 1101  nmbayes    NA                                  Switched from foo to bar

Note, however, that summary_log() does not currently include include any useful summary statistics or diagnostics for NONMEM Bayes models. This will be revisited once model_summary() is implemented (see above).