Issues with the code to create a three compartment model

Banshri · June 15, 2024, 8:41am

Hii dear Pumas team
I have been using the academic version of Pumas software for my dissertation work since March 2024.

I have to create a 3-compartment model Docetaxel for the same. I have been encountering several errors while running the code as listed below:

unable to run using Impute and when I import the package it would show me an error message stating pkg yet to be registered.
unable to use transform function as the error shows “transform not defined”
as the model is a 3-compartment model, I am not able to fit the base model as it shows an error

BoundsError: attempt to access 3-element Vector{ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag, Float64}, Float64, 3}} at index [4]

and not able to move further from there itself.

Kindly help me solve the issue as this is a crucial part of my dissertation thesis work.

vijay · June 15, 2024, 8:42am

Dear Banshri - welcome to the comunity.

Can I please request you to share your code? That way we can debug and see what is happening.

Best,
Vijay

Banshri · June 16, 2024, 1:55pm

Sure!

using Random
using CSV
using Pumas 
using DataFramesMeta 
using NCA
using NCAUtilities, NCA.Unitful 
using Dates 
using DataFrames
 using Chain
 using Pumas.Latexify
 using PumasUtilities 
 using CairoMakie
 using Serialization 
 using AlgebraOfGraphics, MakieCore

# impute col by group when all obs are missing
# for some groups it carries forward from
# the last observation
function locf_if_all_missing(col)
  if all(ismissing, col)
      return col
  else
      return Impute.locf(col)
  end
end

# impute col by group when all obs are missing
# for some groups it carries backward 
# the last observation
function nocb_if_all_missing(col)
  if all(ismissing, col)
      return col
  else
      return Impute.nocb(col)
  end
end

docedata = CSV.read("model_dataset_gedoce.csv", DataFrame; missingstring=["NA","","."])
#fill missing values in dose column with dose administered for that particular subject
filter!(df -> (df.id != 11), docedata)

transform!(
    groupby(docedata, :id),
    [:age, :gender, :bsa, :alb, :bil, :aag] .=>
    Impute.locf .=>
    [:age, :gender, :bsa, :alb, :bil, :aag]
)

transform!(
    groupby(docedata, :id),
    [:age, :sex, :bsa, :albumin, :bilirubin] .=>
    Impute.nocb .=>
    [:age, :sex, :bsa, :albumin, :bilirubin]
)

unique!(docedata, [:id, :time])
sort!(docedata, [:id, :time])

summ_docedf = filter(df -> (df.evid == 1), docedata)
describe(summ_docedf)


#...........................
pop_doce = read_pumas(  docedata;
id = :id,
time = :time,
observations = [:DV_1],
mdv = :mdv,
amt = :amt,
cmt = :cmt,
evid = :evid,
rate = :rate,
covariates = [:age, :gender, :bsa, :alb, :bil, :aag],
addl = :addl, 
ii = :ii,
)

observations_vs_time(pop_doce, loess = false, axis = (yscale = log10, xlabel = "Time (h)", ylabel = "Plasma docetaxel concentration (ug/L)"))
observations_vs_time(pop_doce)

doce_base = @model begin
  @param begin
    tvcl    ∈ RealDomain(lower=0, init = 18)
    tvvc    ∈ RealDomain(lower=0, init = 5.8)
    tvvp    ∈ RealDomain(lower=0, init = 2.65)
    tvvp1   ∈ RealDomain(lower=0, init = 468)
    tvq     ∈ RealDomain(lower=0, init = 6.49)
    tvq1    ∈ RealDomain(lower=0, init = 29.7)
    Ω ∈ PDiagDomain(; init = [0.18,0.25,0.48])
    σ_prop   ∈ RealDomain(lower = 0, init = 0.30)
    end

    @random begin
    η ~ MvNormal(Ω)
    end

    @pre begin
    CL    = tvcl*exp(η[1])
    Vc    = tvvc*exp(η[4])
    Vp    = tvvp
    Q     = tvq
    Vp1    = tvvp1*exp(η[3])
    Q1     = tvq1*exp(η[2])
    end
    
    @dynamics begin
      Central'    = -(CL+Q+Q1)/Vc*Central + Q/Vp*Peripheral + Q1/Vp1*Peripheral1
      Peripheral' = Q/Vc*Central - Q/Vp*Peripheral
      Peripheral1' = Q1/Vc*Central - Q1/Vp1*Peripheral1
      end
    
      @derived begin
        cp    = @. Central/Vc
        DV_1  ~ @. Normal(abs(cp), cp*σ_prop)
        end
end

