Error with SAEM Algorithm

ahmed.salem · September 19, 2022, 1:12pm

Hi,
I wrote the following model to estimate the parameters through FOCEI algorithm


mod = @model begin
        @param   begin
            CLss      ∈ RealDomain(lower = 0.000)
            CLinit      ∈ RealDomain(lower = 0.000)
            ecmo      ∈ RealDomain()
            T12      ∈ RealDomain(lower = 0.000)
            TVVc      ∈ RealDomain(lower  = 0.000) 
            wt_cl     ∈ RealDomain()
            circuit     ∈ RealDomain(lower = 0.000)
            wt_v      ∈ RealDomain()
            Ω         ∈ PDiagDomain(4)
            σ_add     ∈ RealDomain(lower = 0.0000)
            σ_prop    ∈ RealDomain(lower = 0.0000)
        end
        @random begin
            η ~ MvNormal(Ω)
        end
        @covariates WT AGE_D ISMALE SCR GFR AGEYRS TYPE_MODELING IS_BLEEDING IS_CIRCUIT_CHANGE ECMO_DAYS Occassions IS_CIRCUIT_CHANGE_UPDATE IS_BLEEDING_UPDATE
        @pre begin
            CL_SS  = CLss *exp(η[1]) *(WT/50)^wt_cl * circuit^IS_CIRCUIT_CHANGE
            CL_init = CLinit *exp(η[2]) *(WT/50)^wt_cl * circuit^IS_CIRCUIT_CHANGE
            T_12 = T12 * exp(η[3])*(ECMO_DAYS/6)^ecmo   
            Vc     = TVVc * exp(η[4])*(WT/50)^wt_v
            CL     = CL_SS - (CL_SS - CL_init)*exp(-t/T_12) 
        end

        @dynamics  begin
                Central'  = -(CL/Vc) * Central
        end 
      
        @derived begin
            CP = @. Central / Vc
            CONC = @. Normal(abs(CP), sqrt(CP^2*σ_prop+σ_add))
        end
end

I want also to estimate the parameters through SAEM algorithm and compare according to the following code (which I tried to match with the one in the documentation)

df = @chain df begin
     @transform :logwt = log.(:WT/50)
     @transform :logecmo = log.(:ECMO_DAYS/6)
end
saem = @emmodel begin

        @random begin
                CLss   ~ 1  + logwt + IS_CIRCUIT_CHANGE  | LogNormal
                CLinit ~ 1  + logwt + IS_CIRCUIT_CHANGE  | LogNormal
                T12    ~ 1  + logecmo                    | LogNormal
                Vc     ~ 1  + logwt                      | LogNormal
        end
    
        @covariates logwt logecmo WT AGE_D ISMALE SCR GFR AGEYRS TYPE_MODELING IS_BLEEDING IS_CIRCUIT_CHANGE ECMO_DAYS Occassions IS_CIRCUIT_CHANGE_UPDATE IS_BLEEDING_UPDATE
    
        @pre begin
           CL     = CLss - (CLss - CLinit)*exp(-t/T12) 
        end
    
        @dynamics begin
                Central' = -(CL/Vc) * Central
        end
    
        @post begin
            μ = Central/Vc
        end
    
        @error begin
            CONC ~ CombinedNormal(μ)
        end
end

rngv = [MersenneTwister(1941964947i + 1) for i ∈ 1:Threads.nthreads()];
init_saem = init_covariate = (CLss= (1.25, 2/3,0.18), CLinit= (0.53, 2/3,0.18), T12= (1.6, 0.68), Vc = (1.6,1))

@time fit_saem = fit(saem, estimation, init_saem, Pumas.SAEM(), ensemblealg = EnsembleThreads(), rng=rngv)

But I am getting the following error

@time fit_saem = fit(saem, estimation, init_saem, Pumas.SAEM(), ensemblealg = EnsembleThreads(), rng=rngv)
ERROR: MethodError: no method matching distribution_inverse_transform(::Tuple{Float64, Int64}, ::Pumas.Formula{:Vc, (1, :logwt), :LogNormal})
Closest candidates are:
  distribution_inverse_transform(::Tuple{Vararg{T, N}} where {N, T}, ::Pumas.Formula{Y, M, D}) where {Y, M, D} at /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/estimation/saem.jl:939
  distribution_inverse_transform(::StaticArraysCore.SVector, ::Pumas.Formula{Y, M, D}) where {Y, M, D} at /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/estimation/saem.jl:946
