L2 Correlations

parssh.mehta · August 5, 2023, 7:41pm

Hi I am trying to use L2 correlation, since the pk and pd concentrations are taken from the same sample. The simulations work and it also starts the estimation but then fails with an error. Can you please help me fix this?

Attached model code and failed estimation error

true_model_pk_delay_pd_sim = @model begin
    @param begin
        tvcl      ∈ RealDomain(lower=0.0)
        tvvc      ∈ RealDomain(lower=0.0)
        tvq       ∈ RealDomain(lower=0.0)
        tvvp      ∈ RealDomain(lower=0.0)
        tvka      ∈ RealDomain(lower=0.0)
        tvkin     ∈ RealDomain(lower=0.0)
        tvkout    ∈ RealDomain(lower=0.0)
        tvImax    ∈ RealDomain(lower=0.0)
        tvIC50    ∈ RealDomain(lower=0.0)
        Ω         ∈ PDiagDomain(8)
        Σ         ∈ PSDDomain(2)
    end

    @random begin
        η      ~ MvNormal(Ω)
    end

    @pre begin
        CL     = tvcl * exp(η[1])
        Vc     = tvvc * exp(η[2])
        Q      = tvq * exp(η[3])
        Vp     = tvvp * exp(η[4])
        Ka     = tvka * exp(η[5])
        Kin    = tvkin * exp(η[6])
        Kout   = tvkout * exp(η[7])
        Imax   = tvImax
        IC50   = tvIC50 * exp(η[8])

        σ_prop_pk = Σ[1]
        σ_prop_pd = Σ[2]
    end

    @init begin
      Response = Kin/Kout
    end

    @vars begin
      Central_conc = (Central/Vc)*1000
    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
      Response'   =  Kin*(1-((Imax * Central_conc)/(IC50 + Central_conc))) - Kout*Response
    end

    @derived begin
      Cp  = @. Central_conc
      DV  ~ @. Normal(Cp, Cp*σ_prop_pk)

      Resp    = @. Response
      DV_Resp ~ @. Normal(Resp, Resp*σ_prop_pd)
    end
end

##Simulation
final_params = (tvcl      = 11.3,
                tvvc      = 37.6,
                tvq       = 31.4,
                tvvp      = 47.8,
                tvka      = 0.8,
                tvkin     = 20,
                tvkout    = 0.8,
                tvImax    = 1,
                tvIC50    = 2000,
                Ω         = Diagonal([0.0625,0.09,0.0625,0.09,0.09,0.16,0.09,0.09]),
                Σ         = [0.1  0.02
                              0.02 0.2])

N=30
ev0  = DosageRegimen(500, time=0, cmt=1)
protocol_sampling_time = [0.20,1.33,1.66,1.99,2.95,6.96,7.30,7.63,19.22]
pop0 = map(i -> Subject(id=i, events=ev0, time=protocol_sampling_time), 1:N)
pop0_df = DataFrame(pop0)

Random.seed!(123)
sim_base  = simobs(true_model_pk_delay_pd_sim, 
                    pop0,   
                    final_params)
df_base_sim   = DataFrame(sim_base)


##Estimation
initial_params = (  tvcl      = 10.3,
                    tvvc      = 39.6,
                    tvq       = 35.4,
                    tvvp      = 50.8,
                    tvka      = 0.8,
                    tvkin     = 20,
                    tvkout    = 0.7,
                    tvImax    = 0.8,
                    tvIC50    = 2100,
                    Ω         = Diagonal([0.0625,0.09,0.0625,0.09,0.09,0.16,0.09,0.09]),
                    Σ          = [0.1  0.02
                                  0.02 0.2])

read_fit_df = read_pumas(df_base_sim,
                        observations = [:DV, :DV_Resp])

fit_foce = fit(true_model_pk_delay_pd_sim, read_fit_df, initial_params, 
                        FOCE(), ensemblealg=EnsembleThreads())

Error

