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:
modelBuilderUI()
function to launch the Shiny
GUI and create a built-in
model involving covariates and BQLestimatesUI()
function from the Certara.RsNLME
packagemodelTextual()
function to directly edit PML
statements from the Shiny GUI and return the updated model to
RStudioWe 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
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")
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)
.
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.
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
.
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.
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"))
Clicking “SAVE & EXIT” will return the RsNLME model object back to your R environment, ready for execution.
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
= 20tvCl
= 8dVdBW
= 1dCldBW
= 1Click 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
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
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")
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:
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, ))
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