Help with Fitting

alewolf · November 19, 2019, 8:16pm

Makes sense!

Here are the results for NaivePooled

FittedPumasModel

Successful minimization: true

Likelihood approximation: Pumas.NaivePooled
Objective function value: -8.5046991
Total number of observation records: 10
Number of active observation records: 10
Number of subjects: 2

       Estimate

θ₁ 0.01
θ₂ 0.21883
θ₃ 0.049793
θ₄ 0.040644
θ₅ 0.052324
θ₆ 1.2279
σ_prop 0.01

Results for FOCEI are weirs specially theta_6
FittedPumasModel

Successful minimization: true

Likelihood approximation: Pumas.FOCEI
Objective function value: 18.396454
Total number of observation records: 10
Number of active observation records: 10
Number of subjects: 2

       Estimate

θ₁ 0.020736
θ₂ 0.01
θ₃ 0.01
θ₄ 0.01
θ₅ 0.01
θ₆ 3.6819e17
σ_prop 0.36486

My constraints where for all theta as well as sigma_prop lower =0.01

NaivePooled with the 3 samples (res2)
gives

FittedPumasModel

Successful minimization: true

Likelihood approximation: Pumas.NaivePooled
Objective function value: -1.0753737
Total number of observation records: 15
Number of active observation records: 13
Number of subjects: 3

       Estimate

θ₁ 0.01
θ₂ 6.4777
θ₃ 0.010113
θ₄ 0.10051
θ₅ 0.019268
θ₆ 0.010004
σ_prop 0.072009

#but neither 
predict(res2)
#nor
infer(res2)
#nor 
vpc(res2,10) |> plot

work
julia predict(res2) returns

MethodError: no method matching _predict(::PumasModel{ParamSet{NamedTuple{(:θ, :σ_prop),Tuple{VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}},RealDomain{Float64,TransformVariables.Infinity{true},Float64}}}},getfield(Main, Symbol("##353#360")),getfield(Main, Symbol("##354#361")),getfield(Main, Symbol("##355#362")),ODEProblem{Nothing,Tuple{Nothing,Nothing},false,Nothing,ODEFunction{false,getfield(Main, Symbol("##356#363")),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}},DiffEqBase.StandardODEProblem},getfield(Main, Symbol("##358#365")),getfield(Main, Symbol("##359#366"))}, ::Subject{NamedTuple{(:dv,),Tuple{Array{Union{Missing, Float64},1}}},Nothing,Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1},Array{Float64,1}}, ::NamedTuple{(:θ, :σ_prop),Tuple{Array{Float64,1},Float64}}, ::Pumas.NaivePooled, ::Array{Any,1})
Closest candidates are:
_predict(::Any, ::Any, ::Any, !Matched::Pumas.FO, ::Any) at C:\Users\awolf-yadlin.juliapro\JuliaPro_v1.2.0-1\packages\Pumas\6uorK\src\estimation\diagnostics.jl:542
_predict(::Any, ::Any, ::Any, !Matched::Union{Pumas.FOCE, Laplace}, ::Any) at C:\Users\awolf-yadlin.juliapro\JuliaPro_v1.2.0-1\packages\Pumas\6uorK\src\estimation\diagnostics.jl:549
_predict(::Any, ::Any, ::Any, !Matched::Union{Pumas.FOCEI, Pumas.LaplaceI}, ::Any) at C:\Users\awolf-yadlin.juliapro\JuliaPro_v1.2.0-1\packages\Pumas\6uorK\src\estimation\diagnostics.jl:556

julia infer(res2) returns
Calculating: variance-covariance matrix
MethodError: no method matching _orth_empirical_bayes!(::Array{Float64,1}, ::PumasModel{ParamSet{NamedTuple{(:θ, :σ_prop),Tuple{VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}},RealDomain{Float64,TransformVariables.Infinity{true},Float64}}}},getfield(Main, Symbol("##353#360")),getfield(Main, Symbol("##354#361")),getfield(Main, Symbol("##355#362")),ODEProblem{Nothing,Tuple{Nothing,Nothing},false,Nothing,ODEFunction{false,getfield(Main, Symbol("##356#363")),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}},DiffEqBase.StandardODEProblem},getfield(Main, Symbol("##358#365")),getfield(Main, Symbol("##359#366"))}, ::Subject{NamedTuple{(:dv,),Tuple{Array{Union{Missing, Float64},1}}},Nothing,Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1},Array{Float64,1}}, ::NamedTuple{(:θ, :σ_prop),Tuple{Array{Float64,1},Float64}}, ::Pumas.NaivePooled)
Closest candidates are:
_orth_empirical_bayes!(::AbstractArray{T,1} where T, ::PumasModel, ::Subject, ::NamedTuple, !Matched::Union{Pumas.FO, Pumas.FOI, Pumas.HCubeQuad}, !Matched::Any…; kwargs…) at C:\Users\awolf-yadlin.juliapro\JuliaPro_v1.2.0-1\packages\Pumas\6uorK\src\estimation\likelihoods.jl:228
_orth_empirical_bayes!(::AbstractArray{T,1} where T, ::PumasModel, ::Subject, ::NamedTuple, !Matched::Union{Pumas.FOCE, Pumas.FOCEI, Laplace, Pumas.LaplaceI}, !Matched::Any…; reltol, fdtype, fdrelstep, kwargs…) at C:\Users\awolf-yadlin.juliapro\JuliaPro_v1.2.0-1\packages\Pumas\6uorK\src\estimation\likelihoods.jl:245

