# Help with Fitting

**URL:** <https://discourse.pumas.ai/t/help-with-fitting/50>\
**Category:** Estimation\
**Created:** [October 17, 2019, 9:27pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50 "2019-10-17T21:27:10Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![vijay](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/vijay/32/96_2.png) [@vijay](https://discourse.pumas.ai/u/vijay)\
**Post date:** [November 18, 2019, 9:34pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/21 "2019-11-18T21:34:36Z")

</div>

> [@alewolf](#):
>
> Why are you fitting the new model to dv and not to Conc\_Mean?

This is because, I changed the dataset headers so that the name of the concentration column is now `dv`. I posted the dataset in the same message.

> [@alewolf](#):
>
> 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?

Could you provide the full code that you ran. Or actually, the `simobs` that your ran please? From what I see in your plot, it looks like you did not pass the dataset in correctly.

> [@alewolf](#):
>
> 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.

You don’t get the initial conditions from the dataframe. They come from this block below when you use `init_param(PK)` NOT `init(PK)`

```auto

    @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

```

what would be a good set of initial parameter values given the knowledge of your model?

---

<div class="post-metadata">

**Author:** ![alewolf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/alewolf/32/24_2.png) [@alewolf](https://discourse.pumas.ai/u/alewolf)\
**Post date:** [November 18, 2019, 10:08pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/22 "2019-11-18T22:08:44Z")

</div>

This is the exact code I ran:

This is the csv file i used - (from the one you made)

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

```auto

df = CSV.read("Data//FixedData.csv", missingstring="missing", copycols=true)

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

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

```

thats it -

A not so bad set of initial parameters according to literature (1st) and Phoenix fit (second) is:

```auto

param = (θ = [
              2.4 # CL_adc Clearance ADC 2.4 ml/d/kg 
              17 # CLD_adc Clearance distribution ADC 17 ml/d/kg 
              29 # V1_adc Centarl Compartment Distribution Volume 29 ml/kg
              30 # V2_adc Peripheral Compartment Distribution Volume 30 ml/kg
              51 # Vmax Ask Gabby 51 ug/d/kg
              0.93 # Km Ask Gabby 0.93 ug/ml ])

or

param=(θ = [
              5.64445E-05 # CL_adc Clearance ADC 2.4 ml/d/kg 
              143.927 # CLD_adc Clearance distribution ADC 17 ml/d/kg 
              51.0108 # V1_adc Centarl Compartment Distribution Volume 29 ml/kg
              59.4517 # V2_adc Peripheral Compartment Distribution Volume 30 ml/kg
              56.5194 # Vmax Ask Gabby 51 ug/d/kg
              0.93 # Km Ask Gabby 0.93 ug/ml ])

```

Thanks so much! I appreciate all the help!

---

<div class="post-metadata">

**Author:** ![vijay](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/vijay/32/96_2.png) [@vijay](https://discourse.pumas.ai/u/vijay)\
**Post date:** [November 18, 2019, 10:29pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/23 "2019-11-18T22:29:22Z")

</div>

Here are the plots with three different sets of initial parameters

1. `init_param`  

2. set 1 that you provided  

3. set 2 that you provided  

Some quick things to check at your end before moving on to fitting.

1. I don’t see you use `read_pumas` in your code, did you use it to set your `PKdata`?
2. Can you share the literature model where you took the parameters from, it will be nice to ensure that the model is written up correctly
3. What version of JuliaPro are you on and what version of Pumas?

I’ve tried to fit your data with the model you provided, but I am not getting a convergence. What estimation method was used to fit in phoenix? DId you really have 6 eats for 3 groups estimated in phoenix?

---

<div class="post-metadata">

**Author:** ![alewolf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/alewolf/32/24_2.png) [@alewolf](https://discourse.pumas.ai/u/alewolf)\
**Post date:** [November 18, 2019, 10:37pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/24 "2019-11-18T22:37:59Z")

</div>

I am working on getting the Phoenix model from my co-worker that wrote it.  
I have never used Phoenix, so I am hard pressed to answer those details, until she finally ahres it with me.  
My goal is to move people in my company from Phoenix to Pumas - where data can be shared more readily and i can build apps for other people to fit and/or simulate data.

I did do

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

```

before setting up the model

This is my current status

Status `C:\Users\awolf-yadlin\.juliapro\JuliaPro_v1.2.0-1\environments\v1.2\Project.toml`  
[c52e3926] Atom v0.11.2  
[a134a8b2] BlackBoxOptim v0.5.0  
[336ed68f] CSV v0.5.13  
[535c2557] ClinicalTrialUtilities v0.1.14  
[a93c6f00] DataFrames v0.19.4  
[31c24e10] Distributions v0.21.3  
[7073ff75] IJulia v1.20.0  
[e5e0dc1b] Juno v0.7.2  
[91a5bcdd] Plots v0.27.0  
[4f2c3c20] Pumas v0.5.0  
[b7b41870] PumasTutorials v0.0.1  
[f3b207a7] StatsPlots v0.12.0  
[fdbf4ff8] XLSX v0.5.7  
[37e2e46d] LinearAlgebra  
[9a3f8284] Random

As you see julia 1.2 and Pumas 0.5.0 -  
Somehow it seems the simulation is not getting an initial condition for :dv or Conc\_Mean from the table for me.  
My plots just reflect eps()

The literature parameters came from several papers - the model is a simply 2 compartment where the km and vmax are used to simulate an abrupt fall in concentration - I’ll ask my colleague about the papers too, as my part in this project is only to get it to work and estimate in Pumas 🙂

thanks again-

A

---

<div class="post-metadata">

**Author:** ![vijay](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/vijay/32/96_2.png) [@vijay](https://discourse.pumas.ai/u/vijay)\
**Post date:** [November 18, 2019, 10:40pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/25 "2019-11-18T22:40:51Z")

</div>

Can you share the result of  
`sims1 = simobs(PK,PKdata,param1)`

just type sims1 into the console.

Also, can you do `PKdata[1].events` and `PKdata[1].time` and `PKdata[1].observations` in your console and past the output here please

---

<div class="post-metadata">

**Author:** ![alewolf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/alewolf/32/24_2.png) [@alewolf](https://discourse.pumas.ai/u/alewolf)\
**Post date:** [November 18, 2019, 10:43pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/26 "2019-11-18T22:43:26Z")

</div>

so this might explain my zero line

`julia PKdata[1].observations` :(dv = Union{Missing, Float64}[0.0, 0.293033, 0.151067, 0.0420333, missing, missing],)

`julia PKdata[1].observations`: 0-element Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1}  
Shouldn’t there be at least one event occurring?

`PKdata[1].time`

6-element Array{Float64,1}:  
0.0  
0.0034722  
0.25  
1.0  
2.0  
2.99

Thanks

---

<div class="post-metadata">

**Author:** ![vijay](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/vijay/32/96_2.png) [@vijay](https://discourse.pumas.ai/u/vijay)\
**Post date:** [November 18, 2019, 10:50pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/27 "2019-11-18T22:50:07Z")

</div>

> [@alewolf](#):
>
> `julia PKdata[1].observations` : 0-element Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1}  
> Shouldn’t there be at least one event occurring?

yup… somehow you are not reading the events correctly. for me it comes like this…

```auto

julia> PKdata[1].events
1-element Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1}:
 Dose event
  dose amount = 0.03
  dose time = 0.0
  compartment = 1
  instantaneous
  interdose interval = 0.0
  infusion start time = 0.0

```

FWIW, just restart your julia and go through each step…

---

<div class="post-metadata">

**Author:** ![alewolf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/alewolf/32/24_2.png) [@alewolf](https://discourse.pumas.ai/u/alewolf)\
**Post date:** [November 18, 2019, 10:50pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/28 "2019-11-18T22:50:25Z")

</div>

same result on sim1 plot with params 1 ![sim1](https://canada1.discourse-cdn.com/flex030/uploads/pumas/original/1X/c4d496f148f0bbe6537352f8c0a427f98a814a59.png)

the issue I think is ot detecting events form table (at least regarding simulations)

---

<div class="post-metadata">

**Author:** ![alewolf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/alewolf/32/24_2.png) [@alewolf](https://discourse.pumas.ai/u/alewolf)\
**Post date:** [November 18, 2019, 10:51pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/29 "2019-11-18T22:51:06Z")

</div>

I’ll try restarting

A

---

<div class="post-metadata">

**Author:** ![alewolf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/alewolf/32/24_2.png) [@alewolf](https://discourse.pumas.ai/u/alewolf)\
**Post date:** [November 18, 2019, 10:59pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/30 "2019-11-18T22:59:38Z")

</div>

so I discover what happened -  
The csv file column named “evid” is named "evid " -\> I erased the space and now things are working.

Now i need to figure out the fitting :), which is still not working ☹

thanks for so much help and sorry the error was so stupid I couldn’t visually see it, so I didn’t think about it -

Sorry,

A

---

<div class="post-metadata">

**Author:** ![vijay](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/vijay/32/96_2.png) [@vijay](https://discourse.pumas.ai/u/vijay)\
**Post date:** [November 18, 2019, 11:11pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/31 "2019-11-18T23:11:06Z")

</div>

no worries. I am glad you figured it out. These are pretty common errors.  
With regards to the fitting, I will look into it more and get back. Did you put in all those `eta` random effects or did the model you copied to pumas have it. Unusual to have 5 etas for 3 subjects

---

<div class="post-metadata">

**Author:** ![alewolf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/alewolf/32/24_2.png) [@alewolf](https://discourse.pumas.ai/u/alewolf)\
**Post date:** [November 18, 2019, 11:15pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/32 "2019-11-18T23:15:13Z")

</div>

I just added them - as I thought each parameter should have a random effect?  
Is that wrong?  
I am pretty new to PKPD, so I assumed that when we estimate parameters -there should be a random noise to them.  
you mean we would need at least 3 groups/ subjects to have 5 etas?  
Are the etas being calculated by the noise between the 3 groups?

PS- I just runned it, it didn’t crash as badly as before 9using 5 eta still) it gave issues with taking sqrt(negative number) - I assume this could be an issue with starting params?

DomainError with -0.10722130015692335:  
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).

-same issue with 3 etas

A

---

<div class="post-metadata">

**Author:** ![vijay](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/vijay/32/96_2.png) [@vijay](https://discourse.pumas.ai/u/vijay)\
**Post date:** [November 18, 2019, 11:20pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/33 "2019-11-18T23:20:29Z")

</div>

Correct. The etas you added are subject level parameters ( in this case each group is assigned to a subject number). 5 noise parameters to three subjects seems a lot. I would start with testing without any etas using a `Pumas.NaivePooled()` approach, which unfortunately, is not documented yet here [https://docs.pumas.ai/dev/basics/estimation/#Marginal-Likelihood-Approximations-1](https://docs.pumas.ai/dev/basics/estimation/#Marginal-Likelihood-Approximations-1)

will update the doc and try with this and get back to you.

---

<div class="post-metadata">

**Author:** ![vijay](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/vijay/32/96_2.png) [@vijay](https://discourse.pumas.ai/u/vijay)\
**Post date:** [November 18, 2019, 11:21pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/34 "2019-11-18T23:21:39Z")

</div>

yeah, I am getting the same one and am assuming something to do with the error distributions… the second set of parameters you provided are the final values?

---

<div class="post-metadata">

**Author:** ![alewolf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/alewolf/32/24_2.png) [@alewolf](https://discourse.pumas.ai/u/alewolf)\
**Post date:** [November 18, 2019, 11:22pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/35 "2019-11-18T23:22:09Z")

</div>

yes - those are the number that pjoenix came up with

---

<div class="post-metadata">

**Author:** ![vijay](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/vijay/32/96_2.png) [@vijay](https://discourse.pumas.ai/u/vijay)\
**Post date:** [November 18, 2019, 11:22pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/36 "2019-11-18T23:22:35Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![pkofod](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/pkofod/32/59_2.png) [@pkofod](https://discourse.pumas.ai/u/pkofod)\
**Post date:** [November 18, 2019, 11:28pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/37 "2019-11-18T23:28:48Z")

</div>

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

```auto
σ_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. 🙂

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/andreasnoack/32/22_2.png) [@andreasnoack](https://discourse.pumas.ai/u/andreasnoack)\
**Post date:** [November 19, 2019, 12:10pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/38 "2019-11-19T12:10:01Z")

</div>

The immediate error here is because of this issue [https://github.com/JuliaStats/Distributions.jl/issues/1014](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.

---

<div class="post-metadata">

**Author:** ![alewolf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/alewolf/32/24_2.png) [@alewolf](https://discourse.pumas.ai/u/alewolf)\
**Post date:** [November 19, 2019, 7:52pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/39 "2019-11-19T19:52:05Z")

</div>

> [@vijay](#):
>
> Pumas.NaivePooled()

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.

```julia
vpc(res1,10) |> plot

```

returns:

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

```julia
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

```julia
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

---

<div class="post-metadata">

**Author:** ![vijay](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/vijay/32/96_2.png) [@vijay](https://discourse.pumas.ai/u/vijay)\
**Post date:** [November 19, 2019, 7:59pm UTC](https://discourse.pumas.ai/t/help-with-fitting/50/40 "2019-11-19T19:59:43Z")

</div>

regarding `vpc` - [Not getting graph for VPC](https://discourse.pumas.ai/t/not-getting-graph-for-vpc/77/4)

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

[Previous page](https://discourse.pumas.ai/t/help-with-fitting/50.md?page=1)

[Next page](https://discourse.pumas.ai/t/help-with-fitting/50.md?page=3)