Hello!
I’am trying to use OptimalDesign and have a few questions:
- Is-it possible to perform only design evaluation? Currently what I am doing is:
- Defining in the bounds arguments as many intervals as sampling times I want in my design,
- set those intervals to
(t_j..t_j)where t_j are the samplings times, - and constraint having one observation per interval:
obs_times = [0.5,1,2,6,24,36,72,120]
obs_dictionnary = Dict( (t..t) => 1 for t in obs_times )
- Morevoer, I add the argument
max_iter = 0indesign()to only evaluate the inital design
=> is there a better way to do that?
- I wanted to investigate the warfarin example presented in Nyberg, J., Bazzoli, C., Ogungbenro, K., Aliev, A., Leonov, S., Duffull, S., … & Mentré, F. (2015). Methods and software tools for design evaluation in population pharmacokinetics–pharmacodynamics studies. British journal of clinical pharmacology, 79(1), 6-17., nevertheless, with the proportional error model, the execution stops with an error:
PosDefException: matrix is not positive definite; Cholesky factorization failed.while it runs with an additive error model (but still doesn’t give the same SE as in PFIM and PopED for the fixed effects parameters, which suggests there’s a problem with my implementation in Pumas (or with the Pumas source code?)).
Here is my code, can you have a look at it and explore:- why I have an error with the proportional error model
- why the SE with additive error model are not the same as in PFIM and PopED
model_warfarin = @model begin
@param begin
tvcl ∈ RealDomain(; lower = 0.0)
tvvc ∈ RealDomain(; lower = 0.0)
tvKa ∈ RealDomain(; lower = 0.0)
ω²cl ∈ RealDomain(;lower = 0.0)
ω²vc ∈ RealDomain(;lower = 0.0)
ω²Ka ∈ RealDomain(;lower = 0.0)
σ² ∈ RealDomain(; lower = 0.0)
end
@random begin
ηcl ~ Normal(0.0, sqrt(ω²cl))
ηvc ~ Normal(0.0, sqrt(ω²vc))
ηKa ~ Normal(0.0, sqrt(ω²Ka))
end
@pre begin
CL = tvcl * exp(ηcl)
Vc = tvvc * exp(ηvc)
Ka = tvKa * exp(ηKa)
end
@dynamics Depots1Central1
@derived begin
cp = @. Central / Vc
# dv ~ @. Normal(cp, sqrt( cp^2 * σ² )) # Proportional
dv ~ @. Normal(cp, sqrt( σ²)) # Additive
end
end
parameters = (;
tvcl = 0.15,
tvvc = 8,
tvKa = 1,
ω²cl = 0.07,
ω²vc = 0.02,
ω²Ka = 0.6,
σ² = 0.01
)
elementary_design_warfarin = Subject(; events = DosageRegimen(70), observations = (; dv = nothing))
obs_times = [0.5,1,2,6,24,36,72,120]
obs_dictionnary = Dict( (t..t) => 1 for t in obs_times )
# Design evaluation
dec_warfarin = decision(
model_warfarin,
elementary_design_warfarin,
parameters;
type = :observation_times,
bounds = obs_dictionnary,
minimum_offset = 0.25,
init = obs_times,
subject_multiples = 32
);
eval_result_warfarin = design(dec_warfarin; optimality = :doptimal, auto_initialize = false, max_iter = 0 )
SE = sqrt.(diag(inv(eval_result_warfarin.optimalfim)))
paramValues = [0.15, 8, 1, 0.07, 0.02, 0.6, 0.01]
df = DataFrame(Parameter = ["Cl", "V", "Ka", "omega^2_Cl", "omega^2_V", "omega^2_Ka", "b^2" ],
Value = paramValues,
SE = SE,
RSE = 100 * SE./paramValues
)
Thank you very much in advance!!
Best regards
Lucie