and julia vpc(res2,10) |> plot returns
MethodError: Cannot convert an object of type Missing to an object of type Float64
Closest candidates are:
convert(::Type{T<:Real}, !Matched::Quantity) where T<:Real at C:\Users\awolf-yadlin.juliapro\JuliaPro_v1.2.0-1\packages\Unitful\ytsW0\src\conversion.jl:141
convert(::Type{T<:Real}, !Matched::Level) where T<:Real at C:\Users\awolf-yadlin.juliapro\JuliaPro_v1.2.0-1\packages\Unitful\ytsW0\src\logarithm.jl:22
convert(::Type{T<:Real}, !Matched::Gain) where T<:Real at C:\Users\awolf-yadlin.juliapro\JuliaPro_v1.2.0-1\packages\Unitful\ytsW0\src\logarithm.jl:62

which I assume ahve to do with the two missing values in the PKdata, but I am not sure how to deal with them in this context (ie how to remove them).
Moreover, given that infer doesn’t work I assume vpc will not work as it will lack the CIs for the shaded plots even if we didnt have the missing issue correct?

Thanks-

A

PS - I know my posts are long, but I am trying to provide as much info as possible!

vijay · November 19, 2019, 9:13pm

thanks, it helps with the information. did you remove all the etas from your model to run NaivePooled()? and what were the starting parameters you used as the initiat estimates for your fit?

alewolf · November 19, 2019, 9:22pm

Here are the two models:

FOCEI:

PK = @model begin

    @param begin
      θ ∈ VectorDomain(6,lower=[0.01, 0.01, 0.01, 0.01, 0.01, 0.01], init=[0.011, 1000, 1000, 1000, 0.011, 0.011]) 
#       Ω ∈ PSDDomain(5)
       σ_prop ∈ RealDomain(lower=0.01, init=0.04)
    end

    @random begin
      η ~ MvNormal(Matrix{Float64}(0.04I, 5, 5))
#       η ~ MvNormal(Ω)
    end

    @pre begin
        CL_adc        = θ[1]
        CLD_adc       = θ[2]*exp(η[1])
        V1_adc        = θ[3]*exp(η[2])
        V2_adc        = θ[4]*exp(η[3])
        Vmax          = θ[5]*exp(η[2])
        Km            = θ[6]*exp(η[3])
    end


    @dynamics begin
        A1'  = -CL_adc*C1 - CLD_adc*(C1-C2) - Vmax*C1/(Km+C1)
        A2'  =  CLD_adc*(C1-C2)
    end
    
    @vars begin
        C1   :=  A1/V1_adc
        C2   :=  A2/V2_adc
    end
        

    @derived begin
#         a1     = A1
        Conc_Mean     = C1
        dv ~ @. Normal(Conc_Mean, sqrt(Conc_Mean^2*σ_prop)+eps())
    end
    
end

NaivePooled:


PK1 = @model begin

    @param begin
      θ ∈ VectorDomain(6,lower=[0.01, 0.01, 0.01, 0.01, 0.01, 0.01], init=[0.011, 1000, 1000, 1000, 0.011, 0.011]) 
#       Ω ∈ PSDDomain(5)
       σ_prop ∈ RealDomain(lower=0.01, init=0.04)
    end

    @random begin
#       η ~ MvNormal(Matrix{Float64}(0.04I, 5, 5))
#       η ~ MvNormal(Ω)
    end

    @pre begin
        CL_adc        = θ[1]
        CLD_adc       = θ[2]
        V1_adc        = θ[3]
        V2_adc        = θ[4]
        Vmax          = θ[5]
        Km            = θ[6]
    end


    @dynamics begin
        A1'  = -CL_adc*C1 - CLD_adc*(C1-C2) - Vmax*C1/(Km+C1)
        A2'  =  CLD_adc*(C1-C2)
    end
    
    @vars begin
        C1   :=  A1/V1_adc
        C2   :=  A2/V2_adc
    end
        

    @derived begin
