Help with Fitting

Hi I am trying to fit a very simple model and failing - it might be not enough data, bur I’d love to know if I am doing something wrong.

Model

PK = @model begin

    @param begin
      θ ∈ VectorDomain(6,lower=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], init=[1.0, 1.0, 1.0, 1.0, 1.0, 1.0]) 
#       Ω ∈ PSDDomain(5)
#       σ_prop ∈ RealDomain(init=0.1)
    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(η[4])
        Km            = θ[6]*exp(η[5])
    end

    @init begin
    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
    end
    
end

Params:

param = init_param(PK)

Data
PKDataMean
|Group|Time|Conc_Mean|
||Any|Float64|Any|
|1|0.1 mg/kg|0.0034722|2.12803|
|2|0.03 mg/kg|0.0034722|0.293033|
|3|0.3 mg/kg|0.0034722|5.51947|
|4|0.1 mg/kg|0.25|1.07887|
|5|0.03 mg/kg|0.25|0.151067|
|6|0.3 mg/kg|0.25|3.18253|
|7|0.1 mg/kg|1.0|0.7315|
|8|0.03 mg/kg|1.0|0.0420333|
|9|0.3 mg/kg|1.0|2.52407|
|10|0.1 mg/kg|2.0|0.4226|
|11|0.03 mg/kg|2.0|missing|
|12|0.3 mg/kg|2.0|2.0864|
|13|0.1 mg/kg|2.99|0.3219|
|14|0.03 mg/kg|2.99|missing|
|15|0.3 mg/kg|2.99|1.3471|

(I also made a table where I remove Group 0.03 mg/kg, in case the problem was caused by missing data)

Filtering if necesary

PKDataMean=filter(row -> row.Group ∈ ["0.1 mg/kg", "0.3 mg/kg"], PKDataMean)

read_as_Pumas and fit

PKdata = read_pumas(PKDataMean, id=:Group, time=:Time, evid=1, dvs=[:Conc_Mean])
res = fit(PK,PKdata,param,Pumas.FOCEI(), solver=Rodas5())

ERROR

TypeError: in typeassert, expected ForwardDiff.Dual{ForwardDiff.Tag{getfield(Pumas, Symbol("##145#146")){Float64,Base.Iterators.Pairs{Symbol,Rodas5{0,true,DefaultLinSolve,DataType},Tuple{Symbol},NamedTuple{(:solver,),Tuple{Rodas5{0,true,DefaultLinSolve,DataType}}}},PumasModel{ParamSet{NamedTuple{(:θ,),Tuple{VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}}}}},getfield(Main, Symbol("##167#174")),getfield(Main, Symbol("##168#175")),getfield(Main, Symbol("##169#176")),ODEProblem{Nothing,Tuple{Nothing,Nothing},false,Nothing,ODEFunction{false,getfield(Main, Symbol("##170#177")),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},getfield(Main, Symbol("##172#179")),getfield(Main, Symbol("##173#180"))},Subject{NamedTuple{(:Conc_Mean,),Tuple{Array{Union{Missing, Float64},1}}},Nothing,Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1},Array{Float64,1}},NamedTuple{(:θ,),Tuple{Array{Float64,1}}},Pumas.FOCEI,Tuple{}},Float64},Float64,5}, got Float64

if I skip the filtering the error changes to
DimensionMismatch(“vectors must have same length”) -> I assume this is consequence of the missing values.

Thanks and hope this is not to wordy.

You need to specify an error distribution on the measurements

Thanks I’ll try it-

Ale

So maybe, I am not doing the right changes -
I modified the code above, to add the error distribution (I think)
with same model as above


@param begin
      θ ∈ VectorDomain(6,lower=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], init=[1.0, 1.0, 1.0, 1.0, 1.0, 1.0]) 
      Ω ∈ PSDDomain(5)
      σ_prop ∈ RealDomain(init=0.04)
    end

    @derived begin
        Conc_Mean     = C1
        
         # I assume this is how you set up error distribution?
        dv ~ @. Normal(Conc_Mean, sqrt(Conc_Mean^2*σ_prop))
    end


PKdata = read_pumas(PKDataMean, id=:Group, time=:Time, evid=1, dvs=[:Conc_Mean])
param = init_param(PK)
res = fit(PK,PKdata,param,Pumas.FOCEI(), solver=Rodas5())

It reads -
DimensionMismatch(“vectors must have same length”) -> not sure which vectors it is refering to.
Thanks,

