exploring_model_diagnostics.Rmd
The purpose of this vignette is to demonstrate how to:
xpose_data
object from NLME model execution
directory/filesxpose
package to generate
basic GOF plotsxpose
plot functions with
ggplot2
Certara.Xpose.NLME
xpose_data
object and explore
usage of ggcertara
Note: The corresponding R script used in this example can be
found in
RsNLME_Examples/TwoCptIVBolus_FitBaseModel_CovariateSearch_VPC_BootStrapping.R
.
See RsNLME
Example Scripts.
We will be using the built-in pkData
from the
Certara.RsNLME
package for the following examples.
pkData
is a pharmacokinetic dataset containing 16 subjects
with single bolus dose.
library(Certara.RsNLME)
library(dplyr)
library(magrittr)
pkData <- Certara.RsNLME::pkData
glimpse(pkData)
#> Rows: 112
#> Columns: 8
#> $ Subject <dbl> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,…
#> $ Nom_Time <dbl> 0, 1, 1, 2, 4, 8, 24, 0, 1, 1, 2, 4, 8, 24, 0, 1, 1, 2, 4, …
#> $ Act_Time <dbl> 0.00, 0.26, 1.10, 2.10, 4.13, 8.17, 23.78, 0.00, 0.37, 0.83…
#> $ Amount <dbl> 25000, NA, NA, NA, NA, NA, NA, 25000, NA, NA, NA, NA, NA, N…
#> $ Conc <dbl> 2010.0, 1330.0, 565.0, 216.0, 180.0, 120.0, 16.9, 850.0, 79…
#> $ Age <dbl> 22, 22, 22, 22, 22, 22, 22, 46, 46, 46, 46, 46, 46, 46, 52,…
#> $ BodyWeight <dbl> 73, 73, 73, 73, 73, 73, 73, 79, 79, 79, 79, 79, 79, 79, 63,…
#> $ Gender <chr> "male", "male", "male", "male", "male", "male", "male", "fe…
Before exploring usage of Certara.Xpose.NLME
and
evaulating our covariate model diagnostic functions, we will first:
For a detailed explanation on model building and execution in
Certara.RsNLME
visit here.
# Create model
baseModel <- pkmodel(
numCompartments = 2, data = pkData,
ID = "Subject", Time = "Act_Time", A1 = "Amount", CObs = "Conc",
modelName = "TwCpt_IVBolus_FOCE_ELS"
)
# Update parameters
baseModel <- baseModel %>%
structuralParameter(paramName = "V2", hasRandomEffect = FALSE) %>%
fixedEffect(effect = c("tvV", "tvCl", "tvV2", "tvCl2"), value = c(15, 5, 40, 15)) %>%
randomEffect(effect = c("nV", "nCl", "nCl2"), value = rep(0.1, 3)) %>%
residualError(predName = "C", SD = 0.2)
# Fit model
baseFitJob <- fitmodel(baseModel)
#> MPI not found on the system.
#> Using localhost without parallelization.
#> 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-Xpose-NLME/vignettes
covariateModel <- copyModel(baseModel,
acceptAllEffects = TRUE,
modelName = "TwCpt_IVBolus_SelectedCovariateModel_FOCE_ELS"
)
covariateModel <- covariateModel %>%
addCovariate(covariate = c(Sex = "Gender"), effect = c("V", "Cl"), type = "Categorical", levels = c(0, 1), labels = c("female", "male")) %>%
addCovariate(covariate = "Age", effect = c("V", "Cl")) %>%
addCovariate(covariate = c(BW = "BodyWeight"), effect = c("V", "Cl"), center = "Value", centerValue = 70)
covariateJob <- fitmodel(covariateModel)
#> MPI not found on the system.
#> Using localhost without parallelization.
#> 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-Xpose-NLME/vignettes
xpose_data
object from RsNLME model output
We can create our xpose_data
object from our RsNLME
output by supplying the output directory where our model execution files
have been saved. If the model in your R environment has been executed
previously, you may pass the folder path from
model@modelInfo@workingDir
directly to the dir
argument of the xposeNlme()
function.
library(Certara.Xpose.NLME)
library(xpose)
xpdb <- xposeNlme(dir = covariateModel@modelInfo@workingDir)
class(xpdb)
#> [1] "xpose_data" "uneval"
We will use the list_vars()
function from
xpose
to print out useful information about our
xpose_data
object.
list_vars(xpdb)
#>
#> List of available variables for problem no. 1
#> - Subject identifier (id) : ID
#> - Dependent variable (dv) : DV
#> - Independent variable (idv) : IVAR
#> - DV identifier (dvid) : ObsName
#> - Dose amount (amt) : Amount
#> - Event identifier (evid) : EVID
#> - Model typical predictions (pred) : PRED
#> - Model individual predictions (ipred) : IPRED
#> - Model parameter (param) : V, Cl, V2, Cl2
#> - Eta (eta) : nV, nCl, nCl2
#> - Residuals (res) : IRES, IWRES, WRES, CWRES, PCWRES
#> - Categorical covariates (catcov) : Sex
#> - Continuous covariates (contcov) : Age, BW
#> - Not attributed (na) : Subject, Nom_Time, Act_Time, Conc, BodyWeight, Gender, WhichReset, TIME, TAD, PREDSE, Weight, CdfPCWRES, CdfDV, WhichDose
Covariate model diagnostics plot types vary given the covariate type (continuous or categorical). Continuous covariates are plotted in a scatter plot, while categorical covariates in a box plot.
The following “xpose-like” covariate model diagnostics functions are
available in Certara.Xpose.Nlme
:
The above functions follow xpose
conventions, including
usage of the type
argument to specify the type of plot
attributes and usage of the ...
ellipses to allow pass
through geom arguments from ggplot2
.
Note: The returned plot from xpose
and
Certara.Xpose.Nlme
is of class ggplot
,
therefore functions from ggplot2
can be applied to the plot
for further customized styling.
Let’s look inside the data and preview which residuals are available for plotting.
library(dplyr)
library(xpose)
gofData <- get_data(xpdb)
gofData %>%
select(ends_with("RES"))
#> # A tibble: 112 × 6
#> IRES IWRES WRES CWRES PCWRES CdfPCWRES
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 264. 0.951 1.20 1.15 0 0
#> 2 120. 0.625 0.674 0.652 0 0
#> 3 109. 1.51 1.07 1.07 0 0
#> 4 -24.3 -0.637 -1.03 -1.10 0 0
#> 5 19.5 0.767 0.458 0.519 0 0
#> 6 14.5 0.863 0.237 0.267 0 0
#> 7 -4.79 -1.39 -1.67 -1.93 0 0
#> 8 -248. -1.42 -1.59 -1.69 0 0
#> 9 33.5 0.277 0.397 0.402 0 0
#> 10 -140. -1.75 -1.32 -1.27 0 0
#> # … with 102 more rows
library(Certara.Xpose.NLME)
res_vs_cov(xpdb, covariate = "Age", res = "CWRES", type = "ps")
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
library(Certara.Xpose.NLME)
res_vs_cov(xpdb, covariate = "Sex", res = "CWRES", type = "b")
ggplot2
library(Certara.Xpose.NLME)
library(ggplot2)
res_vs_cov(xpdb, covariate = "BW", res = "CWRES", type = "ps") +
geom_hline(yintercept = 0, linetype = "dashed", lwd = 1, color = "red") +
theme_bw2()
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
eta_vs_cov(xpdb, covariate = "Age", type = "ps")
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
prm_vs_cov(xpdb, covariate = "BW", type = "ps")
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
prm_vs_cov(xpdb, covariate = "Sex", type = "b")
Note: We can remove the variable V2
from the
plots, which has a fixed value, by including the argument
drop_fixed = TRUE
.
ggcertara
ggcertara
is an R package to provide a standardized look
for plots employed by pharmacometricians. It provides a ggplot2 theme,
color palette, and a collection of plotting functions for basic
goodness-of-fit diagnostic plots.
Goodness-of-fit (GOF) data can easily be extracted from the
xpose_data
object in order to utilize plot functions from
ggcertara
.
For ggcertara
installation and usage details, visit the
following link.
To learn more about the standardized set of basic goodness-of-fit
diagnostic plots available in ggcertara
, explore the
following vignette.
library(ggcertara)
gofData <- xpose::get_data(xpdb)
colnames(gofData) <- tolower(colnames(gofData)) # ggcertara expects lower case column names
glimpse(gofData)
#> Rows: 112
#> Columns: 37
#> $ subject <int> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,…
#> $ nom_time <int> 0, 1, 1, 2, 4, 8, 24, 0, 1, 1, 2, 4, 8, 24, 0, 1, 1, 2, 4, …
#> $ act_time <dbl> 0.00, 0.26, 1.10, 2.10, 4.13, 8.17, 23.78, 0.00, 0.37, 0.83…
#> $ amount <int> 25000, NA, NA, NA, NA, NA, NA, 25000, NA, NA, NA, NA, NA, N…
#> $ conc <dbl> 2010.0, 1330.0, 565.0, 216.0, 180.0, 120.0, 16.9, 850.0, 79…
#> $ bodyweight <int> 73, 73, 73, 73, 73, 73, 73, 79, 79, 79, 79, 79, 79, 79, 63,…
#> $ gender <chr> "male", "male", "male", "male", "male", "male", "male", "fe…
#> $ id <fct> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,…
#> $ whichreset <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
#> $ time <dbl> 0.00, 0.26, 1.10, 2.10, 4.13, 8.17, 23.78, 0.00, 0.37, 0.83…
#> $ sex <fct> 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1,…
#> $ age <dbl> 22, 22, 22, 22, 22, 22, 22, 46, 46, 46, 46, 46, 46, 46, 52,…
#> $ bw <dbl> 73, 73, 73, 73, 73, 73, 73, 79, 79, 79, 79, 79, 79, 79, 63,…
#> $ v <dbl> 14.31602, 14.31602, 14.31602, 14.31602, 14.31602, 14.31602,…
#> $ cl <dbl> 7.612771, 7.612771, 7.612771, 7.612771, 7.612771, 7.612771,…
#> $ v2 <dbl> 41.44404, 41.44404, 41.44404, 41.44404, 41.44404, 41.44404,…
#> $ cl2 <dbl> 13.19875, 13.19875, 13.19875, 13.19875, 13.19875, 13.19875,…
#> $ ivar <dbl> 0.00, 0.26, 1.10, 2.10, 4.13, 8.17, 23.78, 0.00, 0.37, 0.83…
#> $ tad <dbl> 0.00, 0.26, 1.10, 2.10, 4.13, 8.17, 23.78, 0.00, 0.37, 0.83…
#> $ pred <dbl> 1659.31000, 1165.17000, 461.91900, 255.80400, 175.49500, 11…
#> $ ipred <dbl> 1746.3000, 1209.8100, 455.5610, 240.3010, 160.4620, 105.541…
#> $ dv <dbl> 2010.0, 1330.0, 565.0, 216.0, 180.0, 120.0, 16.9, 850.0, 79…
#> $ ires <dbl> 263.704000, 120.189000, 109.439000, -24.300800, 19.538000, …
#> $ predse <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
#> $ weight <dbl> 3.27917e-07, 6.83226e-07, 4.81845e-06, 1.73177e-05, 3.88379…
#> $ iwres <dbl> 0.9507100, 0.6254520, 1.5124300, -0.6366700, 0.7665790, 0.8…
#> $ wres <dbl> 1.20045000, 0.67415100, 1.07359000, -1.03417000, 0.45761200…
#> $ cwres <dbl> 1.14794000, 0.65225900, 1.06750000, -1.09698000, 0.51872100…
#> $ pcwres <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
#> $ cdfpcwres <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
#> $ cdfdv <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
#> $ obsname <fct> CObs, CObs, CObs, CObs, CObs, CObs, CObs, CObs, CObs, CObs,…
#> $ whichdose <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
#> $ nv <dbl> -0.05109783, -0.05109783, -0.05109783, -0.05109783, -0.0510…
#> $ ncl <dbl> 0.07114151, 0.07114151, 0.07114151, 0.07114151, 0.07114151,…
#> $ ncl2 <dbl> -0.06440154, -0.06440154, -0.06440154, -0.06440154, -0.0644…
#> $ evid <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
Use the gof()
function to plot basic GOF plot with
default panel options.
gof(gofData)
You may select from different combinations of gof plots using the
panels
argument.
gof(gofData, panels = 3:6)
The ggcertara
package contains a useful Shiny
application to explore different color palettes and test how those with
color blindness will perceive your selected colors.
First, create the plot and assign a return value.
myplot <- gof(gofData)
Next, pass the plot to the plotobj
argument in
run_colorexplorer()
.
run_colorexplorer(plotobj = myplot)
To simulate color blindness given the colors chosen in the palette, click the checkbox labeled “Simulate Color Blindness”.