Repository files navigation

Authors: Amine Gassem, Irina Baltcheva, Christian Bartels, SouvikBhattacharya, Thomas Dumortier, Matthew Fidler, Seid Hamzic, IngaLudwig, Ines Paule, Didier Renard, Bruno Bieth

ggPMX is an open-source R package freely available on CRAN since April2019. It generates standard diagnostic plots for mixed effect modelsused in pharmacometric activities. The package builds on the R-packageggplot2 and aims at providing a workflow that isconsistent,reproducible andefficient, resulting inhigh qualitygraphics ready-to-use in submission documents and publications.Intuitive functions and options allow for optimal figure customizationand graphics stratification. ggPMX enables straightforward generation ofPDF, Word or PNG output files that contain all diagnostic plots forkeeping track of modeling results. The package is currently compatiblewith Monolix versions 2016 and later, NONMEM version 7.2 and later andnlmixr.

Using simple syntax, the toolbox produces various goodness-of-fitdiagnostics such as: - residual- and empirical Bayes estimate(EBE)-based plots, - distribution plots, - prediction- andsimulation-based diagnostics (visual predictive checks).

ggPMX creates plot objects which can be further customized and extendedusing the extensiveggplot2syntax.

In addition, shrinkage and summary parameters tables can be alsoproduced. By default, the PDF- or Word-format diagnostic report containsessential goodness-of-fit plots. However, these can be adapted toproduce different sets of diagnostics as desired by the user, and any ofthe plots may be customized individually. The types of supportedcustomizations include modifications of the graphical parameters,labels, and various stratifications by covariates.

# Either install from CRANinstall.packages("ggPMX")# Or the development version from GitHub:# install.packages("devtools")devtools::install_github("ggPMXdevelopment/ggPMX")

Once ggPMX is installed, you can test if everything is working well. Ifsuccessful you should see a plot created after the following build-inexample.

library(ggPMX)ctr<- theophylline()ctr %>% pmx_plot_eta_matrix()

ggPMX is now ready for inputs and enhancements by the pharmacometriccommunity. - Please usepackageissues to fill inyour feedback.

TheggPMX package generates standard diagnostic plots and tables formixed effect models used in Pharmacometric (PMX) activities. The tool isbuilt upon the ggplot2 package and supports models developped withMonolix, NONMEM and nlmixr software. The current release (1.2) supportsmodels fitted with Monolix versions 2016 and later, NONMEM version 7.2and later and nlmixr.

The package aims to provide a workflow that is consistent, efficient andwhich results in high quality graphics ready to use in officialdocuments and reports. The package allows a high degree of flexibilityand customization, yet providing an acceptable default setting. Thepackage also allows to fully automate plots and report generation.

The general context is the analysis of mixed effect models fitted todata. ggPMX was developed in the framework of Pharmacometric activities,in which case the model is a population pharmacokinetic (PK) and/orpharmacodynamic (PD) model and the data is clinical or pre-clinical PKand/or PD data.

In the context of model building, evaluation and qualification, it isgood practice to assess the goodness-of-fit of models by inspecting(qualitatively and quantitatively) a set of graphs that indicate howwell the model describes the data. Several types of diagnostic plotsallow to evaluate a mixed effects model fit, the most common being:

1. residual-based plots
2. prediction-based plots
3. empirical Bayes estimates (EBE)-based plots
4. simulation-based plots.

The following figures are examples of diagnotic plots.

This document introduces the ggPMX functionalities and syntax.

The high level architecture is represented in the figure below. The keycomponents of the package are the following:

• Reader - reads model outputs from different sources (i.e. textfiles containing population parameters, model predictions,individual random effects, simulations and data-related inputs likecovariates) and restructures these outputs into standard formats forinternal processing.
• Generator ??? processes outputs from Reader. It contains Rlanguage code to produce the plots and is factorized into a smallset of flexible key functions. A set of default plots is defined ina configuration file. The configuration file can be adapted, e.g.,to have different configurations for different types of modelingactivities.
• Controller - serves as user interface. The user will callGenerator functions via wrapper functions in the Controller toproduce either all the default plots or selected plots of interest.In addition to editing the configuration, the user has differentoptions to adapt aspects of the plots to specific requirements.Plots may be adapted by setting parameters of the wrapper functionsthat generate the plots; there exist additional wrapper functions tochange aspects of the existing default plots. The plots are, ingeneral, returned as ggplot objects that can be further customizedusing ggplot functionalities.
• Reporter - generates sets of graphs and tables and integratesthem into an output file (Word or PDF) with annotations.

The package is coded using object-oriented programming meaning thatinformation is encoded in objects. The user can change or print thecontent of objects via functions. Such an implementation allows to havecode that is modular and easily customizable.

The typical workflow of ggPMX is as follows:

1. The user creates the Controller using pre-defined configurations(yaml templates) for plot settings.
2. The Controller implicitly calls the Reader that reads and storesmodelling outputs into a standard format. As a result, theController contains all available plots with their default valuesaccording to the configuration used.
3. The Generator allows to print the available plots by calling thecorresponding functions. Plots can be modified by using optionalarguments.
4. A call to the Reporter allows to create a pdf or docx report. Thereport Rmarkdown template can also be personalized.

