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.
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\awolfyadlin.juliapro\JuliaPro_v1.2.01\packages\Pumas\6uorK\src\estimation\diagnostics.jl:542
_predict(::Any, ::Any, ::Any, !Matched::Union{Pumas.FOCE, Laplace}, ::Any) at C:\Users\awolfyadlin.juliapro\JuliaPro_v1.2.01\packages\Pumas\6uorK\src\estimation\diagnostics.jl:549
_predict(::Any, ::Any, ::Any, !Matched::Union{Pumas.FOCEI, Pumas.LaplaceI}, ::Any) at C:\Users\awolfyadlin.juliapro\JuliaPro_v1.2.01\packages\Pumas\6uorK\src\estimation\diagnostics.jl:556
Finally,
neither
infer(res)
infer(res1)
works
Calculating: variancecovariance 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)

 did you get decent parameter estimates for your model?

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\awolfyadlin.juliapro\JuliaPro_v1.2.01\packages\Pumas\6uorK\src\estimation\diagnostics.jl:542
_predict(::Any, ::Any, ::Any, !Matched::Union{Pumas.FOCE, Laplace}, ::Any) at C:\Users\awolfyadlin.juliapro\JuliaPro_v1.2.01\packages\Pumas\6uorK\src\estimation\diagnostics.jl:549
_predict(::Any, ::Any, ::Any, !Matched::Union{Pumas.FOCEI, Pumas.LaplaceI}, ::Any) at C:\Users\awolfyadlin.juliapro\JuliaPro_v1.2.01\packages\Pumas\6uorK\src\estimation\diagnostics.jl:556
julia infer(res2)
returns
Calculating: variancecovariance 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\awolfyadlin.juliapro\JuliaPro_v1.2.01\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\awolfyadlin.juliapro\JuliaPro_v1.2.01\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\awolfyadlin.juliapro\JuliaPro_v1.2.01\packages\Unitful\ytsW0\src\conversion.jl:141
convert(::Type{T<:Real}, !Matched::Level) where T<:Real at C:\Users\awolfyadlin.juliapro\JuliaPro_v1.2.01\packages\Unitful\ytsW0\src\logarithm.jl:22
convert(::Type{T<:Real}, !Matched::Gain) where T<:Real at C:\Users\awolfyadlin.juliapro\JuliaPro_v1.2.01\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*(C1C2)  Vmax*C1/(Km+C1)
A2' = CLD_adc*(C1C2)
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*(C1C2)  Vmax*C1/(Km+C1)
A2' = CLD_adc*(C1C2)
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: variancecovariance 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\awolfyadlin.juliapro\JuliaPro_v1.2.01\packages\Unitful\ytsW0\src\conversion.jl:141
convert(::Type{T<:Real}, !Matched::Level) where T<:Real at C:\Users\awolfyadlin.juliapro\JuliaPro_v1.2.01\packages\Unitful\ytsW0\src\logarithm.jl:22
convert(::Type{T<:Real}, !Matched::Gain) where T<:Real at C:\Users\awolfyadlin.juliapro\JuliaPro_v1.2.01\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 goworker to get back to me, but I’ll start asking today
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
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/#DosageRegimenTerminology1 . 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.
Thanks so much 
This clarifies a lot!
I really appreciate it!
A
I sent the Phoenix files and a nonagregated 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?