Help with Fitting

alewolf · October 17, 2019, 9:27pm

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.

ChrisRackauckas · October 18, 2019, 2:24pm

You need to specify an error distribution on the measurements

alewolf · October 21, 2019, 3:34pm

Thanks I’ll try it-

Ale

alewolf · October 21, 2019, 4:01pm

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

alewolf · October 21, 2019, 6:10pm

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

sai_matcha · November 11, 2019, 6:18pm

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

alewolf · November 11, 2019, 6:30pm

I have not gotten answer yet.

Sorry-

sai_matcha · November 11, 2019, 6:56pm

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…

vijay · November 11, 2019, 8:46pm

@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

alewolf · November 11, 2019, 10:35pm

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

alewolf · November 11, 2019, 10:46pm

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

Thanks,

A

vijay · November 11, 2019, 11:46pm

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?

alewolf · November 12, 2019, 12:17am

Hi Vijay,

Thank you!

I’ll try this tonight

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

vijay · November 12, 2019, 3:24am

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

alewolf · November 12, 2019, 4:02am

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

Ale

alewolf · November 18, 2019, 5:45pm

sorry-
Still waiting from co-worker for Phoenix model

alewolf · November 18, 2019, 8:31pm

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

pkofod · November 18, 2019, 8:40pm

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,

alewolf · November 18, 2019, 8:47pm

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

pkofod · November 18, 2019, 9:00pm

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,

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

Help with Fitting

Related Topics