Help with Fitting

yes - those are the number that pjoenix came up with

thanks. I will take a look and get back to you

It’s always advised to put a lower bound on the variability parameters. Say,

σ_prop ∈ RealDomain(init=0.1, lower=0.001)

Now you have a proportional error model, so conc can still be very small, but even so, this should help because conc can never be negative. The optimizer doesn’t know that one or the other parameter is not meaningful below 0 so you have to tell it. :slight_smile:

1 Like

The immediate error here is because of this issue https://github.com/JuliaStats/Distributions.jl/issues/1014 which will hopefully get fixed soon. However, there might be underlying issues here as well since issue above is triggered when the variance of the error term becomes zero.

Thank you everyone, by adding some constraints - I have been able to run optimization with NaivePooled and FOCEI, by removing subject 1 whose profile was closest to zero.

However - predict, infer and vpc do not work on the model results.

vpc(res1,10) |> plot

returns:

Cannot convert Tuple{Float64,Float64} to series data for plotting

dataframe(predict(res)) # res is teh result of FOCEI optimization  works, but

dataframe(predict(res1)) # res is the result of NaivePooled optimization 

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

Finally,

neither

infer(res)
infer(res1)

works

Calculating: variance-covariance matrix

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

This seems to be related to the issue Andreas commented about having covariance <=0.
This is confusing to me as the model did fit in Phoenix (I am insisting on daily basis to get access to the phoenix code to share with you)

I’ll be happy to share all the code and results if helpful.

A

regarding vpc - Not getting graph for VPC

regarding infer(res) -

  1. did you get decent parameter estimates for your model?
  2. infer is to compute the standard errors and confidence intervals, so I am not that surprised it is unable to compute for just two subjects? (as you removed one?)

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!

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?

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,

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?

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.

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

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 :slight_smile:

A

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

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

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.

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.

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.

3 Likes

Thanks so much -
This clarifies a lot!

I really appreciate it!

A

1 Like

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?