#         a1     = A1
        Conc_Mean     = C1
        dv ~ @. Normal(Conc_Mean, sqrt(Conc_Mean^2*σ_prop)+eps())
    end
    
end

Also FOCEI model with 3 ran under this conditions, but as you see sigma_prop is crazy

FittedPumasModel

Successful minimization: true

Likelihood approximation: Pumas.FOCEI
Objective function value: 345.0777
Total number of observation records: 15
Number of active observation records: 13
Number of subjects: 3

       Estimate

θ₁ 4647.9
θ₂ 0.01
θ₃ 0.01
θ₄ 0.01
θ₅ 0.010978
θ₆ 0.011388
σ_prop 4.5943e39

infer(res3)

Calculating: variance-covariance matrix

PosDefException: matrix is not positive definite; Cholesky factorization failed.

julia preds3 = DataFrame(predict(res3)) works fine

julia vpc(res3,10) fails

MethodError: Cannot convert an object of type Missing to an object of type Float64
Closest candidates are:
convert(::Type{T<:Real}, !Matched::Quantity) where T<:Real at C:\Users\awolf-yadlin.juliapro\JuliaPro_v1.2.0-1\packages\Unitful\ytsW0\src\conversion.jl:141
convert(::Type{T<:Real}, !Matched::Level) where T<:Real at C:\Users\awolf-yadlin.juliapro\JuliaPro_v1.2.0-1\packages\Unitful\ytsW0\src\logarithm.jl:22
convert(::Type{T<:Real}, !Matched::Gain) where T<:Real at C:\Users\awolf-yadlin.juliapro\JuliaPro_v1.2.0-1\packages\Unitful\ytsW0\src\logarithm.jl:62

As you see only prdicts works, but only for FOCEI not NaivePooled. infer or vpc don’t for either.
I’ll try to put upper limit in sigma_prop, that may help
Take care,

vijay · November 19, 2019, 9:30pm

Just to note again - VPC is NOT a supported feature in the version that you have. It will be updated shortly.
The infer issue is something that we know about and should be more stable in the next release. but the main question I have for you is, do these parameter values make sense to you?

alewolf · November 19, 2019, 10:34pm

Yes and no -
The sigma_prop is crazy. The other values make more sense as what we are having is a super quick clearing system. The Vmax and Km were added to account for the quick drop on concentration, while keeping the CL parameter relatively small, as I am not constraining it, the quick disappearance of my analyte can be explained either by this parameters or the Phoenix ones, that I shared with you before.
Overall, I am happy with the results.

pkofod · November 19, 2019, 11:10pm

Is there any chance you could try to simulate data in Pumas and estimate the model in your other software? I’m still curious if something is not matching between the two systems (in terms of syntax or whatever).

alewolf · November 19, 2019, 11:22pm

Do you mean generate data in pumas with some group of parameters and see if I can recover the parameters on Phoenix?

I can do the rims for sure, I’ll then need to ask for the Phoenix runs, it may take a couple of days/weeks for my go-worker to get back to me, but I’ll start asking today

A

pkofod · November 20, 2019, 5:58am

Exactly! Something is obviously strange here, and I just want to verify that the model as written actually does what it’s supposed to

alewolf · November 20, 2019, 5:34pm

Hi -
I finally have the Phoenix model available as well as a larger dataset, and a word file with the textual model. As attachments here are limited to figures. Is there a way for me to share this files as requested with you?

Thanks

alewolf · November 20, 2019, 7:05pm

quick, unrelated, question - when you are doing a dose - jhow do you specify to which of all the variables in your equation system the dose corresponds?
Is it just dosed to the top variable in the dynamic system? I don’t see where to specify that in the regimen.

pkofod · November 20, 2019, 9:39pm

You can chose the cmt https://docs.pumas.ai/dev/basics/doses_subjects_populations/#Dosage-Regimen-Terminology-1 . Edit: but yes, it defaults to the first state/compartment in your dynamics block.

pkofod · November 20, 2019, 9:59pm

Okay, maybe that wasn’t as helpful as I could have been!

Consider this slightly odd model.

using Pumas

mdsl1 = @model begin
    @param begin
        θ ∈ VectorDomain(1, init=[0.5])
        Ω ∈ RealDomain(lower=0.001, init=0.5)
    end

    @random begin
        η ~ Normal(0.0, Ω)
    end

    @pre begin
        CL = θ[1] * exp(η)
        V  = 1.0
    end

    @vars begin
        conc1 = Central1 / V
        conc2 = Central2 / V
    end

    @dynamics begin
      Central1' = -CL/V*Central1
      Central2' = -CL/V*Central2
    end
    @derived begin
        dv1 ~ @. Normal(conc1,0.01)
        dv2 ~ @. Normal(conc2,0.01)
    end
