Defining different error models for data from different sources

Estimation
1

Hello,

I have datasets collected from multiple clinical sites and I would like to define different error model based on the clinical site. I wrote the code as following

mod = @model begin
    @param   begin
        TVCL      ∈ RealDomain(lower= 0.0001, upper = 6)
        tvKTV     ∈ RealDomain(lower = 0.001, upper = 2)
        TVVc      ∈ RealDomain(lower = 0.001, upper = 20) 
        ecmo      ∈ RealDomain(lower = -3,upper = 3)
        cl      ∈ RealDomain(lower = -3,upper = 3)
        v      ∈ RealDomain(lower = -3,upper = 3)
        circuit   ∈ RealDomain(lower = 0.001, upper = 5) 
        Ω         ∈ PDiagDomain(3)
        σ_prop_1    ∈ RealDomain(lower = 0.00001, upper = 1)
        σ_prop_2    ∈ RealDomain(lower = 0.00001, upper = 1)
        σ_add ∈ RealDomain(lower = 1)
    end
    @random begin
        η ~ MvNormal(Ω)
    end
    @covariates WT ECMO_DAYS IS_CIRCUIT_CHANGE site
    @pre begin
        _site = site
        KTV    = tvKTV*exp(η[1])*(ECMO_DAYS/6)^ecmo
        CL     = TVCL * exp(η[2])*(1-exp(-KTV*t))*(WT/10)^cl*circuit^IS_CIRCUIT_CHANGE
        Vc     = TVVc * exp(η[3])*(WT/10)^v
  
    end
    @dynamics  begin 
        Central' = -(CL/Vc)*Central
    end
    @derived begin
        CP = @. Central / Vc
        CONC = if _site ==2  
                  @. Normal(CP, sqrt(CP^2*σ_prop_1 + σ_add)) 
          else 
                 @. Normal(CP, sqrt(CP^2*σ_prop_2))
          end
    end
end

When I run the fit function I get the following error

julia> modfit = fit(mod, estimation, param, Pumas.FOCEI()) 
[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter     Function value   Gradient norm 
     0     1.897660e+04     1.198511e+05
 * time: 3.0040740966796875e-5
ERROR: gradient is exactly zero in σ_prop_1. This indicates that σ_prop_1 isn't identified.
Stacktrace:
  [1] _unidentified_exception(#unused#::Float64, i::Int64, j::Int64, d::Int64, k::Symbol)
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:3484
  [2] _check_zero_gradient(vparam::Vector{Float64}, g::Vector{Float64}, fixedtrf::TransformVariables.TransformTuple{NamedTuple{(:TVCL, :tvKTV, :TVVc, :ecmo, :cl, :v, :circuit, :Ω, :σ_prop_1, :σ_prop_2, :σ_add), Tuple{TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Float64}, Pumas.PDiagTransform, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.Identity, TransformVariables.ShiftedExp{true, Float64}}}})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:3475
  [3] (::Pumas.var"#614#616"{FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:TVCL, :tvKTV, :TVVc, :ecmo, :cl, :v, :circuit, :Ω, :σ_prop_1, :σ_prop_2, :σ_add), Tuple{TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Float64}, Pumas.PDiagTransform, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.Identity, TransformVariables.ShiftedExp{true, Float64}}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}})(state::Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:4047
  [4] #585
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/Td3Jp/src/estimation/likelihoods.jl:3000 [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"#585#587"{Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, Pumas.var"#614#616"{FOCE, Bool, NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, TransformVariables.TransformTuple{NamedTuple{(:TVCL, :tvKTV, :TVVc, :ecmo, :cl, :v, :circuit, :Ω, :σ_prop_1, :σ_prop_2, :σ_add), Tuple{TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Float64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Int64}, TransformVariables.ScaledShiftedLogistic{Float64}, Pumas.PDiagTransform, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.Identity, TransformVariables.ShiftedExp{true, Float64}}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}}}, trace_simplex::Bool)

When I fit each dataset separately, the fit function works fine. I am not sure if the syntax in the @derived block should be modified. Any help would be appreciated.

Thanks

2
  1. What’s your initial estimate for sigma_prop_1?
  2. You have an upper bound of 1 on that parameter, is that intentional?
  3. Could you rewrite the model by constructing the covariate effect in the pre block and then passing it into the derived block? e.g. like this

siteeffect = site==1 ? sigma_prop_1 : sigma_prop_2

Then pass that siteeffect into the derived block

Best,
Vijay

3

Hi @vijay,

Thanks so much for your input, it was very helpful. Here is the modified code if a future user wants to implement it

mod = @model begin
    @param   begin
        TVCL      ∈ RealDomain(lower= 0.0001, upper = 6)
        tvKTV     ∈ RealDomain(lower = 0.001, upper = 2)
        TVVc      ∈ RealDomain(lower = 0.001, upper = 20) 
        ecmo      ∈ RealDomain(lower = -3,upper = 3)
        cl      ∈ RealDomain(lower = -3,upper = 3)
        v      ∈ RealDomain(lower = -3,upper = 3)
        circuit   ∈ RealDomain(lower = 0.001, upper = 5) 
        Ω         ∈ PDiagDomain(3)
        σ_prop_1    ∈ RealDomain(lower = 0.00001)
        σ_prop_2    ∈ RealDomain(lower = 0.00001)
        σ_add ∈ RealDomain(lower = 1)
    end
    @random begin
        η ~ MvNormal(Ω)
    end
    @covariates WT ECMO_DAYS IS_CIRCUIT_CHANGE site
    @pre begin
        siteeffect_1 = if site==2 
            σ_prop_1 
        else
            σ_prop_2
        end 
        siteeffect_2 = site==2 ? σ_add : 0.0
        KTV    = tvKTV*exp(η[1])*(ECMO_DAYS/6)^ecmo
        CL     = TVCL * exp(η[2])*(1-exp(-KTV*t))*(WT/10)^cl*circuit^IS_CIRCUIT_CHANGE
        Vc     = TVVc * exp(η[3])*(WT/10)^v
  
    end
    @dynamics  begin 
        Central' = -(CL/Vc)*Central
    end
    @derived begin
        CP = @. Central / Vc
        CONC = @. Normal(CP, sqrt(CP^2*siteeffect_1 + siteeffect_2)) 

    end
end
1 Like