[ 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     2.810998e+03     2.300546e+04
 * time: 2.5987625122070312e-5
ERROR: MethodError: no method matching iterate(::Nothing)
Closest candidates are:
  iterate(::Union{LinRange, StepRangeLen}) at range.jl:872
  iterate(::Union{LinRange, StepRangeLen}, ::Integer) at range.jl:872
  iterate(::T) where T<:Union{Base.KeySet{<:Any, <:Dict}, Base.ValueIterator{<:Dict}} at dict.jl:712
  ...
Stacktrace:
  [1] indexed_iterate(I::Nothing, i::Int64)
    @ Base ./tuple.jl:91
  [2] _unidentified_exception(A::PDMats.PDMat{Float64, Matrix{Float64}}, i::Int64, j::Int64, d::Int64, k::Symbol)
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:3510
  [3] _check_zero_gradient(vparam::Vector{Float64}, g::Vector{Float64}, fixedtrf::TransformVariables.TransformTuple{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, Pumas.PDiagTransform, Pumas.PSDTransform}}})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:3475
  [4] (::Pumas.var"#619#621"{FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, Pumas.PDiagTransform, Pumas.PSDTransform}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}})(state::Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:4047
  [5] #588
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:3000 [inlined]
  [6] update!(tr::Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}}, iteration::Int64, f_x::Float64, grnorm::Float64, dt::Dict{Any, Any}, store_trace::Bool, show_trace::Bool, show_every::Int64, callback::Pumas.var"#588#590"{Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, Pumas.var"#619#621"{FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, Pumas.PDiagTransform, Pumas.PSDTransform}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}}}, trace_simplex::Bool)
    @ Optim /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/tP8PJ/src/utilities/update.jl:23
  [7] update!
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/tP8PJ/src/utilities/update.jl:11 [inlined]
  [8] trace!(tr::Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}}, d::NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, state::Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}}, iteration::Int64, method::Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, options::Optim.Options{Float64, Pumas.var"#588#590"{Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, Pumas.var"#619#621"{FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, Pumas.PDiagTransform, Pumas.PSDTransform}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}}}}, curr_time::Float64)
    @ Optim /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/tP8PJ/src/multivariate/solvers/first_order/bfgs.jl:191
  [9] optimize(d::NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, method::Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, options::Optim.Options{Float64, Pumas.var"#588#590"{Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, Pumas.var"#619#621"{FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, Pumas.PDiagTransform, Pumas.PSDTransform}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}}}}, state::Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}})
    @ Optim /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/tP8PJ/src/multivariate/optimize/optimize.jl:50
 [10] optimize
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/tP8PJ/src/multivariate/optimize/optimize.jl:36 [inlined]
 [11] DefaultOptimizeFN
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:2993 [inlined]
 [12] _fit(m::PumasModel{(tvcl = 1, tvvc = 1, tvq = 1, tvvp = 1, tvka = 1, tvkin = 1, tvkout = 1, tvImax = 1, tvIC50 = 1, Ω = 8, Σ = 3), 8, (:Depot, :Central, :Peripheral, :Response), ParamSet{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, PSDDomain{Matrix{Float64}}}}}, var"#5#12", var"#6#13", Nothing, var"#7#14", ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, ODEFunction{false, ModelingToolkit.ODEFunctionClosure{var"#8#15", var"#9#16"}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, var"#10#17", var"#11#18", ModelingToolkit.ODESystem}, population::Vector{Subject{NamedTuple{(:DV, :DV_Resp), Tuple{Vector{Union{Missing, Float64}}, Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(), Tuple{}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, param::NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{Float64, Float64, Float64, Float64, Float64, Int64, Float64, Float64, Int64, Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, approx::FOCE, ensemblealg::EnsembleThreads, optimize_fn::Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, fixedparamset::ParamSet{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, PSDDomain{Matrix{Float64}}}}}, fixedparam::NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{Float64, Float64, Float64, Float64, Float64, Int64, Float64, Float64, Int64, Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, checkidentification::Bool, diffeq_options::NamedTuple{(), Tuple{}}, init_vrandeffsorth::Vector{Vector{Float64}}, verbose::Bool)
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:4062
 [13] __fit
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:3804 [inlined]
 [14] fit(m::PumasModel{(tvcl = 1, tvvc = 1, tvq = 1, tvvp = 1, tvka = 1, tvkin = 1, tvkout = 1, tvImax = 1, tvIC50 = 1, Ω = 8, Σ = 3), 8, (:Depot, :Central, :Peripheral, :Response), ParamSet{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, PSDDomain{Matrix{Float64}}}}}, var"#5#12", var"#6#13", Nothing, var"#7#14", ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, ODEFunction{false, ModelingToolkit.ODEFunctionClosure{var"#8#15", var"#9#16"}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, var"#10#17", var"#11#18", ModelingToolkit.ODESystem}, _population::Vector{Subject{NamedTuple{(:DV, :DV_Resp), Tuple{Vector{Union{Missing, Float64}}, Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(), Tuple{}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Symbol}}, Vector{Float64}}}, param::NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :tvka, :tvkin, :tvkout, :tvImax, :tvIC50, :Ω, :Σ), Tuple{Float64, Float64, Float64, Float64, Float64, Int64, Float64, Float64, Int64, Diagonal{Float64, Vector{Float64}}, Matrix{Float64}}}, approx::FOCE; optim_alg::Nothing, optim_options::Nothing, optimize_fn::Nothing, constantcoef::NamedTuple{(), Tuple{}}, omegas::Tuple{}, ensemblealg::EnsembleThreads, checkidentification::Bool, diffeq_options::NamedTuple{(), Tuple{}}, init_randeffs::Nothing, init_vrandeffsorth::Nothing, verbose::Bool)
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/zfm1h/src/estimation/likelihoods.jl:3768
 [15] top-level scope
    @ ~/data/code/Methodology-Project/model/L2 Correlation/l2_correlation.jmd:161

