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?

1 Like

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

sorry-
Still waiting from co-worker for Phoenix model

1 Like

Quick question -
And this might be the base of my confusion.
Why are you fitting the new model to dv and not to Conc_Mean?
Conc_Mean are the meassured concentrations, shouldn’t they be fit to the expected concentrations (Conc_Mean?) and not to dv that I assume are the deviations.

Also I ran your model with julia param =init(PK) - I assume that what yo did as without defining parameters the simulation crashed and got very different plots from what you got. (I copy pasted your code to make sure everything was exactly the same, but I got this very different plot?

sim

I am not sure whats going on. I have Pumas 5.0

Thanks for all the help!

A

Remember, we only have “conc_mean” because we are simulating the data. dv does not stand for deviations, it stands for dependent variable, so it is the observed concentration (including noise).

Best,

thanks-

So, just to be clear this term
julia dv ~ @. Normal(Conc_Mean, sqrt(Conc_Mean^2*σ_prop)+eps())
represents the distribution of my concentrations at any given time?

I assume Conc_mean is just the mean, but I want to fit to the distributions then?

What I am confused now is how do I get initial conditions from the data frame, which I assume is causing the differences between Vijay’s plot and mine.

Thanks

A
Thanks

Yeah, so if you forget that we simulate data, we would normally have some model function f that has a distribution. The model function is often just the mean of the concentration and the distribution is then normal, lognormal, …

You want to fit based on dv because this it what has the distributional assumption encoded.

Can you be specific about “initial conditions”? What is it exactly you’re looking for?

Best,