Design evaluation and optimization in continuous space

Overview

This example is inspired by Sukeishi et al. (2021) (Sukeishi et al. 2022) and illustrates the use of PFIM for a one-compartment PK model with intravenous (IV) infusion and multiple dosing.

Unlike Example 01, where model equations were specified as user-defined ODEs, the PK model here is loaded directly from PFIM’s built-in model library ("Linear1InfusionSingleDose_ClV"). This approach is more concise and avoids the need to manually write sensitivity equations for standard structural models.

Considered original design. A single arm with 150 subjects received multiple IV doses: a loading dose of 400 mg on Day 1, followed by four daily maintenance doses of 200 mg. Each dose was administered as a 1-hour IV infusion. Blood samples were collected at 1, 12, 24, 44, 72, and 120 hours post-first-dose.

Objectives.

  1. Evaluation — Compute and compare three FIM types for the original design: population, individual, and Bayesian FIMs. Assess parameter precision (RSE%) and shrinkage.
  2. Optimization — Find a D-optimal design with only 4 sampling times, using sampling time windows and continuous-space optimization. The number of subjects and the dosing regimen remain unchanged. Three algorithms are compared: PSO (Particle Swarm Optimization), PGBO (Population genetic-based optimization), and Simplex (Nelder–Mead).

Optimisation results are computed by example02_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

PK model from the library of models

PFIM provides a library of pre-implemented standard PK/PD structural models. Using a library model is preferable to a user-defined ODE when the structural form is standard, because it avoids transcription errors and uses pre-computed closed-form sensitivity equations, which are faster and more numerically stable than ODE-based sensitivity computation.

The list modelFromLibrary is passed to Evaluation() or Optimization() in place of modelEquations (the two arguments are mutually exclusive). The key "PKModel" is used for pharmacokinetic models.

The string "Linear1InfusionSingleDose_ClV" decodes as:

Token Meaning
Linear First-order (linear) elimination
1 One-compartment disposition
Infusion IV infusion route; requires Tinf in Administration()
ClV Parameterised by clearance \(Cl\) and volume \(V\)

The predicted concentration during and after infusion follows superposition across repeated doses:

\[C_c(t) = \begin{cases} \dfrac{D}{T_{inf} \cdot Cl}\!\left(1 - e^{-\frac{Cl}{V}t}\right) & 0 \leq t \leq T_{inf} \\[6pt] C_c(T_{inf})\,e^{-\frac{Cl}{V}(t - T_{inf})} & t > T_{inf} \end{cases}\]

Define the PK model Linear1InfusionSingleDose_ClV from the library of models


modelFromLibrary = list("PKModel" = "Linear1InfusionSingleDose_ClV")

Model parameters

The one-compartment model has two structural parameters to estimate, both log-normally distributed:

Parameter Description \(\mu\) \(\omega\)
\(V\) Volume of distribution (L) 50 \(\sqrt{0.26} \approx 0.51\)
\(Cl\) Elimination clearance (L/h) 5 \(\sqrt{0.34} \approx 0.58\)

omega is expressed as the standard deviation on the log scale. Unlike Example 01, no parameters are fixed here (all have inter-individual variability).

Define mu and omega for each parameter

 
modelParameters = list(
  ModelParameter( name = "V",    distribution = LogNormal( mu = 50, omega = sqrt( .26 ) ) ),
  ModelParameter( name = "Cl",   distribution = LogNormal( mu = 5, omega = sqrt( .34 ) ) ) )

Residual error model

A Combined1 (proportional + additive) error model is used for the PK response:

\[\mathrm{SD}(\varepsilon) = \sigma_{\mathrm{inter}} + \sigma_{\mathrm{slope}} \cdot f(\theta, \xi)\]

The additive term dominates at trough concentrations (between doses), while the proportional term prevails near \(C_{\max}\).

Define the error model to the response PK RespPK


errorModelRespPK = Combined1( output = "RespPK", sigmaInter = 0.5, sigmaSlope = sqrt( 0.15 ) )
modelError = list( errorModelRespPK )

Administration parameters

The dosing schedule encodes the full multiple-dose regimen. The attribute Tinf is a numeric vector of infusion durations (hours), one per dose event; it is required for library models of the "Infusion" type.

The loading dose (400 mg = 2× maintenance) rapidly raises \(C_c\) toward the therapeutic target; subsequent 200 mg doses maintain near-steady-state concentrations.