The most important task for the user is the Controller creation. Thisstep requires careful consideration because it involves differentoptions according to the type of model (PK or PKPD) and software(Monolix, NONMEM or nlmixr) used for model fitting. The next sectiondescribes the Controller creation for the different possible cases.

Once the Controller is created, it implicitly calls the Reader andcreates the diagnostic plots. The user can then generate the graphs bycalling pre-defined functions. The same syntax is used independent ofthe model structure (PK or PKPD model) and of the fitting software.

The Reporter creates one report per endpoint, i.e., one report for PKand one for each PD endpoint.

For the sake of this document, three types of datasets are defined.

• Theinput modeling dataset is the one used for model fitting (theactual data). There are no particular requirements regarding thisdataset.
• Theoutput modeling datasets are those output from the fittingtool (Monolix, NONMEM or nlmixr). See the appendix for more detailson software requirements.
• TheggPMX datasets are the ones created within (internal to)ggPMX.

A diagnostic session starts with the creation of a Controller. TheController is the “user interface” of the package and allows to conrolall possible options. It is a container that stores configuration fields(model- and input data-related information), datasets and plots. It canbe used as a reference object in which the user can see the names of theexisitng plots, the names of theggPMX datasets, etc. The syntax ofthe Controller creation differs depending on the software used for modelfitting and on the number of model endpoints (or outputs). This sectionpresents different cases of Controller creation. For simplicity, thecase of models with one single output is presented first, thengeneralized to several outputs. ## Single-endpoint models

In general, models with only one endpoint (or output) are mostly PKmodels, but these could also be k-PD models.

To illustrateggPMX functionalities, the single-endpoint built-inmodel calledtheophylline is used hereafter. Thetheophyllinepopulation PK example has the following characteristics:

• The input modeling data contains PK samples of 2 studies, each with25 individuals who recieved a unique dose of 2000 mg theophylline.
• The model is a simple one-compartmental PK model with first-orderabsorption.
• The following covariates are used: weight (WT0) on volume (V) andclearance (Cl), age (AGE0), sex (SEX) and study (STUD) on clearance.
• Random effects on all three parameters (ka, Cl, V) are included.
• The residual error is proportional.

The input modeling dataset has the following columns:

#>   ID TIME  AMT   Y EVID WT0 AGE0 SEX STUD#> 1  1  0.0 2000   0    1  87   73   1    1#> 2  1  0.5    0 130    0  87   73   1    1#> 3  1  1.0    0 228    0  87   73   1    1#> 4  1  2.0    0 495    0  87   73   1    1#> 5  1  3.0    0 484    0  87   73   1    1#> 6  1  5.0    0 479    0  87   73   1    1

Note that the DVID (or CMT/YTYPE) column is missing, but since this is asingle-endpoint model, it is not necessary in that case.

The Contoller creation is wrapped in a function called “theophylline()”for quick reference:

ctr<- theophylline()

The controller initialization using the Monolix controllerpmx_mlx(),which is a wrapper function forpmx() withsys="mlx" (See AppendixA).

theophylline_path<- file.path(system.file(package="ggPMX"),"testdata","theophylline")work_dir<- file.path(theophylline_path,"Monolix")input_data_path<- file.path(theophylline_path,"data_pk.csv")ctr<- pmx_mlx(directory=work_dir,input=input_data_path,dv="Y")

####pmx_mlxtran()

The controller initialization can be simplified by using the Monolixcontrollerpmx_mlxtran(). This function parses the mlxtran file of aMonolix project and assigns automatically the different fields necessaryto the Controller creation. The only mandatory argument is file_name,the path to the mlxtran file.

mlxtran_path<- file.path(system.file(package="ggPMX"),"testdata","1_popPK_model","project.mlxtran")ctr<- pmx_mlxtran(file_name=mlxtran_path)

The user can verify the content of the Controller and how parameters areassigned by printing it.

The controller initialization using the NONMEM controllerpmx_nm() isbased on reading functions of the xpose package. It is highlyrecommended (but not required) to use the “sdtab, patab, cotab, catab”table naming convention followed by a run number(e.g. sdtab001,cotab001) This will enable automatic recognition ofcovariates. It is also recommended to name the model files accordingly(e.g. run001.lst). In order to generate a VPC a simulation dataset isrequired (see section about VPC)

For controller generation it is recommended to use the model file:

nonmem_dir<- file.path(system.file(package="ggPMX"),"testdata","extdata")ctr<- pmx_nm(directory=nonmem_dir,file="run001.lst" )

or the run number. The standard prefix is “run”, however can bespecified usingprefix

