---
title: "Design evaluation and optimization in discrete space"
classoption: openany
output:
  rmarkdown::html_vignette:
    toc: true
bibliography: references.bib
biblio-style: apalike
link-citations: yes
linkcolor: blue
urlcolor: green
vignette: >
  %\VignetteIndexEntry{Design evaluation and optimization in discrete space}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<style type="text/css">

body {
  font-size: 11pt;
  max-width: 1400px !important;
  width: 92%;
  margin: 0 auto !important;
  padding: 0 1.5rem;
  line-height: 1.45;
}

h1.title { font-size: 28pt; border-bottom: none !important; }
h1, h2, h3, h4, h5, h6 {
  border-bottom: none !important;
  box-shadow: none !important;
}
h1 { font-size: 18pt; margin-top: 1.4em; }
h2 { font-size: 14pt; margin-top: 1.2em; }
h3 { font-size: 12pt; margin-top: 1.0em; }

/* PFIM R script boxes: single black frame (div only; pre has no extra border) */
div.sourceCode {
  border: 1.5px solid #000 !important;
  background: #fafafa !important;
  padding: 0.7em 0.9em !important;
  overflow-x: auto;
  margin: 0.8em 0;
}
div.sourceCode pre,
div.sourceCode pre.sourceCode,
pre.sourceCode {
  border: none !important;
  background: transparent !important;
  padding: 0 !important;
  margin: 0 !important;
  overflow-x: visible;
}
/* Fallback when code is a bare <pre> (no div.sourceCode wrapper) */
pre:not(.sourceCode) {
  border: 1.5px solid #000 !important;
  background: #fafafa !important;
  padding: 0.7em 0.9em !important;
  overflow-x: auto;
  margin: 0.8em 0;
}
code.r,
pre code {
  font-size: 10.5pt;
}

img {
  max-width: 100%;
  height: auto;
  display: block;
  margin: 0.6em auto;
}

table {
  margin: 0.9em auto 1.1em auto;
  border-collapse: collapse;
  font-size: 10.5pt;
}

table caption {
  caption-side: bottom;
  font-style: italic;
  font-size: 10pt;
  padding-top: 0.4em;
  color: #333;
}

th, td {
  padding: 5px 11px;
  vertical-align: middle;
}

table.table > thead > tr > th,
table > thead > tr > th,
thead th {
  background-color: #f3f4f6 !important;
  border-bottom: 1px solid #bbb;
}

table.table > tbody > tr > td,
tbody tr,
tbody tr:nth-child(even),
tbody tr:nth-child(odd) {
  background-color: #ffffff !important;
}


.pfim-show {
  font-size: 11px;
  line-height: 1.25;
  overflow-x: auto;
  background: #fafafa;
  border: 1.5px solid #000;
  padding: 0.6em 0.8em;
}

</style>

```{r global_options, echo = FALSE, include = FALSE}
knitr::opts_knit$set(tangle = FALSE)
backup_options = options()
library(PFIM)
set.seed(42)
options(width = 200)
utils = system.file("vignette-scripts", "pfim-vignette-utils.R", package = "PFIM")
if (!nzchar(utils)) stop("pfim-vignette-utils.R not found.", call. = FALSE)
source(utils, local = knitr::knit_global())
paths = pfimVignetteSetupPaths()
plotOptions = list(unitTime = c("hour"), unitOutcomes = c("mcg/mL", "DI%"))
.pfimVignetteHas = function( name ) {
  exists( name, inherits = TRUE ) && {
    val = get( name, inherits = TRUE )
    !is.null( val ) && ( !is.character( val ) || any( nzchar( val ) ) )
  }
}
knitr::opts_chunk$set(purl = FALSE, collapse = TRUE,
                      comment = "#>", echo = FALSE, warning = FALSE, message = FALSE,
                      cache = FALSE, tidy = FALSE,
                      fig.align = "center", out.width = "100%", dpi = 160,
                      fig.width = 7, fig.height = 4, dev = "png",
                      dev.args = if (isTRUE(capabilities("cairo")))
                        list(png = list(type = "cairo", antialias = "default")) else list())
```

# Overview

This example is based on a landmark PK/PD population study of Tolmetin (a non-steroidal anti-inflammatory drug) in rats [@FloresMurrieta1998].

The model consists of:

- **PK:** one-compartment model with first-order oral absorption.
- **PD:** an indirect response model (Type I) where the drug inhibits the synthesis rate (*Rin*) of an inflammation score (DI) via an Imax function.

## Original experimental design

The original study involved 6 parallel groups of rats (*n* ≥ 6 per group) receiving single oral doses of 1, 3.2, 10, 31.6, 56.2, or 100 mg/kg per os. Blood sampling and drug response (DI score) evaluation were conducted at 0, 15, 30, and 45 min and at 1, 1.25, 1.5, 2, 3 and 4 hours after administration (nine non-zero sampling times in hours: 0.25–4 h).