A

Also a side question -

When I build

PKdata = read_pumas(PKDataMean, id=:Group, time=:Time, evid=1, dvs=[:Conc_Mean])

How can I see the values in PKdata?

thanks

Hi alewolf , did you get the answer. even im getting same error

I have not gotten answer yet.

Sorry-

thank you… one thing found…
i have copy pasted all codes from the working example still i got same error while using the fix code…

@alewolf you are trying to fit a PK model here without a dose into the system. is the grouping column your dose? Pumas is trying to see what goes into the system and currently in read_pumas you specify no amt

Yes the grouping is my dose.
Thanks/

It’s a single injection.

Do I need to make an amnt Column with the group value at t=0 and and 0 everywhere else?

Thanks!

Ale

Also just curious the conc column does not represent amount?. Those are the actual measured PK values

Thanks,

A

hi - I am not sure I understand the experimental design. A dose of 0.1mg/kg was given at time zero and concentrations were measured at regular time intervals. Is that right? If so,
I created a modified datasheet that you can use.

| id | Group      | time      | dv        | amt  | evid |
|----|------------|-----------|-----------|------|------|
| 1  | 0.03 mg/kg | 0         | 0         | 0.03 | 1    |
| 1  | 0.03 mg/kg | 0.0034722 | 0.293033  | 0    | 0    |
| 1  | 0.03 mg/kg | 0.25      | 0.151067  | 0    | 0    |
| 1  | 0.03 mg/kg | 1         | 0.0420333 | 0    | 0    |
| 1  | 0.03 mg/kg | 2         | missing   | 0    | 0    |
| 1  | 0.03 mg/kg | 2.99      | missing   | 0    | 0    |
| 2  | 0.1 mg/kg  | 0         | 0         | 0.1  | 1    |
| 2  | 0.1 mg/kg  | 0.0034722 | 2.12803   | 0    | 0    |
| 2  | 0.1 mg/kg  | 0.25      | 1.07887   | 0    | 0    |
| 2  | 0.1 mg/kg  | 1         | 0.7315    | 0    | 0    |
| 2  | 0.1 mg/kg  | 2         | 0.4226    | 0    | 0    |
| 2  | 0.1 mg/kg  | 2.99      | 0.3219    | 0    | 0    |
| 3  | 0.3 mg/kg  | 0         | 0         | 0.3  | 1    |
| 3  | 0.3 mg/kg  | 0.0034722 | 5.51947   | 0    | 0    |
| 3  | 0.3 mg/kg  | 0.25      | 3.18253   | 0    | 0    |
| 3  | 0.3 mg/kg  | 1         | 2.52407   | 0    | 0    |
| 3  | 0.3 mg/kg  | 2         | 2.0864    | 0    | 0    |
| 3  | 0.3 mg/kg  | 2.99      | 1.3471    | 0    | 0    |
using Pumas
using CSV
using Plots

You can read this dataset from a spreadsheet using

df  = CSV.read("testdata.csv", missingstring="missing")

PKdata = read_pumas(df, id=:id, time=:time, dvs=[:dv])

I slightly modified your model to add the error framework that was discussed earlier. Just for completeness, I am pasting the entire model here.

PK = @model begin

    @param begin
      θ ∈ VectorDomain(6,lower=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0], init=[1.0, 1.0, 1.0, 1.0, 1.0, 1.0]) 
#       Ω ∈ PSDDomain(5)
       σ_prop ∈ RealDomain(init=0.1)
    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(η[4])
        Km            = θ[6]*exp(η[5])
    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

My first check always is to make sure the model is runnable by quickly running a simulation.

sims = simobs(PK,PKdata,param) 
plot(sims)

I know the initial estimates are all 1 and probably don’t make sense, but just tried it out nevertheless

Now, I tried fitting the data

res = fit(PK,PKdata,param,Pumas.FOCEI())

The model does not fit due to some ODE issues, but I am not sure if it mechanstically is correct. Did you try that model elsewhere. If so, we can take a deeper look into why it is not working?

Hi Vijay,

Thank you!

I’ll try this tonight :slight_smile:

The model was run in Phoenix and it did fit, my initial parameters were zero. Other than that. That’s the same ODE model I had.

Thanks!

Ale

share the textual model of phoenix and we can replicate it here.

I’ll get it from my co-worker tomorrow.
Thanks-

Ale