Define the administration parameters of the response PK


administrationRespPK = Administration( outcome = "RespPK",
                                       Tinf = rep( 1, 5 ),
                                       timeDose = seq( 0, 96, 24 ),
                                       dose = c( 400, rep( 200, 4 ) ) )

Sampling times

Six sampling times cover distinct phases of the multi-dose concentration profile:

Time (h) Pharmacokinetic phase
1 End of first infusion — \(C_{\max}\) of dose 1
12 Mid-interdose Day 1 — early elimination
24 Trough before dose 2 — \(C_{\min}\) after dose 1
44 Near \(C_{\max}\) of dose 2 (dose at 24 h + ~20 h)
72 Trough before dose 4 — approaching steady state
120 24 h after the last dose — post-steady-state elimination

Define the sampling times for the response PK RespPK

samplingTimesRespPK = SamplingTimes( outcome = "RespPK", 
                                     samplings = c( 1, 12, 24, 44, 72, 120 ) )

Arm and design

A single arm with all 150 subjects on the same regimen (parallel single-group design). No initialCondition is required: the library model sets \(C_c(0) = 0\) internally for IV infusion models.

Define an arm called arm1 of size 150 with administration administrationRespPK and samplings samplingTimesRespPK


arm1 = Arm( name = "arm1",
            size = 150,
            administrations = list( administrationRespPK ) ,
            samplingTimes   = list( samplingTimesRespPK ) )

Add the arm arm1 to the design design1


design1 = Design( name = "design1", arms = list( arm1 ) )

FIM types: population, individual, and Bayesian

Three Evaluation objects are created, each identical except for fimType. They differ in which parameter subvector enters the FIM and how inter-individual variability (IIV) is handled:

Population FIM (fimType = "population"). The parameter vector includes fixed effects (\(\mu_V\), \(\mu_{Cl}\)) and variance components (\(\omega^2_V\), \(\omega^2_{Cl}\)). This is the standard choice for NLMEM study design.

Individual FIM (fimType = "individual"). Only fixed effects (\(\mu_V\), \(\mu_{Cl}\)) are estimated; \(\omega^2\) are fixed at their specified values. This is relevant when IIV is already well-characterized from a prior study. It tends to give larger (more favorable) FIM values because it ignores the burden of estimating IIV.

Bayesian FIM (fimType = "Bayesian"). The FIM is augmented by a prior precision matrix derived from \(\omega^2\), reflecting the regularization imposed by MAP estimation. This is appropriate for therapeutic drug monitoring or sparse individual PK designs. Shrinkage values are reported for the Bayesian FIM.

Comparing SE and RSE across the three fimType values for the same design quantifies the effect of IIV estimation on parameter precision.

Evaluate the population, individual and Bayesian FIMs


evaluationPop = Evaluation( name = "",
                            modelFromLibrary = modelFromLibrary,
                            modelParameters = modelParameters,
                            modelError = modelError,
                            outputs = list( "RespPK" ),
                            designs = list( design1 ),
                            fimType =  "population",
                            odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) )

evaluationFIMPop = run( evaluationPop )

evaluationInd = Evaluation( name = "",
                            modelFromLibrary = modelFromLibrary,
                            modelParameters = modelParameters,
                            modelError = modelError,
                            outputs = list( "RespPK" ),
                            designs = list( design1 ),
                            fimType = "individual",
                            odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) )

evaluationFIMInd = run( evaluationInd )

evaluationBay = Evaluation( name = "",
                            modelFromLibrary = modelFromLibrary,
                            modelParameters = modelParameters,
                            modelError = modelError,
                            outputs = list( "RespPK" ),
                            designs = list( design1 ),
                            fimType = "Bayesian",
                            odeSolverParameters = list( atol = 1e-8, rtol = 1e-8 ) )

evaluationFIMBay = run( evaluationBay )

Display the results of the design evaluations


show( evaluationFIMPop )
show( evaluationFIMInd )
show( evaluationFIMBay )

fisherMatrix = getFisherMatrix( evaluationFIMPop )
getCorrelationMatrix( evaluationFIMPop )
getSE( evaluationFIMPop )
getRSE( evaluationFIMPop )
getShrinkage( evaluationFIMPop )
getDeterminant( evaluationFIMPop )
getDcriterion( evaluationFIMPop )

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

