This example is taken from Tim Waterhouse’s webinar A Gentle Introduction to Optimal Design for Pharmacometric Models (optimal-design-2023). It shows how PFIM evaluates and optimizes a sparse pharmacokinetic sampling design in pediatric patients.
A population pharmacokinetic model for orally administered Mockdrozaline was developed in adults. The objective here is to assess and improve a sparse pediatric design so that population parameters can be estimated with adequate precision.
The adult model consists of:
Cl,
V, and ka,Clearance and volume enter the oral absorption formula after scaling
with (WT/70)^beta, where WT is body weight in
kg and the exponents are fixed at the adult estimates.
Typical values V, Cl, and ka
refer to a 70 kg adult. Weight exponents beta_V_WT and
beta_Cl_WT are included in the equations but fixed
(fixedMu = TRUE) because they are not re-estimated in the
pediatric study.
modelParameters = list(
ModelParameter(name = "V", distribution = LogNormal(mu = 100, omega = sqrt(0.1))),
ModelParameter(name = "Cl", distribution = LogNormal(mu = 10, omega = sqrt(0.08))),
ModelParameter(name = "ka", distribution = LogNormal(mu = 0.25, omega = sqrt(0.2))),
ModelParameter(
name = "beta_V_WT", distribution = Normal(mu = 1, omega = 0),
fixedMu = TRUE, fixedOmega = TRUE
),
ModelParameter(
name = "beta_Cl_WT", distribution = Normal(mu = 0.75, omega = 0),
fixedMu = TRUE, fixedOmega = TRUE
)
)Concentrations are modelled with combined proportional and additive
residual error (PopED: feps.add.prop). Writing \(f\) for the prediction
(RespPK), both tools use the same variance formula; they
differ in what the numbers mean:
\[ \mathrm{Var}(y \mid f) = \sigma_{\mathrm{add}}^2 + (\sigma_{\mathrm{prop}}\, f)^2 . \]
| Proportional term | Additive term | |
|---|---|---|
PopED / NONMEM (sigma, SIGMA) |
variance 0.05 | variance 1 |
PFIM Combined2 |
SD sigmaSlope = sqrt(0.05) |
SD sigmaInter = 1 |
Tim’s tables sometimes write sigma_prop = 0.05, which
looks like “5% error”. In optimal-design-2023, that value
is still a variance in sigma, not a proportional SD.
Copying 0.05 into sigmaSlope in PFIM would not
reproduce PopED (it would use variance \(0.05^2\) instead of \(0.05\)). Inter-individual variability
follows the same variance → SD rule via omega = sqrt(d) in
the chunk above. Details and a short checklist are in
vignette("Vignette_05_gcse_pkpd", package = "PFIM"),
section PFIM, PopED, and cross-tool equivalence.
Mockdrozaline is administered as 10 mg once daily for 24 weeks
(n = 12 subjects). Doses are expressed in the same units as
in the Tim Waterhouse material (mg × 1000 in the
administration record).
week_h = 24 * 7
lastDose = week_h * 24
dose_times = seq(0, lastDose, by = 24)
administrationRespPK = Administration(
outcome = "RespPK",
timeDose = dose_times,
dose = rep(10 * 1000, times = length(dose_times))
)The initial sparse design includes five samples:
observation_times = c(
5,
week_h * 8 - 1,
week_h * c(12, 24),
lastDose + 168
)
samplingTimesRespPK = SamplingTimes(
outcome = "RespPK",
samplings = observation_times
)
armInitial = Arm(
name = "Initial arm",
size = 12,
administrations = list(administrationRespPK),
samplingTimes = list(samplingTimesRespPK)
)
designInitial = Design(
name = "designInitial",
arms = list(armInitial)
)The population Fisher information matrix (FIM) is computed for the initial design.
Evaluation and PSO runs are computed by
vignette_04_tim_waterhouse_pediatric_execute.R (run once,
then cached as .RDS in data/).
evaluationPopFimInitial = Evaluation(
name = "FIM design initial",
modelParameters = modelParameters,
modelEquations = modelEquations,
modelError = modelError,
designs = list(designInitial),
outputs = list("RespPK"),
fimType = "population"
)
evaluationPopFimInitial = run(evaluationPopFimInitial)#> $designInitial
#> $designInitial$`Initial arm`
#> $designInitial$`Initial arm`$RespPK
show(evaluationPopFimInitial)
#>
#> ***************************************
#> Population Fisher Matrix
#> ***************************************
#>
#> μ_V μ_Cl μ_ka ω²_V ω²_Cl ω²_ka σ_inter_RespPK σ_slope_RespPK
#> μ_V 0.005191013 -0.02212943 -1.264968 0.00000000 0.000000 0.00000000 0.00000000 0.00000
#> μ_Cl -0.022129430 1.24867057 2.833570 0.00000000 0.000000 0.00000000 0.00000000 0.00000
#> μ_ka -1.264968309 2.83356962 317.133094 0.00000000 0.000000 0.00000000 0.00000000 0.00000
#> ω²_V 0.000000000 0.00000000 0.000000 112.27754970 20.404653 41.67043806 0.03434818 42.42871
#> ω²_Cl 0.000000000 0.00000000 0.000000 20.40465277 649.657584 2.09091583 0.12799100 71.00212
#> ω²_ka 0.000000000 0.00000000 0.000000 41.67043806 2.090916 16.36936838 0.01137853 20.54078
#> σ_inter_RespPK 0.000000000 0.00000000 0.000000 0.03434818 0.127991 0.01137853 24.05314802 6.90698
#> σ_slope_RespPK 0.000000000 0.00000000 0.000000 42.42871155 71.002117 20.54078137 6.90697983 965.24067
#>
#> ***************************************
#> Fixed effects (μ)
#> ***************************************
#>
#> μ_V μ_Cl μ_ka
#> μ_V 0.005191013 -0.02212943 -1.264968
#> μ_Cl -0.022129430 1.24867057 2.833570
#> μ_ka -1.264968309 2.83356962 317.133094
#>
#> ***************************************
#> Variance components (ω², γ², σ)
#> ***************************************
#>
#> ω²_V ω²_Cl ω²_ka σ_inter_RespPK σ_slope_RespPK
#> ω²_V 112.27754970 20.404653 41.67043806 0.03434818 42.42871
#> ω²_Cl 20.40465277 649.657584 2.09091583 0.12799100 71.00212
#> ω²_ka 41.67043806 2.090916 16.36936838 0.01137853 20.54078
#> σ_inter_RespPK 0.03434818 0.127991 0.01137853 24.05314802 6.90698
#> σ_slope_RespPK 42.42871155 71.002117 20.54078137 6.90697983 965.24067
#>
#> *********************************************
#> Determinant, condition numbers and D-criterion
#> ***********************************************
#>
#> Determinant: 26095079
#> D-criterion: 8.454148
#> Condition number (fixed effects): 6403571
#> Condition number (variance components): 1353.512
#>
#> ***************************************
#> Parameters estimation
#> ***************************************
#>
#> Parameter Value SE RSE(%)
#> μ_V 100.0000000 142.08529815 142.08530
#> μ_Cl 10.0000000 1.54877798 15.48778
#> μ_ka 0.2500000 0.55839813 223.35925
#> ω²_V 0.1000000 0.41853695 418.53695
#> ω²_Cl 0.0800000 0.04070788 50.88485
#> ω²_ka 0.2000000 1.09951923 549.75962
#> σ_inter_RespPK 1.0000000 0.20411978 20.41198
#> σ_slope_RespPK 0.2236068 0.03318821 14.84222Designs can be compared with the D-criterion; relative standard errors (RSE) are often easier to interpret.
Sampling times can be optimized within protocol windows.
Particle-swarm settings below use a small iteration count so the
vignette renders in reasonable time; increase maxIteration
for production work.
One window covers absorption (0–6 h after the first dose), three windows capture steady-state troughs (weeks 8, 12, and 24), and one window covers 96–168 h after the last dose.
samplingConstraintsRespPK5 = SamplingTimeConstraints(
outcome = "RespPK",
initialSamplings = observation_times,
samplingsWindows = list(
c(0, 6),
c(week_h * 8 - 1, week_h * 8),
c(week_h * 12 - 1, week_h * 12),
c(week_h * 24 - 1, week_h * 24),
c(lastDose + 96, lastDose + 168)
),
numberOfTimesByWindows = c(1, 1, 1, 1, 1),
minSampling = c(1, 1, 1, 1, 1)
)armOptim5 = Arm(
name = "ArmOptim5",
size = 12,
administrations = list(administrationRespPK),
samplingTimes = list(samplingTimesRespPK),
samplingTimesConstraints = list(samplingConstraintsRespPK5)
)
designOptim5 = Design(
name = "designOptim5",
arms = list(armOptim5)
)
optimizationPopFIM5 = Optimization(
name = "Optimization 5 samples",
modelEquations = modelEquations,
modelParameters = modelParameters,
modelError = modelError,
optimizer = "PSOAlgorithm",
optimizerParameters = list(
maxIteration = 10,
populationSize = 50,
personalLearningCoefficient = 2.05,
globalLearningCoefficient = 2.05,
seed = 1234,
showProcess = FALSE
),
designs = list(designOptim5),
fimType = "population",
outputs = list("RespPK")
)
optimizationPopFIM5 = run(optimizationPopFIM5)show(optimizationPopFIM5)
#>
#> --- Optimal design ---
#>
#> Arms name Number of subjects Outcome
#> 1 ArmOptim5 12 RespPK
#> Dose
#> 1 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000
#> Sampling times
#> 1 (0.33, 1344, 2015, 4032, 4128)
#>
#> ***************************************
#> Population Fisher Matrix
#> ***************************************
#>
#> μ_V μ_Cl μ_ka ω²_V ω²_Cl ω²_ka σ_inter_RespPK σ_slope_RespPK
#> μ_V 0.004149819 -0.02856067 -1.327250 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
#> μ_Cl -0.028560666 1.27045005 -1.208882 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
#> μ_ka -1.327250132 -1.20888211 530.944022 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
#> ω²_V 0.000000000 0.00000000 0.000000 71.7541623 33.9879857 45.8748155 0.2105445 14.081909
#> ω²_Cl 0.000000000 0.00000000 0.000000 33.9879857 672.5180543 0.3805719 0.6717039 52.863270
#> ω²_ka 0.000000000 0.00000000 0.000000 45.8748155 0.3805719 45.8824145 0.3673430 23.379602
#> σ_inter_RespPK 0.000000000 0.00000000 0.000000 0.2105445 0.6717039 0.3673430 24.0563894 7.222854
#> σ_slope_RespPK 0.000000000 0.00000000 0.000000 14.0819091 52.8632704 23.3796024 7.2228537 918.086053
#>
#> ***************************************
#> Fixed effects (μ)
#> ***************************************
#>
#> μ_V μ_Cl μ_ka
#> μ_V 0.004149819 -0.02856067 -1.327250
#> μ_Cl -0.028560666 1.27045005 -1.208882
#> μ_ka -1.327250132 -1.20888211 530.944022
#>
#> ***************************************
#> Variance components (ω², γ², σ)
#> ***************************************
#>
#> ω²_V ω²_Cl ω²_ka σ_inter_RespPK σ_slope_RespPK
#> ω²_V 71.7541623 33.9879857 45.8748155 0.2105445 14.081909
#> ω²_Cl 33.9879857 672.5180543 0.3805719 0.6717039 52.863270
#> ω²_ka 45.8748155 0.3805719 45.8824145 0.3673430 23.379602
#> σ_inter_RespPK 0.2105445 0.6717039 0.3673430 24.0563894 7.222854
#> σ_slope_RespPK 14.0819091 52.8632704 23.3796024 7.2228537 918.086053
#>
#> *********************************************
#> Determinant, condition numbers and D-criterion
#> ***********************************************
#>
#> Determinant: 488075482
#> D-criterion: 12.1916
#> Condition number (fixed effects): 11769929
#> Condition number (variance components): 89.6053
#>
#> ***************************************
#> Parameters estimation
#> ***************************************
#>
#> Parameter Value SE RSE(%)
#> μ_V 100.0000000 148.84145795 148.84146
#> μ_Cl 10.0000000 3.81302408 38.13024
#> μ_ka 0.2500000 0.38298792 153.19517
#> ω²_V 0.1000000 0.20391433 203.91433
#> ω²_Cl 0.0800000 0.04002444 50.03055
#> ω²_ka 0.2000000 0.25333384 126.66692
#> σ_inter_RespPK 1.0000000 0.20413127 20.41313
#> σ_slope_RespPK 0.2236068 0.03344128 14.95539The optimized five-sample design improves the D-criterion, but precision on the typical parameters remains limited.
A sixth occasion is added in the washout phase (72 h and 168 h after the last dose), as in the extended Tim Waterhouse scenario.
observation_times_6 = c(
5,
week_h * 8 - 1,
week_h * c(12, 24),
lastDose + 72,
lastDose + 168
)
samplingConstraintsRespPK6 = SamplingTimeConstraints(
outcome = "RespPK",
initialSamplings = observation_times_6,
samplingsWindows = list(
c(0, 6),
c(week_h * 8 - 1, week_h * 8),
c(week_h * 12 - 1, week_h * 12),
c(week_h * 24 - 1, week_h * 24),
c(lastDose, lastDose + 168),
c(lastDose + 167.5, lastDose + 168)
),
numberOfTimesByWindows = c(1, 1, 1, 1, 1, 1),
minSampling = c(1, 1, 1, 1, 1, 1)
)armOptim6 = Arm(
name = "ArmOptim6",
size = 12,
administrations = list(administrationRespPK),
samplingTimes = list(
SamplingTimes(outcome = "RespPK", samplings = observation_times_6)
),
samplingTimesConstraints = list(samplingConstraintsRespPK6)
)
designOptim6 = Design(
name = "designOptim6",
arms = list(armOptim6)
)
optimizationPopFIM6 = Optimization(
name = "Optimization 6 samples",
modelEquations = modelEquations,
modelParameters = modelParameters,
modelError = modelError,
optimizer = "PSOAlgorithm",
optimizerParameters = list(
maxIteration = 10,
populationSize = 50,
personalLearningCoefficient = 2.05,
globalLearningCoefficient = 2.05,
seed = 1234,
showProcess = FALSE
),
designs = list(designOptim6),
fimType = "population",
outputs = list("RespPK")
)
optimizationPopFIM6 = run(optimizationPopFIM6)show(optimizationPopFIM6)
#>
#> --- Optimal design ---
#>
#> Arms name Number of subjects Outcome
#> 1 ArmOptim6 12 RespPK
#> Dose
#> 1 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000
#> Sampling times
#> 1 (0.47, 1343.25, 2015.23, 4031.04, 4072.97, 4199.68)
#>
#> ***************************************
#> Population Fisher Matrix
#> ***************************************
#>
#> μ_V μ_Cl μ_ka ω²_V ω²_Cl ω²_ka σ_inter_RespPK σ_slope_RespPK
#> μ_V 0.006737006 -0.02518152 -0.8183397 0.000000 0.00000000 0.00000000 0.0000000 0.00000
#> μ_Cl -0.025181524 1.27906374 -0.6138614 0.000000 0.00000000 0.00000000 0.0000000 0.00000
#> μ_ka -0.818339690 -0.61386140 631.9699997 0.000000 0.00000000 0.00000000 0.0000000 0.00000
#> ω²_V 0.000000000 0.00000000 0.0000000 189.113534 26.42121406 17.43957937 8.5638376 43.73878
#> ω²_Cl 0.000000000 0.00000000 0.0000000 26.421214 681.66835435 0.09813172 0.4955298 58.22905
#> ω²_ka 0.000000000 0.00000000 0.0000000 17.439579 0.09813172 65.00424488 2.3414104 34.21810
#> σ_inter_RespPK 0.000000000 0.00000000 0.0000000 8.563838 0.49552976 2.34141042 26.5972921 18.50211
#> σ_slope_RespPK 0.000000000 0.00000000 0.0000000 43.738779 58.22905274 34.21809952 18.5021058 991.23066
#>
#> ***************************************
#> Fixed effects (μ)
#> ***************************************
#>
#> μ_V μ_Cl μ_ka
#> μ_V 0.006737006 -0.02518152 -0.8183397
#> μ_Cl -0.025181524 1.27906374 -0.6138614
#> μ_ka -0.818339690 -0.61386140 631.9699997
#>
#> ***************************************
#> Variance components (ω², γ², σ)
#> ***************************************
#>
#> ω²_V ω²_Cl ω²_ka σ_inter_RespPK σ_slope_RespPK
#> ω²_V 189.113534 26.42121406 17.43957937 8.5638376 43.73878
#> ω²_Cl 26.421214 681.66835435 0.09813172 0.4955298 58.22905
#> ω²_ka 17.439579 0.09813172 65.00424488 2.3414104 34.21810
#> σ_inter_RespPK 8.563838 0.49552976 2.34141042 26.5972921 18.50211
#> σ_slope_RespPK 43.738779 58.22905274 34.21809952 18.5021058 991.23066
#>
#> *********************************************
#> Determinant, condition numbers and D-criterion
#> ***********************************************
#>
#> Determinant: 8.4334e+11
#> D-criterion: 30.95639
#> Condition number (fixed effects): 122775.1
#> Condition number (variance components): 38.93445
#>
#> ***************************************
#> Parameters estimation
#> ***************************************
#>
#> Parameter Value SE RSE(%)
#> μ_V 100.0000000 13.93528758 13.935288
#> μ_Cl 10.0000000 0.92863143 9.286314
#> μ_ka 0.2500000 0.04380307 17.521230
#> ω²_V 0.1000000 0.07447811 74.478111
#> ω²_Cl 0.0800000 0.03849369 48.117107
#> ω²_ka 0.2000000 0.12658928 63.294638
#> σ_inter_RespPK 1.0000000 0.19645472 19.645472
#> σ_slope_RespPK 0.2236068 0.03239345 14.486792optArms = prop(optimizationPopFIM6, "optimisationAlgorithmOutputs")$optimalArms
optSamplings = prop(prop(optArms[[1]], "samplingTimes")[[1]], "samplings")
optSamplingsWith these settings, the typical values V,
Cl, and ka are expected to reach RSE below 20%
after optimization (variance components are usually less precisely
estimated).
Pediatric studies allow only a few samples per subject. Starting from a feasible sparse design, PFIM supports identifiability checks, FIM-based precision, D-optimal timing under windows, and comparison of alternative schedules.
This scenario follows Tim Waterhouse’s PopED workshop material
(optimal-design-2023), re-expressed with PFIM objects. For
object mapping, variance versus standard deviation, and cross-checks
with PopED, see
vignette("Vignette_05_gcse_pkpd", package = "PFIM"),
section PFIM, PopED, and cross-tool equivalence.