Using Model Builder in RStudio
modelbuilder_rstudio.Rmd
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 theCertara.RsNLMEpackage - 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 1Launch 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 : CObsBQLFit 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)
#> MPI not found on the system.
#> Using localhost without parallelization.
#>
#> 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
#> TDL5 version: 23.10.1.283
#> Servername: no-net
#> Sentinel License file: C:\PROGRA~1\Certara\NLME_E~1\lservrc
#>
#> License Name: NU_PMLCommandLine_N
#> License Type: Commercial
#> Expiration Date: 31 Jan 2024
#> Current Date: 17 Nov 2023
#> Days until program expires: 75
#> 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.3689 8.19965 0.842822 1.38964 80 640
#> 3 2545.31 26.1404 8.22741 0.560995 1.91161 80 640
#> 4 2494.21 24.3861 8.17291 0.799068 1.92974 80 640
#> 5 2472.67 24.7155 8.08545 0.835674 2.13684 80 640
#> 6 2457.11 24.6638 8.04369 0.781643 2.18438 80 640
#> 7 2455.28 25.0976 8.01678 0.771658 2.17206 80 640
#> 8 2453.26 24.5998 7.90991 0.751189 2.14356 80 640
#> 9 2448.66 24.7020 7.79512 0.741226 2.18738 80 640
#> 10 2441.66 24.9085 7.61638 0.754266 2.21733 80 640
#> 11 2437.08 24.9018 7.33681 0.766579 2.24009 80 640
#> 12 2435.11 24.8010 7.63493 0.766938 2.21619 80 640
#> 13 2434.73 24.7721 7.68259 0.764938 2.22005 80 640
#> 14 2434.71 24.7723 7.71752 0.765022 2.22053 80 640
#> 15 2434.71 24.7724 7.71725 0.765006 2.22040 80 640
#> 16 2434.71 24.7724 7.71714 0.765000 2.22038 80 640
#> 17 2434.71 24.7748 7.71650 0.764909 2.21936 80 640
#> 18 2434.71 24.7784 7.71702 0.764971 2.21954 80 640
#> 19 2434.71 24.7765 7.71832 0.765035 2.21934 80 640
#> 20 2434.71 24.7772 7.72025 0.765100 2.21947 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 1 -1217.353 2434.706 2448.706 2479.936 7 640 80
#> EpsShrinkage Condition
#> 1: 0.10531 26.51068Copy 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))
to
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)
#> MPI not found on the system.
#> Using localhost without parallelization.
#>
#> 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
#> TDL5 version: 23.10.1.283
#> Servername: no-net
#> Sentinel License file: C:\PROGRA~1\Certara\NLME_E~1\lservrc
#>
#> License Name: NU_PMLCommandLine_N
#> License Type: Commercial
#> Expiration Date: 31 Jan 2024
#> Current Date: 17 Nov 2023
#> Days until program expires: 75
#> The model compiled
#>
#>
#> Iteration -2LL tvV tvCl dVdBW dCldBW nSubj nObs
#> 1 2434.72 24.8045 7.72901 0.765058 2.22198 80 640
#> 2 2434.71 24.7987 7.72898 0.765049 2.22193 80 640
#> 3 2434.71 24.7712 7.72856 0.764994 2.22161 80 640
#> 4 2434.71 24.7736 7.72849 0.764997 2.22158 80 640
#> 5 2434.71 24.7738 7.72490 0.764805 2.21975 80 640
#> 6 2434.71 24.7725 7.72284 0.764778 2.21914 80 640
#> 7 2434.71 24.7796 7.72053 0.765017 2.21894 80 640
#> 8 2434.71 24.7765 7.72029 0.765056 2.21907 80 640
#> 9 2434.71 24.7761 7.72035 0.765022 2.21912 80 640
#> 10 2434.71 24.7761 7.72035 0.765021 2.21912 80 640
#> 11 2434.71 24.7761 7.72034 0.765021 2.21912 80 640
#> 12 2434.71 24.7761 7.72034 0.765021 2.21912 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 3 -1217.353 2434.706 2448.706 2479.936 7 640 80
#> EpsShrinkage Condition
#> 1: 0.10507 26.36589