Overview

The purpose of this vignette is to demonstrate how to create an absorption delay model through built-in model library, fit the model, create commonly used diagnostic plots through both command-line and mode results shiny app (from Certara.ModelResults package), perform VPC and create VPC plots through tidyvpc package (command-line) or VPC results shiny app from Certara.VPCResults package.

The model demonstrated is a one-compartment model with 1st-order clearance and the delay time between the administration time of the drug and the time when the drug molecules reach the central compartment is assumed to be gamma distributed. In other words, the model is described as follows: \[\frac{dA1(t)}{dt} = \sum_{i=1}^{m}D_{i} gammaPDF(t - t_{Di}; rateParameter, shapeParameter) - Cl/V * A1(t)\]

where

  • \(A1\) denotes the drug amount at central compartment with \(V\) and \(Cl\) respectively being the central volume distribution and central clearance rate
  • \(m\) denotes the number of bolus dosing events
  • \(D_{i}\) is the amount of dose administered at time \(t_{Di}\) for the ith dosing events
  • gammaPDF denotes the probability density function of a gamma distribution with with \(mean = MeanDelayTime\) and shape parameter being \((ShapeParamMinusOne + 1)\)

(For more information on distribute delays, please visit here).

For the rest of this document, we assume that all the necessary packages are loaded and the directory with NLME Executables is given as an environment variable (INSTALLDIR). They can be loaded and checked with the following commands:

  # loading the package
  library(Certara.RsNLME)
  library(data.table)
  library(dplyr)
  library(xpose)
  library(Certara.Xpose.NLME)
  library(ggplot2)
  library(Certara.ModelResults)
  library(Certara.VPCResults)
  library(tidyvpc)

 # Check the environment variable
  Sys.getenv("INSTALLDIR")
[1] "C:\\Program Files\\Certara\\NLME_Engine"

Create the model

We will use the input dataset OneCpt_GammaDistributedAbsorptionDelay.csv included with the Certara.RsNLME package.

  filename <- system.file("vignettesdata/OneCpt_GammaDistributedAbsorptionDelay.csv",
                          package = "Certara.RsNLME",
                          mustWork = TRUE)
  dt_InputDataSet <- fread(filename)

Next, we will define a one-compartment model with gamma absorption delay and 1st-order clearance and make the following modifications to it:

  • Set \(V = e^{tvlogV + nlogV}\)
  • Set \(Cl = e^{tvlogCl + nlogCl}\)
  • Set \(MeanDelayTime = e^{tvlogMeanDelayTime + nlogMeanDelayTime}\)
  • Set \(ShapeParamMinusOne = e^{tvlogShapeParamMinusOne}\)
  • Change initial values for fixed effects \(tvlogV\) and \(tvlogShapeParamMinusOne\) to be 5 and 2, respectively (the default value is 1)
  • Add secondary parameter \(tvV\), which is set to \(tvV = e^{tvlogV}\)
  • Add secondary parameter \(tvCl\), which is set to \(tvCl = e^{tvlogCl}\)
  • Add secondary parameter \(tvMeanDelayTime\), which is set to \(tvMeanDelayTime = e^{tvlogMeanDelayTime}\)
  • Add secondary parameter \(tvShapeParam\), which is set to \(tvShapeParam = e^{tvlogShapeParamMinusOne} + 1\)
  model <- pkmodel(absorption = "Gamma"
                   , data = dt_InputDataSet
                   , columnMap = FALSE
                   , modelName = "OneCpt_GammaDistributedAbsorptionDelay"
                   ) %>%
  structuralParameter(paramName = "V", 
                      style = "LogNormal2", 
                      fixedEffName = "tvlogV", 
                      randomEffName = "nlogV") %>%
  structuralParameter(paramName = "Cl", 
                      style = "LogNormal2", 
                      fixedEffName = "tvlogCl", 
                      randomEffName = "nlogCl") %>%
  structuralParameter(paramName = "MeanDelayTime", 
                      style = "LogNormal2", 
                      fixedEffName = "tvlogMeanDelayTime", 
                      randomEffName = "nlogMeanDelayTime") %>%
  structuralParameter(paramName = "ShapeParamMinusOne", 
                      style = "LogNormal2", 
                      fixedEffName = "tvlogShapeParamMinusOne", 
                      hasRandomEffect = FALSE) %>%
  fixedEffect(effect = c("tvlogV", "tvlogShapeParamMinusOne"), value = c(5, 2)) %>%
  addSecondary(name = "tvV", definition = "exp(tvlogV)") %>%
  addSecondary(name = "tvCl", definition = "exp(tvlogCl)") %>%
  addSecondary(name = "tvMeanDelayTime", definition = "exp(tvlogMeanDelayTime)") %>%
  addSecondary(name = "tvShapeParam", definition = "exp(tvlogShapeParamMinusOne) + 1") 

