categorical_data.Rmd
The purpose of this vignette is to demonstrate how to:
tidyvpc
(command-line usage) and VPC results Shiny app (in
Certara.VPCResults
package)We assume that all the necessary packages are loaded and the directory with NLME Executables is given as an environment variable (INSTALLDIR).
# 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 data
OneCpt1stOrderAbsorp_Emax_CategoricalModel.csv
and the PML
model file OneCpt1stOrderAbsorp_Emax_CategoricalModel.mdl
distributed with the Certara.RsNLME
package. First, we will
import the data
filename <- system.file("vignettesdata/OneCpt1stOrderAbsorp_Emax_CategoricalModel.csv",
package = "Certara.RsNLME",
mustWork = TRUE)
dt_InputDataSet <- fread(filename)
Next we will locate the PML model file and then create a textual model object for it.
filename <- system.file("vignettesdata/OneCpt1stOrderAbsorp_Emax_CategoricalModel.mdl",
package = "Certara.RsNLME",
mustWork = TRUE)
# Load the PML codes and link it to associated input data to create a model object
model <- textualmodel(modelName = "OneCpt1stOrderAbsorp_Emax_CategoricalModel",
mdl = filename,
data = dt_InputDataSet)
Let’s view the model and its associated column mappings, and then map those un-mapped model variables to their corresponding input data columns:
# View the model and its associated column mappings
print(model)
Model Overview -------------------------------------------
: OneCpt1stOrderAbsorp_Emax_CategoricalModel
Model Name : C:/Users/jcraig/Documents/GitHub/R-RsNLME/vignettes/OneCpt1stOrderAbsorp_Emax_CategoricalModel
Working Directory : Textual
Model Type
PML -------------------------------------------
test(){
# ===============================================================
# PK model: one compartment model with 1st order absorption
# ===============================================================
cfMicro(A1, Cl / V, first = (Aa = Ka))
dosepoint(Aa)
= A1 / V
C
# residual error model
error(CEps = 0.1)
observe(CObs = C * (1 + CEps))
# ----------------------------------------------------------------
# PK model parameters
# ----------------------------------------------------------------
# Structural model parameters
stparm(Ka = exp(tvlogKa + nlogKa))
stparm(V = exp(tvlogV + nlogV))
stparm(Cl = exp(tvlogCl + nlogCl))
# fixed effects
fixef(tvlogKa = c(, -1, ))
fixef(tvlogV = c(, 2, ))
fixef(tvlogCl = c(, 0, ))
# random effects
ranef(diag(nlogV, nlogCl, nlogKa) = c(1, 1, 1))
# ================================================================
# PD model
# ================================================================
= Emax * C / (EC50 + C)
E
## Residual error model
error(EEps = 0.1)
observe(EObs = E * (1 + EEps))
## Categorical model
multi(CategoricalObs, ilogit, -E, -(E + CatParam))
# ----------------------------------------------------------------
# Categorical model parameters
# ----------------------------------------------------------------
# structural model parameters
stparm(EC50 = exp(tvlogEC50 + nlogEC50))
stparm(Emax = exp(tvlogEmax + nlogEmax))
stparm(CatParam = exp(tvlogCatParam + nlogCatParam))
# fixed effects
fixef(tvlogEC50 = c(, 2, ))
fixef(tvlogEmax = c(, -2, ))
fixef(tvlogCatParam = c(, 1, ))
# random effects
ranef(diag(nlogEC50, nlogEmax, nlogCatParam) = c(1, 1, 1))
}
Structural Parameters -------------------------------------------
Ka V Cl EC50 Emax CatParam-------------------------------------------
Column Mappings -------------------------------------------
: Data Column name
Model Variable Name : ?
id : time
time : ?
Aa : CObs
CObs : EObs
EObs : CategoricalObs
CategoricalObs
# Manually map those un-mapped model variables to their corresponding input data columns
<- model %>%
model colMapping(c(id = "SubID", Aa = "dose_Aa"))
Next, we will run the model using the fitmodel
function with default host. We will use the QRPEM method for fitting. We
will also output residuals PCWRES with the number of replicates set to
be 1000 (Note: PCWRES is not outputted by default).
<- fitmodel(model, method = "QRPEM", numRepPCWRES = 1000)
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 2711.294 -5422.588 -5394.588 -5297.671 14 7500 300
EpsShrinkage Condition1: 0.09678 3.46327
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 = "OneCpt1stOrderAbsorp_Emax_CategoricalModel")
## Filter out CategoricalObs
xp <- xp %>%
filter(ObsName != "CategoricalObs")
## observations against population predictions
dv_vs_pred(xp, type = "p", facets = "ObsName", subtitle = "-2LL: @ofv")
## observations against individual predictions
dv_vs_ipred(xp, type = "p", facets = "ObsName", subtitle = "-2LL: @ofv, Eps shrinkage: @epsshk")
## CWRES against population predictions
res_vs_pred(xp, res = "CWRES", type = "ps", facets = "ObsName", 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", facets = "ObsName", subtitle = "-2LL: @ofv")
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
## PCWRES against population predictions
res_vs_pred(xp, res = "PCWRES", type = "ps", facets = "ObsName", subtitle = "-2LL: @ofv")
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
## PCWRES against the independent variable
res_vs_idv(xp, res = "PCWRES", type = "ps", facets = "ObsName", subtitle = "-2LL: @ofv")
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
Alternatively, one can view/customize diagnostic plots as well as
estimation results using the Certara.ModelResults
Shiny
application, which can also be used to generate .R
and/or
.Rmd
code based on operations performed in the GUI. For
installation and usage details, please visit the following link.
Here we only demonstrate how to invoke this Shiny app through either the
NlmePmlModel
object or the xpose_data
object
created above.
We will use the copyModel
function to copy the model into a new object and accept final parameter
estimates from fitting run as initial estimates for VPC simulation:
<- copyModel(model,
modelVPC acceptAllEffects = TRUE,
modelName = "OneCpt1stOrderAbsorp_Emax_CategoricalModel_VPC")
## View model
print(modelVPC)
Model Overview -------------------------------------------
: OneCpt1stOrderAbsorp_Emax_CategoricalModel_VPC
Model Name : C:/Users/jcraig/Documents/GitHub/R-RsNLME/vignettes/OneCpt1stOrderAbsorp_Emax_CategoricalModel_VPC
Working Directory : Textual
Model Type
PML -------------------------------------------
test(){
# ===============================================================
# PK model: one compartment model with 1st order absorption
# ===============================================================
cfMicro(A1, Cl / V, first = (Aa = Ka))
dosepoint(Aa)
= A1 / V
C
# residual error model
error(CEps = 0.121544989414485)
observe(CObs = C * (1 + CEps))
# ----------------------------------------------------------------
# PK model parameters
# ----------------------------------------------------------------
# Structural model parameters
stparm(Ka = exp(tvlogKa + nlogKa))
stparm(V = exp(tvlogV + nlogV))
stparm(Cl = exp(tvlogCl + nlogCl))
# fixed effects
fixef(tvlogKa = c(,-0.356862376082258,))
fixef(tvlogV = c(,1.63701201337159,))
fixef(tvlogCl = c(,-0.224407387640359,))
# random effects
ranef(diag(nlogV, nlogCl, nlogKa) = c(0.086623473, 0.1052053, 0.11503702))
# ================================================================
# PD model
# ================================================================
= Emax * C / (EC50 + C)
E
## Residual error model
error(EEps = 0.179381190657918)
observe(EObs = E * (1 + EEps))
## Categorical model
multi(CategoricalObs, ilogit, -E, -(E + CatParam))
# ----------------------------------------------------------------
# Categorical model parameters
# ----------------------------------------------------------------
# structural model parameters
stparm(EC50 = exp(tvlogEC50 + nlogEC50))
stparm(Emax = exp(tvlogEmax + nlogEmax))
stparm(CatParam = exp(tvlogCatParam + nlogCatParam))
# fixed effects
fixef(tvlogEC50 = c(,2.2899579614783,))
fixef(tvlogEmax = c(,-2.30045875119669,))
fixef(tvlogCatParam = c(,0.418715524594699,))
# random effects
ranef(diag(nlogEC50, nlogEmax, nlogCatParam) = c(0.066093482, 0.10266479, 0.15212821))
}
Structural Parameters -------------------------------------------
Ka V Cl EC50 Emax CatParam-------------------------------------------
Column Mappings -------------------------------------------
: Data Column name
Model Variable Name : SubID
id : time
time : dose_Aa
Aa : CObs
CObs : EObs
EObs : CategoricalObs CategoricalObs
Now, let’s run VPC using the vpcmodel
function with the default host, default values for the relevant NLME
engine arguments, and PRED outputted.
Note: Here we define VPC argument, outputPRED
,
through ellipsis (additional argument). Alternatively, one can define
VPC arguments through vpcParams
argument.
<- vpcmodel(modelVPC, outputPRED = TRUE)
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/OneCpt1stOrderAbsorp_Emax_CategoricalModel_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
## predcheck0 contains observed data for all continuous observed variables
<- vpcJob$predcheck0
dt_ObsData_ContinuousObs
## predcheck0_cat contains observed data for Categorical/count observed variables
<- vpcJob$predcheck0_cat
dt_ObsData_CategoricalObs
## predout contains simulated data for all observed variables
<- 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 contains 3 observed variable with PRED outputted. Hence, to
use this package, we have to do some data preprocessing on both
simulated and observed data to meet the requirements set by the
tidyvpc
package.
First we will process the simulated data output to pass to the
simulated()
function in the tidyvpc
package,
creating a separate data.frame
for each of our
DV
:
## Extract simulated data for observed variable "CObs"
dt_SimData_tidyvpc_CObs <- dt_SimData[OBSNAME == "CObs"]
## Extract simulated data for observed variable "EObs"
dt_SimData_tidyvpc_EObs <- dt_SimData[OBSNAME == "EObs"]
## Extract simulated data for observed variable "CategoricalObs"
dt_SimData_tidyvpc_CategoricalObs <- dt_SimData[OBSNAME == "CategoricalObs"]
Next, we will process the observed data output to pass to the
observed()
function in the tidyvpc
package,
creating a separate data.frame
for each of our
DV
as we did for the simulated data output:
## Extract observed data for observed variable "CObs"
dt_ObsData_ContinuousObs_tidyvpc_CObs <- dt_ObsData_ContinuousObs[ObsName == "CObs"]
## Extract observed data for observed variable "EObs"
dt_ObsData_ContinuousObs_tidyvpc_EObs <- dt_ObsData_ContinuousObs[ObsName == "EObs"]
## Extract observed data for observed variable "CategoricalObs"
dt_ObsData_CategoricalObs_tidyvpc <- dt_ObsData_CategoricalObs[ObsName == "CategoricalObs"]
Finally, we will add the $PRED
column from
REPLICATE == 0
($PRED
may be extracted from
any REPLICATE
) in the simulated data to our observed data
in order to perform a prediction-corrected VPC:
## Add PRED from REPLICATE == 0 of simulated data (CObs) to observed data (CObs)
dt_ObsData_ContinuousObs_tidyvpc_CObs$PRED <- as.numeric(dt_SimData_tidyvpc_CObs[REPLICATE == 0]$PRED)
## Add PRED from REPLICATE == 0 of simulated data (EObs) to observed data (EObs)
dt_ObsData_ContinuousObs_tidyvpc_EObs$PRED <- as.numeric(dt_SimData_tidyvpc_EObs[REPLICATE == 0]$PRED)
Now we can create VPC plots.
### Create a binless VPC plot for CObs
binless_VPC_CObs <- observed(dt_ObsData_ContinuousObs_tidyvpc_CObs, x = IVAR, yobs = DV) %>%
simulated(dt_SimData_tidyvpc_CObs, ysim = DV) %>%
binless() %>%
vpcstats()
plot(binless_VPC_CObs)
### Create a binless pcVPC plot for CObs
binless_pcVPC_CObs <- observed(dt_ObsData_ContinuousObs_tidyvpc_CObs,
x = IVAR,
yobs = DV) %>%
simulated(dt_SimData_tidyvpc_CObs, ysim = DV) %>%
binless(loess.ypc = TRUE) %>%
predcorrect(pred = PRED) %>%
vpcstats()
plot(binless_pcVPC_CObs)
### Create a binless VPC plot for EObs
binless_VPC_EObs <- observed(dt_ObsData_ContinuousObs_tidyvpc_EObs, x = IVAR, yobs = DV) %>%
simulated(dt_SimData_tidyvpc_EObs, ysim = DV) %>%
binless() %>%
vpcstats()
plot(binless_VPC_EObs)
### Create a binless pcVPC plot for EObs
binless_pcVPC_EObs <- observed(dt_ObsData_ContinuousObs_tidyvpc_EObs, x = IVAR, yobs = DV) %>%
simulated(dt_SimData_tidyvpc_EObs, ysim = DV) %>%
binless(loess.ypc = TRUE) %>%
predcorrect(pred = PRED) %>%
vpcstats()
plot(binless_pcVPC_EObs)
### Create a binless VPC plot for CatgoricalObs
binless_VPC_CategoricalObs <- observed(dt_ObsData_CategoricalObs_tidyvpc,
x = IVAR,
yobs = DV) %>%
simulated(dt_SimData_tidyvpc_CategoricalObs, ysim = DV) %>%
binless() %>%
vpcstats(vpc.type = "categorical")
plot(binless_VPC_CategoricalObs
, facet = TRUE
, facet.scales = "fixed"
, legend.position = "bottom"
)
### Create a binning VPC plot for CObs: binning on x-variable itself
binning_VPC_CObs <- observed(dt_ObsData_ContinuousObs_tidyvpc_CObs, x = IVAR, yobs = DV) %>%
simulated(dt_SimData_tidyvpc_CObs, ysim = DV) %>%
binning(bin = IVAR) %>%
vpcstats()
plot(binning_VPC_CObs)
### Create a binning pcVPC plot for CObs: binning on x-variable itself
binning_pcVPC_CObs <- observed(dt_ObsData_ContinuousObs_tidyvpc_CObs,
x = IVAR,
yobs = DV) %>%
simulated(dt_SimData_tidyvpc_CObs, ysim = DV) %>%
binning(bin = IVAR) %>%
predcorrect(pred = PRED) %>%
vpcstats()
plot(binning_pcVPC_CObs)
### Create a binning VPC plot for EObs: binning on x-variable itself
binning_VPC_EObs <- observed(dt_ObsData_ContinuousObs_tidyvpc_EObs, x = IVAR, yobs = DV) %>%
simulated(dt_SimData_tidyvpc_EObs, ysim = DV) %>%
binning(bin = IVAR) %>%
vpcstats()
plot(binning_VPC_EObs)
### Create a binning pcVPC plot for EObs: binning on x-variable itself
binning_pcVPC_EObs <- observed(dt_ObsData_ContinuousObs_tidyvpc_EObs, x = IVAR, yobs = DV) %>%
simulated(dt_SimData_tidyvpc_EObs, ysim = DV) %>%
binning(bin = IVAR) %>%
predcorrect(pred = PRED) %>%
vpcstats()
plot(binning_pcVPC_EObs)
### Create a binning VPC plot for CatgoricalObs: binning on x-variable itself
binning_VPC_CategoricalObs <- observed(dt_ObsData_CategoricalObs_tidyvpc,
x = IVAR,
yobs = DV) %>%
simulated(dt_SimData_tidyvpc_CategoricalObs, ysim = DV) %>%
binning(bin = IVAR) %>%
vpcstats(vpc.type = "categorical")
plot(binning_VPC_CategoricalObs
, facet = TRUE
, facet.scales = "fixed"
, legend.position = "bottom"
)
Alternatively, one can create/customize VPC plots through VPC results
shiny app (in Certara.VPCResults
package), which can also
be used to:
Here we only demonstrate how to invoke this shiny app (Note: The shiny app will automatically preprocess the data as what we did above for tidyvpc package).
library(Certara.VPCResults)
## Invoke VPC results shiny app to create VPC plots for CObs
taggedVPC_CObs <- vpcResultsUI(observed = dt_ObsData_ContinuousObs_tidyvpc_CObs,
simulated = dt_SimData_tidyvpc_CObs)
## Invoke VPC results shiny app to create VPC plots for EObs
taggedVPC_EObs <- vpcResultsUI(observed = dt_ObsData_ContinuousObs_tidyvpc_EObs,
simulated = dt_SimData_tidyvpc_EObs)
## Invoke VPC results shiny app to create VPC plot for categoricalObs
taggedVPC_CategoricalObs <- vpcResultsUI(observed = dt_ObsData_CategoricalObs_tidyvpc,
simulated = dt_SimData_tidyvpc_CategoricalObs,
vpc.type = "categorical")