doce_base_fit = fit(doce_base, pop_doce, init_params(doce_base), FOCE())
  insp_doce_base_fit = inspect(doce_base_fit)
goodness_of_fit(insp_doce_base_fit, loess = false)

icoefdata = icoef(doce_base_fit)
dficoef = DataFrame(icoefdata)
CSV.write("./Docetaxel_Banashri/dficoef.csv", dficoef)

#loglikelihood = 55.08
serialize("./Docetaxel_Banashri/Serialization/basemodel.jls",doce_base_fit)
basem = deserialize("./Docetaxel_Banashri/Serialization/basemodel.jls")

#infer(doce_iivq1_fit, Bootstrap(;samples = 200); level=0.95)

#BSA on CL
doce_bsacl = @model begin
  @param begin
    tvcl    ∈ RealDomain(lower=0, init = 18)
    tvvc    ∈ RealDomain(lower=0, init = 5.8)
    tvvp    ∈ RealDomain(lower=0, init = 2.65)
    tvvp1   ∈ RealDomain(lower=0, init = 468)
    tvq     ∈ RealDomain(lower=0, init = 6.49)
    tvq1    ∈ RealDomain(lower=0, init = 29.7)
    Ω       ∈ PDiagDomain(; init = [0.18,0.25,0.48])
    σ_prop  ∈ RealDomain(; init = 0.30)
    bsa_cl  ∈ RealDomain(; init = 0.99)
    end

    @random begin
    η ~ MvNormal(Ω)
    end

    @covariates bsa

    @pre begin
    CL    = tvcl*((bsa/1.59)^bsa_cl)*exp(η[1])
    Vc    = tvvc
    Vp    = tvvp
    Q     = tvq
    Vp1    = tvvp1*exp(η[3])
    Q1     = tvq1*exp(η[2])
    end
    
    @dynamics begin
      Central'    = -(CL+Q+Q1)/Vc*Central + Q/Vp*Peripheral + Q1/Vp1*Peripheral1
      Peripheral' = Q/Vc*Central - Q/Vp*Peripheral
      Peripheral1' = Q1/Vc*Central - Q1/Vp1*Peripheral1
      end
    
      @derived begin
        cp    = @. Central/Vc
        DV_1    ~ @. Normal(abs(cp), cp*σ_prop)
        end
end

doce_bsacl_fit = fit(doce_bsacl, pop_doce, init_params(doce_bsacl), FOCE())
insp_doce_bsacl_fit = inspect(doce_bsacl_fit)
goodness_of_fit(insp_doce_bsacl_fit, loess = false)

#loglikelihood = 54.79
serialize("./Docetaxel_Banashri/Serialization/bsacl.jls",doce_bsacl_fit)
bsacl_ = deserialize("./Docetaxel_Banashri/Serialization/bsacl.jls")

lrtest(basem,bsacl_)
#p = 1.0, statistic = -0.582
#BSA has no statistically significant effect on CL

#..........................
#age on CL
doce_agecl = @model begin
  @param begin
    tvcl    ∈ RealDomain(lower=0, init = 18)
    tvvc    ∈ RealDomain(lower=0, init = 5.8)
    tvvp    ∈ RealDomain(lower=0, init = 2.65)
    tvvp1   ∈ RealDomain(lower=0, init = 468)
    tvq     ∈ RealDomain(lower=0, init = 6.49)
    tvq1    ∈ RealDomain(lower=0, init = 29.7)
    Ω       ∈ PDiagDomain(; init = [0.18,0.25,0.48])
    σ_prop  ∈ RealDomain(; init = 0.30)
    age_cl  ∈ RealDomain(; init = -0.72)
    end

    @random begin
    η ~ MvNormal(Ω)
    end

    @covariates age

    @pre begin
    CL    = tvcl*((age/65)^age_cl)*exp(η[1])
    Vc    = tvvc
    Vp    = tvvp
    Q     = tvq
    Vp1    = tvvp1*exp(η[3])
    Q1     = tvq1*exp(η[2])
    end
    
    @dynamics begin
      Central'    = -(CL+Q+Q1)/Vc*Central + Q/Vp*Peripheral + Q1/Vp1*Peripheral1
      Peripheral' = Q/Vc*Central - Q/Vp*Peripheral
      Peripheral1' = Q1/Vc*Central - Q1/Vp1*Peripheral1
      end
    
      @derived begin
        cp    = @. Central/Vc
        DV_1    ~ @. Normal(abs(cp), cp*σ_prop)
        end
end

doce_agecl_fit = fit(doce_agecl, pop_doce, init_params(doce_agecl), FOCE())
insp_doce_agecl_fit = inspect(doce_agecl_fit)
goodness_of_fit(insp_doce_agecl_fit, loess = false)