nonmem_dir<- file.path(system.file(package="ggPMX"),"testdata","extdata")ctr<- pmx_nm(directory=nonmem_dir,runno="001"#can be a string or a number)

It is simple to create a ggPMX controller for a nlmixr object usingpmx_nlmixr(). Using the theophylline example with a nlmixr model wehave:

one.cmt<-function() {  ini({## You may label each parameter with a commenttka<-0.45# Log Katcl<-1# Log Cl## This works with interactive models## You may also label the preceding line with label("label text")tv<-3.45; label("log V")## the label("Label name") works with all modelseta.ka~0.6eta.cl~0.3eta.v~0.1add.sd<-0.7  })  model({ka<- exp(tka+eta.ka)cl<- exp(tcl+eta.cl)v<- exp(tv+eta.v)    linCmt()~ add(add.sd)  })}fit<- nlmixr(one.cmt,theo_sd,est="saem",control=list(print=0))

Thefit object is a nlmixr fit; You can read it into the nlmixrcontroller by:

ctr<- pmx_nlmixr(fit,vpc=FALSE## VPC is turned on by default, can turn off)

The following are some optional arguments to the controller function(for details of each option, see the corresponding section or use?pmx_mlx,?pmx_nm,?pmx_nlmixr):

• cats: character vector of categorical covariates
• conts: character vector of continuous covariates
• occ: character occasional covariate variable name (currently notavailable for NONMEM or nlmixr)
• strats: character extra stratification variables
• settings: global pmxSettingsClass (pmx_settings()) shared betweenall plots
• endpoint: pmxEndpointClass (pmx_endpoint()) or integer orcharcater of the endpoint code (depends on the fitting software)
• sim: pmxSimClass object for VPC generation. (syntax for VPCgeneration is depending on the fitting software)

Models with more than one endpoint (or output) are mostly PKPD models,but these could also be, for example, PK binding models in which thereare measurements and predictions of both PK and its target.

ggPMX produces one diagnostics report per endpoint. As a consequence,the endpoint (if more than one) should be set at the time of theController creation in order to filter the observations dataset and tokeep only the values corresponding to the endpoint of interest. Thereare two ways of dealing with endpoints, using pmx_endpoint() (only forMonolix), or a simplified syntax for endpoints which is supported by allsoftware outputs.

To handle this, the user creates an “endpoint” object using the functionpmx_endpoint() having the following attributes:

• code (charcater): how the endpoint is coded in the input(modeling) dataset
• label: can be used in the title of some plots or for the reportfile name
• unit: used in the axis label for some plots
• files (list):list(predictions="predictions1",finegrid ="finegrid1")
• trans: whether the value must be transformed before beingdisplayed and which transformation to use.

To illustrate the Controller creation with multiple-endpoint models, abuilt-in PKPD example is used. The input dataset is called pk_pd.csvand has the following columns.

#>    id time amt  dv dvid   wt sex age#> 1 100  0.0 100   .    3 66.7   1  50#> 2 100  0.5   .   0    3 66.7   1  50#> 3 100  1.0   . 1.9    3 66.7   1  50#> 4 100  2.0   . 3.3    3 66.7   1  50#> 5 100  3.0   . 6.6    3 66.7   1  50#> 6 100  6.0   . 9.1    3 66.7   1  50

The dvid column contains values=3 for PK (first endpoint) and dose and=4 for PD (second endpoint). Monolix2016 outputs are found in folderRESULTS/ which contains predictions1.txt and finegrid1.txt for PKpredictions, and predictions2.txt and finegrid2.txt for PD predictions.The Endpoint and Controller objects are created as follows:

pkpd_path<- file.path(system.file(package="ggPMX"),"testdata","pk_pd")pkpd_work_dir<- file.path(pkpd_path,"RESULTS")pkpd_input_file<- file.path(pkpd_path,"pk_pd.csv")ep<- pmx_endpoint(code="4",label="some_label",unit="some_unit",file.code="2",# will use predictions2.txt and finegrig2.txttrans="log10")ctr<- pmx_mlx(directory=pkpd_work_dir,input=pkpd_input_file,dv="dv",dvid="dvid",endpoint=ep)#> use  predictions2.txt  as  model predictions file .#> use  finegrid2.txt  as  finegrid file .

For NONMEM and nlmixr users, endpoint can be simply specified in thecontroller creation by e.g. endpoint = 1

NONMEM

ctr<- pmx_nm(directory=nonmem_dir,file="run001.lst",endpoint=1## select the first endpointdvid="DVID"## use this column as observation id)

nlmixr

ctr<- pmx_nlmixr(fit,endpoint=1## select the first endpointdvid="DVID"## use this column as observation id)

Also for Monolix users, a simplified syntax for the Endpoint creationexists in the case where the endpoint code matches the files post-fixes(code=1 corresponds to predictions1.txt, code=2 corresponds topredictions2.txt). Instead of passing a pmxEndpoint object as argumentof the Controller, the user can specify the numerical valuecorresponding to the YTYPE/CMT/DVID column.

pmx_mlx(dvid="YTYPE",## use this column as observation idendpoint=1,## select the first endpoint...)## other pmx parameters , config, input,etc..

Internally, a pmxEndpoint object will be created, and observationshaving YTYPE=x will be filtered.

When fitting a multiple-endpoint model, Monolix is known to renameendpoint variables to y1, y2. For controller creation to work it maynecessary to reverse this renaming in the simulation file (e.g. y1 ->LIDV).

Besides the mandatory fields to initialize a Controller, the user canset optional parameters related to covariates. This feature isillustrated below with the Theophylline example.

theophylline_path<- file.path(system.file(package="ggPMX"),"testdata","theophylline")work_dir<- file.path(theophylline_path,"Monolix")input_data_path<- file.path(theophylline_path,"data_pk.csv")ctr<- pmx_mlx(directory=work_dir,input=input_data_path,dv="Y",cats= c("SEX"),conts= c("WT0","AGE0"),strats= c("STUD","SEX"))

Conts are the continuous covariates.Cats are categorical covariatesused in the model, whereasstrats are categorical variables that canbe used for plot stratification, but are not used as covariates in themodel.

The covariates can be accessed using helper functions:

ctr %>% get_cats()#> [1] "SEX"ctr %>% get_conts()#> [1] "WT0"  "AGE0"ctr %>% get_strats()#> [1] "STUD" "SEX"ctr %>% get_covariates()#> [1] "SEX"  "WT0"  "AGE0"

The content of the Controller can be seen by printing it:

ctr#>#> pmx object:#>#>#> |PARAM                |VALUE        |#> |:--------------------|:------------|#> |working directory    |theophylline |#> |Modelling input file |data_pk.csv  |#> |dv                   |Y            |#> |dvid                 |DVID         |#> |cats                 |SEX          |#> |conts                |WT0,AGE0     |#> |strats               |STUD,SEX     |#>#>#> |data_name   |data_file       |data_label                                |#> |:-----------|:---------------|:-----------------------------------------|#> |predictions |predictions.txt |model predictions file                    |#> |estimates   |estimates.txt   |parameter estimates file                  |#> |eta         |indiv_eta.txt   |invidual estimates of random effects file |#> |finegrid    |finegrid.txt    |finegrid file                             |#> |input       |data_pk.csv     |modelling input                           |#>#>#> |plot_name       |plot_type |#> |:---------------|:---------|#> |abs_iwres_ipred |SCATTER   |#> |abs_iwres_time  |SCATTER   |#> |iwres_ipred     |SCATTER   |#> |iwres_time      |SCATTER   |#> |iwres_dens      |PMX_DENS  |#> |iwres_qq        |PMX_QQ    |#> |npde_time       |SCATTER   |#> |npde_pred       |SCATTER   |#> |npde_qq         |PMX_QQ    |#> |dv_pred         |SCATTER   |#> |dv_ipred        |SCATTER   |#> |individual      |IND       |#> |eta_hist        |DIS       |#> |eta_box         |DIS       |#> |eta_matrix      |ETA_PAIRS |#> |eta_cats        |ETA_COV   |#> |eta_conts       |ETA_COV   |#> |eta_qq          |PMX_QQ    |

It contains three tables:

• The first table is the Controller configuration. The user can seethe working directory, the input modeling dataset name, thedependent variable (DV) name and other fields related to the model(e.g., continuous and discrete covariates).

• The second table lists theggPMX datasets. The first column(data_name) of this table contains the ggPMX name of the dataset;the second column (data_file) contains the names of the outputmodeling datasets (for example estimates.txt); in the third column(data_label) contains the dataset description.

• The third table provides the list of available plots in theGenerator. It corresponds to Table .

The Controller is a container that stores all plots. To get the list ofplots, the functionplot_names() is used:

ctr %>% plot_names()#>  [1] "abs_iwres_ipred" "abs_iwres_time"  "iwres_ipred"     "iwres_time"#>  [5] "iwres_dens"      "iwres_qq"        "npde_time"       "npde_pred"#>  [9] "npde_qq"         "dv_pred"         "dv_ipred"        "individual"#> [13] "eta_hist"        "eta_box"         "eta_matrix"      "eta_cats"#> [17] "eta_conts"       "eta_qq"

An alternative way to display the names of the existing plots is byprinting the content of the Controller as done above.

ggPMX provides a specialized function to create and update each plotpmx_plot_xx() wherexx is the plot name from the list above.

Each plot type is a class of similar plots.ggPMX defines thefollowing plot types:

• SCATTER: residual plots
• IND: individual plots display longitudinal (time course) predictionsand data (one panel per individual)
• DIS: distribution of empirical Bayes estimates (EBE) histogram orboxplot
• ETA_PAIRS: random effects (ETA) distributions and correlationsstructure
• ETA_COV: relationships between random effects (ETA) and continuousor categorical covariates
• PMX_QQ: quantile-quantile plot (qq-plot)

The following syntax allows to see which type of plot corresponds towhich plot name:

ctr %>% plots()#>           plot_name plot_type            plot_function#>  1: abs_iwres_ipred   SCATTER pmx_plot_abs_iwres_ipred#>  2:  abs_iwres_time   SCATTER  pmx_plot_abs_iwres_time#>  3:     iwres_ipred   SCATTER     pmx_plot_iwres_ipred#>  4:      iwres_time   SCATTER      pmx_plot_iwres_time#>  5:      iwres_dens  PMX_DENS      pmx_plot_iwres_dens#>  6:        iwres_qq    PMX_QQ        pmx_plot_iwres_qq#>  7:       npde_time   SCATTER       pmx_plot_npde_time#>  8:       npde_pred   SCATTER       pmx_plot_npde_pred#>  9:         npde_qq    PMX_QQ         pmx_plot_npde_qq#> 10:         dv_pred   SCATTER         pmx_plot_dv_pred#> 11:        dv_ipred   SCATTER        pmx_plot_dv_ipred#> 12:      individual       IND      pmx_plot_individual#> 13:        eta_hist       DIS        pmx_plot_eta_hist#> 14:         eta_box       DIS         pmx_plot_eta_box#> 15:      eta_matrix ETA_PAIRS      pmx_plot_eta_matrix#> 16:        eta_cats   ETA_COV        pmx_plot_eta_cats#> 17:       eta_conts   ETA_COV       pmx_plot_eta_conts#> 18:          eta_qq    PMX_QQ          pmx_plot_eta_qq

The diagnostic plots of ggPMX are generated by calling the functionspmx_plot_xx() wherexx is a placeholder for the plot name. The listof names of all available plots can be seen via:

ctr %>% plot_names()#>  [1] "abs_iwres_ipred" "abs_iwres_time"  "iwres_ipred"     "iwres_time"#>  [5] "iwres_dens"      "iwres_qq"        "npde_time"       "npde_pred"#>  [9] "npde_qq"         "dv_pred"         "dv_ipred"        "individual"#> [13] "eta_hist"        "eta_box"         "eta_matrix"      "eta_cats"#> [17] "eta_conts"       "eta_qq"          "pmx_vpc"

As a convention, when plots are described as ???Y vs. X???, it is meantthat Y is plotted on the vertical axis and X on the horizontal axis.

As an example, a plot of individual weighted residuals (IWRES) versustime (with default settings) can be generated using a single-line code:

ctr %>%pmx_plot_iwres_time

The complete list of available plots per plot type (given inparenthesis) is given below:

• Residual plots (SCATTER)

ctr %>%pmx_plot_dv_predctr %>%pmx_plot_dv_ipredctr %>%pmx_plot_iwres_timectr %>%pmx_plot_npde_timectr %>%pmx_plot_iwres_ipredctr %>%pmx_plot_abs_iwres_ipredctr %>%pmx_plot_npde_pred

• Empirical Bayes Estimates (EBE), also called “eta”, histogram andboxplot (DIS)

ctr %>%pmx_plot_eta_histctr %>%pmx_plot_eta_box

• Individual plots (IND)

ctr %>% pmx_plot_individual(which_pages=1)

• QQ-plots (PMX_QQ)

ctr %>%pmx_plot_npde_qqctr %>%pmx_plot_iwres_qq

• Distribution and correlation structure of random effects(ETA_PAIRS)

ctr %>%pmx_plot_eta_matrix

In order to generate VPCs a simulation dataset is requried. Creation ofVPC is slightly different dependening on the fitting software used(Monolix, NONMEM or nlmixr). Make sure to use the same name for the dvcolumn in the simulation file as the one used in the input (modelingdataset).

Monolix users, need to run a simulation with simulx. Here’s an examplecode

## Create simulated object using simulxmysim<- simulx(project=project_dir,nrep=100)### Retrieve simulated dataset (assumed to be in y1)simdata<-mysim$LIDV

For use with ggPMX, it is required that the IDs are reverted to theoriginal IDs as in the modelling dataset for ggPMX.

## Need to revert the original IDs as in modeling dataset for ggPMX## Rename IDs column to same name as in modeling dataset, e.g.## ???id??? in the example belowsimdata<-simdata %>%   mutate(newId= as.numeric(as.character(id))) %>%  left_join(.,mysim$originalId) %>%  mutate(id= as.numeric(as.character(oriId))) %>%  select(-oriId,-newId) %>%data.table::data.table()## It's highly recommended to store your simulation as .csvvpc_file<- write.csv(simdata,file="my_VPC.csv",quote=FALSE,row.names=FALSE)

pmx_sim creates a simulation object. It takes the following arguments:

Argumentss

1. file character path to the simulation file
2. irun character name of the simulation column
3. idv character name of the ind. variable
4. dv character name of the observation variable

Within pmx vpc controller, it is called like :

theoph_path<- file.path(  system.file(package="ggPMX"),"testdata","theophylline")WORK_DIR<- file.path(theoph_path,"Monolix")input_file<- file.path(theoph_path,"data_pk.csv")vpc_file<- file.path(theoph_path,"sim.csv")ctr<- pmx_mlx(directory=WORK_DIR,input=input_file,dv="Y",cats= c("SEX"),conts= c("WT0","AGE0"),strats="STUD",settings= pmx_settings(use.labels=TRUE,cats.labels=list(SEX=c("0"="Male","1"="Female")    )  ),sim= pmx_sim(file=vpc_file,irun="rep",idv="TIME"  ))

It is required to provide simulation tables to generate VPCs.Furthermore, it is highly recommended that simulation tables have a“sim”-suffix and are kept with the same naming convetion as for theprediction tables (e.g sdtab001sim, catab001sim)). In this case they’reautomatically recognized by the runnumber (runno) or by the model fileif specified there. For post-hoc simulation it is possible to include anadditional simfile:

ctr<- pmx_nm(directory=model_dir,file="modelfile.ctl"#or .lstsimfile="simulation_modelfile.ctl"#or .lst)

Important: When simulations are performed post-hoc and the controller isgenerated by run number, the reader might load the wrong .ext file (usedfor parameters). A warning message is displayed.

For nlmixr users, the VPC is generated automatically with the controllercreation and turned on by defaultvpc = TRUE.

ctr<- pmx_nlmixr(fit)## VPC will be generated automatically, vpc = TRUEctr<- pmx_nlmixr(fit,vpc=FALSE## But can be turned off)

The plot options are described in?pmx_plot_vpc function.

ctr %>%pmx_plot_vpc

By default the vpc plot ispercentile ; , but we can plot thescatter type:

ctr %>% pmx_plot_vpc(type="scatter")

ctr %>% pmx_plot_vpc(bin=pmx_vpc_bin(style="kmeans",n=5))

ctr %>% pmx_plot_vpc(labels= c(x="DV axis",y="TIME axis"))

ctr %>% pmx_plot_vpc(strat.facet=~SEX,facets=list(nrow=2))

User can customize the options to get a Monolix-like display.

ctr %>% pmx_plot_vpc(strat.facet=~SEX,facets=list(nrow=2),type="percentile",is.draft=FALSE,pi= pmx_vpc_pi(interval= c(0.1,0.9),median=list(color="green"),extreme=list(color="green")),obs= pmx_vpc_obs(color="blue",shape=18,size=2),ci= pmx_vpc_ci(interval= c(0.1,0.9),median=list(fill="red")))

A report (in pdf and docx format) containing all default diagnosticplots can be created using thepmx_report function. Theoutput cantake three different values:

• “report”: produces a pdf and a docx file (namedname.pdf andname.png specified in argumentname, located insave_dir)with default diagnostic plots
• “plots”: produces a folder named afterplots_subdir parameter(ggpmx_GOF by default, if parameter is not specified) located insave_dir that contains all default diagnostic plots, each in apdf and png file. The different plots are numerated in order to havean unique identifier for each plot (ex: ebe_box-1.pdf). This isnecessary for having correct footnotes that indicated the path tothe source file (for submission reports).
• “all”: is a combination of both options above.

Example:

ctr %>% pmx_report(name='Diagnostic_plots2',save_dir=work_dir,plots_subdir="ggpmx_report",output='all')

Note that running the same command first with the option“output=‘plots’” and then with the option “output=‘report’” willremove theggpmx_GOF folder.

Note also that by default, the report will have the DRAFT label on allplots. The label can be removed by using the settings argument in theController creation.

The user can customize the default report by creating a “template”. Tocreate a template, the user should create first a default report withthe following command:

ctr %>% pmx_report(name='Diagnostic_plots1',save_dir=work_dir,plots_subdir="ggpmx_report",output='report')

The Rmarkdown (.Rmd) file is the “template”. The user can modify theRmarkdown file as desired (ex: changing the size of some figures) andsave the modified file. The new template can be used with the followingcommand:

ctr %>% pmx_report(name='Diagnostic_plots3',save_dir=work_dir,plots_subdir="ggpmx_report",output='report',template=file.path(work_dir,'Diagnostic_plots1.Rmd'))

Any particular plot can be customized in two ways:

• Specifying options in each call of a plot (on the fly, recommended):

ctr %>% pmx_plot_xx(listofoptions)

• Customizing a type of plot for all subsequent calls in the sessionby modifying the Controller:

ctr %>% pmx_update(???xx???,listofoptions)

Help(pmx_gpar) lists some options.

Help(pmx_plot_xx) lists some possible parameters to update aparticular plot.

To obtain an exhaustive list of possible options for a particular plot,use the following:

ctr %>% get_plot_config("xx")

It is possible to visualize BLQ (below the limit of quantification)values in the individual plots. For this,bloq needs to be specifiedusing the pmx_bloq function (see example below with thepmx_mlxtran()function).

ctr %>% pmx_mlxtran(file_name=mlx_file,bloq=pmx_bloq(cens= ???BLQ???,limit= ???LIMIT???))ctr %>% pmx_plot_individual()

Monolix users may want to plot simulated BLQs. This is possible foroutputs with Monolix 2018 and later. An additional dataset is loaded(sim_blq), which will be used for plotting insted the regular“predictions”-dataset.

Thesim_blq statment can be specified within the plot (locally) …

ctr %>% pmx_mlxtran(file_name=mlx_file))ctr %>% pmx_plot_iwres_ipred(sim_blq=TRUE)

… or within the controller (globally). If this statment is used withinthe controller, all corresponding plots will plot simulated BLOQs bydefault.

ctr %>% pmx_mlxtran(file_name=mlx_file,sim_blq=TRUE))ctr %>% pmx_plot_iwres_ipred()

The user can define a Controller with global settings that will beapplied to all plots. For example remove draft annoataion, useabbreviation defintions to define axis labels, etc.

A settings object is defined by using the functionpmx_settings(). Thecreated object is passed as the parameter “settings” topmx(). Bydoing so, the settings are defined globally and apply to all plots. Fora complete list of global settings with their corresponding defaultvalues, please consult the ggPMX Help (?pmx_settings).

## set one or more settingsmy_settings<- pmx_settings(is.draft=FALSE,use.abbrev=TRUE,...)### set other settings parameters herectr<-  pmx_mlx(...,## put here other pmx parametessettings=my_settings  )

By default in the “standing” configuration, a DRAFT label is printed onall plots. In order to switch this label off, the user sets theis.draft option ofpmx_settings() toFALSE.

ctr<- theophylline(settings= pmx_settings(is.draft=FALSE))

The order of factors in a facet plot is determined by the data containedin the predictions data frame. Ordering of facet panels can be performedby changing the factor specification for facet columns in thepredictions data frame.

The example below defines a Controller with changed SEX facet panelsorder: 1 will be plotted before 0 rather than the other way round.

ctr<- theophylline()ctr[["data"]][["predictions"]][["SEX"]]<-factor(ctr[["data"]][["predictions"]][["SEX"]],levels=c("1","0"))pmx_plot_iwres_ipred(ctr,strat.facet=~SEX)

The standing configuration initializes all plots using abbreviations foraxis labels. Each abbreviation has its corresponding definition. To getthe list of abbreviation :

ctr %>%get_abbrev#> NULL

You can update one abbreviation to set a custom label

ctr %>% set_abbrev(TIME="TIME after the first dose")

Usinguse.abbrev flag you can use abbreviation definition to set axislabels:

ctr<- theophylline(settings=pmx_settings(use.abbrev=TRUE))ctr %>% set_abbrev(TIME="Custom TIME axis")ctr %>%pmx_plot_npde_time

within Monolix, user can choose to not use finegrid file even if it ispresent.

ctr<- theophylline()ctr %>% pmx_plot_individual(use.finegrid=FALSE,legend.position="top")#> USE predictions data set

In case of color startfication user can customize the legend. Forexample here using theggplot2::scale_color_manual:

ctr<- theophylline()ctr %>% pmx_plot_npde_time(strat.color="STUD")+ggplot2::scale_color_manual("Study",labels=c("Study 1","Study 2"),values=c("1"="green","2"="blue"))

Another way to do it is to define a globalscales.color parameter thatwill applied to all plots with strat.color :

ctr<- theophylline(settings=    pmx_settings(color.scales=list("Study",labels=c("Study 1","Study 2"),values=c("1"="orange","2"="magenta"))    ))ctr %>% pmx_plot_npde_time(strat.color="STUD")

ctr %>% pmx_plot_eta_box(strat.color="STUD")

In case of faceting by stratification user can redfine categoricallabels to have more human readables strips. Lables are defined withincats.labels argument and user can use them by settinguse.lables toTRUE.

ctr<- theophylline(settings=    pmx_settings(cats.labels=list(SEX=c("0"="M","1"="F"),STUD=c("1"="Study 1","2"="Study 2")      ),use.labels=TRUE    ))ctr %>% pmx_plot_npde_time(strat.facet=~SEX)

ctr<- theophylline(settings=    pmx_settings(cats.labels=list(SEX=c("0"="M","1"="F"),STUD=c("1"="Study 1","2"="Study 2")      ),use.labels=TRUE    ))ctr %>% pmx_plot_npde_time(strat.facet=~SEX)

ctr %>% pmx_plot_eta_box(strat.facet=~SEX)

The type of binning used for VPC plots can be applied as an update tothe controller object, causing the binning to be applied to subsequentVPC plots.

Updating the controller is achieved usingpmx_update in a commandlike:

pmx_update(ctr,"pmx_vpc",bin=pmx_vpc_bin(within_strat=TRUE,style="jenks"))

Subsequent VPC plotting calls will then use “jenks” style binning. Forexample:

ctr %>%pmx_plot_vpc

The functionpmx() is the generic function for greating a Controller.The user needs to specify a set of arguments such as the path to themodel directory, the software used for model fitting (Monolix ornlmixr), the name of a configuration. A list of all existingconfigurations is provided in the Appendix. Allmandatory argumentsofpmx() are listed in Table .

The example below defines a Controller with thestanding (standard)configuration.

theophylline_path<- file.path(system.file(package="ggPMX"),"testdata","theophylline")work_dir<- file.path(theophylline_path,"Monolix")input_data_path<- file.path(theophylline_path,"data_pk.csv")ctr<- pmx(sys="mlx",config="standing",directory=work_dir,input=input_data_path,dv="Y",dvid="DVID")

Note that the column “DVID” of data_pk.csv does not exist; however itis not needed here because there is only one single output of the model.As dvid is a mandatory argument, it still needs to be provided and wasset arbritrarly to “DVID” in the example above.

The input dataset can be provided to ggPMX via its location (as in theexample above) or as a data frame (maybe give an example). The modelingoutput datasets have to be in the location that is indicated as workingdirectory (work_dir in the example above).

ggPMX is compatible with Monolix versions 2016 and later, NONMEM version7.2 and later, and nlmixr.

In order to be able to produce all available diagnostic plots, thefollowing tasks should be executed in Monolix during the model fittingprocedure:

• Population parameters;
• Individual parameters (EBEs);
• Standard errors;
• Plots.

In addition, make sure to export charts data (In Monolix 2018: Settings-> Preferences -> Export -> switch on the Export charts data button).Select at least the following plots to be displayed and saved:individual fits and scatter plot of the residuals.

NONMEM Version 7.2/7.3/7.4 Preferred output tables according ???sdtab,patab, cotab, catab??? convention Simulation are required for creationof VPC (e.g. sdtab1sim)

Fit objects need to be generated by nlmixr and have data attached.Standard errors are required (a successful covariance plot) for fulldiagnostic checks. Can useboostrapFit() to get standard errorestimates if necessary

The main target of ggPMX is to create a report containing the followingplots (see abbreviation list below):

Abbreviations:

• NPDE: normalized prediction distribution errors
• IWRES: individual weighted residuals
• EBE: empirical Bayes estimates
• RE: random effects
• VPC: visual predivtive check

ggPMX implements few functions to generate and manipulate diagnosticplots. (Should we list pmx and pmx_mlx separately and say thedifferences? Or it’s maybe clear from the previous sections.)

(Apparently, it’s not the full list. Add all functions.) The design ofthe package is around the central object: the controller. It canintrospected or piped using the%>% operand.

Note that:

The controller is anR6 object, it behaves like a reference object.Some functions (methods) can have a side effect on the controller andmodify it internally. Technically speaking we talk about chaining notpiping here. However, usingpmx_copy user can work on a copy of thecontroller.

Graphical parameters inggPMX are set internally using thepmx_gparfunction. A number of graphical parameters can be set for the differentplot types.

args(pmx_gpar)#> function (is.title, labels, axis.title, which_pages, print, axis.text,#>     ranges, is.smooth, smooth, is.band, band, is.draft, draft,#>     discrete, is.identity_line, identity_line, smooth_with_bloq,#>     scale_x_log10, scale_y_log10, color.scales, is.legend, legend.position)#> NULL

More information can be found in the help document?pmx_gpar and inthe examples that follow.

For the moment we are mainly using standing configuration. In the nextrelease user can specfiy configuration either by cretaing a custom yamlfile or an R configuration object. Also ggPMX will create many helperfunctions to manipulate the configuration objects.

The shrinkage is a computation within controller data. In general it isused to annotate the plots. Although one can get it independently fromany plot usingpmx_comp_shrink. It is part of thepmx_compt_xxlayer( In the future we will add ,pmx_comp_cor ,pmx_comp_summary,..)

Here the basic call:

ctr %>%pmx_comp_shrink#>    EFFECT   OMEGA    SHRINK       POS FUN#> 1:     Cl 0.22485 0.1125175 0.2934250 var#> 2:      V 0.03939 0.9469996 0.0057915 var#> 3:     ka 0.10024 0.7423478 0.0810850 var

We get the shrinkage for each effect (SHRINK column).

The same values can be shown on distribution plot , for example :

ctr %>%pmx_plot_eta_box

or :

ctr %>%pmx_plot_eta_hist

You can add or remove shrinkage annotation usingis.shrink argument (TRUE by default) :

ctr %>%   pmx_plot_eta_box(is.shrink=FALSE)

You can compute shrinkage by applying either standard deviationsd orvariancevar :

ctr %>% pmx_comp_shrink(fun="var")#>    EFFECT   OMEGA    SHRINK       POS FUN#> 1:     Cl 0.22485 0.1125175 0.2934250 var#> 2:      V 0.03939 0.9469996 0.0057915 var#> 3:     ka 0.10024 0.7423478 0.0810850 var

ctr %>% pmx_plot_eta_box(shrink=pmx_shrink(fun="var"))

Note that for plotting functions which take a shrink parameter, this canbe created using thepmx_shrink function to create apmxShrinkClassobject (or it can be a list which can be converted into such an objectusingpmx_shrink).

Shrinkage can be applied after stratification :

ctr %>% pmx_comp_shrink(strat.facet=~SEX)#>    EFFECT SEX   OMEGA      SHRINK        POS FUN#> 1:     Cl   1 0.22485 -0.08032359 0.29342500 var#> 2:     Cl   0 0.22485  0.51828810 0.12378000 var#> 3:      V   1 0.03939  0.94628054 0.00579150 var#> 4:      V   0 0.03939  0.94818243 0.00437235 var#> 5:     ka   1 0.10024  0.70737008 0.08108500 var#> 6:     ka   0 0.10024  0.80907530 0.03676950 var

or by grouping like :

ctr %>% pmx_comp_shrink(strat.color="SEX")#>    EFFECT SEX   OMEGA      SHRINK        POS FUN#> 1:     Cl   1 0.22485 -0.08032359 0.29342500 var#> 2:     Cl   0 0.22485  0.51828810 0.12378000 var#> 3:      V   1 0.03939  0.94628054 0.00579150 var#> 4:      V   0 0.03939  0.94818243 0.00437235 var#> 5:     ka   1 0.10024  0.70737008 0.08108500 var#> 6:     ka   0 0.10024  0.80907530 0.03676950 var

We can

ctr %>% pmx_plot_eta_hist(is.shrink=TRUE,strat.facet=~SEX,facets=list(scales="free_y"))

or

ctr %>% pmx_plot_eta_box(is.shrink=TRUE,strat.facet="SEX",facets=list(scales="free_y",ncol=2))