## Objectives

1. **Evaluation:** compute the Population Fisher Information Matrix (FIM) for the original design to assess the precision of structural parameters and variance components (RSE%, Shrinkage). A Bayesian FIM is also computed for comparison.
2. **Optimization:** identify the D-optimal design for a total of 30 rats by selecting dose levels and sampling times from a discrete candidate set. The Fedorov-Wynn and Multiplicative algorithms are compared.

Optimisation results are computed by `example01_execute.R` (run once, then cached as `.RDS` in `data/`). HTML reports are written to `results/`. Reports are also available at <https://github.com/packagePFIM>

# Design evaluation

## Model equations

The PKPD model is defined as a system of Ordinary Differential Equations (ODEs) using named character strings. PFIM performs symbolic differentiation on these strings to derive the sensitivity equations required for FIM computation.

**Convention:**

- Prefix `Deriv_`: mandatory; identifies the string as the right-hand side of an ODE.
- Suffix: must match the state variable name (`Cc` or `E`).
- Operators: standard R arithmetic; use `**` for exponentiation.

**Equation 1 (PK)** — one-compartment model with first-order oral absorption:

$$\frac{dC_c}{dt} = \frac{\mathrm{dose_{RespPK}}}{V} \cdot k_a \cdot e^{-k_a t} - \frac{Cl}{V} \cdot C_c$$

| Symbol | Description |
|--------|-------------|
| *V* | volume of distribution (L) |
| *ka* | first-order absorption rate constant (h⁻¹) |
| *Cl* | total clearance (L/h) |
| *Cc* | plasma drug concentration — state variable (mcg/mL) |

**Equation 2 (PD)** — indirect response model, inhibition of production (Type I):

$$\frac{dE}{dt} = R_{in} \left(1 - I_{max} \frac{C_c^\gamma}{C_c^\gamma + IC_{50}^\gamma}\right) - k_{out} \cdot E$$

| Symbol | Description |
|--------|-------------|
| *Rin* | baseline production rate of the inflammation score (h⁻¹) |
| *Imax* | maximum fractional inhibition (dimensionless, 0 < Imax ≤ 1) |
| *IC50* | concentration producing 50% of Imax (mcg/mL) |
| *gamma* | Hill coefficient (dimensionless) |
| *kout* | first-order elimination rate of the effect (h⁻¹) |
| *E* | pharmacodynamic response (DI inflammation score) — state variable |

**Steady-state note:** at *t* = 0 (*Cc* = 0) the system is at equilibrium: *E*(0) = *Rin* / *kout* (here 614 / 6.14 = 100).