storopoli · August 7, 2023, 11:10am

This is probably an indexing issue, since Σ is a matrix. Try:

σ_prop_pk = Σ[1, 1]
σ_prop_pd = Σ[2, 2]

benjaminrich · August 7, 2023, 1:03pm

I don’t think this is the right approach. If you want the errors for DV and DV_Resp to be correlated, you would need to use a multivariate normal distribution for the error vector, something like this:

ϵ ~ @. MvNormal(Σ)

Cp  = @. Central_conc
DV  = @. Cp*(1 + ϵ[1])

Resp    = @. Response
DV_Resp = @. Resp*(1 + ϵ[2])

I’m not sure if this is the exact way to write it in Pumas, maybe @storopoli can help?

storopoli · August 7, 2023, 7:28pm

@andreasnoack and I tried to do a minimum working example of such correlated error models. It would look something like this:

mod = @model begin
  @param begin
    σ_PK in RealDomain(; lower=0)
    σ_PD in RealDomain(; lower=0)
    ρ in RealDomain()
  end

  @derived begin
    DV ~ @. MvNormal([σ_PK^2 σ_PK*σ_PD*ρ; σ_PK*σ_PD*ρ σ_PD^2])
  end
end

Of course you would need to supply your DV as a vector of length-2.

However, this does not currently works in Pumas. We’ll support this type of model in a future Pumas release.

parssh.mehta · August 8, 2023, 1:04am

Thanks @storopoli , that works.

jogagobburu · August 8, 2023, 1:10pm

@parssh.mehta Hope you are doing well. Glad to see the model specification part is taken care of.

I am unable to appreciate the scientific basis for estimating a correlation between the random effects of PK and PD, albeit being bioanalyzed from the same sample. Even if one were to assume PD is say another concentration of some endogenous biomarker- can you please explain the scientific rationale for attempting a correlation of residual errors of PK and PD?
J

parssh.mehta · August 9, 2023, 7:36pm

Hi @jogagobburu I’m doing great, hope you are doing well too.

The current project that I am working on is impact of assessing the influence of systematic errors on Pop PK/PD model. The scientific basis for imputing this correlation is that historically the assay variability contributed the highest to the residual error, but that is not the case these days with improvement in analytical techniques. But still the residual error reported in literature is greater than the assay error and since we are modeling clinical trial perturbations (dose variability, dosing time, sampling time, and bioanalytical errors) assuming both scenarios of PK and PD samples collected at the same time and also collected at different times. If PK and PD are collected at the same time, in the same sample, we would assume that if any errors in dose variability, dosing time, sampling time, and bioanalytical errors can impact both PK and PD and not one. This can contribute that these errors are not independent from PK or PD if they are collected at the same time. Moreover, If we stay at the same sample, any misspecification in PK can propagate to PD and this shows the dependence between these dynamics. Thus we wanted to see if if there is a correlation which can bias our estimation of the PK and PD parameters from the various perturbations.

Happy to hear your thoughts on this any feedback is appreciated.

Parssh

mjaber · August 9, 2023, 7:53pm

Hi @storopoli,

Thanks for investigating this. In simulation exercise, things are promising with this implementation, not yet validated.

model = @model begin
    @param begin
        θ in VectorDomain(6)
        Ω in PDiagDomain(6)
        Σ in PSDDomain(2)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        CL = θ[1] * exp(η[1])
        V = θ[2] * exp(η[2])
        ka = θ[3] * exp(η[3])

        kin = θ[4] * exp(η[4])
        kout = θ[5] * exp(η[5])
      
        IC = θ[6] * exp(η[6])
    end
    @init begin
        y3 =kin/kout
    end

    @dynamics begin
        y1' = - ka * y1 
        y2' = ka * y1 - (CL/V) * y2 
        y3'  = kin * (1 - y2/(y2 + IC)) - kout * y3  #'
    end

    @derived begin
        cp = @. log(y2/V + 0.001) 
        pd = @. log(y3 + 0.01)  

        vec := @. MvLogNormal([cp, pd], Σ)
        dv1 ~ @. vec[1]
        dv2 ~ @. vec[2]
    end
end

dose = 500 
param = (θ = [8.0, 25.0, 0.8, 20.0, 0.2, 1.3], 
        Ω = Diagonal(ones(6) * 0.1), 
        Σ = [0.2 0.1; 0.1 0.3], )

