Using Model Builder in RStudio • Certara.RsNLME.ModelBuilder

The purpose of this vignette is to demonstrate how to utilize the Certara.RsNLME.ModelBuilder R package to build and parameterize a NLME model from the Shiny GUI. In the following examples, you will learn how to:

Use the modelBuilderUI() function to launch the Shiny GUI and create a built-in model involving covariates and BQL
Update values for fixed effects from the Shiny GUI using the estimatesUI() function from the Certara.RsNLME package
Fit the model in RStudio
Copy and update the model with initial values set to be final parameter estimates
Use the modelTextual() function to directly edit PML statements from the Shiny GUI and return the updated model to RStudio
Fit the updated model

Data

We will be using the built-in pkcovbqlData from the Certara.RsNLME package.

library(Certara.RsNLME)

pkcovbqlData <- Certara.RsNLME::pkcovbqlData

head(pkcovbqlData)
#> # A tibble: 6 × 8
#>      ID  Time  Dose  CObs  LLOQ CObsBQL    BW   PMA
#>   <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl>
#> 1     1   0    1200  NA    0.01       0    12     1
#> 2     1   0.5    NA  97.7  0.01       0    12     1
#> 3     1   1      NA  97.4  0.01       0    12     1
#> 4     1   2      NA  79.5  0.01       0    12     1
#> 5     1   4      NA  96.4  0.01       0    12     1
#> 6     1   8      NA  87.2  0.01       0    12     1

Launch Model Builder (built-in model)

We’ll use the modelBuilderUI() function to launch the Shiny GUI and build our model. Note that we will provide the pkcovbqlData data.frame from above to the data argument, and supply a string to the modelName argument, which is used to create the run folder for the associated output files from model execution.

library(Certara.RsNLME.ModelBuilder)

baseModel <- modelBuilderUI(data = pkcovbqlData, modelName = "baseModel")

Specify BQL

We will leave all options in the “PK Structural Model” section as default.

In the “Residual Error Model” selection, we will select the “BQL checkbox”. This will add the bql indicator to line 6 in our PML e.g., observe(CObs=C * ( 1 + CEps), bql).

Add covariates

Use the next button or the hamburger menu in the top left corner of the page to navigate to the “Parameters” step, where you will see options related to “Structural Parameters”, “Covariates”, “Fixed Effects”, and “Random Effects”.

From the “Covariates” tab, we will introduce our two covariates into the model by clicking the “+” button and specifiying the covariate type and whether to introduce the covariate effect on a structural parameter.

For our continuous covariate body weight, denoted as BW, we will center the covariate on a value of 30 (or alternatively, select “mean” or “median” to use as the centering value), and set the covariate effect on our structural parameters, V and Cl.

For our second (continuous) covariate PMA, we will add to our model without explicitly specifying the covariate effect at this time (which will add the fcovariate(PMA) statement into PML code without changing the stparm() statement).

After adding the covariates into our model from the GUI, we can observe the corresponding updates to the PML on lines 7-10.

Update covariance matrices for random effects

Next, from the “Random Effects” tab, select from the available ETA to update one or more covariance matrices. We will specify values of 0.1 for nV and 0.2 for nCl.

Map model variables to columns

We will skip the input options step and navigate to the final step in the Shiny GUI, “Column Mapping”. Here we will map all required model variables to their corresponding columns in the input data. Columns denoted with * are required.

Clicking the “PREVIEW” button (not required) will allow you to inspect the corresponding column definition file used by NLME during execution, ensuring all required model variables have been mapped correctly.

Generate RsNLME code

Whether you have finished building your model, or still in progress, at any time you may navigate to the “RsNLME” tab and click “Generate” to produce the corresponding RsNLME R code given model-building operations you’ve performed in the GUI. This will allow you to recreate the model from the command line in future R sessions, or, share with other stakeholders to reproduce results.

Copy the RsNLME code generated in the Shiny GUI to an .R and/or .Rmd script, where you can recreate and update the model from the R command line.

