--- version: "2024R1" language: "en" --- # MonolixSuite in R ## MonolixSuite in R Applications of MonolixSuite can be controlled via the R package *lixoftConnectors*. This package offers functions that allow users to script all the actions that can be performed via the graphical user interface. The package is not available on CRAN, but can be found in the installation folder of MonolixSuite. The package *Rsmlx* requires *lixoftConnectors* and contains features that are not available in the graphical user interface, such as automated PK model building and calculation of profile likelihood. The package *RsSimulx* is a wrapper for *lixoftConnectors* Simulx functions. It can be used to re-run the scripts made using the deprecated *mlxR* package, and to perform simulations with an alternative syntax. [![174857-20240913-110724.png](https://monolixsuite.slp-software.com/__attachments/a_a6f760e264b12133071d933c47d641667b0910a3b8c34b61e3649cfd81576568/174857-20240913-110724.png?cb=30c453b7f5fbdb09ea0cb42a5dc7a6e5)](https://www.linkedin.com/company/monolix-suite) [![1384060-20240913-110702.png](https://monolixsuite.slp-software.com/__attachments/a_2b503acc99b879b5b1a57e328f6e3a721c7cf9ae36169c0041d03d435e98cbf4/1384060-20240913-110702.png?cb=ca6d67e60f758d352745329b283e8f32)](https://www.youtube.com/@MonolixSuite) [![SLP-Logo-Cube_Full Color-20240827-160842.png](https://monolixsuite.slp-software.com/__attachments/a_e6441758a49ba8f42953aa788cf4fb9498af8950a4abb095992958bbeb595cc5/SLP-Logo-Cube_Full%20Color-20240827-160842.png?cb=c0c95ad1a1c411ba291f3cb92bc598ca)](https://www.simulations-plus.com/software/monolix/) ### Documentation * #### [Package "lixoftConnectors"](https://monolixsuite.slp-software.com/r-functions/2024R1/package-lixoftconnectors.md) * [Installation and initialization](https://monolixsuite.slp-software.com/r-functions/2024R1/installation-and-initialization.md) * [Errors and warnings](https://monolixsuite.slp-software.com/r-functions/2024R1/errors-and-warnings.md) * [Examples](https://monolixsuite.slp-software.com/r-functions/2024R1/examples.md) * [Reference](https://monolixsuite.slp-software.com/r-functions/2024R1/reference.md) * [Documentation of older versions](https://monolixsuite.slp-software.com/r-functions/2024R1/documentation-of-older-versions.md) * [1 more pages](https://monolixsuite.slp-software.com/r-functions/2024R1/package-lixoftconnectors.md) * #### [Package "Rsmlx"](https://monolixsuite.slp-software.com/r-functions/2024R1/package-rsmlx.md) * [Installation instructions](https://monolixsuite.slp-software.com/r-functions/2024R1/installation-instructions.md) * [Release notes](https://monolixsuite.slp-software.com/r-functions/2024R1/release-notes.md) * [Reference](https://monolixsuite.slp-software.com/r-functions/2024R1/reference-1.md) * #### [Package "RsSimulx"](https://monolixsuite.slp-software.com/r-functions/2024R1/package-rssimulx.md) * [Installation instructions](https://monolixsuite.slp-software.com/r-functions/2024R1/installation-instructions-rssimulx.md) * [Release notes](https://monolixsuite.slp-software.com/r-functions/2024R1/release-notes-rssimulx.md) * [Reference](https://monolixsuite.slp-software.com/r-functions/2024R1/reference-2.md) * #### [Package "mlxDesignEval"](https://monolixsuite.slp-software.com/r-functions/2024R1/package-mlxdesigneval.md) * [mlxDesignEval examples](https://monolixsuite.slp-software.com/r-functions/2024R1/mlxdesigneval-examples.md) * [Comparison to popED](https://monolixsuite.slp-software.com/r-functions/2024R1/comparison-to-poped.md) * #### [Package "conc-QTc"](https://monolixsuite.slp-software.com/r-functions/2024R1/package-for-conc-qtc-analysis.md) * [Download and Installation](https://monolixsuite.slp-software.com/r-functions/2024R1/download-and-installation.md) * [Introduction to conc-QTc analyses](https://monolixsuite.slp-software.com/r-functions/2024R1/introduction-to-conc-qtc-analyses.md) * [Conc-QTc modeling and implementation](https://monolixsuite.slp-software.com/r-functions/2024R1/conc-qtc-modeling-and-implementation.md) * [R functions documentation](https://monolixsuite.slp-software.com/r-functions/2024R1/r-functions-documentation.md) * [conc-QTc examples](https://monolixsuite.slp-software.com/r-functions/2024R1/conc-qtc-examples.md) * #### [Package "mlxModelFinder"](https://monolixsuite.slp-software.com/r-functions/2024R1/package-mlxmodelfinder.md) * [Parallel Execution](https://monolixsuite.slp-software.com/r-functions/2024R1/parallel-execution.md) * [Going Beyond the Built-in Model Libraries](https://monolixsuite.slp-software.com/r-functions/2024R1/going-beyond-the-built-in-model-libraries.md) * [Example: Implementing a custom cost function](https://monolixsuite.slp-software.com/r-functions/2024R1/using-custom-cost-functions.md) * [Changelog](https://monolixsuite.slp-software.com/r-functions/2024R1/changelog.md) * #### [Deprecated packages](https://monolixsuite.slp-software.com/r-functions/2024R1/deprecated-packages.md) * [Package "mlxR"](https://monolixsuite.slp-software.com/r-functions/2024R1/package-mlxr.md) --- version: "2024R1" language: "en" --- # addAdditionalCovariate ## \[Monolix - PKanalix\] Add an additional covariate Create an additional covariate for stratification purpose. Notice that these covariates are available only if they are not contant through the dataset. Available column transformations are: |----------------|--------------------------|-------------------------------------------------------------------| | \[continuous\] | 'firstDoseAmount' | (first dose amount) | | \[continuous\] | 'doseNumber' | (dose number) | | \[discrete\] | 'administrationType' | (admninistration type) | | \[discrete\] | 'administrationSequence' | (administration sequence) | | \[discrete\] | 'dosingDesign' | (dose multiplicity) | | \[continuous\] | 'observationNumber' | (observation number per individual, for a given observation type) | ### Usage R addAdditionalCovariate(transformation, base = "", name = "") ### Arguments transformation (character) applied transformation. base (character) \[optional\] base data on which the transformation is applied. name (character) \[optional\] name of the covariate. ### See also [`deleteAdditionalCovariate`](deleteadditionalcovariate) ### Examples R if (FALSE) { addAdditionalCovariate("firstDoseAmount") addAdditionalCovariate(transformation = "observationNumberPerIndividual", headerName = "CONC") } Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # addContinuousTransformedCovariate ## \[Monolix\] Add transformed continuous covariate Create a new continuous covariate by transforming an existing one. Transformed covariates cannot be use to produce new covariates. Call [`getCovariateInformation`](getcovariateinformation) to find out which covariates can be transformed. Those of type `"continuous"` can be used. ### Usage R addContinuousTransformedCovariate(...) ### Arguments ... Comma-separated pairs {transformedCovariateName = (character)"formula"} ### See also [`getCovariateInformation`](getcovariateinformation) get current covariates in the model [`addCategoricalTransformedCovariate`](addcategoricaltransformedcovariate) to add a transformation of a categorical covariate [`addMixture`](addmixture) to add a latent covariate [`removeCovariate`](removecovariate) remove added covariates ### Examples R initializeLixoftConnectors("monolix") #> [INFO] The library lixoftConnectors ("C:\Program Files\Lixoft\MonolixSuite2024R1\lib\lixoftConnectors.dll") is already loaded and initialized for monolix software -> nothing to be done. project_file <- file.path(getDemoPath(), "1.creating_and_using_models", "1.1.libraries_of_models", "theophylline_project.mlxtran") loadProject(project_file) #> [INFO] Results have been successfully loaded addContinuousTransformedCovariate( logtWEIGHT = "log(WEIGHT/70)" ) getCovariateInformation() #> $name #> [1] "WEIGHT" "logtWEIGHT" "SEX" #> #> $type #> WEIGHT logtWEIGHT SEX #> "continuous" "continuoustransformed" "categorical" #> #> $formula #> logtWEIGHT #> "logtWEIGHT = log(WEIGHT/70)" #> #> $range #> $range$WEIGHT #> [1] 54.6 86.4 #> #> $range$logtWEIGHT #> [1] -0.2484614 0.2104924 #> #> #> $categories #> $categories$SEX #> [1] "F" "M" #> #> #> $covariate #> id WEIGHT logtWEIGHT SEX #> 1 1 79.6 0.128518900 M #> 2 2 72.4 0.033711060 M #> 3 3 70.5 0.007117468 M #> 4 4 72.7 0.037846140 M #> 5 5 54.6 -0.248461400 F #> 6 6 80.0 0.133531400 M #> 7 7 64.6 -0.080280830 F #> 8 8 70.5 0.007117468 M #> 9 9 86.4 0.210492400 M #> 10 10 58.2 -0.184609900 F #> 11 11 65.0 -0.074107970 F #> 12 12 60.5 -0.145851900 F #> Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # addCustomNCAParametersFromPreferences ## \[PKanalix\] Add custom parameters from the preferences to the project If a custom NCA parameter exists in preferences, this function can be used to add it to the current project. Available arguments: |---------|-------------------------|------------------------| | "names" | (*character*, required) | Name of the parameter. | ### Usage R addCustomNCAParametersFromPreferences(names = NULL) ### See also [`createCustomNCAParameter`](createcustomncaparameter)`, `[`deleteCustomNCAParameter`](deletecustomncaparameter)`, `[`getCustomNCAParameters`](getcustomncaparameters) Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # addGroup ## \[Simulx\] Add simulation group Add a new simulation group. ### Usage R addGroup(group) ### Arguments group (character) Name of the group to add. ### Details Simulation groups can be added to the simulation [as in Simulx GUI](https://simulx.lixoft.com/simulation/simulation-scenario/). By default, the elements of the newly added group are the same as the first simulation group. To check which elements have been set for this group, please use [`getGroups`](getgroups). To change a group element, use [`setGroupElement`](setgroupelement). Note: when a Simulx project is created, a first group is created by default with the name "simulationGroup1". ### See also [`getGroups`](getgroups) ### Examples R # create two groups with different treatments initializeLixoftConnectors("simulx") #> [INFO] lixoftConnectors package is about to switch from "monolix" mode to "simulx" mode. #> Information relative to the previous "monolix" session won't be accessible anymore and potential unsaved project changes and associated results will be lost. #> Proceed software switch ? [y|N] project_name <- file.path(getDemoPath(), "4.exploration", "PKPD_exploration.smlx") loadProject(project_name) #> [ERROR] 'C:\Users\FranoMihaljevic\lixoft\monolix\monolix2024R1\demos\4.exploration\PKPD_exploration.smlx' is not a regular Monolix project. It must be loaded in Simulx addGroup("simulationGroup2") #> [ERROR] This function relates to "simulx" software. It is not available for "monolix" software. setGroupElement("simulationGroup2", elements = "Dose_4000") #> [ERROR] This function relates to "simulx" software. It is not available for "monolix" software. Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # addIndividualBlock ## \[Simulx\] Add an INDIVIDUAL block to a structural model Add an INDIVIDUAL block to a structural model. The structural model can be either a txt file or a library model, and the resulting model is written to another a file. ### Usage R addIndividualBlock(sourceModelFile, targetModelFile, force = FALSE) ### Arguments sourceModelFile (character) Path to the model file. Can be absolute or relative. Alternatively, the name of a library model (starting with "lib:"); see getLibraryModelName for a list of available models. targetModelFile (character) Path to the resulting file. Can be absolute or relative. force (logical) \[optional\] Should the resulting file be overwritten if already existing? Default: FALSE. ### Details Adding an [INDIVIDUAL block](https://simulx.lixoft.com/definition/model/) enables to add inter-individual variability to all parameters of a structural model. For every parameter theta available in the structural model file, population parameters theta_pop (typical value) and omega_theta (standard deviation of random effects) are added to the model, and lognormal distributions are used by default. If the structural model is a library model, the distributions are set to normal, lognormal or logitnormal depending on the meaning of the parameter. The resulting file can then be used as a starting point to remove IIV from some parameters, change parameter distributions, or add correlations if needed, and the final model can then be used to create a new Simulx project with [`newProject`](newproject). More details on the syntax of the INDIVIDUAL block can be found [here](https://mlxtran.lixoft.com/individual/). ### Examples R if (FALSE) { addIndividualBlock("lib:bolus_1cpt_VCl.txt", "example_output_bolus_1cpt_VCl_with_individual_block.txt") } # Working example: create project from a library model with IIV on all parameters. targetModelFile = tempfile("example_output_bolus_1cpt_VCl_with_individual_block",fileext = ".txt") addIndividualBlock("lib:bolus_1cpt_VCl.txt", targetModelFile) # we add an IIV block to a library model which has parameters V and Cl. #> [ERROR] Only available with simulx. file.show(targetModelFile) # to look at the created file with a text editor #> Warning: file.show(): file 'C:\Users\FRANOM~1\AppData\Local\Temp\RtmpopEtMJ\example_output_bolus_1cpt_VCl_with_individual_block49d056981fd5.txt' does not exist newProject(targetModelFile) # create a new Simulx project with this model #> [ERROR] Unexpected type encountered for the field "data", required. getPopulationElements() # a population parameters element was automatically created with V_pop, omega_V, Cl_pop, omega_Cl. #> [ERROR] This function relates to "simulx" software. It is not available for "monolix" software. #> NULL Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # addMixture ## \[Monolix\] Add latent covariate to the model for a finite mixture model Add a new latent covariate to the current model giving its name and its modality number (how many subpopulations). ### Usage R addMixture(...) ### Arguments ... A list of comma-separated pairs latentCovariateName = modalityNumber, where modalityNumber is an integer ### See also [`getCovariateInformation`](getcovariateinformation) get current covariates in the model [`addContinuousTransformedCovariate`](addcontinuoustransformedcovariate) to add a transformation of a continuous covariate [`addCategoricalTransformedCovariate`](addcategoricaltransformedcovariate) to add a transformation of a categorical covariate [`removeCovariate`](removecovariate) to remove added covariates ### Examples R initializeLixoftConnectors("monolix") #> [INFO] The library lixoftConnectors ("C:\Program Files\Lixoft\MonolixSuite2024R1\lib\lixoftConnectors.dll") is already loaded and initialized for monolix software -> nothing to be done. project_file <- file.path(getDemoPath(), "5.models_for_individual_parameters", "5.3.mixture_of_distributions", "PKgroup_project.mlxtran") loadProject(project_file) addMixture(lcat = 2) getCovariateInformation() #> $modalityNumber #> lcat #> 2 #> #> $name #> [1] "GROUP" "lcat" #> #> $type #> GROUP lcat #> "categorical" "latent" #> #> $categories #> $categories$GROUP #> [1] "0" "1" #> #> #> $covariate #> id GROUP #> 1 1 1 #> 2 2 0 #> 3 3 1 #> 4 4 0 #> 5 5 1 #> 6 6 0 #> 7 7 1 #> 8 8 0 #> 9 9 1 #> 10 10 0 #> 11 11 1 #> 12 12 0 #> 13 13 1 #> 14 14 0 #> 15 15 1 #> 16 16 0 #> 17 17 1 #> 18 18 0 #> 19 19 1 #> 20 20 0 #> 21 21 1 #> 22 22 0 #> 23 23 1 #> 24 24 0 #> 25 25 1 #> 26 26 0 #> 27 27 1 #> 28 28 0 #> 29 29 1 #> 30 30 0 #> 31 31 1 #> 32 32 0 #> 33 33 1 #> 34 34 0 #> 35 35 1 #> 36 36 0 #> 37 37 1 #> 38 38 0 #> 39 39 1 #> 40 40 0 #> 41 41 1 #> 42 42 0 #> 43 43 1 #> 44 44 0 #> 45 45 1 #> 46 46 0 #> 47 47 1 #> 48 48 0 #> 49 49 1 #> 50 50 0 #> 51 51 1 #> 52 52 0 #> 53 53 1 #> 54 54 0 #> 55 55 1 #> 56 56 0 #> 57 57 1 #> 58 58 0 #> 59 59 1 #> 60 60 0 #> 61 61 1 #> 62 62 0 #> 63 63 1 #> 64 64 0 #> 65 65 1 #> 66 66 0 #> 67 67 1 #> 68 68 0 #> 69 69 1 #> 70 70 0 #> 71 71 1 #> 72 72 0 #> 73 73 1 #> 74 74 0 #> 75 75 1 #> 76 76 0 #> 77 77 1 #> 78 78 0 #> 79 79 1 #> 80 80 0 #> 81 81 0 #> 82 82 0 #> 83 83 0 #> 84 84 0 #> 85 85 0 #> 86 86 0 #> 87 87 0 #> 88 88 0 #> 89 89 0 #> 90 90 0 #> 91 91 0 #> 92 92 0 #> 93 93 0 #> 94 94 0 #> 95 95 0 #> 96 96 0 #> 97 97 0 #> 98 98 0 #> 99 99 0 #> 100 100 0 #> Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # Advanced plot manipulation ## Adding horizontal or vertical lines on a plot The plot connectors return ggplot objects which can be further customized by adding other ggplot elements on top. In order to be able to use ggplot functions, the ggplot2 library must be loaded. Below we show several modifications on the prediction distribution plot. The following elements are changed directly in the plot settings: * number of bands for the prediction interval * addition of the observations * modification of the axis labels The following elements are modified using ggplot functions: * theme * position of axis ticks * addition of a horizontal line R library(lixoftConnectors) library(ggplot2) initializeLixoftConnectors() # load project loadProject(paste0(getDemoPath(),"/6.PK_models/6.3.multiple_doses/addl_project.mlxtran")) # run to ensure results are available runScenario() # generate default prediction disribution plot plotPredictionDistribution() # generate customized plot plotPredictionDistribution(settings=list(nbBands=1, obs=T, xlab="Time (hr)", ylab="Concentration (ng/mL)")) + theme_light() + theme(legend.position="none") + scale_x_continuous(breaks=c(0,24,48,72)) + geom_hline(yintercept=4, linetype="dashed", size = 1, color = "#9e2662") The default display when calling `plotPredictionDistribution()` is: ![image-20250110-091102.png](https://monolixsuite.slp-software.com/__attachments/a_3bdcdfad66b6165f68eeb339fd89a8dbb5b3306881994485251ffdeadd47e958/image-20250110-091102.png?cb=b16c70ba062eb733b01c31a40fc30f7a) After customization, the display is the following: ![image-20250110-091032.png](https://monolixsuite.slp-software.com/__attachments/a_fa8197463e6238d47903186f155cfb92b940045069ec7eb1c8efcd7fa18d223d/image-20250110-091032.png?cb=1ff5ced2b8dfc50f4154c595870fca33) ## Changing the order of subplots By default, the plot connectors will sort the output plots lexicographically (alphabetically). In this piece of code, we are using the ggplot2 facet_wrap function to rearrange the subplots, so that the subplot of male subjects comes before the subplot of female subjects (by default, since F is alphabetically before M, the first shown subplot is the subplot of female subjects). R library(lixoftConnectors) library(ggplot2) initializeLixoftConnectors() loadProject( file.path( getDemoPath(), "1.creating_and_using_models", "1.1.libraries_of_models", "theophylline_project.mlxtran" ) ) p <- plotObservedData( settings = list( mean = T, error = T, meanMethod = "geometric", timeAfterLastDose = T, ylog = T, lines = F, scales = "free_x", ncol = 2 ), stratify = list(color = c("SEX"), split = c("SEX")) ) p %+% facet_wrap( ~ factor(split, levels = c("#SEX:M", "#SEX:F")), labeller = p$facet$params$labeller, scales = "free_x" ) ## Renaming legend items In this case, we want to replace the legend items "SEX:M" and "SEX:F" with more meaningful strings "Male subjects" and "Female subjects". We will also reorder the legend items. We are using the ggplot2 function scale_color_manual. R library(lixoftConnectors) initializeLixoftConnectors() loadProject( file.path( getDemoPath(), "1.creating_and_using_models", "1.1.libraries_of_models", "theophylline_project.mlxtran" ) ) p <- plotObservedData( settings = list( mean = T, error = T, meanMethod = "geometric", timeAfterLastDose = T, ylog = T, lines = F, scales = "free_x", legend = T ), stratify = list(color = c("SEX")) ) levels(p$data$color) # get colors used in the legend p + scale_color_manual( name = "Observed data", breaks = c("#4682B4", "#4D4D4D"), values = c("#4682B4", "#4D4D4D"), labels = c("Female subjects", "Male subjects") ) ## Renaming subplots By running this piece of code, we will get the observed data plot with subplot names "#SEX:F" and "#SEX:M" which we would like to change to something more meaningful, like "Female" and "Male": R library(lixoftConnectors) initializeLixoftConnectors() loadProject( file.path( getDemoPath(), "1.creating_and_using_models", "1.1.libraries_of_models", "theophylline_project.mlxtran" ) ) p <- plotObservedData(stratify = list(split = c("SEX"))) p ![Rplot-20250115-160449.svg](https://monolixsuite.slp-software.com/__attachments/a_85acc4730c6642b0292c1e31c6f798ea15d3071c1f5e2ed768aa97d9bc9ab975/Rplot-20250115-160449.svg?cb=958afc7f6578381f5af76d365390e0ef) We can define a list with new labels and change the labeller function to use the strings from the list: new_labels <- list("#SEX:F" = "Female", "#SEX:M" = "Male") p$facet$params$labeller <- function(labels) { labels$split <- new_labels[labels$split] return(labels) } p ![Rplot01-20250115-160624.svg](https://monolixsuite.slp-software.com/__attachments/a_5c4b11e479154d8882bbc3e2272cbc4a2db2f6b0d509a0fed0189b134cd7ccae/Rplot01-20250115-160624.svg?cb=4d380bb48be50aef95d57ae07e121564) ## Dodging error bars By default, when error bars are displayed and several groups are merged on the same subplot, the error bars are not dodged and can be overlapping making them hard to read. R library(lixoftConnectors) initializeLixoftConnectors() library(ggplot2) loadProject(file.path( getDemoPath(), "1.creating_and_using_models", "1.1.libraries_of_models", "warfarinPK_project.mlxtran" ) ) p <- plotObservedData(obsName="y1", settings = list( mean = T, error = T, dots = F, meanMethod = "geometric", ylog = F, lines = F, binsSettings=list(is.fixedNbBins=T,nbBins=12)), stratify = list(split = c("sex"), mergedSplits=T) ) print(p) ![image-20250127-153811.png](https://monolixsuite.slp-software.com/__attachments/a_bcf6844ad84ef6472045c65868891de997129320dbd43787a48287682272cb24/image-20250127-153811.png?cb=a1c61cc9da1d3082b97f0a66bfadb1b1) In order to dodge the error bars, it is possible to edit the corresponding layer of the ggplot object `p`. Note that to use ggplot functions, it is necessary to load the ggplot2 package in addition to the lixoftConnectors. First inspect the layers to find the one that corresponds to the error bars: it is the one containing `ymin = ~errorMin`. Depending on what you have displayed on the plot, it can be the second or third layer for instance. R > p$layers [[1]] mapping: x = ~time, y = ~mean, colour = ~color geom_line: na.rm = FALSE, orientation = NA stat_identity: na.rm = FALSE position_identity [[2]] mapping: x = ~time, ymin = ~errorMin, ymax = ~errorMax, colour = ~color geom_errorbar: na.rm = FALSE, orientation = NA, width = 0.9 stat_identity: na.rm = FALSE position_identity Then change the position argument of that layer using `position_dodge()` and redraw the plot: p$layers[[2]]$position <- position_dodge(width = 2) print(p) ![image-20250127-154854.png](https://monolixsuite.slp-software.com/__attachments/a_86db722bf42f36bb887f6e3fe26dae36a24f6ccd590ae6bcbcdb040f19e0cfe0/image-20250127-154854.png?cb=c8d71e644d992f8d0ab46eedf9a08a60) ## Ajusting axis limits for Individual fits Individuals fits are often displayed using **log scale for the y-axis**. Due to the model prediction being usually zero at time zero, by default the y-axis limits span up to 1e-16, which makes the plot hard to read: R plotIndividualFits( settings = list(cens = TRUE, popCov = TRUE, legend = FALSE, ylog = TRUE, xlab = "Time (hr)", ylab = "Concentration (ng/mL)")) ![image-20250912-161253.png](https://monolixsuite.slp-software.com/__attachments/a_f76f28bc315b210db6fa169cf4d5593b125d0b8a087a145d7acca2f4e4e7f085/image-20250912-161253.png?cb=8b473b96ea7527240f440ba6c4971684) A possible solution is to specify the y-axis limits using the `ylim` argument, but this will apply to all individuals, which may be inconvenient in case of different dose groups and individuals concentration spanning over different order of magnitude. R plotIndividualFits( settings = list(cens = TRUE, popCov = TRUE, legend = FALSE, ylog = TRUE, xlab = "Time (hr)", ylab = "Concentration (ng/mL)", ylim = c(0.1,1000))) ![image-20250912-161644.png](https://monolixsuite.slp-software.com/__attachments/a_ad89ad408821ee0cefd819f75d5ae4448f658c5e11e082028df50f2113cf8b09/image-20250912-161644.png?cb=08f0be09ffde481f1996dbbd8bcf703d) Another option is to leave the upper bound of the `ylim` argument undefined (as `NA`), to let the plot choose automatically the upper limit for each individual based on the data and prediction, while enforcing the lower limit to be a specified value. R plotIndividualFits( settings = list(cens = T, popCov = T, legend = F, ylog = TRUE, xlab = "Time (hr)", ylab = "Concentration (ng/mL)", ylim = c(0.1,NA))) ![image-20250912-161809.png](https://monolixsuite.slp-software.com/__attachments/a_15c726e57cdf09b1d320a1d68ca088798ec8aab158711b005eb18d73bfefe34a/image-20250912-161809.png?cb=5c48b6ea02fe895c076f687b3fe933ca) Last updated: September 12, 2025 --- version: "2024R1" language: "en" --- # AI Generator of Custom Libraries Last updated: March 11, 2026 --- version: "2024R1" language: "en" --- # applyFilter ## \[Monolix - PKanalix\] Apply filter Apply a filter on the current data. ### Usage R applyFilter(filter, name = "") ### Arguments filter (list\< list\< action = "headerName-comparator-value" \> \> or "complement") filter definition. Existing actions are "selectLines", "selectIds", "removeLines" and "removeIds". First vector level is for set unions, the second one for set intersection. It is possible to give only a list of actions if there is only no high-level union. name (character) \[optional\] created data set name. If not defined, the default name is "currentDataSet_filtered". ### Details The possible actions are line selection (selectLines), line removal (removeLines), Ids selection (selectIds) or removal (removeIds). The selection is a string containing the header name, a comparison operator and a value selection = "headerName\*-comparator\*\*-value" (ex: "id=='100'", "WEIGHT\<70", "SEX!='M'") Notice that : - The headerName corresponds to the data set header or one of the header aliases defined in MONOLIX software preferences - The comparator possibilities are "==", "!=" for all types of value and "\<=", "\<", "\>=", "\>" only for numerical types Syntax: \* apply a simple filter: applyFilter( filter = list(act = sel)), e.g. applyFilter( filter = list(removeIds = "WEIGHT\<50")) =\> apply a filter with the action act on the selection sel. In this example, we apply a filter that removes all subjects with a weight less than 50. \* apply a filter with several concurrent conditions, i.e AND condition: applyFilter( list(act1 = sel1, act2 = sel2)), e.g. applyFilter( filter = list(removeIds = "WEIGHT\<50", removeIds = " AGE\<20")) =\> apply a filter with both the action act1 on sel1 AND the action act2 on sel2. In this example, we apply a filter that removes all subjects with a weight less than 50 and an age less than 20. It corresponds to the intersecton of the subjects with a weight less than 50 and the subjects with an age less than 20. \* apply a filter with several non-concurrent conditions, i.e OR condition: applyFilter(filter = list(list(act1 = sel1), list(act2 = sel2)) ), e.g. applyFilter( filter = list(list(removeIds = "WEIGHT\<50"),list(removeIds = " AGE\<20"))) =\> apply a filter with the action act1 on sel1 OR the action act2 on sel2. In this example, we apply a filter that removes all subjects with a weight less than 50 and an age less than 20. It corresponds to the union of the subjects with a weight less than 50 and the subjects with an age less than 20. \* It is possible to have any combination: applyFilter(filter = list(list(act1 = sel1), list(act2 = sel2, act3 = sel3)) ) \<=\> act1,sel1 OR ( act2,sel2 AND act3,sel3 ) \* It is possible to apply the complement of an existing filter: applyFilter(filter = "complement") ### See also [`getAvailableData`](getavailabledata) [`createFilter`](createfilter) [`removeFilter`](removefilter) ### Examples Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # Basic PKanalix workflow In the example below we show how to create a new PKanalix project from a dataset, choose the NCA settings, run the NCA parameter calculation and retrieve the results. In addition, "table1" R package is used to display a nicely formatted table. R library(lixoftConnectors) initializeLixoftConnectors(software="pkanalix") # get path to dataset from demos data_file <- paste0(getDemoPath(),"/2.case_studies/data/Theo_extravascular_singledose.csv") # Create new PKanalix project newProject(data = list(dataFile = data_file, headerTypes = list(ID='id',time='time',DV="observation", Amount="amount", FORM="catcov",AGE="contcov", HT="contcov",SEQ="catcov",Period="occ"))) # Choose settings for NCA task setNCASettings(administrationtype = list("1"="extravascular"), integralMethod = "LinLogTrapLinLogInterp", blqMethodAfterTmax="missing", lambdaRule="adjustedR2", partialAucTime=list(T,c(0,24)), ajdr2AcceptanceCriteria=list(T,0.97), computedNCAParameters=c("AUCINF_obs","Cmax","Tmax")) # Save project saveProject(projectFile = "theo_NCA.pkx") # Compute NCA runNCAEstimation() # Select NCA parameters from results indivParams = getNCAIndividualParameters("AUCINF_obs","Cmax")$parameters print(head(indivParams)) # Post-processing of the results: # summarizing the results split by formulation in a table library(table1) table1::label(indivParams$AUCINF_obs) <- "AUC to infinity" table1::label(indivParams$Cmax) <- "Cmax" table1::table1(~AUCINF_obs + Cmax | FORM, data = indivParams) The generated table is: ![image-20241021-125519.png](https://monolixsuite.slp-software.com/__attachments/a_c2fa2c6b1ae6ffa7106059a755a12387110d41449ed2571206d4969639ad024a/image-20241021-125519.png?cb=c670b8a1821c57e037770d576a3324b8) Last updated: February 21, 2025 --- version: "2024R1" language: "en" --- # Basic Simulx workflow ## Load a project and run the simulation Simple example to load an existing project, run the simulation and look at the results: R # load and initialize the API library(lixoftConnectors) initializeLixoftConnectors(software="simulx") # Get the project project <- paste0(getDemoPath(), "/2.models/longitudinal.smlx") loadProject(projectFile = project) # Run the simulation runSimulation() # The results are accessible through the function getSimulationResults() # The results for the output is the list res and TS is one of the outputs head(getSimulationResults()$res$TS) id time TS 1 1 0 10.00000 2 1 1 11.04939 3 1 2 12.20862 4 1 3 13.48915 5 1 4 14.90359 6 1 5 16.46585 ## Import a Monolix project and simulate a new output variable with the EBEs In this example, the Monolix demo theophylline_project.mlxtran is imported to Simulx. A new variable AUC is added to the model, and it is simulated on a regular grid using the EBEs estimated in the Monolix project. This requires that the tasks "Population parameters" and "EBEs" have run in the Monolix project. R # Import a Monolix project and simulate a new output variable with the EBEs ===== # Load and initialize library ===== library(lixoftConnectors) initializeLixoftConnectors(software = "monolix") #Load project from Monolix Demos and run it to get EBEs estimated ===== MonolixProject <- paste0(getDemoPath(), "/1.creating_and_using_models/1.1.libraries_of_models/theophylline_project.mlxtran") loadProject(MonolixProject) runScenario() initializeLixoftConnectors(software = "simulx", force = TRUE) #IMPORT Monolix project and check that mlx_EBEs element has been created ===== importProject(MonolixProject) getIndividualElements() #SET ADDITIONAL formula in the model to compute the AUC ===== setAddLines("ddt_AUC = Cc") # define a new element output with a regular grid for the variable "AUC" defineOutputElement(name="AUCregular", element = list( data = data.frame( start = 0, interval = 1, final = 120), output = "AUC")) #CHECK the simulation group that already exists ===== getGroups() # => we currently have one group, called "simulationGroup1", which uses the following elements: # - mlx_Pop (population parameter values) => need to be replaced by mlx_EBEs # - mlx_Adm1 (treatment from monolix data set) => OK # - mlx_CONC (output CONC with times as in monolix data set) => need to be replaced by AUC with new grid # - size = 12 (same number of in original data set) => OK #MODIFY the existing simulation group such that it uses EBEs and the new output element===== setGroupElement(group = "simulationGroup1", elements = c("mlx_EBEs", "AUCregular")) # check the new simulation setup getGroups() #Set shared IDs to say the individuals in the treatment element and in the parameter element are the same and should alwyas be sampled together ===== setSharedIds(c("treatment","individual")) # check the simulation setup getGroups() getSharedIds() # - mlx_EBEs => ok # - remaining parameters (error model parameters if residual error is simulated) => not needed as we output model predictions # - mlx_Adm1 => ok # - AUCregular => ok # - size=12 => ok #SAVE and RUN project ===== saveProject("simulx_api_example.smlx") # run simulation runSimulation() #GET the results: simulated values for AUC and sampled individual parameters ===== simulatedParam <- getSimulationResults()$IndividualParameters$simulationGroup1 simulatedOutput <- getSimulationResults()$res$AUC #PLOT the results ===== library(ggplot2) ggplot(data = simulatedOutput, aes(x=time, y=AUC))+ geom_line() + aes(color = factor(as.double(original_id)))+ labs(colour="Original ID") We get the following simulations for the new output variable. "Original_id" is the ID used in the original dataset loaded in the Monolix project. ![2023-03-13_13h36_08-20250221-105339.png](https://monolixsuite.slp-software.com/__attachments/a_23d84267eb640fa5d1dbf43573ab809caceae9579f52d168c50009dd1eb2f0b0/2023-03-13_13h36_08-20250221-105339.png?cb=a12baa4492facf078bf5f36888cdafa2) ## Create a Simulx project from scratch and simulate covariate-dependent treatments This script creates a project similar to the demo project *treatment_weight_and_genotype_based.smlx* , available in the folder 3.1.treatments of Simulx demos. The model file *TMDDmodel.txt* can be downloaded here: [TMDDmodel.txt](https://monolixsuite.slp-software.com/__attachments/a_2882a5090b42732378114060057f71436ac922a7416cceea58fca987e137c6cf/TMDDmodel.txt.md?cb=b23f2259db58c1b05bdab7ad198af59c). While in the interface of Simulx it is possible to directly define distribution laws for covariates and covariate-dependent treatments that are computed at the simulation step, in R it is necessary to sample covariates from distributions and derive the corresponding doses before defining the elements for the simulation. R ############################################################### # Create a Simulx project from scratch and simulate covariate-dependent treatments ############################################################### # load and initialize the API library(lixoftConnectors) initializeLixoftConnectors(software = "simulx") # Create a new project based on the model newProject(modelFile = "TMDDmodel.txt") # Define vector of population parameters definePopulationElement(name="PopParameters", element = data.frame(V_pop=3, omega_V=0.1 ,Cl_pop=0.1, omega_Cl=0.1, Q_pop=1, omega_Q=0.2, V2_pop=3, omega_V2=0.2, beta_V_logtWeight=1, KD_pop=0.01, omega_KD=0.01, R0_pop=0.2, omega_R0=0.3, kint_pop=50, omega_kint=0.2, kon_pop=10, omega_kon=0.01, ksyn_pop=10, omega_ksyn=0.2, beta_Cl_logtWeight=0.75, beta_Q_logtWeight=0.75, beta_V2_logtWeight=0.75, beta_KD_Genotype_Heterozygous=1.3)) # Define covariates: in the interface it is possible to define distribution laws, # but in R we have to sample from the distributions prior to the simulation. # Here we define 3 tables for the 3 simulation groups for (i in 1:3){ covTable <- data.frame(id=1:250, Weight = rlnorm(n=250, mean=log(70), sd=0.3), Genotype=c("Homozygous", "Heterozygous")[sample(1:2, 250, replace=T)]) write.csv(covTable, file = paste0("covTable",i,".csv"), quote=F, row.names=F) defineCovariateElement(name=paste0("covTable",i), element = paste0("covTable",i,".csv")) } # Define common treatment defineTreatmentElement(name="1000nmol", element=list(admID=1, data=data.frame(time=seq(0,by=21,length.out = 5), amount=1000, tInf=0.208))) # Define weight-based treatment: in the interface it can be done directly with a scaling formula, # but in R we have to define a table of individual doses prior to the simulation. # We use covTable2.csv for simulation group 2. covTable <- read.csv("covTable1.csv") trt_14nmolPerKg <- data.frame(id=rep(1:250, each=5), time=rep(seq(0,by=21,length.out = 5), times=250), amount=14*covTable$Weight, tInf = 0.208) write.csv(trt_14nmolPerKg, file="trtTableWeight.csv", quote=F, row.names=F) defineTreatmentElement(name="14nmolPerKg", element=list(admID=1, data="trtTableWeight.csv")) # Define genotype-based treatment: in the interface it can be done directly with a scaling formula, # but in R we have to define a table of individual doses prior to the simulation. # We use covTable3.csv for simulation group 3. covTable <- read.csv("covTable3.csv") trt_genotype <- data.frame(id=rep(1:250, each=5), time=rep(seq(0,by=21,length.out = 5), times=250), amount=ifelse(covTable$Genotype=="Heterozygous", 1000, 800), tInf = 0.208) write.csv(trt_genotype, file="trtTableGenotype.csv", quote=F, row.names=F) defineTreatmentElement(name="1000nmolHomo_800nmolHetero", element=list(admID=1, data="trtTableGenotype.csv")) # Define outputs defineOutputElement(name="Concentration", element = list(output="L", data=data.frame(time=seq(0,96,by=1)))) defineOutputElement(name="TargetOccupancy", element = list(output="TO", data=data.frame(time=seq(0,96,by=1)))) # By default there is a single simulation group named simulationGroup1, # change its name and set all elements for the simulation setGroupSize("simulationGroup1", 250) setGroupElement("simulationGroup1", elements = c("PopParameters","14nmolPerKg","Concentration","TargetOccupancy", "covTable1")) renameGroup("simulationGroup1", "Weight_based") # Define the two other simulation groups with different covariates and treatments. # By default the other elements are the same as the first group addGroup("Flat_dose") setGroupElement("Flat_dose", elements = c("1000nmol","covTable2")) addGroup("Genotype_based") setGroupElement("Genotype_based", elements = c("1000nmolHomo_800nmolHetero","covTable3")) # Save project and run the simulation saveProject(projectFile = "treatment_weight_and_genotype_based.smlx") runSimulation() # Get simulation results: a table of individual parameters for each group, and a table of simulated values for each output sim <- getSimulationResults() names(sim$IndividualParameters) # "Weight_based", "Flat_dose", "Genotype_based" names(sim$res) # "L", "TO" ## Import a Monolix project and simulate with new output times Monolix simulates the original dataset with replicates to generate the VPC. It can be useful to perform the same simulations in Simulx with modified observation times, for example to correct the bias in VPC due to dropout. The example below shows this step, and uses the project PDTTE_dropout. The full example to correct a VPC for dropout is detailed [on this page](https://monolixsuite.slp-software.com/monolix/2024R1/correcting-the-vpc-for-missing-observations.md). In this example the Monolix project with the biased VPC is named PD_TTE.mlxtran, the new regular measurement times are similar to the measurements in the original dataset, but are not impacted by dropout, and the number of replicates is 500 (as the default number of simulations for the VPC in Monolix). The random variable capturing dropouts in the model used in Monolix is called survival. R ############################################################### # Import a Monolix project and simulate with new output times ############################################################### library(lixoftConnectors) initializeLixoftConnectors(software = "simulx") importProject(projectFile = "PDTTE_dropout.mlxtran") defineOutputElement(name="y1_alltimes", element = list(output="y1", data=data.frame(time=seq(-100,1500,by=100)))) setGroupElement("simulationGroup1", elements = c("y1_alltimes","mlx_survival")) setNbReplicates(500) saveProject(projectFile = "PD_TTE_alltimes.smlx") setPreferences(exportsimulationfiles=T) runSimulation() Last updated: February 21, 2025 --- version: "2024R1" language: "en" --- # Bayesian individual dynamics predictions This example shows how to perform Bayesian individual dynamic predictions. Using an existing population PK model with estimated population parameters as prior, samples from the conditional distribution (i.e posterior) are obtained for each individual of a new dataset, to capture the uncertainty of the individual parameters given the observations. ## Step 1: sampling from the conditional distribution We will use the demo Theophylline project as an example, using `loadProject()`. If population parameters have not been estimated yet, the task must be run with `runPopulationParameterEstimation()`. R library("lixoftConnectors") library(dplyr) library(ggplot2) initializeLixoftConnectors(software = "monolix", force=T) # load demo project project <- paste0(getDemoPath(), "/1.creating_and_using_models/1.1.libraries_of_models/theophylline_project.mlxtran") loadProject(projectFile = project) # run parameter estimation (if needed) runPopulationParameterEstimation() In order to use the estimated population parameters as a prior (and not re-estimate them based on the new dataset), we update the initial values based on the previous estimates with `setInitialEstimatesToLastEstimates()` and then fix them using `setPopulationParameterInformation()`. R setInitialEstimatesToLastEstimates() popparams <- getPopulationParameterInformation() popparams$method <- "FIXED" setPopulationParameterInformation(popparams) The next step is to load the new dataset containing the individuals for which we want to perform the Bayesian prediction. In this example, we reuse the original dataset but truncate the time at 10h, to investigate the impact of a shorter trial duration. The observations from this new dataset are also saved in a R object which will be used later to overlay the observations of the prediction plot. R baseDataset <- read.table(getData()$dataFile, header=T) filteredDataset <- baseDataset[baseDataset$TIME<10,] write.csv(filteredDataset, "filteredDataset.csv", row.names = F, quote=F) setData("filteredDataset.csv", getData()$headerTypes, getData()$observationTypes) # save for plotting later obsData <- getObservationInformation()$CONC By default only a few samples from the conditional distribution are saved. In order to obtain a nice prediction interval, we will save 200 samples. This is done by using `nbsimulatedparameters=200` in `setConditionalDistributionSamplingSettings()`. The samples are taken uniformly over an interval of `nbminiterations` iterations. Given that two samples next to each other have a risk to be identical, we usually set nbminiterations = 3 x nbsimulatedparameters to use one sample every 3 iterations. When `nbminiterations` is large, the convergence criteria are harder to achieve. To avoid running for very long, we can give a max number of iterations with `enableMaxIterations=T` and `nbMaxIterations=1200`. The call is done in two steps due avoid being temporarily in a situation where `nbminiterations`\>`nbMaxIterations` (which is forbidden). R setConditionalDistributionSamplingSettings(enableMaxIterations=T, nbMaxIterations=1200) setConditionalDistributionSamplingSettings(nbminiterations=600, nbsimulatedparameters=200) The project can then be saved (to avoid overwriting the original run) and run. Even if the population parameters are fixed, the task Pop param must be run before running the conditional distribution task. We then retrieve the samples from the conditional distribution with `getSimulatedIndividualParameters()`. R saveProject("bayesian_forecasting.mlxtran") runPopulationParameterEstimation() # mandatory before other tasks, but nb of iterations is null as all parameters are fixed runConditionalDistributionSampling() # this is sampling from the posterior distribution for each individual # retrieve the samples from the conditional distribution samplesCondDistrib <- getSimulatedIndividualParameters() The `samplesCondDistrib` contains one line per id and per sample (column rep). R > head(samplesCondDistrib) rep id ka V Cl 1 1 1 1.9582949 0.3939064 0.02888324 2 1 2 2.1925275 0.4230115 0.05076061 3 1 3 1.6761378 0.4503078 0.04192262 4 1 4 1.2335670 0.5261603 0.03801980 5 1 5 1.5879271 0.5561964 0.04986806 6 1 6 0.9479706 0.4571833 0.04743354 ## Step 2: simulation using the samples for each individual To simulate the prediction corresponding to the 200 samples of each individual, we will use Simulx. For this, we will define an Individual parameter element containing the samples. However, this type of element can have an "id" column but not both an "id" and a "rep" column. We will thus work replicate per replicate. Note that it would not be possible to merge the id and rep column into a single one, because then the ids would be renamed and the matching to the id column of the treatment or regressor elements would be broken. The first step is to export from Monolix to Simulx. R exportProject(settings=list(targetSoftware="simulx"), force=T) By default, a simulation is already setup with the same number of individuals (size), treatment (mlx_Adm1), covariates (if applicable) and regressors (if applicable) as in the dataset used in the Monolix project. The population parameters are used by default, and we will change this later to use the sampled individual parameters. R > printSimulationSetup() $simulationGroups simulationGroup1 size 12 parameters mlx_Pop treatments mlx_Adm1 outputs mlx_CONC The simulation is performed into a for loop going over the replicates (i.e samples). Within each iteration of the loop, all individuals are simulated with sample i. To define the individual parameter element, we need to provide a csv file. We thus filter the `samplesCondDistrib`to keep only replicate i and save this as a csv file. The individual parameters element is defined in Simulx using `defineIndividualElement()`. In addition, we can define an output element to control de time grid of the prediction using `defineOutputElement()`. In this example, we use the smooth (without residual error) model output called `Cc`. The two created element, called `samplesCondDistrib` and `outCc` are applied to the simulation using `setGroupElement(group = "simulationGroup1", elements = c("samplesCondDistrib","outCc"))`. Note that `simulationGroup1` is the default name, as visible above when using `printSimulationSetup()`. To ensure that the ids of the individual parameter and other elements (here treatment mainly) are correctly mapped together, we use `setSharedIds(sharedIds = c("output", "treatment", "regressor", "individual"))`. The simulation is then run with `runSimulation()` and the results are retrieved with `getSimulationResults().` The replicate information is added to the result data.frame and the results for all replicates are merged. R sims <- NULL for(i in unique(samplesCondDistrib$rep)){ cat(paste0("Replicate ", i, " out of ", length(unique(samplesCondDistrib$rep)), "\r")) # the samples are filtered for the current replicate and the rep column is removed samplesCondDistrib_repI <- samplesCondDistrib[samplesCondDistrib$rep==i,c(-1)] # save as a csv file (it is not possible to give a data.frame as input to defineIndividualElement) write.csv(samplesCondDistrib_repI, "samples_repI.csv", row.names = F, quote=F) # define an element from an external file containing the samples from the conditional distribution defineIndividualElement(name="samplesCondDistrib", element="samples_repI.csv") # define output time grid defineOutputElement(name="outCc", list(data = data.frame(start = 0, interval = 0.1, final = 24), output = "Cc")) # use the defined element for the simulation setGroupElement(group = "simulationGroup1", elements = c("samplesCondDistrib","outCc")) # make sure the ids of each element are correctly mapped together setSharedIds(sharedIds = c("output", "treatment", "regressor", "individual")) # run the simulation runSimulation() # retrieve the results sim_repI <- getSimulationResults()$res$Cc # add replicate information sim_repI$rep <- i # merge all replicates sims <- rbind(sims, sim_repI) } ## Step 3: plotting To plot the prediction interval of each individual, we calculate the 5th, median and 95th percentile over the replicates for each individual and each time point. The prediction interval is then plotted with ggplot, stratified by id and the data is overlayed. R prctle <- sims %>% group_by(original_id,time) %>% summarise(P05=quantile(Cc,0.05), P50=quantile(Cc,0.5), P95=quantile(Cc,0.95)) %>% ungroup() %>% rename(id=original_id) ggplot(prctle) + geom_ribbon(aes(x=time,ymin=P05,ymax=P95), fill="#ff793f", alpha=0.5) + geom_line(aes(x=time,y=P50), color="#ff793f") + geom_point(data=obsData, aes(x=time,y=CONC)) + facet_wrap(.~id)+ ylab("Conentration")+ xlab('Time (hours)')+theme_bw() The resulting plot is the following with data truncated at **10 hours**. ![image-20260224-145709.png](https://monolixsuite.slp-software.com/__attachments/a_10306e674e319f61822d5e83911634cbca29a9fc727f818965cf257270e6ba7d/image-20260224-145709.png?cb=287ee75d3e4c26dbef4ee5ae1d3214aa) If the data is truncated at **4 hours** instead, we have less data and thus a higher uncertainty in the individual parameters. This leads to wider prediction intervals. For individuals which have no data information the elimination rate, the prediction of the elimination phase is mostly based on the population parameters which were used as a prior. ![image-20260224-145457.png](https://monolixsuite.slp-software.com/__attachments/a_a0719adfec05fe82130e1462925b0a8876ae31c84024650d23ea466e202d5ec6/image-20260224-145457.png?cb=f6cfc07ec19d375c7b3abc3fd995cabc) Last updated: February 24, 2026 --- version: "2024R1" language: "en" --- # buildAll ## Overview ### Description `buildAll` builds the complete statistical model by iteratively calling functions `buildmlx` and `buildVar`. #### Usage buildAll <- function(project=NULL, final.project=NULL, model="all", prior=NULL, weight=NULL, coef.w1=0.5, cv.min=0.001, fError.min=1e-3, paramToUse="all", covToTest="all", covToTransform="none", center.covariate=FALSE, criterion="BICc", linearization=FALSE, ll=T, test=T, direction=NULL, steps=1000, max.iter=20, explor.iter=2, seq.cov=FALSE, seq.corr=TRUE, seq.cov.iter=0, p.max=0.1, p.min=c(0.075, 0.05, 0.1), print=TRUE, nb.model=1, fix.param1=NULL, fix.param0=NULL, remove=T, add=T, delta=c(30,10,5), omega.set=NULL, pop.set1=NULL, pop.set2=NULL) #### Arguments **project** a string: the initial Monolix project **final.project** the final Monolix project (default is the original project)a string: the final Monolix project (default adds "_buildAll" to the original project) **model** the components of the model to optimize c("residualError", "covariate", "correlation"), (default="all") **prior** list of prior probabilities for each component of the model (default=NULL) **weight** list of penalty weights for each component of the model (default=NULL) **coef.w1** multiplicative weight coefficient used for the first iteration only (default=0.5) **cv.min** value of the coefficient of variation below which an individual parameter is considered fixed (default=0.001) **fError.min** minimum fraction of residual variance for combined error model (default = 1e-3) **paramToUse** list of parameters possibly function of covariates (default="all") **covToTest** components of the covariate model that can be modified (default="all") **covToTransform** list of (continuous) covariates to be log-transformed (default="none") **center.covariate** TRUE/{FALSE} center the covariates of the final model (default=FALSE) **criterion** penalization criterion to optimize c("AIC", "BIC", {"BICc"}, gamma) **linearization** TRUE/{FALSE} whether the computation of the likelihood is based on a linearization of the model (default=FALSE) **ll** {TRUE}/FALSE compute the observe likelihood and the criterion to optimize at each iteration (default=TRUE) **test** {TRUE}/FALSE perform additional statistical tests for building the model (default=TRUE) **direction** method for covariate search c({"full"}, "both", "backward", "forward"), (default="full" or "both") **steps** maximum number of iteration for stepAIC (default=1000) **max.iter** maximum number of SAMBA iterations (default=20) **exp.iter** number of iterations during the exploratory phase of SAMBA (default=2) **seq.cov** TRUE/{FALSE} whether the covariate model is built before the correlation model (default=FALSE) **seq.corr** {TRUE}/FALSE whether the correlation model is built iteratively (default=TRUE) **seq.cov.iter** number of iterations before building the correlation model (only when seq.cov=F, default=0) **p.max** maximum p-value used for removing non significant relationships between covariates and individual parameters (default=0.1) **p.min** vector of 3 minimum p-values used for testing the components of a new model (default=c(0.075, 0.05, 0.1)) **print** {TRUE}/FALSE display the results (default=TRUE) **nb.model** number of models to display at each iteration (default=1) **fix.param1** parameters with variability that cannot be removed (default=NULL) **fix.param0** parameters without variability that cannot be added (default=NULL) **remove** {TRUE}/FALSE try to remove random effects (default=TRUE) **add** {TRUE}/FALSE try to add random effects (default=TRUE) **delta** maximum difference in criteria for testing a new model (default=c(30,10,5)) **omega.set** settings to define how a variance varies during iterations of SAEM **pop.set1** Monolix settings 1 **pop.set2** Monolix settings 2 ## Example Using version \<= 4.0.2 of `Rsmlx` requires to load explicitly `lixoftConnectors`. This package is automatically loaded with `Rsmlx 4.0.3`. R library(lixoftConnectors) initializeLixoftConnectors(software="monolix") library(Rsmlx) We select a Monolix project R project <- "projects/simulatedPK2.mlxtran" Example using the defaults settings of `buildmlx` and `buildVar`: R buildAll.res1 <- buildAll(project) Selected model: R loadProject(buildAll.res1$project) getIndividualParameterModel() ## $name ## [1] "ka" "Cl" "V1" "Q" "V2" ## ## $distribution ## ka Cl V1 Q V2 ## "logNormal" "logNormal" "logNormal" "logNormal" "logNormal" ## ## $limits ## named list() ## ## $formula ## [1] "log(ka) = log(ka_pop) + eta_ka\nlog(Cl) = log(Cl_pop) + beta_Cl_X03*X03 + eta_Cl\nlog(V1) = log(V1_pop) + beta_V1_X01*X01 + beta_V1_X02*X02 + eta_V1\nlog(Q) = log(Q_pop) + beta_Q_X01*X01\nlog(V2) = log(V2_pop)\nCorrelations\n\tID : {Cl, V1}\n" ## ## $variability ## $variability$id ## ka Cl V1 Q V2 ## TRUE TRUE TRUE FALSE FALSE ## ## ## $covariateModel ## $covariateModel$ka ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## $covariateModel$Cl ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## $covariateModel$V1 ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## $covariateModel$Q ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## $covariateModel$V2 ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## ## $correlationBlocks ## $correlationBlocks$id ## $correlationBlocks$id[[1]] ## [1] "Cl" "V1" Example using several settings of `buildmlx` and `buildVar`: R buildAll.res2 <- buildAll(project, final.project="projects/buildAll2.mlxtran", paramToUse=c("Cl", "V1", "Q"), covToTest=c("X01", "X02", "X03", "X04"), fix.param0="Q", fix.param1=c("ka", "Cl")) Selected model: R loadProject(buildAll.res2$project) getIndividualParameterModel() ## $name ## [1] "ka" "Cl" "V1" "Q" "V2" ## ## $distribution ## ka Cl V1 Q V2 ## "logNormal" "logNormal" "logNormal" "logNormal" "logNormal" ## ## $limits ## named list() ## ## $formula ## [1] "log(ka) = log(ka_pop) + eta_ka\nlog(Cl) = log(Cl_pop) + beta_Cl_X03*X03 + eta_Cl\nlog(V1) = log(V1_pop) + beta_V1_X01*X01 + beta_V1_X02*X02 + eta_V1\nlog(Q) = log(Q_pop) + beta_Q_X01*X01\nlog(V2) = log(V2_pop)\nCorrelations\n\tID : {Cl, V1}\n" ## ## $variability ## $variability$id ## ka Cl V1 Q V2 ## TRUE TRUE TRUE FALSE FALSE ## ## ## $covariateModel ## $covariateModel$ka ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## $covariateModel$Cl ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## $covariateModel$V1 ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## $covariateModel$Q ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## $covariateModel$V2 ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## ## $correlationBlocks ## $correlationBlocks$id ## $correlationBlocks$id[[1]] ## [1] "Cl" "V1" Following what is done with `buildmlx`and `buildVar`, prior information about the statistical model can be introduced either by defining prior probabilities or introducing a weighting for the penalty term. In this example, relationships between covariates and parameters are strongly penalized (weight=5) while correlations are lightly penalized (weight=0.1). Different weights are used for the variances: variability of Q is privileged while that of V1 is strongly penalized. R buildAll.res3 <- buildAll(project, weight=list(covariate=5, correlation=0.3, variance=c(Q=0.05, V1=20))) R print(buildAll.res3$variability.model) ## NULL R loadProject(buildAll.res3$project) getIndividualParameterModel() ## $name ## [1] "ka" "Cl" "V1" "Q" "V2" ## ## $distribution ## ka Cl V1 Q V2 ## "logNormal" "logNormal" "logNormal" "logNormal" "logNormal" ## ## $limits ## named list() ## ## $formula ## [1] "log(ka) = log(ka_pop) + eta_ka\nlog(Cl) = log(Cl_pop) + beta_Cl_X03*X03 + eta_Cl\nlog(V1) = log(V1_pop) + beta_V1_X01*X01 + beta_V1_X02*X02 + eta_V1\nlog(Q) = log(Q_pop) + beta_Q_X01*X01\nlog(V2) = log(V2_pop)\nCorrelations\n\tID : {Cl, V1}\n" ## ## $variability ## $variability$id ## ka Cl V1 Q V2 ## TRUE TRUE TRUE FALSE FALSE ## ## ## $covariateModel ## $covariateModel$ka ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## $covariateModel$Cl ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## $covariateModel$V1 ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## $covariateModel$Q ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## $covariateModel$V2 ## X01 X02 X03 X04 X05 X06 X07 X08 X09 X10 ## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## ## ## $correlationBlocks ## $correlationBlocks$id ## $correlationBlocks$id[[1]] ## [1] "Cl" "V1" Last updated: July 15, 2025 --- version: "2024R1" language: "en" --- # buildVar ## Overview ### Description `buildVar` is designed to build the best variance model for the random effects by selecting which individual parameters vary and which ones are fixed. Penalization criterion can be either a custom penalization of the form γ\*(number of parameters), AIC (γ=2), BIC (γ=log(N)) or BICc. #### Usage buildVar <- function(project=NULL,final.project=NULL, prior=NULL, weight=NULL, cv.min=0.001, fix.param1=NULL, fix.param0=NULL, criterion="BICc", linearization=F, remove=T, add=T, delta=c(30,10,5), omega.set=NULL, pop.set1=NULL, pop.set2=NULL, print=TRUE) #### Arguments **project** a string: the initial Monolix project **final.project** the final Monolix project (default is the original project)a string: the final Monolix project (default adds "_var" to the original project) **prior** named vector of prior probabilities (default=NULL) **weight** named vector of weights (default=NULL) **cv.min** value of the coefficient of variation below which an individual parameter is considered fixed (default=0.001) **fix.param1** parameters with variability that cannot be removed (default=NULL) **fix.param0** parameters without variability that cannot be added (default=NULL) **criterion** penalization criterion to optimize c("AIC", "BIC", {"BICc"}, gamma) **linearization** TRUE/{FALSE} whether the computation of the likelihood is based on a linearization of the model (default=FALSE) **remove** {TRUE}/FALSE try to remove random effects (default=TRUE) **add** {TRUE}/FALSE try to add random effects (default=TRUE) **delta** maximum difference in criteria for testing a new model (default=c(30,10,5)) **omega.set** settings to define how a variance varies during iterations of SAEM **pop.set1** Monolix settings 1 **pop.set2** Monolix settings 2 **print** {TRUE}/FALSE display the results (default=TRUE) ## Example R library(Rsmlx) We select a Monolix project R project <- "projects/simulatedPK1.mlxtran" The model has 5 parameters: ka, Cl, V1, Q, V2. We will use `buildVar` to distinguish the parameters with and without IIV. R buildVar.res1 <- buildVar(project) ## ## -------------------------------------------------- ## ## Building the variance model ## ## __________________________________________________ ## ## Estimating the population parameters ## ## __________________________________________________ ## Iteration 1 ## ## removing variability... ## ## ----------------------- ## Step 1 ## Parameters without variability: ## Parameters with variability : ka Cl V1 Q V2 ## ## Criterion (linearization): 3843 ## trying to remove omega_V2 : 3839.2 ## trying to remove omega_Q : 3839.3 ## trying to remove omega_V1 : 3887.7 ## trying to remove omega_Cl : 5481.4 ## trying to remove omega_ka : 4011.8 ## ## Criterion: 3860.8 ## fitting the model with no variability on V2 : 3852.6 ## variability on V2 removed ## ## ----------------------- ## Step 2 ## Parameters without variability: V2 ## Parameters with variability : ka Cl V1 Q ## ## Criterion (linearization): 3837.8 ## trying to remove omega_Q : 3834.0 ## ## Criterion: 3852.6 ## fitting the model with no variability on V2 Q : 3849.2 ## variability on Q removed ## ## no more variability can be removed ## _______________________ ## ## adding variability... ## ## ----------------------- ## Step 1 ## Parameters without variability: Q V2 ## Parameters with variability : ka Cl V1 ## ## Criterion (linearization): 3833.9 ## trying to add omega_V2 : 3844.9 ## ## no more variability can be added ## ## __________________________________________________ ## ## Final variance model: ## ## Parameters without variability: Q V2 ## Parameters with variability : ka Cl V1 ## ## Fitting the final model using the original settings... ## ## Estimated criteria (importanceSampling): ## AIC BIC BICc s.e. ## 3795.10 3831.57 3848.96 0.21 ## ## total time: 263.4s R print(buildVar.res1) ## $project ## [1] "projects/simulatedPK1_var.mlxtran" ## ## $niter ## [1] 1 ## ## $change ## [1] TRUE ## ## $variability.model ## ka Cl V1 Q V2 ## TRUE TRUE TRUE FALSE FALSE ## ## $time ## elapsed ## 263.44 It is possible to fix the type of variability of some parameters of the model. In this example, we force Q to have no IIV while ka and Cl vary: R buildVar.res2 <- buildVar(project, final.project="projects/buildVar2.mlxtran", fix.param0="Q", fix.param1=c("ka", "Cl")) R print(buildVar.res2$variability.model) ## ka Cl V1 Q V2 ## TRUE TRUE TRUE FALSE FALSE Rather than forcing certain parameters to vary or to be fixed, one can introduce a priori information to favor the presence or absence of variability. We can either define a prior probability or introduce a weighting for the penalty term. For instance, weight=2 means that, for each parameter, we penalize twice as much the presence of a variance, while weight=0.5 means that the penalty is twice less than when the chosen criterion is used by default (i.e. weight=1). By default, BICc is used. Then, the criteria to minimize is -2 LL + weight x log(N) x (number of variances) In this example, add=F means that the algorithm starts with a full variance model and only tries to remove variances. R buildVar.res3 <- buildVar(project, final.project="projects/buildVar3.mlxtran", weight=3, add=F) R print(buildVar.res3$variability.model) ## ka Cl V1 Q V2 ## TRUE TRUE FALSE FALSE FALSE Different weights can be used for the different parameters. In this example, the variability of Q is privileged while that of V1 is strongly penalized R buildVar.res4 <- buildVar(project, final.project="projects/buildVar4.mlxtran", weight=c(Q=0.05, V1=20)) R print(buildVar.res4$variability.model) ## ka Cl V1 Q V2 ## TRUE TRUE FALSE TRUE FALSE Last updated: July 15, 2025 --- version: "2024R1" language: "en" --- # Changelog ## Version 2.0.1 ### Bugfix * Fix ACO `save_mode = "best"`not reliably saving the best model. In parallel mode, the "improved" flag from workers could become stale, causing a worse model to overwrite the saved best. The flag is now recomputed on the main process. Additionally, ACO now re-runs and saves the best model at the end of the search, matching the behavior of exhaustive search. ## Version 2.0.0 ### New Features #### Custom Model Libraries * **Custom model support** (`library = "custom"`): use your own model library instead of the built-in PK library. Provide a `model_creation_func(library, filters)` that returns a model file path or a Monolix `lib:` reference. Compatible with `"exhaustive_search"` and `"ACO"` algorithms. * **Monolix library integration** : `model_creation_func` can return `lib:` references from `getLibraryModelName()`, enabling search over any Monolix library (PKPD, PD, TMDD, TGI, etc.) without writing model files. * `param_mapping` argument: defines which parameters are available for IIV in each model variant via `common_params` and `filter_params`. #### IIV and Parameter Control * `fixed_iiv` argument: lock IIV on (`TRUE`) or off (`FALSE`) for specific parameters, reducing the search space. Works with all algorithms and with both `iiv = TRUE` and `iiv = FALSE`. * `param_distributions` argument: override the default log-normal distribution for specific parameters (supports `"normal"`, `"logNormal"`, `"logitNormal"` with optional limits). #### Parallel Execution * **Parallel model evaluation** for ACO and exhaustive search using `future` / `future.apply`. Workers run models concurrently across local cores or HPC clusters. #### Other * Specify `table_func` to output additional information about the model. #### Search and Algorithm Improvements * `settings$tasks`: configure which Monolix scenario tasks to run (e.g. disable `standardErrorEstimation` for faster exploration). * `settings$table_func`: custom function to output additional information about each run. Last updated: March 11, 2026 --- version: "2024R1" language: "en" --- # Comparison to popED ## Introduction This page shows a comparison of the results of mlxDesignEval and popED. To do so, we re-use the examples presented in the [popED documentation](https://andrewhooker.github.io/PopED/articles/examples.html). popED calculates the RSE on the variance of the random effects and error model, while mlxDesignEval works with standard deviations by default. To request RSEs on the variances with mlxDesignEval, we will use the option `rse_on_variance=T`. For each example, the RSEs are displayed side by side and the percentage difference is calculated. If the difference is less than 1%, we consider that the results are the same (i.e `is_identical=T`). Note that the OFV and FIM will differ between mlxDesignEval and popED and cannot be compared directly. **The RSE are identical between popED and mlxDesignEval on all examples.** In order to allow an easy comparison of the models syntax, we will use the `inlineModel()` function to create model files. The results would be the same using a Monolix or Simulx project (containing the same model). ## Ex1a: PK,1-comp,oral,MD This example uses the analytical solution of a 1-compartment model with first-order absorption. The bioavailability is set as fixed as it cannot be estimated, and the error model parameters are also fixed. There are two groups with different dose amounts. ### popED ##-- Model: One comp first order absorption ## -- Analytic solution for both mutiple and single dosing ff <- function(model_switch,xt,parameters,poped.db){ with(as.list(parameters),{ y=xt N = floor(xt/TAU)+1 y=(DOSE*Favail/V)*(KA/(KA - CL/V)) * (exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) return(list( y=y,poped.db=poped.db)) }) } ## -- parameter definition function ## -- names match parameters in function ff sfg <- function(x,a,bpop,b,bocc){ parameters=c( V=bpop[1]*exp(b[1]), KA=bpop[2]*exp(b[2]), CL=bpop[3]*exp(b[3]), Favail=bpop[4], DOSE=a[1], TAU=a[2]) return( parameters ) } ## -- Residual unexplained variablity (RUV) function ## -- Additive + Proportional feps <- function(model_switch,xt,parameters,epsi,poped.db){ returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db)) y <- returnArgs[[1]] poped.db <- returnArgs[[2]] y = y*(1+epsi[,1])+epsi[,2] return(list( y= y,poped.db =poped.db )) } ## -- Define design and design space poped.db <- create.poped.database(ff_fun="ff", fg_fun="sfg", fError_fun="feps", bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9), notfixed_bpop=c(1,1,1,0), d=c(V=0.09,KA=0.09,CL=0.25^2), sigma=c(0.04,5e-6), notfixed_sigma=c(0,0), m=2, groupsize=20, xt=c( 1,2,8,240,245), minxt=c(0,0,0,240,240), maxxt=c(10,10,10,248,248), bUseGrouped_xt=1, a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)), maxa=c(DOSE=200,TAU=24), mina=c(DOSE=0,TAU=24)) res1a_popED <- PopED::evaluate_design(poped.db) ### mlxDesignEval With mlxDesignEval: * DOSE and TAU do not need to be specified in the model * the model is defined using the `pkmodel()` macro, which will be replaced by the corresponding analytical solution in the background * the bioavailability `F` has no random effects, it can be set as either `normal` or `logNormal` (same results) * `y = y*(1+epsi[,1])+epsi[,2]` in popED corresponds to a `combined2` error model in mlxtran language model1a <- inlineModel(" [INDIVIDUAL] input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka} DEFINITION: Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl} F = {distribution=normal, typical=F_pop, no-variability} V = {distribution=logNormal, typical=V_pop, sd=omega_V} ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka} [LONGITUDINAL] input = {a, b} input = {F, ka, V, Cl} PK: Cc = pkmodel(ka, V, Cl, p=F) OUTPUT: output = {Cc} DEFINITION: DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}") When defining the `population_parameters`, remember that omega parameters represent a standard deviation while the popED parameter are variances. We thus use the square root of the value given in popED. treatment1<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=20)) treatment2<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40)) g1 <- list(size=20, treatment=treatment1) g2 <- list(size=20, treatment=treatment2) res1a_mlx <- mlxDesignEval::evaluate_design( model_file=model1a, population_parameters = data.frame(F_pop=0.9, ka_pop=0.25, V_pop=72.8, Cl_pop=3.75, omega_ka=sqrt(0.09), omega_V=sqrt(0.09), omega_Cl=0.25, a=sqrt(5e-6), b=sqrt(0.04)), fixed_parameters = c("F_pop","a","b"), group=list(g1,g2), output=list(output="DV", data=data.frame(time=c(1,2,8,240,245))), rse_on_variance = T) ### Comparison #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 V 8.215338 V_pop 8.215339 1.661207e-05 TRUE #> 2 KA 10.090955 ka_pop 10.090937 1.835064e-04 TRUE #> 3 CL 4.400304 Cl_pop 4.400307 6.276719e-05 TRUE #> 4 d_V 39.833230 var_omega_V 39.833208 5.578078e-05 TRUE #> 5 d_KA 60.089601 var_omega_ka 60.089600 1.920909e-06 TRUE #> 6 d_CL 27.391518 var_omega_Cl 27.391561 1.557890e-04 TRUE ## Ex1b: PK,1-comp,oral,MD Same example as before but parametrized with elimination rate instead of clearance. ### popED ## -- names match parameters in function defined in ff_file fg.PK.1.comp.oral.md.param.2 <- function(x,a,bpop,b,bocc){ ## -- parameter definition function parameters=c( V=bpop[1]*exp(b[1]), KA=bpop[2]*exp(b[2]), KE=bpop[3]*exp(b[3]), Favail=bpop[4], DOSE=a[1], TAU=a[2]) return( parameters ) } ## -- Define design and design space poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.md.KE", fg_file="fg.PK.1.comp.oral.md.param.2", fError_file="feps.add.prop", groupsize=20, m=2, sigma=c(0.04,5e-6), bpop=c(V=72.8,KA=0.25,KE=3.75/72.8,Favail=0.9), d=c(V=0.09,KA=0.09,KE=0.25^2), notfixed_bpop=c(1,1,1,0), notfixed_sigma=c(0,0), xt=c( 1,2,8,240,245), minxt=c(0,0,0,240,240), maxxt=c(10,10,10,248,248), bUseGrouped_xt=1, a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)), maxa=c(DOSE=200,TAU=40), mina=c(DOSE=0,TAU=2)) ## evaluate initial design res1b_popED <- PopED::evaluate_design(poped.db) ### mlxDesignEval model1b <- inlineModel(" [INDIVIDUAL] input = {k_pop, omega_k, F_pop, V_pop, omega_V, ka_pop, omega_ka} DEFINITION: k = {distribution=logNormal, typical=k_pop, sd=omega_k} F = {distribution=normal, typical=F_pop, no-variability} V = {distribution=logNormal, typical=V_pop, sd=omega_V} ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka} [LONGITUDINAL] input = {a, b} input = {F, ka, V, k} PK: Cc = pkmodel(ka, V, k, p=F) OUTPUT: output = {Cc} DEFINITION: DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}") treatment1<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=20)) treatment2<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40)) g1 <- list(size=20, treatment=treatment1) g2 <- list(size=20, treatment=treatment2) res1b_mlx <- mlxDesignEval::evaluate_design( model_file = model1b, population_parameters = data.frame(F_pop=0.9, ka_pop=0.25, V_pop=72.8, k_pop=3.75/72.8, omega_ka=sqrt(0.09), omega_V=sqrt(0.09), omega_k=0.25, a=sqrt(5e-6), b=sqrt(0.04)), fixed_parameters = c("F_pop","a","b"), group = list(g1,g2), output = list(output="DV", data=data.frame(time=c( 1,2,8,240,245))), rse_on_variance = T) ### Comparison #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 V 8.215338 V_pop 8.215339 1.980025e-05 TRUE #> 2 KA 10.090955 ka_pop 10.090937 1.838245e-04 TRUE #> 3 KE 7.566975 k_pop 7.566974 4.300219e-06 TRUE #> 4 d_V 31.220520 var_omega_V 31.220506 4.337963e-05 TRUE #> 5 d_KA 44.677836 var_omega_ka 44.677859 5.190823e-05 TRUE #> 6 d_KE 38.005067 var_omega_k 38.005052 3.844144e-05 TRUE ## Ex1c: PK,1-comp,oral,MD Same example as 1a but written with an ODE system instead of an analytical solution. ### popED The ODE system can be solved in R with deSolve (slow) or in C++ using Rcpp (faster). ## define the ODE PK.1.comp.oral.ode <- function(Time, State, Pars){ with(as.list(c(State, Pars)), { dA1 <- -KA*A1 dA2 <- KA*A1 - (CL/V)*A2 return(list(c(dA1, dA2))) }) } ## define the initial conditions and the dosing ff.ode <- function(model_switch, xt, parameters, poped.db){ with(as.list(parameters),{ A_ini <- c(A1=0, A2=0) times_xt <- drop(xt) #xt[,,drop=T] dose_times = seq(from=0,to=max(times_xt),by=TAU) eventdat <- data.frame(var = c("A1"), time = dose_times, value = c(DOSE*Favail), method = c("add")) times <- sort(c(times_xt,dose_times)) out <- ode(A_ini, times, PK.1.comp.oral.ode, parameters, events = list(data = eventdat))#atol=1e-13,rtol=1e-13) y = out[, "A2"]/(V) y=y[match(times_xt,out[,"time"])] y=cbind(y) return(list(y=y,poped.db=poped.db)) }) } ## -- parameter definition function ## -- names match parameters in function ff sfg <- function(x,a,bpop,b,bocc){ parameters=c( V=bpop[1]*exp(b[1]), KA=bpop[2]*exp(b[2]), CL=bpop[3]*exp(b[3]), Favail=bpop[4], DOSE=a[1], TAU=a[2]) return( parameters ) } poped.db <- create.poped.database(ff_fun=ff.ode, fError_fun=feps.add.prop, fg_fun=sfg, groupsize=20, m=2, #number of groups sigma=c(0.04,5e-6), bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9), d=c(V=0.09,KA=0.09,CL=0.25^2), notfixed_bpop=c(1,1,1,0), notfixed_sigma=c(0,0), xt=c( 1,2,8,240,245), minxt=c(0,0,0,240,240), maxxt=c(10,10,10,248,248), discrete_xt = list(0:248), bUseGrouped_xt=1, a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)), maxa=c(DOSE=200,TAU=24), mina=c(DOSE=0,TAU=24)) # calculations are noticeably slower than with the analytic solution res1c_popED_Rode <- PopED::evaluate_design(poped.db) # using Rcpp cppFunction('List one_comp_oral_ode(double Time, NumericVector A, NumericVector Pars) { int n = A.size(); NumericVector dA(n); double CL = Pars[0]; double V = Pars[1]; double KA = Pars[2]; dA[0] = -KA*A[0]; dA[1] = KA*A[0] - (CL/V)*A[1]; return List::create(dA); }') ff.ode.rcpp <- function(model_switch, xt, parameters, poped.db){ with(as.list(parameters),{ A_ini <- c(A1=0, A2=0) times_xt <- drop(xt) #xt[,,drop=T] dose_times = seq(from=0,to=max(times_xt),by=TAU) eventdat <- data.frame(var = c("A1"), time = dose_times, value = c(DOSE*Favail), method = c("add")) times <- sort(c(times_xt,dose_times)) out <- ode(A_ini, times, one_comp_oral_ode, c(CL,V,KA), events = list(data = eventdat)) y = out[, "A2"]/(V) y=y[match(times_xt,out[,"time"])] y=cbind(y) return(list(y=y,poped.db=poped.db)) }) } ## -- Update poped.db with compiled function poped.db.compiled.rcpp <- create.poped.database(poped.db,ff_fun=ff.ode.rcpp) # calculations are much faster than with the pure R solution but slower than .dll compiled solution in desolve res1c_popED_Rcpp <- PopED::evaluate_design(poped.db.compiled.rcpp) ### mlxDesignEval The model in mlxtran language is always and automatically converted to C++ code to run fast. There is nothing to do on the user side. The solver for stiff ODEs has a better precision and is usually preferred. It can be set with `odeType=stiff` in the structural model. model1c <- inlineModel(" [INDIVIDUAL] input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka} DEFINITION: Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl} F = {distribution=logNormal, typical=F_pop, no-variability} V = {distribution=logNormal, typical=V_pop, sd=omega_V} ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka} [LONGITUDINAL] input = {a, b} input = {F, ka, V, Cl} PK: depot(target=Ad) EQUATION: odeType=stiff ddt_Ad = -ka*Ad ddt_Ac = ka*Ad - (Cl/V)*Ac Cc = Ac/V OUTPUT: output = {Cc} DEFINITION: DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}") treatment1<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=20)) treatment2<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40)) g1 <- list(size=20, treatment=treatment1) g2 <- list(size=20, treatment=treatment2) res1c_mlx <- mlxDesignEval::evaluate_design( model_file = model1c, population_parameters = data.frame(F_pop=0.9, ka_pop=0.25, V_pop=72.8, Cl_pop=3.75, omega_ka=sqrt(0.09), omega_V=sqrt(0.09), omega_Cl=0.25, a=sqrt(5e-6), b=sqrt(0.04)), fixed_parameters = c("F_pop","a","b"), group = list(g1,g2), output = list(output="DV", data=data.frame(time=c( 1,2,8,240,245))), rse_on_variance = T) ### Comparison Comparison to solution with deSolve: #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 V 8.215334 V_pop 8.212722 0.031791907 TRUE #> 2 KA 10.090963 ka_pop 10.084245 0.066577753 TRUE #> 3 CL 4.400304 Cl_pop 4.400000 0.006913958 TRUE #> 4 d_V 39.833193 var_omega_V 39.821267 0.029938234 TRUE #> 5 d_KA 60.089706 var_omega_ka 60.030178 0.099065990 TRUE #> 6 d_CL 27.391516 var_omega_Cl 27.387705 0.013911217 TRUE Comparison to solution with Rcpp: #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 V 8.215334 V_pop 8.212722 0.031791907 TRUE #> 2 KA 10.090963 ka_pop 10.084245 0.066577753 TRUE #> 3 CL 4.400304 Cl_pop 4.400000 0.006913958 TRUE #> 4 d_V 39.833193 var_omega_V 39.821267 0.029938234 TRUE #> 5 d_KA 60.089706 var_omega_ka 60.030178 0.099065990 TRUE #> 6 d_CL 27.391516 var_omega_Cl 27.387705 0.013911217 TRUE ## Ex2: warfarin This example uses the Warfarin model used in the design evaluation software comparison [Nyberg et al., "Methods and software tools for design evaluation for population pharmacokinetics-pharmacodynamics studies", Br. J. Clin. Pharm., 2014.](https://bpspubs.onlinelibrary.wiley.com/doi/10.1111/bcp.12352). It uses the analytical solution of a 1-compartment model with first-order absorption with a single dose. ### popED ff <- function(model_switch,xt,parameters,poped.db){ ##-- Model: One comp first order absorption with(as.list(parameters),{ y=xt y=(DOSE*Favail*KA/(V*(KA-CL/V)))*(exp(-CL/V*xt)-exp(-KA*xt)) return(list(y=y,poped.db=poped.db)) }) } sfg <- function(x,a,bpop,b,bocc){ ## -- parameter definition function parameters=c(CL=bpop[1]*exp(b[1]), V=bpop[2]*exp(b[2]), KA=bpop[3]*exp(b[3]), Favail=bpop[4], DOSE=a[1]) return(parameters) } feps <- function(model_switch,xt,parameters,epsi,poped.db){ ## -- Residual Error function ## -- Proportional returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db)) y <- returnArgs[[1]] poped.db <- returnArgs[[2]] y = y*(1+epsi[,1]) return(list(y=y,poped.db=poped.db)) } ## -- Define initial design and design space poped.db <- create.poped.database(ff_file="ff", fg_file="sfg", fError_file="feps", bpop=c(CL=0.15, V=8, KA=1.0, Favail=1), notfixed_bpop=c(1,1,1,0), d=c(CL=0.07, V=0.02, KA=0.6), sigma=0.01, groupsize=32, xt=c( 0.5,1,2,6,24,36,72,120), minxt=0, maxxt=120, a=70) ## evaluate initial design res2a_popED <- PopED::evaluate_design(poped.db) ### mlxDesignEval The model is the same as for the example 1, except that a proportional error model is used. model2 <- inlineModel(" [INDIVIDUAL] input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka} DEFINITION: Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl} F = {distribution=logNormal, typical=F_pop, no-variability} V = {distribution=logNormal, typical=V_pop, sd=omega_V} ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka} [LONGITUDINAL] input = {b} input = {F, ka, V, Cl} PK: Cc = pkmodel(ka, V, Cl, p=F) OUTPUT: output = {Cc} DEFINITION: DV = {distribution=normal, prediction=Cc, errorModel=proportional(b)}") res2a_mlx <- mlxDesignEval::evaluate_design( model_file = model2, population_parameters = data.frame(F_pop=1, ka_pop=1, V_pop=8, Cl_pop=0.15, omega_ka=sqrt(0.6), omega_V=sqrt(0.02), omega_Cl=sqrt(0.07), b=sqrt(0.01)), fixed_parameters = c("F_pop"), group = list(size=32), treatment = list(data=data.frame(time=0, amount=70)), output =list(output="DV", data=data.frame(time=c(0.5,1,2,6,24,36,72,120))), rse_on_variance = T) ### Comparison #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 CL 4.738266 Cl_pop 4.738270 6.837835e-05 TRUE #> 2 V 2.756206 V_pop 2.756208 4.766023e-05 TRUE #> 3 KA 13.925829 ka_pop 13.925839 7.575014e-05 TRUE #> 4 d_CL 25.627205 var_omega_Cl 25.627221 6.124399e-05 TRUE #> 5 d_V 30.344316 var_omega_V 30.344315 3.414725e-06 TRUE #> 6 d_KA 25.777327 var_omega_ka 25.777328 3.518511e-06 TRUE #> 7 SIGMA[1,1] 11.170784 var_b 11.170783 9.912263e-06 TRUE ## Ex3: PKPD,1-comp,oral,MD This example uses a model with two outputs (PK and PD), both using an analytical solution. The PK model has 1 compartment and first order absorption. The PD is an Imax model. ### popED ##-- Model: One comp first order absorption + inhibitory imax ## -- works for both mutiple and single dosing ff <- function(model_switch,xt,parameters,poped.db){ with(as.list(parameters),{ y=xt MS <- model_switch # PK model N = floor(xt/TAU)+1 CONC=(DOSE*Favail/V)*(KA/(KA - CL/V)) * (exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) # PD model EFF = E0*(1 - CONC*IMAX/(IC50 + CONC)) y[MS==1] = CONC[MS==1] y[MS==2] = EFF[MS==2] return(list( y= y,poped.db=poped.db)) }) } ## -- parameter definition function sfg <- function(x,a,bpop,b,bocc){ parameters=c( V=bpop[1]*exp(b[1]), KA=bpop[2]*exp(b[2]), CL=bpop[3]*exp(b[3]), Favail=bpop[4], DOSE=a[1], TAU = a[2], E0=bpop[5]*exp(b[4]), IMAX=bpop[6], IC50=bpop[7]) return( parameters ) } ## -- Residual Error function feps <- function(model_switch,xt,parameters,epsi,poped.db){ returnArgs <- ff(model_switch,xt,parameters,poped.db) y <- returnArgs[[1]] poped.db <- returnArgs[[2]] MS <- model_switch pk.dv <- y*(1+epsi[,1])+epsi[,2] pd.dv <- y*(1+epsi[,3])+epsi[,4] y[MS==1] = pk.dv[MS==1] y[MS==2] = pd.dv[MS==2] return(list( y= y,poped.db =poped.db )) } poped.db <- create.poped.database(ff_fun="ff", fError_fun="feps", fg_fun="sfg", groupsize=20, m=3, bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,E0=1120,IMAX=0.807,IC50=0.0993), notfixed_bpop=c(1,1,1,0,1,1,1), d=c(V=0.09,KA=0.09,CL=0.25^2,E0=0.09), sigma=c(0.04,5e-6,0.09,100), notfixed_sigma=c(0,0,0,0), xt=c( 1,2,8,240,240,1,2,8,240,240), minxt=c(0,0,0,240,240,0,0,0,240,240), maxxt=c(10,10,10,248,248,10,10,10,248,248), discrete_xt = list(0:248), G_xt=c(1,2,3,4,5,1,2,3,4,5), bUseGrouped_xt=1, model_switch=c(1,1,1,1,1,2,2,2,2,2), a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24),c(DOSE=0, TAU=24)), maxa=c(DOSE=200,TAU=40), mina=c(DOSE=0,TAU=2), ourzero=0) ## evaluate initial design res3a_popED <- PopED::evaluate_design(poped.db) ### mlxDesignEval model3 <- inlineModel(" [INDIVIDUAL] input = {V_pop, omega_V, KA_pop, omega_KA, CL_pop, omega_CL, Favail_pop, E0_pop, omega_E0, IMAX_pop, IC50_pop} DEFINITION: V = {distribution=logNormal, typical=V_pop, sd=omega_V} KA = {distribution=logNormal, typical=KA_pop, sd=omega_KA} CL = {distribution=logNormal, typical=CL_pop, sd=omega_CL} Favail = {distribution=logNormal, typical=Favail_pop, no-variability} E0 = {distribution=logNormal, typical=E0_pop, sd=omega_E0} IMAX = {distribution=logNormal, typical=IMAX_pop, no-variability} IC50 = {distribution=logNormal, typical=IC50_pop, no-variability} [LONGITUDINAL] input = {aCONC, bCONC, aEFF, bEFF} input = {V,KA,CL,Favail,E0,IMAX,IC50} EQUATION: Cc = pkmodel(ka=KA,V,Cl=CL,p=Favail) Eff = E0*(1 - Cc*IMAX/(IC50 + Cc)) OUTPUT: output = {Cc,Eff} DEFINITION: yCONC = {distribution=normal, prediction=Cc, errorModel=combined2(aCONC, bCONC)} yEFF = {distribution=normal, prediction=Eff, errorModel=combined2(aEFF, bEFF)}") The popED example defines two PK measurements at the same time t=240. If the same is done with mlxDesignEval, the two times are merged into a single one, and the information that two measurements are done at the same time is lost. To capture the double measurements, we give two slightly different times: 240 and 240.001. treatment1 <- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=20)) treatment2 <- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40)) treatment3 <- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=0)) g1 <- list(size=20, treatment=treatment1) g2 <- list(size=20, treatment=treatment2) g3 <- list(size=20, treatment=treatment3) # it is not possible to give twice the same time (will be merged into one by simulx) so giving slightly different value outPK <- list(output="yCONC", data=data.frame(time=c(1,2,8,240,240.001))) outPD <- list(output="yEFF", data=data.frame(time=c(1,2,8,240,240.001))) res3a_mlx_option1 <- mlxDesignEval::evaluate_design( model_file = model3, population_parameters = data.frame(Favail_pop=0.9, KA_pop=0.25, V_pop=72.8, CL_pop=3.75, E0_pop=1120,IMAX_pop=0.807,IC50_pop=0.0993, omega_KA=sqrt(0.09), omega_V=sqrt(0.09), omega_CL=sqrt(0.25^2),omega_E0=sqrt(0.09), bCONC=sqrt(0.04),aCONC=sqrt(5e-6),bEFF=sqrt(0.09),aEFF=sqrt(100)), fixed_parameters = c("Favail_pop","bCONC","aCONC","bEFF","aEFF"), group = list(g1,g2,g3), output = list(outPK,outPD), rse_on_variance = T) ### Comparison #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 V 8.119842 V_pop 8.119667 2.153595e-03 TRUE #> 2 KA 9.968612 KA_pop 9.968355 2.580165e-03 TRUE #> 3 CL 4.304635 CL_pop 4.304765 3.011526e-03 TRUE #> 4 E0 7.076883 E0_pop 7.076857 3.718358e-04 TRUE #> 5 IMAX 9.895340 IMAX_pop 9.894347 1.003892e-02 TRUE #> 6 IC50 39.478269 IC50_pop 39.477335 2.365973e-03 TRUE #> 7 d_V 38.960998 var_omega_V 38.960436 1.443299e-03 TRUE #> 8 d_KA 58.523188 var_omega_KA 58.521487 2.905635e-03 TRUE #> 9 d_CL 25.832775 var_omega_CL 25.833807 3.996265e-03 TRUE #> 10 d_E0 22.036110 var_omega_E0 22.036132 9.915972e-05 TRUE ## Ex4: PKPD,1-comp,Emax This example uses a model with two outputs (PK and PD), both using an analytical solution. The PK model has 1 compartment and bolus administration. The PD is an Emax model. ### popED ff <- function(model_switch,xt,parameters,poped.db){ with(as.list(parameters),{ y=xt MS <- model_switch # PK model CONC = DOSE/V*exp(-CL/V*xt) # PD model EFF = E0 + CONC*EMAX/(EC50 + CONC) y[MS==1] = CONC[MS==1] y[MS==2] = EFF[MS==2] return(list( y= y,poped.db=poped.db)) }) } ## -- parameter definition function sfg <- function(x,a,bpop,b,bocc){ parameters=c( CL=bpop[1]*exp(b[1]) , V=bpop[2]*exp(b[2]) , E0=bpop[3]*exp(b[3]) , EMAX=bpop[4], EC50=bpop[5]*exp(b[4]) , DOSE=a[1] ) return( parameters ) } ## -- Residual Error function ## -- Proportional PK + additive PD feps <- function(model_switch,xt,parameters,epsi,poped.db){ returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db)) y <- returnArgs[[1]] poped.db <- returnArgs[[2]] MS <- model_switch prop.err <- y*(1+epsi[,1]) add.err <- y+epsi[,2] y[MS==1] = prop.err[MS==1] y[MS==2] = add.err[MS==2] return(list( y= y,poped.db =poped.db )) } poped.db <- create.poped.database(ff_fun=ff, fError_fun=feps, fg_fun=sfg, groupsize=20, m=3, sigma=diag(c(0.15,0.015)), bpop=c(CL=0.5,V=0.2,E0=1,EMAX=1,EC50=1), d=c(CL=0.09,V=0.09,E0=0.04,EC50=0.09), xt=c( 0.33,0.66,0.9,5,0.1,1,2,5), bUseGrouped_xt=1, model_switch=c( 1,1,1,1,2,2,2,2), minxt=0, maxxt=5, ourzero = 0, a=list(c(DOSE=0),c(DOSE=1),c(DOSE=2)), maxa=c(DOSE=10), mina=c(DOSE=0)) ## evaluate initial design res4_popED <- PopED::evaluate_design(poped.db) #> Problems inverting the matrix. Results could be misleading. ### mlxDesignEval model4 <- inlineModel(" [INDIVIDUAL] input = {V_pop, omega_V, CL_pop, omega_CL, E0_pop, omega_E0, EC50_pop, omega_EC50, EMAX_pop} DEFINITION: V = {distribution=logNormal, typical=V_pop, sd=omega_V} CL = {distribution=logNormal, typical=CL_pop, sd=omega_CL} E0 = {distribution=logNormal, typical=E0_pop, sd=omega_E0} EC50 = {distribution=logNormal, typical=EC50_pop, sd=omega_EC50} EMAX = {distribution=logNormal, typical=EMAX_pop, no-variability} [LONGITUDINAL] input = {bCONC, aEFF} input = {V,CL,E0,EMAX,EC50} EQUATION: Cc = pkmodel(V=V,Cl=CL) Eff = E0 + Cc*EMAX/(EC50 + Cc) OUTPUT: output = {Cc,Eff} DEFINITION: yCONC = {distribution=normal, prediction=Cc, errorModel=proportional(bCONC)} yEFF = {distribution=normal, prediction=Eff, errorModel=constant(aEFF)}") treatment1<- list(data=data.frame(time=0, amount=0)) treatment2<- list(data=data.frame(time=0, amount=1)) treatment3<- list(data=data.frame(time=0, amount=2)) g1 <- list(size=20, treatment=treatment1) g2 <- list(size=20, treatment=treatment2) g3 <- list(size=20, treatment=treatment3) outPK <- list(output="yCONC", data=data.frame(time=c(0.33,0.66,0.9,5))) outPD <- list(output="yEFF", data=data.frame(time=c(0.1,1,2,5))) res4_mlx <- mlxDesignEval::evaluate_design( model_file = model4, population_parameters = data.frame(CL_pop=0.5,V_pop=0.2,E0_pop=1,EMAX_pop=1,EC50_pop=1, omega_V=sqrt(0.09), omega_CL=sqrt(0.09),omega_E0=sqrt(0.04),omega_EC50=sqrt(0.09), bCONC=sqrt(0.15),aEFF=sqrt(0.015)), group = list(g1,g2,g3), output = list(outPK,outPD), rse_on_variance = T) ### Comparison #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 CL 6.032896 CL_pop 6.032893 4.680235e-05 TRUE #> 2 V 6.288412 V_pop 6.288410 2.768071e-05 TRUE #> 3 E0 2.802303 E0_pop 2.802303 1.030050e-05 TRUE #> 4 EMAX 3.555580 EMAX_pop 3.555587 1.938459e-04 TRUE #> 5 EC50 14.508492 EC50_pop 14.508480 7.788633e-05 TRUE #> 6 d_CL 30.322296 var_omega_CL 30.322269 8.930289e-05 TRUE #> 7 d_V 32.964398 var_omega_V 32.964358 1.224145e-04 TRUE #> 8 d_E0 20.298320 var_omega_E0 20.298325 2.401653e-05 TRUE #> 9 d_EC50 142.920289 var_omega_EC50 142.920335 3.266862e-05 TRUE #> 10 SIGMA[1,1] 14.728239 var_bCONC 14.728240 1.091985e-05 TRUE #> 11 SIGMA[2,2] 11.723581 var_aEFF 11.723583 2.183795e-05 TRUE ## Ex5: PD,Emax-Hill This example uses a model with PD only, with an analytical solution of an Emax model with sigmoidicity depending on the dose. ### popED Note that the time `xt` is used as a surrogate for the dose (`DOSE=xt`). ff.emax.hill <- function(model_switch,xt,parameters,poped.db){ with(as.list(parameters),{ y=xt DOSE = xt y=BASE + EMAX*DOSE^(GAMMA)/(ED50^(GAMMA) + DOSE^(GAMMA)) return(list( y= y,poped.db=poped.db)) }) } ## -- parameter definition function sfg.emax.hill <- function(x,a,bpop,b,bocc){ parameters=c( EMAX=bpop[1]*exp(b[1]), ED50=bpop[2]*exp(b[2]), GAMMA=bpop[3], BASE=bpop[4]+b[3]) return( parameters ) } poped.db <- create.poped.database(ff_fun=ff.emax.hill, fError_fun=feps.add.prop, fg_fun=sfg.emax.hill, groupsize=100, m=1, bpop=c(EMAX=100,ED50=20,GAMMA=4.5,BASE=1), d=c(EMAX=0.0625,ED50=0.0625,BASE=0.0625), sigma=diag(c(0.01,.1)), xt=seq(0,50,length.out=8), minxt=0, maxxt=50, ourzero=0) ## evaluate initial design res5_popED <- PopED::evaluate_design(poped.db) ### mlxDesignEval `Emax` and `EC50` have a lognormal distribution while `Base` has a normal distribution. `Gamma` has no variability and can be set as either logNormal or normal. In order to consider the dose on the x axis instead of time, we rename the internal time variable `t` as `DOSE` (in the same way as done with popED). model5 <- inlineModel(" [INDIVIDUAL] input = {EMAX_pop, BASE_pop, omega_BASE, ED50_pop, omega_ED50, GAMMA_pop, omega_EMAX} DEFINITION: EMAX = {distribution=logNormal, typical=EMAX_pop, sd=omega_EMAX} BASE = {distribution=normal, typical=BASE_pop, sd=omega_BASE} ED50 = {distribution=logNormal, typical=ED50_pop, sd=omega_ED50} GAMMA = {distribution=normal, typical=GAMMA_pop, no-variability} [LONGITUDINAL] input = {aCONC, bCONC} input = {BASE,EMAX,ED50,GAMMA} EQUATION: DOSE=t Pred=BASE + EMAX*DOSE^(GAMMA)/(ED50^(GAMMA) + DOSE^(GAMMA)) OUTPUT: output = {Pred} DEFINITION: yDV = {distribution=normal, prediction=Pred, errorModel=combined2(aCONC, bCONC)}") The different doses tested are defined in the `output` element, which data frame column should stay `time` (even if it now represents the dose). out <- list(output="yDV", data=data.frame(time=seq(0,50,length.out=8))) res5_mlx <- mlxDesignEval::evaluate_design( model_file = model5, population_parameters = data.frame(EMAX_pop=100,ED50_pop=20,GAMMA_pop=4.5,BASE_pop=1, omega_EMAX=sqrt(0.0625),omega_ED50=sqrt(0.0625),omega_BASE=sqrt(0.0625), bCONC=sqrt(0.01),aCONC=sqrt(0.1)), group = list(size=100), output = out, rse_on_variance = T) ### Comparison #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 EMAX 2.588804 EMAX_pop 2.588805 1.742681e-05 TRUE #> 2 ED50 2.554668 ED50_pop 2.554668 9.715213e-06 TRUE #> 3 GAMMA 1.112444 GAMMA_pop 1.112444 1.946455e-05 TRUE #> 4 BASE 3.922208 BASE_pop 3.922206 5.208619e-05 TRUE #> 5 d_EMAX 14.825482 var_omega_EMAX 14.825492 6.249867e-05 TRUE #> 6 d_ED50 14.375072 var_omega_ED50 14.375067 3.504967e-05 TRUE #> 7 d_BASE 33.257104 var_omega_BASE 33.257108 1.172602e-05 TRUE #> 8 SIGMA[1,1] 7.072694 var_bCONC 7.072680 1.946146e-04 TRUE #> 9 SIGMA[2,2] 18.487445 var_aCONC 18.487444 1.021925e-06 TRUE ## Ex6: PK,1-comp,oral,SD This example uses a simple PK model with 1 compartment and first-order absorption, receiving a single dose. The ODE system solution and analytical solution are compared. ### popED ##-- Model: One comp first order absorption, analytic solution PK.1.comp.oral.sd.ff <- function(model_switch,xt,parameters,poped.db){ with(as.list(parameters),{ y=xt y=(DOSE*Favail*KA/(V*(KA-KE)))*(exp(-KE*xt)-exp(-KA*xt)) return(list( y= y,poped.db=poped.db)) }) } # ODE solution PK.1.comp.oral.ode <- function(Time, State, Pars){ with(as.list(c(State, Pars)), { dA1 <- -KA*A1 dA2 <- KA*A1 - KE*A2 return(list(c(dA1, dA2))) }) } PK.1.comp.oral.sd.ff.ode <- function(model_switch,xt,parameters,poped.db){ with(as.list(parameters),{ A_ini <- c(A1 = DOSE, A2 = 0) times <- drop(xt) times <- sort(times) times <- c(0,times) # add extra time for start of integration out <- ode(A_ini, times, PK.1.comp.oral.ode, parameters) #,atol=1e-13,rtol=1e-13) y = out[,"A2"]/(V/Favail) y=y[-1] # remove initial time for start of integration y = cbind(y) ## must be a row vector return(list( y= y,poped.db=poped.db)) }) } ## -- parameter definition function PK.1.comp.oral.sd.fg.param.1 <- function(x,a,bpop,b,bocc){ parameters=c( V=bpop[1]*exp(b[1]), KA=bpop[2]*exp(b[2]), CL=bpop[3]*exp(b[3]), Favail=bpop[4], DOSE=a[1]) parameters["KE"]=parameters["CL"]/parameters["V"] return( parameters ) } PK.1.comp.oral.sd.fg.param.2 <- function(x,a,bpop,b,bocc){ parameters=c( V=bpop[1]*exp(b[1]), KA=bpop[2]*exp(b[2]), KE=bpop[3]*exp(b[3]), Favail=bpop[4], DOSE=a[1]) return( parameters ) } poped.db.1 <- create.poped.database(ff_fun="PK.1.comp.oral.sd.ff", fg_fun="PK.1.comp.oral.sd.fg.param.1", fError_fun="feps.add.prop", groupsize=32, m=1, sigma=diag(c(0.01,0.25)), bpop=c(V=8,KA=1,CL=0.15,Favail=1), d=c(V=0.02,KA=0.6,CL=0.07), notfixed_bpop=c(1,1,1,0), xt=c(1,2,3,6,24,36,72,120), minxt=0, maxxt=c(25,25,25,120,120,120,120,120), discrete_xt = list(1:120), a=cbind(c(70)), bUseGrouped_xt=1, maxa=c(200), mina=c(0)) poped.db.2 <- create.poped.database(poped.db.1, fg_fun="PK.1.comp.oral.sd.fg.param.2", bpop=c(V=8,KA=1,KE=0.15/8,Favail=1), d=c(V=0.02,KA=0.6,KE=0.07)) poped.db.3 <- create.poped.database(poped.db.1, ff_fun="PK.1.comp.oral.sd.ff.ode") ## evaluate initial designs res6a_popED <- PopED::evaluate_design(poped.db.1) # analytical solution with Cl res6b_popED <- PopED::evaluate_design(poped.db.2) # analytical solution with KE res6c_popED <- PopED::evaluate_design(poped.db.3) # ODEs ### mlxDesignEval # analytical solution with pkmodel() parametrized with Cl model6a <- inlineModel(" [INDIVIDUAL] input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka} DEFINITION: Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl} F = {distribution=normal, typical=F_pop, no-variability} V = {distribution=logNormal, typical=V_pop, sd=omega_V} ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka} [LONGITUDINAL] input = {a, b} input = {F, ka, V, Cl} PK: Cc = pkmodel(ka, V, Cl, p=F) OUTPUT: output = {Cc} DEFINITION: DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}") # analytical solution with pkmodel() parametrized with k model6b <- inlineModel(" [INDIVIDUAL] input = {k_pop, omega_k, F_pop, V_pop, omega_V, ka_pop, omega_ka} DEFINITION: k = {distribution=logNormal, typical=k_pop, sd=omega_k} F = {distribution=normal, typical=F_pop, no-variability} V = {distribution=logNormal, typical=V_pop, sd=omega_V} ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka} [LONGITUDINAL] input = {a, b} input = {F, ka, V, k} PK: Cc = pkmodel(ka, V, k, p=F) OUTPUT: output = {Cc} DEFINITION: DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}") # ODE model parametrized with Cl model6c <- inlineModel(" [INDIVIDUAL] input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka} DEFINITION: Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl} F = {distribution=normal, typical=F_pop, no-variability} V = {distribution=logNormal, typical=V_pop, sd=omega_V} ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka} [LONGITUDINAL] input = {a, b} input = {F, ka, V, Cl} PK: depot(target=Ad) EQUATION: odeType=stiff ddt_Ad = -ka*Ad ddt_Ac = ka*Ad - (Cl/V)*Ac Cc = Ac/V OUTPUT: output = {Cc} DEFINITION: DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}") The same design is evaluated for the three models. # analytical solution with Cl res6a_mlx <- mlxDesignEval::evaluate_design( model_file = model6a, population_parameters = data.frame(F_pop=1, ka_pop=1, V_pop=8, Cl_pop=0.15, omega_ka=sqrt(0.6), omega_V=sqrt(0.02), omega_Cl=sqrt(0.07), a=sqrt(0.25),b=sqrt(0.01)), fixed_parameters = c("F_pop"), group = list(size=32), treatment = list(data=data.frame(time=0,amount=70)), output = list(output="DV", data=data.frame(time=c(1,2,3,6,24,36,72,120))), rse_on_variance = T) # analytical solution with k res6b_mlx <- mlxDesignEval::evaluate_design( model_file = model6b, population_parameters = data.frame(F_pop=1, ka_pop=1, V_pop=8, k_pop=0.15/8, omega_ka=sqrt(0.6), omega_V=sqrt(0.02), omega_k=sqrt(0.07), a=sqrt(0.25),b=sqrt(0.01)), fixed_parameters = c("F_pop"), group = list(size=32), treatment = list(data=data.frame(time=0,amount=70)), output = list(output="DV", data=data.frame(time=c(1,2,3,6,24,36,72,120))), rse_on_variance = T) # ODE system res6c_mlx <- mlxDesignEval::evaluate_design( model_file = model6c, population_parameters = data.frame(F_pop=1, ka_pop=1, V_pop=8, Cl_pop=0.15, omega_ka=sqrt(0.6), omega_V=sqrt(0.02), omega_Cl=sqrt(0.07), a=sqrt(0.25),b=sqrt(0.01)), fixed_parameters = c("F_pop"), group = list(size=32), treatment = list(data=data.frame(time=0,amount=70)), output = list(output="DV", data=data.frame(time=c(1,2,3,6,24,36,72,120))), rse_on_variance = T) ### Comparison Analytical solution with Cl: #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 V 2.946875 V_pop 2.946876 4.858601e-05 TRUE #> 2 KA 14.572054 ka_pop 14.572062 5.500982e-05 TRUE #> 3 CL 5.096568 Cl_pop 5.096567 1.770528e-05 TRUE #> 4 d_V 34.369463 var_omega_V 34.369453 3.158243e-05 TRUE #> 5 d_KA 27.846026 var_omega_ka 27.846004 7.837298e-05 TRUE #> 6 d_CL 29.782073 var_omega_Cl 29.782069 1.165045e-05 TRUE #> 7 SIGMA[1,1] 26.863004 var_b 26.862986 6.770186e-05 TRUE #> 8 SIGMA[2,2] 25.642019 var_a 25.642028 3.589427e-05 TRUE Analytical solution with k: #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 V 2.946875 V_pop 2.946876 5.191786e-05 TRUE #> 2 KA 14.572054 ka_pop 14.572062 5.524263e-05 TRUE #> 3 KE 5.437768 k_pop 5.437766 3.526612e-05 TRUE #> 4 d_V 33.258184 var_omega_V 33.258192 2.252980e-05 TRUE #> 5 d_KA 27.727936 var_omega_ka 27.727963 9.708115e-05 TRUE #> 6 d_KE 32.722277 var_omega_k 32.722299 6.570502e-05 TRUE #> 7 SIGMA[1,1] 26.509409 var_b 26.509395 5.598068e-05 TRUE #> 8 SIGMA[2,2] 25.306924 var_a 25.306916 3.139526e-05 TRUE ODE system with Cl: #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 V 2.947023 V_pop 2.946876 4.989732e-03 TRUE #> 2 KA 14.572062 ka_pop 14.572062 5.305491e-06 TRUE #> 3 CL 5.097250 Cl_pop 5.096567 1.339638e-02 TRUE #> 4 d_V 34.368092 var_omega_V 34.369453 3.957485e-03 TRUE #> 5 d_KA 27.845897 var_omega_ka 27.846004 3.848681e-04 TRUE #> 6 d_CL 29.780671 var_omega_Cl 29.782069 4.694171e-03 TRUE #> 7 SIGMA[1,1] 26.864337 var_b 26.862986 5.030488e-03 TRUE #> 8 SIGMA[2,2] 25.645051 var_a 25.641997 1.190959e-02 TRUE ## Ex7: PK,1-comp,maturation This example uses a 1-comp PK model with weight and post-menstrual age (PMA) as a covariate on Cl and V. The weight is not an input, instead it is calculated from the PMA value with different formulas for males and females. Four groups with different PMA values are defined in the design. ### popED The function to calculate the weight from the post-menstrual age is defined as a separate function which is then called in the model. # typical WT prediction from Sumpter & Holford 2011 # PMA in years get_typical_weight <- function(PMA,SEX=2){ WTmax1 <- 2.76 TM50wt1 <- 38.5 HILLwt1 <- 12.9 HILL2wt1 <- 2.74 PMA1 <- PMA*52 WTmax2 <- 16.4 TM50wt2 <- 2.1 HILLwt2 <- 2.04 HILL2wt2 <- 1 WTmax3 <- 40.2 TM50wt3 <- 12.4 HILLwt3 <- 2.87 HILL2wt3 <- 0 TLAGwt4 <- 12.4 WTmax4 <- 33.6 THALFwt4 <- 3.61 FFEM <- 1 if(SEX==1) FFEM <- 0.884 wt1 <- FFEM*WTmax1/(1+(TM50wt1/PMA1)^((PMA1=TM50wt1)*HILL2wt1)) wt2 <- WTmax2/(1+(TM50wt2/PMA)^((PMA=TM50wt2)*HILL2wt2)) wt3 <- WTmax3/(1+(TM50wt3/PMA)^((PMA=TM50wt3)*HILL2wt3)) wt4 <- (PMA > TLAGwt4)*(FFEM*WTmax4*(1-exp(-log(2)/THALFwt4*(PMA-TLAGwt4)))) return(wt1 + wt2 + wt3 + wt4) } PK.1.comp.maturation.ff <- function(model_switch,xt,parameters,poped.db){ with(as.list(parameters),{ y=xt WT <- get_typical_weight(PMA,SEX=SEX) CL=CL*(WT/70)^(3/4)*(PMA^HILL)/(TM50^HILL+PMA^HILL) V=V*(WT/70) DOSE=1000*(WT/70) y = DOSE/V*exp(-CL/V*xt) return(list( y= y,poped.db=poped.db)) }) } PK.1.comp.maturation.fg <- function(x,a,bpop,b,bocc){ parameters=c( CL=bpop[1]*exp(b[1]), V=bpop[2]*exp(b[2]), TM50=bpop[3]*exp(b[3]), HILL=bpop[4], PMA=a[1], SEX=a[2]) return( parameters ) } poped.db <- create.poped.database(ff_file="PK.1.comp.maturation.ff", fError_file="feps.add.prop", fg_file="PK.1.comp.maturation.fg", groupsize=rbind(50,20,20,20), m=4, sigma=c(0.015,0.0015), notfixed_sigma = c(1,0), bpop=c(CL=3.8,V=20,TM50=60/52,HILL=3), d=c(CL=0.05,V=0.05,TM50=0.05), xt=c( 1,2,4,6,8,24), minxt=0, maxxt=24, bUseGrouped_xt=1, a=list(c(PMA=25,SEX=2), c(PMA=15,SEX=2), c(PMA=10,SEX=2), c(PMA=5,SEX=2)), maxa=list(c(PMA=30,SEX=2)), mina=list(c(PMA=1,SEX=2))) ## evaluate initial design res7_popED <- PopED::evaluate_design(poped.db) ### mlxDesignEval PMA and SEX are passed as regressors such that they can be handled in a flexible way in the structural model. The calculation of the typical weight for each PMA is done directly in the model. In this example, the dose is defined directly in the model and the explicit analytical formula is used. Note that it would also be possible to use the `pkmodel()` macro and define the treatment as part of the `evaluate_design` input. In that case, the scaling by weight could be done via the `p=` argument of `pkmodel`. model7 <- inlineModel(" [INDIVIDUAL] input = {CL_pop, omega_CL, V_pop, omega_V, TM50_pop, omega_TM50, HILL_pop} DEFINITION: CL = {distribution=logNormal, typical=CL_pop, sd=omega_CL} V = {distribution=logNormal, typical=V_pop, sd=omega_V} TM50 = {distribution=logNormal, typical=TM50_pop, sd=omega_TM50} HILL = {distribution=logNormal, typical=HILL_pop, no-variability} [LONGITUDINAL] input = {a, b} input = {CL,V,TM50,HILL,PMA,SEX} PMA = {use=regressor} SEX = {use=regressor} EQUATION: WTmax1 = 2.76 TM50wt1 = 38.5 HILLwt1 = 12.9 HILL2wt1 = 2.74 PMA1 = PMA*52 WTmax2 = 16.4 TM50wt2 = 2.1 HILLwt2 = 2.04 HILL2wt2 = 1 WTmax3 = 40.2 TM50wt3 = 12.4 HILLwt3 = 2.87 HILL2wt3 = 0 TLAGwt4 = 12.4 WTmax4 = 33.6 THALFwt4 = 3.61 if(SEX==1) FFEM = 0.884 else FFEM = 1 end if (PMA1=TM50wt1) wt1 = FFEM*WTmax1/(1+(TM50wt1/PMA1)^(HILLwt1+HILL2wt1)) else wt1 = FFEM*WTmax1/(1+(TM50wt1/PMA1)^(HILLwt1)) end else if (PMA1>=TM50wt1) wt1 = FFEM*WTmax1/(1+(TM50wt1/PMA1)^(HILL2wt1)) else wt1 = FFEM*WTmax1/(1+(TM50wt1/PMA1)^(0)) end end if (PMA=TM50wt2) wt2 = WTmax2/(1+(TM50wt2/PMA)^(HILLwt2+HILL2wt2)) else wt2 = WTmax2/(1+(TM50wt2/PMA)^(HILLwt2)) end else if (PMA>=TM50wt2) wt2 = WTmax2/(1+(TM50wt2/PMA)^(HILL2wt2)) else wt2 = WTmax2/(1+(TM50wt2/PMA)^(0)) end end if (PMA=TM50wt3) wt3 = WTmax3/(1+(TM50wt3/PMA)^(HILLwt3+HILL2wt3)) else wt3 = WTmax3/(1+(TM50wt3/PMA)^(HILLwt3)) end else if (PMA>=TM50wt2) wt3 = WTmax3/(1+(TM50wt3/PMA)^(HILL2wt3)) else wt3 = WTmax3/(1+(TM50wt3/PMA)^(0)) end end if (PMA > TLAGwt4) wt4 = (FFEM*WTmax4*(1-exp(-log(2)/THALFwt4*(PMA-TLAGwt4)))) else wt4 = 0*(FFEM*WTmax4*(1-exp(-log(2)/THALFwt4*(PMA-TLAGwt4)))) end WT = wt1 + wt2 + wt3 + wt4 CL1 = CL*(WT/70)^(3/4)*(PMA^HILL)/(TM50^HILL+PMA^HILL) V1 = V *(WT/70) DOSE=1000*(WT/70) Cc = DOSE/V1*exp(-CL1/V1*t) OUTPUT: output = {Cc} DEFINITION: y1 = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}") Each group has its own regressor element defining the PMA and SEX. There is no treatment definition in `evaluate_design` as the dose is defined as part of the structural model. reg1 <- data.frame(time=0,PMA=25, SEX=2) reg2 <- data.frame(time=0,PMA=15, SEX=2) reg3 <- data.frame(time=0,PMA=10, SEX=2) reg4 <- data.frame(time=0,PMA=5, SEX=2) g1 <-list(size=50, regressor=reg1) g2 <-list(size=20, regressor=reg2) g3 <-list(size=20, regressor=reg3) g4 <-list(size=20, regressor=reg4) res7_mlx <- mlxDesignEval::evaluate_design( model_file = model7, population_parameters = data.frame(CL_pop=3.8,V_pop=20,TM50_pop=60/52,HILL_pop=3, omega_CL=sqrt(0.05), omega_V=sqrt(0.05), omega_TM50=sqrt(0.05), a=sqrt(0.0015),b=sqrt(0.015)), fixed_parameters = c("a"), group = list(g1,g2,g3,g4), output = list(output="y1", data=data.frame(time=c(1,2,4,6,8,24))), rse_on_variance = T) ### Comparison #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 CL 3.724808 CL_pop 3.724809 1.139376e-05 TRUE #> 2 V 2.251799 V_pop 2.251799 9.947249e-06 TRUE #> 3 TM50 3048.507973 TM50_pop 3048.507832 4.620080e-06 TRUE #> 4 HILL 2117.356641 HILL_pop 2117.356843 9.528602e-06 TRUE #> 5 d_CL 15.751133 var_omega_CL 15.751139 3.949789e-05 TRUE #> 6 d_V 14.997327 var_omega_V 14.997333 3.759456e-05 TRUE #> 7 d_TM50 27968.366291 var_omega_TM50 27968.367565 4.556077e-06 TRUE #> 8 SIGMA[1,1] 6.950869 var_b 6.950875 7.528790e-05 TRUE ## Ex8: TMDD This example uses a TMDD model with quasi steady-state (QSS) approximation, encoded as an ODE system. ### popED The system is solved using deSolve in R (results with C++ using Rcpp are identical). The results are quite sensitive to the ODE solver tolerance. We fave set the tolerance in the popED script to be the same as the default Monolix tolerances. Starting from version 2024, the tolerances can also be changed in Monolix. tmdd_qss_one_target_model_ode <- function(Time,State,Pars){ with(as.list(c(State, Pars)), { RTOT = A4 CTOT= A2/V1 CFREE = 0.5*((CTOT-RTOT-KSSS)+sqrt((CTOT-RTOT-KSSS)^2+4*KSSS*CTOT)) dA1 = -KA*A1 dA2 = FAVAIL*KA*A1+(Q/V2)*A3-(CL/V1+Q/V1)*CFREE*V1-RTOT*KINT*CFREE*V1/(KSSS+CFREE) dA3 = (Q/V1)*CFREE*V1 - (Q/V2)*A3 dA4 = R0*KDEG - KDEG*RTOT - (KINT-KDEG)*(RTOT*CFREE/(KSSS+CFREE)) return(list(c(dA1,dA2,dA3,dA4))) }) } sfg <- function(x,a,bpop,b,bocc){ parameters=c( CL=bpop[1]*exp(b[1]) , V1=bpop[2]*exp(b[2]) , Q=bpop[3]*exp(b[3]) , V2=bpop[4]*exp(b[4]) , FAVAIL=bpop[5]*exp(b[5]) , KA=bpop[6]*exp(b[6]) , VMAX=bpop[7]*exp(b[7]) , KMSS=bpop[8]*exp(b[8]) , R0=bpop[9]*exp(b[9]) , KSSS=bpop[10]*exp(b[10]) , KDEG=bpop[11]*exp(b[11]) , KINT=bpop[12]*exp(b[12]) , DOSE=a[1] , SC_FLAG=a[2]) return(parameters) } tmdd_qss_one_target_model <- function(model_switch,xt,parameters,poped.db){ with(as.list(parameters),{ y=xt #The initialization vector for the compartment A_ini <- c(A1=DOSE*SC_FLAG, A2=DOSE*(1-SC_FLAG), A3=0, A4=R0) #Set up time points for the ODE times_xt <- drop(xt) times <- sort(times_xt) times <- c(0,times) ## add extra time for start of integration # solve the ODE out <- ode(A_ini, times, tmdd_qss_one_target_model_ode, parameters, atol=1e-9,rtol=1e-6) # as monolix default # extract the time points of the observations out = out[match(times_xt,out[,"time"]),] # Match ODE output to measurements RTOT = out[,"A4"] CTOT = out[,"A2"]/V1 CFREE = 0.5*((CTOT-RTOT-KSSS)+sqrt((CTOT-RTOT-KSSS)^2+4*KSSS*CTOT)) COMPLEX=((RTOT*CFREE)/(KSSS+CFREE)) RFREE= RTOT-COMPLEX y[model_switch==1]= RTOT[model_switch==1] y[model_switch==2]= CFREE[model_switch==2] #y[model_switch==3]=RFREE[model_switch==3] return(list( y=y,poped.db=poped.db)) }) } tmdd_qss_one_target_model_ruv <- function(model_switch,xt,parameters,epsi,poped.db){ returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db)) y <- returnArgs[[1]] poped.db <- returnArgs[[2]] y[model_switch==1] = log(y[model_switch==1])+epsi[,1] y[model_switch==2] = log(y[model_switch==2])+epsi[,2] #y[model_switch==3] = log(y[model_switch==3])+epsi[,3] return(list(y=y,poped.db=poped.db)) } ################################################# # for study 1 in gibiansky,JPKPD,2012 table 2 ################################################# # for study 1 in gibiansky,JPKPD,2012 table 2 poped.db.1 <- create.poped.database(ff_fun=tmdd_qss_one_target_model, fError_fun=tmdd_qss_one_target_model_ruv, fg_fun=sfg, groupsize=6, m=4, #number of groups sigma=c(0.04,0.0225), bpop=c(CL=0.3,V1=3,Q=0.2,V2=3,FAVAIL=0.7,KA=0.5,VMAX=0, KMSS=0,R0=0.1,KSSS=0.015,KDEG=10,KINT=0.05), d=c(CL=0.09,V1=0.09,Q=0.04,V2=0.04,FAVAIL=0.04,KA=0.16,VMAX=0, KMSS=0,R0=0.09,KSSS=0.09,KDEG=0.04,KINT=0.04), notfixed_bpop=c( 1,1,1,1,1,1,0,0,1,1,1,1), notfixed_d=c( 1,1,1,1,1,1,0,0,1,1,1,1), xt=c(0.0417,0.25,0.5,1,3,7,14,21,28,35,42,49,56, 0.0417,0.25,0.5,1,3,7,14,21,28,35,42,49,56), model_switch=c(1,1,1,1,1,1,1,1,1,1,1,1,1, 2,2,2,2,2,2,2,2,2,2,2,2,2), bUseGrouped_xt=1, G_xt=c(1,2,3,4,5,6,7,8,9,10,11,12,13, 1,2,3,4,5,6,7,8,9,10,11,12,13), a=list(c(DOSE=100, SC_FLAG=0), c(DOSE=300, SC_FLAG=0), c(DOSE=600, SC_FLAG=0), c(DOSE=1000, SC_FLAG=1)), discrete_a = list(DOSE=seq(100,1000,by=100), SC_FLAG=c(0,1))) res8a_popED <- PopED::evaluate_design(poped.db.1) ### mlxDesignEval The model in mlxtran language is always and automatically converted to C++ code to run fast. There is nothing to do on the user side. The solver for stiff ODEs has a better precision and is usually preferred. It can be set with `odeType=stiff` in the structural model. IV ans sub-cutaneous (SC) administrations are distinguished by the `adm` keyword. model8 <- inlineModel(" [INDIVIDUAL] input = {Q_pop, omega_Q, R0_pop, omega_R0, V2_pop, omega_V2, CL_pop, omega_CL, FAVAIL_pop, omega_FAVAIL, KA_pop, omega_KA, KDEG_pop, omega_KDEG, KINT_pop, omega_KINT, KSSS_pop, omega_KSSS, V1_pop, omega_V1} DEFINITION: Q = {distribution=logNormal, typical=Q_pop, sd=omega_Q} R0 = {distribution=logNormal, typical=R0_pop, sd=omega_R0} V2 = {distribution=logNormal, typical=V2_pop, sd=omega_V2} CL = {distribution=logNormal, typical=CL_pop, sd=omega_CL} FAVAIL = {distribution=logNormal, typical=FAVAIL_pop, sd=omega_FAVAIL} KA = {distribution=logNormal, typical=KA_pop, sd=omega_KA} KDEG = {distribution=logNormal, typical=KDEG_pop, sd=omega_KDEG} KINT = {distribution=logNormal, typical=KINT_pop, sd=omega_KINT} KSSS = {distribution=logNormal, typical=KSSS_pop, sd=omega_KSSS} V1 = {distribution=logNormal, typical=V1_pop, sd=omega_V1} [LONGITUDINAL] input = {aCfree, aRtot} input = {KA, FAVAIL, V1, KINT, KSSS, KDEG, R0, CL, Q, V2} PK: ; IV depot(target = A2, adm = 1) ; SC depot(target = A1, adm = 2) EQUATION: odeType = stiff A1_0 = 0 A2_0 = 0 A3_0 = 0 A4_0 = R0 ; Parameter transformations RTOT = A4 CTOT = A2/V1 CFREE = 0.5*((CTOT-RTOT-KSSS)+sqrt((CTOT-RTOT-KSSS)^2+4*KSSS*CTOT)) ddt_A1 = -KA*A1 ddt_A2 = FAVAIL*KA*A1+(Q/V2)*A3-(CL/V1+Q/V1)*CFREE*V1-RTOT*KINT*CFREE*V1/(KSSS+CFREE) ddt_A3 = (Q/V1)*CFREE*V1 - (Q/V2)*A3 ddt_A4 = R0*KDEG - KDEG*RTOT - (KINT-KDEG)*(RTOT*CFREE/(KSSS+CFREE)) OUTPUT: output = {CFREE, RTOT} DEFINITION: yCfree = {distribution=logNormal, prediction=CFREE, errorModel=constant(aCfree)} yRtot = {distribution=logNormal, prediction=RTOT, errorModel=constant(aRtot)}") We define 3 IV groups with ascending doses and one SC group. trt1 <- list(data=data.frame(time=0, amount=100), admID=1) trt2 <- list(data=data.frame(time=0, amount=300), admID=1) trt3 <- list(data=data.frame(time=0, amount=600), admID=1) trt4 <- list(data=data.frame(time=0, amount=1000), admID=2) g1 <- list(size=6, treatment=trt1) g2 <- list(size=6, treatment=trt2) g3 <- list(size=6, treatment=trt3) g4 <- list(size=6, treatment=trt4) outCfree <- list(output="yCfree", data=data.frame(time=c(0.0417,0.25,0.5,1,3,7,14,21,28,35,42,49,56))) outRtot <- list(output="yRtot", data=data.frame(time=c(0.0417,0.25,0.5,1,3,7,14,21,28,35,42,49,56))) res8_mlx <- mlxDesignEval::evaluate_design( model_file = model8, population_parameters = data.frame(KA_pop=0.5, V1_pop=3, CL_pop=0.3, FAVAIL_pop=0.7, KINT_pop=0.05, KSSS_pop=0.015, Q_pop=0.2, V2_pop=3, R0_pop=0.1, KDEG_pop=10, omega_KA=sqrt(0.16), omega_V1=sqrt(0.09), omega_CL=sqrt(0.09), omega_FAVAIL=sqrt(0.04),omega_KINT=sqrt(0.04), omega_KSSS=sqrt(0.09), omega_Q=sqrt(0.04), omega_V2=sqrt(0.04), omega_R0=sqrt(0.09), omega_KDEG=sqrt(0.04), aCfree=sqrt(0.0225), aRtot=sqrt(0.04)), group = list(g1,g2,g3,g4), output = list(outCfree,outRtot), rse_on_variance = T) ### Comparison #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 CL 7.125817 CL_pop 7.125686 1.844242e-03 TRUE #> 2 V1 6.866975 V1_pop 6.867008 4.873102e-04 TRUE #> 3 Q 7.409078 Q_pop 7.408239 1.132901e-02 TRUE #> 4 V2 10.435157 V2_pop 10.436235 1.032650e-02 TRUE #> 5 FAVAIL 11.471389 FAVAIL_pop 11.471443 4.658107e-04 TRUE #> 6 KA 19.592448 KA_pop 19.592448 3.111509e-06 TRUE #> 7 R0 8.408518 R0_pop 8.408258 3.089987e-03 TRUE #> 8 KSSS 11.031286 KSSS_pop 11.028962 2.107175e-02 TRUE #> 9 KDEG 8.176590 KDEG_pop 8.175696 1.094110e-02 TRUE #> 10 KINT 7.344912 KINT_pop 7.343644 1.726866e-02 TRUE #> 11 d_CL 32.695336 var_omega_CL 32.695158 5.452809e-04 TRUE #> 12 d_V1 33.548081 var_omega_V1 33.548108 7.942589e-05 TRUE #> 13 d_Q 74.871951 var_omega_Q 74.864811 9.535884e-03 TRUE #> 14 d_V2 84.900937 var_omega_V2 84.901826 1.046666e-03 TRUE #> 15 d_FAVAIL 98.434249 var_omega_FAVAIL 98.434293 4.446352e-05 TRUE #> 16 d_KA 78.750482 var_omega_KA 78.750556 9.337863e-05 TRUE #> 17 d_R0 36.358020 var_omega_R0 36.359746 4.746907e-03 TRUE #> 18 d_KSSS 49.621583 var_omega_KSSS 49.607974 2.742569e-02 TRUE #> 19 d_KDEG 66.473682 var_omega_KDEG 66.478869 7.803029e-03 TRUE #> 20 d_KINT 49.673985 var_omega_KINT 49.675044 2.131655e-03 TRUE #> 21 SIGMA[1,1] 8.935594 var_aRtot 8.935620 2.905946e-04 TRUE #> 22 SIGMA[2,2] 9.676886 var_aCfree 9.676859 2.780330e-04 TRUE ## Ex9: PK,2-comp,oral,ODE This is a two compartment model with oral absorption using ODEs. ### popED The system is solved using deSolve in R. PK.2.comp.oral.ode <- function(Time, State, Pars){ with(as.list(c(State, Pars)), { dA1 <- -KA*A1 dA2 <- KA*A1 + A3* Q/V2 -A2*(CL/V1+Q/V1) dA3 <- A2* Q/V1-A3* Q/V2 return(list(c(dA1, dA2, dA3))) }) } #' define the initial conditions and the dosing ff.PK.2.comp.oral.md.ode <- function(model_switch, xt, parameters, poped.db){ with(as.list(parameters),{ A_ini <- c(A1=0, A2=0, A3=0) times_xt <- drop(xt) dose_times = seq(from=0,to=max(times_xt),by=TAU) eventdat <- data.frame(var = c("A1"), time = dose_times, value = c(DOSE), method = c("add")) times <- sort(c(times_xt,dose_times)) out <- ode(A_ini, times, PK.2.comp.oral.ode, parameters, events = list(data = eventdat))#atol=1e-13,rtol=1e-13) y = out[, "A2"]/(V1/Favail) y=y[match(times_xt,out[,"time"])] y=cbind(y) return(list(y=y,poped.db=poped.db)) }) } #' parameter definition function #' names match parameters in function ff fg <- function(x,a,bpop,b,bocc){ parameters=c( CL=bpop[1]*exp(b[1]), V1=bpop[2], KA=bpop[3]*exp(b[2]), Q=bpop[4], V2=bpop[5], Favail=bpop[6], DOSE=a[1], TAU=a[2]) return( parameters ) } #' create poped database poped.db <- create.poped.database(ff_fun="ff.PK.2.comp.oral.md.ode", fError_fun="feps.add.prop", fg_fun="fg", groupsize=20, m=1, #number of groups sigma=c(prop=0.1^2,add=0.05^2), bpop=c(CL=10,V1=100,KA=1,Q= 3.0, V2= 40.0, Favail=1), d=c(CL=0.15^2,KA=0.25^2), notfixed_bpop=c(1,1,1,1,1,0), xt=c( 48,50,55,65,70,85,90,120), minxt=0, maxxt=144, discrete_xt = list(0:144), a=c(DOSE=100,TAU=24), maxa=c(DOSE=1000,TAU=24), mina=c(DOSE=0,TAU=8), discrete_a = list(DOSE=seq(0,1000,by=100),TAU=8:24)) res9_popED <- PopED::evaluate_design(poped.db) ### mlxDesignEval model9 <- inlineModel(" [INDIVIDUAL] input = {Cl_pop, omega_Cl, Q_pop, V1_pop, V2_pop, ka_pop, omega_ka} DEFINITION: Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl} Q = {distribution=normal, typical=Q_pop, no-variability} V1 = {distribution=normal, typical=V1_pop, no-variability} V2 = {distribution=normal, typical=V2_pop, no-variability} ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka} [LONGITUDINAL] input = {a, b} input = {ka, Cl, V1, Q, V2} PK: depot(target=Ad) EQUATION: odeType=stiff ; Parameter transformations V = V1 k12 = Q/V1 k21 = Q/V2 k = Cl/V1 ddt_Ad = -ka * Ad ddt_Ac = ka*Ad - k*Ac - k12*Ac + k21*Ap ddt_Ap = k12*Ac - k21*Ap Cc=Ac/V OUTPUT: output = {Cc} DEFINITION: DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}") res9_mlx <- mlxDesignEval::evaluate_design( model_file = model9, population_parameters = data.frame(ka_pop=1, V1_pop=100, Cl_pop=10, Q_pop=3, V2_pop =40, omega_ka=0.25, omega_Cl=0.15, a=0.05, b=0.1), group = list(size=20), treatment = list(data=data.frame(start=0, interval=24, nbDoses=6, amount=100)), output = list(output="DV", data=data.frame(time=c(48,50,55,65,70,85,90,120))), rse_on_variance = T) ### Comparison #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 CL 4.409665 Cl_pop 4.40966 1.204776e-04 TRUE #> 2 V1 12.071097 V1_pop 12.07108 1.341462e-04 TRUE #> 3 KA 32.825506 ka_pop 32.82545 1.773338e-04 TRUE #> 4 Q 65.765415 Q_pop 65.76515 3.965393e-04 TRUE #> 5 V2 77.649537 V2_pop 77.64977 3.058298e-04 TRUE #> 6 d_CL 36.561733 var_omega_Cl 36.56174 2.154027e-05 TRUE #> 7 d_KA 191.947572 var_omega_ka 191.94783 1.324427e-04 TRUE #> 8 sig_prop 64.576067 var_b 64.57616 1.385793e-04 TRUE #> 9 sig_add 20.574939 var_a 20.57494 1.310480e-05 TRUE ## Ex10: PKPD,HCV This example uses the hepatitis C virus (HCV) model from the publication [Nyberg et al., "Methods and software tools for design evaluation for population pharmacokinetics-pharmacodynamics studies", Br. J. Clin. Pharm., 2014.](https://bpspubs.onlinelibrary.wiley.com/doi/10.1111/bcp.12352). ### popED The popED example uses the compiled code. It is therefore not possible to see the equations. sfg <- function(x,a,bpop,b,bocc){ ## -- parameter definition function parameters=c(p=bpop[1], d=bpop[2], e=bpop[3], s=bpop[4], KA=bpop[5] + b[1], KE=bpop[6] + b[2], VD=bpop[7] + b[3], EC50=bpop[8] + b[4], n=bpop[9] + b[5], delta=bpop[10] + b[6], c=bpop[11] + b[7], DOSE=a[1], TINF=a[2], TAU=a[3]) return(parameters) } ff_ODE_compiled <- function(model_switch,xt,parameters,poped.db){ parameters[5:11] <- exp(parameters[5:11]) with(as.list(parameters),{ A_ini <- c(A1 = 0, A2 = 0, A3=c*delta/(p*e), A4=(s*e*p-d*c*delta)/(p*delta*e), A5=(s*e*p-d*c*delta)/(c*delta*e)) #Set up time points for the ODE times_xt <- drop(xt) times <- c(0,times_xt) ## add extra time for start of integration times <- sort(times) times <- unique(times) # remove duplicates # compute values from ODEs out <- ode(A_ini, times, func = "derivs", parms = parameters, #method="daspk", #jacfunc = "jac", dllname = "HCV_ode", initfunc = "initmod", #nout = 1, outnames = "Sum", atol=1e-12,rtol=1e-12) # grab timepoint values out = out[match(times_xt,out[,"time"]),] y <- xt*0 pk <- out[,"A2"]/VD pd <- out[,"A5"] y[model_switch==1] <- pk[model_switch==1] y[model_switch==2] <- log10(pd[model_switch==2]) y = cbind(y) # must be a column matrix return(list( y= y,poped.db=poped.db)) }) } feps_ODE_compiled <- function(model_switch,xt,parameters,epsi,poped.db){ ## -- Residual Error function MS<-model_switch y <- ff_ODE_compiled(model_switch,xt,parameters,poped.db)[[1]] pk.dv <- y + epsi[,1] pd.dv <- y + epsi[,2] y[MS==1] = pk.dv[MS==1] y[MS==2] = pd.dv[MS==2] return(list(y=y,poped.db=poped.db)) } ## -- Define initial design and design space poped_db_compiled <- create.poped.database(ff_file="ff_ODE_compiled", fg_file="sfg", fError_file="feps_ODE_compiled", bpop=c(p=100, d=0.001, e=1e-7, s=20000, KA=log(0.8), KE=log(0.15), VD=log(100),#VD=log(100000), EC50=log(0.12), #EC50=log(0.00012), n=log(2), delta=log(0.2), c=log(7)), notfixed_bpop=c(0,0,0,0,1,1,1,1,1,1,1), d=c(KA=0.25, KE=0.25, VD=0.25, EC50=0.25, n=0.25, delta=0.25, c=0.25), sigma=c(0.04,0.04), groupsize=30, xt=c(0,0.25,0.5,1,2,3,4,7,10,14,21,28, 0,0.25,0.5,1,2,3,4,7,10,14,21,28), model_switch=c(rep(1,12),rep(2,12)), a=c(180,1,7)) res10_popED <- PopED::evaluate_design(poped_db_compiled) ### mlxDesignEval model10 <- inlineModel(" [INDIVIDUAL] input = {lEC50_pop, omega_lEC50, lVd_pop, omega_lVd, lc_pop, omega_lc, ldelta_pop, omega_ldelta, lka_pop, omega_lka, lke_pop, omega_lke, ln_pop, omega_ln} DEFINITION: lEC50 = {distribution=normal, typical=lEC50_pop, sd=omega_lEC50} lVd = {distribution=normal, typical=lVd_pop, sd=omega_lVd} lc = {distribution=normal, typical=lc_pop, sd=omega_lc} ldelta = {distribution=normal, typical=ldelta_pop,sd=omega_ldelta} lka = {distribution=normal, typical=lka_pop, sd=omega_lka} lke = {distribution=normal, typical=lke_pop, sd=omega_lke} ln = {distribution=normal, typical=ln_pop, sd=omega_ln} [LONGITUDINAL] input = {aC, aW} input = {lka,lke,lVd,lEC50,ln,ldelta,lc} PK: depot(target = X, p = 1, adm = 1) EQUATION: odeType = stiff ka=exp(lka) ke=exp(lke) Vd=exp(lVd) EC50=exp(lEC50) n=exp(ln) delta=exp(ldelta) c=exp(lc) q=100 d=0.001 e=0.0000001 s=20000 ; Initial conditions t_0 = 0 X_0 = 0 T_0 = (c*delta)/(q*e) A_0 = 0 I_0 = (s*e*q-d*c*delta)/(q*delta*e) W_0 = (s*e*q-d*c*delta)/(c*delta*e) ; Ordinary Differential Equations C=A/Vd ddt_X = -ka*X ddt_A = ka*X-ke*A ddt_T = s-T*(e*W+d) ddt_I = e*W*T-delta*I ddt_W = q*(1-(max(C,0)^n)/(max(C,0)^n+EC50^n))*I-c*W log10W = log10(W) OUTPUT: output = {C,log10W} DEFINITION: yC = {distribution=normal, prediction=C, errorModel=constant(aC)} yW = {distribution=normal, prediction=log10W, errorModel=constant(aW)}") outPK <- list(output="yC", data=data.frame(time=c(0,0.25,0.5,1,2,3,4,7,10,14,21,28))) outPD <- list(output="yW", data=data.frame(time=c(0,0.25,0.5,1,2,3,4,7,10,14,21,28))) res10_mlx <- mlxDesignEval::evaluate_design( model_file = model10, population_parameters = data.frame(lka_pop=log(0.8), lke_pop=log(0.15), lVd_pop=log(100), lEC50_pop=log(0.12), ln_pop=log(2), ldelta_pop=log(0.2), lc_pop=log(7), omega_lka=sqrt(0.25), omega_lke=sqrt(0.25), omega_lVd=sqrt(0.25), omega_lEC50=sqrt(0.25), omega_ln=sqrt(0.25), omega_ldelta=sqrt(0.25), omega_lc=sqrt(0.25), aC=sqrt(0.04), aW=sqrt(0.04)), group = list(size=30), output = list(outPK,outPD), treatment = list(data=data.frame(start=0, interval=7, nbDoses=6, amount=180, tinf=1)), rse_on_variance = T) ### Comparison #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 KA 54.103880 lka_pop 54.103938 1.079865e-04 TRUE #> 2 KE 5.520538 lke_pop 5.520530 1.466914e-04 TRUE #> 3 VD 2.163313 lVd_pop 2.163313 9.106420e-06 TRUE #> 4 EC50 7.432080 lEC50_pop 7.432087 8.158977e-05 TRUE #> 5 n 15.079243 ln_pop 15.079235 5.244492e-05 TRUE #> 6 delta 5.833940 ldelta_pop 5.833941 1.394196e-05 TRUE #> 7 c 5.656054 lc_pop 5.656058 5.658128e-05 TRUE #> 8 d_KA 39.856357 var_omega_lka 39.856362 1.173282e-05 TRUE #> 9 d_KE 30.738289 var_omega_lke 30.738302 4.195786e-05 TRUE #> 10 d_VD 28.610608 var_omega_lVd 28.610599 3.133398e-05 TRUE #> 11 d_EC50 60.518673 var_omega_lEC50 60.518691 3.001940e-05 TRUE #> 12 d_n 28.979133 var_omega_ln 28.979137 1.595303e-05 TRUE #> 13 d_delta 27.179856 var_omega_ldelta 27.179845 3.854860e-05 TRUE #> 14 d_c 32.751496 var_omega_lc 32.751501 1.418497e-05 TRUE #> 15 SIGMA[1,1] 8.373173 var_aC 8.373171 1.446238e-05 TRUE #> 16 SIGMA[2,2] 9.164741 var_aW 9.164742 8.385020e-06 TRUE ## Ex11: PK, prior FIM This example shows the usage of a prior FIM. ### popED # This example shows how to include a prior FIM into the design evaluation. # We look at PK assessment in pediatrics, where we are mainly interested # in assessing if there is a substantial difference of more than 20% in # clearance (CL) between children and adults. # First we setup the general model, then the designs for adults and pediatrics, # and finally we can evaluate the separate designs and the pooled data. ##-- Model: One comp first order absorption ## -- Analytic solution for both multiple and single dosing ff <- function(model_switch,xt,parameters,poped.db){ with(as.list(parameters),{ y=xt N = floor(xt/TAU)+1 y=(DOSE*Favail/V)*(KA/(KA - CL/V)) * (exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) return(list( y=y,poped.db=poped.db)) }) } ## -- parameter definition function ## -- names match parameters in function ff ## -- note, covariate on clearance for pediatrics (using isPediatric x[1]) sfg <- function(x,a,bpop,b,bocc){ parameters=c( V=bpop[1]*exp(b[1]), KA=bpop[2]*exp(b[2]), CL=bpop[3]*exp(b[3])*bpop[5]^a[3], # add covariate for pediatrics Favail=bpop[4], isPediatric = a[3], DOSE=a[1], TAU=a[2]) return( parameters ) } ## -- Residual unexplained variablity (RUV) function ## -- Additive + Proportional feps <- function(model_switch,xt,parameters,epsi,poped.db){ returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db)) y <- returnArgs[[1]] poped.db <- returnArgs[[2]] y = y*(1+epsi[,1])+epsi[,2] return(list( y= y,poped.db =poped.db )) } ## -- Define design and design space for adults (isPediatric = 0) ## Two arms, 5 time points poped.db <- create.poped.database(ff_fun="ff", fg_fun="sfg", fError_fun="feps", bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,pedCL=0.8), notfixed_bpop=c(1,1,1,0,1), d=c(V=0.09,KA=0.09,CL=0.25^2), sigma=c(0.04,5e-6), notfixed_sigma=c(0,0), m=2, groupsize=20, xt=c( 1,8,10,240,245), bUseGrouped_xt=1, a=list(c(DOSE=20,TAU=24,isPediatric = 0), c(DOSE=40, TAU=24,isPediatric = 0))) # Note, to be able to use the adults FIM to combine with the pediatrics, # both have to have the parameter "pedCL" defined and set notfixed_bpop to 1. ## Define pediatric model/design (isPediatric = 1) ## One arm, 4 time points only poped.db.ped <- create.poped.database(ff_fun="ff", fg_fun="sfg", fError_fun="feps", bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,pedCL=0.8), notfixed_bpop=c(1,1,1,0,1), d=c(V=0.09,KA=0.09,CL=0.25^2), sigma=c(0.04,5e-6), notfixed_sigma=c(0,0), m=1, groupsize=6, xt=c( 1,2,6,240), bUseGrouped_xt=1, a=list(c(DOSE=40,TAU=24,isPediatric = 1))) # We can evaluate the adult design without warning, by setting the pedCL # parameter to be fixed (i.e., not estimated) res11_popED_adult <- PopED::evaluate_design(create.poped.database(poped.db, notfixed_bpop=c(1,1,1,0,0))) # One obtains good estimates for all parameters for adults (<60% RSE for all). ## evaluate design of pediatrics only - insufficient # Similarly as before with only pediatrics we cannot estimate the covariate effect, # so we fix it. res11_popED_ped <- PopED::evaluate_design(create.poped.database(poped.db.ped, notfixed_bpop=c(1,1,1,0,0))) # Due to having less subjects, less samples per subject, and only one dose level # the variability in pediatrics cannot be estimated well. ## Add adult prior # Now we combined the two sutdies, where we assume that we only assess a difference # between adults and pediatrics in CL. # We can set the prior FIM to the adult one: outAdult <- res11_popED_adult poped.db.all <- create.poped.database( poped.db.ped, prior_fim = outAdult$fim, notfixed_bpop=c(1,1,1,0,0) ) ## evaluate design using prior FIM from adults res11_popED_both <- PopED::evaluate_design(poped.db.all) # Obviously, the pooled data leads to much higher precision in parameter estimates # compared to the pediatrics only. ### mlxDesignEval The isPediatric indicator is passed as a regressor to adjust the clearance value between adults and children. model11 <- inlineModel(" [INDIVIDUAL] input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka, pedCL_pop} DEFINITION: Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl} F = {distribution=logNormal, typical=F_pop, no-variability} V = {distribution=logNormal, typical=V_pop, sd=omega_V} ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka} pedCL = {distribution=logNormal, typical=pedCL_pop, no-variability} [LONGITUDINAL] input = {a, b} input = {F, ka, V, Cl, pedCL, isPediatric} isPediatric = {use=regressor} PK: Clnew = Cl * pedCL^isPediatric Cc = pkmodel(ka, V, Cl=Clnew, p=F) OUTPUT: output = {Cc} DEFINITION: DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}") treatment1<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=20)) treatment2<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40)) g1 <- list(size=20, treatment=treatment1) g2 <- list(size=20, treatment=treatment2) res11_mlx_adult <- mlxDesignEval::evaluate_design( model_file = model11, population_parameters = data.frame(F_pop=0.9, ka_pop=0.25, V_pop=72.8, Cl_pop=3.75, pedCL_pop=0.8, omega_ka=sqrt(0.09), omega_V=sqrt(0.09), omega_Cl=0.25, a=sqrt(5e-6), b=sqrt(0.04)), fixed_parameters = c("F_pop","pedCL_pop","a","b"), group = list(g1,g2), regressor = data.frame(time=0,isPediatric=0), output = list(output="DV", data=data.frame(time=c( 1,8,10,240,245))), rse_on_variance = T) res11_mlx_ped <- mlxDesignEval::evaluate_design( model_file = model11, population_parameters = data.frame(F_pop=0.9, ka_pop=0.25, V_pop=72.8, Cl_pop=3.75, pedCL_pop=0.8, omega_ka=sqrt(0.09), omega_V=sqrt(0.09), omega_Cl=0.25, a=sqrt(5e-6), b=sqrt(0.04)), fixed_parameters = c("F_pop","pedCL_pop","a","b"), group = list(size=6), treatment = list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40)), regressor = data.frame(time=0,isPediatric=1), output = list(output="DV", data=data.frame(time=c( 1,2,6,240))), rse_on_variance = T) res11_mlx_both <- mlxDesignEval::evaluate_design( model_file = model11, population_parameters = data.frame(F_pop=0.9, ka_pop=0.25, V_pop=72.8, Cl_pop=3.75, pedCL_pop=0.8, omega_ka=sqrt(0.09), omega_V=sqrt(0.09), omega_Cl=0.25, a=sqrt(5e-6), b=sqrt(0.04)), fixed_parameters = c("F_pop","pedCL_pop","a","b"), group = list(size=6), treatment = list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40)), regressor = data.frame(time=0,isPediatric=1), output = list(output="DV", data=data.frame(time=c( 1,2,6,240))), rse_on_variance = T, prior_FIM = res11_mlx_adult$FIM) ### Comparison Adult trial: #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 V 6.634931 V_pop 6.634930 3.545338e-06 TRUE #> 2 KA 8.587203 ka_pop 8.587206 3.117019e-05 TRUE #> 3 CL 4.354792 Cl_pop 4.354790 2.426530e-05 TRUE #> 4 d_V 33.243601 var_omega_V 33.243613 3.433354e-05 TRUE #> 5 d_KA 55.689432 var_omega_ka 55.689417 2.726942e-05 TRUE #> 6 d_CL 27.133255 var_omega_Cl 27.133212 1.586636e-04 TRUE Pediatric trial: #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 V 24.72088 V_pop 24.72088 5.977031e-06 TRUE #> 2 KA 30.84953 ka_pop 30.84954 1.908690e-05 TRUE #> 3 CL 11.94767 Cl_pop 11.94768 8.210502e-05 TRUE #> 4 d_V 116.23095 var_omega_V 116.23100 4.371910e-05 TRUE #> 5 d_KA 181.19778 var_omega_ka 181.19774 2.271892e-05 TRUE #> 6 d_CL 77.29188 var_omega_Cl 77.29188 2.266350e-06 TRUE Pediatric trial taking into account knowledge from adult trial: #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 V 6.379075 V_pop 6.379074 3.848623e-06 TRUE #> 2 KA 8.219929 ka_pop 8.219932 2.814061e-05 TRUE #> 3 CL 4.086855 Cl_pop 4.086855 1.122750e-05 TRUE #> 4 d_V 31.808871 var_omega_V 31.808884 3.972277e-05 TRUE #> 5 d_KA 52.858399 var_omega_ka 52.858388 2.094961e-05 TRUE #> 6 d_CL 25.601551 var_omega_Cl 25.601515 1.415043e-04 TRUE ## Ex12: Covariates This example shows the usage of covariate distributions. The model includes the effect of weight on Cl and V. ### popED ## Introduction # Lets assume that we have a model with a covariate included in the model description. # Here we define a one-compartment PK model that has weight on both clearance and volume of distribution. mod_1 <- function(model_switch,xt,parameters,poped.db){ with(as.list(parameters),{ y=xt CL=CL*(WT/70)^(WT_CL) V=V*(WT/70)^(WT_V) DOSE=1000*(WT/70) y = DOSE/V*exp(-CL/V*xt) return(list( y= y,poped.db=poped.db)) }) } par_1 <- function(x,a,bpop,b,bocc){ parameters=c( CL=bpop[1]*exp(b[1]), V=bpop[2]*exp(b[2]), WT_CL=bpop[3], WT_V=bpop[4], WT=a[1]) return( parameters ) } # distribution of covariates: We set #`groupsize=1`, the number of groups to be 50 (`m=50`) and assume that WT is #sampled from a normal distribution with mean=70 and sd=10 set.seed(123456) weigths <- rnorm(50,mean = 70,sd=10) poped_db_2 <- create.poped.database(ff_fun=mod_1, fg_fun=par_1, fError_fun=feps.add.prop, groupsize=1, m=50, sigma=c(0.015,0.0015), notfixed_sigma = c(1,0), bpop=c(CL=3.8,V=20,WT_CL=0.75,WT_V=1), d=c(CL=0.05,V=0.05), xt=c( 1,2,4,6,8,24), minxt=0, maxxt=24, bUseGrouped_xt=1, a=as.list(weigths) ) res12b_popED <- PopED::evaluate_design(poped_db_2) ### mlxDesignEval The covariate effect is defined in the \[INDIVIDUAL\] block with a covariate transformation in the \[COVARIATE\] block to obtain a power-law relationship between the parameter and the covariate. Less standard parameter-covariate relationships can also be coded in the structural model using regressors, see our [dedicated documentation](https://monolixsuite.slp-software.com/monolix/?contextKey=complex-cov-param-relationships). model12 <- inlineModel(" [COVARIATE] input = WT EQUATION: tWT = log(WT/70) [INDIVIDUAL] input = {Cl_pop, omega_Cl, V_pop, omega_V, tWT, beta_V_tWT, beta_Cl_tWT} DEFINITION: Cl = {distribution=logNormal, typical=Cl_pop, covariate=tWT, coefficient=beta_Cl_tWT, sd=omega_Cl} V = {distribution=logNormal, typical=V_pop, covariate=tWT, coefficient=beta_V_tWT, sd=omega_V} [LONGITUDINAL] input = {a, b} input = {V, Cl} PK: Cc = pkmodel(V, Cl) OUTPUT: output = {Cc} DEFINITION: DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}") In R seed is fixed to use the same weights table in popED and mlxDesignEval. Instead of using one group per covariate value (as done in popED), we can simply define a data frame of weight values as `covariate` argument. # table of weight values and weight-based dosing set.seed(123456) covweigths <- rnorm(50,mean = 70,sd=10) res12b_mlx <- mlxDesignEval::evaluate_design( model_file = model12, population_parameters = data.frame(V_pop=20, Cl_pop=3.8, beta_V_tWT=1,beta_Cl_tWT=0.75, omega_V=sqrt(0.05), omega_Cl=sqrt(0.05), a=sqrt(0.0015),b=sqrt(0.015)), fixed_parameters = c("a"), group = list(size=50), covariate = data.frame(id=1:50, WT=covweigths), treatment = list(data=data.frame(time=0, amount=1000/70), scale=list(covariate="WT", intercept=0)), output = list(output="DV", data=data.frame(time=c(1,2,4,6,8,24))), rse_on_variance = T) ### Comparison #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 CL 3.249976 Cl_pop 3.249969 1.930655e-04 TRUE #> 2 V 3.319989 V_pop 3.319985 1.253870e-04 TRUE #> 3 WT_CL 29.001417 beta_Cl_tWT 29.001425 2.897006e-05 TRUE #> 4 WT_V 22.228796 beta_V_tWT 22.228788 3.603646e-05 TRUE #> 5 d_CL 21.026128 var_omega_Cl 21.026136 3.565315e-05 TRUE #> 6 d_V 21.953885 var_omega_V 21.953879 2.690191e-05 TRUE #> 7 SIGMA[1,1] 10.064637 var_b 10.064643 5.809050e-05 TRUE ## Ex13: Shrinkage It is not possible to calculate shrinkage with mlxDesignEval. ## Ex14: PK, IOV This example shows the usage of inter-occasion variability (IOV) on the clearance parameter. ### popED With popED, we need to define a different clearance parameter for each occasion and choose to use the first or second one using an if/else statement. cppFunction( 'List one_comp_oral_ode(double Time, NumericVector A, NumericVector Pars) { int n = A.size(); NumericVector dA(n); double CL_OCC_1 = Pars[0]; double CL_OCC_2 = Pars[1]; double V = Pars[2]; double KA = Pars[3]; double TAU = Pars[4]; double N,CL; N = floor(Time/TAU)+1; CL = CL_OCC_1; if(N>6) CL = CL_OCC_2; dA[0] = -KA*A[0]; dA[1] = KA*A[0] - (CL/V)*A[1]; return List::create(dA); }' ) ff.ode.rcpp <- function(model_switch, xt, parameters, poped.db){ with(as.list(parameters),{ A_ini <- c(A1=0, A2=0) times_xt <- drop(xt) #xt[,,drop=T] dose_times = seq(from=0,to=max(times_xt),by=TAU) eventdat <- data.frame(var = c("A1"), time = dose_times, value = c(DOSE), method = c("add")) times <- sort(c(times_xt,dose_times)) out <- ode(A_ini, times, one_comp_oral_ode, c(CL_OCC_1,CL_OCC_2,V,KA,TAU), events = list(data = eventdat))#atol=1e-13,rtol=1e-13) y = out[, "A2"]/(V) y=y[match(times_xt,out[,"time"])] y=cbind(y) return(list(y=y,poped.db=poped.db)) }) } ## -- parameter definition function ## -- names match parameters in function ff sfg <- function(x,a,bpop,b,bocc){ parameters=c( CL_OCC_1=bpop[1]*exp(b[1]+bocc[1,1]), CL_OCC_2=bpop[1]*exp(b[1]+bocc[1,2]), V=bpop[2]*exp(b[2]), KA=bpop[3]*exp(b[3]), DOSE=a[1], TAU=a[2]) return( parameters ) } ## -- Define design and design space poped.db <- create.poped.database(ff_fun=ff.ode.rcpp, fError_fun=feps.add.prop, fg_fun=sfg, bpop=c(CL=3.75,V=72.8,KA=0.25), d=c(CL=0.25^2,V=0.09,KA=0.09), sigma=c(0.04,5e-6), notfixed_sigma=c(0,0), docc = matrix(c(0,0.09,0),nrow = 1), m=2, groupsize=20, xt=c( 1,2,8,240,245), minxt=c(0,0,0,240,240), maxxt=c(10,10,10,248,248), bUseGrouped_xt=1, a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)), maxa=c(DOSE=200,TAU=24), mina=c(DOSE=0,TAU=24)) ## evaluate initial design res14_popED <- PopED::evaluate_design(poped.db) ### mlxDesignEval The IOV is defined in the \[INDIVIDUAL\] block. Omega_Cl represents the IIV and gamma_Cl the IOV. Given the lognormal distribution, the formula for Cl is the same as in popED. model14 <- inlineModel(" [INDIVIDUAL] input = {Cl_pop, omega_Cl, V_pop, omega_V, ka_pop, omega_ka, gamma_Cl} DEFINITION: Cl = {distribution=logNormal, typical=Cl_pop, varlevel={id, id*occ}, sd={omega_Cl, gamma_Cl}} V = {distribution=logNormal, typical=V_pop, sd=omega_V} ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka} [LONGITUDINAL] input = {a, b} input = {ka, V, Cl} PK: Cc = pkmodel(ka, V, Cl) OUTPUT: output = {Cc} DEFINITION: DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}") The occasions are defined in the `occasion` argument and repeated in the other elements as well. treatment1<- list(data=data.frame(occ=c(1,1,1,1,1,1,2,2,2,2,2) , time=c(0,24,48,72,96,120,144,168,192,216,240), amount=20)) treatment2<- list(data=data.frame(occ=c(1,1,1,1,1,1,2,2,2,2,2) , time=c(0,24,48,72,96,120,144,168,192,216,240), amount=40)) g1 <- list(size=20, treatment=treatment1) g2 <- list(size=20, treatment=treatment2) res14_mlx <- mlxDesignEval::evaluate_design( model_file = model14, population_parameters = data.frame(ka_pop=0.25, V_pop=72.8, Cl_pop=3.75, omega_ka=sqrt(0.09), omega_V=sqrt(0.09), omega_Cl=0.25, gamma_Cl=sqrt(0.09), a=sqrt(5e-6), b=sqrt(0.04)), fixed_parameters = c("a","b"), group = list(g1,g2), occasion = data.frame(time=c(0,144),occ=c(1,2)), output = list(output="DV", data=data.frame(occ=c(1,1,1,2,2),time=c(1,2,8,240,245))), rse_on_variance = T) ### Comparison The RSEs are identical, except a small difference for the IOV variance. #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 CL 6.225029 Cl_pop 6.244430 0.31166092 TRUE #> 2 V 8.817098 V_pop 8.826273 0.10406706 TRUE #> 3 KA 10.672172 ka_pop 10.680824 0.08106658 TRUE #> 4 d_CL 106.438935 var_omega_Cl 106.309736 0.12138269 TRUE #> 5 d_V 42.915831 var_omega_V 42.942908 0.06309316 TRUE #> 6 d_KA 62.884047 var_omega_ka 62.911048 0.04293874 TRUE #> 7 D.occ[1,1] 79.419819 var_gamma_Cl 78.406349 1.27609228 FALSE ## Ex15: full covariance matrix This example shows the inclusion of a full covariance matrix, i.e correlations between the random effects. ### popED With popED, the covariance terms are defined in the `covd` argument of the database. ff <- function(model_switch,xt,parameters,poped.db){ ##-- Model: One comp first order absorption with(as.list(parameters),{ y=xt y=(DOSE*Favail*KA/(V*(KA-CL/V)))*(exp(-CL/V*xt)-exp(-KA*xt)) return(list(y=y,poped.db=poped.db)) }) } sfg <- function(x,a,bpop,b,bocc){ ## -- parameter definition function parameters=c(CL=bpop[1]*exp(b[1]), V=bpop[2]*exp(b[2]), KA=bpop[3]*exp(b[3]), Favail=bpop[4], DOSE=a[1]) return(parameters) } feps <- function(model_switch,xt,parameters,epsi,poped.db){ ## -- Residual Error function ## -- Proportional returnArgs <- ff(model_switch,xt,parameters,poped.db) y <- returnArgs[[1]] poped.db <- returnArgs[[2]] y = y*(1+epsi[,1]) return(list(y=y,poped.db=poped.db)) } ## -- With covariances poped.db_with <- create.poped.database(ff_file="ff", fg_file="sfg", fError_file="feps", bpop=c(CL=0.15, V=8, KA=1.0, Favail=1), notfixed_bpop=c(1,1,1,0), d=c(CL=0.07, V=0.02, KA=0.6), covd = c(.03,.1,.09), sigma=0.01, groupsize=32, xt=c( 0.5,1,2,6,24,36,72,120), minxt=0, maxxt=120, a=70) res15_popED <- PopED::evaluate_design(poped.db_with) ### mlxDesignEval With mlxDesignEval, we define correlations instead of covariances, and this is done in the \[INDIVIDUAL\] block. model15 <- inlineModel(" [INDIVIDUAL] input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka, corr_ka_V, corr_V_Cl, corr_ka_Cl} DEFINITION: Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl} F = {distribution=logNormal, typical=F_pop, no-variability} V = {distribution=logNormal, typical=V_pop, sd=omega_V} ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka} correlation = {level=id, r(V, Cl)=corr_V_Cl, r(ka, Cl)=corr_ka_Cl, r(ka, V)=corr_ka_V} [LONGITUDINAL] input = {b} input = {F, ka, V, Cl} PK: Cc = pkmodel(ka, V, Cl, p=F) OUTPUT: output = {Cc} DEFINITION: DV = {distribution=normal, prediction=Cc, errorModel=proportional(b)}") To calculate the correlation value corresponding to the covariance value in popED, use: corr_i_j = covar_i_j/(sqrt(var_i) \* sqrt(var_j)) res15_mlx <- mlxDesignEval::evaluate_design( model_file = model15, population_parameters = data.frame(F_pop=1, Cl_pop=0.15,V_pop=8, ka_pop=1, omega_Cl=sqrt(0.07),omega_V=sqrt(0.02),omega_ka=sqrt(0.6), corr_V_Cl=0.8017837,corr_ka_Cl=0.4879500, corr_ka_V=0.8215838, b=sqrt(0.01)), fixed_parameters = c("F_pop"), group = list(size=32), treatment = list(data=data.frame(time=0, amount=70)), output = list(output="DV", data=data.frame(time=c(0.5,1,2,6,24,36,72,120))), rse_on_variance = T) ### Comparison #> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical #> 1 CL 4.738266 Cl_pop 4.738270 6.837835e-05 TRUE #> 2 V 2.756206 V_pop 2.756208 4.766023e-05 TRUE #> 3 KA 13.925829 ka_pop 13.925839 7.575014e-05 TRUE #> 4 d_CL 25.660169 var_omega_Cl 25.660200 1.211581e-04 TRUE #> 5 d_V 30.482032 var_omega_V 30.482018 4.601698e-05 TRUE #> 6 d_KA 25.860002 var_omega_ka 25.860007 1.764386e-05 TRUE #> 7 D[2,1] 30.823749 covar_V_Cl 30.823745 1.409248e-05 TRUE #> 8 D[3,1] 41.469557 covar_ka_Cl 41.469561 1.042334e-05 TRUE #> 9 D[3,2] 30.732069 covar_ka_V 30.732062 2.328111e-05 TRUE #> 10 SIGMA[1,1] 11.180340 var_b 11.180340 5.004789e-12 TRUE Last updated: October 15, 2024 --- version: "2024R1" language: "en" --- # computeBins ## Compute bins Compute bins values, middles, and data repartition over bins. Available options are: |---------------|---------------|----------------------------------------------------------------------------------------------| | "criteria" | (*character*) | Bins criteria: "equalwidth", "equalsize" or "leastsquare" (default). | | "useFixedNb" | (*logical*) | TRUE to fix the number of bins, FALSE (default) to estimate it. | | "fixedNb" | (*integer*) | Fixed number of bins (default: 10). | | "estimatedNb" | (*pair*) | Bounds for bins number estimation (default: \[5,30\]). | | "nbBinData" | (*pair*) | Minimum and maximum number of data per bin for bins number estimation (default: \[10,200\]). | ### Usage R computeBins(data, options = list()) ### Arguments data (vector) Input data. options (list) \[optional\] Computation options. ### Value A list bins values ("values"), middles ("middles") and the actual number of data per bin ("repartition"). ### Examples R if (FALSE) { computeBins(data = c(1, 1.25, 2.5, 5, 5.5, 5.75, 7.5, 15, 16.5), options = list(nbBinData = c(1,10))) } Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # computeChartsData ## \[Monolix - PKanalix - Simulx\] Compute the charts data Compute (if needed) and export the charts data of a given plot or, if not specified, all the available project plots. ### Usage R computeChartsData(plot = NULL, output = NULL, exportVPCSimulations = NULL) ### Arguments plot (character) \[optional\]\[Monolix\] Plot type. If not specified, all the available project plots will be considered. Available plots: bivariatedataviewer, covariateviewer, outputplot, indfits, obspred, residualsscatter, residualsdistribution, vpc, npc, predictiondistribution, parameterdistribution, randomeffects, covariancemodeldiagnosis, covariatemodeldiagnosis, likelihoodcontribution, fisher, saemresults, condmeanresults, likelihoodresults. output (character) \[optional\]\[Monolix\] Plotted output (depending on the software, it can represent an observation, a simulation output, ...). By default, all available outputs are considered. exportVPCSimulations (logical) \[optional\]\[Monolix\] Should VPC simulations be exported if available. Equals FALSE by default. NOTE: If 'plot" argument is not provided, 'output' and "task' arguments are ignored. ### Details computeChartsData can be used to compute and export the charts data for plots available in the graphical user interface as in [Monolix](https://monolix.lixoft.com/graphics/), [PKanalix](https://pkanalix.lixoft.com/non-compartmental-analysis/#ncaplots) or [Simulx](https://simulx.lixoft.com/simulation/plots/), when you export \> export charts data. The exported charts data is saved as txt files in the result folder, in the ChartsData subfolder. Notice that it does not impact the current scenario. To get a ggplot equivalent to the plot in the GUI, but customizable in R with the ggplot2 library, better use one of the plot... functions available in the connectors for Monolix and PKanalix (not available for Simulx). To get the charts data for one of these plot functions as a dataframe, you can use [`getChartsData`](getchartsdata). ### See also [`getChartsData`](getchartsdata) ### Examples R if (FALSE) { computeChartsData() # Monolix - PKanalix - Simulx computeChartsData(plot = "vpc", output = "y1") # Monolix } Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # computePredictions ## \[Monolix\] Compute predictions from the structural model Compute predicted observation model values on observation times for each individual in a set of individuals. By default, the predictions are computed for all individuals. If the parameter `individualIds` is specified, the `individualParameters` must contain only the parameters for those individuals. That is, the number of rows in `individualParameters` must be the same as the length of `individualIds`. ### Usage R computePredictions(individualParameters, individualIds = NULL) ### Arguments individualParameters Individual parameter values for each parameter present in the model, either for all individuals or for the set of individuals specified in individualIds". This should be a data.frame with a column for individual id and a column for each parameter. individualIds \[optional\] vector Ids of the individuals for which observation models should be computed. By default, all the individuals are used. ### Value A list of predictions, where each prediction is a vector giving the computed prediction at observation times for each individual ### Details For each prediction, all individual values are returned as a single vector. Thus the number of observations per individual can be used to separate the predictions per individual. ### See also [`getIndividualParameterModel`](getindividualparametermodel) to get the individual parameter model used for the prediction [`getEstimatedIndividualParameters`](getestimatedindividualparameters) to get individual parameters to use in the prediction ### Examples R initializeLixoftConnectors("monolix") project_file <- file.path(getDemoPath(), "1.creating_and_using_models", "1.1.libraries_of_models", "theophylline_project.mlxtran") loadProject(project_file) runScenario() ids <- c(1,4) individualParamsForAllIndiv <- getEstimatedIndividualParameters()$saem predictions <- computePredictions( individualParameters = individualParamsForAllIndiv[ids,], individualIds = ids ) obsInfo = getObservationInformation() allIds = unique(obsInfo$y$id) obsPred = cbind(obsInfo$y[obsInfo$y$id %in% allIds[ids]], y_pred = predictions$Cc) Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # conc-QTc examples Here we present in more details several examples of concentration-QTc analysis using the provided R package. All example datasets and R script are available in the conc-QTc R package zip. * [Quinidine: crossover with placebo and drug periods](https://monolixsuite.slp-software.com/r-functions/2024R1/quinidine-crossover-with-placebo-and-drug-periods.md) * [Vanoxerine: identification of a delay](https://monolixsuite.slp-software.com/r-functions/2024R1/vanoxerine-identification-of-a-delay.md) * [Ranolazine: parallel design](https://monolixsuite.slp-software.com/r-functions/2024R1/ranolazine-parallel-design.md) * [No placebo and time matched baseline](https://monolixsuite.slp-software.com/r-functions/2024R1/no-placebo-and-time-matched-baseline.md) * [Multi-Day QTc Studies](https://monolixsuite.slp-software.com/r-functions/2024R1/multi-day-qtc-studies.md) Last updated: February 19, 2026 --- version: "2024R1" language: "en" --- # Conc-QTc modeling and implementation ## Conc-QTc Modeling ### Pre-specified linear model The [Garnett et al. White Paper](https://link.springer.com/article/10.1007/s10928-017-9558-5) suggests to use a linear model with concentration effect, treatment effect, time effect and centered baseline effect: with *i* the individual, *j* the treatment and *k* the time. The *ϴ* parameters represent fixed effects and the *η* terms indicate random effects. It is assumed the random effects are normally distributed with mean 0 and an unstructured covariance matrix. The residual error are normally distributed with mean 0 and a constant variance. ### Variations to the pre-specified model Specific situations require adaptations of the pre-specified model, as listed in Table 4 of [Garnett et al](https://doi.org/10.1007/s10928-017-9558-5). #### ΔΔQTc as dependent variable In crossover studies where each individual receives both placebo and active treatment, it is possible to compute ΔΔQTc for each individual by subtracting the ΔQTc for placebo from the ΔQTcF for each treatment arm. In this case, ΔΔQTc can be modeled directly, leading to a simpler model. With the intercept (ϴ~0,pop~+η~0,i~) and time point terms, ϴ~3,T1~\*I~(TIME=T1)~ + ϴ~3,T2~\*I~(TIME=T2)~ + ..., being identical for active and placebo treatment, these terms cancel out. The parameter ϴ~1~ denoting the treatment effect for ΔQTc becomes the intercept for ΔΔQTc (without random effects). A new covariate appears, the difference between centered baselines of the treatment and the placebo group. The model reads: with #### No placebo data If by design or ethical reasons a concurrent placebo arm was not included, it is not possible to estimate the ϴ~1~term associated to the placebo vs active treatment. The [Garnett et al](https://doi.org/10.1007/s10928-017-9558-5) white paper also suggests to remove the ϴ~3~ term associated to the time effect. When time-matched baseline record are used, diurnal variations are already corrected for and the ϴ~3~ term is indeed not needed (see next section). However, when datasets without placebo data use predose baselines, [Orihashi et al](https://link.springer.com/article/10.1007/s10928-021-09737-0) have shown that including categorical time effects in the model leads to less biased estimates and better controlled false negative rates. The model thus reads: #### Time-matched baseline Time-matched baseline adjustments use the corresponding predose QTc measurements collected prior to drug administration on Day -1. A time-matched baseline adjustment may be used when placebo data are not available, in order to minimize the effect of diurnal variation in QTc. In case of time-matched baseline, diurnal variations are already corrected for and the ϴ~3~ term is not needed. The model reads: Note that if there is no placebo data, the ϴ~1~term also drops out (see section above). #### Data over multiple days If QTc measurements are collected over multiple days, the TIME parameter can be reduced to time after dose on same day and DAY is included in model as an additional factor. The model reads: #### Non-linear concentration effect If exploratory plots indicate a non-linear effect of the concentration on ΔQTc, alternative models can be used. Typical non-linear models are Emax, Emax with sigmoidicity (Hill) and loglinear relationships. Emax model: Emax with sigmoidicity (Hill equation): Loglinear model: ### Implementation in Monolix Monolix is a tool for parameter estimation and model diagnostic for linear and non-linear mixed-effects models. It can thus be used to implement and use the models presented above. The implementation in Monolix requires to slightly reorganize the model into two distinct pieces (still corresponding to the same final equation). The implementation of the model will be done automatically in the R package. This section is for information purpose but does not need to be applied to the user himself. #### Separation into structural and statistical model In a pharmacometric framework, the pre-specified linear model can be rewritten as a structural model which described the drug effect on ΔQTc: and a statistical model, which describes the interindividual variability and covariate effects: This formulation allows for an easy extension to non-linear relationships between ΔQTc and the concentration, by defining alternative structural models. The implementation of the other variations to the pre-specified model are easily implemented by including less or different covariates. For ϴ~0,i~, we would like to consider random effects varying from individual to individual (*η* ~0,i~) and covariates which varies from period to period (TRT, BLQTc_cent) or from time point to time point (TIME). Due to a limitation in Monolix, this forces us to separate the fixed effects and the random effects into two different terms of the structural model. The structural model we will use is thus: and the individual model is: The term ϴ~0,i~ now contains only fixed effects: the intercept ϴ~0,pop~ and the covariate effects ϴ~1,~ϴ~3,k~, and ϴ~4~. The term *η* ~0,i~ represents the random effects for the intercept. We will fix its typical value *η* ~0,pop~ to zero such that the term contains only random effects. The slope remains as before with a random and a fixed component. #### Structural model in mlxtran language In Monolix, the structural model is defined via a txt file using mlxtran language. The concentration is passed as a regressor (to be read directly from the dataset) and its effect is defined in the structural model. The content of the mlxtran file for the linear model for ΔQTc is: [LONGITUDINAL] input = {Cc, dQTc0, eta_dQTc0, slope} Cc = {use=regressor} EQUATION: dQTc = dQTc0 + eta_dQTc0 + slope*Cc OUTPUT: output = dQTc with `dQTc0` corresponding to the ϴ~0,i~individual parameter (with fixed effect and covariate effects), `eta_dQTc0` corresponding to *η* ~0,i~ (random effects only) and `slope` corresponding to the ϴ~2,i~individual parameter (with fixed effect and random effects). #### Statistical model in the GUI The statistical model is defined in the graphical user interface (or via the lixoftConnectors). For ΔQTc, we will set: * `dQTc0`: normal distribution to ensure that covariate are added linearly into the model, no random effects, and covariates TRT (categorical), TIME (categorical) and BLcent (continuous), * `eta_dQTc`: normal distribution, with random effects, and typical population value fixed to zero, * `slope`: normal distribution with random effects. ![image-20260216-144621.png](https://monolixsuite.slp-software.com/__attachments/a_713fa1ef3c542461c51da76defb8a4436c67590a696c1f61cc32f4fde55ddadd/image-20260216-144621.png?cb=5d46774b731aeb35734878bcc8d213f8) The equations corresponding to this setup can be displayed with the "equation" icon: ![image-20260216-104130.png](https://monolixsuite.slp-software.com/__attachments/a_b696574948c8370a49bafa32da64f7f59c21845a7ecbbd5960b69cc1489cd78f/image-20260216-104130.png?cb=55ee7125d25faa5155f19eb009802f00) The typical value for the eta_dQTc0 parameter is fixed to zero in the Initial Estimates tab. ![image-20260216-104528.png](https://monolixsuite.slp-software.com/__attachments/a_f6c0d02f0d9a4278b61bede33b4f7c2f3478957f5fe99f761039e1e74fc92f6a/image-20260216-104528.png?cb=72cda293adf0b746067a138489006c4a) We consider correlations between the random effects (which corresponds to an unstructured covariance structure) and a constant residual error model, as suggested in Garnett et al. #### Parameter estimation Running the "population parameters" task will estimate the model parameters using the SAEM algorithm. The standard errors task will estimation the Fisher information matrix (FIM), the variance-covariance matrix of the estimates and the corresponding relative standard errors (RSE). They represent the uncertainty of the estimated parameters and will be used in the next section to compute the confidence interval. ## Model-derived ΔΔQTc at concentration(s) of interest ### Derivation of the ΔΔQTc equation The conc-ΔQTc model is used to compute the ΔΔQTc at concentrations of interest. The concentrations of interest are typically concentrations at the therapeutic dose in a patient population taking into account possible concentration increase due to drug interactions, impaired hepatic or renal function, and polymorphisms in CYP enzymes. The concentrations of interest should be inside the range of observed concentrations used to generate the model, such that the model is used for interpolation but not extrapolation. The PD metric of interest is ΔΔQTc, the change from baseline QTc adjusted for placebo. It is calculated as the difference between the model-derived ΔQTc at the concentrations of interest and model-derived ΔQTc for placebo with concentration = 0. Because the intercept (ϴ~0,pop~+η~0,i~) and time point terms, ϴ~3,T1~\*I~(TIME=T1)~ + ϴ~3,T2~\*I~(TIME=T2)~ + ..., are identical for the two terms, they cancel out. The parameter ϴ~1~ (treatment effect) which is present only for the first term remains. The BLQTc_cent is assumed to be null on average and is ignored. The random effects for the slope *η* ~2,i~ are also averaged by the mean. We thus obtain: with *Est* denoting the estimate of the true value. To calculate the confidence interval, an analytical formula can be derived when the linear model is used (see Garnett et al.), based on the standard errors of the estimates. However, this formula cannot be generalized when non-linear models are used. An alternative would be to use non-parametric bootstrap, but this would require many runs. In Monolix, we can use the variance-covariance matrix of the estimates (from which the standard errors are derived) to sample sets of (ϴ~1~,ϴ~2~) parameters which capture their uncertainty, instead of using bootstrap replicates. Each set is used to generate a ΔΔQTc prediction and the 90% confidence interval is calculated based on these predictions, as described in Garnett. et al. The resulting plot overlays the ΔΔQTc model prediction (black line), its 90% prediction interval capturing the uncertainty (blue band), the 10 ms increase in ΔΔQTc threshold (dashed line) and the upper limit of the prediction interval at the concentration of interest (blue arrow). ![image-20260216-144118.png](https://monolixsuite.slp-software.com/__attachments/a_06869c3481a5990d4eb31e21dcbad95ea8560a5d54ef1507722432186e9a3881/image-20260216-144118.png?cb=77201b9afa0c41f5335ab1d8d7ba0f61) ### Implementation in Simulx In order to reuse directly the developed model for ΔQTc, we will calculate the ΔΔQTc as the difference between the model-derived ΔQTc at the concentrations of interest and model-derived ΔQTc for placebo with concentration = 0. To perform the model prediction, the Monolix project is exported to Simulx. We will consider two simulation groups, one for the active treatment and one for the placebo. The simulation setup is the following: * **parameters**: the random effects are ignored and the fixed effects are sampled from the uncertainty distribution. This is done by selecting the element "mlx_TypicalUncertainLin". * **covariates**: the TRT covariates is set to 1 for the active group and 0 for the placebo. The values for the BLQTc_cent and TIME covariates are set to 0 and will anyway cancel out. * **concentration**: a fine grid of concentrations of interest is defined, via the regressor element. For the placebo group, the concentration is set to zero. * **replicates**: 500 replicates are used, each corresponding to a different set of population parameters sampled from the covariance matrix of the estimates (uncertainty distribution) * **group size**: for each replicate, a single simulation is done (as the random effects are ignored, all simulations would be identical) Outside of Simulx, we compute ΔΔQTc as the difference between the predictions of the two groups, for each replicates. The 5th, median and 95th percentiles over replicates is finally computed. This procedure for both linear and non-linear models and any included covariates. Last updated: February 26, 2026 --- version: "2024R1" language: "en" --- # confintmlx ## Overview ### Description Compute confidence intervals for the population parameters estimated by Monolix. The method used for computing the confidence intervals can be either based on the standard errors derived from an estimation of the Fisher Information Matrix ("fim"), on the profile likelihood ("proflike") or on nonparametric bootstrap estimate ("bootstrap"). is used by default. When method="fim", the FIM can be either estimated using a linearization of the model or a stochastic approximation. When method="proflike", the observed likelihood can be either estimated using a linearization of the model or an importance sampling Monte Carlo procedure. When method="bootstrap", the bootstrap estimates are obtained using the bootmlx function #### Usage r <- confintmlx(project, parameters="all", method="fim", level=0.90, linearization=TRUE, Nboot=100, settings=NULL) #### Arguments **project** a Monolix project **parameters** list of parameters for which confidence intervals are computed (default="all") **level** confidence level, a real number between 0 and 1 (default=0.90) **linearization** {TRUE}/FALSE whether the calculation of the standard errors or the profile likelihood is based on a linearization of the model (default=TRUE) **Nboot** number of bootstrat replicates (default=100, used when method="bootstrap") **settings** a list of optional settings * `max.iter` : maximum number of iterations to find the solution (default=10), * `tol.LL` : absolute tolerance for -2LL (default=0.001), * `tol.param` : relative tolerance for the parameter (default=0.01), * `print` : {TRUE}/FALSE display the results (default=TRUE) ## Example Compute confidence intervals using the standard errors derived from the estimated Fisher Information Matrix R library(Rsmlx) project <- "projects/warfarinPK1.mlxtran" r.fim <- confintmlx(project) print(r.fim) ## $confint ## estimate lower upper ## ka_pop 0.55941719 0.38420502 0.8145328 ## V_pop 7.77238591 7.39471923 8.1693410 ## beta_V_lw70 0.89284510 0.64956631 1.1361239 ## Cl_pop 0.13436970 0.12338053 0.1463377 ## omega_ka 0.74765531 0.50507055 1.1067532 ## omega_V 0.12073213 0.08098506 0.1799869 ## omega_Cl 0.27892262 0.22251102 0.3496358 ## a1 0.54294872 0.39471520 0.6911822 ## b1 0.07352272 0.04625716 0.1007883 ## ## $level ## [1] 0.9 ## ## $method ## [1] "fim" Compute confidence intervals using the profile likelihood method: R r.prl <- confintmlx(project, method="proflike", parameters = c("V_pop", "beta_V_lw70", "omega_V")) ## /**********************************************************************/ ## LL search on V_pop ## Upper bound search / Iteration 2 / V_pop = 9.493 ## Upper bound search / Iteration 3 / V_pop = 8.333 ## Upper bound search / Iteration 4 / V_pop = 8.34 ## Lower bound search / Iteration 2 / V_pop = 6.363 ## Lower bound search / Iteration 3 / V_pop = 7.34 ## Lower bound search / Iteration 4 / V_pop = 7.371 ## parameter V_pop ## Value 7.772 ## CI = [7.371 , 8.333] ## diff. = [-0.401 , 0.56] ## rel. diff. = [-5.157 , 7.216] ## /**********************************************************************/ ## LL search on beta_V_lw70 ## Upper bound search / Iteration 2 / beta_V_lw70 = 1.092 ## Upper bound search / Iteration 3 / beta_V_lw70 = 1.176 ## Lower bound search / Iteration 2 / beta_V_lw70 = 0.692 ## Lower bound search / Iteration 3 / beta_V_lw70 = 0.628 ## Lower bound search / Iteration 4 / beta_V_lw70 = 0.651 ## parameter beta_V_lw70 ## Value 0.892 ## CI = [0.651 , 1.176] ## diff. = [-0.242 , 0.283] ## rel. diff. = [-27.032 , 31.801] ## /**********************************************************************/ ## LL search on omega_V ## Upper bound search / Iteration 2 / omega_V = 0.147 ## Upper bound search / Iteration 3 / omega_V = 0.191 ## Upper bound search / Iteration 4 / omega_V = 0.171 ## Upper bound search / Iteration 5 / omega_V = 0.175 ## Lower bound search / Iteration 2 / omega_V = 0.098 ## Lower bound search / Iteration 3 / omega_V = 0.023 ## Lower bound search / Iteration 4 / omega_V = 0.077 ## Lower bound search / Iteration 5 / omega_V = 0.074 ## Lower bound search / Iteration 6 / omega_V = 0.07 ## Lower bound search / Iteration 7 / omega_V = 0.042 ## Lower bound search / Iteration 8 / omega_V = 0.069 ## Lower bound search / Iteration 9 / omega_V = 0.07 ![Untitled-20250715-125113.png](https://monolixsuite.slp-software.com/__attachments/a_6c6523658a08ecc89704ac3343905c23f06249d3258a4cfa0e8a67cd4ee5e6e7/Untitled-20250715-125113.png?cb=1e3b49c73f0f035e0f017ed69102add0) ## parameter omega_V ## Value 0.12 ## CI = [0.07 , 0.175] ## diff. = [-0.051 , 0.054] ## rel. diff. = [-41.675 , 45.187] ## /**********************************************************************/ ## parameter V_pop ## Value 7.772 ## CI = [7.371 , 8.333] ## diff. = [-0.401 , 0.56] ## rel. diff. = [-5.157 , 7.216] ## /**********************************************************************/ ## parameter beta_V_lw70 ## Value 0.892 ## CI = [0.651 , 1.176] ## diff. = [-0.242 , 0.283] ## rel. diff. = [-27.032 , 31.801] ## /**********************************************************************/ ## parameter omega_V ## Value 0.12 ## CI = [0.07 , 0.175] ## diff. = [-0.051 , 0.054] ## rel. diff. = [-41.675 , 45.187] R print(r.prl) ## $confint ## estimate lower upper ## V_pop 7.7723900 7.37158778 8.3333183 ## beta_V_lw70 0.8928451 0.65149574 1.1767846 ## omega_V 0.1207321 0.07041813 0.1752875 ## ## $proflike ## $proflike[[1]] ## param paramInit -2LL name thresh tol.param tol.LL ## 1 6.363495 7.77239 852.2823 V_pop 2.705543 0.01 0.1 ## 2 7.340131 7.77239 824.9214 V_pop 2.705543 0.01 0.1 ## 3 7.371588 7.77239 824.2444 V_pop 2.705543 0.01 0.1 ## 4 7.772390 7.77239 819.2299 V_pop 2.705543 0.01 0.1 ## 5 8.333318 7.77239 821.3445 V_pop 2.705543 0.01 0.1 ## 6 8.340940 7.77239 821.3185 V_pop 2.705543 0.01 0.1 ## 7 9.493219 7.77239 841.5167 V_pop 2.705543 0.01 0.1 ## useLinearization ## 1 TRUE ## 2 TRUE ## 3 TRUE ## 4 TRUE ## 5 TRUE ## 6 TRUE ## 7 TRUE ## ## $proflike[[2]] ## param paramInit -2LL name thresh tol.param tol.LL ## 1 0.6289709 0.8928451 822.4290 beta_V_lw70 2.705543 0.01 0.1 ## 2 0.6514957 0.8928451 821.8525 beta_V_lw70 2.705543 0.01 0.1 ## 3 0.6928451 0.8928451 820.7842 beta_V_lw70 2.705543 0.01 0.1 ## 4 0.8928451 0.8928451 819.2299 beta_V_lw70 2.705543 0.01 0.1 ## 5 1.0928451 0.8928451 820.5722 beta_V_lw70 2.705543 0.01 0.1 ## 6 1.1767846 0.8928451 821.9819 beta_V_lw70 2.705543 0.01 0.1 ## useLinearization ## 1 TRUE ## 2 TRUE ## 3 TRUE ## 4 TRUE ## 5 TRUE ## 6 TRUE ## ## $proflike[[3]] ## param paramInit -2LL name thresh tol.param tol.LL ## 1 0.02330622 0.1207321 828.9542 omega_V 2.705543 0.01 0.1 ## 2 0.04240121 0.1207321 826.0104 omega_V 2.705543 0.01 0.1 ## 3 0.06921732 0.1207321 822.3420 omega_V 2.705543 0.01 0.1 ## 4 0.07004615 0.1207321 822.2433 omega_V 2.705543 0.01 0.1 ## 5 0.07041813 0.1207321 821.7530 omega_V 2.705543 0.01 0.1 ## 6 0.07463491 0.1207321 821.3072 omega_V 2.705543 0.01 0.1 ## 7 0.07723433 0.1207321 821.5800 omega_V 2.705543 0.01 0.1 ## 8 0.09884708 0.1207321 819.1226 omega_V 2.705543 0.01 0.1 ## 9 0.12073210 0.1207321 819.2299 omega_V 2.705543 0.01 0.1 ## 10 0.14746252 0.1207321 819.7394 omega_V 2.705543 0.01 0.1 ## 11 0.17149316 0.1207321 821.5071 omega_V 2.705543 0.01 0.1 ## 12 0.17528745 0.1207321 821.8978 omega_V 2.705543 0.01 0.1 ## 13 0.19141786 0.1207321 823.7563 omega_V 2.705543 0.01 0.1 ## useLinearization ## 1 TRUE ## 2 TRUE ## 3 TRUE ## 4 TRUE ## 5 TRUE ## 6 TRUE ## 7 TRUE ## 8 TRUE ## 9 TRUE ## 10 TRUE ## 11 TRUE ## 12 TRUE ## 13 TRUE ## ## ## $level ## [1] 0.9 ## ## $method ## [1] "proflike" ![Untitled-1-20250715-125206.png](https://monolixsuite.slp-software.com/__attachments/a_60401a6b79e86761f0ef8077ded8ed5079e67783b8fbedd3ca286481307b3373/Untitled-1-20250715-125206.png?cb=1e3b49c73f0f035e0f017ed69102add0) Compute confidence intervals using the bootstrap method: R r.boot <- confintmlx(project, method="bootstrap", nboot=20) ## Generating data sets with initial data set resampling... ## Generating projects with bootstrap data sets... ## Project 1/20 => Population parameters already estimated ## Project 2/20 => Population parameters already estimated ## Project 3/20 => Population parameters already estimated ## Project 4/20 => Population parameters already estimated ## Project 5/20 => Population parameters already estimated ## Project 6/20 => Estimating the population parameters ## Project 7/20 => Estimating the population parameters ## Project 8/20 => Estimating the population parameters ## Project 9/20 => Estimating the population parameters ## Project 10/20 => Estimating the population parameters ## Project 11/20 => Estimating the population parameters ## Project 12/20 => Estimating the population parameters ## Project 13/20 => Estimating the population parameters ## Project 14/20 => Estimating the population parameters ## Project 15/20 => Estimating the population parameters ## Project 16/20 => Estimating the population parameters ## Project 17/20 => Estimating the population parameters ## Project 18/20 => Estimating the population parameters ## Project 19/20 => Estimating the population parameters ## Project 20/20 => Estimating the population parameters print(r.boot) ## $confint ## estimate lower upper ## ka_pop 0.55941719 0.43381532 0.9713527 ## V_pop 7.77238591 7.32236653 8.2168901 ## beta_V_lw70 0.89284510 0.65441378 1.2730923 ## Cl_pop 0.13436970 0.12447216 0.1446367 ## omega_ka 0.74765531 0.56088014 0.8638393 ## omega_V 0.12073213 0.07689146 0.1436114 ## omega_Cl 0.27892262 0.21887904 0.3196775 ## a1 0.54294872 0.31344991 0.6793549 ## b1 0.07352272 0.03661033 0.1084073 ## ## $level ## [1] 0.9 ## ## $method ## [1] "bootstrap" Last updated: July 15, 2025 --- version: "2024R1" language: "en" --- # Convergence assessment ## Using lixoftConnectors The convergence assessment can be run from the command line using the corresponding connector. As in the GUI, the following options can be chosen via the settings: the number of runs, the estimation of the standard errors and likelihood, linearization method and bounds for the sampling of the initial values. R # load and initialize the API library(lixoftConnectors) initializeLixoftConnectors(software="monolix") # load a project from the demos project <- paste0(getDemoPath(), "/1.creating_and_using_models/1.1.libraries_of_models/theophylline_project.mlxtran") loadProject(projectFile = project) # get the default settings and set new values set <- getAssessmentSettings() set$nbRuns <- 10 set$extendedEstimation <- T set$useLin <- T set$initialParameters <- data.frame(parameters = c("ka_pop","V_pop","Cl_pop"), fixed = FALSE, min=c(0.1, 0.1, 0.01), max=c(10,10,1)) # run the convergence assessment runAssessment(settings = set) # retrieve the results getAssessmentResults() When using the runAssessment() connector, its is possible to sample the initial value only for the fixed effects (except betas of covariate effects) and the sampling is uniform over an interval. If other strategies are needed, it is possible to implement a custom convergence assessment, as shown below. ## Custom convergence assessment The example below shows the functions to build a project from scratch using one of the demo data sets and a model from the libraries, and run a convergence assessment to evaluate the robustness of the convergence. Compared to the built-in convergence assessment, the strategy below samples new initial values for all parameters instead of fixed effects only, and it does not change the random seed: R # load and initialize the API library(lixoftConnectors) initializeLixoftConnectors(software="monolix") # get folder containing demo datasets demoPath = paste0(getDemoPath(), '/1.creating_and_using_models/1.1.libraries_of_models/') # get name of model from library model <- getLibraryModelName(library="pk", filters = list(administration="oral", delay="lagTime", absorption = "firstOrder", distribution = "1compartment", elimination = "linear", parametrization = "clearance")) # create a new project by setting a data set and a structural model newProject(data = list(dataFile = paste0(demoPath,'data/warfarin_data.csv'), headerTypes =c("id", "time", "amount", "observation", "obsid", "contcov", "catcov", "ignore")), mapping = list(list(obsId="1", observationName="y1", modelOutput="Cc")), modelFile = model) # set tasks in scenario scenario <- getScenario() scenario$tasks = c(populationParameterEstimation = T, conditionalModeEstimation = T, conditionalDistributionSampling = T, standardErrorEstimation=T, logLikelihoodEstimation=T) scenario$linearization = TRUE setScenario(scenario) # ---------------------------------------------------------------------------- # convergence assessment: run 5 estimations with different initial estimates, # store the results in tabestimates # ---------------------------------------------------------------------------- popparams <- getPopulationParameterInformation() tabestimates <- NULL; tabiters <- NULL for(i in 1:5){ # sample new initial estimates popini <- sapply(1:nrow(popparams), function(j){runif(n=1, min=popparams$initialValue[j]/2, max=popparams$initialValue[j]*2)}) # set sampled values as new initial estimates newpopparams <- popparams newpopparams$initialValue <- popini setPopulationParameterInformation(newpopparams) # run the estimation runScenario() # store the estimates and s.e. in a table estimates <- as.data.frame(getEstimatedPopulationParameters()) names(estimates) <- "estimate" rses <- getEstimatedStandardErrors()$linearization$rse names(rses) <- getEstimatedStandardErrors()$linearization$parameter rses <- as.data.frame(rses) estimates <- merge(estimates, rses, by = "row.names") estimates$run <- i names(estimates)[names(estimates) == "Row.names"] <- "param" tabestimates <- rbind(tabestimates, estimates) # store the iterations iters <- getChartsData("plotSaem") iters$run <- i tabiters <- rbind(tabiters, iters) } # load plotting libraries library(ggplot2) library(gridExtra) # plot SAEM iterations plotList <- list() i <- 1 for (param in popparams$name) { if (popparams[popparams$name == param, ]$method == "FIXED") next changePhase <- tabiters$iteration[which(diff(tabiters$phase) == 1) + 1] plotList[[i]] <- ggplot(tabiters, aes_string(x = "iteration", y = param)) + geom_line(aes(group = run, color = factor(run))) + theme(legend.position = "none", plot.title = element_text(hjust = .5)) + geom_vline(xintercept = changePhase, color = 1:length(changePhase)) + labs(title = param, x = NULL, y = NULL) i <- i + 1 } grid.arrange(grobs = plotList, ncol = 3) # plot population parameters plotList <- list() i <- 1 for (param in popparams$name) { if (popparams[popparams$name == param, ]$method == "FIXED") next estimates <- tabestimates[tabestimates$param == param, ] plotList[[i]] <- ggplot(estimates, aes(x = run, y = estimate)) + geom_point(aes(color = factor(run))) + geom_errorbar(aes(ymax = estimate * (1+rses/100), ymin = estimate * (1-rses/100), color = factor(run))) + theme(legend.position = "none", plot.title = element_text(hjust = .5)) + labs(title = param, x = NULL, y = NULL) i <- i + 1 } grid.arrange(grobs = plotList, ncol = 3) ![iters-1-20240829-123153.svg](https://monolixsuite.slp-software.com/__attachments/a_41bf8586678fbd77a40616aa739adb08135d0f39056bc9fda099afa145fe48fa/iters-1-20240829-123153.svg?cb=e685297d21d88ea40bf0a0079ec3f274) ![params-20240829-123202.svg](https://monolixsuite.slp-software.com/__attachments/a_a770ad18285cf72893c5c43789c2e6b01c6edb986a0491ef18496800070dec59/params-20240829-123202.svg?cb=8755d8cd8f5169154a3576efdbefb904) Last updated: October 21, 2024 --- version: "2024R1" language: "en" --- # Covariate search The automatic covariate search procedures can be launched using the lixoftConnectors. Note that when reloading the project (in the GUI or via the lixoftConnectors) with the 2020R1, the covariate search results are displayed but the settings are lost. From the 2021R1 version on, the settings will be saved and reloaded. In the example below, we do the following steps: * load a base monolix project, * add covariate transformations, * get the default settings of the covariate search and modify them, * specify which covariates and parameters to test along with relationships which are lock-in or out, * run the covariate search, * get the results. As an example, we are using the demo 5.2 warfarin_covariate1_project.mlxtran as base project. R library(lixoftConnectors) initializeLixoftConnectors(software="monolix") # makes the link to the MonolixSuite installation # load base project (from the demos) loadProject("warfarin_covariate1_project.mlxtran") # wt and sex are defined as covariates. Add logtWT=log(wt/70) addContinuousTransformedCovariate(logtWt = "log(wt/70)" ) # save the project with the additional covariate under a new name saveProject("warfarin_covariate1_project_covsearch.mlxtran") # get the default covariate search settings set <- getModelBuildingSettings() # modify the settings set$strategy <- 'cossac' # "cossac" or "scm" set$useSambaBeforeCossac <- FALSE # if true, corresponds to "covSAMBA-COSSAC" in the GUI set$criterion <- "LRT" # "BIC" or "LRT" set$threshold$lrt[1] <- 0.01 # forward p-value threshold for LRT set$threshold$lrt[2] <- 0.001 # backward p-value threshold for LRT # set logtWt and sex as covariates to test on all parameters except Tlag set$covariates <- c("sex","logtWt") set$parameters <- c("ka","V","Cl") # in addition lock-out (never test it) "sex on ka" and lock-in (have it in all tested models) "logtWt on Cl" set$relationships[1,] <- c("ka", "sex", FALSE) # sex never included on ka set$relationships[2,] <- c("Cl", "logtWt", TRUE) # logtWt always included on Cl # all other relationships corresponding to set$covariates and set$parameters will be tested # run the covariate search runModelBuilding(settings=set) # get the results. The best model is indicated with $bestModel = T getModelBuildingResults() Last updated: August 29, 2024 --- version: "2024R1" language: "en" --- # createCustomNCAParameter ## \[PKanalix\] Create a new NCA parameter as a formula of existing parameters Create a new NCA parameter as a formula of existing parameters and covariates. Available arguments: |--------------------|-------------------------|-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------| | "name" | (*character*, required) | Name of the parameter. | | "formula" | (*character*, required) | Formula used to calculate the parameter. This mathematical expression can contain PKanalix names of parameters (including partial AUCs, but excluding other custom parameters), names of covariates, operators +, -, /, \* and \^, parentheses and functions abs, sqrt, exp, log, log10, logit, invlogit, sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, floor, ceil, factorial, min, max, atan2, rem. Partial AUC parameters with non-integer times need to be specifically formatted (e.g., AUC_0_24.5 becomes AUC_0_24d5 and AUC_0_1e-07 becomes AUC_0_1em07). | | "alias" | (*character*, optional) | Stylized name of the parameter that will be displayed in results. Can contain tags and for subscript and superscript. | | "unit" | (*character*, optional) | Units of the parameter. It can contain substrings TIME, AMOUNT, VOLUME, CONC, GRADING to use data set units of time, dose amount, volume, concentration and normalization units. Different units should be combined with the character ".", and "\^-1" can be used. | | "previousName" | (*character*, optional) | Should be used only when editing an existing custom parameter. In that case the previous parameter is replaced with the new parameter. | | "addToPreferences" | (*logical*, optional) | Whether to add the custom NCA parameter to PKanalix preferences (default=FALSE). | ### Usage R createCustomNCAParameter( name, formula, alias = "", unit = "", addToPreferences = FALSE, previousName = "" ) ### See also [`deleteCustomNCAParameter`](deletecustomncaparameter)`, `[`getCustomNCAParameters`](getcustomncaparameters)`, `[`addCustomNCAParametersFromPreferences`](addcustomncaparametersfrompreferences) ### Examples R initializeLixoftConnectors("pkanalix") loadProject(paste0(getDemoPath(),"/1.basic_examples/project_covariates.pkx")) createCustomNCAParameter(name="CLperkg", formula = "Cl_F_obs/WEIGHT", alias = "CL/kg", unit = "VOLUME.TIME^-1.KG^1", addToPreferences = FALSE) Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # createFilter ## \[Monolix - PKanalix\] Create filter Create a new filtered data set by applying a filter on an existing one and/or complementing it. ### Usage R createFilter(filter, name = "", origin = "") ### Arguments filter (list\< list\< action = "headerName-comparator-value" \> \> or "complement") \[optional\] filter definition. Existing actions are "selectLines", "selectIds", "removeLines" and "removeIds". First vector level is for set unions, the second one for set intersection. It is possible to give only a list of actions if there is only no high-level union. name (character) \[optional\] created data set name. If not defined, the default name is "currentDataSet_filtered". origin (character) \[optional\] name of the data set to be filtered. The current one is used by default. ### Details The possible actions are line selection (selectLines), line removal (removeLines), Ids selection (selectIds) or removal (removeIds). The selection is a string containing the header name, a comparison operator and a value selection = "headerName\*-comparator\*\*-value" (ex: `"id=='100'"`, `"WEIGHT<70"`, `"SEX!='M'"`) Notice that : - The headerName corresponds to the data set header or one of the header aliases defined in MONOLIX software preferences - The comparator possibilities are "==", "!=" for all types of value and "\<=", "\<", "\>=", "\>" only for numerical types Syntax: \* create a simple filter: createFilter( filter = list(act = sel)), e.g. createFilter( filter = list(removeIds = "WEIGHT\<50")) =\> create a filter with the action act on the selection sel. In this example, we create a filter that removes all subjects with a weight less than 50. \* create a filter with several concurrent conditions, i.e AND condition: createFilter( list(act1 = sel1, act2 = sel2)), e.g. createFilter( filter = list(removeIds = "WEIGHT\<50", removeIds = " AGE\<20")) =\> create a filter with both the action act1 on sel1 AND the action act2 on sel2. In this example, we create a filter that removes all subjects with a weight less than 50 and an age less than 20. It corresponds to the intersecton of the subjects with a weight less than 50 and the subjects with an age less than 20. \* create a filter with several non-concurrent conditions, i.e OR condition: createFilter(filter = list(list(act1 = sel1), list(act2 = sel2)) ), e.g. createFilter( filter = list(list(removeIds = "WEIGHT\<50"),list(removeIds = " AGE\<20"))) =\> create a filter with the action act1 on sel1 OR the action act2 on sel2. In this example, we create a filter that removes all subjects with a weight less than 50 and an age less than 20. It corresponds to the union of the subjects with a weight less than 50 and the subjects with an age less than 20. \* It is possible to have any combinaison: createFilter(filter = list(list(act1 = sel1), list(act2 = sel2, act3 = sel3)) ) \<=\> act1,sel1 OR ( act2,sel2 AND act3,sel3 ) \* It is possible to create the complement of an existing filter: createFilter(filter = "complement") ### See also [`applyFilter`](applyfilter) ### Examples Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # createNCARatio ## \[PKanalix\] Create a ratio of NCA parameters across occasions. To use this, the dataset must have occasions. The ratios of parameters are calculated for each individual, across occasions. The ratio corresponds to the parameter value for the test modality divided by the parameter value for the reference modality. If certain subjects have multiple values of the test or reference value (for example, one occasion of R drug and two occasions of T drug), then the arithmetic mean of values of the subject-occasion results is taken. ### Usage R createNCARatio(name, variable, reference, test, parameter) ### Arguments name (character) Name of the new parameter defined as ratio of existing parameters. variable (character) Occasion or categorical covariate column. reference (character) Reference modality. test (character) Test modality. parameter (character) Name of the NCA parameter used in the ratio. ### See also [`deleteNCARatio`](deletencaratio)`, `[`getNCARatios`](getncaratios)`, `[`getNCAIndividualRatios`](getncaindividualratios)`, `[`getNCARatioStatistics`](getncaratiostatistics) ### Examples R initializeLixoftConnectors("pkanalix") loadProject(paste0(getDemoPath(),"/2.case_studies/project_Theo_extravasc_SD.pkx")) createNCARatio(name="AUCratio", variable="FORM", reference = "ref", test = "test", parameter = "AUCINF_obs") Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # Creating forest plots The process of creating forest plots based on a Monolix project is explained in this video: Below you can find an example R function that creates forest plots for visualization of covariate effects on a typical individual based on a Simulx project, as described in the video above. The function *plotCovariateEffects* takes two arguments: * *project* -- path to the Simulx project, * *outcome* -- name of the outcome that calculates the exposure parameter. Certain elements of created forest plots will be based on the names of simulation groups and outcomes in the used Simulx project. Here is an image describing the origin of those elements: ![image-20241011-123214.png](https://monolixsuite.slp-software.com/__attachments/a_b0a9f332cdfd013b589e873cc9946161950159a6593fab20fefc528a268a8b88/image-20241011-123214.png?cb=ef253e070b32eb3fc3735e6d78070dc5) R library(ggplot2) library(dplyr) library(lixoftConnectors) plotCovariateEffects <- function(project, outcome) { # Load a project initializeLixoftConnectors("simulx", force = TRUE) loadProject(project) # Get and format results results <- getEndpointsResults()$outcomes summary <- results[[outcome]] %>% group_by(group) %>% summarise( mid = quantile(.data[[outcome]], 0.5), lower = quantile(.data[[outcome]], 0.05), upper = quantile(.data[[outcome]], 0.95) ) reference_value <- summary[summary$group == "Reference", ]$mid summary <- summary %>% mutate(across(mid:upper, .fns = ~.x/reference_value)) %>% mutate( LABEL = paste0( format(round(mid, 2), nsmall = 2), " [", format(round(lower, 2), nsmall = 2), "-", format(round(upper, 2), nsmall = 2), "]" ) ) summary$covname <- gsub("_.*", "", summary$group) summary$label <- gsub("^(.*?)_", "", summary$group) summary$label <- gsub("_", " ", summary$label) %>% paste("\n", summary$LABEL) # Plot the results print( ggplot(data = summary[summary$covname != "Reference",], aes_string( y = "label", x = "mid", xmin = "lower", xmax = "upper" )) + geom_pointrange( aes(color = "90% CI\nCovariate Effects"), size = 1, alpha = 1 ) + annotate("rect", xmin = min(0.8), xmax = max(1.25), ymin = -Inf, ymax = Inf, fill = "gray", alpha=0.1) + geom_vline(aes(xintercept = 1, linetype = "Reference"), linetype = "dashed") + facet_grid(covname ~ ., scales = "free_y", switch = "y") + labs(y = "", x = paste(outcome, "Relative to Reference Value"), colour = "", linetype = "") + theme_bw() + scale_color_manual(values = c("blue")) ) } It is possible to automate the process of creating the Simulx project using *lixoftConnectors* . Here you can download two examples that start from a Monolix project and create a Simulx project and forest plots completely in R using *lixoftConnectors* : [download](https://simplus.sharefile.com/public/share/web-sbea0ea3b28c44ef7b7080be148355e43). Example 1 creates forest plots that visualize effect of covariates on a typical individual, while example 2 creates forest plots that illustrate individual exposure. The explanation of the process is described in this video: Last updated: February 21, 2025 --- version: "2024R1" language: "en" --- # defineCovariateElement ## \[Simulx\] Define covariate element Define or edit a covariate element. Covariate elements are defined and used for simulation [as in Simulx GUI](https://simulx.lixoft.com/definition/covariates/). As for other elements, the covariate elements can be defined or imported, and they are saved with the Simulx project if calling [`saveProject`](saveproject)`. Once a covariate element is defined, it needs to be added to a simulation group with `[`setGroupElement`](setgroupelement) to be used in simulation. ### Usage R defineCovariateElement(name, element) ### Arguments name (character) Element name. element (character or dataFrame or list) Element definition from external file path or data frame with covariates as columns, or list to select the sheet of an excel file:file (character) Path to the population file. sheet (character) Name of the sheet in xlsx/xls file. ### Details Covariate elements can be defined only if the model used in the Simulx project contains [a block \[COVARIATE\]](https://simulx.lixoft.com/definition/model/). A covariate element can be defined as an external file (csv, xlsx, xlsx, sas7bdat, xpt or txt) or as a data frame. In any case, it can contain columns occasions (optional), and it should contain one column per covariate (mandatory). Covariate names and categorical covariate values must correspond to covariates and categories defined in the model (block \[COVARIATE\]). The occasion headers must correspond to the occasion names defined in the occasion element. A data frame can be used only to define covariate elements of type 'common', i.e the same for all individuals (potentially occasion-wise). If you want to define subject-specific covariates, use an external file with an "id" column. Covariate definition with distributions is only possible in the GUI (in R, please sample from the desired distribution to generate an external file). An external file can be used in all cases (common or subject-specific). It can contain a column id (optional) in addition to occasions (optional), and should contain one column per covariate (mandatory). When id and occasion columns are present, then they must be the first columns. When the id column is not present, the covariate is considered common. If the project has a subject-specific occasion structure (defined with an external file with an ID column (see [`defineOccasionElement`](defineoccasionelement))), occasion-wise common elements are not allowed. Covariate elements must be either common over all subjects and all occasions, or can be defined with subject-specific occasion-wise values as an external table, with the same occasion structure. ### See also [`getCovariateElements`](getcovariateelements) ### Examples R if (FALSE) { defineCovariateElement(name = "name", element = "file/path") defineCovariateElement(name = "name", element = list(file = "file/path", sheet = "sheet_name")) defineCovariateElement(name = "name", element = data.frame(wt = 70, sex = 1, age = 35)) } # Create a manual and a subject-specific element initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "1.overview", "importFromMonolix_clinicalTrial.smlx") loadProject(project_name) defineCovariateElement(name = "wt_typical", element = data.frame(wt = 70, sex = 1, age = 35)) # manual samples <- data.frame(id = 1:10, sex = sample(0:1, 10, replace = TRUE), age = rnorm(10, 30, 10)) samples$wt <- rnorm(10, mean = ifelse(samples$sex == 0, 62, 75), sd = 10) # mean weight dependent on sex file_name <- tempfile("cov", fileext = ".csv") write.csv(samples, file_name, row.names = FALSE) defineCovariateElement(name = "wt_distribution", element = file_name) # subject-specific # Create an element with common occasions initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "3.definition", "3.7.occasions", "occasions_common.smlx") loadProject(project_name) defineCovariateElement(name = "Fasted_Fed", element = data.frame(occ = c(1, 2), FOOD = c("Fasted", "Fed"))) # Create an element with subject-specific occasions initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "3.definition", "3.7.occasions", "occasions_external.smlx") loadProject(project_name) occasions <- getOccasionElements() covariates <- data.frame(id = occasions$id, occ = unlist(occasions$occasions), FOOD = rep(c("Fasted", "Fed", "Fasted"), 9)) file_name <- tempfile("cov", fileext = ".csv") write.csv(covariates, file_name, row.names = FALSE) defineCovariateElement(name = "cov_external", element = file_name) Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # defineEndpoint ## \[Simulx\] Define endpoint element Define or edit an endpoint. Endpoints summarize the outcome values over all individuals, for each simulation group and each replicate. Endpoints are defined [as in Simulx GUI](https://simulx.lixoft.com/simulation/outcomes-and-endpoints/). ### Usage R defineEndpoint(name, element) ### Arguments name (character) (required) Endpoint name. element (list) (required) List with the endpoint settings:outcome (character or list) (required) - Outcome on which the endpoint is based on. If one outcome, use a string containing its name. Use list to combine outcomes with:names (vector of character) - Vector of outcome names groupName (character) - Name you want to give to the outcome combination operator (character) - Way in which output should be combined. One of "and"/"or" (in case of boolean outcomes) or "min"/"max" (in other cases) metric (character) (optional) - Calculation method for the endpoint. One of "arithmeticMean" (default), "geometricMean" or "median" if value-based outcome. In case of event-based outcomes, "kaplanMeier" (median survival) will be used and in case of boolean outcomes, "percentTrue" will be used by default. groupComparison (optional) (list) - Group comparison settings. List of:type (character) (optional) - one of "directComparison", "statisticalTest" (default). h1 (list) (optional) - a list containing hypothesis information:operator (character) - one of "!=" (default), "\>" or "\<", threshold (double) - a real number indicating the threshold for difference/oddsRatio (0 by default) pvalue (double) - a real number indicating the p-value (if type is "statisticalTest", 0.05 by default) ### Details To compute the defined endpoints, use [`runEndpoints`](runendpoints)` and get the results with `[`getEndpointsResults`](getendpointsresults). To specify if endpoints should be compared across simulation groups, use [`setGroupComparisonSettings`](setgroupcomparisonsettings). If group comparison is relevant, the way the comparison will be done for each endpoint (eg if statistical test and which p-value) should be defined in the endpoint element. ### See also [`getEndpoints`](getendpoints) ### Examples R # Endpoint with group comparison initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "6.outcome_endpoints", "6.1.outcome_endpoints", "OutcomeEndpoint_PKPD_changeFromBaseline.smlx") loadProject(project_name) defineEndpoint(name = "comparison", element = list(outcome = "changeFromBaseline", metric = "geometricMean", groupComparison = list(type = "statisticalTest", operator = "!=", threshold = "0", pvalue = 0.05))) # Combine multiple outcomes initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "6.outcome_endpoints", "6.1.outcome_endpoints", "OutcomeEndpoint_PDTTE_survival_NADIR_timeToNADIR.smlx") loadProject(project_name) defineEndpoint(name = "combined_endpoint", element = list(outcome = list(names = c("Survival", "TimeToNADIR_AsEvent"), groupName = "combined", operator = "min"))) Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # defineIndividualElement ## \[Simulx\] Define individual parameters element Define or edit an element of individual parameters. Individual parameter elements are defined and used for simulation [as in Simulx GUI](https://simulx.lixoft.com/definition/individual-parameters/). As for other elements, the individual parameters elements can be defined or imported, and they are saved with the Simulx project if calling [`saveProject`](saveproject)`. Once an individual parameters element is defined, it needs to be added to a simulation group with `[`setGroupElement`](setgroupelement) to be used in simulation. ### Usage R defineIndividualElement(name, element) ### Arguments name (character) Element name. element (character or dataFrame or list) Element definition from external file path or data frame with individual parameters as columns, or list to select the sheet of an excel file:file (character) Path to the individual parameters file. sheet (character) Name of the sheet in xlsx/xls file. ### Details Individual parameters to be defined depend on the model of the simulx project. If only the \[LONGITUDINAL\] block is present in the model: all parameters of the input list, except those defined as regressors. If both the \[LONGITUDINAL\] and \[INDIVIDUAL\] blocks are present: parameters defined in the DEFINITION section of the \[INDIVIDUAL\] block. An individual parameters element can be defined as an external file (csv, xlsx, xlsx, sas7bdat, xpt or txt) or as a data frame. In any case, it can contain columns occasions (optional), and it should contain one column per individual parameter (mandatory). The parameter headers must correspond to the parameter names used in the model. The occasion headers must correspond to the occasion names defined in the occasion element. A data frame can be used only to define individual parameter elements of type 'common', i.e the same for all individuals (but potentially occasion-wise). If you want to define subject-specific individual parameters, use an external file with an "id" column. An external file can be used in all cases (common or subject-specific). It can contain a column id (optional) in addition to occasions (optional), and should contain one column per parameter (mandatory). When id and occasion columns are present, then they must be the first columns. When the id column is not present, the parameter element is considered common. If the project has a subject-specific occasion structure (defined with an external file with an ID column (see [`defineOccasionElement`](defineoccasionelement))), occasion-wise common elements are not allowed. In this case, individual parameters elements have to be be either common over all subjects and all occasions, or can be defined with subject-specific occasion-wise values as an external table, with the same occasion structure. ### See also [`getIndividualElements`](getindividualelements) ### Examples R if (FALSE) { defineIndividualElement(name = "name", element = "file/path") defineIndividualElement(name = "name", element = list(file = "file/path", sheet = "sheet_name")) defineIndividualElement(name = "name", element = data.frame(Tlag = 0.5, ka = 0.25, V = 70, Cl = 12, F = 0.7)) } # Defining elements with one and multiple sets of indiv params initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "2.models", "longitudinal_individual.smlx") loadProject(project_name) defineIndividualElement(name = "custom_params", element = data.frame(Tlag = 0.5, ka = 0.25, V = 70, Cl = 12, F = 0.7)) # one set params <- data.frame(id = c(1, 2, 3), Tlag = 0.5, ka = 0.25, V = 70, Cl = 12, F = c(0.6, 0.7, 0.8)) file_name <- tempfile("cov", fileext = ".csv") write.csv(params, file_name, row.names = FALSE) defineIndividualElement(name = "different_F", element = file_name) # multiple sets # Common occasions initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "3.definition", "3.7.occasions", "occasions_common.smlx") loadProject(project_name) defineIndividualElement(name = "params_per_occ", element = data.frame(occ = c(1, 2), ka = c(0.2, 0.4), V = 10, Cl = 5)) # Subject-specific occasions initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "3.definition", "3.7.occasions", "occasions_external.smlx") loadProject(project_name) occasions <- getOccasionElements() params <- data.frame(ID = occasions$id, occ = unlist(occasions$occasions)) params$ka <- rlnorm(27, log(0.2), 0.1 + 0.1) params$V <- rlnorm(27, log(5), 0.2) params$Cl <- rlnorm(27, log(0.3), 0.15) file_name <- tempfile("cov", fileext = ".csv") write.csv(params, file_name, row.names = FALSE) defineIndividualElement(name = "params_external", element = file_name) Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # defineOccasionElement ## \[Simulx\] Define occasion element Define the occasion structure of the project. ### Usage R defineOccasionElement(element) ### Arguments element (character or dataFrame or list) Element definition from external file path or data frame with time and occasion levels as columns, or list to select the sheet of an excel file:file (character) Path to the occasion file. sheet (character) Name of the sheet in xlsx/xls file. ### Details The occasion structure of a project is defined and used for simulation [as in Simulx GUI](https://simulx.lixoft.com/definition/occasions/). The occasion element impacts the definition of other elements and the simulation. As for other elements, the occasion element can be defined or imported, and it is saved with the Simulx project if calling [`saveProject`](saveproject). If can be defined as an external file (csv, xlsx, xlsx, sas7bdat, xpt or txt) or as a data.frame. In any case, it should contain a column time and one column per occasion level. The header names for these occasion levels are free and used by Simulx. A data frame can be used only to define a common structure of occasions applied to all simulated subjects. An external file can be used in all cases (common or subject-specific structure). It can contain a column id (optional) in addition to other columns. When the id column is not present, the occasion structure is considered common. After defining a subject-specific occasion structure, other elements (parameters, covariates, treatments, outputs and regressors) must be either common over all subjects and all occasions, or can be defined with subject-specific occasion-wise values as an external table, with the same occasion structure. ### See also [`getOccasionElements`](getoccasionelements) ### Examples R if (FALSE) { defineOccasionElement(element = "file/path") defineOccasionElement(element = list(file = "file/path", sheet = "sheet_name")) defineOccasionElement(element = data.frame(time = c(0, 0.5, 2), occ1 = c(1, 1, 2), occ2 = c(1, 2, 3))) } # Example: define subject-specific occasions through external file initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "5.simulation", "replicates.smlx") loadProject(project_name) occasions <- data.frame(id = c(1, 1, 2, 2), time = c(0, 24, 0, 36), occ = c(1, 2, 1, 2)) file_name <- tempfile("cov", fileext = ".csv") write.csv(occasions, file_name, row.names = FALSE) defineOccasionElement(element = file_name) # Example: define common occasions through data.frame initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "5.simulation", "replicates.smlx") loadProject(project_name) defineOccasionElement(element = data.frame(time = c(0, 0.5, 2), occ1 = c(1, 1, 2), occ2 = c(1, 2, 3))) Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # defineOutcome ## \[Simulx\] Define outcome element Define or edit an outcome. Outcomes represent a post-processing of the simulation outputs done for each individual. Outcomes are defined [as in Simulx GUI](https://simulx.lixoft.com/simulation/outcomes-and-endpoints/). Outcomes can only be defined by the user (no outcome is imported), and they are saved with the Simulx project if calling [`saveProject`](saveproject). Contrary to the GUI, defining an outcome with the connectors does not automatically generate an endpoint. To compute the defined outcomes, use them in endpoints with [`defineEndpoint`](defineendpoint), run [`runEndpoints`](runendpoints) and get the results with [`getEndpointsResults`](getendpointsresults). ### Usage R defineOutcome(name, element) ### Arguments name (character) Outcome name. element (list) List with the outcome settings:output(character)Name of the output element on which the outcome is based.perOccasion(logical)(optional) If occasions are present in the simulation, it indicates if the outcome is calculated for each id or each occasion of each id. TRUE by default.Then, if the outcome is based on a continuous or categorical output:relativeTo(list)(optional) List of elements that define reference settings: reference - one of "baseline", "min", "max", "minCurrentTime", "maxCurrentTime" or "customValue", type - "ratio" or "difference", value - if reference is "customValue". processing(list)(required) List of elements that define how the output will be processed: operator (required) - one of "avg", "min", "max", "first", "last", "durationBelow", "durationAbove", "durationBetween", "timePoint", or "none" (if the output has a single time point and does not require processing). type (optional)if operator is "min" or "max", one of "value" (default), "timeContinuous" or "timeEvent", if operator is "durationBelow", "durationAbove" or "durationBetween", one of "cumulativeTime" (default), "percentTime", "nbObs", "firstOccurenceContinuous" or "firstOccurenceEvent", value (optional) - vector of boundaries if operator is "durationBelow", "durationAbove" or "durationBetween", or time point if operator is "timePoint" (0 by default). applyThreshold(list)(optional) List of elements: operator - one of "==", "!=", "\>=", "\>", "\<=" or "\<", value - a real number indicating the threshold value. Or if the outcome is based on a TTE output:event(list)(required) List of arguments that define event settings: type (required) - one of "timeOfEvents" (in case of single and repeated TTE), "hasAnEvent", "hasNoEvent" (in case of single TTE) or "numberOf" (in case of repeated TTE), rank (optional) - rank of the event of which time is the outcome (if type is "timeOfEvents" and repeated TTE). ### See also [`getOutcomes`](getoutcomes) ### Examples R # Define an outcome to calculate Cmax initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "6.outcome_endpoints", "6.1.outcome_endpoints", "OutcomeEndpoint_PKPD_Cmax_targetInhibition.smlx") loadProject(project_name) defineOutcome(name = "Cmax_outcome", element = list(output = "Plasma_concentration", processing = list(operator = "max", type = "value"))) # Define time above MIC as percentage initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "6.outcome_endpoints", "6.1.outcome_endpoints", "OutcomesEndpoints_antibiotics_TaboveMIC.smlx") loadProject(project_name) defineOutcome(name = "TimeAboveMIC", element = list(output = "mlx_Cc", processing = list(operator = "durationAbove", type = "percentTime", value = 0.5))) # Relative outcome with threshold per individual and occasion initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "6.outcome_endpoints", "6.1.outcome_endpoints", "OutcomeEndpoint_PKPD_LastPerOccasion.smlx") loadProject(project_name) defineOutcome(name = "custom_outcome", element = list(output = "prediction_Cc_per_id", perOccasion = TRUE, relativeTo = list(reference = "customValue", type = "difference", "value" = 2), processing = list(operator = "avg"), applyThreshold = list(operator = "<=", threshold = 2))) # Event based outcome initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "6.outcome_endpoints", "6.1.outcome_endpoints", "OutcomeEndpoint_PDTTE_survival_NADIR_timeToNADIR.smlx") loadProject(project_name) defineOutcome(name = "survival_outcome", element = list(output = "Death", event = list(type = "timeOfEvents"))) Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # defineOutputElement ## \[Simulx\] Define output element Define or edit an output element. ### Usage R defineOutputElement(name, element) ### Arguments name (character) Element name. element List with:data (character or dataFrame or list): data frame or path to external file (csv, xlsx, xlsx, sas7bdat, xpt or txt), or list to select the sheet of an excel file:file (character) Path to the output file. sheet (character) Name of the sheet in xlsx/xls file. output (character): name of any variable from the \[LONGITUDINAL\] block of the model (variable in EQUATION, smooth prediction listed under OUTPUT or noisy observation listed under DEFINITION). ### Details Output elements are defined and used for simulation [as in Simulx GUI](https://simulx.lixoft.com/definition/outputs/). As for other elements, the output elements can be defined or imported, and they are saved with the Simulx project if calling [`saveProject`](saveproject)`. Once an output element is defined, it needs to be added to a simulation group with `[`setGroupElement`](setgroupelement) to be used in simulation. To define an output element, in addition to the element name, you need to provide the time grid for simulation and the output to simulate. The field data can be specified with a data frame or an external file (csv, xlsx, xlsx, sas7bdat, xpt or txt). To define a regular output, common to all individuals, you can use a data frame, with column headers start, interval and final. All time points from "start" to "final" by steps of "interval" will be used for simulation. If the project has a common occasion structure, this data frame can contain a column occasion and several lines to define the regular treatment occasion-wise. To define an output by giving a specific list of times, both data frames and external files (csv, xlsx, xlsx, sas7bdat, xpt or txt) can be used, with a column time. They can contain columns occasions (optional). The occasion headers must correspond to the occasion names defined in the occasion element. Data frames can be used only to define output elements of type 'common', i.e the same for all individuals (potentially occasion-wise). If you want to define subject-specific output elements, you have to use an external file with an "id" column. An external file can be used in all cases (common or subject-specific). It can contain a column id (optional) in addition to occasions (optional), and should contain one column time (mandatory). When id and occasion columns are present, then they must be the first columns. When the id column is not present, the covariate is considered common. If the project has a subject-specific occasion structure (defined with an external file with an ID column (see [`defineOccasionElement`](defineoccasionelement))), occasion-wise common elements are not allowed. Output elements must be either common over all subjects and all occasions, or can be defined with subject-specific occasion-wise values as an external table, with the same occasion structure. ### See also [`getOutputElements`](getoutputelements) ### Examples R if (FALSE) { defineOutputElement(name = "name", element = list(data ="file/path")) defineOutputElement(name = "name", element = list(data = list(file = "file/path", sheet = "sheet_name"))) defineOutputElement(name = "name", element = list(data = data.frame(time = 24), output = "AUC")) } # Define subject-specific outputs using an external file (saved in tmp directory) initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "5.simulation", "replicates.smlx") loadProject(project_name) samplings <- data.frame(id = c(1, 1, 2, 2, 3, 3), time = c(0, 24, 0, 48, 0, 72)) file_name <- tempfile("cov", fileext = ".csv") write.csv(samplings, file_name, row.names = FALSE) defineOutputElement(name = "external_output", element = list(data = file_name, output = "Cc")) # Define manual and regular output initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "1.overview", "importFromMonolix_clinicalTrial.smlx") loadProject(project_name) defineOutputElement(name = "AUC_0_24", element = list(data = data.frame(time = 24), output = "AUC")) defineOutputElement(name = "Cc_7days", element = list(data = data.frame(start = 0, interval = 1, final = 168), output = "Cc")) # Define manual and regular occasion-wise output initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "3.definition", "3.7.occasions", "occasions_two_levels.smlx") loadProject(project_name) defineOutputElement(name = "manualOcc", element = list(data = data.frame(time = c(0, 2, 24, 4, 36), occ1 = c(1, 1, 1, 2, 2), occ2 = c(1, 1, 2, 1, 2)), output = "Y")) defineOutputElement(name = "regularOcc", element = list(data = data.frame(start = c(0, 0, 24, 24), interval = c(1, 2, 1, 2), final = c(24, 48, 48, 72), occ1 = c(1, 2, 1, 2), occ2 = c(1, 1, 2, 2)), output = "Cc")) Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # definePopulationElement ## \[Simulx\] Define population element Define or edit an element of population parameters. Population parameter elements are defined and used for simulation [as in Simulx GUI](https://simulx.lixoft.com/definition/population-parameters/). As for other elements, the population parameters elements can be defined or imported, and they are saved with the Simulx project if calling [`saveProject`](saveproject). Once a population parameters element is defined, it needs to be added to a simulation group with [`setGroupElement`](setgroupelement) to be used in simulation. ### Usage R definePopulationElement(name, element) ### Arguments name (character) Element name. element (character or dataFrame or list) Element definition from external file path or data frame with population parameters as columns, or list to select the sheet of an excel file:file (character) Path to the population file. sheet (character) Name of the sheet in xlsx/xls file. ### Details Definition of population parameters as simulation elements allows to simulate individual parameters from probability distributions. It is possible only if the model includes [a block \[INDIVIDUAL\]](https://simulx.lixoft.com/definition/model/) to consider the inter-individual variability. A population parameters element can be defined as an external file (csv, xlsx, xlsx, sas7bdat, xpt or txt) or as a data frame. In any case, it should contain one column per population parameter (mandatory). To check exactly which column names to use, you can use [`getPopulationElements`](getpopulationelements) and view the population parameters element that was automatically created after defining the model (if the model has an \[INDIVIDUAL\] block). To define a population parameters element with several lines, with several sets to be used in different replicate simulations, an external file is required. In this case, you should also set the number of replicates for your simulation with [`setNbReplicates`](setnbreplicates), otherwise only the first set of population parameters will be taken into account. Each replicate uses one set of population parameters with the order of the appearance in the table. Note: It is not possible to define population elements with distributions with the connectors as in the GUI. To do this, please sample values from distributions in R and create the element with different rows as in the last example below. ### See also [`getPopulationElements`](getpopulationelements), [`setNbReplicates`](setnbreplicates) ### Examples R if (FALSE) { definePopulationElement(name = "name", element = "file/path") definePopulationElement(name = "name", element = list(file = "file/path", sheet = "sheet_name")) definePopulationElement(name = "name", element = data.frame(Cl_pop = 1, V_pop = 2.5)) } # Define a pop param element with a data frame loadProject(file.path(getDemoPath(),"3.definition","3.3.population_parameters","pop_parameters_manual.smlx")) # get the existing pop param element as an example ExamplePopParamData = getPopulationElements()$manual_parameters$data ExamplePopParamData[] = rep(1,11) #set the desired values, for simplicity we use all param =1 definePopulationElement(name = "PopParam", element = ExamplePopParamData) # Check impact of varying ka with replicates (external file required) loadProject(file.path(getDemoPath(),"3.definition","3.3.population_parameters","pop_parameters_manual.smlx")) # get the existing pop param element as an example: ExamplePopParamData = getPopulationElements()$manual_parameters$data # add lines to the data frame to have different values for ka: PopParamReplicates = ExamplePopParamData[rep(1,10),] PopParamReplicates$ka_pop = rnorm(10,mean = 0.8, sd = 0.1) # write the table to an external file (required because it has several lines): file_name = tempfile("PopParamReplicates.csv") write.csv(PopParamReplicates, file_name, row.names = FALSE) # define the new element and use it in simulation: definePopulationElement(name = "PopParamReplicates", element = file_name) setGroupElement(group = "simulationGroup1", elements = "PopParamReplicates") setNbReplicates(nb = 10) # to simulate the project 10x, each time with a different population parameter element runSimulation() if (FALSE) getSimulationResults() Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # defineRegressorElement ## \[Simulx\] Define regressor element Define or edit a regressor element. Regressor elements are defined and used for simulation [as in Simulx GUI](https://simulx.lixoft.com/definition/regressors/). Once a regressor element is defined, it needs to be added to a simulation group with [`setGroupElement`](setgroupelement) to be used in simulation. Regressor elements can be defined only if regressors are defined in the model loaded in the Simulx project. ### Usage R defineRegressorElement(name, element) ### Arguments name (character) Element name. element (character or dataFrame or list) Element definition from external file path or data frame with time and regressors as columns, or list to select the sheet of an excel file:file (character) Path to the regressors file. sheet (character) Name of the sheet in xlsx/xls file. ### Details To define a regressor element, in addition to the element name, you need to provide in the field data the time points and values of the regressor at each time point. To simulate the model for points outside of this grid, last value carried forward interpolation is used. The field data can be specified with a data frame or an external file (csv, xlsx, xlsx, sas7bdat, xpt or txt). They should contain a column time and a column for each regressor variable. They can contain columns occasions (optional). The occasion headers must correspond to the occasion names defined in the occasion element. Data frames can be used only to define regressor elements of type 'common', i.e the same for all individuals (potentially occasion-wise). If you want to define subject-specific regressor elements, you have to use an external file with an additional "id" column. An external file can be used in all cases (common or subject-specific). It can contain a column id (optional) in addition to occasions (optional), time (mandatory) and one column per regressor defined in the model (mandatory). When id and occasion columns are present, then they must be the first columns. When the id column is not present, the covariate is considered common. If the project has a subject-specific occasion structure (defined with an external file with an ID column (see [`defineOccasionElement`](defineoccasionelement))), occasion-wise common elements are not allowed. In this case, regressors must be either common over all subjects and all occasions, or defined with subject-specific occasion-wise values as an external table, with the same occasion structure. Note: It is not possible to define regressor elements with distributions with the connectors as in the GUI. To do this, please sample values from distributions in R and create the element with different rows as in the examples below. ### See also [`getRegressorElements`](getregressorelements) ### Examples R if (FALSE) { defineRegressorElement(name = "name", element = "file/path") defineRegressorElement(name = "name", element = list(file = "file/path", sheet = "sheet_name")) defineRegressorElement(name = "name", element = data.frame(time = c(0, 0.5, 2), PCA = c(1, 2, 5.25))) } # Define subject-specific regressors using an external file (saved in tmp directory) initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "3.definition", "3.6.regressors", "regressor_manual_paramCovRelationship.smlx") loadProject(project_name) samplings <- data.frame(id = c(1, 1, 2, 2, 3, 3), time = c(0, 24, 0, 48, 0, 72), PCA = c(9, 15, 5, 20, 3, 14)) file_name <- tempfile("cov", fileext = ".csv") write.csv(samplings, file_name, row.names = FALSE) defineRegressorElement(name = "reg_external", element = file_name) # Define manual regressor element initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "3.definition", "3.6.regressors", "regressor_manual_paramCovRelationship.smlx") loadProject(project_name) defineRegressorElement(name = "reg_manual", element = data.frame(time = c(0, 0.5, 2), PCA = c(1, 2, 5.25))) # Define manual occasion-wise regressors initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "3.definition", "3.6.regressors", "regressor_manual_paramCovRelationship.smlx") loadProject(project_name) defineOccasionElement(element = data.frame(time = c(0, 0, 0, 0), occ1 = c(1, 2, 1, 2), occ2 = c(1, 1, 2, 2))) defineRegressorElement(name = "name", element = data.frame(time = c(0, 0.5, 2, 5, 6), PCA = c(1, 2, 5.25, 6, 7), occ1 = c(1, 1, 1, 2, 2), occ2 = c(1, 1, 2, 1, 2))) Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # defineTreatmentElement ## \[Simulx\] Define treatment element Define or edit a treatment element. ### Usage R defineTreatmentElement(name, element) ### Arguments name (character) Element name. element (list) List with the treatment settings:data (mandatory)(data frame OR path to external file (csv, xlsx, xlsx, sas7bdat, xpt or txt) OR list to select the sheet of an excel file: list( file = path_to_excel, sheet = sheet_name )Column headers:\[only data frame\] for a regular treatment common to all ids:occ (optional, if common occasion structure, same header as in the occasion element), start, interval, nbDoses, amount \[only data frame\] for a manual treatment common to all ids: occ (optional, if common occasion structure, same header as in the occasion element), time, amount, washout (optional, to add a washout just before the dose, otherwise washout = FALSE by default) \[only external file\] for a manual treatment common or specific to each id: id (optional), occ (optional), time, amount, washout (optional) \[data frame or external\] in case of infusion: tInf (duration) OR rate (mutually exclusive). admID (optional)(integer) same as integer in the model as administration typeprobaMissDose (optional)(double) probability to miss a dose (number in \[0 1\])repeats (optional)(vector) to repeat the specified treatment after a specific duration.Elements:cycleDuration duration after which the treatment will be repeated (can be longer than the treatment duration) NumberOfRepetitions number of times the treatment will be repeated scale (optional)(list) to scale the dose amount by a covariate. The scaled amount will be administered instead of amount.covariate (character) covariate to use for scaling (same name as in the \[COVARIATE\] block of the model) intercept (double, required for continuous covariate): intercept to use in the scaling formula: scaledAmount = amount\*cov + intercept modalities (list, required for categorical covariate): list of lists with, for each modality, the name of the modality (eg "male"), the factor, and the intercept to use in the scaling formula: scaledAmount = \[cov = modality\]\*factor\*amt + intercept (no scaling if factor =1 and intercept = 0) scaleDuration (optional, logical) if TRUE (default), infusion duration will be scaled, otherwise it will be rate ### Details Treatment elements are defined and used for simulation [as in Simulx GUI](https://simulx.lixoft.com/definition/treatments/). As for other elements, the treatment elements can be defined or created at import, and they are saved with the Simulx project if calling [`saveProject`](saveproject)`. Once an output element is defined, it needs to be added to a simulation group with `[`setGroupElement`](setgroupelement) to be used in simulation. Several treatment elements can be added to the same simulation group and they will be both administered to every individual in the group. To define a treatment element, in addition to the element name, you need to provide a list with at list one field "data" containing the dose information. The field data can be specified with a data frame or an external file (csv, xlsx, xlsx, sas7bdat, xpt or txt). To define a regular treatment, common to all individuals, you can use a data frame, with column headers start, interval, nbDoses and amount. You can include an optional column tInf or rate to define an infusion. If the project has a common occasion structure (i.e. same occasions for all individuals), this data frame can contain a column occasion to define the regular treatment occasion-wise. The occasion headers must correspond to the occasion names defined in the occasion element. To define a treatment by giving a specific list of times and amounts, both data frames and external files (csv, xlsx, xlsx, sas7bdat, xpt or txt) can be used, with a column time. They can contain columns id and occasions (optional). The occasion headers must correspond to the occasion names defined in the occasion element. Data frames can be used only to define output elements of type 'common', i.e the same for all individuals (potentially occasion-wise). If you want to define subject-specific treatment elements, you have to use an external file with an "id" column. An external file can be used in all cases (common or subject-specific). It can contain a column id (optional) in addition to occasions (optional), and should contain one column time (mandatory) and one column amount (mandatory). When id and occasion columns are present, then they must be the first columns. When the id column is not present, the covariate is considered common. ### Note To define a regular schedule, it is advised to use a regular treatment without repeats, rather than a manual treatment with repeats. Repeats are useful to create more complex schedules in addition to a manual or regular definition, such as dosing regimen 3 weeks ON, 1 week OFF. To see the impact of a treatment until the end of a dosing regimen, you should set an output element that spans the duration of the treatment to the same simulation group. ### See also [`getTreatmentElements`](gettreatmentelements) ### Examples R if (FALSE) { defineTreatmentElement(name = "name", element = list(data = "file/path")) defineTreatmentElement(name = "name", element = list(data = list(file = "file/path", sheet = "sheetname"))) defineTreatmentElement(name = "name", element = list(probaMissDose=0, admID=1, repeats=c(cycleDuration = 1, NumberOfRepetitions=12), data=data.frame(time=c(1,2), amount=c(1,2), tInf=c(0, 1), washout=c(TRUE, FALSE)))) defineTreatmentElement(name = "name", element = list(probaMissDose=0, admID=1, repeats=c(cycleDuration = 1, NumberOfRepetitions=12), data=data.frame(time=c(10,10,10,10), amount=c(10,20,30,40), occ1=c(1,1,2,2), occ2=c(1,2,1,2)))) defineTreatmentElement(name = "name", element = list(probaMissDose=0, admID=1, repeats=c(cycleDuration = 1, NumberOfRepetitions=12), data=data.frame(start=1, interval=2, nbDoses=10, amount=1))) defineTreatmentElement(name = "name", element = list(admID=1, scale=list(covariate="age", intercept=12), data=data.frame(start=1, interval=2, nbDoses=10, amount=1))) defineTreatmentElement(name = "name", element = list(admID=1, scale=list(covariate="sex", modalities=list(list(name="0", factor=1, intercept=10), list(name="1", factor=1.5, intercept=10))), data=data.frame(start=1, interval=2, nbDoses=10, amount=1))) } ##### Working example with treatment combinations ##### # In this demo, the first group receives only the chemotherapy, while the second group receives both the chemotherapy and the anti-angiogenic therapy. # Note that the chemotherapy treatment uses adm=1 to be applied to compartment 1 via the macro iv(adm=1, cmt=1) in the model representing the concentration of the chemo drug. # The anti-angiogenic treatment is defined with adm=2 which is applied via the macro iv(adm=2, cmt=2) to compartment 2 representing the concentration of anti-angiogenic drug. initializeLixoftConnectors("simulx") loadProject(paste0(getDemoPath(), "/3.definition/3.1.treatments/treatment_combinations.smlx")) # to see how the structural model is defined: file.show(getStructuralModel()) defineTreatmentElement(name = "Chemotherapy", element = list(data=data.frame(start=10, interval=14, nbDoses=10, amount=1))) defineTreatmentElement(name = "AntiAngionenic_treatment", element = list(admID = 2, data=data.frame(start=10, interval=7, nbDoses=20, amount=1))) setGroupElement("simulationGroup1","Chemotherapy") setGroupElement("simulationGroup2",c("Chemotherapy","AntiAngionenic_treatment")) runSimulation() # use ggplot or export to Monolix/PKanalix to plot trajectories exportProject(settings = list(targetSoftware = "monolix"),force = TRUE) if (FALSE) { plotObservedData( settings = list(dots = FALSE, ylab = "Target Occupancy", legend = TRUE), stratify = list(colorGroup = list(name = "group")), preferences = list(obs = list(lineWidth = 0.5))) } ##### Working example with a treatment scaled by weight and based on genotype ##### # In this demo, a weight-based dose is defined by indicating the dose per unit weight in the amount box (14 nmol/kg) and using the "Scale amount by a covariate" option with "Weight" selected as covariate. # The "intercept" could be used to define a offset common to all weights (e.g 14nmol/kg + 10nmol). # When an infusion duration or rate has been defined, the user can choose if the infusion duration or the infusion rate is scaled by the covariate. # For categorical covariates, such as the genotype, a scaling factor and an intercept can be defined for each category. # In this demo, the scaling for Homozygous is 1 meaning that they receive the dose defined in the amount box. # For heterozygous, the scaling is 0.8, meaning that they receive 0.8 times the amount in the amount box. initializeLixoftConnectors("simulx") loadProject(paste0(getDemoPath(), "/3.definition/3.1.treatments/treatment_weight_and_genotype_based.smlx")) defineTreatmentElement(name = "14nmolPerKg", element = list(data=data.frame(start=0, interval=21, nbDoses=5, amount=14, tInf = 0.208), scale=list(covariate="Weight", intercept = 0, scaleDuration = FALSE))) defineTreatmentElement(name = "1000nmol", element = list(data=data.frame(start=0, interval=21, nbDoses=5, amount=1000, tInf = 0.208))) defineTreatmentElement(name = "1000nmolHomo_800nmolHetero", element = list(data=data.frame(start=0, interval=21, nbDoses=5, amount=1000, tInf = 0.208), scale=list(covariate="Genotype", modalities=list(list(name="Homozygous", factor=1, intercept=0), list(name="Heterozygous", factor=0.8, intercept=0)), scaleDuration = FALSE))) setGroupElement("Weight_based","14nmolPerKg") setGroupSize("Weight_based",20) setGroupElement("Flat_dose","1000nmol") setGroupSize("Flat_dose",20) setGroupElement("Genotype_based","1000nmolHomo_800nmolHetero") setGroupSize("Genotype_based",20) runSimulation() # use ggplot or export to Monolix/PKanalix to plot trajectories exportProject(settings = list(targetSoftware = "monolix"),force = TRUE) plotObservedData(obsName = "yTO", settings = list(dots = FALSE, ylab = "Target Occupancy"), stratify = list(splitGroup = list(name = "group")), preferences = list(obs = list(lineWidth = 0.5))) #> Warning: No shared levels found between `names(values)` of the manual scale and the #> data's linetype values. ![defineTreatmentElement-1.png](https://monolixsuite.slp-software.com/__attachments/a_0f9dacb5983e62c7bdbd1c0b12172bce8982c6d82bd03d37f72321a7d28dc71b/defineTreatmentElement-1.png?cb=e36de9ce1694531d696da4e30c8055f5) R ##### Working example with a probability to miss a dose ##### # In this demo, the second treatment is defined with a probability to miss a dose of 0.1, meaning that on average 10% of the doses will not be taken. The missed doses are random. initializeLixoftConnectors("simulx") loadProject(paste0(getDemoPath(), "/3.definition/3.1.treatments/treatment_non_adherence.smlx")) defineTreatmentElement(name = "OncePerDay_full_compliance", element = list(data=data.frame(start=0, interval=1, nbDoses=112, amount=100))) defineTreatmentElement(name = "OncePerDay_partial_compliance", element = list(data=data.frame(start=0, interval=1, nbDoses=112, amount=100),probaMissDose = 0.1)) setGroupElement("simulationGroup1","OncePerDay_full_compliance") renameGroup("simulationGroup1","FullCompliance") setGroupElement("simulationGroup2","OncePerDay_partial_compliance") renameGroup("simulationGroup2","NonAdherence") setGroupSize("FullCompliance",20) setGroupSize("NonAdherence",20) runSimulation() # use ggplot or export to Monolix/PKanalix to plot trajectories exportProject(settings = list(targetSoftware = "monolix"),force = TRUE) plotObservedData( settings = list(dots = FALSE, ylab = "Target Occupancy"), stratify = list(splitGroup = list(name = "group")), preferences = list(obs = list(lineWidth = 0.5))) #> Warning: No shared levels found between `names(values)` of the manual scale and the #> data's linetype values. ![defineTreatmentElement-2.png](https://monolixsuite.slp-software.com/__attachments/a_42c67dcf9b5f4370d8144852c6c1ba869fbae65b44d1163082934b31749d908c/defineTreatmentElement-2.png?cb=aeb7c6063c505c5e8942409792bcb70a) R ##### Working example with an external file ##### # Demo: use an external file to define a dose by age group: 1-2 years 12.5 mg, 3-6 years 18.75 mg and 7-15 years 25 mg. # The age also appears as covariate in the model and the covariate element is defined via an external file. # To make sure the covariates are sampled from the covariate external file and the doses sampled from the treatment external file are consistent (i.e correspond to the same id and thus the same age), the option "shared id" is selected between covariate and treatment elements. initializeLixoftConnectors("simulx") loadProject(paste0(getDemoPath(), "/3.definition/3.1.treatments/treatment_external_byAgeGroup.smlx")) tableAge = getCovariateElements()$External_AGE_values$data AmtByAgeGroups = (tableAge$AGE < 3)*12.5 + ((tableAge$AGE >=3) & (tableAge$AGE < 7))*18.75 + (tableAge$AGE >= 7)*25 Nid = length(AmtByAgeGroups) dataAmtByAgeGroups = data.frame(id = tableAge$ID, time = rep(0,Nid), amount = AmtByAgeGroups) file_name <- tempfile("trt", fileext = ".csv") write.csv(dataAmtByAgeGroups, file_name, row.names = FALSE) defineTreatmentElement(name = "doseByAgeGroup", element = list(data = file_name)) setGroupElement("simulationGroup1",c("doseByAgeGroup","External_AGE_values","regularCc")) setSharedIds(c("covariate", "treatment")) runSimulation() # use ggplot or export to Monolix/PKanalix to plot trajectories exportProject(settings = list(targetSoftware = "monolix"),force = TRUE) plotObservedData(settings = list(dots = FALSE, ylab = "Cc",legend = TRUE, ylim = c(0,13)), stratify = list(splitGroup = list(name = "AGE", breaks = c(2,7)), colorGroup = list(name = "ID")), preferences = list(obs = list(lineWidth = 0.5))) #> [ERROR] - Invalid stratify inputs: #> - Invalid covariate input for color: 'ID' #> [ERROR] Error during observed data (outputplot) chart creation: #> - Invalid stratify inputs: #> - Invalid covariate input for color: 'ID' ##### Working example with washout ##### # In this demo, two different formulations are given. # The reference formulation is given at time zero. # The test formulation is given after a long washout period. # In order not to simulate this washout period, the test dose is defined at time 48 and a washout is applied just before the test dose to reset the model to its initial state. initializeLixoftConnectors("simulx") loadProject(paste0(getDemoPath(), "/3.definition/3.1.treatments/treatment_washout.smlx")) defineTreatmentElement(name = "ReferenceFormulation_atTime0", element = list(data=data.frame(time=0, amount=600))) defineTreatmentElement(name = "TestFormulation_atTime48", element = list(admID = 2, data=data.frame(time=0, amount=600, washout = TRUE))) setGroupElement("simulationGroup1",c("ReferenceFormulation_atTime0","ReferenceFormulation_atTime0")) ##### Working example with a regular treatment and repeats ##### # The "repeat" option allows to repeat a base pattern several times with a defined periodicity. # In this demo, the first treatment is defined as one dose per day during 112 days. # The second treatment is defined as one dose per day during 14 days and this is repeated every 28 days leading to a 2 weeks on / 2 weeks off pattern. initializeLixoftConnectors("simulx") loadProject(paste0(getDemoPath(), "/3.definition/3.1.treatments/treatment_regular_cycles.smlx")) defineTreatmentElement(name = "OncePerDay_4weeksOn", element = list(data=data.frame(start=0, interval=1, nbDoses=112, amount=100))) defineTreatmentElement(name = "OncePerDay_2weeksOn2weeksOff", element = list(repeats=c(cycleDuration = 28, NumberOfRepetitions=4), data=data.frame(start=0, interval=1, nbDoses=14, amount=100))) setGroupElement("simulationGroup1","OncePerDay_4weeksOn") renameGroup("simulationGroup1","4weeksOn") setGroupElement("simulationGroup2","OncePerDay_2weeksOn2weeksOff") renameGroup("simulationGroup2","2weeksOn2weeksOff") runSimulation() # use ggplot or export to Monolix/PKanalix to plot trajectories exportProject(settings = list(targetSoftware = "monolix"),force = TRUE) plotObservedData(settings = list(dots = FALSE, ylab = "Cc",legend = TRUE), stratify = list(color = list(name = "group")), preferences = list(obs = list(lineWidth = 0.5))) #> [ERROR] Unexpected type encountered for 'state' : 'color' field must be a vector of string. ![defineTreatmentElement-3.png](https://monolixsuite.slp-software.com/__attachments/a_9265bcd5f631d3140e70d91016ad544d6ae21d889c3ce97e0b214c2ac52240e3/defineTreatmentElement-3.png?cb=d7becdd527c77e8a8216b001bb36867b) R ``` ``` Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # deleteAdditionalCovariate ## \[Monolix - PKanalix\] Delete additional covariate Delete a created additinal covariate. ### Usage R deleteAdditionalCovariate(name) ### Arguments name (character) name of the covariate. ### See also [`addAdditionalCovariate`](addadditionalcovariate) ### Examples Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # deleteCustomNCAParameter ## \[PKanalix\] Delete a custom NCA parameter Remove a previously created custom parameter from the current project. Available arguments: |--------|-------------------------|------------------------| | "name" | (*character*, required) | Name of the parameter. | ### Usage R deleteCustomNCAParameter(name) ### See also [`createCustomNCAParameter`](createcustomncaparameter)`, `[`getCustomNCAParameters`](getcustomncaparameters)`, `[`addCustomNCAParametersFromPreferences`](addcustomncaparametersfrompreferences) Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # deleteElement ## \[Simulx\] Delete element Delete an element of any type. ### Usage R deleteElement(name) ### Arguments name (character) Element name. ### Details Elements defined are created in the background and saved with the Simulx project if calling [`saveProject`](saveproject). To check which elements of a certain type have been defined so far, please use one of the "get..Elements" connectors: [`getCovariateElements`](getcovariateelements), [`getPopulationElements`](getpopulationelements), [`getIndividualElements`](getindividualelements), [`getTreatmentElements`](gettreatmentelements), [`getOccasionElements`](getoccasionelements), [`getRegressorElements`](getregressorelements). Elements cannot be deleted if they are used for the simulation. To remove an element from the simulation, use [`removeGroupElement`](removegroupelement). ### Examples R initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "1.overview", "importFromMonolix_clinicalTrial.smlx") loadProject(project_name) deleteElement(name = "mlx_CovDist") Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # deleteEndpoint ## \[Simulx\] Delete an endpoint Delete an endpoint. ### Usage R deleteEndpoint(name) ### Arguments name (character) Endpoint name ### Details Endopints defined are created in the background and saved with the Simulx project if calling [`saveProject`](saveproject). To check which endpoints have been defined, please use [`getEndpoints`](getendpoints). ### See also [`deleteOutcome`](deleteoutcome) ### Examples R initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "6.outcome_endpoints", "6.1.outcome_endpoints", "OutcomeEndpoint_PDTTE_survival_NADIR_timeToNADIR.smlx") loadProject(project_name) deleteEndpoint("mean_NADIR") Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # deleteFilter ## \[Monolix - PKanalix\] Delete filter Delete a data set. Only filtered data set which are not active and whose children are not active either can be deleted. ### Usage R deleteFilter(name) ### Arguments name (character) data set name. ### See also [`createFilter`](createfilter) ### Examples R if (FALSE) { deleteFilter(name = "filter2") } Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # deleteNCARatio ## \[PKanalix\] Delete an NCA ratio. \[PKanalix\] Delete an NCA ratio. ### Usage R deleteNCARatio(name) ### Arguments name (character) Name of the parameter defined as ratio of existing NCA parameters. ### See also [`createNCARatio`](createncaratio)`, `[`getNCARatios`](getncaratios) ### Examples R initializeLixoftConnectors("pkanalix") loadProject(file.path(getDemoPath(), "1.basic_examples", "project_accumulationRatio.pkx")) deleteNCARatio(name = "accumulationRatio") getNCARatios() #> named list() Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # deleteOccasionElement ## \[Simulx\] Delete occasion element Delete the occasion element. ### Usage R deleteOccasionElement() ### Details The occasion element impacts the definition of other elements and the simulation. As for other elements, the occasion element can be defined or imported, and it is saved with the Simulx project if calling [`saveProject`](saveproject). To check if an occasion element has been defined, please use [`getOccasionElements`](getoccasionelements). The occasion element may impact the definition of other elements. When deleting the occasion element, all other elements that depend on occasions are also deleted. ### Examples R initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "3.definition", "3.7.occasions", "occasions_common.smlx") loadProject(project_name) deleteOccasionElement() Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # deleteOutcome ## \[Simulx\] Delete an outcome Delete an outcome. ### Usage R deleteOutcome(name) ### Arguments name (character) Outcome name ### Details Outcomes defined are created in the background and saved with the Simulx project if calling [`saveProject`](saveproject). To check which outcomes have been defined, please use [`getOutcomes`](getoutcomes). An outcome used in an endpoint cannot be deleted. The related endpoint must be deleted first with [`deleteEndpoint`](deleteendpoint). ### See also [`deleteEndpoint`](deleteendpoint) ### Examples R initializeLixoftConnectors("simulx") project_name <- file.path(getDemoPath(), "6.outcome_endpoints", "6.1.outcome_endpoints", "OutcomeEndpoint_PDTTE_survival_NADIR_timeToNADIR.smlx") loadProject(project_name) deleteEndpoint("mean_NADIR") deleteOutcome("NADIR") Last updated: August 08, 2024 --- version: "2024R1" language: "en" --- # Deprecated packages * [Package "mlxR"](https://monolixsuite.slp-software.com/r-functions/2024R1/package-mlxr.md) Last updated: February 19, 2025 --- version: "2024R1" language: "en" --- # Documentation of older versions Here you can download the documentation for the older versions of lixoftConnectors: * [lixoftConnectors_2023R1.pdf](https://monolixsuite.slp-software.com/__attachments/a_90c14f5fadd3a2ababb425abba2e64224d83525d8ef0a9ae946b50f15d8f34ba/lixoftConnectors_2023R1.pdf.md?cb=d5793b54edd430ef32d7f3ba7ced23d7) * [lixoftConnectors_2021R2.pdf](https://monolixsuite.slp-software.com/__attachments/a_bd7616c823f65b6243ba6f8dca9c9397950cdccf8ff71420bddfdbde570066fa/lixoftConnectors_2021R2.pdf.md?cb=460d1b6079b46f6ed56535ebb7472772) * [lixoftConnectors_2021R1.pdf](https://monolixsuite.slp-software.com/__attachments/a_497c039c898986623df1fc0305624b5185eb5c5f2117ecbb67388c3a4717c510/lixoftConnectors_2021R1.pdf.md?cb=3395c2759dd90398e5250108463e3694) * [lixoftConnectors_2020R1.pdf](https://monolixsuite.slp-software.com/__attachments/a_f8629b9aebd6f915b9c6d7808905333951aee0c066940339fada4a8832649497/lixoftConnectors_2020R1.pdf.md?cb=2b834cd746aa92d133bdd62fd7998eb1) * [lixoftConnectors_2019R2.pdf](https://monolixsuite.slp-software.com/__attachments/a_00ac952b5c99aa2f1d3539963751e15e4ab6682fa85710864a4b9b3afa9403c1/lixoftConnectors_2019R2.pdf.md?cb=d3a4aea539a3e5d153d9e8531bed494c) Last updated: October 11, 2024 --- version: "2024R1" language: "en" --- # Dose linearity Dose linearity can be described when the systemic exposure (e.g. Cmax, AUC) is a straight-line function of the administered dose amount. But unlike dose proportionality, that line need not pass through the origin. It is a **weaker condition than** [**dose proportionality**](https://monolixsuite.slp-software.com/r-functions/2024R1/dose-proportionality-assessment): it permits a constant, dose-independent offset/intercept in exposure, so a fold change in dose does not necessarily produce the same fold change in exposure. This offset is occasionally relevant in early pharmacokinetic analysis, e.g. for **endogenous compounds** , where measurable exposure is present even in the absence of drug. The simple **linear-model** assessment described here **complements the power-model dose-proportionality** analysis: the power model characterises the *shape* of the dose--exposure relationship, while the **linear model isolates its** ***intercept*** \[[Vogel HG, Maas J, Gebauer A, editors. *Drug Discovery and Evaluation: Methods in Clinical Pharmacology.* Berlin, Heidelberg: Springer; 2011.](https://e-library.nu.edu.sd/book/76/view)\]. ## Dose linearity linear pharamcokinetics Dose linearity refers to the *shape of the* ***dose--exposure relationship*** , not to linear first-order pharmacokinetics (where rate of elimination ke concentration). A compound can show dose-linear exposure without its clearance and volume being strictly dose-independent and vice versa. The two terms are easily confused but describe different things. ## Concept and mathematics The assessment fits a **linear model** by ordinary least squares, which relates a PK parameter to dose through and focuses on the **intercept** . Dose proportionality is the special case . The whole assessment reduces to a single question: *Is the intercept* *different from zero?* This can be answered by a linear regression model along the intercept **p-value** : if were really 0, how often would chance alone produce an estimate this far from 0? A **large** p-value means the intercept cannot be told apart from 0; a **small** one (below the significance level, here 10%) means it is judged to differ from 0. * **and p-value** ***not*** **significant** : the intercept cannot be told apart from 0. The line passes through the origin → ***consistent with dose proportionality***. * **and p-value significant:** real positive offset linear but not proportional; plausible for **endogenous compounds** or an assay baseline. * **and p-value signifcant:** impossible negative exposure at dose 0 → red flag the linear assumption is wrong The **90% confidence interval** of is an equivalent check and gives the same insight: it **contains 0** when the p-value is large and **excludes 0** when it is small. Only the interval's position relative to **0** is informative. **Note:** ++Linearity is assumed, not tested.++ Fitting a line does not verify that the relationship *is* a line. A saturating curve can still yield a respectable and a plausible-looking intercept. ## Function `doseLinearity()` The above explaned methodology is implemented as a single function that takes the data for one NCA parameter, retrieved by `getNCAIndividualParameters()$parameters`,and returns the assessment: R doseLinearity <- function(ds, # data set (long format) parNm, # parameter name parVal = "Value", # column holding the parameter values doseVal = "Dose", # column holding the dose values level = 0.90) # confidence level for the intercept / slope CIs { # Keep only finite rows (the linear model tolerates a zero dose, but drop NA/Inf). ok <- is.finite(ds[[doseVal]]) & is.finite(ds[[parVal]]) if (any(!ok)) warning(sprintf("doseLinearity: dropping %d row(s) with non-finite values.", sum(!ok))) ds <- ds[ok, , drop = FALSE] if (length(unique(ds[[doseVal]])) < 2) stop("doseLinearity: need at least two distinct dose levels.") # Clean two-column frame so the fitted model carries tidy variable names. fit_df <- data.frame(param = ds[[parVal]], dose = ds[[doseVal]]) # Linear dose-linearity model: C = alpha0 + alpha * dose. modlin <- lm(param ~ dose, data = fit_df) smry <- summary(modlin) a0.est <- unname(coef(modlin)[1]) # intercept (alpha0) a1.est <- unname(coef(modlin)[2]) # slope (alpha) ci <- confint(modlin, level = level) a0.int <- ci[1, ] # CI of the intercept a1.int <- ci[2, ] # CI of the slope # a0.t <- unname(smry$coefficients[1, "t value"]) # t = estimate / SE, i.e. how many SEs alpha0 sits from 0 a0.p <- unname(smry$coefficients[1, "Pr(>|t|)"]) # p-value of H0: alpha0 = 0 (two-sided) r2 <- smry$r.squared # goodness of the linear fit # Outcome of the intercept test (H0: alpha0 = 0), read off the CI / p-value. # a0 not significant -> line effectively through origin (consistent with, not proof of, proportionality); # a0 significant -> genuine dose-independent offset. if (a0.int[1] <= 0 & a0.int[2] >= 0) { test.res <- "a0 not signif. diff. from 0" } else { test.res <- "a0 signif. diff. from 0" } result.tab <- data.frame(parameter = parNm, `Dosing range` = paste0(min(ds[[doseVal]]), "--", max(ds[[doseVal]])), `Intercept (a0)` = round(a0.est, 3), `Intercept CI` = paste0("(", round(a0.int[1], 3), "; ", round(a0.int[2], 3), ")"), # `t-value (a0)` = round(a0.t, 3), `p-value (a0)` = signif(a0.p, 3), `Slope (a)` = round(a1.est, 3), `Slope CI` = paste0("(", round(a1.int[1], 3), "; ", round(a1.int[2], 3), ")"), `R2` = round(r2, 3), Conclusion = test.res, check.names = FALSE) # Return the tidy table; attach data, model and numeric stats as attributes. attr(result.tab, "data") <- fit_df attr(result.tab, "modlin") <- modlin attr(result.tab, "stats") <- data.frame(a0 = a0.est, a0_low = unname(a0.int[1]), a0_high = unname(a0.int[2]), a0_p = a0.p, # a0_t = a0.t, a1 = a1.est, a1_low = unname(a1.int[1]), a1_high = unname(a1.int[2]), r2 = r2) result.tab } The core of the function is `lm(param ~ dose)`: R's `lm()` fits a linear regression by ordinary least squares, i.e. it finds the intercept and slope of the straight line that best matches exposure against dose. `confint()` and `summary()` then read the 90% confidence interval and the p-value of the intercept off that fitted model. ## Input data format The function requires **long-format** data with one row per subject, containing (at least) a dose column and a parameter-value column. Only these two columns are used; extra columns are ignored. | id | Dose | parameter | Value | |-----|------|-----------|-------| | 1 | 150 | Cmax | 12.3 | | 2 | 150 | Cmax | 10.8 | | ... | ... | ... | ... | Point the function at the relevant columns via `parVal`/`doseVal` (defaults `"Value"`/`"Dose"`), and set the confidence level via `level` (default `0.90`). Call it **once per parameter**. ## Example Using the `aPCSK9_SAD.pkx` PKanalix demo project, we assess dose proportionality for `Cmax`, `AUClast` and `AUCINF_obs`. Via `lixoftConnectors` function `loadProject()` load the project, restrict the computed parameters to `Dose`, `Cmax`, `AUClast` and `AUCINF_obs` with `setNCASettings()`, run the analysis with `runNCAEstimation()`, and retrieve the individual NCA parameters with `getNCAIndividualParameters()$parameters`. R library(dplyr) library(ggplot2) path.software <- "C:/Program Files/Lixoft/MonolixSuite2024R1" library(lixoftConnectors) initializeLixoftConnectors(path.software, software = "pkanalix") library(flextable) ############################################################################## # helper to prepare the data prepDataset <- function(ds,input_param){ param <- sym(input_param) dt <- ds %>% dplyr::rename(Value=!!param)%>% mutate(parameter=input_param, Value =as.numeric(Value))%>% select(id,Dose,parameter,Value) return(dt) } prj <- paste0(getDemoPath(),"/2.case_studies/project_aPCSK9_SAD.pkx") loadProject(prj) setNCASettings(computedNCAParameters=c("Dose","Cmax","AUClast", "AUCINF_obs")) runNCAEstimation() df_nca <- getNCAIndividualParameters()$parameters df_nca <- rbind(prepDataset(df_nca, "Cmax"), prepDataset(df_nca, "AUClast"), prepDataset(df_nca, "AUCINF_obs")) The retrieved table is in wide format (one column per parameter), so `prepDataset()` reshapes it into **long format** : one row per subject (`id`) and parameter, each row carrying the corresponding `Dose`. Once the dataset is in this format, `doseLinearity()` can be applied (once per parameter). R ############################################################################## # Dose linearity analysis # Assess each parameter with the linear model; intercept judged against 0. DL_params <- unique(df_nca$parameter) tb.res <- data.frame() ci.res <- data.frame() stat.res <- data.frame() for (i in DL_params){ nca.i <- df_nca %>% filter(parameter == i) res.i <- doseLinearity(ds = nca.i, parNm = i, parVal = "Value", doseVal = "Dose", level = 0.90) tb.res <- rbind(tb.res, res.i) # intercept + slope + 90% CIs + R2 st.i <- attr(res.i, "stats"); st.i$parameter <- i stat.res <- rbind(stat.res, st.i) # fitted line + 90% confidence / prediction bands (for the linear-scale plot) mdl.i <- attr(res.i, "modlin") grid.i <- data.frame(dose = seq(min(nca.i$Dose), max(nca.i$Dose), length.out = 100)) ci.i <- predict(mdl.i, newdata = grid.i, interval = "confidence", level = 0.90) pi.i <- predict(mdl.i, newdata = grid.i, interval = "prediction", level = 0.90) ci.res <- rbind(ci.res, data.frame(parameter = i, Dose = grid.i$dose, fit = ci.i[, "fit"], ci_lwr = ci.i[, "lwr"], ci_upr = ci.i[, "upr"], pi_lwr = pi.i[, "lwr"], pi_upr = pi.i[, "upr"])) } flextable(tb.res) %>% autofit() ![image-20260720-100156.png](https://monolixsuite.slp-software.com/__attachments/a_92356ae0d00dd8671b2dcbefa9cc7a569eefaaf88da472ce00548e5721c9d58d/image-20260720-100156.png?cb=c13693d7a2ff31f039b3b2eecabdf3da) ++**Interpretation:**++ Over the 150--800 mg dose range all three parameters return `a0 not significant`: although the intercept point estimates look different from 0 (e.g. roughly -100 for `AUClast`), their confidence intervals are wide and comfortably contain 0, so the intercept cannot be distinguished from zero. The result is *consistent with* dose proportionality, in agreement with the power-model assessment. ### Exposure vs. dose The linear model fit (solid), its 90% confidence interval (dashed) and 90% prediction interval (dotted) for each parameter: R ggplot(df_nca, aes(Dose, Value)) + geom_point(alpha = 0.6, size = 2.5) + geom_line(data = ci.res, aes(Dose, fit, linetype = "Linear-model fit"), colour = "#CD2626", linewidth = 1.2) + geom_line(data = ci.res, aes(Dose, ci_lwr, linetype = "90% confidence interval"), colour = "#1874CD", linewidth = 1.1) + geom_line(data = ci.res, aes(Dose, ci_upr, linetype = "90% confidence interval"), colour = "#1874CD", linewidth = 1.1) + geom_line(data = ci.res, aes(Dose, pi_lwr, linetype = "90% prediction interval"), colour = "#008B45", linewidth = 1.0) + geom_line(data = ci.res, aes(Dose, pi_upr, linetype = "90% prediction interval"), colour = "#008B45", linewidth = 1.0) + scale_linetype_manual(name = NULL, breaks = c("Linear-model fit", "90% confidence interval", "90% prediction interval"), values = c("Linear-model fit" = "solid", "90% confidence interval" = "dashed", "90% prediction interval" = "dotted")) + expand_limits(x = 0) + guides(linetype = guide_legend()) + facet_wrap(~ parameter, scales = "free_y") + labs(x = "Dose", y = "Parameter value") + theme_bw() ![DoseLin-20260715-124606.svg](https://monolixsuite.slp-software.com/__attachments/a_8a2749887dade9d9901869d6534c8cf3f0f143aec79ab7c15adb9e8249f832fa/DoseLin-20260715-124606.svg?cb=5ded6508840752c27c063ecd4ba4e12a) Over the 150--800 mg range the observations fall close to the fitted line with a good fit, and when the line is extrapolated back to dose 0 it meets the axis approximatively at the zero-exposure line for `Cmax` and only below it for `AUClast` and `AUCINF_obs`. ### Intercept vs. zero (decision plot) The intercept estimate (point estimate) and its 90% confidence interval R intercept_a0 <- ggplot(stat.res, aes(x = a0, y = factor(parameter, levels = DL_params))) + geom_vline(xintercept = 0, linetype = "dotted", linewidth = 0.6) + geom_errorbar(aes(xmin = a0_low, xmax = a0_high), width = 0.15) + geom_point(size = 2.5) + annotate("text", x = 0, y = Inf, label = "Intercept = 0", vjust = 1.4, fontface = "bold", size = 3.2) + labs(x = expression(alpha[0] ~ "(intercept) with 90% CI"), y = NULL) + theme_bw() + theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5, face = "bold"), plot.margin = margin(t = 18, r = 12, b = 6, l = 6)) ![Intercept_a0-20260715-095501.svg](https://monolixsuite.slp-software.com/__attachments/a_5681d0d585efc73b66da9f30319a822596f985faf0b6d055398f7f1551d7a976/Intercept_a0-20260715-095501.svg?cb=22ac0a5b53e0fb60cded6d520f9028ea) Each point estimates confidence interval contains 0 (crosses the "Intercept = 0" line). For all three parameters the intercept cannot be distinguished from zero consistent with dose proportionality. Last updated: July 20, 2026 --- version: "2024R1" language: "en" --- # Dose proportionality assessment Dose proportionality can be described when the systemic exposure (e.g. Cmax, AUC) increases in direct proportion to the administered dose amount, meaning doubling the dose doubles the exposure. Characterising how exposure changes with dose is a routine part of early pharmacokinetic analysis. When exposure is dose proportional, concentrations at an untested dose could be approximated by simple proportional scaling; departures from proportionality (supra- or sub-proportional increases) point to non-linear pharmacokinetics that require closer attention. The **power-model** assessment described here is a common exploratory tool for this purpose \[[Smith *et al.* 2000](https://physiologie.envt.fr/wp-content/uploads/2011/04/smithConfidence_Interval_Criteria_for_assessment_dose_proportionality_smith_Pharm_research.pdf); [Hummel *et al.*2008](https://onlinelibrary.wiley.com/doi/10.1002/pst.326)\]. ## Concept and mathematics The assessment implemented here uses the **power model** , which relates a PK parameter to dose through The model is fit via a linear regression of on ; the slope is the estimate of interest: * **exact dose proportionality** (exposure scales 1:1 with dose), * exposure increases more than proportionally * exposure increases less than proportionally Because a point estimate of exactly is never observed, dose proportionality is assessed with an **equivalence criterion** (Smith *et al.* 2000; Hummel *et al.* 2008). Given the ratio of the highest to the lowest dose and acceptance bounds for the ratio of dose-normalised exposures over that dose range, proportionality is concluded when the **90% confidence interval of** lies entirely within that critical region ## Choice of acceptance bounds The bounds are set on the **ratio of dose-normalised exposure** between the extremes of the dose range, **not** on the slope directly. Two conventions are used in the literature: | **Bounds** | **Origin** | **When appropriate** | |--------------|--------------------------------------------------------------------------------------------------------------------------------------------|----------------------------------------------------------------------------------------------------------------| | (0.80, 1.25) | *Confidence Interval Criteria for Assessment of Dose Proportionality (Smith et al. 2000)* | Strict; reasonable when doses are only \~2-fold apart. Becomes very hard to satisfy over a wide dose range. | | (0.50, 2.00) | *Exploratory assessment of dose proportionality: review of current approaches and proposal for a practical criterion (Hummel et al. 2008)* | Allows a -fold change in dose-normalised exposure; intended for exploratory assessment over a wide dose range. | Both are symmetric on the log scale ( ). The bounds should be **pre-specified** in the analysis plan and justified by what fold-change in exposure is clinically acceptable. ## Interpretation The conclusion follows from the position of the **90% CI of the slope** relative to the **critical region** where and | **Position of the 90% CI of** | **Conclusion** | |---------------------------------------|------------------| | Entirely **inside** | Proportional | | Entirely **outside** (below or above) | Not proportional | | Overlapping the interval | Inconclusive | ## Function `doseProportionality()` The above explaned methodology can be implemented as follows: it requires the individual NCA parameters retrieved by `getNCAIndividualParameters()$parameters` and returns the assessment: R doseProportionality <- function(ds, # data set (longitudinal format) parNm, # parameter name (label in the output table) parVal = "Value", # column holding the parameter values doseVal = "Dose", # column holding the dose values thetaL = 0.5, # lower acceptance bound for the dose-normalised ratio thetaH = 2.0) # upper acceptance bound for the dose-normalised ratio { # Keep only rows usable on the log-log scale (positive, finite dose and value). ok <- is.finite(ds[[doseVal]]) & is.finite(ds[[parVal]]) & ds[[doseVal]] > 0 & ds[[parVal]] > 0 if (any(!ok)) warning(sprintf("doseProportionality: dropping %d row(s) with non-positive/non-finite values.", sum(!ok))) ds <- ds[ok, , drop = FALSE] if (length(unique(ds[[doseVal]])) < 2) stop("doseProportionality: need at least two distinct dose levels.") # Clean two-column frame so the fitted models carry tidy variable names. fit_df <- data.frame(param = ds[[parVal]], dose = ds[[doseVal]]) rr <- max(fit_df$dose) / min(fit_df$dose) # highest / lowest dose parfit <- lm(log(param) ~ log(dose), data = fit_df) # power fit beta.est <- unname(coef(parfit)[2]) # slope estimate beta.int <- confint(parfit, level = 0.9)[2, ] # 90% CI of the slope llim <- 1 + log(thetaL) / log(rr) # critical region, lower limit ulim <- 1 + log(thetaH) / log(rr) # critical region, upper limit # CI inside -> Proportional; straddling -> Inconclusive; outside -> Not Proportional. if (beta.int[1] > llim & beta.int[2] < ulim) { test.res <- "Proportional" } else if ((beta.int[1] <= llim & (beta.int[2] > llim & beta.int[2] < ulim)) | ((beta.int[1] > llim & beta.int[1] < ulim) & beta.int[2] >= ulim) | ( beta.int[1] <= llim & beta.int[2] >= ulim)) { test.res <- "Inconclusive" } else { test.res <- "Not proportional" } result.tab <- data.frame(parameter = parNm, `Dosing range` = paste0(min(ds[[doseVal]]), "--", max(ds[[doseVal]])), `Point estimate` = round(beta.est, 2), `Criteria interval` = paste0("(", round(llim, 3), "; ", round(ulim, 3), ")"), `Confidence interval` = paste0("(", round(beta.int[1], 3), "; ", round(beta.int[2], 3), ")"), Conclusion = test.res, check.names = FALSE) parfitlin <- lm(param ~ dose - 1, data = fit_df) # through-origin fit, for plotting # Return the tidy table; attach data, models and numeric stats as attributes. attr(result.tab, "data") <- fit_df attr(result.tab, "modlog") <- parfit attr(result.tab, "modlin") <- parfitlin attr(result.tab, "stats") <- data.frame(beta = beta.est, ci_low = unname(beta.int[1]), ci_high = unname(beta.int[2]), llim = llim, ulim = ulim) result.tab } The core of the function is `lm(log(param) ~ log(dose))`: R's `lm()` fits a linear regression by ordinary least squares on the **log-transformed** data, so the fitted line is the power model on the log--log scale and its slope *is* the exponent . `confint()` then reads the 90% confidence interval of that slope off the fitted model. (The second fit, `lm(param ~ dose - 1)`, is a plain through-origin line used only for plotting.) ### Input data format The function requires **long-format** data with one row per subject, containing (at least) a dose column and a parameter-value column. Only these two columns are used; extra columns are ignored. | id | Dose | parameter | Value | |-----|------|-----------|-------| | 1 | 150 | Cmax | 12.3 | | 2 | 150 | Cmax | 10.8 | | ... | ... | ... | ... | Point the function `doseProportionality()` at the relevant columns via `parVal`/`doseVal` (defaults `"Value"`/`"Dose"`), and set the acceptance bounds `thetaL`/`thetaH `. Call it **once per parameter**. ## Example Using the `aPCSK9_SAD.pkx` PKanalix demo project, we assess dose proportionality for `Cmax`, `AUClast` and `AUCINF_obs`. Use the raw exposure parameters, not dose-normalised (`_D`). Pass the **original** `Cmax` / `AUC` parameters to `doseProportionality()`, not dose-normalised (`_D`) parameters. Dose-normalisation is already built in: under the dose-normalised exposure ratio is , so the bounds on map exactly onto the critical region for . A `_D` parameter would give a slope of (proportionality at 0, not 1) and normalise the data twice, invalidating the criterion. Via the `lixoftConnectors` we load the project with `loadProject()`, restrict the computed parameters to `Dose`, `Cmax`, `AUClast` and `AUCINF_obs` with `setNCASettings()`, run the analysis with `runNCAEstimation()`, and retrieve the individual NCA parameters with `getNCAIndividualParameters()$parameters`. R library(dplyr) library(ggplot2) path.software <- "C:/Program Files/Lixoft/MonolixSuite2024R1" library(lixoftConnectors) initializeLixoftConnectors(path.software, software = "pkanalix") library(flextable) ############################################################################## # helper to prepare the data prepDataset <- function(ds,input_param){ param <- sym(input_param) dt <- ds %>% dplyr::rename(Value=!!param)%>% mutate(parameter=input_param, Value =as.numeric(Value))%>% select(id,Dose,parameter,Value) return(dt) } prj <- paste0(getDemoPath(),"/2.case_studies/project_aPCSK9_SAD.pkx") loadProject(prj) setNCASettings(computedNCAParameters=c("Dose","Cmax","AUClast", "AUCINF_obs")) runNCAEstimation() df_nca <- getNCAIndividualParameters()$parameters df_nca <- rbind(prepDataset(df_nca, "Cmax"), prepDataset(df_nca, "AUClast"), prepDataset(df_nca, "AUCINF_obs")) The retrieved table is in wide format (one column per parameter), so `prepDataset()` reshapes it into **long format** : one row per subject (`id`) and parameter, each row carrying the corresponding `Dose`. Once the dataset is in this format, `doseProportionality()` can be applied (once per parameter). R ############################################################################## # Dose proportionality analysis # Assess each parameter; theta bounds set to the (0.5, 2.0) criterion DP_params <- unique(df_nca$parameter) tb.res <- data.frame() ci.res <- data.frame() stat.res <- data.frame() for (i in DP_params) { nca.i <- df_nca %>% filter(parameter == i) res.i <- doseProportionality(ds = nca.i, parNm = i, parVal = "Value", doseVal = "Dose", thetaL = 0.5, thetaH = 2.0) tb.res <- rbind(tb.res, res.i) # slope + 90% CI + critical region st.i <- attr(res.i, "stats"); st.i$parameter <- i stat.res <- rbind(stat.res, st.i) # fitted line + 90% confidence / prediction bands (for the log-log plot) mdl.i <- attr(res.i, "modlog") grid.i <- data.frame(dose = 10^seq(log10(min(nca.i$Dose)), log10(max(nca.i$Dose)), length.out = 100)) ci.i <- exp(predict(mdl.i, newdata = grid.i, interval = "confidence", level = 0.90)) pi.i <- exp(predict(mdl.i, newdata = grid.i, interval = "prediction", level = 0.90)) ci.res <- rbind(ci.res, data.frame(parameter = i, Dose = grid.i$dose, fit = ci.i[, "fit"], ci_lwr = ci.i[, "lwr"], ci_upr = ci.i[, "upr"], pi_lwr = pi.i[, "lwr"], pi_upr = pi.i[, "upr"])) } flextable(tb.res) %>% autofit() ![image-20260710-090543.png](https://monolixsuite.slp-software.com/__attachments/a_d6c81dcbca6ca9c23f8555e5167a21b69ae0dd1c71dee622418acaeec7dca6e4/image-20260710-090543.png?cb=468b3b489d1ceff4a0f9004cbc594ed0) ++**Interpretation:**++ Over the 150--800 mg dose range, **Cmax is dose‑proportional** (β = 1.15; 90% CI 0.974--1.318, entirely within the acceptance region 0.586--1.414). For **AUClast and AUCINF_obs the assessment is inconclusive** (β = 1.27; 90% CI 1.067--1.474 and 1.074--1.474), since the upper CI bounds slightly exceeded the upper limit (1.414), which can indicate a possible tendency toward more‑than‑proportional total exposure. ### Exposure vs. dose (log--log) The power-model fit (solid), its 90% confidence interval (dashed) and 90% prediction interval (dotted) for each parameter: R ggplot(df_nca, aes(Dose, Value)) + geom_point(alpha = 0.6) + geom_line(data = ci.res, aes(Dose, fit, linetype = "Power-model fit"), colour = "#CD2626", linewidth = 1.2) + geom_line(data = ci.res, aes(Dose, ci_lwr, linetype = "90% confidence interval"), colour = "#1874CD", linewidth = 1.1) + geom_line(data = ci.res, aes(Dose, ci_upr, linetype = "90% confidence interval"), colour = "#1874CD", linewidth = 1.1) + geom_line(data = ci.res, aes(Dose, pi_lwr, linetype = "90% prediction interval"), colour = "#008B45", linewidth = 1.0) + geom_line(data = ci.res, aes(Dose, pi_upr, linetype = "90% prediction interval"), colour = "#008B45", linewidth = 1.0) + scale_linetype_manual(name = NULL, breaks = c("Power-model fit", "90% confidence interval", "90% prediction interval"), values = c("Power-model fit" = "solid", "90% confidence interval" = "dashed", "90% prediction interval" = "dotted")) + guides(linetype = guide_legend()) + # green legend keys scale_x_log10() + scale_y_log10() + annotation_logticks(sides = "bl") + facet_wrap(~ parameter, scales = "free_y") + labs(x = "log(Dose)", y = "log(Parameter value)") + theme_bw() ![DP_PCSK9-20260709-164046.svg](https://monolixsuite.slp-software.com/__attachments/a_794574ae3129244d27627801dc48eb9eaf86acb1dddc7195629129614f38a802/DP_PCSK9-20260709-164046.svg?cb=d81db7b7dfe7a321f7655dd4f40b7f0f) On the log--log scale, Cmax, AUClast and AUCINF_obs increased approximately linearly with dose (power‑model slopes near 1, all points within the 90% prediction interval). ### Slope vs. critical region (decision plot) The slope estimate (point estimate) and its 90% CI (bar) for each parameter, shown vs the critical-region limits (dotted lines). A parameter is **proportional**when its CI lies entirely between the two limits. The limits are drawn as single vertical lines because all parameters share the same dose range. The critical region is identical for each parameter. R llim_val <- stat.res$llim[1] ulim_val <- stat.res$ulim[1] ggplot(stat.res, aes(x = beta, y = factor(parameter, levels = DP_params))) + geom_vline(xintercept = c(llim_val, ulim_val), linetype = "dotted", linewidth = 0.6) + geom_errorbarh(aes(xmin = ci_low, xmax = ci_high), height = 0.15) + geom_point(size = 2.5) + annotate("text", x = llim_val, y = Inf, label = "Lower limit", vjust = 1.4, fontface = "bold", size = 3.2) + annotate("text", x = ulim_val, y = Inf, label = "Upper limit", vjust = 1.4, fontface = "bold", size = 3.2) + scale_x_continuous(limits = c(0, max(2, max(stat.res$ci_high, ulim_val) * 1.05))) + coord_cartesian(clip = "off") + labs(title = "Dose proportionality evaluation", x = expression(beta ~ "with 90% CI"), y = NULL) + theme_bw() + theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5, face = "bold"), plot.margin = margin(t = 18, r = 12, b = 6, l = 6)) ![SlopeVsCriticalRegion-20260710-120310.svg](https://monolixsuite.slp-software.com/__attachments/a_6ea7fa839b9f94a744780f2a481ed7e8a3671f70d95b76d73a178747de72cc77/SlopeVsCriticalRegion-20260710-120310.svg?cb=6ccf6b9187b7b0aa375aef1460a63c87) The power‑model slope with its 90% CI relative to the acceptance limits (0.586--1.414) shows `Cmax `lying fully within the region (dose‑proportional), whereas `AUClast `and `AUCINF_obs `intervals exceed the upper limit (inconclusive). Last updated: July 15, 2026 [Next Page](https://monolixsuite.slp-software.com/llms-full.txt/1)