Let’s view the model and its associated column mappings, and then map those un-mapped model variables to their corresponding input data columns:

  print(model)

 Model Overview 
 ------------------------------------------- 
Is population     :  TRUE
Model Type        :  PK

 PK 
 ------------------------------------------- 
Parameterization  :  Clearance
Absorption        :  Gamma
Num Compartments  :  1
Dose Tlag?        :  FALSE
Elimination Comp ?:  FALSE
Infusion Allowed ?:  FALSE
Sequential        :  FALSE
Freeze PK         :  FALSE

 PML 
 ------------------------------------------- 
test(){
    dosepoint(A1)
    C = A1 / V
    delayInfCpt(A1, MeanDelayTime, ShapeParamMinusOne, out =  -  Cl * C, dist = Gamma)
    error(CEps=0.1)
    observe(CObs=C * ( 1 + CEps))
    stparm(MeanDelayTime = exp(tvlogMeanDelayTime + nlogMeanDelayTime))
    stparm(ShapeParamMinusOne = exp(tvlogShapeParamMinusOne ))
    stparm(V = exp(tvlogV + nlogV))
    stparm(Cl = exp(tvlogCl + nlogCl))
    fixef( tvlogMeanDelayTime = c(,1,))
    fixef( tvlogShapeParamMinusOne = c(,2,))
    fixef( tvlogV = c(,5,))
    fixef( tvlogCl = c(,1,))
    ranef(diag(nlogMeanDelayTime,nlogV,nlogCl) =  c(1,1,1))
    secondary(tvV=exp(tvlogV))
    secondary(tvCl=exp(tvlogCl))
    secondary(tvMeanDelayTime=exp(tvlogMeanDelayTime))
    secondary(tvShapeParam=exp(tvlogShapeParamMinusOne) + 1)
}

 Structural Parameters 
 ------------------------------------------- 
 MeanDelayTime ShapeParamMinusOne V Cl
 ------------------------------------------- 
Observations:
Observation Name :  CObs
Effect Name      :  C
Epsilon Name     :  CEps
Epsilon Type     :  Multiplicative
Epsilon frozen   :  FALSE
is BQL           :  FALSE
 ------------------------------------------- 
 Column Mappings 
 ------------------------------------------- 
Model Variable Name : Data Column name
id                  : ID
time                : Time
A1                  : ?
CObs                : CObs

  ## Manually map those un-mapped model variables to their corresponding input data columns
  model <- model %>% colMapping(c(A1 = "Dose")) 

Model Fitting

Let’s run the model using the fitmodel function with default host and engine method set to be QRPEM.

  job <- fitmodel(model, method = "QRPEM")

Host argument is not given. 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/vignettes

 Compiling 1 of 1 NLME models
