In MikeKSmith/rspeaksnonmem: Read, Modify and run NONMEM from within R

knitr::opts_chunk$set(fig.width=6, fig.height=6)

What does rspeaksnonmem do?

rspeaksnonmem is designed to allow the user to craft workflows based on a given NONMEM model. It is not intended to replace Perl speaks Nonmem (PsN) - in fact rspeaksnonmem can call and use PsN functionality from an R script. In fact rspeaksnonmem could equally be called "rspeaksPsN". rspeaksnonmem helps the user by allow them to write complex workflows using NONMEM, PsN and R. For example performing data checkout (exploratory analysis), using a "template" model as a starting point for model refinement, specifying and running a sequence of tasks for every model, collating results across models for model comparison.

After importing and parsing a control stream to an R object using the RNMImport package, rspeaksnonmem allows the user to easily change initial estimates, change data attributes or change task properties (estimation or simulation settings) without having to change the model. This saves a lot of textual searching by the user.

To best use rspeaksnonmem it is useful to think of the NMTRAN expression of the model as a combination of: Data specification Parameter specification Model specification (structural, covariate and stochastic components) Task specification

With rspeaksnonmem it is easy to alter data, parameter and task information, but altering the model specification has knock-on effects across many other aspects. We recommend that the user sets up the model in such a way that it is easy to change the model simply by fixing or estimating certain parameters i.e. fixing population parameters and/or altering the OMEGA and SIGMA parameter specification.