```{r, echo = TRUE, eval = FALSE, 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*). PFIM assumes a log-normal distribution for all parameters, guaranteeing strict positivity.

The parameter vector estimated by the population FIM is:

$$\theta = \{\mu_V, \mu_{Cl}, \mu_{kout}, \mu_{Imax}, \mu_{IC50}, \mu_{gamma}, \omega^2_V, \omega^2_{Cl}, \omega^2_{kout}, \omega^2_{Imax}, \omega^2_{IC50}, \omega^2_{gamma}\}$$

Parameters with `fixedMu = TRUE` are considered known constants and are excluded from the FIM. Parameters with *omega* = 0 have no IIV component and their variance is not estimated.

| Name | Description | *mu* | *omega* | fixedMu |
|------|-------------|------|---------|---------|
| V | Volume of distribution (L) | 0.74 | 0.316 | FALSE |
| Cl | Total clearance (L/h) | 0.28 | 0.456 | FALSE |
| ka | Absorption rate constant (h⁻¹) | 10 | 0 | TRUE |
| kout | Effect elimination rate (h⁻¹) | 6.14 | 0.947 | FALSE |
| Rin | Baseline production rate (h⁻¹) | 614 | 0 | TRUE |
| Imax | Maximum inhibition (-) | 0.76 | 0.439 | FALSE |
| IC50 | Potency (mcg/mL) | 9.22 | 0.452 | FALSE |
| gamma | Hill coefficient (-) | 2.77 | 1.761 | FALSE |

**Rationale for fixed parameters:**

- *ka*: very rapid oral absorption makes this parameter weakly identifiable from sparse PK data; fixing avoids numerical instability in the FIM.
- *Rin*: entirely determined by the steady-state constraint *E*(0) = *Rin* / *kout*; estimating it jointly with *kout* would create a structural redundancy.

```{r, echo = TRUE, eval = FALSE, 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 distinct residual error models are specified, one per response.

`Combined1(output, sigmaInter, sigmaSlope)` — combined additive + proportional model:

$$\mathrm{SD}(\epsilon) = \sigma_{inter} + \sigma_{slope} \cdot f(\theta, \xi)$$

Setting `sigmaInter = 0` reduces it to a pure proportional model: SD(ε_PK) = 0.21 × *Cc* (21% proportional error). This is appropriate for plasma concentrations where measurement error scales with the signal magnitude across the dynamic range.

`Constant(output, sigmaInter)` — additive error model: SD(ε_PD) = 9.6 DI units. This is appropriate for bounded inflammation scores whose measurement precision does not depend on the response level.

```{r, echo = TRUE, eval = FALSE, comment=''}

errorModelRespPK = Combined1(output = "RespPK", sigmaInter = 0,   sigmaSlope = 0.21)
errorModelRespPD = Constant( output = "RespPD", sigmaInter = 9.6)

modelError = list(errorModelRespPK, errorModelRespPD)
```

## Sampling times

Nine observation times (hours) are used for both responses, spanning:

- Absorption and peak phase: 0.25–1.00 h
- Post-peak / distribution: 1.25–2.00 h
- Elimination / offset phase: 3.00–4.00 h

Time 0 is omitted: *Cc*(0) = 0 and *E*(0) = 100 are fixed initial conditions that carry no information about model parameters (zero sensitivity).

```{r, echo = TRUE, eval = FALSE, 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 (dose groups)

Six arms correspond to the original dose levels, converted from mg/kg to absolute doses for a 200 g rat (dose_mg = dose_mg/kg × 0.200 kg).

| 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 |

Design summary: 6 arms × 6 subjects = 36 subjects total.

Initial conditions:

- *Cc*(0) = 0: no drug present at baseline.
- *E*(0) = 100: inflammation score at steady-state (= *Rin* / *kout* = 614 / 6.14).

```{r, echo = TRUE, eval = FALSE, comment=''}

administrationRespPK1 = Administration(
  outcome  = "RespPK",
  timeDose = c(0),
  dose     = c(0.2)
)

arm1 = Arm(
  name             = "0.2mg Arm",
  size             = 6,
  administrations  = list(administrationRespPK1),
  samplingTimes    = list(samplingTimesRespPK, samplingTimesRespPD),
  initialCondition = list("Cc" = 0, "E" = 100)
)

administrationRespPK2 = Administration(outcome = "RespPK", timeDose = c(0), dose = c(0.64))

arm2 = Arm(
  name             = "0.64mg Arm",
  size             = 6,
  administrations  = list(administrationRespPK2),
  samplingTimes    = list(samplingTimesRespPK, samplingTimesRespPD),
  initialCondition = list("Cc" = 0, "E" = 100)
)

administrationRespPK3 = Administration(outcome = "RespPK", timeDose = c(0), dose = c(2))

arm3 = Arm(
  name             = "2mg Arm",
  size             = 6,
  administrations  = list(administrationRespPK3),
  samplingTimes    = list(samplingTimesRespPK, samplingTimesRespPD),
  initialCondition = list("Cc" = 0, "E" = 100)
)

administrationRespPK4 = Administration(outcome = "RespPK", timeDose = c(0), dose = c(6.32))

arm4 = Arm(
  name             = "6.32mg Arm",
  size             = 6,
  administrations  = list(administrationRespPK4),
  samplingTimes    = list(samplingTimesRespPK, samplingTimesRespPD),
  initialCondition = list("Cc" = 0, "E" = 100)
)

administrationRespPK5 = Administration(outcome = "RespPK", timeDose = c(0), dose = c(11.24))

arm5 = Arm(
  name             = "11.24mg Arm",
  size             = 6,
  administrations  = list(administrationRespPK5),
  samplingTimes    = list(samplingTimesRespPK, samplingTimesRespPD),
  initialCondition = list("Cc" = 0, "E" = 100)
)

administrationRespPK6 = Administration(outcome = "RespPK", timeDose = c(0), dose = c(20))

arm6 = Arm(
  name             = "20mg Arm",
  size             = 6,
  administrations  = list(administrationRespPK6),
  samplingTimes    = list(samplingTimesRespPK, samplingTimesRespPD),
  initialCondition = list("Cc" = 0, "E" = 100)
)
```

## Design assembly

`Design()` aggregates all arms into a single experimental design object.

```{r, echo = TRUE, eval = FALSE, comment=''}

design1 = Design(
  name = "design1",
  arms = list(arm1, arm2, arm3, arm4, arm5, arm6)
)
```

## Population and Bayesian FIM evaluation

The `Evaluation()` constructor specifies the full statistical model:

| Argument | Role |
|:--------:|:-----|
| `modelEquations` | user-defined ODE system |
| `modelParameters` | fixed effects + IIV |
| `modelError` | intra-individual error |
| `outputs` | named list mapping outcome labels to ODE state variables (`"RespPK"` → `Cc`, `"RespPD"` → `E`) |
| `designs` | list of `Design` objects to evaluate |
| `fimType` | `"population"` estimates θ = {*mu*, ω²}; `"Bayesian"` is the individual FIM regularized by a prior on ω² |
| `odeSolverParameters` | passed to `deSolve::lsoda`; tight tolerances (1e-8) are required because the PD sub-model is moderately stiff (*kout* = 6.14 h⁻¹ implies rapid equilibration of the effect) |

PFIM applies a First-Order (FO) linearization, expanding the individual model in a first-order Taylor series around the typical values *mu*. The population FIM has dimension *p* × *p*, where *p* = number of estimable parameters (here *p* = 12: 6 fixed effects + 6 variance components).

### Population FIM

```{r, echo = TRUE, eval = FALSE, comment=''}

evaluationPop = Evaluation(
  name                = "evaluation",
  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)
)

