Exploring Model Diagnostics

The purpose of this vignette is to demonstrate how to:

Create an xpose_data object from NLME model execution directory/files
Utilize functions from the xpose package to generate basic GOF plots
Extend xpose plot functions with ggplot2
Explore application of covariate model diagnostics functions, only available in Certara.Xpose.NLME
Extract GOF data from the 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.

Data

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…

Model Building and Execution

Before exploring usage of Certara.Xpose.NLME and evaulating our covariate model diagnostic functions, we will first:

Build and execute our base model
Copy the base model and accept parameter estimates
Add continuous and categorical covariates to the new model
Execute the new covariate model

For a detailed explanation on model building and execution in Certara.RsNLME visit here.

Build and execute base model

# 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
#> TDL5 version: 23.10.1.283
#> Servername: no-net
#> Sentinel License file: C:\PROGRA~1\Certara\NLME_E~1\lservrc
#> 
#> License Name: NU_PMLCommandLine_N
#> License Type: Commercial
#> Expiration Date: 31 Jan 2024
#> Current Date: 17 Nov 2023
#> Days until program expires: 75
#> The model compiled
#> 
#> Trying to generate job results...
#> Done generating job results.

Copy model and accept parameter estimates


covariateModel <- copyModel(baseModel,
  acceptAllEffects = TRUE,
  modelName = "TwCpt_IVBolus_SelectedCovariateModel_FOCE_ELS"
)

Add covariates to new model and fit

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
#> TDL5 version: 23.10.1.283
#> Servername: no-net
#> Sentinel License file: C:\PROGRA~1\Certara\NLME_E~1\lservrc
#> 
#> License Name: NU_PMLCommandLine_N
#> License Type: Commercial
#> Expiration Date: 31 Jan 2024
#> Current Date: 17 Nov 2023
#> Days until program expires: 75
#> The model compiled
#> 
#> Trying to generate job results...
#> Done generating job results.

Create `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"

List available variables

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

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:

res_vs_cov()
prm_vs_cov()
eta_vs_cov()

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.

Residuals vs Covariate

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  262.    0.943  1.18   1.13       0         0
#>  2  120.    0.622  0.670  0.649      0         0
#>  3  110.    1.52   1.08   1.08       0         0
#>  4  -24.0  -0.630 -1.02  -1.08       0         0
#>  5   19.7   0.775  0.488  0.546      0         0
#>  6   14.6   0.873  0.257  0.285      0         0
#>  7   -4.73 -1.38  -1.67  -1.92       0         0
#>  8 -250.   -1.43  -1.59  -1.69       0         0
#>  9   32.2   0.266  0.395  0.398      0         0
#> 10 -141.   -1.76  -1.32  -1.26       0         0
#> # ℹ 102 more rows

CWRES vs Continuous Covariate

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'

CWRES vs Categorical Covariate

library(Certara.Xpose.NLME)

res_vs_cov(xpdb, covariate = "Sex", res = "CWRES", type = "b")

Extending with `ggplot2`

Add reference line

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'

ETAs vs Covariates

ETAs vs Continuous Covariate


eta_vs_cov(xpdb, covariate = "Age", type = "ps")
#> `geom_smooth()` using formula = 'y ~ x'
#> `geom_smooth()` using formula = 'y ~ x'

ETAs vs Categorical Covariate


eta_vs_cov(xpdb, covariate = "Sex", type = "b")

Structural Parameters vs Covariates

Structural Parameters vs Continuous Covariate


prm_vs_cov(xpdb, covariate = "BW", type = "ps")
#> `geom_smooth()` using formula = 'y ~ x'
#> `geom_smooth()` using formula = 'y ~ x'

Structural Parameters vs Categorical Covariate


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.

Using `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.30108, 14.30108, 14.30108, 14.30108, 14.30108, 14.30108,…
#> $ cl         <dbl> 7.621016, 7.621016, 7.621016, 7.621016, 7.621016, 7.621016,…
#> $ v2         <dbl> 41.44576, 41.44576, 41.44576, 41.44576, 41.44576, 41.44576,…
#> $ cl2        <dbl> 13.19662, 13.19662, 13.19662, 13.19662, 13.19662, 13.19662,…
#> $ 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> 1666.63000, 1167.58000, 460.35100, 254.62900, 174.81800, 11…
#> $ ipred      <dbl> 1748.1200, 1210.4700, 455.2260, 240.0120, 160.2790, 105.389…
#> $ dv         <dbl> 2010.0, 1330.0, 565.0, 216.0, 180.0, 120.0, 16.9, 850.0, 79…
#> $ ires       <dbl> 261.881000, 119.526000, 109.774000, -24.011500, 19.721200, …
#> $ 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.27234e-07, 6.82479e-07, 4.82554e-06, 1.73594e-05, 3.89267…
#> $ iwres      <dbl> 0.943138, 0.621656, 1.518150, -0.629839, 0.774638, 0.872849…
#> $ wres       <dbl> 1.1756900, 0.6703270, 1.0830500, -1.0173700, 0.4876680, 0.2…
#> $ cwres      <dbl> 1.12728000, 0.64910300, 1.07547000, -1.07688000, 0.54619700…
#> $ 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.04773594, -0.04773594, -0.04773594, -0.04773594, -0.0477…
#> $ ncl        <dbl> 0.06744269, 0.06744269, 0.06744269, 0.06744269, 0.06744269,…
#> $ ncl2       <dbl> -0.06546828, -0.06546828, -0.06546828, -0.06546828, -0.0654…
#> $ evid       <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…

Basic GOF

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)

Color Explorer

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”.