Multiple errors while fitting model

Hi! first and foremost I’d like to thank everyone at pumas.ai
this is my first time using pumas for a project so forgive my amateur mistakes :sweat_smile:
I’ve tried to fit my data set with 4 slightly different models just to get different error msgs such as :

ERROR: DomainError with x > shift must hold. Got
x => 0
shift => 0.0:

ERROR: DimensionMismatch(“array could not be broadcast to match destination”)

ERROR: BoundsError: attempt to access 2-element Vector{Int64} at index [3

I’m not sure what to do at this point but I’m guessing it could be because of my dataset? I changed all the observations from time zero to time 0.001 to avoid errors with no luck.
I’ll upload my dataset and I’m happy to share any of the codes I’ve tried out.

Could you please post the model code with respect to this error?

Also notice that in most of the datasets the dosing events are at time 0.0 and have a specific EVID, AMT and CMT.

I went back and checked the data . I added an evid column but no luck . still seeing errors.
one of the codes I used :

using CSV

using DataFramesMeta

using Dates

using Pumas

using Pumas.Latexify

using PumasUtilities

using Random

using CairoMakie

using Serialization

using LinearAlgebra

using Optim

pkdata = CSV.read(“/home/jrun/data/code/Lithium evid10.csv”, DataFrame ; missingstrings=[“NA”,“.”]);

pkdata_selected = @rsubset pkdata :OCC == 1

pop = read_pumas(pkdata_selected,

              id=:id,

              time=:time,

              amt=:amount,

              covariates = [:weight,:age,:sex],

              observations = [:observations],

              cmt = :cmt,

              evid = :evid)

pk_onecomp = @model begin

@metadata begin

desc = "One Compartment Model with lag time"

timeu = u"minute"

end

@param begin

"Absorption rate constant (1/min)"

tvka     ∈ RealDomain(lower=0)

"Elimination rate constant (1/min)"

tvkel    ∈ RealDomain(lower=0)

"Volume (L)"

tvvc     ∈ RealDomain(lower=0)

"Lag Time (mins)"

tvlag    ∈ RealDomain(lower=0)

Ω        ∈ PDiagDomain(4)

"Proportional RUV"

σ²_prop  ∈ RealDomain(lower=0)

end

@random begin

η        ~ MvNormal(Ω)

end

@pre begin

Ka       = tvka * exp(η[1])

Kel      = tvkel * exp(η[2])

Vc       = tvvc * exp(η[3])

end

@dosecontrol begin

lags     = (Depot=tvlag * exp(η[4]),)

end

@dynamics begin

Depot'   = -Ka*Depot

Central' =  Ka*Depot - Kel*Central

end

@derived begin

"""

PK02 Concentration (ug/L)

"""

cp       = @. Central/Vc

"""

PK02 Concentration (ug/L)

"""

observations   ~ @. Normal(cp, sqrt(cp^2*σ²_prop))

end

end

params_pk_onecomp = (tvka = 0.013,

tvkel = 0.013,

tvvc = 32,

tvlag = 0,

Ω = Diagonal([0.0,0.0,0.0,0.0]),

σ²_prop = 0.015)

ee_pk_onecomp = explore_estimates(pk_onecomp,

                             pop,

                             params_pk_onecomp)

                             latexify(pk_onecomp, :dynamics)

render(latexify(pk_onecomp, :dynamics))

Update coefficients based on app

params_pk_onecomp = coef(ee_pk_onecomp)

Create a DataFrame from the exploration

DataFrame(ee_pk_onecomp)

Evaluate initial loglikelihood

pkfit_pk_onecomp = loglikelihood(pk_onecomp,

pop,

params_pk_onecomp,

Pumas.FOCE())

explore influential individuals

pkfit_pk_onecomp = findinfluential(pk_onecomp,

pop,

params_pk_onecomp,

Pumas.FOCE())

Run model fitting

pkfit_pk_onecomp = fit(pk_onecomp,

pop,

params_pk_onecomp,

Pumas.FOCE())

Blockquote

which lead to this error :
ERROR: DomainError with x > shift must hold. Got
x => 0
shift => 0.0:

I’d like to upload my dataset as well but I can’t seem to figure out how to do it here.

another one i tried

using CSV

using DataFramesMeta

using Dates

using Pumas

using Pumas.Latexify

using PumasUtilities

using Random

using CairoMakie

using Serialization

using LinearAlgebra

using Optim

pkdata = CSV.read(“/home/jrun/data/code/Lithium evid10.csv”, DataFrame ; missingstrings=[“NA”,“.”]);

pkdata_selected = @rsubset pkdata :OCC == 1

pop = read_pumas(pkdata_selected,

              id=:id,

              time=:time,

              amt=:amount,

              covariates = [:weight,:age,:sex],

              observations = [:observations],

              cmt = :cmt,

              evid = :evid)

             

              pk_onecomp  = @model begin

               

                    @param begin

                      tvcl ∈ RealDomain(lower=0)

                      tvvc ∈ RealDomain(lower=0)

                      tvka ∈ RealDomain(lower=0)

                      Ω    ∈ PDiagDomain(3)

                    end

                 

                    @random begin

                      η ~ MvNormal(Ω)

                    end

                 

                    @covariates weight age sex

                 

                    @pre begin

                      CL = tvcl*(weight/70)^0.75*exp(η[1])

                      Vc = tvvc*(weight/70)*exp(η[2])

                      Ka = tvka*exp(η[3])

                    end

                 

                    @dynamics begin

                      Depot'   = -Ka*Depot

                      Central' =  Ka*Depot - (CL/Vc)*Central

                    end

                 

                    @observed begin

                      observations = @. Central/Vc

                    end

                  end

                  params_pk_onecomp = (

tvcl = 5, tvvc = 20, tvka = 1,

Ω = Diagonal([0.04, 0.04, 0.04]))

                    ee_pk_onecomp = explore_estimates(pk_onecomp,

                             pop,

                             params_pk_onecomp)

                             latexify(pk_onecomp, :dynamics)

render(latexify(pk_onecomp, :dynamics))

Update coefficients based on app

params_pk_onecomp = coef(ee_pk_onecomp)

Create a DataFrame from the exploration

DataFrame(ee_pk_onecomp)

Evaluate initial loglikelihood

pkfit_pk_onecomp = loglikelihood(pk_onecomp,

pop,

params_pk_onecomp,

Pumas.FOCE())

explore influential individuals

pkfit_pk_onecomp = findinfluential(pk_onecomp,

pop,

params_pk_onecomp,

Pumas.FOCE())

Run model fitting

pkfit_pk_onecomp = fit(pk_onecomp,

pop,

params_pk_onecomp,

Pumas.FOCE())

it gave this error :
ERROR: type NamedTuple has no field observations

i also used this model :

using CSV

using DataFramesMeta

using Dates

using Pumas

using Pumas.Latexify

using PumasUtilities

using Random

using CairoMakie

using Serialization

using LinearAlgebra

using Optim

pkdata = CSV.read(“/home/jrun/data/code/Lithium evid10.csv”, DataFrame ; missingstrings=[“NA”,“.”]);

pkdata_selected = @rsubset pkdata :OCC == 1

pop = read_pumas(pkdata_selected,

              id=:id,

              time=:time,

              amt=:amount,

              covariates = [:weight,:age,:sex],

              observations = [:observations],

              cmt = :cmt,

              evid = :evid)

              pk_onecomp  = @model begin

                @metadata begin

                  desc = "One Compartment Model "

                  timeu = u"hr"

                end

                @param begin

                    "Absorption rate constant (1/min)"

                    tvka    ∈ RealDomain(lower=0)

                    "Clearance (L/hr)"

                    tvcl   ∈ RealDomain(lower = 0.0001)

                    "Volume (L)"

                    tvvc   ∈ RealDomain(lower = 0.0001)

                    """

                    - ΩCL

                    - ΩVc

                    """

                    Ω      ∈ PDiagDomain(2) #covariance matrix

                     #σ_add  ∈ RealDomain(lower = 0.0001) # standard deviation

                     #σ_prop ∈ RealDomain(lower = 0.0001) # standard deviation

                    "Additive RUV"

                    σ²_add ∈ RealDomain(lower = 0.0001) # variance

                    "Proportional RUV"

                    σ²_prop ∈ RealDomain(lower = 0.0001) # variance  

                  end

                  @random begin

                    η ~ MvNormal(Ω)

                  end

                  @covariates weight age sex

                  @pre begin

                    CL = tvcl * exp(η[1])

                    Vc = tvvc * exp(η[2])

                    Ka       = tvka * exp(η[1])

                  end

               

                    @dynamics Central1



                    @derived begin

                      CONC := @. Central/Vc #suppressing CONC from being output using ":="

                      # additive error model

                      #DV ~ @. Normal(CONC, σ_add) # using standard deviation

                      #DV ~ @. Normal(CONC, sqrt(σ²_add)) # using variance

                      # proportional error model

                      #DV ~ @. Normal(CONC, abs(CONC)*σ_prop) # using standard deviation

                      #DV ~ @. Normal(CONC, sqrt(CONC^2*σ²_prop))   # using variance

                      # combination error model

                      """

                      DrugY Concentration (ng/mL)

                      """

                      observations ~ @. Normal(CONC, √(abs(CONC)^2*σ²_prop + σ²_add)) # using variance

                    end

                    @options begin

                      checklinear=false

                    end  

                end

                params_pk_onecomp = (

                  tvcl = 1.48, tvvc = 34, tvka = 0.64,

                  σ²_prop = 0.015 , σ²_add = 0.015,

                  Ω = Diagonal([0.04, 0.04, 0.04]))

                  ee_pk_onecomp = explore_estimates(pk_onecomp,

                           pop,

                           params_pk_onecomp)

                           latexify(pk_onecomp, :dynamics)

render(latexify(pk_onecomp, :dynamics))

Update coefficients based on app

params_pk_onecomp = coef(ee_pk_onecomp)

Create a DataFrame from the exploration

DataFrame(ee_pk_onecomp)

Evaluate initial loglikelihood

pkfit_pk_onecomp = loglikelihood(pk_onecomp,

pop,

params_pk_onecomp,

Pumas.FOCE())

explore influential individuals

pkfit_pk_onecomp = findinfluential(pk_onecomp,

pop,

params_pk_onecomp,

Pumas.FOCE())

and I got this error :
ERROR: DimensionMismatch(“array could not be broadcast to match destination”)

Your omega matrix is only two random effects while you are specifying three random effects in @pre block. you need to change

Ω      ∈ PDiagDomain(2) #covariance matrix     to      Ω      ∈ PDiagDomain(3) #covariance matrix

and

 @pre begin

                    CL = tvcl * exp(η[1])

                    Vc = tvvc * exp(η[2])

                    Ka       = tvka * exp(η[1])

                  end

to 
 @pre begin

                    CL = tvcl * exp(η[1])

                    Vc = tvvc * exp(η[2])

                    Ka       = tvka * exp(η[3])

                  end

Try to make lower bounds for PK parameters and residual error as small values such as 0.0001 instead of zero

thank you, I fixed that bit.
but I got an error which reads as ERROR: gradient is exactly zero in tvka. This indicates that tvka isn’t identified.
I made lower bounds and I don’t have zero values but rather 0.0001 . still I don’t get why it won’t work