evaluationPop = run(evaluationPop)
```

### Bayesian FIM

The Bayesian FIM augments the individual FIM with the inverse prior covariance matrix (Ω⁻¹), acting as a regularization term. This is equivalent to a Maximum A Posteriori (MAP) estimation framework and is relevant when prior information on ω² is available. All other arguments are identical to the population evaluation above.

```{r, echo = TRUE, eval = FALSE, comment=''}

evaluationBay = Evaluation(
  name                = "evaluation",
  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)
)

evaluationBay = run(evaluationBay)
```

```{r ex01_run, include = FALSE}
script = system.file("vignette-scripts", "example01_execute.R", package = "PFIM")
if (!nzchar(script) || !file.exists(script)) {
  stop("example01_execute.R not found; reinstall PFIM or rebuild vignettes.")
}
source(script, local = knitr::knit_global())
```

## Results: display and export

Accessor functions for retrieving results from `Evaluation` / `Optimization` objects:

| Function | Description |
|:--------:|:------------|
| `show(x)` | formatted summary of all statistical metrics |
| `getFisherMatrix(x)` | retrieves the FIM |
| `getCorrelationMatrix(x)` | normalized FIM to identify parameter collinearity |
| `getSE(x)` | asymptotic Standard Errors (SE) |
| `getRSE(x)` | Relative Standard Errors (RSE%) |
| `getShrinkage(x)` | Shrinkage (%) for random effects |
| `getDeterminant(x)` | determinant of the FIM |
| `getDcriterion(x)` | D-optimality criterion of the FIM |

```{r, echo = TRUE, eval = FALSE, comment=''}

show(evaluationPop)
writeLines(capture.output(show(evaluationPop)),
           file.path(paths$outputs, "vignette1_evaluation_populationFIM_show.txt"))
fisherMatrix = getFisherMatrix(evaluationPop)
getCorrelationMatrix(evaluationPop)
getSE(evaluationPop)
getRSE(evaluationPop)
getShrinkage(evaluationPop)
getDeterminant(evaluationPop)
getDcriterion(evaluationPop)

show(evaluationBay)
writeLines(capture.output(show(evaluationBay)),
           file.path(paths$outputs, "vignette1_evaluation_BayesianFIM_show.txt"))
fisherMatrix = getFisherMatrix(evaluationBay)
getCorrelationMatrix(evaluationBay)
getSE(evaluationBay)
getRSE(evaluationBay)
getShrinkage(evaluationBay)
getDeterminant(evaluationBay)
getDcriterion(evaluationBay)
```

```{r ex01_show_pop, echo = FALSE, results = "asis"}
cat("<pre class=\"pfim-show\">", showOutputEvaluationPop, "</pre>\n", sep = "")
```

```{r ex01_show_bay, echo = FALSE, results = "asis"}
cat("<pre class=\"pfim-show\">", showOutputEvaluationBay, "</pre>\n", sep = "")
```

## Diagnostic plots

`plotEvaluation()`, `plotSensitivityIndices()`, `plotSE()`, and `plotRSE()` are the PFIM entry points for evaluation graphics:

- `plotEvaluation()` — model predictions vs sampling times, nested as `result$designName$armName$outcomeName`
- `plotSensitivityIndices()` — sensitivity index curves, nested as `result$design$arm$outcome$param`
- `plotSE()` / `plotRSE()` — ggplot2 bar charts of Standard and Relative Standard Errors

`plotOptions` controls axis labels in all PFIM graphics.

```{r, echo = TRUE, eval = FALSE, comment=''}

plotOptions = list(unitTime = c("hour"), unitOutcomes = c("mcg/mL", "DI%"))