plotsEval2_eval = plotEvaluation( evaluationFIMPop, plotOptions )
plotsEval2_si   = plotSensitivityIndices( evaluationFIMPop, plotOptions )

plotOutcomesEvaluationRespPK = plotsEval2_eval$design1$arm1$RespPK
plotSensitivityIndice_RespPK_V = plotsEval2_si$design1$arm1$RespPK$V
plotSensitivityIndice_RespPK_Cl = plotsEval2_si$design1$arm1$RespPK$Cl
plotEval_SE = PFIM::plotSE( evaluationFIMPop )
plotEval_RSE = PFIM::plotRSE( evaluationFIMPop )
*************************************** 
  Population Fisher Matrix 
*************************************** 

                      μ_V       μ_Cl      ω²_V     ω²_Cl σ_inter_RespPK σ_slope_RespPK
μ_V             0.1393880 -0.2899109   0.00000   0.00000        0.00000         0.0000
μ_Cl           -0.2899109 14.3204999   0.00000   0.00000        0.00000         0.0000
ω²_V            0.0000000  0.0000000 404.77110  17.51007       51.06227       259.6216
ω²_Cl           0.0000000  0.0000000  17.51007 427.24316       54.71243        80.1080
σ_inter_RespPK  0.0000000  0.0000000  51.06227  54.71243     1982.21804      1361.3658
σ_slope_RespPK  0.0000000  0.0000000 259.62165  80.10800     1361.36582      1648.5030

*************************************** 
  Fixed effects (μ) 
*************************************** 

            μ_V       μ_Cl
μ_V   0.1393880 -0.2899109
μ_Cl -0.2899109 14.3204999

*************************************** 
  Variance components (ω², γ², σ) 
*************************************** 

                    ω²_V     ω²_Cl σ_inter_RespPK σ_slope_RespPK
ω²_V           404.77110  17.51007       51.06227       259.6216
ω²_Cl           17.51007 427.24316       54.71243        80.1080
σ_inter_RespPK  51.06227  54.71243     1982.21804      1361.3658
σ_slope_RespPK 259.62165  80.10800     1361.36582      1648.5030

********************************************* 
  Determinant, condition numbers and D-criterion 
 *********************************************** 

Determinant: 380849659428 
D-criterion: 85.13842 
Condition number (fixed effects): 107.3432 
Condition number (variance components): 12.64203 

*************************************** 
  Parameters estimation 
*************************************** 

Parameter               Value           SE     RSE(%)
μ_V               50.0000000   2.73670901   5.473418
μ_Cl               5.0000000   0.26999904   5.399981
ω²_V              0.2600000   0.05481857  21.084067
ω²_Cl             0.3400000   0.04861157  14.297522
σ_inter_RespPK      0.5000000   0.03570392   7.140784
σ_slope_RespPK      0.3872983   0.04131312  10.667003

[1] 85.13842
*************************************** 
 Individual Fisher Matrix 
*************************************** 

                        μ_V        μ_Cl σ_inter_RespPK σ_slope_RespPK
μ_V             0.003831386 -0.03578704        0.00000        0.00000
μ_Cl           -0.035787035  0.74495930        0.00000        0.00000
σ_inter_RespPK  0.000000000  0.00000000       16.79391       12.15818
σ_slope_RespPK  0.000000000  0.00000000       12.15818       20.61786

*************************************** 
 Fixed effects (μ) 
*************************************** 

              μ_V        μ_Cl
μ_V   0.003831386 -0.03578704
μ_Cl -0.035787035  0.74495930

*************************************** 
 Variance components (σ) 
*************************************** 

               σ_inter_RespPK σ_slope_RespPK
σ_inter_RespPK       16.79391       12.15818
σ_slope_RespPK       12.15818       20.61786

*********************************************** 
  Determinant, condition numbers and D-criterion 
 *********************************************** 

Determinant: 0.3122376 
D-criterion: 0.7475174 
Condition number (fixed effects): 354.3252 
Condition number (variance effects): 4.847154 

*************************************** 
 Parameters estimation 
*************************************** 

Parameter               Value           SE     RSE(%)
μ_V               50.0000000   21.7585945   43.51719
μ_Cl               5.0000000    1.5604237   31.20847
σ_inter_RespPK      0.5000000    0.3223403   64.46806
σ_slope_RespPK      0.3872983    0.2909168   75.11439

*************************************** 
 Bayesian Fisher Matrix 