The model compiled

 Iteration    -2LL tvlogMeanDelayTime tvlogShapeParamMinusOne  tvlogV tvlogCl
         1 -537.53           0.561276                 1.96611 4.64304 1.07284
 nSubj nObs
    80  720
         2 -1227.75           0.571959                 1.93745 4.64224 1.12164
         3 -1267.67           0.580982                 1.90653 4.63884 1.14450
         4 -1289.29           0.590877                 1.88101 4.63452 1.15754
         5 -1304.35           0.599273                 1.85722 4.63041 1.16602
         6 -1317.09           0.607051                 1.83942 4.62675 1.17206
         7 -1326.00           0.613089                 1.82105 4.62363 1.17653
         8 -1334.32           0.619103                 1.80210 4.62091 1.18000
         9 -1342.02           0.625291                 1.78799 4.61828 1.18351
        10 -1347.54           0.629984                 1.76981 4.61597 1.18647
        11 -1353.35           0.635822                 1.75386 4.61368 1.18932
        12 -1358.09           0.641057                 1.74270 4.61142 1.19220
        13 -1360.84           0.644762                 1.73192 4.60960 1.19441
        14 -1363.42           0.648390                 1.72428 4.60790 1.19656
        15 -1364.65           0.650958                 1.71791 4.60662 1.19805
        16 -1365.92           0.653090                 1.71199 4.60563 1.19938
        17 -1366.78           0.655109                 1.70778 4.60464 1.20055
        18 -1367.49           0.656452                 1.70789 4.60400 1.20139
        19 -1367.20           0.656579                 1.69600 4.60362 1.20177
        20 -1368.94           0.660179                 1.68416 4.60278 1.20290
        21 -1369.63           0.663853                 1.68110 4.60157 1.20452
        22 -1369.79           0.665087                 1.67546 4.60064 1.20569
        23 -1370.17           0.666959                 1.67232 4.59983 1.20695
        24 -1370.03           0.668003                 1.66526 4.59923 1.20765
        25 -1370.50           0.670198                 1.66564 4.59849 1.20855
        26 -1370.24           0.670209                 1.67387 4.59821 1.20904
        27 -1370.04           0.667842                 1.67255 4.59851 1.20860
        28 -1370.06           0.668154                 1.67032 4.59869 1.20848
        29 -1370.20           0.668781                 1.67232 4.59858 1.20858
        30 -1370.43           0.668258                 1.67091 4.59861 1.20837
        31 -1370.21           0.668602                 1.67171 4.59858 1.20834
        32 -1370.25           0.668413                 1.67050 4.59860 1.20830
        33 -1370.34           0.668780                 1.66523 4.59852 1.20842
        34 -1370.21           0.670339                 1.66350 4.59816 1.20897
        35 -1370.19           0.670946                 1.66944 4.59788 1.20930
        36 -1370.38           0.669240                 1.67088 4.59812 1.20918
        37 -1370.31           0.668720                 1.66754 4.59831 1.20897
        38 -1370.23           0.669679                 1.66857 4.59819 1.20905
        39 -1370.10           0.669420                 1.66819 4.59823 1.20899
        40 -1370.44           0.669555                 1.66492 4.59822 1.20898
        41 -1370.21           0.670489                 1.66428 4.59799 1.20906
        42 -1370.14           0.670733                 1.66445 4.59782 1.20943
        43 -1370.56           0.670744                 1.66610 4.59775 1.20946
        44 -1370.11           0.670210                 1.66092 4.59783 1.20952
        45 -1370.12           0.671723                 1.66811 4.59756 1.20963
        46 -1370.48           0.669657                 1.66483 4.59783 1.20934
        47 -1370.15           0.670543                 1.66347 4.59785 1.20939
        48 -1370.58           0.670975                 1.66481 4.59772 1.20970
        49 -1370.38           0.670630                 1.66379 4.59774 1.20963
        50 -1370.19           0.670926                 1.66768 4.59768 1.20962
        51 -1370.34           0.669787                 1.66858 4.59790 1.20935
        52 -1370.27           0.669475                 1.66765 4.59805 1.20943
        53 -1370.48           0.669689                 1.66586 4.59809 1.20934
        54 -1370.17           0.670231                 1.66503 4.59798 1.20939
        55 -1370.26           0.670518                 1.66480 4.59784 1.20958
        56 -1370.21           0.670589                 1.66570 4.59779 1.20956
        57 -1370.47           0.670336                 1.66893 4.59783 1.20947
        58 -1370.21           0.669411                 1.67058 4.59804 1.20914
        59 -1370.31           0.668834                 1.67131 4.59831 1.20874
        60 -1370.12           0.668565                 1.67338 4.59846 1.20862

