Estimation failing: parameter not identified

EmReg · March 29, 2022, 12:26pm

Hello,

First of all, many thanks to Pumas team for this package as well as the environment provided to work with it.

I have started using Pumas recently to deal with NLME models. I followed some tutorials and I am now trying to make modifications to meet my needs.
Briefly, I have tried to implement a simple model for the fluorescence of individuals (e.g. cells) and to keep it simple fluorescence level will increase or decrease depending on the state of an activation function (see the @dynamicspart of the code).

So far, I have been able to simulate and estimate parameters (logκ,logγ,σ²κ,σ²γ,σ²f) of the model when giving/fixing an initial state for the fluorescence (f=0 in the @init part and commenting all f0 related lines) and it worked great.
However, I will be confronted to data where the starting value of the fluorescence is unknown, so I am now trying to take that into account in my model definition.
In the code below I am able to simulate data with, as far as I understand, random values for f0. The issue arrises when I am using the fit function.
First I get some warnings of dt set as NaN and then when the minimization starts an error (see below code). One of the parameter, logκ, is not identified.
Unfortunately, I do not really know what I can do to change that, it seems to me that I should do something to prevent the dt to be set to NaN but I do not know how (nor if) I can do that.
As already told, this estimation problem worked when considering only parameters κ and γ and it also works when estimating κ, γ and f0 when there is no activation function (commenting it in the @dynamics macro). Also, in all cases I was able to simulate data with the simobs function.

I am thus a bit puzzled and was wondering if they were things to check, first to know whether my problem is theoretically identifiable, and to make the estimation work in Pumas (I have tried other algorithms from Optim). Am I doing something wrong in the model macro?

Also, I might have to re-write this as an @emmodel later, is there more chances for this kind of estimation to work then?

Thanks for the help/advise you can give!

Working example:

using Pumas
using Random
using PumasUtilities
using CairoMakie


## Defining Pumas model with macro

fluo_model = @model begin

@param begin
logκ ∈ RealDomain()
logγ ∈ RealDomain()
logf0 ∈ RealDomain() #
σ²κ ∈ RealDomain(lower=0.0001)
σ²γ ∈ RealDomain(lower=0.0001)
σ²f0 ∈ RealDomain(lower=0.0) #
σ²f ∈ RealDomain(lower=0.0001, upper=20.0)
end

@random begin
ηκ ~ Normal(0.0,sqrt(σ²κ))
ηγ ~ Normal(0.0,sqrt(σ²γ))
ηf0 ~ Normal(0.0,sqrt(σ²f0)) #
end

@pre begin
log_κ = logκ + ηκ
log_γ = logγ + ηγ
log_f0 = logf0 + ηf0 #
end

@init begin
f = exp(log_f0) # f = 0
end

@dynamics begin
f' = - exp(log_γ)*f + exp(log_κ)*1/(1+exp(-(mod(t,40)-30)))
end

@derived begin
Fluo = @. Normal(f,sqrt(σ²f))
end

end


params = (logκ=log(0.31),logγ=log(0.21),logf0=0.01,σ²κ=0.1^2,σ²γ=0.1^2,σ²f=0.1^2,σ²f0=0.1^2)

#params = (logκ=log(0.31),logγ=log(0.21),σ²κ=0.1^2,σ²γ=0.1^2,σ²f=0.1^2)

## For simulations
reg = DosageRegimen(0.01,cmt=1,time=0)
# subject = Subject(id=1,events=reg)
pop = Population(map(i -> Subject(id=i,events=reg),1:5))
obstime = [i for i=1:300];
sim = simobs(fluo_model,pop,params,obstimes=obstime)
sim_plot(fluo_model,sim,observations=[:Fluo])

### Estimates from simulated data
result = fit(fluo_model,Subject.(sim),params, Pumas.FOCEI())

┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq /builds/PumasAI/PumasSystemImages-jl/.julia/packages/OrdinaryDiffEq/UcJTX/src/solve.jl:510
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/LK7N9/src/integrator_interface.jl:325

(The warning messages appear several time)

Iter     Function value   Gradient norm 
     0              Inf     0.000000e+00
 * time: 2.002716064453125e-5