*************************************** 

           μ_V      μ_Cl
μ_V  13.424620 -8.946759
μ_Cl -8.946759 21.565159

*************************************** 
 Fixed effects 
*************************************** 

           μ_V      μ_Cl
μ_V  13.424620 -8.946759
μ_Cl -8.946759 21.565159

*********************************************** 
 Determinant, condition numbers and D-criterion 
*********************************************** 

Determinant: 209.459560208699 
D-criterion: 14.4727177892993 
Condition number of the fixed effects: 3.56441811487241 

*************************************** 
 Shrinkage 
*************************************** 

               μ_V    μ_Cl
Shrinkage 39.59854 18.8505

*************************************** 
 Parameters estimation 
*************************************** 

Parameter             Value           SE     RSE(%)
μ_V                     50    16.043394   32.08679
μ_Cl                     5     1.265817   25.31634

Create and save the report for the design evaluation


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

outputFile = "Example02_EvaluationPopFIM.html"
Report( evaluationFIMPop, paths$reports, outputFile, plotOptions )

outputFile = "Example02_EvaluationIndFIM.html"
Report( evaluationFIMInd, paths$reports, outputFile, plotOptions )

outputFile = "Example02_EvaluationBayFIM.html"
Report( evaluationFIMBay, paths$reports, outputFile, plotOptions )

Design optimization

Once the pharmacokinetic model has been evaluated using the specified design, the next step is to optimize the sampling times to obtain precise parameter estimates. The number of sampling times per individual is reduced from six in the initial design to four, which helps minimize the burden on patients and biomedical staff, as well as the study cost.

Continuous design space

In contrast to Example 01 (discrete space: finite set of candidate times and doses), the optimization here operates in a continuous space: sampling times are real-valued variables constrained to lie within specified time windows. The algorithms navigate a continuous, non-convex objective function surface rather than enumerate a finite set of candidates.

Unlike the multiplicative and Fedorov–Wynn algorithms (discrete optimization), PSO, PGBO, and Simplex preserve the distribution of individuals across the initial arms and focus solely on optimizing observation times within continuous intervals. In continuous optimization with PFIM, dosing regimens are not optimized; the administration defined for evaluation remains unchanged.

The sampling time constraints involve selecting 4 sampling times according to the following criteria:

Create sampling times for the initial design

We create sampling times that will be used in the initial design for comparison during the optimization process. As continuous optimizers do not optimize dose regimens, we keep the same administration.

samplingTimesRespPK = SamplingTimes( outcome = "RespPK", samplings = c( 1, 48, 72, 120  ) )

Define design constraints

The attributes to be defined are: outcome, initialSamplings, numberOfTimesByWindows (number of samples per window), samplingsWindows (lower and upper bounds of each window), and minSampling (minimum gap between consecutive samples). The sum of numberOfTimesByWindows must equal the total number of sampling times.


samplingConstraintsRespPK  = SamplingTimeConstraints( outcome = "RespPK",
                                                      initialSamplings = c( 1, 48, 72, 120 ),
                                                      numberOfTimesByWindows = c( 2, 2 ),
                                                      samplingsWindows = list( c( 1, 48 ), c( 72, 120 ) ),
                                                      minSampling = 5 )

Create the constraint arm and the associated design


arm2 = Arm( name = "arm2",
            size = 150,
            administrations = list( administrationRespPK ) ,
            samplingTimes   = list( samplingTimesRespPK ),
            samplingTimesConstraints = list( samplingConstraintsRespPK ) )

design2 = Design( name = "design2", arms = list( arm2 ), numberOfArms = 150 )

PSO algorithm (Particle Swarm Optimization)

The Particle Swarm Optimization (PSO) algorithm is a metaheuristic optimization technique inspired by social behaviors observed in birds, fish, and bees. PSO begins with the initialization of a swarm of particles at random positions, each representing a potential solution (a vector of sampling times). Particles iteratively explore the solution space, guided by both their individual experiences (personal best) and the collective knowledge of the swarm (global best).

Key optimizerParameters for PSOAlgorithm:

Set the parameters of the PSO algorithm


optimizationPSO = Optimization( name = "",
                             
                             modelFromLibrary = modelFromLibrary,
                             modelParameters = modelParameters,
                             modelError = modelError,
                             
                             optimizer = "PSOAlgorithm",
                             
                             optimizerParameters = list(
                               maxIteration = 100,
                               populationSize = 50,
                               personalLearningCoefficient = 2.05,
                               globalLearningCoefficient = 2.05,
                               seed = 42,
                               showProcess = FALSE  ),
                             
                             designs = list( design2 ),
                             
                             fimType = "population",
                             
                             outputs = list( "RespPK") )