end

There are two central compartments that act exactly the same, but they are disconnected. So how can we dose into Central1?

dose = DosageRegimen(43, cmt=1, time=0, ii=1, addl=0, rate=0)
pop   = Population(map(i -> Subject(id=i, evs=dose),1:4))
_simobs = simobs(mdsl1, pop, init_param(mdsl1); obstime=range(0,10;step=0.1))

and you can plot.

using Plots
plot(_simobs)

And intro Central2?

dose = DosageRegimen(43, cmt=2, time=0, ii=1, addl=0, rate=0)
pop   = Population(map(i -> Subject(id=i, evs=dose),1:4))
_simobs = simobs(mdsl1, pop, init_param(mdsl1); obstime=range(0,10;step=0.1))
plot(_simobs)

What about in Central1 at time 0 and then in Central2 at time 2?

dose1 = DosageRegimen(43, cmt=1, time=0, ii=1, addl=0, rate=0)
dose2 = DosageRegimen(43, cmt=2, time=2, ii=1, addl=0, rate=0)
dose = DosageRegimen(dose1, dose2)
pop   = Population(map(i -> Subject(id=i, evs=dose),1:4))
_simobs = simobs(mdsl1, pop, init_param(mdsl1); obstime=range(0,10;step=0.1))
plot(_simobs)

But really, I prefer to use the symbols, say :Central2. But you can combine them:

dose1 = DosageRegimen(43, cmt=1, time=0, ii=1, addl=0, rate=0)
dose2 = DosageRegimen(43, cmt=:Central2, time=2, ii=1, addl=0, rate=0)
dose = DosageRegimen(dose1, dose2)
pop   = Population(map(i -> Subject(id=i, evs=dose),1:4))
_simobs = simobs(mdsl1, pop, init_param(mdsl1); obstime=range(0,10;step=0.1))
plot(_simobs)

So if th cmt keyword is an integer it’s the order of specification in dynamics and if it’s a symbol, it’s the name of the state / compartment.

alewolf · November 20, 2019, 10:02pm

Thanks so much -
This clarifies a lot!

I really appreciate it!

A

alewolf · November 21, 2019, 6:55pm

I sent the Phoenix files and a non-agregated dataset.
The parameter estimation should be similar to the data that I already shared as NaivePooled was also used to estimate.

Thanks for all the help?

Topic		Replies	Views
How to estimate parameters for data maximum-likelihood	10	1090	November 12, 2019
Constraining Random effects? Estimation	2	493	December 10, 2019
Fixing parameters using ConstDomain Estimation	3	776	November 14, 2019
Esimated parameters is equal to initial value & inspect() function didn't work Estimation	12	424	November 18, 2021
Defining different error models for data from different sources Estimation	2	183	February 24, 2023

Help with Fitting

θ₁ 0.01 θ₂ 0.21883 θ₃ 0.049793 θ₄ 0.040644 θ₅ 0.052324 θ₆ 1.2279 σ_prop 0.01

θ₁ 0.020736 θ₂ 0.01 θ₃ 0.01 θ₄ 0.01 θ₅ 0.01 θ₆ 3.6819e17 σ_prop 0.36486

θ₁ 0.01 θ₂ 6.4777 θ₃ 0.010113 θ₄ 0.10051 θ₅ 0.019268 θ₆ 0.010004 σ_prop 0.072009

θ₁ 4647.9 θ₂ 0.01 θ₃ 0.01 θ₄ 0.01 θ₅ 0.010978 θ₆ 0.011388 σ_prop 4.5943e39

Related Topics

θ₁ 0.01
θ₂ 0.21883
θ₃ 0.049793
θ₄ 0.040644
θ₅ 0.052324
θ₆ 1.2279
σ_prop 0.01

θ₁ 0.020736
θ₂ 0.01
θ₃ 0.01
θ₄ 0.01
θ₅ 0.01
θ₆ 3.6819e17
σ_prop 0.36486

θ₁ 0.01
θ₂ 6.4777
θ₃ 0.010113
θ₄ 0.10051
θ₅ 0.019268
θ₆ 0.010004
σ_prop 0.072009

θ₁ 4647.9
θ₂ 0.01
θ₃ 0.01
θ₄ 0.01
θ₅ 0.010978
θ₆ 0.011388
σ_prop 4.5943e39