Error while running simobs using Bayesian fit

Estimation bayesian
1

Hi I ran the following model

BAYESIAN_PKPD = @model begin
    @param begin
        tvKTV1 ~ LogNormal(log(0.4), 0.5)
        tvKTV2 ~ LogNormal(log(0.1), 0.5)
        tvecmo ~ Normal(-0.7, 0.5)
        tvCL1 ~ LogNormal(log(0.9), 0.5)
        tvCL2 ~ LogNormal(log(0.8), 0.5)
        tvcircuit ~ LogNormal(log(1.2), 0.5)
        tvVc1 ~ LogNormal(log(1.8), 0.25)
        tvVc2 ~ LogNormal(log(3.0), 0.25)
        σ_prop_1 ~ truncated(Cauchy(0, 0.4), 0, 1.0)
        σ_prop_2 ~ truncated(Cauchy(0, 0.2), 0, 1.0)
        σ_prop_3 ~ truncated(Cauchy(0, 0.2), 0, 1.0)
        σ_add ~ truncated(Cauchy(0, 80), 0, Inf)
        ω ∈ Constrained(
            MvNormal(zeros(3), Diagonal([1.0,0.2,1.0])),
            lower=zeros(3),
            upper=fill(Inf, 3),
            init=ones(3),
        )
    end

    @random begin
        ηstd ~ MvNormal(I(3))
    end
    @covariates WT ECMO_DAYS IS_CIRCUIT_CHANGE site T
    @pre begin
        # compute the η from the ηstd
        # using lower Cholesky triangular matrix
        η = ω .* ηstd

        # PK parameters

        siteeffect_1 = if site==1 
            σ_prop_1^2 
        elseif site==2
            σ_prop_2^2
        else 
            σ_prop_3^2
        end 
        siteeffect_2 = site==2 ? σ_add^2 : 0.0

        tvKTV = site == 1 ? tvKTV1 : tvKTV2*(ECMO_DAYS/6)^tvecmo
        KTV    = tvKTV*exp(η[1])
        TVCL = site == 1 ? tvCL1 * (WT/10)^0.75 : tvCL2* (WT/10)^0.75  * tvcircuit^IS_CIRCUIT_CHANGE  
        CL     = TVCL * exp(η[2])*(1-exp(-KTV*t))
        TVVc = site == 1 ? tvVc1 * (WT/10)^1 : tvVc2 * (WT/10)^1
        Vc     = TVVc * exp(η[3])
    end

    @dynamics begin
        Central' = -(CL/Vc) * Central
    end

    @derived begin
        CP = @. Central / Vc
        CONC = @. Normal(CP, sqrt(CP^2*siteeffect_1 + siteeffect_2))
    end
end


  
  param = (
  tvKTV1=0.4,
  tvKTV2 = 0.1,
  tvecmo = -0.7,
  tvCL1 = 1.0,
  tvCL2 = 1.0,
  tvcircuit = 1.2,
  tvVc1 = 1.8,
  tvVc2 = 3.6,
  ω=[1.0, 0.2,1.0],
  σ_prop_1=0.4,
  σ_prop_2 = 0.2,
  σ_prop_3 = 0.2,
  σ_add = 80
)


result_PKPD   = @time fit(BAYESIAN_PKPD, estimation, param, Pumas.BayesMCMC(nsamples=2000, nadapts=1000,nchains =4,target_accept = 0.6,ensemblealg = EnsembleThreads()))

when I try to run the following code

posterior_sims = simobs(
    result_PKPD;
    samples=200,
    simulate_error=true
)

I get the following error