Run the PSO algorithm for the optimization with a population FIM


optimizationPSO = run(optimizationPSO)
saveRDS(optimizationPSO,
        file.path(paths$data, "vignette2_optimization_PSO_populationFIM.RDS"))

Display the results of the design optimization


show( optimizationPSO )

fisherMatrix = getFisherMatrix( optimizationPSO )
getCorrelationMatrix( optimizationPSO )
getSE( optimizationPSO )
getRSE( optimizationPSO )
getShrinkage( optimizationPSO )
getDeterminant( optimizationPSO )
getDcriterion( optimizationPSO )

plotPSO_SE = PFIM::plotSE( optimizationPSO )
plotPSO_RSE = PFIM::plotRSE( optimizationPSO )
--- Optimal design ---

  Arms name Number of subjects Outcome                    Dose   Sampling times
1      arm2                150  RespPK 400, 200, 200, 200, 200 (1, 24, 97, 120)

*************************************** 
  Population Fisher Matrix 
*************************************** 

                      μ_V       μ_Cl       ω²_V      ω²_Cl σ_inter_RespPK σ_slope_RespPK
μ_V             0.1553906 -0.1846977   0.000000   0.000000        0.00000        0.00000
μ_Cl           -0.1846977 12.4964508   0.000000   0.000000        0.00000        0.00000
ω²_V            0.0000000  0.0000000 503.046381   7.106922       52.27905      247.66451
ω²_Cl           0.0000000  0.0000000   7.106922 325.336003      101.02165       95.47256
σ_inter_RespPK  0.0000000  0.0000000  52.279048 101.021645      699.53147      623.60788
σ_slope_RespPK  0.0000000  0.0000000 247.664509  95.472563      623.60788     1641.23454

*************************************** 
  Fixed effects (μ) 
*************************************** 

            μ_V       μ_Cl
μ_V   0.1553906 -0.1846977
μ_Cl -0.1846977 12.4964508

*************************************** 
  Variance components (ω², γ², σ) 
*************************************** 

                     ω²_V      ω²_Cl σ_inter_RespPK σ_slope_RespPK
ω²_V           503.046381   7.106922       52.27905      247.66451
ω²_Cl            7.106922 325.336003      101.02165       95.47256
σ_inter_RespPK  52.279048 101.021645      699.53147      623.60788
σ_slope_RespPK 247.664509  95.472563      623.60788     1641.23454

********************************************* 
  Determinant, condition numbers and D-criterion 
 *********************************************** 

Determinant: 207863614655 
D-criterion: 76.96556 
Condition number (fixed effects): 81.89388 
Condition number (variance components): 6.914956 

*************************************** 
  Parameters estimation 
*************************************** 

Parameter               Value           SE     RSE(%)
μ_V               50.0000000   2.55938922   5.118778
μ_Cl               5.0000000   0.28540088   5.708018
ω²_V              0.2600000   0.04652999  17.896150
ω²_Cl             0.3400000   0.05673074  16.685513
σ_inter_RespPK      0.5000000   0.04735254   9.470508
σ_slope_RespPK      0.3872983   0.03155626   8.147792

[1] 76.96556

Create and save the report for the design optimization


outputFile = "Example02_OptimizationPSOPopFIM.html"

Report( optimizationPSO, paths$reports, outputFile, plotOptions )

PGBO algorithm (Population genetic-based optimization)

The Population genetic-based optimization (PGBO) algorithm is an evolutionary algorithm that optimizes parameters using selection, mutations, and periodic purge. A mutant is a feasible solution (a vector of sampling times); its fitness is the inverse of the D-criterion for the associated design. Depending on the comparison between mutant and resident fitness, the mutant either replaces the resident or is eliminated. PGBO maintains a balance between converging to the optimum and escaping local minima.

Key optimizerParameters for PGBOAlgorithm:

Set the parameters of the PGBO algorithm


