In Open-Systems-Pharmacology/OSPSuite.ReportingEngine: R package to create reports for the Open Systems Pharmacology models

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
source(system.file("extdata", "vignette-helpers.R", package = "ospsuite.reportingengine"))

require(ospsuite.reportingengine)

This vignette focuses on time profiles and residuals plots (also called goodness of fit plots) generated using either mean model or population workflows.

Overview

The time profiles and residuals task first aims at generating time profile plots of specified simulation paths. The task can account for observed data, in which case residuals are calculated and plotted. For population workflows, aggregation functions are used, by default mean, median, 5th and 95th percentiles are computed and plotted in the time profiles.

Inputs

Results obtained from the simulate task and stored as csv files within the subdirectory SimulationResults from the workflowFolder directory are required in order to perform the goodness of fit (plotTimeProfilesAndResiduals) task. As a consequence, if the workflow output folder does not include the appropriate simulations, the task simulate needs to be active.

The objects SimulationSet (or PopulationSimulationSet for population workflows) and Output define which and how the simulated and, if available, observed data will be plotted.

Outputs

The time profiles and residuals task generates for each simulation set plots of:

• the simulated time profile
• in linear and logarithmic scale
• if available with data
• for total simulated time range.

Multiple applications

For simulations with multiple applications, up to 6 time profile plots can be generated corresponding to simulated time profiles in linear and logarithmic scale for:

• the total simulated time range
• the time range of the first application
• the time range of the last application

The input applicationRanges from SimulationSet objects can be used to select which ranges to include in the evaluation.

Residuals

If observed data are available in the range selected in the output definition and selected in dataSelection, they will be added to the time profile plots.

Furthermore, residuals will be calculated according to the residual scale defined by the Output objects included within the simulation sets and 6 additional plots will be created figures:

• Predicted vs observed in linear scale
• Predicted vs observed in log scale
• Residuals vs time
• Residuals vs simulated output
• Histogram of residuals
• qq-plot of residuals

The residual scale corresponds to one element from the enum ResidualScales:

ResidualScales

The default uses "Logarithmic" and calculate residuals as: $Residuals=log(Observed)-log(Simulated)$

If the option "Linear" is used instead, the residuals are calculated as: $Residuals=Observed-Simulated$

A histogram and a qq-plot of residuals across all the simulation sets will also be added at the end of the analysis.

For population workflows, residuals are calculated from the aggregation of the Simulated output. By default, $median(Simulated)$ is used instead of $Simulated$. It is possible to modify this aggregation function by updating the fields of the object AggregationConfiguration.

Template examples

The following sections will introduce template scripts including goodness of fit tasks as well as their associated reports. Table 1 shows the features tested in each template script.

examplesTable <- data.frame(
  Example = seq(1, 7),
  Workflow = c(rep("Mean model", 5), rep("Population model", 2)),
  `Simulation Sets` = c(1, 1, 1, 1, 2, 1, 1),
  Outputs = c(1, 1, 1, 2, 1, 1, 1),
  `Observed Data` = c("No", "Yes", "Yes", "Yes", "Yes", "No", "Yes"),
  check.names = FALSE
)
knitr::kable(examplesTable, caption = "Table 1: Features tested by each template script")

Example 1

Example 1 shows one of the simplest workflows including a goodness of fit. In this example, only a simulation time profile is performed since there is no observed data. Furthermore, only one simulation set and one output are requested.

Code

# Get the pkml simulation file: "MiniModel2.pkml"
simulationFile <- system.file("extdata", "MiniModel2.pkml",
  package = "ospsuite.reportingengine"
)

# Define the input parameters
outputA <- Output$new(
  path = "Organism|A|Concentration in container",
  displayName = "Concentration of A",
  displayUnit = "nmol/ml"
)

setA <- SimulationSet$new(
  simulationSetName = "A",
  simulationFile = simulationFile,
  outputs = outputA
)

