The core nm object is created via the new_nm(). Subsequent child objects are created via the child() function. All interactions with NONMEM occur through the this object. It contains metadata about the NONMEM run and contain the contents of what will written as the control file (also called model file).
NONMEM control files are only written to the Models sub-directory just before running NONMEM via the run_nm() command
Let’s create our first NONMEM object with new_nm(). This will be a parent run to all other runs and thus requires more set up than other runs. Subsequent child objects created with child() will inherits characteristics of the parent.
Three arguments are required to create the parent, a run identifier, run_id, a control file it’s based on, based_on, and a (relative) dataset path data_path:
m1 <- new_nm(run_id = "m1",
based_on = "staging/Models/ADVAN2.mod",
data_path = "DerivedData/data.csv")
m1
#> List of 1
#> $ execute.Models/m1
#> ..$ type : chr "execute"
#> ..$ run_id : chr "m1"
#> ..$ run_in : chr "Models"
#> ..$ executed : logi FALSE
#> ..$ ctl_contents: chr collapsed - view with ctl_contents()
#> ..$ ctl_orig : chr collapsed - view with ctl_orig()
#> ..$ data_path : chr "DerivedData/data.csv"
#> ..$ cmd : chr "execute runm1.mod -dir=m1"
#> ..$ cores : int 1
#> ..$ parafile : chr "path/to/parafile.pnm"
#> ..$ run_dir : chr "m1"
#> ..$ ctl_name : chr "runm1.mod"
#> ..$ results_dir : chr "Results"
#> ..$ unique_id : chr "execute.Models/m1"
#> ..$ lst_path : chr "m1/NM_run1/psn.lst"
#> ..- attr(*, "class")= chr [1:3] "nm_execute" "nm_generic" "list"
#> - attr(*, "class")= chr [1:2] "nm_list" "list"Display the object by typing m1 in the console. Notice that the run_in field is point to Models. This can be changed by piping the function of the same name e.g.
m1 <- new_nm(run_id = "m1",
based_on = "staging/Models/ADVAN2.mod",
data_path = "DerivedData/data.csv") %>%
run_in("NONMEM/base_model")will run all models in a subdirectory of a subdirectory NONMEM/base_model/ instead (the $DATA (relative) path to the dataset will be automatically be updated to reflect the new location of the control file). Any field can by modified using a similar heuristic.
Piping is encouraged because it enables you to read the sequence of events in the order that they occur. Here we just have a single pipe, but you’ll see that we will frequently pipe longer chains of commands together. Each step is applying a transformation to the core nm object. The chain can be partially run in the console to see what each transformation is doing. The rest of this document is full of pipes so if you are unfamiliar with pipes, please consult the magrittr documentation.
To extract a field (rather than set a field) use the same function without additional arguments:
run_in(m1)
#> [1] "NONMEM/base_model"For now though, we’ll remove the last piped command and stay with the default run_in location.
NOTE: the field ctl_name refers to the name of the control that will be created. This will only be created when it the model is run with the run_nm() function (described later). For now, the control file contents reside inside the object. To view these you can use “show model/ctl file” RStudio ‘Addin’ , or text(m1).
A few automatic edits from the staged control file and a compact representation of these changes can be shown by highlighting the above code and selecting the “nm_diff” RStudio ‘Addin’ which show what has been changed.
Learning how to read to diffs will be an important skill in NMproject you will pick up over time. Notice how the $DATA has been updated to refer to the new location.
The default cmd() field of the object is execute {ctl_name} -dir={run_dir}. The braces are referred to as glue fields using the glue package. These refer to field names of the object that will be substituted in. For completeness on the next step we will explicitly set this to ensure our model development is easy to read.
m1 <- new_nm(run_id = "m1",
based_on = "staging/Models/ADVAN2.mod",
data_path = "DerivedData/data.csv") %>%
cmd("execute {ctl_name} -dir={run_dir}")The final steps to gets the NONMEM model ready is to fill in the remaining blanks in the template. They are the $INPUT and $THETA, $OMEGA. For this we will use the fill_input() and init_theta() and init_omega().
m1 <- new_nm(run_id = "m1",
based_on = "staging/Models/ADVAN2.mod",
data_path = "DerivedData/data.csv") %>%
cmd("execute {ctl_name} -dir={run_dir}") %>%
fill_input() %>%
init_theta(init = c(-2, 0.5, 1)) %>%
init_sigma(init = c(0.1, 0.1))Thus far, we have not executed NONMEM nor saved the control file to the file system. To execute, we simply pipe into the run_nm() function which will often form the last step of the chain and then run the command. For more information about executing NONMEM, performing NMTRAN checks and monitoring convergence see the Execution article.
m1 <- new_nm(run_id = "m1",
based_on = "staging/Models/ADVAN2.mod",
data_path = "DerivedData/data.csv") %>%
cmd("execute {ctl_name} -dir={run_dir}") %>%
fill_input() %>%
init_theta(init = c(-2, 0.5, 1)) %>%
init_sigma(init = c(0.1, 0.1)) %>%
run_nm()Often there will not be the functions to do the control file manipulation you want. Although it is preferable to stick to automatic control file manipulation functions, you can do fully tracked manual edits via the “manual edit” RStudio ‘Addins’ menu. Again, just highlight the code, click the ‘Addin’ and follow the instructions:
NMproject contains several functions for automatic control file edits. We have already seen fill_input() and init_theta() etc. There are higher order functions which make multiple changes to your control stream, one of which is the subroutine(). If we already have a parent run m1 using ADVAN2 TRANS1, we can create a child() run that uses TRANS2 using:
View exactly what’s been changed by highlighting the above code in RStudio and clicking the "nm_diff RStudio ‘Addin’ to see what’s been changed before running. Here, changes will be in the $SUB, $PK and $THETA.
To add a covariate using PsN coding conventions use add_cov():
We saw a little at how the functions init_theta, init_omega and init_sigma can be used to set initial estimates earlier. The functions actually allow manipulation of initial estimates, parameter bounds, names, units, etc.
## display R representaton of $THETA
it <- m1 %>% init_theta() %>% dplyr::first()
it
#> name parameter lower init upper theta FIX unit trans line pos orig_line
#> 2 KA THETA1 NA 0.5 NA 1 FALSE h-1 LOG 2 1 2
#> 3 K THETA2 NA -2.5 NA 2 FALSE h-1 LOG 3 1 3
#> 4 V THETA3 NA -0.5 NA 3 FALSE L LOG 4 1 4Note that the use of dplyr::first() is because init_theta returns a list of a data.frame and since we want to manipulate io, it’s easier if it’s a data.frame or tibble as then we can use R extensive data.frame manipulation functions. We use the same init_theta with the modified tibble to update the nm object like so:
## this is the slower method - only for illustration purposes
it$init <- c(0, 1, 2)
m1 <- m1 %>% init_theta(it)
m1 %>% dollar("THETA")
#> $`execute.Models/m1`
#> 1| $THETA
#> 2| 0 ; KA ; h-1 ; LOG
#> 3| 1 ; K ; h-1 ; LOG
#> 4| 2 ; V ; L ; LOGThis is quite inconvenient though as it requires 3 lines of R code. Knowing the columns of the ‘tibble’ though, we can manipulate values directly like so:
## Reset in one line to what we set it to earlier
m1 <- m1 %>% init_theta(init = c(-2, 0.5, 1))
m1 %>% dollar("THETA")
#> $`execute.Models/m1`
#> 1| $THETA
#> 2| -2 ; KA ; h-1 ; LOG
#> 3| 0.5 ; K ; h-1 ; LOG
#> 4| 1 ; V ; L ; LOG
## This however requires knowing the order of parameters
## here we supply a named vector in a different order
m1 <- m1 %>% init_theta(init = c(KA = -2, V = 1))
m1 %>% dollar("THETA")
#> $`execute.Models/m1`
#> 1| $THETA
#> 2| -2 ; KA ; h-1 ; LOG
#> 3| 0.5 ; K ; h-1 ; LOG
#> 4| 1 ; V ; L ; LOG
## can also manipulate other aspects (like the FIX column) similarly
m1 <- m1 %>% init_theta(init = c(KA = -2, V = 1),
FIX = c(KA = TRUE))
m1 %>% dollar("THETA")
#> $`execute.Models/m1`
#> 1| $THETA
#> 2| -2 FIX ; KA ; h-1 ; LOG
#> 3| 0.5 ; K ; h-1 ; LOG
#> 4| 1 ; V ; L ; LOGIt works similarly for $OMEGA and $SIGMA with init_omega and init_sigma, respectively.
To modify initial estimates, we’ll use the mutate like behaviour of init_* functions. We will modify the init by referencing itself. We’ll modified all our fixed effects (log transformed) by 30%
m1 <- m1 %>% init_theta(init = rnorm(init, mean = init, sd = 0.3))
m1 %>% dollar("THETA")
#> $`execute.Models/m1`
#> 1| $THETA
#> 2| -2.3115 FIX ; KA ; h-1 ; LOG
#> 3| 0.6613 ; K ; h-1 ; LOG
#> 4| 1.2536 ; V ; L ; LOG
m1 <- m1 %>% init_omega(init = runif(init, min = init/2, max = init*2))
m1 %>% dollar("OMEGA")
#> $`execute.Models/m1`
#> 1| $OMEGA
#> 2| 0.071963 ; IIV_KA ; LOG
#> 3| 0.19997 ; IIV_K ; LOG
#> 4| 0.18206 ; IIV_V ; LOGWe can include and remove $OMEGA blocks with the functions block and unblock (to create $OMEGA BLOCKS for correlated random effects).
io <- m1 %>% init_omega() ## note we dont need dplyr::first as block()/unblock() also work on lists of tibbles.
io
#> $`execute.Models/m1`
#> name omega1 omega2 lower init upper block mblock FIX unit trans line
#> 1 <NA> NA NA NA NA NA 1 NA FALSE <NA> <NA> 1
#> 2 IIV_KA 1 1 NA 0.071963 NA 1 1 FALSE <NA> LOG 2
#> 3 IIV_K 2 2 NA 0.199970 NA 2 1 FALSE <NA> LOG 3
#> 4 IIV_V 3 3 NA 0.182060 NA 3 1 FALSE <NA> LOG 4
#> pos orig_line orig_pos
#> 1 1 1 1
#> 2 1 2 1
#> 3 1 3 1
#> 4 1 4 1
io <- io %>% block(c(2,3)) ## make block out ETA 2 and 3
## put modified io wit
m1 <- m1 %>% init_omega(io)
m1 %>% dollar("OMEGA")
#> $`execute.Models/m1`
#> 1| $OMEGA
#> 2| 0.071963 ; IIV_KA ; LOG
#> 3| $OMEGA BLOCK (2)
#> 4| 0.19997 ; IIV_K ; LOG
#> 5| 0.01 0.18206 ; IIV_V ; LOG
## for demo purposes we'll reverse the process with unblock()
io <- m1 %>% init_omega()
io <- io %>% unblock(c(2,3))
m1 <- m1 %>% init_omega(io)
m1 %>% dollar("OMEGA")
#> $`execute.Models/m1`
#> 1| $OMEGA
#> 2| 0.071963 ; IIV_KA ; LOG
#> 3| 0.19997 ; IIV_K ; LOG
#> 4| 0.18206 ; IIV_V ; LOG