ev = DosageRegimen(dose)       
sub = [Subject(;id=i, events=ev) for i in 1:50]
sim = simobs(model, sub, param, obstimes=0:24) |> DataFrame

However, as you mentioned read_pumas need to be modified to handle this type of data. I also want to ask about your first post, in pre chunk. What does pumas do under the hood? How it can handle the covariance or calculate the correlation since it sounds like you are using standard deviation rather than variance term?

This implementation can be helpful not only in PK/PD systems, but with parent/metabolite and plasma/urine when we do not assume independence (or at least independence is not valid assumption).

storopoli · August 10, 2023, 12:52pm

We don’t know how this would be implemented.
But there are some issues with the code that you’ve presented.

Every ~ in the derived block is some sort of “likelihood statement”. Hence you would need to use ~ with any distribution. Hence it is analogous to a probabilistic statement (since you are familiar to PPLs).
= assignments are only deterministic statements, thus, cannot have any distribution on the right hand side.

so this would probably become something like this:

DVs ~ @ MvLogNormal([cp, pd], Σ)
dv1 = @. vec[1]
dv2 = @. vec[2]

And in the read_pumas you would need to pass a column of length-2 vector elements to the observations kwarg:

pop = read_pumas(
  df;
  observations=[:DVs]
  ...
)

mjaber · August 11, 2023, 3:44am

Thanks Jose! Does it matter when we conduct simulations?

jogagobburu · August 11, 2023, 10:53am

@parssh.mehta Thank you. I would be curious to see what your research teaches us. Still wrapping my head around what practical problem such a correlation of random effects will solve (ignoring the pain of simultaneous PKPD modeling in the real world). One potential application would be for performing more realistic clinical trial simulations (CTS) - however handling residual error in CTS is by itself a challenging topic. Once you have some findings, we can discuss if the juice is worth the squeeze!
On a personal note, thank you for using Pumas for such innovative research - it is built for power users such as yourself!
J

farmacyst · October 31, 2024, 9:36pm

@storopoli I’m facing the same issue with my estimation fit(). I also get the MethodError: no method matching iterate(::Nothing)

Here is my model file

using Pumas
using NCA
using Random
using PumasUtilities
using NCAUtilities
using Distributions
using DataFramesMeta
using Dates
using AlgebraOfGraphics

Random.seed!(1234)

brooks_model = @model begin

    @metadata begin
        desc = "BROOKS MODEL FOR PEDIATRIC FLUDARABINE"
        timeu = u"hr"                                   # time in hours
        amtu = u"mg"                                    # amount in miligrams
    end
    
    @param begin
        θ1 in RealDomain(lower = 0, upper = 4)          # 1 CL
        θ2 in RealDomain(lower = 0, upper = 0.002)      # 2 CRCLEFF
        θ3 in RealDomain(lower = 0, upper = 12)         # 3 V
        θ4 in RealDomain(lower = 0, upper = 2)          # 4 Q
        θ5 in RealDomain(lower = 0, upper = 9)          # 5 V2
        θ6 in RealDomain(lower = 0, upper = 0.006)      # 6 KIN
        θ7 in RealDomain(lower = -0.5 , upper = 0)      # 7 KIN TIME EFFECT
        θ8 in RealDomain(lower = 0, upper = 0.1)        # 8 KOUT
        θ11 in RealDomain(lower = 0, upper = 3.7)       # FMAT shape

        Ω in PSDDomain(2)                               # For building a correlation matric
        ΩKin in RealDomain(lower = 0)                   # Kin is not correlated to V & CL unlike in Ivaturi Model

        σ1 in RealDomain(lower = 0, upper = 1)          # F-ara-A concentration in Plasma
        σ2 in RealDomain(lower = 0, upper = 1)          # F-ara-ATP amount in Intracellular Compartment
    end

    @random begin
        ηCL ~ Normal(0, Ω[1])
        ηV ~ Normal(0,Ω[2])
        ηKin ~ Normal(0,ΩKin)
    end

    @covariates WT HT AGE CRCL SEX

    @pre begin
        HTM = HT/100
        BMI = WT/(HTM^2)
        BSA = (WT^0.425 * HT^0.725) * 0.007184

        FFM = 
        if SEX == 1
            (0.88 + ((1-0.88)/(1+(AGE/13.4)^(-12.7)))) * 
            ((9270.0 * WT)/(6680.0 + (216.0 * BMI)))
        else
            (1.11 + ((1 - 1.11)/(1 + (AGE / 7.1)^(-1.1)))) * 
            ((9270.0 * WT)/(8780.0 + (244.0 * BMI)))
        end

        FMAT = (1-exp(-AGE * θ11))
        BIO = 0.781
        CRCLunits = CRCL * 16.666667 * 1.73 / BSA
        
        CRCLEFF = 
        if CRCLunits > 150
            1 + θ2 * (150-100)
        else
            1 + θ2 * (CRCLunits-100)
        end
        
        TVCL = θ1 * (FFM/12)^0.75 * CRCLEFF * FMAT
        CL = (TVCL/BIO) * exp(ηCL)
        TVV = θ3 * (FFM/12)
        V = (TVV/BIO) * exp(ηV)
        Q = (θ4/BIO) * (FFM/12)^0.75
        V2 = (θ5/BIO) * (FFM/12)
        TVKIN = θ6
        KIN = TVKIN * exp(ηKin)

        KINt = 
        if t > 0
            KIN * (t + 0.001)^(θ7)
        else
            KIN
        end

        KOUT = θ8
    end
    
    @init begin                                        # Since we are taking log error model, we do not want conc at t=0 to be Inf
        TRIP = 1
        CENT = 1
    end

    @vars begin
        Cc = CENT/(V/1000)
        trip = TRIP * 1000
    end

    @dynamics begin
        CENT' = -(CL/V)*CENT + (Q/V2)*PERI - (Q/V)*CENT 
        PERI' =              - (Q/V2)*PERI + (Q/V)*CENT
        TRIP' = KINt*(CENT/V) - KOUT*TRIP
    end

    @derived begin
        cp = @. Cc
        int_cp = @. trip
        log_dv ~ @. Normal(log(abs(cp)), σ1)
        log_int_dv ~ @. Normal(log(abs(int_cp)), σ2)
    end