# Create the workflow instance
workflow1 <-
  MeanModelWorkflow$new(
    simulationSets = setA,
    workflowFolder = "Example-1"
  )

# Set the workflow tasks to be run
workflow1$activateTasks(c("simulate", "plotTimeProfilesAndResiduals"))

# Run the workflow
workflow1$runWorkflow()

The output report for Example 1 is shown below. There are only 2 figures in the report (linear and log scales) showing the time profile of the output path 'Organism|A|Concentration in container' displayed as 'Concentration of A'.

cat(includeReportFromWorkflow(workflow1))

Example 2

Example 2 shows a workflow that used the simulation results performed in Example 1. This example shows 3 new features:

• The presence of observed data compared to the simulation
• The use of results obtained from a previous workflow
• The change of display unit without re-running the simulate task

Code

# Get the pkml simulation file: "MiniModel2.pkml"
simulationFile <- system.file("extdata", "MiniModel2.pkml",
  package = "ospsuite.reportingengine"
)

# Get the observation file and its dictionary
observationFile <- system.file(
  "extdata", "SimpleData.nmdat",
  package = "ospsuite.reportingengine"
)
dictionaryFile <- system.file(
  "extdata", "tpDictionary.csv",
  package = "ospsuite.reportingengine"
)
dataSource <- DataSource$new(
  dataFile = observationFile,
  metaDataFile = dictionaryFile
)

# Define the input parameters
outputA <- Output$new(
  path = "Organism|A|Concentration in container",
  displayName = "Simulated concentration of A",
  displayUnit = "µmol/ml",
  dataSelection = "MDV==0",
  dataDisplayName = "Observed concentration of A"
)

setA <- SimulationSet$new(
  simulationSetName = "A",
  simulationFile = simulationFile,
  outputs = outputA,
  dataSource = dataSource
)

# Create the workflow instance
workflow2 <-
  MeanModelWorkflow$new(
    simulationSets = setA,
    workflowFolder = "Example-1"
  )

# Set the workflow tasks to be run
workflow2$activateTasks("plotTimeProfilesAndResiduals")
workflow2$inactivateTasks("simulate")

# Run the workflow
workflow2$runWorkflow()

It can be seen from the run results that the task simulate was not run.

In this example, the content of the observed data file and its dictionary are respectively described Table 2 and 3.

obsTableExample2 <- readObservedDataFile(observationFile)
knitr::kable(obsTableExample2, caption = "Table 2: Observed data content")

dictTableExample2 <- readObservedDataFile(dictionaryFile)
knitr::kable(dictTableExample2, caption = "Table 3: Dictionary content")

The output report for Example 2 is shown below. There are now many more figures in the report showing not only time profiles in µmol/ml instead of nmol/ml but also residuals.

cat(includeReportFromWorkflow(workflow2))

Example 3

In example 3, the dictionary is changed and indicates that the lower limit of quantitation (LLOQ) is included in the observed data. Like the previous example, there is no need to re-run the simulate task.

Code

# Get the pkml simulation file: "MiniModel2.pkml"
simulationFile <- system.file("extdata", "MiniModel2.pkml",
  package = "ospsuite.reportingengine"
)

# Get the observation file and its dictionary
dataSource <- DataSource$new(
  dataFile = observationFile,
  metaDataFile = system.file(
    "extdata", "tpDictionaryLoq.csv",
    package = "ospsuite.reportingengine"
  )
)

# Define the input parameters
outputA <- Output$new(
  path = "Organism|A|Concentration in container",
  displayName = "Simulated concentration of A",
  displayUnit = "µmol/ml",
  dataSelection = "MDV==0",
  dataDisplayName = "Observed concentration of A"
)

setA <- SimulationSet$new(
  simulationSetName = "A",
  simulationFile = simulationFile,
  outputs = outputA,
  dataSource = dataSource
)