plotsEval1_eval = plotEvaluation(evaluationPop, plotOptions)
plotsEval1_si   = plotSensitivityIndices(evaluationPop, plotOptions)

plotOutcomesEvaluationRespPK = plotsEval1_eval$design1$`20mg Arm`$RespPK
plotOutcomesEvaluationRespPD = plotsEval1_eval$design1$`20mg Arm`$RespPD
plotSensitivityIndice_RespPK_Cl = plotsEval1_si$design1$`20mg Arm`$RespPK$Cl
plotSensitivityIndice_RespPK_V = plotsEval1_si$design1$`20mg Arm`$RespPK$V
plotEval_SE = PFIM::plotSE(evaluationPop)
plotEval_RSE = PFIM::plotRSE(evaluationPop)

ggsave(file.path(paths$figures, "vignette1_evaluation_populationFim_design1_arm20mg_RespPK.pdf"),
       plotOutcomesEvaluationRespPK, width = 8, height = 5)
```

```{r ex01_plot_eval_resppk, echo = FALSE, fig.width = 10, fig.height = 5.5, out.width = "100%"}
plotOutcomesEvaluationRespPK
```

```{r ex01_plot_eval_resppd, echo = FALSE, fig.width = 10, fig.height = 5.5, out.width = "100%"}
plotOutcomesEvaluationRespPD
```

```{r ex01_plot_si_cl, echo = FALSE, fig.width = 10, fig.height = 5.5, out.width = "100%"}
plotSensitivityIndice_RespPK_Cl
```

```{r ex01_plot_si_v, echo = FALSE, fig.width = 10, fig.height = 5.5, out.width = "100%"}
plotSensitivityIndice_RespPK_V
```

```{r ex01_plot_se, echo = FALSE, fig.width = 11, fig.height = 4.4, out.width = "100%"}
plotEval_SE
```

```{r ex01_plot_rse, echo = FALSE, fig.width = 11, fig.height = 4.4, out.width = "100%"}
plotEval_RSE
```

## HTML report

```{r, echo = TRUE, eval = FALSE, comment=''}

outputFile = "Example01_EvaluationPopFIM.html"
Report(evaluationPop, paths$reports, outputFile, plotOptions)
```

# Design optimization

## Objectives and constraints

Building on the evaluation above, we now seek an optimal design for a future study under practical constraints:

- Total sample size fixed at 30 subjects (a single arm, subdivided into elementary protocols by the optimizer)
- Doses restricted to the discrete set: {0.2, 0.64, 2, 6.32, 11.24, 20} mg
- For RespPK: 4 sampling times in total, including 2 fixed times; the remaining 2 are chosen from {0.75, 1, 1.5, 2, 6}
- For RespPD: 4 sampling times in total, including 2 fixed times; the remaining 2 are chosen from {0.25, 0.75, 1.5, 3, 8, 12}

## Algorithms compared

Both algorithms operate in the discrete candidate space and maximize the D-criterion of the population FIM.

**Fedorov-Wynn (FW):** an exact exchange algorithm that iteratively adds the elementary protocol (dose × sampling schedule pair) that most increases the FIM determinant, then removes the least contributing one. It converges to a D-optimal design on the discrete support of candidate protocols.

**Multiplicative Algorithm (MA):** a continuous relaxation method that assigns and iteratively updates weights to all candidate protocols. Weights below a threshold are zeroed at convergence, yielding a sparse approximate D-optimal design. Unlike FW, MA does not require initial elementary protocols and explores the full candidate space simultaneously.

## Runtime note

Both algorithms require approximately 10–15 minutes to run on a standard workstation. During vignette rendering, `example01_execute.R` runs each optimization once (with `showProcess = FALSE`) and saves the result to `data/`; subsequent renders load the cached `.RDS` files.

## Initial administration and candidate sampling grids

Starting dose for the constrained arm. The optimizer will reassign doses to subjects from the discrete set defined in the administration constraints below.

```{r, echo = TRUE, eval = FALSE, comment=''}

administrationRespPK = Administration(outcome = "RespPK", timeDose = c(0), dose = c(6.32))
```

These are the full sets of candidate time points (hours) from which the optimizer will select the most informative subset, subject to the constraints defined in the next sections.

```{r, echo = TRUE, eval = FALSE, comment=''}

samplingTimesRespPK = SamplingTimes(
  outcome   = "RespPK",
  samplings = c(0.25, 0.75, 1, 1.5, 2, 4, 6)
)