#loglikelihood = 55.42
serialize("./Docetaxel_Banashri/Serialization/agecl.jls",doce_agecl_fit)
agecl_ = deserialize("./Docetaxel_Banashri/Serialization/agecl.jls")

lrtest(basem,agecl_)
#p = 0.412 (not significant), statistic = 0.673

#.............................
#alb on CL
doce_albcl = @model begin
  @param begin
    tvcl    ∈ RealDomain(lower=0, init = 18)
    tvvc    ∈ RealDomain(lower=0, init = 5.8)
    tvvp    ∈ RealDomain(lower=0, init = 2.65)
    tvvp1   ∈ RealDomain(lower=0, init = 468)
    tvq     ∈ RealDomain(lower=0, init = 6.49)
    tvq1    ∈ RealDomain(lower=0, init = 29.7)
    Ω       ∈ PDiagDomain(; init = [0.18,0.25,0.48])
    σ_prop  ∈ RealDomain(; init = 0.30)
    alb_cl  ∈ RealDomain(; init = 1.0)
    end

    @random begin
    η ~ MvNormal(Ω)
    end

    @covariates albumin

    @pre begin
    CL    = tvcl*((albumin/2)^alb_cl)*exp(η[1])
    Vc    = tvvc
    Vp    = tvvp
    Q     = tvq
    Vp1    = tvvp1*exp(η[3])
    Q1     = tvq1*exp(η[2])
    end
    
    @dynamics begin
      Central'    = -(CL+Q+Q1)/Vc*Central + Q/Vp*Peripheral + Q1/Vp1*Peripheral1
      Peripheral' = Q/Vc*Central - Q/Vp*Peripheral
      Peripheral1' = Q1/Vc*Central - Q1/Vp1*Peripheral1
      end
    
      @derived begin
        cp    = @. Central/Vc
        DV_1    ~ @. Normal(abs(cp), cp*σ_prop)
        end
end

doce_albcl_fit = fit(doce_albcl, pop_doce, init_params(doce_albcl), FOCE())
insp_doce_albcl_fit = inspect(doce_albcl_fit)
goodness_of_fit(insp_doce_albcl_fit, loess = false)

#loglikelihood = 48.51
serialize("./Docetaxel_Banashri/Serialization/albcl.jls",doce_albcl_fit)
albcl_ = deserialize("./Docetaxel_Banashri/Serialization/albcl.jls")

lrtest(basem,albcl_)
#p = 1.0 (No Significant effect of alb on CL)

#...............................
#bilirubin on CL

doce_bilcl = @model begin
  @param begin
    tvcl    ∈ RealDomain(lower=0, init = 18)
    tvvc    ∈ RealDomain(lower=0, init = 5.8)
    tvvp    ∈ RealDomain(lower=0, init = 2.65)
    tvvp1   ∈ RealDomain(lower=0, init = 468)
    tvq     ∈ RealDomain(lower=0, init = 6.49)
    tvq1    ∈ RealDomain(lower=0, init = 29.7)
    Ω       ∈ PDiagDomain(; init = [0.18,0.25,0.48])
    σ_prop  ∈ RealDomain(; init = 0.30)
    bil_cl  ∈ RealDomain(; init = 1.0)
    end

    @random begin
    η ~ MvNormal(Ω)
    end

    @covariates bilirubin

    @pre begin
    CL    = tvcl*((bilirubin/0.82)^bil_cl)*exp(η[1])
    Vc    = tvvc
    Vp    = tvvp
    Q     = tvq
    Vp1    = tvvp1*exp(η[3])
    Q1     = tvq1*exp(η[2])
    end
    
    @dynamics begin
      Central'     = -(CL+Q+Q1)/Vc*Central + Q/Vp*Peripheral + Q1/Vp1*Peripheral1
      Peripheral'  = Q/Vc*Central - Q/Vp*Peripheral
      Peripheral1' = Q1/Vc*Central - Q1/Vp1*Peripheral1
      end
    
      @derived begin
        cp    = @. Central/Vc
        DV_1    ~ @. Normal(abs(cp), cp*σ_prop)
        end
end

doce_bilcl_fit = fit(doce_bilcl, pop_doce, init_params(doce_bilcl), FOCE())
insp_doce_bilcl_fit = inspect(doce_bilcl_fit)
goodness_of_fit(insp_doce_bilcl_fit, loess = false)