# Create the workflow instance
workflow3 <-
  MeanModelWorkflow$new(
    simulationSets = setA,
    workflowFolder = "Example-1"
  )

# Set the workflow tasks to be run
workflow3$activateTasks("plotTimeProfilesAndResiduals")
workflow3$inactivateTasks("simulate")

# Run the workflow
workflow3$runWorkflow()

The new dictionary content is described in Table 4.

dictTableExample2 <- readObservedDataFile(dictionaryFile)
knitr::kable(dictTableExample2, caption = "Table 4: Dictionary content")

The output report for Example 3 is shown below. The time profiles now indicate the LLOQ.

cat(includeReportFromWorkflow(workflow3))

# Remove the workflow folders
unlink(workflow3$workflowFolder, recursive = TRUE)

Example 4

In example 4, the simulation set now includes 2 output paths. Output path Organism|A|Concentration in container has observed data to be compared to whereas output path Organism|B|Concentration in container does not. Since output path Organism|B|Concentration in container was not requested by the examples 1, 2 and 3, it is not included in the simulation results. As a consequence, the simulate task needs to be re-run.

Code

# Get the pkml simulation file: "MiniModel2.pkml"
simulationFile <- system.file("extdata", "MiniModel2.pkml",
  package = "ospsuite.reportingengine"
)

# Get the observation file and its dictionary
dataSource <- DataSource$new(
  dataFile = observationFile,
  metaDataFile = dictionaryFile
)

# Define the input parameters
outputA <- Output$new(
  path = "Organism|A|Concentration in container",
  displayName = "Simulated concentration of A",
  displayUnit = "µmol/ml",
  dataSelection = "MDV==0",
  dataDisplayName = "Observed concentration of A"
)
outputB <- Output$new(
  path = "Organism|B|Concentration in container",
  displayName = "Simulated concentration of B",
  displayUnit = "µmol/ml"
)

setAB <- SimulationSet$new(
  simulationSetName = "A and B",
  simulationFile = simulationFile,
  outputs = c(outputA, outputB),
  dataSource = dataSource
)

# Create the workflow instance
workflow4 <-
  MeanModelWorkflow$new(
    simulationSets = setAB,
    workflowFolder = "Example-4"
  )

# Set the workflow tasks to be run
workflow4$activateTasks(c("simulate", "plotTimeProfilesAndResiduals"))

# Run the workflow
workflow4$runWorkflow()

The output report for Example 4 is shown below.

cat(includeReportFromWorkflow(workflow4))

# Remove the workflow folders
unlink(workflow4$workflowFolder, recursive = TRUE)

Example 5

In example 5, there are 2 simulation sets, each one includes one output path. Simulation set A includes output path Organism|A|Concentration in container and has observed data. Simulation set B includes output path Organism|B|Concentration in container and has no observed data. The simulate task needs to be re-run also in this case since the simulation results are saved in files that use the simulation set names.

Code

# Get the pkml simulation file: "MiniModel1.pkml"
simulationFileA <- system.file("extdata", "MiniModel1.pkml",
  package = "ospsuite.reportingengine"
)

# Get the pkml simulation file: "MiniModel2.pkml"
simulationFileB <- system.file("extdata", "MiniModel2.pkml",
  package = "ospsuite.reportingengine"
)

# Get the observation file and its dictionary
dataSource <- DataSource$new(
  dataFile = observationFile,
  metaDataFile = dictionaryFile
)

# Define the input parameters
outputA <- Output$new(
  path = "Organism|A|Concentration in container",
  displayName = "Simulated concentration of A",
  displayUnit = "µmol/ml",
  dataSelection = "MDV==0",
  dataDisplayName = "Observed concentration of A"
)
outputB <- Output$new(
  path = "Organism|B|Concentration in container",
  displayName = "Simulated concentration of B",
  displayUnit = "µmol/ml"
)

