How translate the individual model equation obtained in Monolix into Pumas

Estimation maximum-likelihood
1

Hi!
Would anyone help me translate the individual model equation obtained in Monolix into Pumas language? FC is my categorical covariable and the dose was given by iv.

[individual model – Model monolix]
log(Cl) = log(Cl_pop) + beta_Cl_FC_2*[FC = 2] + eta_Cl
log(V1) = log(V1_pop) + beta_V1_FC_2*[FC = 2] + eta_V1
log(Q) = log(Q_pop) + beta_Q_FC_2*[FC = 2] + eta_Q
log(V2) = log(V2_pop)

I tried to do this using the IF/ELSE command, based on pumas documentation, but some Gofs didn’t look good. In this case dFC2dcl sounds like beta values in monolix, right?

[Pumas model, 2 cmt model]

mdl_2cmp_prop = @ model begin
        @ param begin
            # model parameters
            tvcl ∈ RealDomain(lower = 0.001)
            tvv  ∈ RealDomain(lower = 0.001)
            tvvp ∈ RealDomain(lower = 0.001)
            tvq  ∈ RealDomain(lower = 0.001)
            Ω    ∈ PDiagDomain(3)
            σ ∈ RealDomain(lower = 0.001)   
            dFC1dcl ∈ RealDomain(lower = 0.001) # FC effect on parameter
            dFC2dcl ∈ RealDomain(lower = 0.001)
            dFC1dv  ∈ RealDomain(lower = 0.001)
            dFC2dv  ∈ RealDomain(lower = 0.001)
            dFC1dq  ∈ RealDomain(lower = 0.001)
            dFC2dq  ∈ RealDomain(lower = 0.001)

        end
      @ pre begin
#https://tutorials.pumas.ai/html/introduction/covariate.html#categorical-covariates
            CL = if FC == 1
                tvcl * exp(η[1]) * dFC1dcl 
            elseif FC == 2
                tvcl * exp(η[1]) * dFC2dcl 
            end
    
            Vc = if FC == 1
                tvv * exp(η[2]) * dFC1dv 
            elseif FC == 2
                tvv * exp(η[2]) * dFC2dv 
            end
    
            Q = if FC == 1
                tvq * exp(η[3]) * dFC1dq 
            elseif FC == 2
               tvq * exp(η[3]) * dFC2dq  
            end
            Vp = tvvp

        end

I appreciate it!

2

Dear @anafunguetto - welcome to the Pumas community. Glad to help you.

While your approach should work, you could start with a one to one translation of your monolix code and use something like this below.


@param begin
    ....
    beta_Cl_FC_2 in RealDomain()
    beta_V1_FC_2 in RealDomain()
    beta_Q_FC_2 in RealDomain()
    ....
end
@pre begin
    Cl = Cl_pop * beta_Cl_FC_2*(FC == 2) * exp(eta_Cl)
    V1 = V1_pop * beta_V1_FC_2*(FC == 2) * exp(eta_V1)
    Q = Q_pop   * beta_Q_FC_2*(FC == 2) * exp(eta_Q)
    V2 = V2_pop
end

This means that Cl_pop is the population value of clearance for FC==1, V1_pop is the population value of volume for FC==1 etc. Correspondingly beta_Cl_FC_2 is the multiplicative difference in population clearance between FC==1 and FC==2

Let me know if that is clear.

3

Thank you for your time! I did it, but now I have this error message in maximum LL estimation:

"ERROR: DomainError with Inf:
The initial parameter values cause a negative log likelihood of Inf. This can be due to model mis-specification or bad initial parameter values.