Stacktrace:
 [1] __inverse_formula_transform
   @ /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/estimation/saem.jl:950 [inlined]
 [2] _inverse_formula_transform
   @ /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/estimation/saem.jl:962 [inlined]
 [3] _inverse_formula_transform (repeats 3 times)
   @ /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/estimation/saem.jl:956 [inlined]
 [4] inverse_formula_transform(formulas::Tuple{Pumas.Formula{:CLss, (1, :logwt, :IS_CIRCUIT_CHANGE), :LogNormal}, Pumas.Formula{:CLinit, (1, :logwt, :IS_CIRCUIT_CHANGE), :LogNormal}, Pumas.Formula{:T12, (1, :logecmo), :LogNormal}, Pumas.Formula{:Vc, (1, :logwt), :LogNormal}}, nt::NamedTuple{(:CLss, :CLinit, :T12, :Vc), Tuple{Tuple{Float64, Float64, Float64}, Tuple{Float64, Float64, Float64}, Tuple{Float64, Float64}, Tuple{Float64, Int64}}})
   @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/estimation/saem.jl:967
 [5] inverse_formula_transform(m::Pumas.PumasEMModel{(1, 1, 1, 1), 0, (:logwt, :logecmo, :WT, :AGE_D, :ISMALE, :SCR, :GFR, :AGEYRS, :TYPE_MODELING, :IS_BLEEDING, :IS_CIRCUIT_CHANGE, :ECMO_DAYS, :Occassions, :IS_CIRCUIT_CHANGE_UPDATE, :IS_BLEEDING_UPDATE), (:CL, Symbol("##A##"), Symbol("##b##")), (), (:Central,), (:μ,), Nothing, Tuple{Pumas.Formula{:CLss, (1, :logwt, :IS_CIRCUIT_CHANGE), :LogNormal}, Pumas.Formula{:CLinit, (1, :logwt, :IS_CIRCUIT_CHANGE), :LogNormal}, Pumas.Formula{:T12, (1, :logecmo), :LogNormal}, Pumas.Formula{:Vc, (1, :logwt), :LogNormal}}, var"#14#17", var"#13#16", Pumas.LinearODE, var"#15#18", Tuple{Pumas.Error{:CONC, :μ, :CombinedNormal}}, ModelingToolkit.ODESystem}, nt::NamedTuple{(:CLss, :CLinit, :T12, :Vc), Tuple{Tuple{Float64, Float64, Float64}, Tuple{Float64, Float64, Float64}, Tuple{Float64, Float64}, Tuple{Float64, Int64}}})
   @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/estimation/saem.jl:968
 [6] fit(m::Pumas.PumasEMModel{(1, 1, 1, 1), 0, (:logwt, :logecmo, :WT, :AGE_D, :ISMALE, :SCR, :GFR, :AGEYRS, :TYPE_MODELING, :IS_BLEEDING, :IS_CIRCUIT_CHANGE, :ECMO_DAYS, :Occassions, :IS_CIRCUIT_CHANGE_UPDATE, :IS_BLEEDING_UPDATE), (:CL, Symbol("##A##"), Symbol("##b##")), (), (:Central,), (:μ,), Nothing, Tuple{Pumas.Formula{:CLss, (1, :logwt, :IS_CIRCUIT_CHANGE), :LogNormal}, Pumas.Formula{:CLinit, (1, :logwt, :IS_CIRCUIT_CHANGE), :LogNormal}, Pumas.Formula{:T12, (1, :logecmo), :LogNormal}, Pumas.Formula{:Vc, (1, :logwt), :LogNormal}}, var"#14#17", var"#13#16", Pumas.LinearODE, var"#15#18", Tuple{Pumas.Error{:CONC, :μ, :CombinedNormal}}, ModelingToolkit.ODESystem}, pop::Vector{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:logwt, :logecmo, :WT, :AGE_D, :ISMALE, :SCR, :GFR, :AGEYRS, :TYPE_MODELING, :IS_BLEEDING, :IS_CIRCUIT_CHANGE, :ECMO_DAYS, :Occassions, :IS_CIRCUIT_CHANGE_UPDATE, :IS_BLEEDING_UPDATE), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, μnt::NamedTuple{(:CLss, :CLinit, :T12, :Vc), Tuple{Tuple{Float64, Float64, Float64}, Tuple{Float64, Float64, Float64}, Tuple{Float64, Float64}, Tuple{Float64, Int64}}}, saem::Pumas.SAEM; ensemblealg::EnsembleThreads, diffeq_options::NamedTuple{(:alg,), Tuple{CompositeAlgorithm{Tuple{Vern7, Rodas5{0, true, DefaultLinSolve, Val{:forward}, true}}, AutoSwitch{Vern7, Rodas5{0, true, DefaultLinSolve, Val{:forward}, true}, Rational{Int64}, Int64}}}}, rng::Vector{MersenneTwister})
   @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/estimation/saem.jl:5630
 [7] top-level scope
   @ ./timing.jl:220 [inlined]
 [8] top-level scope
   @ ~/data/code/UFH_ECMO_UTAH/UFH_UTAH/ECMO_UTAH_NEW.jl:0

