distr_delay.Rmd
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
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" [
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:
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 -------------------------------------------
: OneCpt_GammaDistributedAbsorptionDelay
Model Name : C:/Users/jcraig/Documents/GitHub/R-RsNLME/vignettes/OneCpt_GammaDistributedAbsorptionDelay
Working Directory : TRUE
Is population : PK
Model Type
PK -------------------------------------------
: Clearance
Parameterization : Gamma
Absorption : 1
Num Compartments : FALSE
Dose Tlag? : FALSE
Elimination Comp ?: FALSE
Infusion Allowed ?: FALSE
Sequential : FALSE
Freeze PK
PML -------------------------------------------
test(){
dosepoint(A1)
= A1 / V
C 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: CObs
Observation Name : C
Effect Name : CEps
Epsilon Name : Multiplicative
Epsilon Type : FALSE
Epsilon frozen : FALSE
is BQL -------------------------------------------
Column Mappings -------------------------------------------
: Data Column name
Model Variable Name : ID
id : Time
time : ?
A1 : CObs
CObs
## Manually map those un-mapped model variables to their corresponding input data columns
<- model %>% colMapping(c(A1 = "Dose")) model
Let’s run the model using the fitmodel
function with default host and engine method set to be QRPEM.
<- fitmodel(model, method = "QRPEM")
job 4 cores
Using MPI host with in the host class instance.
sharedDirectory is not given :
Valid NLME_ROOT_DIRECTORY is not given, using current working directory:/Users/jcraig/Documents/GitHub/R-RsNLME/vignettes C
print(job$Overall)
-2LL AIC BIC nParm nObs nSub
Scenario RetCode LogLik 1: WorkFlow 1 685.0575 -1370.115 -1354.115 -1317.481 8 720 80
EpsShrinkage Condition1: 0.14701 2.65545
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")
eta_distrib(xp)
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.
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:
<- vpcmodel(modelVPC)
vpcJob
Using localhost without parallelization.
NLME Jobin the host class instance.
sharedDirectory is not given :
Valid NLME_ROOT_DIRECTORY is not given, using current working directory:/Users/jcraig/Documents/GitHub/R-RsNLME/vignettes
C
1 of 1 NLME models
Compiling
The model compiled
Trying to generate job results...
Copying tables predcheck_bql.csv, predcheck0.csv, predcheck0_cat.csv, predcheck1.csv, predcheck1_cat.csv, predcheck2.csv, predcheck2_cat.csv, predout.csv
Finished summarizing results. Transferring data and loading the results...
Done generating job results.
/Simulation results are ready in C:/Users/jcraig/Documents/GitHub/R-RsNLME/vignettes/OneCpt_GammaDistributedAbsorptionDelay_VPC
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
<- vpcJob$predcheck0
dt_ObsData
## Simulated data
<- vpcJob$predout dt_SimData
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:
Here we only demonstrate how to invoke this shiny app:
vpcResultsUI(observed = dt_ObsData, simulated = dt_SimData)