Using a template model and a user-specified workflow function e.g. a function combining execute_PsN + sumo_PsN + basic_GOF the user could specify a number of combinations of THETA, OMEGA and SIGMA that would test a range of plausible models, run these models and then compare the output (OFV / AIC and model diagnostics) to find the best model fit. This is suggested not as a way of finding the best final model, but the best base model for further refinement. There is no substitute for the knowledge and skill of an analyst in building and assessing models. We refer the reader to Schmidt and Radivojevic (2014), JPKPD (http://www.ncbi.nlm.nih.gov/pubmed/25056507).

Installation

Dependencies

rspeaksnonmem relies on the package RNMImport function importNmMod which reads and parses the NONMEM control stream. rspeaksnonmem then works with the data, parameter values, and task information separately from the model.

Before installation of rspeaksnonmem, the package RNMImport needs to be installed.

Install RNMIMport from Github using the devtools package. rspeaksnonmem requires the version of RNMImport on MikeKSmith's Github repository which includes handling of OMEGA and SIGMA blocks and identifying which elements are fixed for these parameters.

  #devtools::install_github("MikeKSmith/RNMImport")

Install rspeaksnonmem

Eventually, rspeaksnonmem will be released to CRAN, but while still in development rspeaksnonmem can most easily be installed from GitHub using the devtools package:

    #devtools::install_github("MikeKSmith/rspeaksnonmem")

Alternatively, if the source files are available locally, you can load the functions directly.

 devtools::load_all(pkg = ".")
 devtools::load_all(pkg = "C:\\Users\\smith_mk\\Documents\\Working documents\\RNMImport")

#library(rspeaksnonmem)

Copy an example dataset and model to a directory of your choice

getwd()

file.copy(from = file.path(system.file("exdata", package = "rspeaksnonmem"),
                           "warfarin_conc_pca.csv"),
            to = getwd(), overwrite = T )
file.copy(from = file.path(system.file("exdata", package = "rspeaksnonmem"),
                           "warfarin_bootstrap50.dta"),
            to = getwd(), overwrite = T )
file.copy(from = file.path(system.file("exdata", package = "rspeaksnonmem"),
                           "warfarin.ctl"),
            to = getwd(), overwrite = T )

The initial model within the workflow should act as a "template" for modifications. The best practice is to use a model where all possible parameters are defined (including OMEGAs and covariances / correlations) but where it is possible to fix parameters to zero or some null value. We can then run and test a wide variety of models simply by updating the $THETA, $OMEGA and $SIGMA parameters to allow estimation.

Read the control stream using RNMImport

First, we need to read the NONMEM control stream into R using the importNmMod function of RNMImport.

warfModel <- importNmMod("warfarin.ctl")
class(warfModel)
names(warfModel)

This function reads NONMEM control streams and parses the code. The initial R object contains the raw code as a vector of character strings, any comments from the NONMEM code (anything after ";" in the NMTRAN), the name of the control file in controlFile and the parsed code within the list problemContents. RNMImport allows for cases where there is more than one NONMEM $PROBLEM statement within the code.

We can select the parsed control stream as the basis for modification. Here we take the first / only $PROBLEM statement content.

warfTemplateModel <- warfModel$problemContents[[1]]
str(warfTemplateModel)

Now with the parsed content, it is much easier to change elements within the NONMEM control stream since we're only manipulating text strings. But since the code has been parsed it's more transparent to a third party exactly what is being changed.

First let's change the dataset. Here I've sampled 50 subjects from the original warfarin dataset (warfarin.csv) into a dataset called "warfarin_bootstrap50.dta". Since the names and uses of the data columns is unchanged, we do not need to change anything in the $INPUT statement of NMTRAN. So we create a new object called newData containing the original code, then update the file name.

newData <- warfModel$problemContents[[1]]$Data
newData[,"File"] <- "warfarin_bootstrap50.dta"

In a similar way, we can then update any element of the model using this parsed set of commands. rspeaksnonmem provides some additional functions to extract certain elements of the parsed control stream.

"Object" | rspeaksnonmem Function | NMTRAN blocks ---------|---------------------------|----------------------------- Data | getNMDataObjects | $DATA, $INPUT Parameter| getNMParameterObjects | $THETA, $OMEGA, $SIGMA Task | getNMTaskPropertiesObjects| $EST,$COV, $TAB Model | getNMModelObjects | everything else

For example: getNMParameterObjects returns the $THETA, $OMEGA and $SIGMA records. Since the $THETA parameters are represented as a data frame, we can easily update the initial values, lower or upper bounds, fix or estimate (unfix) THETA parameters. With each object returned by the various getNM<...> functions we get not only the parsed objects but the raw control stream as well.

params <- getNMParameterObjects(warfModel)
str(params)

Together the Data, Parameters, Model and Task "objects" form a "Modelling Object Group" or MOG which is used for a specific estimation task.

In our example, the template model has POP_TLAG fixed to zero. To "unfix" this, we need to provide lower, initial (Est) and upper bounds for POP_TLAG. Note that params$THETA has named rows using the comments after each THETA line in the NMTRAN control stream.

We can create multiple "objects" corresponding to different models that we wish to explore. So below we create an object called thetaNoLag where the $THETA for the LAG term is fixed to zero. We also create an object called "thetaLag" where we set FIX = FALSE and give a reasonable initial value for the Lag.

FIX_POP_TLAG <- list(Lower = 0, Est = 0, Upper = 0, FIX = TRUE, comments = "POP_TLAG")
EST_POP_TLAG <- list(Lower = 0, Est = 0.75, Upper = 1.5, FIX = FALSE, comments = "POP_TLAG")

thetaNoLag <- warfTemplateModel$Theta
thetaNoLag["POP_TLAG",] <- FIX_POP_TLAG

thetaLag <- warfTemplateModel$Theta
thetaLag["POP_TLAG",] <- EST_POP_TLAG

The template model has a combined proportional and additive residual error structure. This has been parameterised as W = THETA(5)+THETA(6)*CONC in the NMTRAN code. W then multiplies a standard Normal variate to give the residual error: W * SIGMA where SIGMA ~ N(0,1).

To be able to estimate different residual error models, we can then simply fix one or other component of the residual error model. Fixing THETA(5) to zero will estimate with only proportional error, while fixing THETA(6) to zero will estimate with only additive error. Estimating both provides the combined additive and proportional residual error.

Similarly to the above, we can then create a new object which estimates THETA(5) which we will give the name "RUV_ADD" for Residual Unexplained Variability ADDitive component.

ruvPropOnly <- list(Est = 0, FIX = TRUE, comments = "RUV_ADD")

Again, in a similar way we can update the between subject OMEGA parameter specification to examine a variety of models - where Clearance (CL) and Volume of Distribution (V) are independent / uncorrelated, or where they are correlated. We can also easily fix or estimate between subject variability on parameters such as KA by specifying initial values and ensuring FIX = FALSE.

The template model specifies an OMEGA block (covariance) between ETA_CL and ETA_V.

CLVBlock_noTLAG <- params$OMEGA
diagPPV_noTLAG <- list(data.frame(values=c(0.1,0.1,0.1), 
                            FIX=c(F,F,F), 
                            comments=c("; PPV_CL", "; PPV_V", "; PPV_KA")),
                       list(values = 0, FIX = TRUE))

Bringing all of these different components together, we can examine the fit of a number of basic models, starting only with the template model. We use the updateModel function to update the template model with objects that we have defined or altered.

  run1 <- updateModel(warfTemplateModel, 
                      theta = thetaNoLag, 
                      omega = diagPPV_noTLAG, 
                      data = newData, 
                      runno = 1) ## First model

  run2 <- updateModel(warfTemplateModel, 
                      theta = thetaLag, 
                      omega = diagPPV_noTLAG , 
                      data = newData,   
                      runno = 2)  ## Add POP_TLAG

  run3 <- updateModel(warfTemplateModel, 
                      theta = thetaLag, 
                      omega = CLVBlock_noTLAG, 
                      data = newData,    
                      runno = 3)  ## Change to block CL, V

rspeaksnonmem provides function that enable us to run NONMEM or PsN without leaving R. By combining these, we can define a sequence of steps that we may wish to perform in a script, or define a function that combines workflow steps that we would wish to perform on each model.

myPopPKWorkflow <- function(model){
  controlFile <- paste(deparse(substitute(model)),".ctl",sep = "")
  lstFile <- paste(deparse(substitute(model)),".lst",sep = "")

  workingDir <- deparse(substitute(model))

  writeNMControlStream(templateModel = warfModel$Raw, 
                       parsedControl = model,
                       outputFile = controlFile)

 execute_PsN(installPath = "c:/strawberry/perl",
              version = "4.7.0",
              modelFile = controlFile, 
              lstFile = lstFile,
              clean = 1,
              working.dir = workingDir)

  sumo_PsN(tool = "sumo-3.5.4",
           installPath = "c:/strawberry/perl",
           version = "4.7.0",
           lstFile = lstFile)

  runno <- as.numeric(gsub("[a-z]", "", deparse(substitute(model))))

  ## ----createXpdb----------------------------------------------------------
  xpdb <- xpose4::xpose.data(runno = runno, quiet = T)
  # save(base.xpdb, file='Xpose database.RData')

  ## ----xposeGOF------------------------------------------------------------
  print(xpose4::dv.vs.pred.ipred(xpdb))
  print(xpose4::pred.vs.idv(xpdb))
  print(xpose4::ipred.vs.idv(xpdb))
  print(xpose4::wres.vs.idv(xpdb))
  print(xpose4::wres.vs.pred(xpdb))
  print(xpose4::ranpar.hist(xpdb))
  print(xpose4::ind.plots(xpdb, layout = c(4, 4)))  
}

We can then run the workflow for each model in turn or "apply" the function across a list of models to be evaluated.

myPopPKWorkflow(run1)
myPopPKWorkflow(run2)
myPopPKWorkflow(run3)

MikeKSmith/rspeaksnonmem documentation built on March 12, 2023, 3:25 p.m.