OptimalDesign for design evaluation + warfarin example

How-to
1

Hello!

I’am trying to use OptimalDesign and have a few questions:

  1. 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 = 0 in design() to only evaluate the inital design

=> is there a better way to do that?

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

2

Hi, Lucie!

Sorry to hear you are running into problems. Assuming you are using Pumas 2.6, there are some utility functions that can help you. To evaluate the design, you can use the function evaluate_design():

ev = evaluate_design(
  [dec_warfarin],
  [elementary_design_warfarin],
  [obs_times],
  :doptimal
)

It returns a DesignEvaluation object, containing the FIM, objective value, parameter weights, and the type of objective.

Also, to determine the Cramer-Rao lower bounds (CRLBs) on the standard errors, here I would usually recommend that you use crlb() function, instead of SE = sqrt.(diag(inv(eval_result_warfarin.optimalfim))).

SE = crlb(
  model_warfarin,
  parameters,
  eval_result_warfarin
)

It returns a DataFrame with the CRLBs and their normalized values for each parameter.

However, by adding using the verbose keyword argument in eval_result_warfarin = design(dec_warfarin; optimality = :doptimal, auto_initialize = false, max_iter = 0, verbose = true), the optimizer outputs “EXIT: Problem has inconsistent variable bounds or constraint sides.”, which is caused by defining time windows without length (i.e. (t..t)).

So, to avoid using the results of a problematic optimization, you can determine the CRLBs through SE = sqrt.(diag(inv(ev.FIM))), where ev was returned by evaluate_design() as above.

And, when it comes to the error and matching your results against the study, I’ll need to look further into it.

2 Likes
3

Hi!

Thank you very much for your help!

I didn’t know about the evaluate_design() and crlb() functions, thanks!

I still obtain different results from PopED and PFIM for RSE for the additive error model and not positive definite matrix for the proportional error model…

4

Hi, Lucie.

Through simulations, it was verified that cp = @. Central / V is sometimes null, which creates numerical problems. Changing the residual error definition to dv ~ @. Normal(cp, sqrt( cp * σ² + 1e-5)) avoids that.

And regarding the mismatching results, could you please provide the PFIM and PopED scripts you used to reproduce the first example in the paper you mentioned?

5

Indeed, thank you!

of course, here’s the PopED code:

library(PopED)

ff <- function(model_switch,xt,parameters,poped.db){
  ##-- Model: One comp first order absorption
  with(as.list(parameters),{
    y=xt
    y=(DOSE*Favail*KA/(V*(KA-CL/V)))*(exp(-CL/V*xt)-exp(-KA*xt))
    return(list(y=y,poped.db=poped.db))
  })
}

sfg <- function(x,a,bpop,b,bocc){
  ## -- parameter definition function
  parameters=c(CL=bpop[1]*exp(b[1]),
               V=bpop[2]*exp(b[2]),
               KA=bpop[3]*exp(b[3]),
               Favail=bpop[4],
               DOSE=a[1])
  return(parameters)
}

feps <- function(model_switch,xt,parameters,epsi,poped.db){
  ## -- Residual Error function
  ## -- Proportional
  returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db))
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]
  y = y*(1+epsi[,1])
  
  return(list(y=y,poped.db=poped.db))
}


## -- Define initial design  and design space
poped.db <- create.poped.database(ff_file="ff",
                                  fg_file="sfg",
                                  fError_file="feps",
                                  bpop=c(CL=0.15, V=8, KA=1.0, Favail=1),
                                  notfixed_bpop=c(1,1,1,0),
                                  d=c(CL=0.07, V=0.02, KA=0.6),
                                  sigma=0.01,
                                  groupsize=32,
                                  xt=c( 0.5,1,2,6,24,36,72,120),
                                  minxt=0,
                                  maxxt=120,
                                  a=70)

## evaluate initial design
resWarfarin_PopED <- PopED::evaluate_design(poped.db)
resWarfarin_PopED$rse

and with PFIM, RSE are the same

6

Hi!
Do you have any news about with issue with design evaluation?
Thank you very much in advance!

7

Hi Lucie,

Happy new year! We haven’t yet identified the source of discrepancy. To do that, we will need to dig into the PopED source code a bit more to understand the differences in implementation. I am hoping to get some time to dig into the PopED and Pumas implementations if not this week then next week. We appreciate your patience while we try to figure this out.

Regards,
Mohamed

8

Hi Mohamed,

Happy new year to you too!

Thank you for the update. Please let me know if there is any additional information or support I can provide to assist you.

Looking forward to your findings!

Best regards

1 Like