library(Certara.RsNLME)
library(magrittr)

baseModel <- pkmodel(isPopulation = TRUE, parameterization = "Clearance", absorption = "Intravenous", numCompartments = 1, isClosedForm = FALSE, isTlag = FALSE, hasEliminationComp = FALSE, isFractionExcreted = FALSE, isSaturating = FALSE, infusionAllowed = FALSE, isDuration = FALSE, columnMap = FALSE, modelName = "baseModel") %>%
  residualError(predName = "C", errorType = "Multiplicative", SD = 0.1, isFrozen = FALSE, isBQL = TRUE) %>%
  structuralParameter(paramName = "Cl", fixedEffName = "tvCl", randomEffName = "nCl", style = "LogNormal", hasRandomEffect = TRUE) %>%
  structuralParameter(paramName = "V", fixedEffName = "tvV", randomEffName = "nV", style = "LogNormal", hasRandomEffect = TRUE) %>%
  addCovariate(covariate = "BW", effect = c("V", "Cl"), type = "Continuous", direction = "Forward", center = "Value", centerValue = 30L) %>%
  addCovariate(covariate = "PMA", effect = NULL, type = "Continuous", direction = "Forward", center = "Mean") %>%
  randomEffect(effect = c("nV", "nCl"), value = c(0.1, 0.2), isDiagonal = TRUE, isFrozen = FALSE) %>%
  dataMapping(pkcovbqlData) %>%
  colMapping(c(id = "ID", time = "Time", CObs = "CObs", CObsBQL = "CObsBQL", BW = "BW", PMA = "PMA", A1 = "Dose"))

Exit Model Builder

Clicking “SAVE & EXIT” will return the RsNLME model object back to your R environment, ready for execution.

Update Initial Values for Fixed Effects

You may update initial values for fixed effects from the “Fixed Effects” tab in the modelBuilderUI() GUI, or alternatively, use the estimatesUI() function from Certara.RsNLME to visually determine a set of reasonable initial values for fixed effects.

library(Certara.RsNLME)

baseModel <- estimatesUI(baseModel)

We will set the following values for fixed effects given evaluation of fit in corresponding plot.

tvV = 20
tvCl = 8
dVdBW = 1
dCldBW = 1

Click the “EXIT” button then check “Update fixed effects in the model” checkbox to accept the initial estimates and update the model object in your R environment, or, exit without saving to return to R without updating initial estimates for your model.

We can print() the model and see the updated values for our fixed effects.

print(baseModel)
#> 
#>  Model Overview 
#>  ------------------------------------------- 
#> Model Name        :  baseModel
#> Working Directory :  C:/Users/jcraig/Documents/GitHub/R-RsNLME-model-builder/vignettes/baseModel
#> Is population     :  TRUE
#> Model Type        :  PK
#> 
#>  PK 
#>  ------------------------------------------- 
#> Parameterization  :  Clearance
#> Absorption        :  Intravenous
#> Num Compartments  :  1
#> Dose Tlag?        :  FALSE
#> Elimination Comp ?:  FALSE
#> Infusion Allowed ?:  FALSE
#> Sequential        :  FALSE
#> Freeze PK         :  FALSE
#> 
#>  PML 
#>  ------------------------------------------- 
#> test(){
#>     dosepoint(A1)
#>     C = A1 / V
#>     deriv(A1 =  -  Cl * C)
#>     error(CEps=0.1)
#>     observe(CObs=C * ( 1 + CEps), bql)
#>     stparm(V = tvV * ((BW/30)^dVdBW)   * exp(nV))
#>     stparm(Cl = tvCl * ((BW/30)^dCldBW)   * exp(nCl))
#>     fcovariate(BW)
#>     fcovariate(PMA)
#>     fixef( tvV = c(,20,))
#>     fixef( tvCl = c(,8,))
#>     fixef( dVdBW(enable=c(0)) = c(,1,))
#>     fixef( dCldBW(enable=c(1)) = c(,1,))
#>     ranef(diag(nV,nCl) = c(0.1,0.2))
#> }
#> 
#>  Structural Parameters 
#>  ------------------------------------------- 
#>  V Cl
#>  ------------------------------------------- 
#> Observations:
#> Observation Name :  CObs
#> Effect Name      :  C
#> Epsilon Name     :  CEps
#> Epsilon Type     :  Multiplicative
#> Epsilon frozen   :  FALSE
#> is BQL           :  TRUE
#>  ------------------------------------------- 
#>  Column Mappings 
#>  ------------------------------------------- 
#> Model Variable Name : Data Column name
#> id                  : ID
#> time                : Time
#> A1                  : Dose
#> BW                  : BW
#> PMA                 : PMA
#> CObs                : CObs
#> CObsBQL             : CObsBQL