samplingTimesRespPD = SamplingTimes(
  outcome   = "RespPD",
  samplings = c(0.25, 0.75, 1.5, 2, 3, 6, 8, 12)
)
```

## Sampling time constraints

`SamplingTimeConstraints(outcome, initialSamplings, fixedTimes, numberOfsamplingsOptimisable, ...)` restricts which time points can be assigned; one object per outcome in the arm.

**RespPK constraints:**

- `fixedTimes`: 0.25 h captures the rising absorption phase and Cmax region; 4.0 h anchors the late elimination phase for reliable Cl/V estimation.
- `numberOfsamplingsOptimisable = 4`: total times per protocol for this outcome (the 2 fixed times plus 2 free times chosen from {0.75, 1, 1.5, 2, 6}).

**RespPD constraints:**

- `fixedTimes`: 2 h is near the expected peak effect (captures Emax region and IC50 estimation); 6 h is the mid-recovery phase (informative for *kout* estimation).
- `numberOfsamplingsOptimisable = 4`: total times per protocol (2 fixed plus 2 free from {0.25, 0.75, 1.5, 3, 8, 12}).

```{r, echo = TRUE, eval = FALSE, comment=''}

samplingConstraintsRespPK = SamplingTimeConstraints(
  outcome                      = "RespPK",
  initialSamplings             = c(0.25, 0.75, 1, 1.5, 2, 4, 6),
  fixedTimes                   = c(0.25, 4),
  numberOfsamplingsOptimisable = 4
)

samplingConstraintsRespPD = SamplingTimeConstraints(
  outcome                      = "RespPD",
  initialSamplings             = c(0.25, 0.75, 1.5, 2, 3, 6, 8, 12),
  fixedTimes                   = c(2, 6),
  numberOfsamplingsOptimisable = 4
)
```

## Initial elementary protocols for Fedorov-Wynn

An elementary protocol is a (dose, sampling schedule) combination for a homogeneous subgroup of subjects. In a PK/PD arm the sampling schedule is **multi-outcome**: PK times and PD times are concatenated on the FW candidate grid.

The initial support used here is **one** protocol (with `proportionsOfSubjects` of length 1):

- RespPK times: {0.25, 0.75, 1, 4} h  
- RespPD times: {1.5, 2, 6, 12} h  

Pass it as a nested list (one element = one support point), or equivalently as a single flat vector of length 8. A vignette-style `list(pk, pd)` with one proportion is also accepted.

`initialElementaryProtocols` is passed to `optimizerParameters$elementaryProtocols` inside the `Optimization()` call; it is only used by `FedorovWynnAlgorithm`.

```{r, echo = TRUE, eval = FALSE, comment=''}

initialElementaryProtocols = list(
  list(
    c(0.25, 0.75, 1, 4),
    c(1.5, 2, 6, 12)
  )
)
```

## Dose constraints

`AdministrationConstraints(outcome, doses)` restricts which dose values the optimizer can assign to each elementary protocol. The discrete set (in mg) corresponds to the 6 dose levels of the original study.

```{r, echo = TRUE, eval = FALSE, comment=''}

administrationConstraintsRespPK = AdministrationConstraints(
  outcome = "RespPK",
  doses   = list(0.2, 0.64, 2, 6.32, 11.24, 20)
)
```

## Constrained arm and design

The arm encodes all constraints simultaneously and serves as the template that both optimization algorithms will operate on.

- `administrationsConstraints`: list of `AdministrationConstraints` objects; restricts which doses can be assigned.
- `samplingTimesConstraints`: list of `SamplingTimeConstraints` objects; one per outcome.

The initial condition for *E* is specified as `"Rin/kout"` (a formula string) rather than the numeric value 100. PFIM evaluates this expression at the typical parameter values, giving *E*â‚€ = 614 / 6.14 = 100. This ensures that the initial condition remains consistent with the model structure if parameter estimates are updated.

`numberOfArms` is the upper bound on the number of distinct elementary protocols the optimizer can create.

```{r, echo = TRUE, eval = FALSE, comment=''}

armConstraint = Arm(
  name                       = "armConstraint",
  size                       = 30,
  administrations            = list(administrationRespPK),
  samplingTimes              = list(samplingTimesRespPK, samplingTimesRespPD),
  administrationsConstraints = list(administrationConstraintsRespPK),
  samplingTimesConstraints   = list(samplingConstraintsRespPK, samplingConstraintsRespPD),
  initialCondition           = list("Cc" = 0, "E" = "Rin/kout")
)

designConstraint = Design(
  name         = "designConstraint",
  arms         = list(armConstraint),
  numberOfArms = 30
)

numberOfSubjects      = c(30)
proportionsOfSubjects = c(30) / 30
```

For very large dose × sampling grids, cap FIM evaluations (deterministic subsample):

```{r, echo = TRUE, eval = FALSE, comment=''}

