--- title: "Design Evaluation and Optimization in Discrete Space" date: "" output: rmarkdown::html_vignette: toc: true toc_depth: 3 bibliography: references.bib biblio-style: apalike link-citations: yes vignette: > %\VignetteIndexEntry{Design Evaluation and Optimization in Discrete Space} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} backup_options = options() options(width = 120) knitr::opts_chunk$set( echo = TRUE, eval = FALSE, # chunks do NOT execute by default (CRAN compliance) warning = FALSE, message = FALSE, comment = "#>", cache = FALSE, fig.width = 8, fig.height = 5, fig.align = "center", out.width = "90%" ) library(PFIM) # paths$data : pre-computed RDS files shipped in inst/extdata/ # paths$outputs : text outputs -> tempdir() (CRAN-compliant) # paths$reports : HTML reports -> tempdir() paths = list( data = system.file("extdata", package = "PFIM"), outputs = file.path(tempdir(), "pfim_out"), reports = file.path(tempdir(), "pfim_rep") ) dir.create(paths$outputs, showWarnings = FALSE, recursive = TRUE) dir.create(paths$reports, showWarnings = FALSE, recursive = TRUE) # Shared plot options plotOptions = list(unitTime = "hour", unitOutcomes = c("mcg/mL", "DI%")) ``` --- # Scientific Background This example is based on a landmark PK/PD population study of **Tolmetin** (a non-steroidal anti-inflammatory drug) in rats [@FloresMurrieta1998]. The integrated model couples a one-compartment PK model with a PD model to describe nociception, through the dysfunction index. The antinociceptive effect of Tolmetin was characterized using an indirect response model. **Original experimental design.** Six parallel dose groups of rats (n=6 per group, total N=36) received single oral doses from 1 to 100 mg/kg. Blood samples and inflammation scores (DI) were collected at 10 time points: 0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 2, 3, and 4 h. **Objectives:** 1. **Evaluation** -- Compute the Population and Bayesian Fisher Information Matrices (FIM) for the original design. Assess parameter precision via RSE% and Shrinkage. 2. **Optimization** -- Find the D-optimal design for 30 rats by selecting doses and sampling times from a discrete candidate set. Two algorithms are compared: Fedorov-Wynn and Multiplicative. --- # Design Evaluation ## Model Equations ### PK model: one-compartment, first-order oral absorption $$\frac{dC_c}{dt} = \frac{D}{V}\,k_a\,e^{-k_a t} - \frac{Cl}{V}\,C_c$$ | Symbol | Description | Unit | |:------:|:------------|:----:| | $C_c$ | Plasma drug concentration (state variable) | mcg/mL | | $V$ | Volume of distribution | L | | $k_a$ | First-order absorption rate constant | h$^{-1}$ | | $Cl$ | Total clearance | L/h | | $D$ | Administered dose (`dose_RespPK`) | mg | ### PD model: indirect response, inhibition of production $$\frac{dE}{dt} = R_{in}\!\left(1 - \frac{I_{max}\,C_c^{\gamma}}{C_c^{\gamma} + IC_{50}^{\gamma}}\right) - k_{out}\,E$$ | Symbol | Description | Unit | |:------:|:------------|:----:| | $E$ | Dysfunction index $%$ | -- | | $R_{in}$ | Zero-order input rate | $%/h$ | | $I_{max}$| Maximum fractional inhibition | -- | | $IC_{50}$| Concentration at 50% of $I_{max}$ | $mg/mL$ | | $\gamma$ | Hill coefficient | -- | | $k_{out}$| First-order dissipation rate | $/h$ | > **Steady-state initial condition.** At $t = 0$ (i.e. $C_c = 0$), the PD system > is at equilibrium: $E(0) = R_{in}/k_{out}. **Naming conventions for ODE models:** - The prefix `Deriv_` is mandatory; the suffix must match the state variable name (`Cc` or `E`). - The dose is specified as `dose_` followed by the name of the PK response (e.g. `dose_RespPK`). - Use `**` for exponentiation. - Model outcomes are represented as a list associating each response with its state variable: `RespPK` $\leftrightarrow$ `Cc`, `RespPD` $\leftrightarrow$ `E`. ```{r model-equations, echo=TRUE, comment=''} modelEquations = list( Deriv_Cc = "dose_RespPK/V*ka*exp(-ka*t) - Cl/V*Cc", Deriv_E = "Rin*(1-Imax*(Cc**gamma)/(Cc**gamma + IC50**gamma)) - kout*E" ) ``` ## Model Parameters Parameters are specified via their population typical value (fixed effect $\mu$) and inter-individual variability (IIV $\omega$). We assume a **log-normal** distribution for all parameters. The estimable parameter vector has dimension $p = 14$ (6 fixed effects + 8 variance components ) $$\theta = \bigl\{\mu_V,\,\mu_{Cl},\,\mu_{k_{out}},\,\mu_{I_{max}},\, \mu_{IC_{50}},\,\mu_{\gamma},\; \omega^2_V,\,\omega^2_{Cl},\,\omega^2_{k_{out}},\,\omega^2_{I_{max}},\, \omega^2_{IC_{50}},\,\omega^2_{\gamma}, \sigma_{slope_{RespPK}}, \sigma_{inter_{RespPD}}\bigr\}$$ Setting `fixedMu = TRUE` excludes the fixed effect from the FIM. Setting `omega = 0` implies no IIV; the corresponding variance is not estimated (`fixedOmega = TRUE` is set automatically). | Parameter | Description | $\mu$ | $\omega$ | Fixed $\mu$ | Fixed $\omega$ | |:----------:|:------------------------------------|--------:|---------:|:-----------:|:--------------:| | $V$ | Volume of distribution (L) | 0.74 | 0.316 | No | No | | $Cl$ | Total clearance (L/h) | 0.28 | 0.456 | No | No | | $k_a$ | Absorption rate constant (h$^{-1}$) | 10 | 0 | Yes | Yes | | $k_{out}$ | Effect elimination rate (h$^{-1}$) | 6.14 | 0.947 | No | No | | $R_{in}$ | Baseline production rate (h$^{-1}$) | 614 | 0 | Yes | Yes | | $I_{max}$ | Maximum inhibition (--) | 0.76 | 0.439 | No | No | | $IC_{50}$ | Potency (mcg/mL) | 9.22 | 0.452 | No | No | | $\gamma$ | Hill coefficient (--) | 2.77 | 1.761 | No | No | ```{r model-parameters, echo=TRUE, comment=''} modelParameters = list( ModelParameter(name = "V", distribution = LogNormal(mu = 0.74, omega = 0.316)), ModelParameter(name = "Cl", distribution = LogNormal(mu = 0.28, omega = 0.456)), ModelParameter(name = "ka", distribution = LogNormal(mu = 10, omega = 0), fixedMu = TRUE), ModelParameter(name = "kout", distribution = LogNormal(mu = 6.14, omega = 0.947)), ModelParameter(name = "Rin", distribution = LogNormal(mu = 614, omega = 0), fixedMu = TRUE), ModelParameter(name = "Imax", distribution = LogNormal(mu = 0.76, omega = 0.439)), ModelParameter(name = "IC50", distribution = LogNormal(mu = 9.22, omega = 0.452)), ModelParameter(name = "gamma", distribution = LogNormal(mu = 2.77, omega = 1.761)) ) ``` ## Residual Error Models Two residual error structures are used, one per response. A proportional error model is considered for the PK (`Combined1` with $\sigma_{inter}$= 0, which is equivalent to `Proportional`). An additive error model is considered for the PD (Constant). | Response | Error model | $\sigma_{inter}$ | $\sigma_{slope}$ | |:--------:|:-----------:|:----------------:|:----------------:| | PK | `Combined1` | 0.0 | 0.21 | | PD | `Constant` | 9.6 | -- | ```{r error-models, echo=TRUE, comment=''} errorModelRespPK = Combined1(output = "RespPK", sigmaInter = 0, sigmaSlope = 0.21) errorModelRespPD = Constant(output = "RespPD", sigmaInter = 9.6) modelError = list(errorModelRespPK, errorModelRespPD) ``` ## Sampling Times ```{r sampling-times, echo = TRUE, comment=''} samplingTimesRespPK = SamplingTimes( outcome = "RespPK", samplings = c(0.25, 0.5, 0.75, 1, 1.25, 1.5, 2, 3, 4) ) samplingTimesRespPD = SamplingTimes( outcome = "RespPD", samplings = c(0.25, 0.5, 0.75, 1, 1.25, 1.5, 2, 3, 4) ) ``` ## Arms Six arms correspond to the original dose levels converted to absolute doses for a 200 g rat: $D_\text{mg} = D_\text{mg/kg} \times 0.200\,\text{kg}$. **Total: 36 subjects** (6 arms x 6 subjects each). | Arm | mg/kg | Absolute dose | Subjects | |:---:|------:|:-------------:|:--------:| | 1 | 1.0 | 0.20 mg | 6 | | 2 | 3.2 | 0.64 mg | 6 | | 3 | 10.0 | 2.00 mg | 6 | | 4 | 31.6 | 6.32 mg | 6 | | 5 | 56.2 | 11.24 mg | 6 | | 6 | 100.0 | 20.00 mg | 6 | Initial conditions: $C_c(0) = 0$ and $E(0) = 100 = R_{in}/k_{out}$ ```{r arms, echo = TRUE, comment=''} # Helper to avoid repeating the Arm() constructor six times makeArm = function(armName, dose, size = 6) { admin = Administration(outcome = "RespPK", timeDose = 0, dose = dose) Arm( name = armName, size = size, administrations = list(admin), samplingTimes = list(samplingTimesRespPK, samplingTimesRespPD), initialCondition = list(Cc = 0, E = 100) ) } arm1 = makeArm("0.20mg Arm", 0.20) arm2 = makeArm("0.64mg Arm", 0.64) arm3 = makeArm("2.00mg Arm", 2.00) arm4 = makeArm("6.32mg Arm", 6.32) arm5 = makeArm("11.24mg Arm", 11.24) arm6 = makeArm("20.00mg Arm", 20.00) ``` ## 1.6 Design Assembly ```{r design1, echo = TRUE, comment=''} design1 = Design( name = "design1", arms = list( arm1, arm2, arm3, arm4, arm5, arm6 ) ) ``` # Evaluation of the population and Bayesian FIM ## Population FIM Evaluation ```{r evaluation-objects, echo = TRUE, comment=''} evaluationPop = Evaluation( name = "evaluationPop", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, outputs = list("RespPK" = "Cc", "RespPD" = "E"), designs = list(design1), fimType = "population", odeSolverParameters = list(atol = 1e-8, rtol = 1e-8) ) evaluationPopResults = run (evaluationPop ) ``` ## Bayesian FIM Evaluation ```{r evaluation-run, echo=TRUE, comment=''} evaluationBayesian = Evaluation( name = "evaluationBayesian", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, outputs = list("RespPK" = "Cc", "RespPD" = "E"), designs = list(design1), fimType = "Bayesian", odeSolverParameters = list(atol = 1e-8, rtol = 1e-8) ) evaluationBayesianResults = run( evaluationBayesian ) ``` ## Results ### Population FIM ```{r results-pop-show, echo = TRUE, eval= FALSE, comment=''} show(evaluationPopResults) ``` ```{r results-pop-show_from_outputs, echo = FALSE, eval=TRUE, comment=''} lines = readLines("outputs/vignette1_evaluation_populationFIM_show.txt") cat(paste(lines, collapse = "\n")) ``` All these elements can also be accessed using the following methods: ```{r results-pop-detail} cat("Fisher Information Matrix") print(getFisherMatrix(evaluationPopResults)) cat("Correlation Matrix") print(getCorrelationMatrix(evaluationPopResults)) cat("Standard Errors (SE)") print(getSE(evaluationPopResults)) cat("Relative Standard Errors") print(getRSE(evaluationPopResults)) cat("Shrinkage (%)") print(getShrinkage(evaluationPopResults)) cat("Determinant") print(getDeterminant(evaluationPopResults)) cat("D-Criterion") print(getDcriterion(evaluationPopResults)) ``` ### Bayesian FIM ```{r results-bayesian-show, echo = TRUE, eval= FALSE, comment=''} show(evaluationBayesianResults) ``` ```{r results-bayesian_from_outputs, echo = FALSE, eval=TRUE, comment=''} lines = readLines("outputs/vignette1_evaluation_BayesianFIM_show.txt") cat(paste(lines, collapse = "\n")) ``` ```{r results-bayesian-detail} cat("Fisher Information Matrix") print(getFisherMatrix(evaluationBayesianResults)) cat("Correlation Matrix") print(getCorrelationMatrix(evaluationBayesianResults)) cat("Standard Errors (SE)") print(getSE(evaluationBayesianResults)) cat("Relative Standard Errorn") print(getRSE(evaluationBayesianResults)) cat("Shrinkage (%)") print(getShrinkage(evaluationBayesianResults)) cat("Determinant") print(getDeterminant(evaluationBayesianResults)) cat("D-Criterion") print(getDcriterion(evaluationBayesianResults)) ``` ## Diagnostic Plots ### Response profiles ```{r plot-eval-pk, echo = TRUE, eval=FALSE, comment=''} # plot() is the unified OO entry point — dispatches on class, returns a named list: # $evaluation -> nested [["design"]][["arm"]][["outcome"]] # $sensitivityIndices -> nested [["design"]][["arm"]][["outcome"]][["param"]] # $SE / $RSE -> ggplot2 bar charts # which = c(...) selects a subset; omitting it computes all plots for the class. evalPlots = plot(evaluationPopResults, plotOptions = plotOptions, which = c("evaluation", "sensitivityIndices", "SE", "RSE")) print(evalPlots$evaluation[["design1"]][["20.00mg Arm"]][["RespPK"]]) ``` ```{r plot-eval-pk_from_figures, echo=FALSE, eval=TRUE, out.width="50%", fig.align="center", fig.cap="PK response profile -- 20.00 mg arm (population FIM)"} knitr::include_graphics("figures/vignette1_evaluation_populationFim_design1_arm20mg_RespPK.png") ``` ```{r plot-eval-pd, echo = TRUE, eval=FALSE, comment=''} print(evalPlots$evaluation[["design1"]][["20.00mg Arm"]][["RespPD"]]) ``` ```{r plot-eval-pd_from_figures, echo = FALSE, eval=TRUE, out.width="50%",comment='', fig.align = "center", fig.cap="PD response profile -- 20.00 mg arm (population FIM)"} knitr::include_graphics("figures/vignette1_evaluation_populationFim_design1_arm20mg_RespPD.png") ``` ### Sensitivity indices ```{r plot-si-cl, echo = TRUE, eval= FALSE, comment=''} # $sensitivityIndices is computed in the same plot() call above print(evalPlots$sensitivityIndices[["design1"]][["20.00mg Arm"]][["RespPK"]][["Cl"]]) ``` ```{r plot-si-cl_from_figures, echo = FALSE, eval=TRUE, out.width="50%", comment='', fig.align = "center", fig.cap="Sensitivity index for Cl -- RespPK, 20.00 mg arm"} knitr::include_graphics("figures/vignette1_evaluation_populationFim_design1_SI_RespPK_Cl.png") ``` ### Standard errors and Relative Standard errors ```{r plot-se, echo = TRUE, eval= FALSE, comment=''} print(evalPlots$SE) ``` ```{r plot-se_figures, out.width="33%", echo = FALSE, eval= TRUE, comment='', fig.align = "center",fig.cap="Standard Errors (SE)"} knitr::include_graphics("figures/vignette1_evaluation_populationFim_design1_SE.png") ``` ```{r plot-rse, fig.width=5, fig.width=15, echo = TRUE, eval= FALSE, comment=''} print(evalPlots$RSE) ``` ```{r plot-rse_figures, out.width="33%", echo = FALSE, eval= TRUE, comment='', fig.align = "center",fig.cap="Relative Standard Errors (RSE %)"} knitr::include_graphics("figures/vignette1_evaluation_populationFim_design1_RSE.png") ``` ## HTML Report ```{r report-pop, eval=FALSE, echo=TRUE, comment=''} #Define your path to save your report: pathsReports ="C:/..." Report(evaluationPopResults, pathsReports, "vignette1_evaluation_popFIM.html", plotOptions) ``` --- # Design Optimization **Objective.** Find a D-optimal design for a future study under practical constraints: - Total sample size fixed at **30 subjects**. - Doses restricted to the discrete set: {0.2, 0.64, 2, 6.32, 11.24, 20} mg. - **RespPK:** 2 mandatory times + 4 optimisable times (from 7 candidates). - **RespPD:** 2 mandatory times + 4 optimisable times (from 8 candidates). ## Shared Optimization Setup ```{r optim-setup, echo = TRUE, eval = FALSE, comment=''} # Initial administration # Starting dose; optimizer will reassign from the discrete set in Section 2.6. adminRespPK = Administration(outcome = "RespPK", timeDose = 0, dose = 6.32) # Candidate sampling grids sampTimesOptPK = SamplingTimes( outcome = "RespPK", samplings = c(0.25, 0.75, 1, 1.5, 2, 4, 6) ) sampTimesOptPD = SamplingTimes( outcome = "RespPD", samplings = c(0.25, 0.75, 1.5, 2, 3, 6, 8, 12) ) # Sampling constraints -- RespPK # Fixed: 0.25 h (absorption/Cmax) and 4 h (late elimination). # Optimisable: 4 times chosen from {0.75, 1, 1.5, 2, 6}. sampConstraintsPK = SamplingTimeConstraints( outcome = "RespPK", initialSamplings = c(0.25, 0.75, 1, 1.5, 2, 4, 6), fixedTimes = c(0.25, 4), numberOfsamplingsOptimisable = 4 ) # Sampling constraints -- RespPD # Fixed: 2 h (near peak effect / IC50) and 6 h (mid-recovery / kout). # Optimisable: 4 times chosen from {0.25, 0.75, 1.5, 3, 8, 12}. sampConstraintsPD = SamplingTimeConstraints( outcome = "RespPD", initialSamplings = c(0.25, 0.75, 1.5, 2, 3, 6, 8, 12), fixedTimes = c(2, 6), numberOfsamplingsOptimisable = 4 ) # 2.5 Initial elementary protocols (Fedorov-Wynn only) initialElementaryProtocols = list( c(0.25, 0.75, 1, 4), # PK-focused: absorption + elimination c(1.5, 2, 6, 12) # PD-focused: near-peak + late recovery ) # Dose constraints -- discrete set matching the original study adminConstraintsPK = AdministrationConstraints( outcome = "RespPK", doses = list(0.2, 0.64, 2, 6.32, 11.24, 20) ) # Constrained arm and design # E(0) = "Rin/kout" as a formula string: PFIM evaluates this at typical values, # giving E0 = 614/6.14 = 100. This is preferable to hardcoding the numeric # value because it stays consistent if parameter estimates are updated. armConstraint = Arm( name = "armConstraint", size = 30, administrations = list(adminRespPK), samplingTimes = list(sampTimesOptPK, sampTimesOptPD), administrationsConstraints = list(adminConstraintsPK), samplingTimesConstraints = list(sampConstraintsPK, sampConstraintsPD), initialCondition = list(Cc = 0, E = "Rin/kout") ) # numberOfArms: upper bound on distinct elementary protocols designConstraint = Design( name = "designConstraints", arms = list(armConstraint), numberOfArms = 30 ) numberOfSubjects = 30 proportionsOfSubjects = 1 # single group at init; optimizer redistributes ``` ## Fedorov-Wynn Algorithm ```{r fw-object, echo = TRUE, eval=FALSE, comment=''} optimizationFW = Optimization( name = "FedorovWynn", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, optimizer = "FedorovWynnAlgorithm", optimizerParameters = list( elementaryProtocols = initialElementaryProtocols, numberOfSubjects = numberOfSubjects, proportionsOfSubjects = proportionsOfSubjects, showProcess = TRUE ), designs = list(designConstraint), fimType = "population", outputs = list("RespPK" = "Cc", "RespPD" = "E"), odeSolverParameters = list(atol = 1e-8, rtol = 1e-8) ) ``` ```{r fw-run, echo=TRUE, eval=FALSE, comment=''} optimizationFWResults = run(optimizationFW) ``` ### Results ```{r fw-show, echo=TRUE, eval=FALSE, comment=''} show(optimizationFWResults) ``` ```{r results-optiFW-popFIM_from_outputs, echo = FALSE, eval=TRUE, comment=''} lines = readLines("outputs/vignette1_optimization_FedorovWynn_populationFIM_show.txt") cat(paste(lines, collapse = "\n")) ``` ```{r fw-detail, echo=TRUE, eval=FALSE, comment=''} cat("Fisher Information Matrix\n"); print(getFisherMatrix(optimizationFWResults)) cat("\nCorrelation Matrix\n"); print(getCorrelationMatrix(optimizationFWResults)) cat("\nStandard Errors (SE)\n"); print(getSE(optimizationFWResults)) cat("\nRelative Standard Errors\n"); print(getRSE(optimizationFWResults)) cat("\nShrinkage (%)\n"); print(getShrinkage(optimizationFWResults)) cat("\nDeterminant\n"); print(getDeterminant(optimizationFWResults)) cat("\nD-Criterion\n"); print(getDcriterion(optimizationFWResults)) ``` ### HTML Report ```{r fw-report, eval=FALSE, echo=TRUE, comment=''} #Define your path to save your report: pathsReports ="C:/..." Report(optimizationFWResults, pathsReports, "vignette1_optimization_FedorovWynn_populationFIM.html", plotOptions) ``` ## Multiplicative Algorithm Unlike Fedorov-Wynn algorithm, the Multiplicative Algorithm requires no initial protocols and may converge to a different local optimum -- comparing both confirms robustness. ```{r mult-object, echo=TRUE, eval=FALSE, comment=''} optimizationMult = Optimization( name = "Multiplicative", modelEquations = modelEquations, modelParameters = modelParameters, modelError = modelError, optimizer = "MultiplicativeAlgorithm", optimizerParameters = list( lambda = 0.99, numberOfIterations = 1000, weightThreshold = 0.01, delta = 1e-4, showProcess = TRUE ), designs = list( designConstraint), fimType = "population", outputs = list( "RespPK" = "Cc", "RespPD" = "E" ), odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) ) ``` ```{r mult-run, echo=TRUE, eval=FALSE, comment=''} optimizationMultResults = run( optimizationMult ) ``` ### Results ```{r mult-show, echo=TRUE, eval=FALSE, comment=''} show(optimizationMultResults) ``` ```{r results-optiAlgoMult-popFIM_from_outputs, echo = FALSE, eval=TRUE, comment=''} lines = readLines("outputs/vignette1_optimization_MultiplicativeAlgorithm_populationFIM_show.txt") cat(paste(lines, collapse = "\n")) ``` ```{r mult-detail} cat("Fisher Information Matrix") print(getFisherMatrix(optimizationMultResults)) cat("Correlation Matrix") print(getCorrelationMatrix(optimizationMultResults)) cat("Standard Errors (SE)") print(getSE(optimizationMultResults)) cat("Relative Standard Errors") print(getRSE(optimizationMultResults)) cat("Shrinkage (%)") print(getShrinkage(optimizationMultResults)) cat("Determinant") print(getDeterminant(optimizationMultResults)) cat("D-Criterion") print(getDcriterion(optimizationMultResults)) ``` ### Arm weights ```{r algoMul-show, echo=TRUE, eval=FALSE, comment=''} # plot() on an Optimization object -- $weights is specific to MultiplicativeAlgorithm: # final weight per candidate protocol; non-zero bars define the design support. plotsMult = plot(optimizationMultResults, plotOptions = plotOptions, which = c("evaluation", "SE", "RSE", "weights")) print( plotsMult$weights ) ``` ```{r mult-plot, out.width="50%", echo = FALSE, eval= TRUE, comment='', fig.align = "center",fig.cap="Weights"} knitr::include_graphics("figures/vignette1_optimization_MultiplicativeAlgorithm_populationFIM_weights.png") ``` ### HTML Report ```{r mult-report, eval=FALSE, echo=TRUE, comment=''} #Define your path to save your report: pathsReports ="C:/... Report(optimizationMultResults,pathsReports ,"vignette1_optimization_Multiplicative_populationFIM.html",plotOptions) ``` --- # References ```{r teardown, echo=FALSE, include=FALSE} options(backup_options) ```