Fit the initial model

Now that we’ve defined appropriate values for our fixed effects, we will fit the model in RStudio using the function fitmodel() from Certara.RsNLME.

library(Certara.RsNLME)

baseJob <- fitmodel(baseModel)
#> Using MPI host with 4 cores
#> 
#> NLME Job
#> sharedDirectory is not given in the host class instance.
#> Valid NLME_ROOT_DIRECTORY is not given, using current working directory:
#> C:/Users/jcraig/Documents/GitHub/R-RsNLME-model-builder/vignettes
#> 
#> Compiling 1 of 1 NLME models
#> The model compiled
#> 
#>  Iteration    -2LL     tvV    tvCl    dVdBW  dCldBW nSubj nObs
#>          1 2642.96 23.8689 8.18440 0.902166 1.28015    80  640
#>          2 2615.21 24.3696 8.19992 0.842820 1.38963    80  640
#>          3 2545.36 26.1442 8.22825 0.560945 1.91162    80  640
#>          4 2494.24 24.3840 8.17353 0.798957 1.92962    80  640
#>          5 2472.69 24.7145 8.08533 0.835758 2.13667    80  640
#>          6 2457.11 24.6638 8.04323 0.781674 2.18437    80  640
#>          7 2455.27 25.0973 8.01639 0.771684 2.17212    80  640
#>          8 2453.25 24.5995 7.91023 0.751087 2.14399    80  640
#>          9 2448.65 24.7030 7.79540 0.741252 2.18770    80  640
#>         10 2441.62 24.9100 7.61455 0.754401 2.21735    80  640
#>         11 2437.08 24.9012 7.33535 0.766621 2.23979    80  640
#>         12 2435.12 24.8008 7.63446 0.766946 2.21612    80  640
#>         13 2434.73 24.7722 7.68220 0.764927 2.22003    80  640
#>         14 2434.71 24.7723 7.71726 0.765024 2.22054    80  640
#>         15 2434.71 24.7724 7.71702 0.765006 2.22041    80  640
#>         16 2434.71 24.7753 7.71617 0.764904 2.21941    80  640
#>         17 2434.71 24.7751 7.71617 0.764903 2.21939    80  640
#> 
#> Trying to generate job results...
#> 
#> Generating Overall.csv
#> Generating EtaEta.csv
#> Generating EtaCov.csv
#>  Generating EtaCovariate.csv
#>  Generating EtaCovariateCat.csv
#> Generating StrCovariate.csv
#> Generating Eta.csv
#>  Generating EtaStacked.csv
#>  Generating bluptable.dat
#> Generating ConvergenceData.csv
#> Generating initest.csv
#> Generating doses.csv
#> Generating omega.csv
#> Generating omega_stderr.csv
#> Generating theta.csv
#>  Generating thetaCorrelation.csv
#>  Generating thetaCovariance.csv
#>  Generating Covariance.csv
#> Generating Residuals.csv
#> Generating posthoc.csv
#> 
#> Finished summarizing results. Transferring data and loading the results...
#> Done generating job results.