pfim_set_option(constraints.maxTasks = 500)
```

## Fedorov-Wynn algorithm

The Fedorov-Wynn algorithm is an exact combinatorial exchange method for finding D-optimal designs in a discrete candidate space. It proceeds as:

1. Start from the user-supplied initial elementary protocols.
2. Evaluate the directional derivative of the D-criterion for every candidate protocol not currently in the design support.
3. Add the protocol with the highest derivative (greedy ascent step).
4. Remove the protocol that contributes least to the current FIM.
5. Repeat steps 2–4 until the improvement in |FIM| falls below a tolerance.

`FedorovWynnAlgorithm` `optimizerParameters`:

| Parameter | Description |
|:---------:|:------------|
| `elementaryProtocols` | list of protocols; each protocol is a flat numeric vector (full grid row) or a list of per-outcome vectors (e.g. PK then PD) |
| `numberOfSubjects` | integer vector; total *N* to distribute across protocols |
| `proportionsOfSubjects` | numeric vector summing to 1; initial allocation fractions |
| `showProcess` | logical; if `TRUE`, prints per-iteration D-criterion values |

```{r, echo = TRUE, eval = FALSE, comment=''}

optimizationFWPopFIM = Optimization(
  name                = "PKPD_ODE_multi_doses_populationFIM",
  modelEquations      = modelEquations,
  modelParameters     = modelParameters,
  modelError          = modelError,
  optimizer           = "FedorovWynnAlgorithm",
  optimizerParameters = list(
    elementaryProtocols   = initialElementaryProtocols,
    numberOfSubjects      = numberOfSubjects,
    proportionsOfSubjects = proportionsOfSubjects,
    showProcess           = FALSE
  ),
  designs             = list(designConstraint),
  fimType             = "population",
  outputs             = list("RespPK" = "Cc", "RespPD" = "E"),
  odeSolverParameters = list(atol = 1e-8, rtol = 1e-8)
)

optimizationFWPopFIM = run(optimizationFWPopFIM)
saveRDS(optimizationFWPopFIM,
        file.path(paths$data, "vignette_1_optimization_FedorovWynn_populationFIM.RDS"))
```

### Display and plot Fedorov-Wynn results

Comparing RSE and D-criterion between `evaluationPop` and `optimizationFWPopFIM` quantifies the information gain from optimization with 30 vs 36 subjects.

For `FedorovWynnAlgorithm`, `plotFrequencies()` returns a bar chart showing how the 30 subjects are distributed across the selected elementary protocols. The number of non-zero bars is the support size of the D-optimal design.

```{r, echo = TRUE, eval = FALSE, comment=''}

show(optimizationFWPopFIM)
writeLines(capture.output(show(optimizationFWPopFIM)),
           file.path(paths$outputs, "vignette1_optimization_FedorovWynn_populationFIM_show.txt"))
fisherMatrix = getFisherMatrix(optimizationFWPopFIM)
getCorrelationMatrix(optimizationFWPopFIM)
getSE(optimizationFWPopFIM)
getRSE(optimizationFWPopFIM)
getShrinkage(optimizationFWPopFIM)
getDeterminant(optimizationFWPopFIM)
getDcriterion(optimizationFWPopFIM)

plotFWFrequencies = PFIM::plotFrequencies(optimizationFWPopFIM)
plotFW_SE = PFIM::plotSE(optimizationFWPopFIM)
plotFW_RSE = PFIM::plotRSE(optimizationFWPopFIM)
plotFWFrequencies
```

```{r ex01_show_fw, echo = FALSE, results = "asis", eval = .pfimVignetteHas("showOutputFW")}
cat("<pre class=\"pfim-show\">", showOutputFW, "</pre>\n", sep = "")
```

```{r ex01_plot_fw_freq, echo = FALSE, fig.width = 6.5, fig.height = 3.8, out.width = "60%", eval = .pfimVignetteHas("plotFWFrequencies")}
plotFWFrequencies
```

```{r ex01_plot_fw_se, echo = FALSE, fig.width = 11, fig.height = 4.4, out.width = "100%", eval = .pfimVignetteHas("plotFW_SE")}
plotFW_SE
```

```{r ex01_plot_fw_rse, echo = FALSE, fig.width = 11, fig.height = 4.4, out.width = "100%", eval = .pfimVignetteHas("plotFW_RSE")}
plotFW_RSE
```

```{r, echo = TRUE, eval = FALSE, comment=''}