Trying to generate job results...

 Generating Overall.csv
 Generating StatusWindow.txt
 Generating EtaEta.csv
 Generating Secondary.csv
 Generating ConvergenceData.csv
 Generating initest.csv
 Generating doses.csv
 Generating omega.csv
 Generating omega_stderr.csv
 Generating thetas.csv
 Generating Eta.csv and EtaStacked.csv
 Generating Residuals.csv
 Generating posthoc.csv and table0n.csv
 
Finished summarizing results. Transferring data and loading the results...
Done generating job results.
  print(job)
$Overall
   Scenario RetCode   LogLik      -2LL       AIC       BIC nParm nObs nSub
1: WorkFlow       1 685.0575 -1370.115 -1354.115 -1317.481     8  720   80
   EpsShrinkage Condition
1:      0.14751   2.65545

$theta
   #Scenario               Parameter  Estimate Units      Stderr       CV%
1:  WorkFlow      tvlogMeanDelayTime 0.6685651    NA 0.023438366 3.5057718
2:  WorkFlow tvlogShapeParamMinusOne 1.6733772    NA 0.033084909 1.9771340
3:  WorkFlow                  tvlogV 4.5984615    NA 0.014494277 0.3151984
4:  WorkFlow                 tvlogCl 1.2086231    NA 0.036531295 3.0225548
5:  WorkFlow                    CEps 0.1212948    NA 0.003399118 2.8023597
      2.5%CI   97.5%CI Var.Inf.factor
1: 0.6225485 0.7145817             NA
2: 1.6084215 1.7383328             NA
3: 4.5700049 4.6269182             NA
4: 1.1369011 1.2803450             NA
5: 0.1146213 0.1279683             NA

$omega
               Label nlogMeanDelayTime       nlogV     nlogCl
1: nlogMeanDelayTime         0.0392779 0.000000000 0.00000000
2:             nlogV         0.0000000 0.009778162 0.00000000
3:            nlogCl         0.0000000 0.000000000 0.08210199

$omega_Correlation
               Label nlogMeanDelayTime nlogV nlogCl
1: nlogMeanDelayTime                 1     0      0
2:             nlogV                 0     1      0
3:            nlogCl                 0     0      1

$Eta_Shrinkage
       Label nlogMeanDelayTime     nlogV     nlogCl
1: Shrinkage        0.01704205 0.2077418 0.08547095

$omega_stderr
               Label nlogMeanDelayTime      nlogV     nlogCl
1: nlogMeanDelayTime       0.006427425 0.00000000 0.00000000
2:             nlogV       0.000000000 0.00248433 0.00000000
3:            nlogCl       0.000000000 0.00000000 0.01557804

$Secondary
   Scenario       Secondary  Estimate Units     Stderr      CV%    2.5%CI
1: WorkFlow             tvV 99.331381    NA 1.43973730 1.449428 96.504743
2: WorkFlow            tvCl  3.348870    NA 0.12233863 3.653131  3.108683
3: WorkFlow tvMeanDelayTime  1.951435    NA 0.04573847 2.343838  1.861637
4: WorkFlow    tvShapeParam  6.330138    NA 0.17634723 2.785835  5.983916
      97.5%CI Var.Inf.factor