print(baseJob$Overall)
#>    Scenario RetCode    LogLik     -2LL      AIC      BIC nParm nObs nSub
#> 1: WorkFlow       2 -1217.353 2434.707 2448.707 2479.937     7  640   80
#>    EpsShrinkage Condition
#> 1:      0.10496   26.4688

Copy model and accept parameter estimates

After fitting our initial model, we’ll create a new model using the copyModel() function and set the argument acceptAllEffects = TRUE to accept the parameter estimates from our model run.

finalModel <- copyModel(baseModel, acceptAllEffects = TRUE, modelName = "finalModel")

Edit/Update Model

Using the modelTextualUI() function from the Certara.RsNLME.ModelBuilder package, we can make changes to the PML code for our model from the Shiny GUI.

Functionality includes:

Real time syntax check for PML code
Add/edit column mapping specifications
Add additional input options such as ADDL and steady state

library(Certara.RsNLME.ModelBuilder)

model <- modelTextualUI(finalModel)

To incorporate the effect of “PMA” on the structural parameter “Cl”, we will change the following statement:

stparm(Cl = tvCl * ((BW/30)^dCldBW) * exp(nCl))

stparm(Cl = tvCl * (BW/30)^dCldBW * (PMA^Gam/(PMA^Gam + PMA50^Gam)) * exp(nCl))

Then add the following statements (right before the ranef statement) to define the newly introduced fixed effects.

fixef(PMA50 = c(, 5, ))

fixef(Gam = c(, 1, ))

Fit final model

After introducing the covariate effect of PMA to our structural parameter, Cl, and adding our additional fixed effects, we will fit the final model in RStudio once more using the function fitmodel() from Certara.RsNLME.

library(Certara.RsNLME)

finalJob <- fitmodel(finalModel)
#> Using MPI host with 4 cores
#> 
#> NLME Job
#> sharedDirectory is not given in the host class instance.
#> Valid NLME_ROOT_DIRECTORY is not given, using current working directory:
#> C:/Users/jcraig/Documents/GitHub/R-RsNLME-model-builder/vignettes
#> 
#> Compiling 1 of 1 NLME models
#> The model compiled
#> 
#>  Iteration    -2LL     tvV    tvCl    dVdBW  dCldBW nSubj nObs
#>          1 2434.71 24.8021 7.72478 0.764872 2.22186    80  640
#>          2 2434.71 24.7966 7.72478 0.764872 2.22181    80  640
#>          3 2434.71 24.7719 7.72457 0.764870 2.22153    80  640
#>          4 2434.71 24.7735 7.72452 0.764875 2.22150    80  640
#>          5 2434.71 24.7739 7.72098 0.764880 2.21893    80  640
#>          6 2434.71 24.7738 7.72095 0.764879 2.21893    80  640
#> 
#> Trying to generate job results...
#> 
#> Generating Overall.csv
#> Generating EtaEta.csv
#> Generating EtaCov.csv
#>  Generating EtaCovariate.csv
#>  Generating EtaCovariateCat.csv
#> Generating StrCovariate.csv
#> Generating Eta.csv
#>  Generating EtaStacked.csv
#>  Generating bluptable.dat
#> Generating ConvergenceData.csv
#> Generating initest.csv
#> Generating doses.csv
#> Generating omega.csv
#> Generating omega_stderr.csv
#> Generating theta.csv
#>  Generating thetaCorrelation.csv
#>  Generating thetaCovariance.csv
#>  Generating Covariance.csv
#> Generating Residuals.csv
#> Generating posthoc.csv
#> 
#> Finished summarizing results. Transferring data and loading the results...
#> Done generating job results.

print(finalJob$Overall)
#>    Scenario RetCode    LogLik     -2LL      AIC      BIC nParm nObs nSub
#> 1: WorkFlow       2 -1217.353 2434.706 2448.706 2479.937     7  640   80
#>    EpsShrinkage Condition
#> 1:      0.10487  26.44691