setA <- SimulationSet$new(
  simulationSetName = "A",
  simulationFile = simulationFileA,
  outputs = outputA,
  dataSource = dataSource
)
setB <- SimulationSet$new(
  simulationSetName = "B",
  simulationFile = simulationFileB,
  outputs = outputB
)

# Create the workflow instance
workflow5 <-
  MeanModelWorkflow$new(
    simulationSets = c(setA, setB),
    workflowFolder = "Example-5"
  )

# Set the workflow tasks to be run
workflow5$activateTasks(c("simulate", "plotTimeProfilesAndResiduals"))

# Run the workflow
workflow5$runWorkflow()

The output report for Example 5 is shown below.

cat(includeReportFromWorkflow(workflow5))

# Remove the workflow folders
unlink(workflow5$workflowFolder, recursive = TRUE)

Example 6

Example 6 shows a goodness of fit for a population workflow. In this example, there is no observed data, only one simulation set and one output are requested.

Code

# Get the pkml simulation file: "MiniModel2.pkml"
simulationFile <- system.file("extdata", "MiniModel2.pkml",
  package = "ospsuite.reportingengine"
)

populationFile <- system.file("extdata", "Pop500_p1p2p3.csv",
  package = "ospsuite.reportingengine"
)

# Define the input parameters
outputA <- Output$new(
  path = "Organism|A|Concentration in container",
  displayName = "Concentration of A",
  displayUnit = "µmol/ml"
)

setA <- PopulationSimulationSet$new(
  simulationSetName = "A",
  simulationFile = simulationFile,
  populationFile = populationFile,
  outputs = outputA
)

# Create the workflow instance
workflow6 <-
  PopulationWorkflow$new(
    workflowType = PopulationWorkflowTypes$parallelComparison,
    simulationSets = setA,
    workflowFolder = "Example-6"
  )

# Set the workflow tasks to be run
workflow6$activateTasks(c("simulate", "plotTimeProfilesAndResiduals"))

# Run the workflow
workflow6$runWorkflow()

The output report for Example 6 is shown below. The 2 figures of the report (linear and log scales) include median, mean 5th and 95th percentiles of the population.

cat(includeReportFromWorkflow(workflow6))

Example 7

Example 7 add observed data to the previous workflow.

Code

# Get the pkml simulation file: "MiniModel2.pkml"
simulationFile <- system.file("extdata", "MiniModel2.pkml",
  package = "ospsuite.reportingengine"
)
populationFile <- system.file("extdata", "Pop500_p1p2p3.csv",
  package = "ospsuite.reportingengine"
)

# Get the observation file and its dictionary
dataSource <- DataSource$new(
  dataFile = observationFile,
  metaDataFile = dictionaryFile
)

# Define the input parameters
outputA <- Output$new(
  path = "Organism|A|Concentration in container",
  displayName = "Simulated concentration of A",
  displayUnit = "µmol/ml",
  dataSelection = "MDV==0",
  dataDisplayName = "Observed concentration of A"
)

setA <- PopulationSimulationSet$new(
  simulationSetName = "A",
  simulationFile = simulationFile,
  populationFile = populationFile,
  outputs = outputA,
  dataSource = dataSource
)

# Create the workflow instance
workflow7 <-
  PopulationWorkflow$new(
    workflowType = PopulationWorkflowTypes$parallelComparison,
    simulationSets = setA,
    workflowFolder = "Example-6"
  )

# Set the workflow tasks to be run
workflow7$activateTasks(StandardPlotTasks$plotTimeProfilesAndResiduals)
workflow7$inactivateTasks(StandardSimulationTasks$simulate)

# Run the workflow
workflow7$runWorkflow()

The output report for Example 7 is shown below. The residuals are calculated on the median simulated values.

cat(includeReportFromWorkflow(workflow7))

# Remove the workflow folders
unlink(workflow7$workflowFolder, recursive = TRUE)

Open-Systems-Pharmacology/OSPSuite.ReportingEngine documentation built on May 1, 2024, 12:27 p.m.