outputFile = "Example01_OptimizationFWPopFIM.html"
Report(optimizationFWPopFIM, paths$reports, outputFile, plotOptions)
```

## Multiplicative algorithm

The Multiplicative algorithm (cocktail / multiplicative weights update) is an iterative continuous relaxation approach to D-optimal design in a discrete candidate space. It proceeds as:

1. Initialize uniform weights *w*â‚– = 1/*K* for each of *K* candidate protocols.
2. At each iteration, update weights multiplicatively: *w*ₖ(*t*+1) ∝ *w*ₖ(*t*) × *d*ₖ(ξₜ), where *d*ₖ(ξₜ) is the normalized directional derivative of the D-criterion for protocol *k* under the current design ξₜ.
3. After convergence (change in D-criterion < *delta*), set to zero all weights below `weightThreshold`, yielding a sparse approximate D-optimal design.

Unlike the Fedorov-Wynn algorithm, the MA explores the full candidate space simultaneously without requiring initial elementary protocols, and may converge to a different local optimum.

`MultiplicativeAlgorithm` `optimizerParameters`:

| Parameter | Description |
|:---------:|:------------|
| `lambda` | step-size dampening factor (0 < λ < 1); values close to 1 slow convergence but reduce oscillations |
| `numberOfIterations` | maximum number of multiplicative update cycles |
| `weightThreshold` | protocols with weight < threshold at convergence are zeroed (here 0.01 = < 1% of subjects) |
| `delta` | convergence tolerance on the relative D-criterion change (here 1e-4: stop when improvement < 0.01%) |
| `showProcess` | logical; if `TRUE`, prints per-iteration D-criterion values |

```{r, echo = TRUE, eval = FALSE, comment=''}

optimizationMultPopFIM = Optimization(
  name                = "PKPD_ODE_multi_doses_populationFIM",
  modelEquations      = modelEquations,
  modelParameters     = modelParameters,
  modelError          = modelError,
  optimizer           = "MultiplicativeAlgorithm",
  optimizerParameters = list(
    lambda             = 0.99,
    numberOfIterations = 1000,
    weightThreshold    = 0.01,
    delta              = 1e-04,
    showProcess        = FALSE
  ),
  designs             = list(designConstraint),
  fimType             = "population",
  outputs             = list("RespPK" = "Cc", "RespPD" = "E"),
  odeSolverParameters = list(atol = 1e-8, rtol = 1e-8)
)

optimizationMultPopFIM = run(optimizationMultPopFIM)
saveRDS(optimizationMultPopFIM,
        file.path(paths$data, "vignette_1_optimization_multiplicativeAlgorithm_populationFIM.RDS"))
```

### Display and plot Multiplicative algorithm results

For `MultiplicativeAlgorithm`, `plotWeights()` returns a bar chart of the final weight of each candidate protocol after convergence. Non-zero weights define the design support; comparing this plot to `plotFrequencies()` from FW reveals whether both algorithms converge to the same support, confirming robustness.

```{r, echo = TRUE, eval = FALSE, comment=''}

show(optimizationMultPopFIM)
writeLines(capture.output(show(optimizationMultPopFIM)),
           file.path(paths$outputs, "vignette1_optimization_MultiplicativeAlgorithm_populationFIM_show.txt"))
fisherMatrix = getFisherMatrix(optimizationMultPopFIM)
getCorrelationMatrix(optimizationMultPopFIM)
getSE(optimizationMultPopFIM)
getRSE(optimizationMultPopFIM)
getShrinkage(optimizationMultPopFIM)
getDeterminant(optimizationMultPopFIM)
getDcriterion(optimizationMultPopFIM)

plotMultWeights = PFIM::plotWeights(optimizationMultPopFIM)
plotMult_SE = PFIM::plotSE(optimizationMultPopFIM)
plotMult_RSE = PFIM::plotRSE(optimizationMultPopFIM)
ggsave(file.path(paths$figures, "vignette1_optimization_MultiplicativeAlgorithm_populationFIM_weights.pdf"),
       plotMultWeights, width = 8, height = 5)
plotMultWeights
```

```{r ex01_show_mult, echo = FALSE, results = "asis", eval = .pfimVignetteHas("showOutputMult")}
cat("<pre class=\"pfim-show\">", showOutputMult, "</pre>\n", sep = "")
```

```{r ex01_plot_mult_weights, echo = FALSE, fig.width = 6.5, fig.height = 3.8, out.width = "60%", eval = .pfimVignetteHas("plotMultWeights")}
plotMultWeights
```

```{r ex01_plot_mult_se, echo = FALSE, fig.width = 11, fig.height = 4.4, out.width = "100%", eval = .pfimVignetteHas("plotMult_SE")}
plotMult_SE
```

```{r ex01_plot_mult_rse, echo = FALSE, fig.height = 4.2, fig.width = 9.5, out.width = "100%", eval = .pfimVignetteHas("plotMult_RSE")}
plotMult_RSE
```

```{r, echo = TRUE, eval = FALSE, comment=''}

outputFile = "Example01_OptimizationMultPopFIM.html"
Report(optimizationMultPopFIM, paths$reports, outputFile, plotOptions)
```

# References

```{r global_options_end, echo = FALSE, include = FALSE}
options(backup_options)
```