Basic Goodness-of-fit Diagnostic Plots

Certara PMxO R Tools Workstream

17-Nov-2020

Setup

Before starting, load the ggplot2 and ggcertara packages, and set the default ggplot theme to theme_certara():

suppressPackageStartupMessages({
  library(ggplot2)
  library(ggcertara)
  library(kableExtra) #used for the vignette not for ggcertara demonstration
})
theme_set(theme_certara())

Overview

The ggcertara package and theme_certara theme are meant to provide a standardized look to typical gof plots employed by pharmacometricians. The syntax is based on ggplot2. The ggcertara also utilizes the patchwork package and syntax. These will be touched upon later in this vignette.

Data:

Data Expectations

The gof() function uses standard column names in lower case. The following column names are recognized by the functions in this package:

Not all of these need to be present, only those that are needed for the disired plots (for instance, if no plots involving iwres are requested, then it does not need to be present). The column names must be in lower case.

Reading:

The gof_read_data() function simplifies the data loading process by:

  1. loading the data
  2. converting all names to lower case
  3. checking for the presence of an mdv column and, if present, filtering the rows to keep only those with mdv = 0

For this vignette we will be using some sample data that is provided with the package. The data is contained in a file named nmtable.csv, that was actually generated by NONMEM using a $TABLE statement (i.e., $TABLE FILE=nmtable.csv FORMAT=,1PE15.8 ...). It’s contents look like this:

ID TIME DV PRED IPRED RES WRES CWRES CIWRES IRES IWRES MDV EVID TAD CL V KA ETA1 ETA2 ETA3
1 0.0000000 0.00000 0.00000 0.00000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1 1 0.0000000 3.311365 10.26801 2.860697 -0.144312 -0.1010664 0.0747987
1 0.1107401 16.84466 14.39002 16.98482 2.4546404 0.0413047 0.0306011 -0.1043802 -0.1401575 -0.1038259 0 0 0.1107401 3.311365 10.26801 2.860697 -0.144312 -0.1010664 0.0747987
1 0.2493163 33.22216 26.62746 31.08233 6.5946926 1.2657456 1.1250718 0.8657021 2.1398227 0.8661927 0 0 0.2493163 3.311365 10.26801 2.860697 -0.144312 -0.1010664 0.0747987
1 0.4983335 43.48446 38.21232 43.90453 5.2721396 0.2099566 0.1983909 -0.1206951 -0.4200702 -0.1203824 0 0 0.4983335 3.311365 10.26801 2.860697 -0.144312 -0.1010664 0.0747987
1 0.8483932 48.27355 42.64485 48.29827 5.6286966 0.3379642 0.3135719 -0.0066225 -0.0247160 -0.0064387 0 0 0.8483932 3.311365 10.26801 2.860697 -0.144312 -0.1010664 0.0747987
1 2.1580771 32.07656 31.68700 35.66761 0.3895546 -1.0379582 -0.9256054 -1.2668550 -3.5910498 -1.2667701 0 0 2.1580771 3.311365 10.26801 2.860697 -0.144312 -0.1010664 0.0747987

In a normal project workflow, the location of the appropriate file would be clear. For the vignette, the file can be located as follows:

Now we load the file:

Examining the data, we can see that the names are lower case and the rows have been filtered such that mdv = 0:

#>   id      time       dv     pred    ipred       res        wres      cwres       ciwres        ires        iwres mdv evid       tad
#> 2  1 0.1107401 16.84466 14.39002 16.98482 2.4546404  0.04130471  0.0306011 -0.104380162 -0.14015753 -0.103825922   0    0 0.1107401
#> 3  1 0.2493163 33.22216 26.62746 31.08233 6.5946926  1.26574562  1.1250718  0.865702143  2.13982267  0.866192653   0    0 0.2493163
#> 4  1 0.4983335 43.48446 38.21232 43.90453 5.2721396  0.20995661  0.1983909 -0.120695079 -0.42007023 -0.120382378   0    0 0.4983335
#> 5  1 0.8483932 48.27355 42.64485 48.29827 5.6286966  0.33796419  0.3135719 -0.006622542 -0.02471603 -0.006438689   0    0 0.8483932
#> 6  1 2.1580771 32.07656 31.68700 35.66761 0.3895546 -1.03795820 -0.9256054 -1.266855000 -3.59104976 -1.266770110   0    0 2.1580771
#> 7  1 3.9157063 22.78889 17.64876 20.31905 5.1401243  2.03836621  1.8503260  1.529135290  2.46983532  1.529378070   0    0 3.9157063
#>         cl        v       ka      eta1       eta2       eta3
#> 2 3.311365 10.26801 2.860697 -0.144312 -0.1010664 0.07479871
#> 3 3.311365 10.26801 2.860697 -0.144312 -0.1010664 0.07479871
#> 4 3.311365 10.26801 2.860697 -0.144312 -0.1010664 0.07479871
#> 5 3.311365 10.26801 2.860697 -0.144312 -0.1010664 0.07479871
#> 6 3.311365 10.26801 2.860697 -0.144312 -0.1010664 0.07479871
#> 7 3.311365 10.26801 2.860697 -0.144312 -0.1010664 0.07479871
#> 
#>   0 
#> 831

Basic usage

The basic functionality is accessed through the gof() function:

gof(dat) #Note: output aspect ratio is set to 3:2, in this case width=9, height=6

By default, the gof() function produces a figure with 6 panels arranged in a 2-by-3 grid.

Column_1 Column_2 Column_3
DV vs. IPRED CWRES vs. TIME CWRES vs. TAD
DV vs. PRED CWRES vs. PRED |IWRES| vs. IPRED

Note that all axes are nicely labelled (though they lack units at this point). With this simple function call, we have produced a report ready figure!

Customizing the output

Selecting panels

The gof() function can easily produce a variety of panels which can be referenced as either numbered panels, named functions, or as strings.

Number Function Name Description
1 gof_cwres_histogram() cwres_histogram Histogram of CWRES
2 gof_cwres_qqplot() cwres_qqplot QQ-plot of CWRES
3 gof_dv_vs_ipred() dv_vs_ipred DV vs. IPRED
4 gof_dv_vs_pred() dv_vs_pred DV vs. PRED
5 gof_cwres_vs_pred() cwres_vs_pred CWRES vs. PRED
6 gof_cwres_vs_time() cwres_vs_time CWRES vs. TIME
7 gof_cwres_vs_tad() cwres_vs_tad CWRES vs. TAD
8 gof_absiwres_vs_ipred() absiwres_vs_ipred |IWRES| vs. IPRED
9 gof_absiwres_vs_time() absiwres_vs_time |IWRES| vs. TIME
10 gof_absiwres_vs_tad() absiwres_vs_tad |IWRES| vs. TAD
11, -3 gof_dv_vs_ipred(log_xy=T) dv_vs_ipred_log DV vs. IPRED log scale
12, -4 gof_dv_vs_pred(log_xy=T) dv_vs_pred_log DV vs. PRED log scale
13, -5 gof_cwres_vs_pred(log_x=T) cwres_vs_pred_log CWRES vs. PRED log scale
14, -6 gof_cwres_vs_time(log_x=T) cwres_vs_time_log CWRES vs. TIME log scale
15, -7 gof_cwres_vs_tad(log_x=T) cwres_vs_tad_log CWRES vs. TAD log scale
16, -8 gof_absiwres_vs_ipred(log_x=T) absiwres_vs_ipred_log |IWRES| vs. IPRED log scale
17, -9 gof_absiwres_vs_time(log_x=T) absiwres_vs_time_log |IWRES| vs. TIME log scale
18, -10 gof_absiwres_vs_tad(log_x=T) absiwres_vs_tad_log |IWRES| vs. TAD log scale

Note that the standard 6 plots in the gof() function correspond to the numbered panels 3, 6, 7, 4, 5, 8 from this list.

Calling the panels option within gof()

Using the panels argument, select panels can be selected using either the respective panel numbers or names (as strings). For example, to prduce both a histogram and QQ-plot of CWRES, select panels 1 and 2 or callthem by name cwres_histogram and cwres_qqplot:

For log-scale plots (panels 11-18 above), we can call their panel number or add a - to their linear scale counterpart. For example, to obtain a plot of DV vs. IPRED on log-log scale, we can set panel=11 or panel=-3:

We can use the layout option to specify the number of rows and columns in a multi-panel output.

Named function:

We can also call specific panels using the corresponding named function. Log variants of the plot can be obtained by specifying either log_x=T or log_xy=T depending on the function.

The gof_list object

Individual panels can also be with the gof_list() function. This returns a numbered list of the 6 standard panels by default, or all panels if all=T is specified:

#> List of gof plots:
#> 
#>   1. cwres_histogram
#>   2. cwres_qqplot
#>   3. dv_vs_ipred
#>   4. dv_vs_pred
#>   5. cwres_vs_pred
#>   6. cwres_vs_time
#>   7. cwres_vs_tad
#>   8. absiwres_vs_ipred
#>   9. absiwres_vs_time
#>   10. absiwres_vs_tad
#>   11. dv_vs_ipred_log
#>   12. dv_vs_pred_log
#>   13. cwres_vs_pred_log
#>   14. cwres_vs_time_log
#>   15. cwres_vs_tad_log
#>   16. absiwres_vs_ipred_log
#>   17. absiwres_vs_time_log
#>   18. absiwres_vs_tad_log

Individual panels can be extracted from the list:

gof_list can also be used to build a new list of plots, by specifying empty=T:

A composite plot can be obtained from a list using the gof_layout() function:

Labels and units

There are several ways to customize the plot axis labels/units. When a recognized variable name is used in a plot (e.g. ipred, cwres, etc.), the functions will find a corresponding axis label in a list object which by default corresponds to a global option:

Option Default
gof.label.dv Observed Value
gof.label.pred Population Prediction
gof.label.ipred Individual Prediction
gof.label.cwres Conditional Weighted Residual
gof.label.iwres Individual Weighted Residual
gof.label.time Time
gof.label.tad Time After Dose
gof.units.dv NA
gof.units.time NA
gof.units.tad NA

Every function has a labels argument to override the default, should it be necessary. Also note when updating labels/units, it is not required to specify all labels/units, only the labels/units that are being updated need be specified.

Global label/units update

We can update labels and units at the global level by calling the desired option within the options function. Units for concentration (dv, pred, ipred), time, and tad can also be set using global options and will be appended to the label:

Update labels for a single plot

We can update a default label for a single plot by calling the label option as follows.