#loglikelihood = 62.86
serialize("./Docetaxel_Banashri/Serialization/bilcl.jls",doce_bilcl_fit)
bilcl_ = deserialize("./Docetaxel_Banashri/Serialization/bilcl.jls")

lrtest(basem,bilcl_)
#statistic = 15.6, p value = 0.0
#p < 0.05 (bilirubin has a significant effect on CL)
#now doce_bilcl is our base model

#...........................
#sex on CL

doce_sexcl = @model begin
  @param begin
    tvcl    ∈ RealDomain(lower=0, init = 18)
    tvvc    ∈ RealDomain(lower=0, init = 5.8)
    tvvp    ∈ RealDomain(lower=0, init = 2.65)
    tvvp1   ∈ RealDomain(lower=0, init = 468)
    tvq     ∈ RealDomain(lower=0, init = 6.49)
    tvq1    ∈ RealDomain(lower=0, init = 29.7)
    Ω       ∈ PDiagDomain(; init = [0.18,0.25,0.48])
    σ_prop  ∈ RealDomain(; init = 0.30)
    bil_cl  ∈ RealDomain(; init = 1.0)
    tvsex_cl  ∈ RealDomain(; init = 1.0)
    end

    @random begin
    η ~ MvNormal(Ω)
    end

    @covariates bsa albumin bilirubin sex

    @pre begin
    sex_cl  = sex == 2 ? tvsex_cl : one(tvsex_cl)  #Sex=0, male , if 1-female
    CL    = tvcl*((bilirubin/3.22)^bil_cl)*sex_cl*exp(η[1])
    Vc    = tvvc
    Vp    = tvvp
    Q     = tvq
    Vp1    = tvvp1*exp(η[3])
    Q1     = tvq1*exp(η[2])
    end
    
    @dynamics begin
      Central'     = -(CL+Q+Q1)/Vc*Central + Q/Vp*Peripheral + Q1/Vp1*Peripheral1
      Peripheral'  = Q/Vc*Central - Q/Vp*Peripheral
      Peripheral1' = Q1/Vc*Central - Q1/Vp1*Peripheral1
      end
    
      @derived begin
        cp    = @. Central/Vc
        DV_1    ~ @. Normal(abs(cp), cp*σ_prop)
        end
end

doce_sexcl_fit = fit(doce_sexcl, pop_doce, init_params(doce_sexcl), FOCE())
insp_doce_sexcl_fit = inspect(doce_sexcl_fit)
goodness_of_fit(insp_doce_sexcl_fit, loess = false)
#loglikelihood = 62.88
serialize("./Docetaxel_Banashri/Serialization/sexcl.jls",doce_sexcl_fit)
sexcl_ = deserialize("./Docetaxel_Banashri/Serialization/sexcl.jls")

lrtest(bilcl_,sexcl_)
#p = 0.84, statistic = 0.04

#...........................
#inference table
infer(doce_bilcl_fit, Bootstrap(;samples = 200); level=0.95)

#creating a VPC plot
vpc_new = vpc(doce_bilcl_fit)
vpcplot = vpc_plot(vpc_new, observations = true, markersize = 8, simquantile_medians = true, simulated_linestyle = "-", simulated_color = :black, facet = true,
axis = (; xlabel = "Time (hr)", ylabel = "Predicted Concentration (ug/mL)"),)

#rich predictions
pred_param = (
  tvcl = 8.47,
  tvvc = 22.65,
  tvvp = 1.97,
  tvvp1 = 1105.6,
  tvq = 6.96,
  tvq1 = 42.78,
  Ω = Diagonal([0.11,0.54,0.80]),
  σ_prop = 0.45,
  bil_cl = -0.60
)

doce_pred = predict(doce_bilcl, pop_doce, pred_param, obstimes = 0.1:0.1:2016)
figures = subject_fits(
  doce_pred;
    separate=true,
    paginate=true,
    facet=(; combinelabels=true),
)

figures[1]
figures[2]
figures[3]
figures[4]
figures[5]
figures[6]

Banshri · June 17, 2024, 1:27pm

Hello Pumas community!
As requested I have posted my code.
Kindly let me know what changes are to be made.

Regards
Banshri

Topic		Replies	Views
Multiple errors while fitting model Estimation	6	351	June 4, 2022
Converting dataframe to pumas data frame Data read	3	620	August 14, 2020
Infer() error due to variance-covariance matrix calculating failed Estimation	3	248	March 7, 2023
Pumas version 0.12.0 (Alpha 12) released Announcements	0	618	May 14, 2020
Pumas version 0.14.0 (Pumas Beta v1.0) released Announcements	0	545	July 14, 2020

Issues with the code to create a three compartment model

Related topics