┌ Warning: Unrecognized keyword arguments found. Future versions will error.
│ The only allowed keyword arguments to `solve` are:
│ (:dense, :saveat, :save_idxs, :tstops, :tspan, :d_discontinuities, :save_everystep, :save_on, :save_start, :save_end, :initialize_save, :adaptive, :abstol, :reltol, :dt, :dtmax, :dtmin, :force_dtmin, :internalnorm, :controller, :gamma, :beta1, :beta2, :qmax, :qmin, :qsteady_min, :qsteady_max, :qoldinit, :failfactor, :calck, :alias_u0, :maxiters, :callback, :isoutofdomain, :unstable_check, :verbose, :merge_callbacks, :progress, :progress_steps, :progress_name, :progress_message, :timeseries_errors, :dense_errors, :weak_timeseries_errors, :weak_dense_errors, :calculate_errors, :initializealg, :alg, :save_noise, :delta, :seed, :alg_hints, :kwargshandle, :trajectories, :batch_size, :sensealg, :advance_to_tstop, :stop_at_next_tstop, :default_set, :second_time, :prob_choice)
│ 
│ See https://diffeq.sciml.ai/stable/basics/common_solver_opts/ for more details.
│ 
│ Set kwargshandle=KeywordArgError for an error message.
│ Set kwargshandle=KeywordArgSilent to ignore this message.
└ @ DiffEqBase /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/DiffEqBase/VN57T/src/solve.jl:835

Could you please let me know what does this error mean ?

2

Is this the full error? Can you post the complete stack trace?

3

This is what is shown but it repeats infinite number of time through the terminal and I have to stop it usint ctrl+C

4

Just one thing to clarify is that it appears only when I use differential equation while using t in

CL     = TVCL * exp(η[2])*(1-exp(-KTV*t))

since Pumas does not allow closed form solution if t is recalled in @pre block but If I used t as a covariate T like the following:

df = @chain df begin
    @transform :T = :TIME
end

estimation = read_pumas(df,
    observations=[:CONC], 
    id=:ID,
    time=:TIME,
    evid=:evid,
    mdv = :mdv,
    amt=:AMT,
    rate = :RATE,
    cmt = :cmt,
    covariates = [:WT,:ECMO_DAYS,:IS_CIRCUIT_CHANGE, :site,:T]
)

Thin I modify the equation in @pre block to

    @covariates WT IS_CIRCUIT_CHANGE ECMO_DAYS site T
@pre begin
CL     = TVCL * exp(η[2])*(1-exp(-KTV*T))
end

So I can run the closed form solution. If I did so I do not get the error above and simobs runs ok.

5

So, it seems the error is related to using the differential equation.

FYI, both differential equation and closed form provide the same estimates at the end. But I am uncertain if using time as covariate to circumvent Pumas requirement for using differential equation if t is recalled in @pre block is mathematically correct or not.

6

If your ODE parameters depend nonlinearly on time, there is no general closed form solution for the ODE (only some special cases) so you are probably better off relying on an ODE solver. Covariates in Pumas are interpolated using a piece-wise constant function so using time as a covariate is not mathematically correct and can introduce significant numerical inaccuracy. So I would suggest using an ODE for that problem instead of treating time as a covariate. Cool hack though!

Regarding the particular error, I haven’t seen it before but I can look into it and get back.

7

@ahmed.salem I just tried your model with a simple subject and simulated with t and I did not get any errors. Note that I used the initial parameters for simobs


julia> df = Subject(
           id = 1,
           events = DosageRegimen(100),
           covariates = (WT = 70, ECMO_DAYS = 2, IS_CIRCUIT_CHANGE = 1, site = 2),
           observations = (CONC = nothing,),
       )
Subject
  ID: 1
  Events: 1
  Observations: CONC: (n=0)
  Covariates: WT, ECMO_DAYS, IS_CIRCUIT_CHANGE, site

julia> posterior_sims = simobs(
           BAYESIAN_PKPD, 
           df, 
           param
       )                            
SimulatedObservations
  Simulated variables: CP, CONC
  Time: 0.0:1.0:24.0
1 Like
8

@ahmed.salem I could reproduce the warning and root caused it. It’s a mild warning that can be ignored and should disappear in the next release.

1 Like
9

@mohamed82008 alright, Thanks for looking into it !