optimizationPGBO = Optimization( name = "",
                                 modelFromLibrary = modelFromLibrary,
                                 modelParameters = modelParameters,
                                 modelError = modelError,
                                 
                                 optimizer = "PGBOAlgorithm",
                                 
                                 optimizerParameters = list(
                                   N = 30,
                                   muteEffect = 0.65,
                                   maxIteration = 1000,
                                   purgeIteration = 200,
                                   seed = 42,
                                   showProcess = FALSE  ),
                                 
                                 designs = list( design2 ),
                                 
                                 fimType = "population",
                                 
                                 outputs = list( "RespPK") )

Run the PGBO algorithm for the optimization with a population FIM


optimizationPGBO = run(optimizationPGBO)
saveRDS(optimizationPGBO,
        file.path(paths$data, "vignette2_optimization_PGBO_populationFIM.RDS"))

Display the results of the design optimization


show( optimizationPGBO )

fisherMatrix = getFisherMatrix( optimizationPGBO )
getCorrelationMatrix( optimizationPGBO )
getSE( optimizationPGBO )
getRSE( optimizationPGBO )
getShrinkage( optimizationPGBO )
getDeterminant( optimizationPGBO )
getDcriterion( optimizationPGBO )

plotPGBO_SE = PFIM::plotSE( optimizationPGBO )
plotPGBO_RSE = PFIM::plotRSE( optimizationPGBO )
--- Optimal design ---

  Arms name Number of subjects Outcome                    Dose               Sampling times
1      arm2                150  RespPK 400, 200, 200, 200, 200 (41.71, 47.62, 81.97, 97.42)

*************************************** 
  Population Fisher Matrix 
*************************************** 

                      μ_V       μ_Cl      ω²_V     ω²_Cl σ_inter_RespPK σ_slope_RespPK
μ_V             0.1041418 -0.2792929   0.00000   0.00000        0.00000        0.00000
μ_Cl           -0.2792929 13.4185743   0.00000   0.00000        0.00000        0.00000
ω²_V            0.0000000  0.0000000 225.94815  16.25095       90.26602      223.80104
ω²_Cl           0.0000000  0.0000000  16.25095 375.12112       72.09968       92.64606
σ_inter_RespPK  0.0000000  0.0000000  90.26602  72.09968      854.87188      826.06775
σ_slope_RespPK  0.0000000  0.0000000 223.80104  92.64606      826.06775     1235.94838

*************************************** 
  Fixed effects (μ) 
*************************************** 

            μ_V       μ_Cl
μ_V   0.1041418 -0.2792929
μ_Cl -0.2792929 13.4185743

*************************************** 
  Variance components (ω², γ², σ) 
*************************************** 

                    ω²_V     ω²_Cl σ_inter_RespPK σ_slope_RespPK
ω²_V           225.94815  16.25095       90.26602      223.80104
ω²_Cl           16.25095 375.12112       72.09968       92.64606
σ_inter_RespPK  90.26602  72.09968      854.87188      826.06775
σ_slope_RespPK 223.80104  92.64606      826.06775     1235.94838

********************************************* 
  Determinant, condition numbers and D-criterion 
 *********************************************** 

Determinant: 31562433262 
D-criterion: 56.21623 
Condition number (fixed effects): 136.5858 
Condition number (variance components): 15.11498 

*************************************** 
  Parameters estimation 
*************************************** 

Parameter               Value           SE     RSE(%)
μ_V               50.0000000   3.18904067   6.378081
μ_Cl               5.0000000   0.28094375   5.618875
ω²_V              0.2600000   0.07585423  29.174706
ω²_Cl             0.3400000   0.05214123  15.335655
σ_inter_RespPK      0.5000000   0.05939051  11.878102
σ_slope_RespPK      0.3872983   0.05339959  13.787716

[1] 56.21623

Create and save the report for the design optimization


outputFile = "Example02_OptimizationPGBOPopFIM.html"

Report( optimizationPGBO, paths$reports, outputFile, plotOptions )

Simplex algorithm (Nelder–Mead)

The Nelder–Mead simplex algorithm is a deterministic derivative-free local optimizer. It operates on a simplex of \(d+1\) vertices in \(\mathbb{R}^d\) (here \(d = 4\) sampling times). At each step, the worst vertex is reflected, expanded, or contracted through the centroid of the remaining vertices until convergence. In this comparison, Simplex serves as a local reference: convergence of PSO, PGBO, and Simplex to the same D-criterion strongly suggests a global optimum.

Key optimizerParameters for SimplexAlgorithm:

Set the parameters of the Simplex algorithm