Is there anything wrong with how the model was written ?
thanks

chriselrod · September 19, 2022, 2:57pm

It should work if you replace:

(CLss= (1.25, 2/3,0.18), CLinit= (0.53, 2/3,0.18), T12= (1.6, 0.68), Vc = (1.6,1))

with

(CLss= (1.25, 2/3,0.18), CLinit= (0.53, 2/3,0.18), T12= (1.6, 0.68), Vc = (1.6,1.0))

it seems it currently doesn’t handle the 1 being an integer.

We’ll fix this in an upcoming release.

ahmed.salem · September 20, 2022, 12:11am

@chriselrod Great, it works now! I have a follow up question. If I do not want to add the effect of covariate weight in the random block but rather in the pre block since I want to estimate one value for weight effect on CL and not on CLss and CLinit separately as the following code :

saem = @emmodel begin

        @random begin
                CLss   ~ 1     | LogNormal
                CLinit ~ 1     | LogNormal
                T12    ~ 1  + logecmo | LogNormal
                Vc     ~ 1  + logwt   | LogNormal
        end
    
        @covariates logwt logecmo WT AGE_D ISMALE SCR GFR AGEYRS TYPE_MODELING IS_BLEEDING IS_CIRCUIT_CHANGE ECMO_DAYS Occassions IS_CIRCUIT_CHANGE_UPDATE IS_BLEEDING_UPDATE
    
        @pre begin
           CL     = CLss - (CLss - CLinit)*exp(-t/T12) * (WT/50)^wt
        end
    
        @dynamics begin
                Central' = -(CL/Vc) * Central
        end
    
        @post begin
            μ = Central/Vc
        end
    
        @error begin
            CONC ~ CombinedNormal(μ)
        end
end

rngv = [MersenneTwister(1941964947i + 1) for i ∈ 1:Threads.nthreads()];
init_saem = init_covariate = (CLss= (1.25, 2/3), CLinit= (0.53, 2/3), T12= (1.6, 0.68), Vc = (1.6,1.0), wt = (2/3,))

@time fit_saem = fit(saem, estimation, init_saem, Pumas.SAEM(), ensemblealg = EnsembleThreads(), rng=rngv)

but I am getting an error with the above code. So how to specify the wt exponent in that case and how to include it in the init_covariate NamedTuple ? thanks

vijay · September 20, 2022, 9:28am

Since wt is a fixed effect covariate parameter that you want to estimate (without a random effect), you can include a @param block where you state

@param begin
    wt ~ 1 | LogNormal
end

then you can include wt in the NamedTuple

ahmed.salem · September 20, 2022, 2:47pm

Ok thanks. I tried to use ηshrinkage(fit_saem) & infer(fit_saem, Pumas.Bootstrap(samples=500, ensemblealg = EnsembleThreads())) functions but I am getting the following errors

