using OptimalDesign, Pumas, Random, StableRNGs
Random.seed!(Pumas.default_rng(), 500)
rng = StableRNG(20);
Two-compartment PK optimal design
1 Case study description
In this case study, the problem of designing the sample collection times of a phase II clinical trial for an adult population is considered. Phase I data is assumed to be available and the population PK (popPK) model developed in phase I is assumed to be a two-compartment model with random effects on all the PK parameters. In phase I, the safety of the drug had been established using a wide dose range administered to healthy volunteers. Regular blood draws enabled the characterization of the popPK model and some population parameter estimates \(\hat{\theta}\). The goal of the phase II study is to quantify the drug’s efficacy for multiple dosing regimens. The study is assumed to be randomized and double-blinded.
1.1 Model definition and parameter values
First, the environment is setup for the experiment as below:
The optimal design process begins by defining the model and parameter estimates obtained using phase I data. The Pumas PK model and assumed parameter estimates are shown below.
model = @model begin
@param begin
tvcl ∈ RealDomain(lower = 0.001)
tvvc ∈ RealDomain(lower = 0.001)
tvka ∈ RealDomain(lower = 0.001)
tvq ∈ RealDomain(lower = 0.001)
tvvp ∈ RealDomain(lower = 0.001)
σ_prop ∈ RealDomain(lower = 0.001)
σ_add ∈ RealDomain(lower = 0.001)
end
@random η ~ MvNormal(Diagonal([0.25, 0.25, 0.25, 0.25, 0.49]))
@pre begin
CL = tvcl * exp(η[1])
Vc = tvvc * exp(η[2])
Q = tvq * exp(η[3])
Vp = tvvp * exp(η[4])
Ka = tvka * exp(η[5])
end
@vars begin
Cp = Central / Vc
end
@dynamics begin
Depot' = -Ka * Depot
Central' =
Ka * Depot - (CL / Vc) * Central - (Q / Vc) * Central + (Q / Vp) * Peripheral
Peripheral' = (Q / Vc) * Central - (Q / Vp) * Peripheral
end
@derived begin
DV ~ @. CombinedNormal(Cp, σ_add, σ_prop)
end
end
params = (tvcl = 25, tvvc = 50, tvka = 0.5, tvq = 31, tvvp = 48, σ_prop = 0.1, σ_add = 0.2);Note that the inter-individual variability parameter \(\Omega\) has been hard-coded into the model. The model includes the following population parameters (@param block):
- tvcl: typical value of clearance.
- tvvc: typical value of the volume of distribution in the central compartment.
- tvka: typical value of the absorption rate constant.
- tvq: typical value of the inter-compartmental clearance.
- tvvp: typical value of the volume of distribution in the peripheral compartment.
- \(\sigma_\text{prop}\): proportional error rate.
- \(\sigma_\text{add}\): additive error.
The estimates of these parameters, fitted using Phase I data, are defined as a named tuple after the model definition.
Furthermore, in the @random block, 5 random effects are defined, one for each PK parameter, with mean 0 and a diagonal covariance matrix. In the @pre block, the individual PK parameters are defined by combining the typical values and random effects. In the @dynamics block, the differential equations of the standard 2 compartment PK model are defined.
Finally, a combined additive and proportional residual error model is use in the @derived block. The model predicts the individual drug concentration \(C_p = \text{Central}/\text{Vc}\), where \(\text{Central}\) and \(\text{Vc}\) are the amount of drug in the plasma and volume of the central compartment, respectively.
1.2 Defining a population
To perform optimal design, a population needs to be defined to contain:
- The covariates of the target population of phase II, e.g. sex and weight.
- The dosing regimens to be administered to the subjects in phase II.
The covariates are simulated from an assumed distribution of covariates in the target population. The number of dosing regimens and dose values were determined based on the safety analysis and simulation results from phase I. Subjects were randomly assigned an arm in the study.
The code below creates a population of 250 subjects evenly split between 5 dose levels (doses). This done by means of a vector of integers, dose_arm. The doses are all orally administered at the beginning of the study, time = 0.
doses = [25, 50, 100, 200, 400] .* 1e3
dose_arm = vcat([fill(i, 50) for i = 1:length(doses)]...)
nSubjects = length(dose_arm)
population = map(1:nSubjects) do i
dose = doses[dose_arm[i]]
return Subject(
id = i,
events = DosageRegimen(dose, time = 0),
observations = (; DV = fill(missing, 25)),
time = 0:24.0,
)
end;1.3 Time windows
Next, we specify the time window in which the sample times are to be selected. The code below sets the start of the time window at the beginning of the experiment (t0 = 0.0) and the end at 24h (t_end = 24.0). Time window definitions can also be provided through date-time structs (example in the appendix). Next, we specify the number of measurements (5) to be taken inside each time window for each subject.
t0 = 0.0
t_end = 24.0
time_window = Dict((t0 .. t_end) => 5);1.4 Optimization
Next, we group all of the components of an optimal design problem in a decision_info variable. This variable groups all of the information relevant to the design optimization, such as model, population, parameter estimates, time windows, initial design and a number of other options, e.g. a minimum offset of 0.25 between every 2 consecutive samples is used below to prevent samples from being too close to each other. Initial sample times are provided through init to initialize the optimization; here the same five sample times are used for every subject. If init is not specified, the sample times are initialized randomly within the provided time windows, and the random times need not satisfy the minimum offset.
decision_info = decision(
model,
population,
params;
type = :observation_times,
bounds = time_window,
minimum_offset = 0.25,
init = [[0.5, 2.0, 6.0, 12.0, 24.0] for _ in population],
)Population decision - sample times
Population size: 250
Samples per subject: 25
Number of time windows: 1An important verification is that of identifiability (Van Noort et al. 2024). This refers to the ability of the data to enable identification of the parameters. If a parameter is not identifiable, it cannot be estimated (Shivva et al. 2013), which is indicated by a very large variance after fitting the model. In the code below, the function analyze_identifiability computes the FIM, and verifies if any parameters are locally non-identifiable. The parameters of the model used in this problem are all identifiable.
id = analyze_identifiability(model, population, params, nsol = 5);[ Info: Computing FIM. [ Info: Running identifiability analysis. [ Info: Unable to prove local non-identifiability. The model may be identifiable.
Next, we perform the design optimization using the design function as shown below. The only mandatory inputs are the decision_info variable and the type of optimality criteria to use, e.g. :doptimal, :aoptimal, or :toptimal. Here, the :aoptimal optimality criteria is used. An optional limit of 10 iterations is also set. Additionally, the optional keyword argument param_weights is used to scale the diagonal of the inverse FIM by the inverse of the squares of the parameter values. This prevents the optimizer from prioritizing a certain parameter just because it is numerically larger than others. This parameter scaling method assumes that none of the parameters is exactly 0.
itrf = Pumas.toidentitytransform(model.param)
paramVec = Pumas.TransformVariables.inverse(itrf, params)
objectiveType = :aoptimal
result = design(
decision_info;
param_weights = [1 / p^2 for p in paramVec],
optimality = objectiveType,
max_iter = 1,
)
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.9.0.
Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 2000
Number of nonzeros in Lagrangian Hessian.............: 0
Total number of variables............................: 1250
variables with only lower bounds: 0
variables with lower and upper bounds: 1250
variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 1000
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 1000
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 6.5650440e+00 0.00e+00 6.59e-01 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 6.3666646e+00 0.00e+00 3.60e-02 -1.1 1.21e-01 - 9.91e-01 1.00e+00f 1
Number of Iterations....: 1
(scaled) (unscaled)
Objective...............: 6.3666646373590492e+00 6.3666646373590492e+00
Dual infeasibility......: 3.6026909669323898e-02 3.6026909669323898e-02
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 3.5905996769330195e-01 3.5905996769330195e-01
Overall NLP error.......: 3.5905996769330195e-01 3.5905996769330195e-01
Number of objective function evaluations = 2
Number of objective gradient evaluations = 2
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 2
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations = 0
Total seconds in IPOPT = 1553.514
EXIT: Maximum Number of Iterations Exceeded.Optimal schedule:
Time points per subject: 5
Subjects: 250
Initial objective: 5.863282615242788
Optimal objective: 6.366664637359049
Size of FIM: (7, 7)1.5 Getting the optimal design and tweaking it
After the optimization, another utility function optimaltimes can be used to extract the final design as a vector, as in the listing below. Each component of the vector is itself a vector of the sample times associated with a specific subject. These values can be manually altered according to a context, as is done on the second line, which rounds all values to integers. Furthermore, the initial time values used in the optimization, which is stored in the decision_info variable, can be directly accessed, as demonstrated in the third line.
opt_times = optimaltimes(result)
rounded_times = map(t -> round.(t), opt_times)
init_times = [dec.t.t for dec in decision_info];1.6 Post-processing
Using the decision_info variable, the utility function check_feasibility can be used to check if a certain design is feasible or not according to the constraints in decision_info. Feasibility checking checks that the sample times are:
- In the pre-defined time windows.
- Unique.
- Non-negative.
- Sorted within subjects.
- Satisfying the optional minimum offset between any two consecutive samples.
The output of the check_feasibility function is a boolean, true (feasible) or false (not feasible). This function is useful if manual adjustments are made to the design and the user would like to confirm that the new design is feasible.
feasible = check_feasibility(decision_info, init_times);Another utility function, named_fim can be used to calculate the FIM as a matrix with columns and rows labeled by the parameter names. The observation time points provided are used to calculate the FIM.
nf = named_fim(model, population, params)┌ 7×7 DimArray{Float64, 2} ┐ ├──────────────────────────┴───────────────────────── dims ┐ ↓ X Categorical{String} ["tvcl", …, "σ_add"] Unordered, → Y Categorical{String} ["tvcl", …, "σ_add"] Unordered └──────────────────────────────────────────────────────────┘ ↓ → "tvcl" "tvvc" … "σ_prop" "σ_add" "tvcl" 171.956 -30.5941 -7917.01 -132.283 "tvvc" -30.5941 6.41924 958.382 21.1257 "tvka" 1483.53 -339.151 -90123.6 -1453.92 "tvq" 23.155 -4.73128 -147.672 -9.21629 "tvvp" -53.5414 10.4892 … 1171.89 32.4158 "σ_prop" -7917.01 958.382 9.68344e5 9731.65 "σ_add" -132.283 21.1257 9731.65 14023.3
The CRLBs at the optimal design can be directly retrieved by the function crlb as in the code block below. The output is a table with the parameters’ names, absolute CRLBs, and relative CRLBs.
lowBounds = crlb(model, params, result)7×3 DataFrame
| Row | Parameter | crlb | relative_crlb |
|---|---|---|---|
| Symbol | Float64 | Float64 | |
| 1 | tvcl | 0.427526 | 0.017101 |
| 2 | tvvc | 1.85351 | 0.0370702 |
| 3 | tvka | 0.0128201 | 0.0256403 |
| 4 | tvq | 1.54108 | 0.0497124 |
| 5 | tvvp | 1.24942 | 0.0260295 |
| 6 | σ_prop | 0.00475344 | 0.0475344 |
| 7 | σ_add | 0.504338 | 2.52169 |
To evaluate a design, one can use the utility function evaluate_design which returns an object containing the FIM and the objective function value. The listing below includes an example using the initial random sample times, and the design after optimization. In this example, the optimization reduced the objective function value from 26.9 to 23.6.
initialPerf = evaluate_design(decision_info, population, init_times, objectiveType)
println(initialPerf)
optimPerf = evaluate_design(decision_info, population, opt_times, objectiveType)FIM: 7x7 matrix
Objective value: 7.594
Optimality criterion: average-variance (aoptimal)
Parameter weights: nothing
FIM: 7x7 matrix
Objective value: 7.809
Optimality criterion: average-variance (aoptimal)
Parameter weights: nothingThe optimal design specifies the sample times for each subject. Depending on the size of the population, that amount of information may be overwhelming to communicate to staff and patients. Therefore, in practice, it is common to summarize the design, looking for a few schedules (or even only one) that are applied to all subjects. This simplifies the final design, ideally without sacrificing too much estimation precision.
As an example, in the code below, medianSummary summarizes the designs by calculating the medians of each sample across subjects. It is therefore a vector with the 5 time instants of the samples.
# Organize times in matrix
opt_times_mat = reduce(hcat, opt_times)
# calculate median of each sample time across subjects
medianSummary = median.(eachrow(opt_times_mat))
# round medians to the nearest multiple of 0.5
halfHourMeds = round.(Int, medianSummary * 2) / 2;halfHourMeds then also simplifies the global schedule obtained by rounding the values to the nearest multiple of 0.5. The reason for this is that the optimization process is performed with 64-bit floating point numbers. This can lead to time definitions that are unrealistically precise, considering that this information will define instructions for patients and medical staff. These procedures result in the global schedule of \([3, 7.5, 12, 16.5, 21]\) h.
Another way to summarize the optimal design is to visualize it by plotting either the exact sample times or some summary statistics of each sample. The first sub-plot in Figure 1 shows the distribution of each of the 5 samples per subject in this case study using a box plot per sample. The second sub-plot shows the median sample times overlaid on a plot of simulated concentrations using the model and design population.
In the code below, the optimal design is used to simulate data and re-estimate the parameters and infer the standard errors to double check the output of the optimal design procedure. Simulations are generated using the simobs function. The model is then fitted to the data simulated. Additionally, after fitting, inference is performed, using the infer function, to obtain asymptotic estimates of the standard errors and 95% confidence intervals of the population parameter estimates. This information can also be summarized in a table format using the coeftable function.
# Simulation observations
simulations =
simobs(model, population, params, obstimes = 0:0.5:24; rng, simulate_error = false)
simPop = Subject.(simulations)
# Fit model to simulated data
fitModel = fit(model, simPop, params, FOCE();)
# Inference
inferModel = infer(fitModel)
coeftable(inferModel)[ Info: Checking the initial parameter values. [ Info: The initial negative log-likelihood and its gradient are finite. Check passed. Iter Function value Gradient norm 0 3.974384e+04 9.872397e+03 * time: 0.0338749885559082 1 3.028313e+04 8.597023e+03 * time: 14.146250009536743 2 2.391304e+04 9.935850e+03 * time: 21.24119997024536 3 2.275421e+04 4.897548e+03 * time: 25.815321922302246 4 2.127510e+04 1.773770e+03 * time: 30.394495964050293 5 2.100858e+04 1.401892e+03 * time: 34.95059108734131 6 2.089575e+04 7.138312e+02 * time: 39.4462411403656 7 2.071221e+04 5.453079e+02 * time: 44.03213596343994 8 2.058151e+04 5.918361e+02 * time: 48.65134000778198 9 2.053281e+04 2.280995e+02 * time: 53.281095027923584 10 2.052773e+04 2.438685e+01 * time: 57.81851005554199 11 2.052758e+04 2.357119e+01 * time: 62.25895309448242 12 2.052741e+04 1.933372e+01 * time: 66.74560904502869 13 2.052723e+04 1.024736e+01 * time: 71.24774193763733 14 2.052714e+04 4.203283e+00 * time: 75.75803399085999 15 2.052712e+04 2.809475e+00 * time: 80.2187750339508 16 2.052711e+04 2.667765e+00 * time: 84.5896589756012 17 2.052711e+04 2.440767e+00 * time: 88.96452713012695 18 2.052711e+04 1.716868e+00 * time: 93.39512991905212 19 2.052711e+04 7.812909e-01 * time: 97.80201005935669 20 2.052710e+04 4.526789e-01 * time: 102.26627993583679 21 2.052711e+04 3.585443e-01 * time: 106.61963605880737 22 2.052710e+04 2.878888e-01 * time: 110.95225191116333 23 2.052711e+04 1.976883e-01 * time: 115.29017400741577 24 2.052711e+04 1.230349e-01 * time: 119.62473797798157 25 2.052710e+04 6.122840e-02 * time: 123.97920298576355 26 2.052710e+04 3.180907e-02 * time: 128.28963899612427 27 2.052711e+04 1.985414e-02 * time: 132.5398781299591 28 2.052711e+04 1.631333e-02 * time: 136.77780103683472 29 2.052711e+04 1.319303e-02 * time: 141.01375913619995 30 2.052711e+04 6.803166e-03 * time: 145.24089694023132 31 2.052711e+04 1.779449e-03 * time: 149.46376991271973 32 2.052711e+04 8.307552e-03 * time: 153.6277940273285 33 2.052711e+04 1.553246e-01 * time: 157.95451593399048 34 2.052710e+04 1.442339e-03 * time: 162.29864597320557 35 2.052711e+04 6.888507e-04 * time: 166.4821491241455 [ Info: Calculating: variance-covariance matrix. [ Info: Done.
7×8 DataFrame
| Row | parameter | description | constant | estimate | se | relative_se | ci_lower | ci_upper |
|---|---|---|---|---|---|---|---|---|
| String | String | Bool | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | tvcl | tvcl | false | 24.9019 | 0.842176 | 0.0338198 | 23.2512 | 26.5525 |
| 2 | tvvc | tvvc | false | 49.7912 | 1.94898 | 0.039143 | 45.9712 | 53.6111 |
| 3 | tvka | tvka | false | 0.480199 | 0.0202717 | 0.0422152 | 0.440468 | 0.519931 |
| 4 | tvq | tvq | false | 22.4786 | 2.28078 | 0.101465 | 18.0083 | 26.9488 |
| 5 | tvvp | tvvp | false | 62.3912 | 3.04567 | 0.0488157 | 56.4218 | 68.3606 |
| 6 | σ_prop | σ_prop | false | 0.0100056 | 0.0016413 | 0.164039 | 0.00678868 | 0.0132225 |
| 7 | σ_add | σ_add | false | 0.0378932 | 0.0180722 | 0.476924 | 0.0024724 | 0.073314 |
Often the optimized design, or its summarized version, is used only as a starting point, and manual adjustments are made, for example to adjust for constraints that were difficult to incorporate in the optimization process. Let the manually adjusted design be \([6, 9, 12, 16.5, 19]\) h. Two designs (summarized and manual) can be compared in a table as shown below.
# Altered summarized design
altSchedule = [6, 9, 12, 16.5, 19]
# Designs
halfHourDesign = fill(halfHourMeds, nSubjects)
altDesign = fill(altSchedule, nSubjects)
# Design evaluations
performances = map([halfHourDesign, altDesign]) do des
evaluate_design(decision_info, population, des, objectiveType)
end
# CRLBs
CRLBs = map(p -> p.FIM |> inv |> diag .|> sqrt, performances)
# Normalizations
itrf = Pumas.toidentitytransform(model.param)
paramVec = Pumas.TransformVariables.inverse(itrf, params)
normCRLBs = map(lb -> lb ./ paramVec, CRLBs)
# Percent change
percentChangeCRLB = @. (normCRLBs[2] / normCRLBs[1] - 1) * 100
df = DataFrame(
Parameter = Pumas._coeftable(params).key,
Design_1_CRLBs = normCRLBs[1],
Design_2_CRLBs = normCRLBs[2],
Percent_change = round.(percentChangeCRLB; digits = 2),
)7×4 DataFrame
| Row | Parameter | Design_1_CRLBs | Design_2_CRLBs | Percent_change |
|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | |
| 1 | tvcl | 0.0170999 | 0.0209747 | 22.66 |
| 2 | tvvc | 0.0366766 | 0.125333 | 241.72 |
| 3 | tvka | 0.0257016 | 0.0467235 | 81.79 |
| 4 | tvq | 0.04947 | 0.127874 | 158.49 |
| 5 | tvvp | 0.0257444 | 0.0851005 | 230.56 |
| 6 | σ_prop | 0.047909 | 0.0441754 | -7.79 |
| 7 | σ_add | 2.72039 | 1.51197 | -44.42 |
The resulting table includes the relative CRLBs for the parameters estimated (5 fixed effects as well as additive and proportional components of the residual error). The CRLBs are divided by the respective parameter values to represent lower bounds on the relative standard error. Intuitively, the exact CRLBs indicate the highest achievable estimation precision (lowest standard error) from the data acquired during an experiment1. The manual adaptation of the summarized design worsened the objective value from 49.9 to 72. Finally, the poorer theoretical performance of the altered design is noticed by a median 21.2% increase in the relative CRLB.
1.7 Simulation experiment
When selecting the number of subjects in a study, it is useful to understand how the CRLB and the distribution of simulation-based standard error estimates depend on the population size. In this case study, we calculate the CRLB and simulation-based standard error estimates for multiple population sizes, starting from a base population of 50 subjects (10 per dose levels). The base population is then gradually expanded, emulating the participation of more individuals in a study. For each population size, the CRLB is estimated and the simulation-based analysis is performed 100 times to get a distribution for the simulated-based standard error estimates. The following steps detail the exact procedure used to generate Figure 2 and Figure 3. The script used is included in the appendix.
- Starting from the base population, sequentially create populations of increasing size, always appending 20 new subjects randomly sampled from the original population. Repeat the remaining steps for each population size.
- Evaluate the design (using the
evaluate_designfunction) to estimate the CRLBs for the given population size. - Perform 100 simulations of observations for the current population using the reference model and parameter values.
- Fit the model to each of the 100 simulated datasets.
- For each of the 100 fits, obtain standard error estimates of the population parameters.
- Plot the relative CRLB estimates and distribution of simulation-based relative standard error estimates for each population size.
Figure 2 and Figure 3 show a clear asymptotic decrease of the approximate relative CRLBs and relative standard errors, as expected from statistical theory. A narrowing of the range of errors across simulations is also observed with more subjects. Additionally with smaller populations, the lower bound property of CRLBs appears to be violated for some parameters. This is possible due to the approximation error when calculating the FIM using the first-order method. With more subjects, such violations occur less often. These plots help to visualize how the precision of parameter estimates changes as the population size increase. This can be helpful when selecting a population size to achieve a target precision level.
2 Appendix A - time windows with DateTime
The following listing includes an example of time windows definition using the DateTime format. The first window goes from 8 am of January 1st, 2024, to 16 pm of that day. And the second window spans the same interval of the next day.
using Dates
t0 = DateTime(2024, 1, 1, 8)
t_end = DateTime(2024, 1, 1, 16)
time_window = Dict((t0 .. t_end) => 3, (t0 .. t_end) + Day(1) => 3)3 Appendix B - Convergence simulation
using OptimalDesign, Pumas, Random, StableRNGs
Random.seed!(Pumas.default_rng(), 500)
rng = StableRNG(20)
### Initial definitions
# Model
model = @model begin
@param begin
tvcl ∈ RealDomain(lower = 0.001)
tvvc ∈ RealDomain(lower = 0.001)
tvka ∈ RealDomain(lower = 0.001)
tvq ∈ RealDomain(lower = 0.001)
tvvp ∈ RealDomain(lower = 0.001)
σ_prop ∈ RealDomain(lower = 0.001)
σ_add ∈ RealDomain(lower = 0.001)
end
@random η ~ MvNormal(Diagonal([0.25, 0.25, 0.25, 0.25, 0.49]))
@pre begin
CL = tvcl * exp(η[1])
Vc = tvvc * exp(η[2])
Q = tvq * exp(η[3])
Vp = tvvp * exp(η[4])
Ka = tvka * exp(η[5])
end
@vars begin
Cp = Central / Vc
end
@dynamics begin
Depot' = -Ka * Depot
Central' =
Ka * Depot - (CL / Vc) * Central - (Q / Vc) * Central + (Q / Vp) * Peripheral
Peripheral' = (Q / Vc) * Central - (Q / Vp) * Peripheral
end
@derived begin
DV ~ @. CombinedNormal(Cp, σ_add, σ_prop)
end
end
# Estimates
params = (tvcl = 25, tvvc = 50, tvka = 0.5, tvq = 31, tvvp = 48, σ_prop = 0.1, σ_add = 0.2)
# Population
nSubjects = 50
doses = [10, 30, 100, 300, 850, 2500] .* 1e3
dose_arm = [ceil(Int, i / 9) for i = 1:nSubjects]
population = map(1:nSubjects) do i
_dose = doses[dose_arm[i]]
return Subject(
id = i,
events = DosageRegimen(_dose, time = 0),
observations = (; DV = fill(missing, 25)),
time = 0:24.0,
)
end
# Time windows
t0 = 0.0
t_end = 24.0
time_window = Dict((t0 .. t_end) => 5)
# Decisions
decisions = decision(
model,
population,
params;
type = :observation_times,
bounds = time_window,
minimum_offset = 0.25,
init = [[0.5, 2.0, 6.0, 12.0, 24.0] for _ in population],
)
### Optimization
itrf = Pumas.toidentitytransform(model.param)
paramVec = Pumas.TransformVariables.inverse(itrf, parameters)
objectiveType = :aoptimal
result = design(
decisions;
param_weights = [1 / p^2 for p in paramVec],
optimality = objectiveType,
max_iter = 10,
)
# Optimal design
opt_times = optimaltimes(result)
### Summarize and simplify design
# Round medians of samples to nearest multiple of 0.5
samplesPerSubject = 5
medianSummary = [median(getindex.(opt_times, i)) for i = 1:samplesPerSubject]
halfHourMeds = round.(Int, medianSummary * 2) / 2
halfHourDesign = [halfHourMeds for _ = 1:nSubjects]
### Run analysis
# Setup for batch of iterations
iters = 0:1
indices = 1:length(iters)
popIncrement = 20
repeatSim = 10
# Create progressively larger populations
allPopulations = [population]
allDecisions = [decisions]
for _ in iters
randomNewIDs = rand(1:nSubjects, popIncrement)
push!(allPopulations, vcat(allPopulations[end], population[randomNewIDs]))
push!(allDecisions, vcat(allDecisions[end], decisions[randomNewIDs]))
end
popsAndDecisions = Dict(:pops => allPopulations, :decisions => allDecisions)
# Convergence analysis
_crlb = [params |> length |> zeros for _ in indices]
df = [[DataFrame() for _ = 1:repeatSim] for _ in indices]
for (index, convIter) in enumerate(iters)
# Population for current iteration. Samples are always on halfHourMeds
newPopulation = map(
currentSubs -> OD.set_observation_times(currentSubs, halfHourMeds),
popsAndDecisions[:pops][convIter+1],
)
newDecisions = popsAndDecisions[:decisions][convIter+1]
# CRLB for expanded design
_crlb[index] .=
evaluate_design(
newDecisions,
newPopulation,
[halfHourMeds for _ = 1:length(newPopulation)],
objectiveType,
).FIM |>
inv |>
diag .|>
sqrt
# Standard errors through inference from multiple simulations
dataFrames = map(1:repeatSim) do s
simPop = simobs(model, newPopulation, params; rng) .|> Subject
# Fit model to simulated data
fitModel =
fit(model, simPop, params, FOCE(); optim_options = (show_trace = false,))
# Inference
inferModel = infer(fitModel)
# Estimate standard errors from table
return coeftable(inferModel)
end
df[index] .= dataFrames
end3.1 References
Shivva, V., J. Korell, I. Tucker, and S. Duffull. 2013. “An Approach for Identifiability of Population Pharmacokinetic–Pharmacodynamic Models.” CPT: Pharmacometrics & Systems Pharmacology 2 (6): 1–9. https://doi.org/10.1038/psp.2013.25.
Van Noort, M., M. Ruppert, J. DeJongh, et al. 2024. “Navigating the Landscape of Parameter Identifiability Methods: A Workflow Recommendation for Model Development.” CPT: Pharmacometrics & Systems Pharmacology 13 (7): 1170–79. https://doi.org/10.1002/psp4.13148.
Footnotes
However, the CRLBs calculated here are only approximate given that the FIM was approximated using the first-order method. Even then they are useful for design comparison↩︎