1: 102.158019             NA
2:   3.589058             NA
3:   2.041234             NA
4:   6.676361             NA
                                                                                     `

Diagnostic Plots

We will use the xposeNlme function from the Certara.Xpose.NLME package to import estimation results to xpose database to create some commonly used diagnostic plots. All the functions provided in the xpose package can be used. Here we only demonstrate several of these functions.

  ## Imports results of an NLME run into xpose database to create commonly used diagnostic plots
  xp <- xposeNlme(dir = model@modelInfo@workingDir, modelName = "OneCpt_GammaDistributedAbsorptionDelay")
  
  ## observations against population predictions 
  dv_vs_pred(xp, type = "p", subtitle = "-2LL: @ofv")


  ## observations against individual predictions 
  dv_vs_ipred(xp, type = "p", subtitle = "-2LL: @ofv, Eps shrinkage: @epsshk")

  
  ## CWRES against population predictions 
  res_vs_pred(xp, res = "CWRES", type = "ps", subtitle = "-2LL: @ofv")
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'


  ## CWRES against the independent variable
  res_vs_idv(xp, res = "CWRES", type = "ps", subtitle = "-2LL: @ofv")
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'

  
  ## Observations, individual predictions and population predictions plotted against the independent variable for every individual   
  ind_plots(xp, subtitle = "-2LL: @ofv, Eps shrinkage: @epsshk")

Alternatively, one can view/customize diagnostic plots as well as estimation results through model results shiny app (from the Certara.ModelResults package), which can also be used to generate R script and report as well as the associated R script and report as well as the associated R markdown. For details on this app as well as how to use it, please visit the following link. Here we only demonstrate how to invoke this shiny app through either model object or the xpose data bases created above.


  ## Invoke model results shiny app through model object defined above 
  resultsUI(model)

  ## Alternatively, one can invoke model results shiny app through xpose data bases created above 
  resultsUI(xpdb = xp)

VPC

We will use the copyModel function to copy model into a new object and accept final parameter estimates from fitting run as initial estimates for VPC simulation:

  modelVPC <- copyModel(model, acceptAllEffects = TRUE, modelName = "OneCpt_GammaDistributedAbsorptionDelay_VPC")

Now, let’s run VPC using the vpcmodel function with the default host, default values for the relevant NLME engine arguments, and default values for VPC arguments:

  vpcJob <- vpcmodel(modelVPC)

Host argument is not given. 
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/vignettes

 Compiling 1 of 1 NLME models
The model compiled

Trying to generate job results...

 
Finished summarizing results. Transferring data and loading the results...
Done generating job results.

VPC/Simulation results are ready in  C:/Users/jcraig/Documents/GitHub/R-RsNLME/vignettes/OneCpt_GammaDistributedAbsorptionDelay_VPC
Loading the results
Loading  predcheck_bql.csv 
Loading  predcheck0.csv 
Loading  predcheck0_cat.csv 
Loading  predcheck1.csv 
Loading  predcheck1_cat.csv 
Loading  predcheck2.csv 
Loading  predcheck2_cat.csv 
Loading  predout.csv 

  ## Observed data
  dt_ObsData <- vpcJob$predcheck0

  ## Simulated data
  dt_SimData <- vpcJob$predout

Next, we will create VPC plots through tidyvpc package. The tidyvpc package provides support for both continuous and categorical VPC using both binning and binless methods. For details on this package, please visit the following link. Note that this example only contains one observed variable and both simulated and observed data meet the requirements set by tidyvpc package. Hence, there is no need to do any data preprocessing for this example.

  ## Create a binless VPC plot 

  plot_binless_vpc <- observed(dt_ObsData, x = IVAR, yobs = DV) %>%
    simulated(dt_SimData, ysim = DV) %>%
    binless() %>%
    vpcstats()
  plot(plot_binless_vpc)

  
  ## Create a binning VPC plot: binning on x-variable itself 
  plot_binning_vpc <- observed(dt_ObsData, x = IVAR, yobs = DV) %>%
    simulated(dt_SimData, ysim = DV) %>%
    binning(bin = IVAR) %>%
    vpcstats()
  plot(plot_binning_vpc)

Alternatively, one can create/customize VPC plots through VPC results shiny app (from Certara.VPCResults package), which can also be used to:

  • generate corresponding tidyvpc code to reproduce the VPC ouput from R command line
  • generate report as well as the associated R markdown.

Here we only demonstrate how to invoke this shiny app:

  vpcResultsUI(observed = dt_ObsData, simulated = dt_SimData)