ηshrinkage(fit_saem)
ERROR: BoundsError
Stacktrace:
 [1] getindex
   @ ./number.jl:98 [inlined]
 [2] #1026
   @ /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/estimation/diagnostics.jl:1154 [inlined]
 [3] macro expansion
   @ ./ntuple.jl:74 [inlined]
 [4] ntuple
   @ ./ntuple.jl:69 [inlined]
 [5] _ηshrinkage(m::Pumas.PumasEMModel{(1, 1, 1, 1, 1), 1, (:logwt, :logecmo, :WT, :AGE_D, :ISMALE, :SCR, :GFR, :AGEYRS, :TYPE_MODELING, :IS_BLEEDING, :IS_CIRCUIT_CHANGE, :ECMO_DAYS, :Occassions, :IS_CIRCUIT_CHANGE_UPDATE, :IS_BLEEDING_UPDATE), (:_CL, :CL, Symbol("##A##"), Symbol("##b##")), (), (:Central,), (:μ,), Nothing, Tuple{Pumas.Formula{:wt, (1,), :LogNormal}, Pumas.Formula{:CLss, (1,), :LogNormal}, Pumas.Formula{:CLinit, (1,), :LogNormal}, Pumas.Formula{:T12, (1, :logecmo), :LogNormal}, Pumas.Formula{:Vc, (1, :logwt), :LogNormal}}, var"#32#35", var"#31#34", Pumas.LinearODE, var"#33#36", Tuple{Pumas.Error{:CONC, :μ, :CombinedNormal}}, ModelingToolkit.ODESystem}, param::NamedTuple{(:μ, :ω, :σ), Tuple{NamedTuple{(:wt, :CLss, :CLinit, :T12, :Vc), Tuple{Float64, Float64, Float64, StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SVector{2, Float64}}}, NTuple{4, Float64}, Tuple{NamedTuple{(:ₐ, :ₚ), Tuple{Float64, Float64}}}}}, η::Array{Float64, 3})
   @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/estimation/diagnostics.jl:1154
 [6] ηshrinkage(fpm::Pumas.FittedPumasEMModel{NamedTuple{(:wt, :CLss, :CLinit, :T12, :Vc), Tuple{Float64, Float64, Float64, StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SVector{2, Float64}}}, NTuple{4, Float64}, Tuple{NamedTuple{(:ₐ, :ₚ), Tuple{Float64, Float64}}}, Pumas.SubjectSetup{Pumas.LinearODE, StaticArraysCore.SVector{1, Float64}, Tuple{Float64, Float64}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}, Subject{NamedTuple{(:CONC,), Tuple{Vector{Float64}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:logwt, :logecmo, :WT, :AGE_D, :ISMALE, :SCR, :GFR, :AGEYRS, :TYPE_MODELING, :IS_BLEEDING, :IS_CIRCUIT_CHANGE, :ECMO_DAYS, :Occassions, :IS_CIRCUIT_CHANGE_UPDATE, :IS_BLEEDING_UPDATE), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, Pumas.SAEM, Array{Float64, 3}, Pumas.PumasEMModel{(1, 1, 1, 1, 1), 1, (:logwt, :logecmo, :WT, :AGE_D, :ISMALE, :SCR, :GFR, :AGEYRS, :TYPE_MODELING, :IS_BLEEDING, :IS_CIRCUIT_CHANGE, :ECMO_DAYS, :Occassions, :IS_CIRCUIT_CHANGE_UPDATE, :IS_BLEEDING_UPDATE), (:_CL, :CL, Symbol("##A##"), Symbol("##b##")), (), (:Central,), (:μ,), Nothing, Tuple{Pumas.Formula{:wt, (1,), :LogNormal}, Pumas.Formula{:CLss, (1,), :LogNormal}, Pumas.Formula{:CLinit, (1,), :LogNormal}, Pumas.Formula{:T12, (1, :logecmo), :LogNormal}, Pumas.Formula{:Vc, (1, :logwt), :LogNormal}}, var"#32#35", var"#31#34", Pumas.LinearODE, var"#33#36", Tuple{Pumas.Error{:CONC, :μ, :CombinedNormal}}, ModelingToolkit.ODESystem}, Vector{Tuple{NTuple{7, Float64}, NTuple{5, Float64}, Tuple{Tuple{Float64, Float64}}, Float64}}, NamedTuple{(:ensemblealg, :diffeq_options), Tuple{EnsembleThreads, NamedTuple{(:alg,), Tuple{CompositeAlgorithm{Tuple{Vern7, Rodas5{0, true, DefaultLinSolve, Val{:forward}, true}}, AutoSwitch{Vern7, Rodas5{0, true, DefaultLinSolve, Val{:forward}, true}, Rational{Int64}, Int64}}}}}}})
   @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/estimation/diagnostics.jl:1114
 [7] top-level scope
   @ ~/data/code/UFH_ECMO_UTAH/UFH_UTAH/ECMO_UTAH_NEW.jl:421

bs_saem = infer(fit_saem, Pumas.Bootstrap(samples=500, ensemblealg = EnsembleThreads()))
ERROR: MethodError: no method matching infer(::Pumas.FittedPumasEMModel{NamedTuple{(:wt, :CLss, :CLinit, :T12, :Vc), Tuple{Float64, Float64, Float64, StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SVector{2, Float64}}}, NTuple{4, Float64}, Tuple{NamedTuple{(:ₐ, :ₚ), Tuple{Float64, Float64}}}, Pumas.SubjectSetup{Pumas.LinearODE, StaticArraysCore.SVector{1, Float64}, Tuple{Float64, Float64}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}, Subject{NamedTuple{(:CONC,), Tuple{Vector{Float64}}}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, NamedTuple{(:logwt, :logecmo, :WT, :AGE_D, :ISMALE, :SCR, :GFR, :AGEYRS, :TYPE_MODELING, :IS_BLEEDING, :IS_CIRCUIT_CHANGE, :ECMO_DAYS, :Occassions, :IS_CIRCUIT_CHANGE_UPDATE, :IS_BLEEDING_UPDATE), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}}}, Symbol}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, Pumas.SAEM, Array{Float64, 3}, Pumas.PumasEMModel{(1, 1, 1, 1, 1), 1, (:logwt, :logecmo, :WT, :AGE_D, :ISMALE, :SCR, :GFR, :AGEYRS, :TYPE_MODELING, :IS_BLEEDING, :IS_CIRCUIT_CHANGE, :ECMO_DAYS, :Occassions, :IS_CIRCUIT_CHANGE_UPDATE, :IS_BLEEDING_UPDATE), (:_CL, :CL, Symbol("##A##"), Symbol("##b##")), (), (:Central,), (:μ,), Nothing, Tuple{Pumas.Formula{:wt, (1,), :LogNormal}, Pumas.Formula{:CLss, (1,), :LogNormal}, Pumas.Formula{:CLinit, (1,), :LogNormal}, Pumas.Formula{:T12, (1, :logecmo), :LogNormal}, Pumas.Formula{:Vc, (1, :logwt), :LogNormal}}, var"#32#35", var"#31#34", Pumas.LinearODE, var"#33#36", Tuple{Pumas.Error{:CONC, :μ, :CombinedNormal}}, ModelingToolkit.ODESystem}, Vector{Tuple{NTuple{7, Float64}, NTuple{5, Float64}, Tuple{Tuple{Float64, Float64}}, Float64}}, NamedTuple{(:ensemblealg, :diffeq_options), Tuple{EnsembleThreads, NamedTuple{(:alg,), Tuple{CompositeAlgorithm{Tuple{Vern7, Rodas5{0, true, DefaultLinSolve, Val{:forward}, true}}, AutoSwitch{Vern7, Rodas5{0, true, DefaultLinSolve, Val{:forward}, true}, Rational{Int64}, Int64}}}}}}}, ::Pumas.Bootstrap)
Closest candidates are:
  infer(::Pumas.FittedPumasModel, ::Pumas.Bootstrap; level) at /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/estimation/inference.jl:107
  infer(::Pumas.AbstractFittedPumasModel; level, rethrow_error, sandwich_estimator) at /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/estimation/inference.jl:59
Stacktrace:
 [1] top-level scope
   @ ~/data/code/UFH_ECMO_UTAH/UFH_UTAH/ECMO_UTAH_NEW.jl:427

Is the syntax is different from when using FOCEI ?

chriselrod · September 21, 2022, 1:37pm

Hi, sorry for the poor error messages. They are currently not supported for SAEM fits, but we intend to start supporting both methods in the Pumas 2.4 release.

Topic		Replies	Views
SAEM in Pumas: ERROR: TaskFailedException Estimation	1	95	May 10, 2024
SAEM algorithm parameters saem	26	681	April 20, 2022
Covariates in SAEM	2	239	May 4, 2022
Defining different error models for data from different sources Estimation	2	245	February 24, 2023
M3 method with SAEM Estimation	2	132	February 9, 2024

Error with SAEM Algorithm

Related Topics