optimizationSimplex = Optimization( name = "",
                                   modelFromLibrary = modelFromLibrary,
                                   modelParameters = modelParameters,
                                   modelError = modelError,
                                   
                                   optimizer = "SimplexAlgorithm",
                                   
                                   optimizerParameters = list( pctInitialSimplexBuilding = 10,
                                                               maxIteration = 1000,
                                                               tolerance = 1e-10,
                                                               showProcess = FALSE  ),
                                   
                                   designs = list( design2 ),
                                   
                                   fimType = "population",
                                   
                                   outputs = list( "RespPK") )

Run the Simplex algorithm for the optimization with a population FIM


optimizationSimplex = run(optimizationSimplex)
saveRDS(optimizationSimplex,
        file.path(paths$data, "vignette2_optimization_Simplex_populationFIM.RDS"))

Display the results of the design optimization


show( optimizationSimplex )

fisherMatrix = getFisherMatrix( optimizationSimplex )
getCorrelationMatrix( optimizationSimplex )
getSE( optimizationSimplex )
getRSE( optimizationSimplex )
getShrinkage( optimizationSimplex )
getDeterminant( optimizationSimplex )
getDcriterion( optimizationSimplex )

plotSimplex_SE = PFIM::plotSE( optimizationSimplex )
plotSimplex_RSE = PFIM::plotRSE( optimizationSimplex )
--- Optimal design ---

  Arms name Number of subjects Outcome                    Dose   Sampling times
1      arm2                150  RespPK 400, 200, 200, 200, 200 (1, 48, 72, 108)

*************************************** 
  Population Fisher Matrix 
*************************************** 

                      μ_V       μ_Cl      ω²_V     ω²_Cl σ_inter_RespPK σ_slope_RespPK
μ_V             0.1348140 -0.2668512   0.00000   0.00000        0.00000        0.00000
μ_Cl           -0.2668512 12.7511475   0.00000   0.00000        0.00000        0.00000
ω²_V            0.0000000  0.0000000 378.64174  14.83533       51.14649      269.76601
ω²_Cl           0.0000000  0.0000000  14.83533 338.73284       90.52798       91.56249
σ_inter_RespPK  0.0000000  0.0000000  51.14649  90.52798     1082.56615      725.21291
σ_slope_RespPK  0.0000000  0.0000000 269.76601  91.56249      725.21291      892.90492

*************************************** 
  Fixed effects (μ) 
*************************************** 

            μ_V       μ_Cl
μ_V   0.1348140 -0.2668512
μ_Cl -0.2668512 12.7511475

*************************************** 
  Variance components (ω², γ², σ) 
*************************************** 

                    ω²_V     ω²_Cl σ_inter_RespPK σ_slope_RespPK
ω²_V           378.64174  14.83533       51.14649      269.76601
ω²_Cl           14.83533 338.73284       90.52798       91.56249
σ_inter_RespPK  51.14649  90.52798     1082.56615      725.21291
σ_slope_RespPK 269.76601  91.56249      725.21291      892.90492

********************************************* 
  Determinant, condition numbers and D-criterion 
 *********************************************** 

Determinant: 57264569164 
D-criterion: 62.08419 
Condition number (fixed effects): 98.75797 
Condition number (variance components): 13.88599 

*************************************** 
  Parameters estimation 
*************************************** 

Parameter               Value           SE     RSE(%)
μ_V               50.0000000   2.78175798   5.563516
μ_Cl               5.0000000   0.28603036   5.720607
ω²_V              0.2600000   0.06457386  24.836098
ω²_Cl             0.3400000   0.05516612  16.225328
σ_inter_RespPK      0.5000000   0.05010192  10.020383
σ_slope_RespPK      0.3872983   0.06227156  16.078448

[1] 62.08419

Create and save the report for the design optimization


outputFile = "Example02_OptimizationSimplexPopFIM.html"

Report( optimizationSimplex, paths$reports, outputFile, plotOptions )

References

Sukeishi, Asami, Kotaro Itohara, Atsushi Yonezawa, Yuki Sato, Katsuyuki Matsumura, Yoshiki Katada, Takayuki Nakagawa, et al. 2022. “Population Pharmacokinetic Modeling of GS-441524, the Active Metabolite of Remdesivir, in Japanese COVID-19 Patients with Renal Dysfunction.” CPT: Pharmacometrics & Systems Pharmacology 11 (1): 94–103.