Stacktrace:

  1. _optim_check_initial_values(costf::Pumas.var"#621#625"{PumasModel{(tvcl = 1, tvv = 1, tvvp = 1, tvq = 1, Ω = 3, σ = 1, beta_Cl_FC_2 = 1, beta_V1_FC_2 = 1, beta_Q_FC_2 = 1), 3, (:Central, :Peripheral), (Symbol("##A##"), Symbol("##b##")), ParamSet{NamedTuple{(:tvcl, :tvv, :tvvp, :tvq, :Ω, :σ, :beta_Cl_FC_2, :beta_V1_FC_2, :beta_Q_FC_2), Tuple{RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}, RealDomain{Int64, TransformVariables.Infinity{true}, Int64}}}}, var"#81#88", Pumas.TimeDispatcher{var"#82#89", var"#83#90"}, Nothing, var"#85#92", Pumas.LinearODE, var"#86#93", var"#87#94", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Vector{Subject{NamedTuple{(:DV,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:FC,), Tuple{Int64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}, EnsembleThreads, NamedTuple{(), Tuple{}}, TransformVariables.TransformTuple{NamedTuple{(:tvcl, :tvv, :tvvp, :tvq, :Ω, :σ, :beta_Cl_FC_2, :beta_V1_FC_2, :beta_Q_FC_2), Tuple{TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, TransformVariables.ShiftedExp{true, Float64}, Pumas.PDiagTransform, Vararg{TransformVariables.ShiftedExp{true, Float64}, 4}}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Vector{Float64}}, vparam::Vector{Float64}, verbose::Bool) at likelihoods.jl"

I guess It is because of my initial values, so a changed but it didn’t work.
I’m working on this code now:

mdl_2cmp_prop = @ model begin
@ param begin
# Parâmetros do modelo
tvcl ∈ RealDomain(lower = 0.001)
tvv ∈ RealDomain(lower = 0.001)
tvvp ∈ RealDomain(lower = 0.001)
tvq ∈ RealDomain(lower = 0.001)
Ω ∈ PDiagDomain(3)
σ ∈ RealDomain(lower = 0.001)
beta_Cl_FC_2 in RealDomain(lower = -100)
beta_V1_FC_2 in RealDomain(lower = -100)
beta_Q_FC_2 in RealDomain(lower = -100)

@ random begin
η ~ MvNormal(Ω)
end

    @ covariates FC

    @ pre begin
        # https://tutorials.pumas.ai/html/introduction/covariate.html#categorical-covariates
        CL = tvcl * beta_Cl_FC_2*(FC == 2) * exp(η[1])
        Vc = tvv * beta_V1_FC_2*(FC == 2) * exp(η[2])
        Q = tvq   * beta_Q_FC_2*(FC == 2) * exp(η[3])
        Vp = tvvp

    end

    @ dynamics begin
        # https://docs.pumas.ai/stable/model_components/dynamical_types/
        Central'    = -(CL+Q)/Vc*Central + Q/Vp*Peripheral
        Peripheral' = Q/Vc*Central - Q/Vp*Peripheral

    end

    @ derived begin
        # https://docs.pumas.ai/stable/model_components/error_models/
        cp := @. Central / Vc
        DV ~ @. Normal(cp, cp * σ)
    end 
end

#initial estimates values
params_2cmp_prop = (
tvcl = 7119.54,
tvv = 1795.2,
tvq = 3471.5,
Ω = Diagonal([0.1, 0.1, 0.1]),
tvvp =7124.9,
σ = 0.1,
beta_Cl_FC_2 = 1,
beta_V1_FC_2 = 1,
beta_Q_FC_2 = 1,
)

maximum likelihood estimation

pkfit_2cmp_prop= fit(mdl_2cmp_prop, pop1, params_2cmp_prop, Pumas.FOCE())

thank you!

4

Can you try it with this instead?

    Cl = Cl_pop * exp(beta_Cl_FC_2*(FC == 2)) * exp(eta_Cl)
    V1 = V1_pop * exp(beta_V1_FC_2*(FC == 2)) * exp(eta_V1)
    Q = Q_pop   * exp(beta_Q_FC_2*(FC == 2)) * exp(eta_Q)