ERROR: LoadError: gradient is exactly zero in logκ. This indicates that logκ isn't identified.
Stacktrace:
  [1] _unindentified_exception(#unused#::Float64, i::Int64, j::Int64, d::Int64, k::Symbol)
    @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/estimation/likelihoods.jl:2242
  [2] _check_zero_gradient(vparam::Vector{Float64}, g::Vector{Float64}, fixedtrf::TransformVariables.TransformTuple{NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), Tuple{TransformVariables.Identity, TransformVariables.Identity, TransformVariables.Identity, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ScaledShiftedLogistic{Float64}}}})
    @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/estimation/likelihoods.jl:2233
  [3] (::Pumas.var"#442#447"{Pumas.FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), Tuple{TransformVariables.Identity, TransformVariables.Identity, TransformVariables.Identity, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ScaledShiftedLogistic{Float64}}}}})(state::Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}})
    @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/estimation/likelihoods.jl:2531
  [4] #426
    @ /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/estimation/likelihoods.jl:1980 [inlined]
  [5] update!(tr::Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}}, iteration::Int64, f_x::Float64, grnorm::Float64, dt::Dict{Any, Any}, store_trace::Bool, show_trace::Bool, show_every::Int64, callback::Pumas.var"#426#428"{Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, Pumas.var"#442#447"{Pumas.FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), Tuple{TransformVariables.Identity, TransformVariables.Identity, TransformVariables.Identity, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ScaledShiftedLogistic{Float64}}}}}}, trace_simplex::Bool)
    @ Optim /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Optim/uwNqi/src/utilities/update.jl:23
  [6] update!
    @ /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Optim/uwNqi/src/utilities/update.jl:11 [inlined]
  [7] trace!(tr::Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}}, d::NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, state::Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}}, iteration::Int64, method::Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, options::Optim.Options{Float64, Pumas.var"#426#428"{Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, Pumas.var"#442#447"{Pumas.FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), Tuple{TransformVariables.Identity, TransformVariables.Identity, TransformVariables.Identity, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ScaledShiftedLogistic{Float64}}}}}}}, curr_time::Float64)
    @ Optim /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Optim/uwNqi/src/multivariate/solvers/first_order/bfgs.jl:172
  [8] optimize(d::NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, method::Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, options::Optim.Options{Float64, Pumas.var"#426#428"{Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, Pumas.var"#442#447"{Pumas.FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), Tuple{TransformVariables.Identity, TransformVariables.Identity, TransformVariables.Identity, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ScaledShiftedLogistic{Float64}}}}}}}, state::Optim.BFGSState{Vector{Float64}, Matrix{Float64}, Float64, Vector{Float64}})
    @ Optim /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Optim/uwNqi/src/multivariate/optimize/optimize.jl:53
  [9] optimize
    @ /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Optim/uwNqi/src/multivariate/optimize/optimize.jl:35 [inlined]
 [10] (::Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}})(cost::NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, p::Vector{Float64}, callback::Pumas.var"#442#447"{Pumas.FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), Tuple{TransformVariables.Identity, TransformVariables.Identity, TransformVariables.Identity, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ScaledShiftedLogistic{Float64}}}}})
    @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/estimation/likelihoods.jl:1973
 [11] _fit(m::PumasModel{ParamSet{NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), Tuple{RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, Float64, Float64}}}}, var"#39#55", var"#40#56", var"#42#58", var"#44#60", ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, ODEFunction{false, ModelingToolkit.ODEFunctionClosure{var"#45#61", var"#46#62"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, var"#47#63", var"#51#67", ModelingToolkit.ODESystem}, population::Vector{Subject{NamedTuple{(:Fluo,), Tuple{Vector{Float64}}}, Pumas.ConstantCovar{NamedTuple{(), Tuple{}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Int64}}}, param::NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), NTuple{7, Float64}}, approx::Pumas.FOCE, ensemblealg::EnsembleThreads, optimize_fn::Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, fixedparamset::ParamSet{NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), Tuple{RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, Float64, Float64}}}}, fixedparam::NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), NTuple{7, Float64}}, checkidentification::Bool, diffeq_options::NamedTuple{(), Tuple{}})
    @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/estimation/likelihoods.jl:2548
 [12] fit(m::PumasModel{ParamSet{NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), Tuple{RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, Float64, Float64}}}}, var"#39#55", var"#40#56", var"#42#58", var"#44#60", ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, ODEFunction{false, ModelingToolkit.ODEFunctionClosure{var"#45#61", var"#46#62"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, var"#47#63", var"#51#67", ModelingToolkit.ODESystem}, population::Vector{Subject{NamedTuple{(:Fluo,), Tuple{Vector{Float64}}}, Pumas.ConstantCovar{NamedTuple{(), Tuple{}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Int64}}}, param::NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), NTuple{7, Float64}}, approx::Pumas.FOCE; optimize_fn::Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, constantcoef::NamedTuple{(), Tuple{}}, omegas::Tuple{}, ensemblealg::EnsembleThreads, checkidentification::Bool, diffeq_options::NamedTuple{(), Tuple{}})
    @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/estimation/likelihoods.jl:2409
 [13] fit(m::PumasModel{ParamSet{NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), Tuple{RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, Float64, Float64}}}}, var"#39#55", var"#40#56", var"#42#58", var"#44#60", ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, ODEFunction{false, ModelingToolkit.ODEFunctionClosure{var"#45#61", var"#46#62"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, var"#47#63", var"#51#67", ModelingToolkit.ODESystem}, population::Vector{Subject{NamedTuple{(:Fluo,), Tuple{Vector{Float64}}}, Pumas.ConstantCovar{NamedTuple{(), Tuple{}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Int64}}}, param::NamedTuple{(:logκ, :logγ, :logf0, :σ²κ, :σ²γ, :σ²f0, :σ²f), NTuple{7, Float64}}, approx::Pumas.FOCE)
    @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/estimation/likelihoods.jl:2392
 [14] top-level scope
    @ ~/data/code/testing nlme/Test_define_model.jl:72
in expression starting at /home/jrun/data/code/testing nlme/Test_define_model.jl:72

patrick · March 29, 2022, 2:58pm

Dear,

Could you try to add

@options begin
    checklinear=false
end

at the top? I’m not 100% sure I understand the issue. Can you also tell me how you’re accessing Pumas?

I get

FittedPumasModel

Successful minimization:                     false

Likelihood approximation:               Pumas.FOCE
Log-likelihood value:                    1301.9996
Number of subjects:                              5
Number of parameters:         Fixed      Optimized
                                  0              7
Observation records:         Active        Missing
    Fluo:                      1500              0
    Total:                     1500              0

---------------------
          Estimate
---------------------
logκ      -1.1416
logγ      -1.7244
logf0      0.034739
σ²κ        0.013475
σ²γ        0.0071028
σ²f0       1.4477e-6
σ²f        0.010048
---------------------

EmReg · March 29, 2022, 3:22pm

Thank you Patrick for your answer!

I added the @options block before the @param one in the model but it does not seem to solve my problem, I get the same error (and warnings).

I am using Pumas app on JuliaHub, I didn’t find (nor received) a link to download it. Could the “set-up” have an impact?
Also sorry if the issue is not clear, let me know what I could clarify.
Did you also have the same error without the @options or was it not reproducible?

patrick · March 29, 2022, 7:25pm

Could you tell me what URL you use to access pumas. Is it http://juliahub.com ? Are you an academic user?

EmReg · March 29, 2022, 10:15pm

Yes, I am an academic user and I use it from juliahub.com (ui/Applications)

korsbo · March 30, 2022, 8:53am

Hi,

I don’t think you’ve done anything wrong at all. Your code works fine on our machines but we’re running a later version of Pumas than you are. I would think that the problem lies in how mod(t, 40) is handled and that machinery has recently undergone a pretty big update. I think the best real fix for your problem is for us to simply expedite the upgrade of the Pumas version that’s on juliahub.com.

In the meantime though, I think I see a little hack if you want to get things running.

In your example, it seems to me that the dose is not actually doing anything, or at least that it’s too small for any real effect. By hijacking the dosing mechanism and adding a new dynamical variable with a constant derivative, we could replace mod.

## Defining Pumas model with macro

fluo_model = @model begin

@param begin
logκ ∈ RealDomain()
logγ ∈ RealDomain()
logf0 ∈ RealDomain() #
σ²κ ∈ RealDomain(lower=0.0001)
σ²γ ∈ RealDomain(lower=0.0001)
σ²f0 ∈ RealDomain(lower=0.0) #
σ²f ∈ RealDomain(lower=0.0001, upper=20.0)
end

@random begin
ηκ ~ Normal(0.0,sqrt(σ²κ))
ηγ ~ Normal(0.0,sqrt(σ²γ))
ηf0 ~ Normal(0.0,sqrt(σ²f0)) #
end

@pre begin
log_κ = logκ + ηκ
log_γ = logγ + ηγ
log_f0 = logf0 + ηf0 #
end

@init begin
f = exp(log_f0) # f = 0
  ramp = 40
end

@dynamics begin
ramp' = -1.
f' = - exp(log_γ)*f + exp(log_κ)*1/(1+exp(-((40-ramp)-30)))
end

@derived begin
Fluo = @. Normal(f,sqrt(σ²f))
end

end


params = (logκ=log(0.31),logγ=log(0.21),logf0=0.01,σ²κ=0.1^2,σ²γ=0.1^2,σ²f=0.1^2,σ²f0=0.1^2)

#params = (logκ=log(0.31),logγ=log(0.21),σ²κ=0.1^2,σ²γ=0.1^2,σ²f=0.1^2)

## For simulations
reg = DosageRegimen(40,cmt=1,time=40, ii=40, addl=10)
# subject = Subject(id=1,events=reg)
pop = Population(map(i -> Subject(id=i,events=reg),1:5))
obstime = [i for i=1:300];
sim = simobs(fluo_model,pop,params,obstimes=obstime)
sim_plot(fluo_model,sim,observations=[:Fluo])


### Estimates from simulated data
result = fit(fluo_model,Subject.(sim),params, Pumas.FOCEI())

Here, ramp ends up with the same temporal dynamics of ramp(t) = 40 - mod(t, 40) so we can replace mod.

So, if you’re in a hurry or if you’re just itching to see if your model works then you could try this. If you’re more patient than that, you could wait for the updated Pumas to come around and your original code should work.

Best regards,
Niklas

EmReg · March 30, 2022, 9:40am

Hello korsbo and thank you for your answer!

I tried what you suggested and it worked! As I might not be able to do it like that with my data I also tried something based on your suggestion. Simply replacing the variable t by ramp (hence having ramp'= 1 and in @init ramp = 0) seems to also resolve my problem. So I guess the issue might be about how t variable itself is handled in the version I am using (mod seems to do fine)?

There is still a difference between applying your answer and mine, in my case the minimization was not successful, so I will work on that now.

Thank you again for your time!

korsbo · March 30, 2022, 10:20am

hmm, that’s a clever thing to try. I’m trying it myself and there is a distinct difference in the fit which is now struggling. I suspect that the core issue is whether the ODE solver is able to deal with the discontinuous derivatives that mod is causing. I’m not sure why it would work with t (on my machine) but not with ramp' = 1. Now that I think about it though, I wonder whether the handling of these discontinuities might have been your original problem.

Could you try your original code but with the addition of a keyword argument to simobs and fit where you force the ODE solver to exactly hit the discontinuous points?

result = fit(fluo_model,Subject.(sim),params, Pumas.FOCEI(); diffeq_options=(; tstops=0:40:300))

could you try this with your original explicit use of t (because your ramp is subject to roundoff errors that might minutely displace the time of the discontinuity)?

EmReg · March 31, 2022, 12:07pm

I tried it, but got the exact same warnings and errors than what I had in the first place.

When I was doing the fit with f0 fixed, and using the mod function already (and t variable), the estimation was working, I am not entirely sure then it is solely a discontinuity issue (at least on the version I am running).

EmReg · March 31, 2022, 12:36pm

I am also trying alternatives to the mod function (mainly because I will need to have activation at non regular points in time), I am now trying to use isless function.

f' = - exp(log_γ)*f + exp(log_κ)*1*(isless(t,38.0)*isless(-t,-30.0)+isless(t,68.0)*isless(-t,-60.0)+isless(t,98.0)*isless(-t,-90.0)+isless(t,128.0)*isless(-t,-120.0)+isless(t,158.0)*isless(-t,-150.0)+isless(t,188.0)*isless(-t,-180.0)+isless(t,218.0)*isless(-t,-210.0)+isless(t,249.91)*isless(-t,-241.91)+isless(t,298.44)*isless(-t,-290.44))

Using this monstruosity in the @dynamics block seems to also work without causing any issue during the fit.
I also tried to define an activ function outside of the @model returning the activation function but then I get the error:

ERROR: LoadError: LoadError: UndefVarError: activ not defined
Stacktrace:
  [1] getproperty
    @ ./Base.jl:26 [inlined]
  [2] convert_rhs_to_Expression(ex::Expr, bvars::Expr, dvars::Vector{SymbolicUtils.Term}, params::Vector{SymbolicUtils.Sym}, t::SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}})
    @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/dsl/model_macro.jl:840
  [3] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
  [5] getindex
    @ ./broadcast.jl:575 [inlined]
  [6] copyto_nonleaf!(dest::Vector{SymbolicUtils.Term{Real, Nothing}}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Pumas.convert_rhs_to_Expression), Tuple{Base.Broadcast.Extruded{Vector{Any}, Tuple{Bool}, Tuple{Int64}}, Tuple{Expr}, Tuple{Vector{SymbolicUtils.Term}}, Tuple{Vector{SymbolicUtils.Sym}}, Tuple{SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}}}}, iter::Base.OneTo{Int64}, state::Int64, count::Int64)
    @ Base.Broadcast ./broadcast.jl:1078
  [7] copy
    @ ./broadcast.jl:930 [inlined]
  [8] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(Pumas.convert_rhs_to_Expression), Tuple{Vector{Any}, Tuple{Expr}, Tuple{Vector{SymbolicUtils.Term}}, Tuple{Vector{SymbolicUtils.Sym}}, Tuple{SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}}}})
    @ Base.Broadcast ./broadcast.jl:883
  [9] convert_rhs_to_Expression(ex::Expr, bvars::Expr, dvars::Vector{SymbolicUtils.Term}, params::Vector{SymbolicUtils.Sym}, t::SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}})
    @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/dsl/model_macro.jl:841
 [10] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
 [11] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
 [12] getindex
    @ ./broadcast.jl:575 [inlined]
 [13] copyto_nonleaf!(dest::Vector{SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing}}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Pumas.convert_rhs_to_Expression), Tuple{Base.Broadcast.Extruded{Vector{Any}, Tuple{Bool}, Tuple{Int64}}, Tuple{Expr}, Tuple{Vector{SymbolicUtils.Term}}, Tuple{Vector{SymbolicUtils.Sym}}, Tuple{SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}}}}, iter::Base.OneTo{Int64}, state::Int64, count::Int64)
    @ Base.Broadcast ./broadcast.jl:1078
 [14] copy
    @ ./broadcast.jl:930 [inlined]
 [15] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(Pumas.convert_rhs_to_Expression), Tuple{Vector{Any}, Tuple{Expr}, Tuple{Vector{SymbolicUtils.Term}}, Tuple{Vector{SymbolicUtils.Sym}}, Tuple{SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}}}})
    @ Base.Broadcast ./broadcast.jl:883
 [16] convert_rhs_to_Expression(ex::Expr, bvars::Expr, dvars::Vector{SymbolicUtils.Term}, params::Vector{SymbolicUtils.Sym}, t::SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}})
    @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/dsl/model_macro.jl:841
 [17] build_mtk_system(pre::OrderedCollections.OrderedSet{Symbol}, odevars::OrderedCollections.OrderedSet{Symbol}, callvars::OrderedCollections.OrderedSet{Symbol}, bvars::Expr, eqs::Expr)
    @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/dsl/model_macro.jl:697
 [18] dynamics_obj(odeexpr::Expr, preexpr::Expr, pre::OrderedCollections.OrderedSet{Symbol}, odevars::OrderedCollections.OrderedSet{Symbol}, callvars::OrderedCollections.OrderedSet{Symbol}, bvars::Expr, eqs::Expr, options::NamedTuple{(:inplace, :subject_t0, :checklinear), Tuple{Bool, Bool, Bool}})
    @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/dsl/model_macro.jl:718
 [19] var"@model"(__source__::LineNumberNode, __module__::Module, expr::Any)
    @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/XHa5R/src/dsl/model_macro.jl:1312
in expression starting at /home/jrun/data/code/clean_test.jl:12

Topic		Replies	Views
Parameter Estimation of DelayDiffEq (and other Non-ODE) Problems Estimation	5	215	October 25, 2024
Bayes estimation functions not working Estimation	4	786	December 24, 2019
Estimation with stochastic differential equations Estimation	6	662	August 4, 2020
Error with fitting using simulated data Estimation	2	271	March 15, 2022
Pumas Model fitting: Error Estimation	5	130	May 10, 2024

Estimation failing: parameter not identified

Related Topics