end

params = (;
    θ1 = 3.14,
    θ2 = 0.00185,
    θ3 = 11.5,
    θ4 = 1.46,
    θ5 = 8.48,
    θ6 = 0.00599,
    θ7 = -0.415,
    θ8 = 0.0949,
    θ11 = 3.64,
    Ω = [
        0.013689 0.019881;
        0.019881 0.091204
        ],
    ΩKin = 0.44997264,
    σ1 = 0.1296,
    σ2 = 0.512656,
)

# Covariates for Population
function choose_covariates()
    WT =rand((TruncatedNormal(23.6, 12, 7.3, 135)))
    HT = rand(TruncatedNormal(121.7, 32, 67, 186)) 
    AGE = rand(TruncatedNormal(8, 7, 2, 22))
    CRCL = rand(TruncatedNormal(137, 47, 57, 231))
    SEX = rand(Binomial(1,0.7))
    
    return (; WT, HT, AGE, CRCL, SEX)
end

# Calculation BSA because traditional dosing of Fludarabine is based on BSA
function bsa(WT,HT)
    return (WT^0.425 * HT^0.725 * 0.007184)
end

population = map(
    i -> Subject(;
        id = i,
        covariates = choose_covariates(),
        events = DosageRegimen(bsa(choose_covariates().WT,choose_covariates().HT)*40; 
            ii = 24, addl = 3, rate = 28, route = NCA.IVInfusion
        ),
        observations = (; dv = nothing),
    ),
    1:10,
)
pop = DataFrame(population)

@time sim = simobs(brooks_model, population, params; obstimes = 0:0.5:23.9)
simdf = DataFrame(sim)
CSV.write("data.csv", simdf)

# Fitting Model
fit_pumas = read_pumas(
    simdf,
    observations = [:log_dv],
    id = :id,
    time = :time,
    amt = :amt,
    evid = :evid,
    covariates = [:WT,:AGE,:HT,:CRCL,:SEX],
    rate = :rate,
    cmt = :cmt
)

fit(brooks_model,
    fit_pumas,params,
    FOCEI();
    constantcoef=(θ6=0.00599,θ7=-0.415,θ8=0.0949,σ1=0.1296,σ2=0.512656)
)

Error:

[ 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.044205e+02     3.446869e+01
 * time: 0.0011739730834960938
ERROR: MethodError: no method matching iterate(::Nothing)

Closest candidates are:
  iterate(::Union{LinRange, StepRangeLen})
   @ Base range.jl:880
  iterate(::Union{LinRange, StepRangeLen}, ::Integer)
   @ Base range.jl:880
  iterate(::T) where T<:Union{Base.KeySet{<:Any, <:Dict}, Base.ValueIterator{<:Dict}}
   @ Base dict.jl:698
  ...

Stacktrace:
  [1] indexed_iterate(I::Nothing, i::Int64)
    @ Base /Applications/Pumas-2.5.1.app/Contents/Resources/julia/lib/julia/sys.pumas.dylib:-1
  [2] _unidentified_exception(A::PDMats.PDMat{Float64, Matrix{Float64}}, i::Int64, j::Int64, d::Int64, k::Symbol)
    @ Pumas /Users/runner/work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Zl9W3/src/estimation/likelihoods.jl:3677
  [3] _check_zero_gradient(vparam::Vector{Float64}, g::Vector{Float64}, fixedtrf::TransformVariables.TransformTuple{NamedTuple{(:θ1, :θ2, :θ3, :θ4, :θ5, :θ6, :θ7, :θ8, :θ11, :Ω, :ΩKin, :σ1, :σ2), Tuple{TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}, TransformVariables.ScaledShiftedLogistic{Float64}, Pumas.PSDTransform, TransformVariables.ShiftedExp{true, Float64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}}}})
    @ Pumas /Users/runner/work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Zl9W3/src/estimation/likelihoods.jl:3642
  [4] (::Pumas.var"#622#626"{FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:θ1, :θ2, :θ3, :θ4, :θ5, :θ6, :θ7, :θ8, :θ11, :Ω, :ΩKin, :σ1, :σ2), Tuple{TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}, TransformVariables.ScaledShiftedLogistic{Float64}, Pumas.PSDTransform, TransformVariables.ShiftedExp{true, Float64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}})(state::Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}})
    @ Pumas /Users/runner/work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Zl9W3/src/estimation/likelihoods.jl:4290
  [5] #591
    @ /Users/runner/work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Zl9W3/src/estimation/likelihoods.jl:3086 [inlined]
  [6] update!(tr::Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}}, iteration::Int64, f_x::Float64, grnorm::Float64, dt::Dict{Any, Any}, store_trace::Bool, show_trace::Bool, show_every::Int64, callback::Pumas.var"#591#593"{Pumas.DefaultOptimizeFN{Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, Pumas.var"#622#626"{FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:θ1, :θ2, :θ3, :θ4, :θ5, :θ6, :θ7, :θ8, :θ11, :Ω, :ΩKin, :σ1, :σ2), Tuple{TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}, TransformVariables.ScaledShiftedLogistic{Float64}, Pumas.PSDTransform, TransformVariables.ShiftedExp{true, Float64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}}}, trace_simplex::Bool)
    @ Optim /Users/runner/work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/V8ZEC/src/utilities/update.jl:23
  [7] update!
    @ /Users/runner/work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/V8ZEC/src/utilities/update.jl:11 [inlined]
  [8] trace!(tr::Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}}, d::NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, state::Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}}, iteration::Int64, method::Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, options::Optim.Options{Float64, Pumas.var"#591#593"{Pumas.DefaultOptimizeFN{Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, Pumas.var"#622#626"{FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:θ1, :θ2, :θ3, :θ4, :θ5, :θ6, :θ7, :θ8, :θ11, :Ω, :ΩKin, :σ1, :σ2), Tuple{TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}, TransformVariables.ScaledShiftedLogistic{Float64}, Pumas.PSDTransform, TransformVariables.ShiftedExp{true, Float64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}}}}, curr_time::Float64)
    @ Optim /Users/runner/work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/V8ZEC/src/multivariate/solvers/first_order/bfgs.jl:191
  [9] optimize(d::NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, method::Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, options::Optim.Options{Float64, Pumas.var"#591#593"{Pumas.DefaultOptimizeFN{Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, Pumas.var"#622#626"{FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:θ1, :θ2, :θ3, :θ4, :θ5, :θ6, :θ7, :θ8, :θ11, :Ω, :ΩKin, :σ1, :σ2), Tuple{TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}, TransformVariables.ScaledShiftedLogistic{Float64}, Pumas.PSDTransform, TransformVariables.ShiftedExp{true, Float64}, Pumas.ConstantTransform{Float64}, Pumas.ConstantTransform{Float64}}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}}}}, state::Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}})
    @ Optim /Users/runner/work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Optim/V8ZEC/src/multivariate/optimize/optimize.jl:50
 [10] DefaultOptimizeFN
    @ /Users/runner/work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Zl9W3/src/estimation/likelihoods.jl:3096 [inlined]
 [11] _fit(m::PumasModel{(θ1 = 1, θ2 = 1, θ3 = 1, θ4 = 1, θ5 = 1, θ6 = 1, θ7 = 1, θ8 = 1, θ11 = 1, Ω = 3, ΩKin = 1, σ1 = 1, σ2 = 1), 3, (:CENT, :PERI, :TRIP), (:HTM, :BMI, :BSA, :FFM, :FMAT, :BIO, :CRCLunits, :CRCLEFF, :TVCL, :CL, :TVV, :V, :Q, :V2, :TVKIN, :KIN, :KINt, :KOUT), ParamSet{NamedTuple{(:θ1, :θ2, :θ3, :θ4, :θ5, :θ6, :θ7, :θ8, :θ11, :Ω, :ΩKin, :σ1, :σ2), Tuple{RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Float64, Float64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Float64, Float64}, RealDomain{Float64, Int64, Float64}, RealDomain{Int64, Float64, Float64}, RealDomain{Int64, Float64, Float64}, PSDDomain{PDMats.PDMat{Float64, Matrix{Float64}}}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Int64, Float64}}}}, var"#705#715", Pumas.TimeDispatcher{var"#706#716", var"#708#718"}, Nothing, var"#710#720", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, SciMLBase.NullParameters, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.ODEFunctionClosure{var"#711#721", var"#712#722"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, var"#713#723", var"#714#724", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, population::Vector{Subject{NamedTuple{(:log_dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:WT, :AGE, :HT, :CRCL, :SEX), Tuple{Float64, Float64, Float64, Float64, Int64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, param::NamedTuple{(:θ1, :θ2, :θ3, :θ4, :θ5, :θ6, :θ7, :θ8, :θ11, :Ω, :ΩKin, :σ1, :σ2), Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Matrix{Float64}, Float64, Float64, Float64}}, approx::FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}, ensemblealg::EnsembleThreads, optimize_fn::Pumas.DefaultOptimizeFN{Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, fixedparamset::ParamSet{NamedTuple{(:θ1, :θ2, :θ3, :θ4, :θ5, :θ6, :θ7, :θ8, :θ11, :Ω, :ΩKin, :σ1, :σ2), Tuple{RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Float64, Float64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Int64, Float64}, Pumas.ConstDomain{Float64}, Pumas.ConstDomain{Float64}, Pumas.ConstDomain{Float64}, RealDomain{Int64, Float64, Float64}, PSDDomain{PDMats.PDMat{Float64, Matrix{Float64}}}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, Pumas.ConstDomain{Float64}, Pumas.ConstDomain{Float64}}}}, fixedparam::NamedTuple{(:θ1, :θ2, :θ3, :θ4, :θ5, :θ6, :θ7, :θ8, :θ11, :Ω, :ΩKin, :σ1, :σ2), Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Matrix{Float64}, Float64, Float64, Float64}}, checkidentification::Bool, diffeq_options::NamedTuple{(), Tuple{}}, init_vrandeffsorth::Vector{Vector{Float64}}, verbose::Bool, optim_state::Nothing)
    @ Pumas /Users/runner/work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Zl9W3/src/estimation/likelihoods.jl:4305
 [12] __fit(m::PumasModel{(θ1 = 1, θ2 = 1, θ3 = 1, θ4 = 1, θ5 = 1, θ6 = 1, θ7 = 1, θ8 = 1, θ11 = 1, Ω = 3, ΩKin = 1, σ1 = 1, σ2 = 1), 3, (:CENT, :PERI, :TRIP), (:HTM, :BMI, :BSA, :FFM, :FMAT, :BIO, :CRCLunits, :CRCLEFF, :TVCL, :CL, :TVV, :V, :Q, :V2, :TVKIN, :KIN, :KINt, :KOUT), ParamSet{NamedTuple{(:θ1, :θ2, :θ3, :θ4, :θ5, :θ6, :θ7, :θ8, :θ11, :Ω, :ΩKin, :σ1, :σ2), Tuple{RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Float64, Float64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Float64, Float64}, RealDomain{Float64, Int64, Float64}, RealDomain{Int64, Float64, Float64}, RealDomain{Int64, Float64, Float64}, PSDDomain{PDMats.PDMat{Float64, Matrix{Float64}}}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Int64, Float64}}}}, var"#705#715", Pumas.TimeDispatcher{var"#706#716", var"#708#718"}, Nothing, var"#710#720", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, SciMLBase.NullParameters, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.ODEFunctionClosure{var"#711#721", var"#712#722"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, var"#713#723", var"#714#724", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, population::Vector{Subject{NamedTuple{(:log_dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:WT, :AGE, :HT, :CRCL, :SEX), Tuple{Float64, Float64, Float64, Float64, Int64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, param::NamedTuple{(:θ1, :θ2, :θ3, :θ4, :θ5, :θ6, :θ7, :θ8, :θ11, :Ω, :ΩKin, :σ1, :σ2), Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Matrix{Float64}, Float64, Float64, Float64}}, approx::FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}, ensemblealg::EnsembleThreads, optimize_fn::Pumas.DefaultOptimizeFN{Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, constantcoef::NTuple{5, Symbol}, omegas::Tuple{}, checkidentification::Bool, diffeq_options::NamedTuple{(), Tuple{}}, init_vrandeffsorth::Vector{Vector{Float64}}, verbose::Bool, optim_state::Nothing)
    @ Pumas /Users/runner/work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Zl9W3/src/estimation/likelihoods.jl:3992
 [13] __fit (repeats 2 times)
    @ /Users/runner/work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Zl9W3/src/estimation/likelihoods.jl:4045 [inlined]
 [14] fit(m::PumasModel{(θ1 = 1, θ2 = 1, θ3 = 1, θ4 = 1, θ5 = 1, θ6 = 1, θ7 = 1, θ8 = 1, θ11 = 1, Ω = 3, ΩKin = 1, σ1 = 1, σ2 = 1), 3, (:CENT, :PERI, :TRIP), (:HTM, :BMI, :BSA, :FFM, :FMAT, :BIO, :CRCLunits, :CRCLEFF, :TVCL, :CL, :TVV, :V, :Q, :V2, :TVKIN, :KIN, :KINt, :KOUT), ParamSet{NamedTuple{(:θ1, :θ2, :θ3, :θ4, :θ5, :θ6, :θ7, :θ8, :θ11, :Ω, :ΩKin, :σ1, :σ2), Tuple{RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Float64, Float64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Float64, Float64}, RealDomain{Float64, Int64, Float64}, RealDomain{Int64, Float64, Float64}, RealDomain{Int64, Float64, Float64}, PSDDomain{PDMats.PDMat{Float64, Matrix{Float64}}}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, Int64, Float64}, RealDomain{Int64, Int64, Float64}}}}, var"#705#715", Pumas.TimeDispatcher{var"#706#716", var"#708#718"}, Nothing, var"#710#720", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, SciMLBase.NullParameters, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.ODEFunctionClosure{var"#711#721", var"#712#722"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, var"#713#723", var"#714#724", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, _population::Vector{Subject{NamedTuple{(:log_dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:WT, :AGE, :HT, :CRCL, :SEX), Tuple{Float64, Float64, Float64, Float64, Int64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Symbol}}, Vector{Float64}}}, param::NamedTuple{(:θ1, :θ2, :θ3, :θ4, :θ5, :θ6, :θ7, :θ8, :θ11, :Ω, :ΩKin, :σ1, :σ2), Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Matrix{Float64}, Float64, Float64, Float64}}, approx::FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}; optim_alg::Nothing, optim_options::Nothing, optimize_fn::Nothing, constantcoef::NamedTuple{(:θ6, :θ7, :θ8, :σ1, :σ2), NTuple{5, Float64}}, omegas::Tuple{}, ensemblealg::EnsembleThreads, checkidentification::Bool, diffeq_options::NamedTuple{(), Tuple{}}, init_randeffs::Nothing, init_vrandeffsorth::Nothing, verbose::Bool)
    @ Pumas /Users/runner/work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Zl9W3/src/estimation/likelihoods.jl:3948
 [15] top-level scope
    @ ~/Desktop/codespaces/pmx/brooks-model.jl:181

andreasnoack · November 4, 2024, 1:11pm

You are hitting an issue with how we show an error message which hides the problem with your model. The error that isn’t shown correctly is trying to tell you that your Ω[2,2] element is not identified. When setting

Ω in PSDDomain(2)

you declare Ω to be a 2x2 positive definite matrix. However, you then set

        ηCL ~ Normal(0, Ω[1])
        ηV ~ Normal(0,Ω[2])

When indexing Ω you extract the (1,1) and (2,1) element of Ω. This is called linear indexing in Julia. However, you don’t use the (2,2) element at all. The fit function detects that and tries to tell you but fails for a separate reason. We will fix that in a future release. The fix is either to keep the correlation between the two random effects and then replace

ηCL ~ Normal(0, Ω[1])
ηV ~ Normal(0,Ω[2])

and

CL = (TVCL/BIO) * exp(ηCL)
...
V = (TVV/BIO) * exp(ηV)

with

ηCLV ~ MvNormal(Ω)

and

CL = (TVCL/BIO) * exp(ηCLV[1])
...
V = (TVV/BIO) * exp(ηCLV[2])

Alternative, you can drop the correlation and replace

Ω in PSDDomain(2)

and

ηCL ~ Normal(0, Ω[1])
ηV ~ Normal(0,Ω[2])

with

ωCL in RealDomain(lower=0.0)
ωV in RealDomain(lower=0.0)

and

ηCL ~ Normal(0, ωCL)
ηV ~ Normal(0, ωV)

Topic		Replies	Views
Error: Gradient is exactly zero in theta Estimation	2	143	March 14, 2024
Turnover PK model simulation Simulation	8	343	March 14, 2022
Non linear kinetics flow dependant Concepts	11	1292	April 20, 2020
Target Mediated Drug Disposition (TMDD)	4	899	April 18, 2020
Error while using infer(res) Estimation	8	1067	March 12, 2020

L2 Correlations

Here is my model file

Error:

Related topics