Covariate Models

Authors

Jose Storopoli

Joel Owen 

using Dates
using Pumas
using PumasUtilities
using CairoMakie
using DataFramesMeta
using CSV
using PharmaDatasets
Caution

Some functions in this tutorial are only available after you load the PumasUtilities package.

In this tutorial we’ll cover a workflow around covariate model building. You’ll learn how to:

  1. include covariates in your model
  2. parse data with covariates
  3. understand the difference between constant and time-varying covariates
  4. handle continuous and categorical covariates
  5. deal with missing data in your covariates
  6. deal with categorical covariates

1 nlme_sample Dataset

For this tutorial we’ll use the nlme_sample dataset from PharmaDatasets.jl:

pkfile = dataset("nlme_sample", String)
pkdata = CSV.read(pkfile, DataFrame; missingstring = ["NA", ""])
first(pkdata, 5)
5×15 DataFrame
Row ID TIME DV AMT EVID CMT RATE WT AGE SEX CRCL GROUP ROUTE DURATION OCC
Int64 Float64 Float64? Int64? Int64 Int64? Int64 Int64 Int64 String1 Int64 String7 Float64 Int64? Int64
1 1 0.0 missing 1000 1 1 500 90 47 M 75 1000 mg Inf 2 1
2 1 0.001 0.0667231 missing 0 missing 0 90 47 M 75 1000 mg Inf missing 1
3 1 1.0 112.817 missing 0 missing 0 90 47 M 75 1000 mg Inf missing 1
4 1 2.0 224.087 missing 0 missing 0 90 47 M 75 1000 mg Inf missing 1
5 1 4.0 220.047 missing 0 missing 0 90 47 M 75 1000 mg Inf missing 1
Note

The nlme_sample dataset has different missing values as the standard datasets in the PharmaDatasets.jl. That’s why we are first getting a String representation of the dataset as a CSV file as pkfile variable. Then, we use it to customize the missingstring keyword argument inside CSV.read function in order to have a working DataFrame for the nlme_sample dataset.

If you want to know more about data wrangling and how to read and write data in different formats, please check out the Data Wrangling Tutorials at tutorials.pumas.ai.

As you can see, the nlme_sample dataset has the standard PK dataset columns such as :ID, :TIME, :DV, :AMT, :EVID and :CMT. The dataset also contains the following list of covariates:

  • :WT: subject weight in kilograms
  • :SEX: subject sex, either "F" or "M"
  • :CRCL: subject creatinine clearance
  • :GROUP: subject dosing group, either "500 mg", "750 mg", or "1000 mg"

And we’ll learn how to include them in our Pumas modeling workflows.

describe(pkdata, :mean, :std, :nunique, :first, :last, :eltype)
15×7 DataFrame
Row variable mean std nunique first last eltype
Symbol Union… Union… Union… Any Any Type
1 ID 15.5 8.661 1 30 Int64
2 TIME 82.6527 63.2187 0.0 168.0 Float64
3 DV 157.315 110.393 missing missing Union{Missing, Float64}
4 AMT 750.0 204.551 1000 500 Union{Missing, Int64}
5 EVID 0.307692 0.461835 1 1 Int64
6 CMT 1.0 0.0 1 1 Union{Missing, Int64}
7 RATE 115.385 182.218 500 250 Int64
8 WT 81.6 11.6051 90 96 Int64
9 AGE 40.0333 11.6479 47 56 Int64
10 SEX 2 M F String1
11 CRCL 72.5667 26.6212 75 90 Int64
12 GROUP 3 1000 mg 500 mg String7
13 ROUTE Inf NaN Inf Inf Float64
14 DURATION 2.0 0.0 2 2 Union{Missing, Int64}
15 OCC 4.15385 2.62836 1 8 Int64
Tip

As you can see, all these covariates are constant. That means, they do not vary over time. We’ll also cover time-varying covariates later in this tutorial.

2 Step 1 - Parse Data into a Population

The first step in our covariate model building workflow is to parse data into a Population. This is accomplished with the read_pumas function. Here we are to use the covariates keyword argument to pass a vector of column names to be parsed as covariates:

pop = read_pumas(
    pkdata;
    id = :ID,
    time = :TIME,
    amt = :AMT,
    covariates = [:WT, :AGE, :SEX, :CRCL, :GROUP],
    observations = [:DV],
    cmt = :CMT,
    evid = :EVID,
    rate = :RATE,
)
Population
  Subjects: 30
  Covariates: WT, AGE, SEX, CRCL, GROUP
  Observations: DV

Due to Pumas’ dynamic workflow capabilities, we only need to define our Population once. That is, we tell read_pumas to parse all the covariates available, even if we will be using none or a subset of those in our models.

This is a feature that highly increases workflow efficiency in developing and fitting models.

3 Step 2 - Base Model

The second step of our covariate model building workflow is to develop a base model, i.e., a model without any covariate effects on its parameters. This represents the null model against which covariate models can be tested after checking if covariate inclusion is helpful in our model.

Let’s create a combined-error simple 2-compartment base model:

base_model = @model begin
    @param begin
        tvcl ∈ RealDomain(; lower = 0)
        tvvc ∈ RealDomain(; lower = 0)
        tvq ∈ RealDomain(; lower = 0)
        tvvp ∈ RealDomain(; lower = 0)
        Ω ∈ PDiagDomain(2)
        σ_add ∈ RealDomain(; lower = 0)
        σ_prop ∈ RealDomain(; lower = 0)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @pre begin
        CL = tvcl * exp(η[1])
        Vc = tvvc * exp(η[2])
        Q = tvq
        Vp = tvvp
    end

    @dynamics Central1Periph1

    @derived begin
        cp := @. Central / Vc
        DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
    end
end
PumasModel
  Parameters: tvcl, tvvc, tvq, tvvp, Ω, σ_add, σ_prop
  Random effects: η
  Covariates:
  Dynamical system variables: Central, Peripheral
  Dynamical system type: Closed form
  Derived: DV
  Observed: DV

And fit it accordingly:

iparams_base_model = (;
    tvvc = 5,
    tvcl = 0.02,
    tvq = 0.01,
    tvvp = 10,
    Ω = Diagonal([0.01, 0.01]),
    σ_add = 0.1,
    σ_prop = 0.1,
)
(tvvc = 5,
 tvcl = 0.02,
 tvq = 0.01,
 tvvp = 10,
 Ω = [0.01 0.0; 0.0 0.01],
 σ_add = 0.1,
 σ_prop = 0.1,)
fit_base_model = fit(base_model, pop, iparams_base_model, FOCE())
[ 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     8.300164e+03     4.360770e+03
 * time: 0.025621891021728516
     1     3.110315e+03     9.706222e+02
 * time: 1.0510759353637695
     2     2.831659e+03     7.817006e+02
 * time: 1.2130820751190186
     3     2.405281e+03     2.923696e+02
 * time: 1.2819879055023193
     4     2.370406e+03     3.032286e+02
 * time: 1.3345019817352295
     5     2.313631e+03     3.126188e+02
 * time: 1.3888940811157227
     6     2.263986e+03     2.946697e+02
 * time: 1.437803030014038
     7     2.160182e+03     1.917599e+02
 * time: 1.5090949535369873
     8     2.112467e+03     1.588288e+02
 * time: 1.5737171173095703
     9     2.090339e+03     1.108334e+02
 * time: 1.618088960647583
    10     2.078171e+03     8.108027e+01
 * time: 1.6669321060180664
    11     2.074517e+03     7.813928e+01
 * time: 1.7105250358581543
    12     2.066270e+03     7.044632e+01
 * time: 1.7542519569396973
    13     2.049660e+03     1.062584e+02
 * time: 1.8558480739593506
    14     2.021965e+03     1.130570e+02
 * time: 1.8961169719696045
    15     1.994936e+03     7.825801e+01
 * time: 1.935425043106079
    16     1.979337e+03     5.318263e+01
 * time: 1.9743189811706543
    17     1.972141e+03     6.807046e+01
 * time: 2.0127758979797363
    18     1.967973e+03     7.896361e+01
 * time: 2.050477981567383
    19     1.962237e+03     8.343757e+01
 * time: 2.0900189876556396
    20     1.952791e+03     5.565304e+01
 * time: 2.130244016647339
    21     1.935857e+03     3.923284e+01
 * time: 2.171376943588257
    22     1.926254e+03     5.749643e+01
 * time: 2.211822032928467
    23     1.922144e+03     4.306225e+01
 * time: 2.2494490146636963
    24     1.911566e+03     4.810496e+01
 * time: 2.2914700508117676
    25     1.906464e+03     4.324267e+01
 * time: 2.336946964263916
    26     1.905339e+03     1.207954e+01
 * time: 2.398267984390259
    27     1.905092e+03     1.093479e+01
 * time: 2.43369197845459
    28     1.904957e+03     1.057034e+01
 * time: 2.4687340259552
    29     1.904875e+03     1.060882e+01
 * time: 2.5022690296173096
    30     1.904459e+03     1.031525e+01
 * time: 2.5377860069274902
    31     1.903886e+03     1.041179e+01
 * time: 2.5754880905151367
    32     1.903313e+03     1.135672e+01
 * time: 2.6124370098114014
    33     1.903057e+03     1.075683e+01
 * time: 2.647097110748291
    34     1.902950e+03     1.091295e+01
 * time: 2.682292938232422
    35     1.902887e+03     1.042409e+01
 * time: 2.716481924057007
    36     1.902640e+03     9.213043e+00
 * time: 2.7520220279693604
    37     1.902364e+03     9.519321e+00
 * time: 2.7880680561065674
    38     1.902156e+03     5.590984e+00
 * time: 2.8246970176696777
    39     1.902094e+03     1.811898e+00
 * time: 2.866909980773926
    40     1.902086e+03     1.644722e+00
 * time: 2.9269349575042725
    41     1.902084e+03     1.653520e+00
 * time: 2.960542917251587
    42     1.902074e+03     1.720184e+00
 * time: 2.994633913040161
    43     1.902055e+03     2.619061e+00
 * time: 3.0286879539489746
    44     1.902015e+03     3.519729e+00
 * time: 3.0628750324249268
    45     1.901962e+03     3.403372e+00
 * time: 3.098323106765747
    46     1.901924e+03     1.945644e+00
 * time: 3.131861925125122
    47     1.901914e+03     6.273342e-01
 * time: 3.1654140949249268
    48     1.901913e+03     5.374557e-01
 * time: 3.198422908782959
    49     1.901913e+03     5.710007e-01
 * time: 3.230139970779419
    50     1.901913e+03     6.091390e-01
 * time: 3.262113094329834
    51     1.901912e+03     6.692417e-01
 * time: 3.2949130535125732
    52     1.901909e+03     7.579620e-01
 * time: 3.330374002456665
    53     1.901903e+03     8.798211e-01
 * time: 3.366763114929199
    54     1.901889e+03     1.002981e+00
 * time: 3.4236209392547607
    55     1.901864e+03     1.495512e+00
 * time: 3.4584569931030273
    56     1.901837e+03     1.664621e+00
 * time: 3.491044044494629
    57     1.901819e+03     8.601118e-01
 * time: 3.5245680809020996
    58     1.901815e+03     4.525470e-02
 * time: 3.5580921173095703
    59     1.901815e+03     1.294280e-02
 * time: 3.590919017791748
    60     1.901815e+03     2.876568e-03
 * time: 3.620661973953247
    61     1.901815e+03     2.876568e-03
 * time: 3.705904960632324
    62     1.901815e+03     2.876568e-03
 * time: 3.7934229373931885
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
Likelihood Optimizer:                         BFGS
Dynamical system type:                 Closed form

Log-likelihood value:                    -1901.815
Number of subjects:                             30
Number of parameters:         Fixed      Optimized
                                  0              8
Observation records:         Active        Missing
    DV:                         540              0
    Total:                      540              0

---------------------
           Estimate
---------------------
tvcl        0.1542
tvvc        4.5856
tvq         0.042341
tvvp        3.7082
Ω₁,₁        0.26467
Ω₂,₂        0.10627
σ_add       0.032183
σ_prop      0.061005
---------------------

4 Step 3 - Covariate Model

The third step of our covariate model building workflow is to actually develop one or more covariate models. Let’s develop three covariate models:

  1. allometric scaling based on weight
  2. clearance effect based on creatinine clearance
  3. separated error model based on sex

To include covariates in a Pumas model we need to first include them in the @covariates block. Then, we are free to use them inside the @pre block

Here’s our covariate model with allometric scaling based on weight:

Tip

When building covariate models, unlike in NONMEM, it is highly recommended to derive covariates or perform any required transformations or scaling as a data wrangling step and pass the derived/transformed into read_pumas as a covariate. In this particular case, it is easy to create two columns in the original data as:

@rtransform! pkdata :ASWT_CL = (:WT / 70)^0.75
@rtransform! pkdata :ASWT_Vc = (:WT / 70)

Then, you bring these covariates into read_pumas and use them directly on CL and Vc. This will avoid computations of your transformations during the model fit.

covariate_model_wt = @model begin
    @param begin
        tvcl ∈ RealDomain(; lower = 0)
        tvvc ∈ RealDomain(; lower = 0)
        tvq ∈ RealDomain(; lower = 0)
        tvvp ∈ RealDomain(; lower = 0)
        Ω ∈ PDiagDomain(2)
        σ_add ∈ RealDomain(; lower = 0)
        σ_prop ∈ RealDomain(; lower = 0)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates begin
        WT
    end

    @pre begin
        CL = tvcl * (WT / 70)^0.75 * exp(η[1])
        Vc = tvvc * (WT / 70) * exp(η[2])
        Q = tvq
        Vp = tvvp
    end

    @dynamics Central1Periph1

    @derived begin
        cp := @. Central / Vc
        DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
    end
end
PumasModel
  Parameters: tvcl, tvvc, tvq, tvvp, Ω, σ_add, σ_prop
  Random effects: η
  Covariates: WT
  Dynamical system variables: Central, Peripheral
  Dynamical system type: Closed form
  Derived: DV
  Observed: DV

Once we included the WT covariate in the @covariates block we can use the WT values inside the @pre block. For both clearance (CL) and volume of the central compartment (Vc), we are allometric scaling by the WT value by the mean weight 70 and, in the case of CL using an allometric exponent with value 0.75.

Let’s fit our covariate_model_wt. Notice that we have not added any new parameters inside the @param block, thus, we can use the same iparams_base_model initial parameters values’ list:

fit_covariate_model_wt = fit(covariate_model_wt, pop, iparams_base_model, FOCE())
[ 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     7.695401e+03     4.898919e+03
 * time: 6.890296936035156e-5
     1     2.846050e+03     1.128657e+03
 * time: 0.08052682876586914
     2     2.472982e+03     7.008264e+02
 * time: 0.13461995124816895
     3     2.180690e+03     2.312704e+02
 * time: 0.18978190422058105
     4     2.125801e+03     1.862929e+02
 * time: 0.23441481590270996
     5     2.103173e+03     1.320946e+02
 * time: 0.349808931350708
     6     2.091136e+03     1.103035e+02
 * time: 0.38991284370422363
     7     2.081443e+03     1.091133e+02
 * time: 0.42706894874572754
     8     2.071793e+03     7.197675e+01
 * time: 0.46520495414733887
     9     2.062706e+03     7.623310e+01
 * time: 0.5034699440002441
    10     2.057515e+03     6.885476e+01
 * time: 0.541003942489624
    11     2.051133e+03     6.368504e+01
 * time: 0.5771958827972412
    12     2.038626e+03     7.730243e+01
 * time: 0.6136679649353027
    13     2.019352e+03     1.136864e+02
 * time: 0.6523489952087402
    14     1.997136e+03     1.005748e+02
 * time: 0.6909389495849609
    15     1.983023e+03     6.831478e+01
 * time: 0.7352168560028076
    16     1.977700e+03     5.649783e+01
 * time: 0.7798378467559814
    17     1.974583e+03     6.357642e+01
 * time: 0.8248069286346436
    18     1.967292e+03     7.658974e+01
 * time: 0.9084019660949707
    19     1.951045e+03     6.130573e+01
 * time: 0.959183931350708
    20     1.935868e+03     4.845839e+01
 * time: 1.025730848312378
    21     1.929356e+03     6.325782e+01
 * time: 1.0863349437713623
    22     1.925187e+03     3.142245e+01
 * time: 1.1253979206085205
    23     1.923733e+03     4.623400e+01
 * time: 1.1632418632507324
    24     1.918498e+03     5.347738e+01
 * time: 1.2014379501342773
    25     1.912383e+03     5.849125e+01
 * time: 1.244764804840088
    26     1.905510e+03     3.254038e+01
 * time: 1.2863948345184326
    27     1.903629e+03     2.905618e+01
 * time: 1.323462963104248
    28     1.902833e+03     2.907696e+01
 * time: 1.3623178005218506
    29     1.902447e+03     2.746037e+01
 * time: 1.3999559879302979
    30     1.899399e+03     1.930949e+01
 * time: 1.4423038959503174
    31     1.898705e+03     1.186361e+01
 * time: 1.5081398487091064
    32     1.898505e+03     1.050402e+01
 * time: 1.5439999103546143
    33     1.898474e+03     1.042186e+01
 * time: 1.5754528045654297
    34     1.897992e+03     1.238729e+01
 * time: 1.6078779697418213
    35     1.897498e+03     1.729368e+01
 * time: 1.6418910026550293
    36     1.896954e+03     1.472554e+01
 * time: 1.6753458976745605
    37     1.896744e+03     5.852825e+00
 * time: 1.7081198692321777
    38     1.896705e+03     1.171353e+00
 * time: 1.737980842590332
    39     1.896704e+03     1.216117e+00
 * time: 1.7688119411468506
    40     1.896703e+03     1.230336e+00
 * time: 1.7994458675384521
    41     1.896698e+03     1.250715e+00
 * time: 1.830901861190796
    42     1.896688e+03     1.449552e+00
 * time: 1.8623628616333008
    43     1.896666e+03     2.533300e+00
 * time: 1.8960678577423096
    44     1.896631e+03     3.075536e+00
 * time: 1.9325978755950928
    45     1.896599e+03     2.139797e+00
 * time: 1.9917538166046143
    46     1.896587e+03     6.636031e-01
 * time: 2.0249388217926025
    47     1.896585e+03     6.303517e-01
 * time: 2.0566539764404297
    48     1.896585e+03     5.995265e-01
 * time: 2.086967945098877
    49     1.896584e+03     5.844159e-01
 * time: 2.1187939643859863
    50     1.896583e+03     6.083858e-01
 * time: 2.151124954223633
    51     1.896579e+03     8.145326e-01
 * time: 2.1826817989349365
    52     1.896570e+03     1.293490e+00
 * time: 2.2141737937927246
    53     1.896549e+03     1.877705e+00
 * time: 2.2457098960876465
    54     1.896513e+03     2.217391e+00
 * time: 2.277435779571533
    55     1.896482e+03     1.658147e+00
 * time: 2.3091259002685547
    56     1.896466e+03     5.207215e-01
 * time: 2.3410727977752686
    57     1.896463e+03     1.177468e-01
 * time: 2.373464822769165
    58     1.896463e+03     1.678937e-02
 * time: 2.405090808868408
    59     1.896463e+03     2.666635e-03
 * time: 2.4585328102111816
    60     1.896463e+03     2.666635e-03
 * time: 2.5200469493865967
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
Likelihood Optimizer:                         BFGS
Dynamical system type:                 Closed form

Log-likelihood value:                   -1896.4632
Number of subjects:                             30
Number of parameters:         Fixed      Optimized
                                  0              8
Observation records:         Active        Missing
    DV:                         540              0
    Total:                      540              0

---------------------
           Estimate
---------------------
tvcl        0.13915
tvvc        3.9754
tvq         0.041988
tvvp        3.5722
Ω₁,₁        0.23874
Ω₂,₂        0.081371
σ_add       0.032174
σ_prop      0.061012
---------------------

We can definitely see that, despite not increasing the parameters of the model, our loglikelihood is higher (implying lower objective function value). Furthermore, our additive and combined error values remain almost the same while our Ωs decreased for CL and Vc. This implies that the WT covariate is definitely assisting in a better model fit by capturing the variability in CL and Vc. We’ll compare models in the next step.

Let’s now try to incorporate into this model creatinine clearance (CRCL) effect on clearance (CL):

Tip

Like the tip above, you can derive covariates or perform any required transformations or scaling as a data wrangling step and pass the derived/transformed into read_pumas as a covariate. In this particular case, it is easy to create three more columns in the original data as:

@rtransform! pkdata :CRCL_CL = :CRCL / 100
@rtransform! pkdata :ASWT_CL = (:WT / 70)^0.75
@rtransform! pkdata :ASWT_Vc = (:WT / 70)

Then, you bring these covariates into read_pumas and use them directly on renCL, CL and Vc. This will avoid computations of your transformations during the model fit.

covariate_model_wt_crcl = @model begin
    @param begin
        tvvc ∈ RealDomain(; lower = 0)
        tvq ∈ RealDomain(; lower = 0)
        tvvp ∈ RealDomain(; lower = 0)
        tvcl_hep ∈ RealDomain(; lower = 0)
        tvcl_ren ∈ RealDomain(; lower = 0)
        Ω ∈ PDiagDomain(2)
        σ_add ∈ RealDomain(; lower = 0)
        σ_prop ∈ RealDomain(; lower = 0)
        dCRCL ∈ RealDomain()
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates begin
        WT
        CRCL
    end

    @pre begin
        hepCL = tvcl_hep * (WT / 70)^0.75
        renCL = tvcl_ren * (CRCL / 100)^dCRCL
        CL = (hepCL + renCL) * exp(η[1])
        Vc = tvvc * (WT / 70) * exp(η[2])
        Q = tvq
        Vp = tvvp
    end

    @dynamics Central1Periph1

    @derived begin
        cp := @. Central / Vc
        DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
    end
end
PumasModel
  Parameters: tvvc, tvq, tvvp, tvcl_hep, tvcl_ren, Ω, σ_add, σ_prop, dCRCL
  Random effects: η
  Covariates: WT, CRCL
  Dynamical system variables: Central, Peripheral
  Dynamical system type: Closed form
  Derived: DV
  Observed: DV

In the covariate_model_wt_crcl model we are keeping our allometric scaling on WT from before. But we are also adding a new covariate creatinine clearance (CRCL), dividing clearance (CL) into hepatic (hepCL) and renal clearance (renCL), along with a new parameter dCRCL.

dCRCL is the exponent of the power function for the effect of creatinine clearance on renal clearance. In some models this parameter is fixed, however we’ll allow the model to estimate it.

This is a good example on how to add covariate coefficients such as dCRCL in any Pumas covariate model. Now, let’s fit this model. Note that we need a new initial parameters values’ list since the previous one we used doesn’t include dCRCL, tvcl_hep or tvcl_ren:

iparams_covariate_model_wt_crcl = (;
    tvvc = 5,
    tvcl_hep = 0.01,
    tvcl_ren = 0.01,
    tvq = 0.01,
    tvvp = 10,
    Ω = Diagonal([0.01, 0.01]),
    σ_add = 0.1,
    σ_prop = 0.1,
    dCRCL = 0.9,
)
(tvvc = 5,
 tvcl_hep = 0.01,
 tvcl_ren = 0.01,
 tvq = 0.01,
 tvvp = 10,
 Ω = [0.01 0.0; 0.0 0.01],
 σ_add = 0.1,
 σ_prop = 0.1,
 dCRCL = 0.9,)
fit_covariate_model_wt_crcl =
    fit(covariate_model_wt_crcl, pop, iparams_covariate_model_wt_crcl, FOCE())
[ 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     8.554152e+03     5.650279e+03
 * time: 6.699562072753906e-5
     1     3.572050e+03     1.302046e+03
 * time: 0.2336139678955078
     2     3.266947e+03     5.384244e+02
 * time: 0.29645705223083496
     3     3.150635e+03     1.918079e+02
 * time: 0.3518071174621582
     4     3.108122e+03     1.277799e+02
 * time: 0.40232205390930176
     5     3.091358e+03     8.838080e+01
 * time: 0.4514639377593994
     6     3.082997e+03     8.321689e+01
 * time: 0.500093936920166
     7     3.076379e+03     8.167702e+01
 * time: 0.5494511127471924
     8     3.067422e+03     1.177822e+02
 * time: 0.6019189357757568
     9     3.048580e+03     2.526969e+02
 * time: 0.661107063293457
    10     2.993096e+03     6.325396e+02
 * time: 0.7535440921783447
    11     2.965744e+03     7.634718e+02
 * time: 1.0046489238739014
    12     2.921212e+03     9.704020e+02
 * time: 1.1059620380401611
    13     2.553649e+03     6.642510e+02
 * time: 1.1847500801086426
    14     2.319495e+03     3.204711e+02
 * time: 1.3167660236358643
    15     2.183040e+03     2.174905e+02
 * time: 1.421375036239624
    16     2.157621e+03     3.150983e+02
 * time: 1.4783430099487305
    17     2.132395e+03     2.847614e+02
 * time: 1.5323691368103027
    18     2.084799e+03     1.563370e+02
 * time: 1.5849881172180176
    19     2.071497e+03     1.006429e+02
 * time: 1.6786539554595947
    20     2.064983e+03     9.753313e+01
 * time: 1.727445125579834
    21     2.048289e+03     9.230405e+01
 * time: 1.7754809856414795
    22     2.020938e+03     1.292359e+02
 * time: 1.8222119808197021
    23     1.983888e+03     2.284990e+02
 * time: 1.8713979721069336
    24     1.962132e+03     1.220188e+02
 * time: 1.9189250469207764
    25     1.945947e+03     1.035894e+02
 * time: 1.9631850719451904
    26     1.917782e+03     8.246698e+01
 * time: 2.0145421028137207
    27     1.905967e+03     5.408054e+01
 * time: 2.0671679973602295
    28     1.898569e+03     2.172222e+01
 * time: 2.1201059818267822
    29     1.897473e+03     1.689350e+01
 * time: 2.171347141265869
    30     1.897019e+03     1.666689e+01
 * time: 2.2229740619659424
    31     1.896796e+03     1.699751e+01
 * time: 2.274718999862671
    32     1.896559e+03     1.645900e+01
 * time: 2.3637280464172363
    33     1.896223e+03     1.415504e+01
 * time: 2.407750129699707
    34     1.895773e+03     1.630081e+01
 * time: 2.452207088470459
    35     1.895309e+03     1.723930e+01
 * time: 2.497921943664551
    36     1.895004e+03     1.229983e+01
 * time: 2.5436911582946777
    37     1.894871e+03     5.385102e+00
 * time: 2.588304042816162
    38     1.894827e+03     3.465463e+00
 * time: 2.6330151557922363
    39     1.894816e+03     3.387474e+00
 * time: 2.6799120903015137
    40     1.894807e+03     3.295388e+00
 * time: 2.728188991546631
    41     1.894786e+03     3.089194e+00
 * time: 2.777161121368408
    42     1.894737e+03     2.928080e+00
 * time: 2.826712131500244
    43     1.894624e+03     3.088723e+00
 * time: 2.876559019088745
    44     1.894413e+03     3.493791e+00
 * time: 2.927276134490967
    45     1.894127e+03     3.142865e+00
 * time: 3.015350103378296
    46     1.893933e+03     2.145253e+00
 * time: 3.0607969760894775
    47     1.893888e+03     2.172800e+00
 * time: 3.1054229736328125
    48     1.893880e+03     2.180509e+00
 * time: 3.1487860679626465
    49     1.893876e+03     2.134257e+00
 * time: 3.1914479732513428
    50     1.893868e+03     2.032354e+00
 * time: 3.2356481552124023
    51     1.893846e+03     1.760874e+00
 * time: 3.2791149616241455
    52     1.893796e+03     1.779016e+00
 * time: 3.3269240856170654
    53     1.893694e+03     2.018956e+00
 * time: 3.3759100437164307
    54     1.893559e+03     2.366854e+00
 * time: 3.4247360229492188
    55     1.893474e+03     3.690142e+00
 * time: 3.47455096244812
    56     1.893446e+03     3.675109e+00
 * time: 3.523712158203125
    57     1.893439e+03     3.426419e+00
 * time: 3.5726230144500732
    58     1.893429e+03     3.183164e+00
 * time: 3.6607561111450195
    59     1.893398e+03     2.695171e+00
 * time: 3.7047131061553955
    60     1.893328e+03     2.753548e+00
 * time: 3.757009983062744
    61     1.893169e+03     3.589748e+00
 * time: 3.83332896232605
    62     1.892920e+03     3.680718e+00
 * time: 3.8928189277648926
    63     1.892667e+03     2.568107e+00
 * time: 3.93630313873291
    64     1.892514e+03     1.087910e+00
 * time: 3.9801180362701416
    65     1.892493e+03     3.287296e-01
 * time: 4.026593923568726
    66     1.892492e+03     2.967465e-01
 * time: 4.073759078979492
    67     1.892492e+03     3.020682e-01
 * time: 4.119004964828491
    68     1.892491e+03     3.034704e-01
 * time: 4.161945104598999
    69     1.892491e+03     3.091846e-01
 * time: 4.205899953842163
    70     1.892491e+03     3.224170e-01
 * time: 4.250786066055298
    71     1.892490e+03     6.494197e-01
 * time: 4.29577898979187
    72     1.892488e+03     1.115188e+00
 * time: 4.376759052276611
    73     1.892483e+03     1.838833e+00
 * time: 4.4189629554748535
    74     1.892472e+03     2.765371e+00
 * time: 4.461066961288452
    75     1.892452e+03     3.463807e+00
 * time: 4.503055095672607
    76     1.892431e+03     2.805270e+00
 * time: 4.544783115386963
    77     1.892411e+03     5.758916e-01
 * time: 4.587821006774902
    78     1.892410e+03     1.434041e-01
 * time: 4.6302330493927
    79     1.892409e+03     1.639246e-01
 * time: 4.677207946777344
    80     1.892409e+03     1.145856e-01
 * time: 4.724751949310303
    81     1.892409e+03     3.966861e-02
 * time: 4.771382093429565
    82     1.892409e+03     3.550808e-02
 * time: 4.816441059112549
    83     1.892409e+03     3.456241e-02
 * time: 4.863898992538452
    84     1.892409e+03     3.114018e-02
 * time: 4.908630132675171
    85     1.892409e+03     4.080806e-02
 * time: 5.006025075912476
    86     1.892409e+03     6.722726e-02
 * time: 5.047101974487305
    87     1.892409e+03     1.006791e-01
 * time: 5.087952136993408
    88     1.892409e+03     1.303988e-01
 * time: 5.12930703163147
    89     1.892409e+03     1.228919e-01
 * time: 5.169641971588135
    90     1.892409e+03     6.433813e-02
 * time: 5.211593151092529
    91     1.892409e+03     1.314164e-02
 * time: 5.253737926483154
    92     1.892409e+03     4.929931e-04
 * time: 5.293926954269409
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
Likelihood Optimizer:                         BFGS
Dynamical system type:                 Closed form

Log-likelihood value:                    -1892.409
Number of subjects:                             30
Number of parameters:         Fixed      Optimized
                                  0             10
Observation records:         Active        Missing
    DV:                         540              0
    Total:                      540              0

-----------------------
             Estimate
-----------------------
tvvc          3.9757
tvq           0.042177
tvvp          3.6434
tvcl_hep      0.058572
tvcl_ren      0.1337
Ω₁,₁          0.18299
Ω₂,₂          0.081353
σ_add         0.032174
σ_prop        0.06101
dCRCL         1.1331
-----------------------

As before, our loglikelihood is higher (implying lower objective function value). Furthermore, our additive and combined error values remain almost the same while our Ω on CL, Ω₁,₁, decreased. This implies that the CRCL covariate with an estimated exponent dCRCL is definitely assisting in a better model fit.

Finally let’s include a separated CL model based on sex as a third covariate model:

covariate_model_wt_crcl_sex = @model begin
    @param begin
        tvvc ∈ RealDomain(; lower = 0)
        tvq ∈ RealDomain(; lower = 0)
        tvvp ∈ RealDomain(; lower = 0)
        tvcl_hep_M ∈ RealDomain(; lower = 0)
        tvcl_hep_F ∈ RealDomain(; lower = 0)
        tvcl_ren_M ∈ RealDomain(; lower = 0)
        tvcl_ren_F ∈ RealDomain(; lower = 0)
        Ω ∈ PDiagDomain(2)
        σ_add ∈ RealDomain(; lower = 0)
        σ_prop ∈ RealDomain(; lower = 0)
        dCRCL_M ∈ RealDomain()
        dCRCL_F ∈ RealDomain()
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates begin
        WT
        CRCL
        SEX
    end

    @pre begin
        tvcl_hep = ifelse(SEX == "M", tvcl_hep_M, tvcl_hep_F)
        tvcl_ren = ifelse(SEX == "M", tvcl_ren_M, tvcl_ren_F)
        dCRCL = ifelse(SEX == "M", dCRCL_M, dCRCL_F)
        hepCL = tvcl_hep * (WT / 70)^0.75
        renCL = tvcl_ren * (CRCL / 100)^dCRCL
        CL = (hepCL + renCL) * exp(η[1])
        Vc = tvvc * (WT / 70) * exp(η[2])
        Q = tvq
        Vp = tvvp
    end

    @dynamics Central1Periph1

    @derived begin
        cp := @. Central / Vc
        DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
    end
end
PumasModel
  Parameters: tvvc, tvq, tvvp, tvcl_hep_M, tvcl_hep_F, tvcl_ren_M, tvcl_ren_F, Ω, σ_add, σ_prop, dCRCL_M, dCRCL_F
  Random effects: η
  Covariates: WT, CRCL, SEX
  Dynamical system variables: Central, Peripheral
  Dynamical system type: Closed form
  Derived: DV
  Observed: DV

In the covariate_model_wt_crcl_sex model we are keeping our allometric scaling on WT, the CRCL effect on renCL, and the breakdown of CL into hepCL and renCL from before. However we are separating them with different values by sex. Hence, we have a new covariate SEX and six new parameters in the @param block by expanding tvcl_hep, tvcl_ren, and dCRCL into male (suffix M) and female (suffix F).

This is a good example on how to add create binary values based on covariate values such as SEX inside the @pre block with the ifelse function. Now, let’s fit this model. Note that we need a new initial parameters values’ list since the previous one we used had a single tvcl_hep, tvcl_ren, and dCRCL:

iparams_covariate_model_wt_crcl_sex = (;
    tvvc = 5,
    tvcl_hep_M = 0.01,
    tvcl_hep_F = 0.01,
    tvcl_ren_M = 0.01,
    tvcl_ren_F = 0.01,
    tvq = 0.01,
    tvvp = 10,
    Ω = Diagonal([0.01, 0.01]),
    σ_add = 0.1,
    σ_prop = 0.1,
    dCRCL_M = 0.9,
    dCRCL_F = 0.9,
)
(tvvc = 5,
 tvcl_hep_M = 0.01,
 tvcl_hep_F = 0.01,
 tvcl_ren_M = 0.01,
 tvcl_ren_F = 0.01,
 tvq = 0.01,
 tvvp = 10,
 Ω = [0.01 0.0; 0.0 0.01],
 σ_add = 0.1,
 σ_prop = 0.1,
 dCRCL_M = 0.9,
 dCRCL_F = 0.9,)
fit_covariate_model_wt_crcl_sex =
    fit(covariate_model_wt_crcl_sex, pop, iparams_covariate_model_wt_crcl_sex, FOCE())
[ 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     8.554152e+03     5.650279e+03
 * time: 7.295608520507812e-5
     1     3.641387e+03     1.432080e+03
 * time: 0.09355998039245605
     2     3.290450e+03     5.274782e+02
 * time: 0.15543007850646973
     3     3.185512e+03     2.173676e+02
 * time: 0.2118821144104004
     4     3.143009e+03     1.479653e+02
 * time: 0.2698330879211426
     5     3.128511e+03     8.980031e+01
 * time: 0.32782912254333496
     6     3.123188e+03     5.033125e+01
 * time: 0.39075517654418945
     7     3.120794e+03     4.279722e+01
 * time: 0.45375609397888184
     8     3.118627e+03     3.971051e+01
 * time: 0.5186901092529297
     9     3.115300e+03     8.456587e+01
 * time: 0.5837070941925049
    10     3.109353e+03     1.350354e+02
 * time: 2.1211180686950684
    11     3.095894e+03     1.998258e+02
 * time: 2.17923903465271
    12     2.988214e+03     4.366433e+02
 * time: 2.267376184463501
    13     2.896081e+03     5.505943e+02
 * time: 2.4788661003112793
    14     2.652467e+03     7.300323e+02
 * time: 3.8669991493225098
    15     2.560937e+03     6.973661e+02
 * time: 4.256755113601685
    16     2.254941e+03     2.740033e+02
 * time: 4.370934009552002
    17     2.222509e+03     2.034303e+02
 * time: 4.498098134994507
    18     2.171255e+03     2.449580e+02
 * time: 4.602251052856445
    19     2.024532e+03     1.121511e+02
 * time: 4.7019171714782715
    20     1.993723e+03     1.042814e+02
 * time: 4.792510986328125
    21     1.985113e+03     8.079014e+01
 * time: 4.888077974319458
    22     1.976757e+03     7.054196e+01
 * time: 4.974152088165283
    23     1.969970e+03     6.070322e+01
 * time: 5.072671175003052
    24     1.961095e+03     6.810782e+01
 * time: 5.159670114517212
    25     1.947983e+03     8.116920e+01
 * time: 5.2527430057525635
    26     1.930371e+03     8.530051e+01
 * time: 5.3479321002960205
    27     1.910209e+03     6.993170e+01
 * time: 5.442656993865967
    28     1.899107e+03     3.362640e+01
 * time: 5.543466091156006
    29     1.898022e+03     2.642220e+01
 * time: 5.611865997314453
    30     1.897055e+03     1.213144e+01
 * time: 5.680430173873901
    31     1.896596e+03     7.773239e+00
 * time: 5.749536037445068
    32     1.896538e+03     7.997039e+00
 * time: 5.815093994140625
    33     1.896451e+03     8.160909e+00
 * time: 5.879706144332886
    34     1.896283e+03     8.237721e+00
 * time: 5.966552019119263
    35     1.895903e+03     1.520219e+01
 * time: 6.053897142410278
    36     1.895272e+03     2.358916e+01
 * time: 6.12993597984314
    37     1.894536e+03     2.461296e+01
 * time: 6.218669176101685
    38     1.893995e+03     1.546128e+01
 * time: 6.306293964385986
    39     1.893858e+03     6.976137e+00
 * time: 6.398015022277832
    40     1.893833e+03     6.019466e+00
 * time: 6.471220016479492
    41     1.893786e+03     3.827201e+00
 * time: 6.542675971984863
    42     1.893714e+03     3.323412e+00
 * time: 6.610147953033447
    43     1.893592e+03     3.215150e+00
 * time: 6.656937122344971
    44     1.893435e+03     6.534965e+00
 * time: 6.703886032104492
    45     1.893286e+03     7.424154e+00
 * time: 6.766718149185181
    46     1.893190e+03     5.552627e+00
 * time: 6.814779043197632
    47     1.893139e+03     3.222316e+00
 * time: 6.861797094345093
    48     1.893120e+03     3.015339e+00
 * time: 6.907350063323975
    49     1.893107e+03     3.244809e+00
 * time: 6.951519966125488
    50     1.893080e+03     6.163100e+00
 * time: 6.997014045715332
    51     1.893027e+03     9.824713e+00
 * time: 7.058634042739868
    52     1.892912e+03     1.390100e+01
 * time: 7.10608696937561
    53     1.892734e+03     1.510937e+01
 * time: 7.153138160705566
    54     1.892561e+03     1.008563e+01
 * time: 7.199366092681885
    55     1.892485e+03     3.730668e+00
 * time: 7.245110034942627
    56     1.892471e+03     3.380261e+00
 * time: 7.2892069816589355
    57     1.892463e+03     3.167904e+00
 * time: 7.349027156829834
    58     1.892441e+03     4.152065e+00
 * time: 7.394813060760498
    59     1.892391e+03     7.355996e+00
 * time: 7.4400999546051025
    60     1.892268e+03     1.195397e+01
 * time: 7.486729145050049
    61     1.892026e+03     1.640783e+01
 * time: 7.532876968383789
    62     1.891735e+03     1.593576e+01
 * time: 7.595043182373047
    63     1.891569e+03     8.316423e+00
 * time: 7.643293142318726
    64     1.891494e+03     3.948212e+00
 * time: 7.6889731884002686
    65     1.891481e+03     3.911593e+00
 * time: 7.733180999755859
    66     1.891457e+03     3.875559e+00
 * time: 7.777536153793335
    67     1.891405e+03     3.811247e+00
 * time: 7.822023153305054
    68     1.891262e+03     3.657045e+00
 * time: 7.88151216506958
    69     1.890930e+03     4.957405e+00
 * time: 7.928992033004761
    70     1.890317e+03     6.657726e+00
 * time: 7.978309154510498
    71     1.889660e+03     6.086302e+00
 * time: 8.026012182235718
    72     1.889303e+03     2.270929e+00
 * time: 8.07204008102417
    73     1.889253e+03     7.695301e-01
 * time: 8.117138147354126
    74     1.889252e+03     7.382144e-01
 * time: 8.176329135894775
    75     1.889251e+03     7.187898e-01
 * time: 8.22130799293518
    76     1.889251e+03     7.215047e-01
 * time: 8.264449119567871
    77     1.889250e+03     7.235155e-01
 * time: 8.307260990142822
    78     1.889249e+03     7.246818e-01
 * time: 8.350273132324219
    79     1.889244e+03     7.257796e-01
 * time: 8.394158124923706
    80     1.889233e+03     7.198190e-01
 * time: 8.454494953155518
    81     1.889204e+03     1.089029e+00
 * time: 8.502465963363647
    82     1.889142e+03     1.801601e+00
 * time: 8.549183130264282
    83     1.889043e+03     2.967917e+00
 * time: 8.59532904624939
    84     1.888889e+03     2.965856e+00
 * time: 8.645277976989746
    85     1.888705e+03     5.933557e-01
 * time: 8.69391417503357
    86     1.888655e+03     9.577696e-01
 * time: 8.775105953216553
    87     1.888582e+03     1.498494e+00
 * time: 8.821997165679932
    88     1.888533e+03     1.502753e+00
 * time: 8.868581056594849
    89     1.888490e+03     1.184664e+00
 * time: 8.916514158248901
    90     1.888480e+03     6.684517e-01
 * time: 8.963972091674805
    91     1.888476e+03     3.680034e-01
 * time: 9.02361512184143
    92     1.888476e+03     4.720040e-01
 * time: 9.084444999694824
    93     1.888476e+03     4.768646e-01
 * time: 9.138576984405518
    94     1.888475e+03     4.736673e-01
 * time: 9.18110203742981
    95     1.888475e+03     4.552764e-01
 * time: 9.223320007324219
    96     1.888474e+03     5.193733e-01
 * time: 9.265405178070068
    97     1.888473e+03     8.850112e-01
 * time: 9.3227219581604
    98     1.888468e+03     1.461600e+00
 * time: 9.367063999176025
    99     1.888458e+03     2.209127e+00
 * time: 9.41123104095459
   100     1.888437e+03     2.961239e+00
 * time: 9.455619096755981
   101     1.888407e+03     2.978463e+00
 * time: 9.501255989074707
   102     1.888384e+03     1.707201e+00
 * time: 9.545618057250977
   103     1.888381e+03     6.199354e-01
 * time: 9.605647087097168
   104     1.888380e+03     5.170909e-01
 * time: 9.651757001876831
   105     1.888378e+03     1.037408e-01
 * time: 9.69723105430603
   106     1.888378e+03     8.473247e-02
 * time: 9.738769054412842
   107     1.888378e+03     8.364978e-02
 * time: 9.78000807762146
   108     1.888378e+03     8.080446e-02
 * time: 9.821596145629883
   109     1.888378e+03     7.873905e-02
 * time: 9.87653398513794
   110     1.888378e+03     7.798398e-02
 * time: 9.91849398612976
   111     1.888378e+03     7.788170e-02
 * time: 9.9601731300354
   112     1.888378e+03     7.776461e-02
 * time: 10.000357151031494
   113     1.888378e+03     9.023662e-02
 * time: 10.040759086608887
   114     1.888378e+03     1.631390e-01
 * time: 10.081798076629639
   115     1.888378e+03     2.768731e-01
 * time: 10.138453960418701
   116     1.888377e+03     4.462386e-01
 * time: 10.182181119918823
   117     1.888377e+03     6.643297e-01
 * time: 10.226317167282104
   118     1.888375e+03     8.433397e-01
 * time: 10.269840955734253
   119     1.888374e+03     7.596814e-01
 * time: 10.31320309638977
   120     1.888373e+03     3.638119e-01
 * time: 10.356076002120972
   121     1.888372e+03     8.306034e-02
 * time: 10.41214895248413
   122     1.888372e+03     2.084513e-02
 * time: 10.455090999603271
   123     1.888372e+03     2.056419e-02
 * time: 10.498250961303711
   124     1.888372e+03     2.044080e-02
 * time: 10.537997961044312
   125     1.888372e+03     2.035196e-02
 * time: 10.5762460231781
   126     1.888372e+03     2.021264e-02
 * time: 10.614962100982666
   127     1.888372e+03     1.998163e-02
 * time: 10.65455412864685
   128     1.888372e+03     3.161638e-02
 * time: 10.710819005966187
   129     1.888372e+03     5.509218e-02
 * time: 10.75266408920288
   130     1.888372e+03     9.275848e-02
 * time: 10.793850183486938
   131     1.888372e+03     1.528749e-01
 * time: 10.834737062454224
   132     1.888372e+03     2.461766e-01
 * time: 10.876012086868286
   133     1.888372e+03     3.799362e-01
 * time: 10.917564153671265
   134     1.888371e+03     5.311665e-01
 * time: 10.976363182067871
   135     1.888369e+03     6.019039e-01
 * time: 11.020753145217896
   136     1.888367e+03     4.664896e-01
 * time: 11.064316034317017
   137     1.888366e+03     1.404934e-01
 * time: 11.10679006576538
   138     1.888365e+03     8.513331e-02
 * time: 11.149557113647461
   139     1.888364e+03     1.244077e-01
 * time: 11.191015005111694
   140     1.888364e+03     1.028085e-01
 * time: 11.248507022857666
   141     1.888364e+03     5.162231e-02
 * time: 11.292840957641602
   142     1.888364e+03     5.149630e-02
 * time: 11.335594177246094
   143     1.888364e+03     3.147284e-02
 * time: 11.377110958099365
   144     1.888364e+03     2.104766e-02
 * time: 11.417720079421997
   145     1.888364e+03     6.539151e-03
 * time: 11.45687198638916
   146     1.888364e+03     2.540196e-03
 * time: 11.51254916191101
   147     1.888364e+03     4.362661e-03
 * time: 11.553251028060913
   148     1.888364e+03     3.034416e-03
 * time: 11.59301209449768
   149     1.888364e+03     5.953892e-04
 * time: 11.631872177124023
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
Likelihood Optimizer:                         BFGS
Dynamical system type:                 Closed form

Log-likelihood value:                   -1888.3638
Number of subjects:                             30
Number of parameters:         Fixed      Optimized
                                  0             13
Observation records:         Active        Missing
    DV:                         540              0
    Total:                      540              0

--------------------------
               Estimate
--------------------------
tvvc            3.976
tvq             0.04239
tvvp            3.7249
tvcl_hep_M      1.7175e-7
tvcl_hep_F      0.13348
tvcl_ren_M      0.19378
tvcl_ren_F      0.042211
Ω₁,₁            0.14046
Ω₂,₂            0.081349
σ_add           0.032171
σ_prop          0.061007
dCRCL_M         0.94821
dCRCL_F         1.9405
--------------------------

As before, our loglikelihood is higher (implying lower objective function value). This is expected since we also added six new parameters to the model.

5 Step 4 - Model Comparison

Now that we’ve fitted all of our models we need to compare them and choose one for our final model.

We begin by analyzing the model metrics. This can be done with the metrics_table function:

metrics_table(fit_base_model)
WARNING: using deprecated binding Distributions.MatrixReshaped in Pumas.
, use Distributions.ReshapedDistribution{2, S, D} where D<:Distributions.Distribution{Distributions.ArrayLikeVariate{1}, S} where S<:Distributions.ValueSupport instead.
17×2 DataFrame
Row Metric Value
String Any
1 Successful true
2 Estimation Time 3.794
3 Subjects 30
4 Fixed Parameters 0
5 Optimized Parameters 8
6 DV Active Observations 540
7 DV Missing Observations 0
8 Total Active Observations 540
9 Total Missing Observations 0
10 Likelihood Approximation Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}
11 LogLikelihood (LL) -1901.82
12 -2LL 3803.63
13 AIC 3819.63
14 BIC 3853.96
15 (η-shrinkage) η₁ -0.015
16 (η-shrinkage) η₂ -0.013
17 (ϵ-shrinkage) DV 0.056 
metrics_table(fit_covariate_model_wt)
17×2 DataFrame
Row Metric Value
String Any
1 Successful true
2 Estimation Time 2.52
3 Subjects 30
4 Fixed Parameters 0
5 Optimized Parameters 8
6 DV Active Observations 540
7 DV Missing Observations 0
8 Total Active Observations 540
9 Total Missing Observations 0
10 Likelihood Approximation Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}
11 LogLikelihood (LL) -1896.46
12 -2LL 3792.93
13 AIC 3808.93
14 BIC 3843.26
15 (η-shrinkage) η₁ -0.014
16 (η-shrinkage) η₂ -0.012
17 (ϵ-shrinkage) DV 0.056 
metrics_table(fit_covariate_model_wt_crcl)
17×2 DataFrame
Row Metric Value
String Any
1 Successful true
2 Estimation Time 5.294
3 Subjects 30
4 Fixed Parameters 0
5 Optimized Parameters 10
6 DV Active Observations 540
7 DV Missing Observations 0
8 Total Active Observations 540
9 Total Missing Observations 0
10 Likelihood Approximation Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}
11 LogLikelihood (LL) -1892.41
12 -2LL 3784.82
13 AIC 3804.82
14 BIC 3847.73
15 (η-shrinkage) η₁ -0.014
16 (η-shrinkage) η₂ -0.012
17 (ϵ-shrinkage) DV 0.056 
metrics_table(fit_covariate_model_wt_crcl_sex)
17×2 DataFrame
Row Metric Value
String Any
1 Successful true
2 Estimation Time 11.632
3 Subjects 30
4 Fixed Parameters 0
5 Optimized Parameters 13
6 DV Active Observations 540
7 DV Missing Observations 0
8 Total Active Observations 540
9 Total Missing Observations 0
10 Likelihood Approximation Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}
11 LogLikelihood (LL) -1888.36
12 -2LL 3776.73
13 AIC 3802.73
14 BIC 3858.52
15 (η-shrinkage) η₁ -0.013
16 (η-shrinkage) η₂ -0.012
17 (ϵ-shrinkage) DV 0.056

metrics_table outputs all of the model metrics we might be interested with respect to a certain model. That includes metadata such as estimation time, number of subjects, how many parameters were optimized and fixed, and number of observations. It also includes common model metrics like AIC, BIC, objective function value with constant (-2 loglikelihood), and so on.

We can also do an innerjoin (check our Data Wrangling Tutorials) to get all metrics into a single DataFrame:

all_metrics = innerjoin(
    metrics_table(fit_base_model),
    metrics_table(fit_covariate_model_wt),
    metrics_table(fit_covariate_model_wt_crcl),
    metrics_table(fit_covariate_model_wt_crcl_sex);
    on = :Metric,
    makeunique = true,
);
rename!(
    all_metrics,
    :Value => :Base_Model,
    :Value_1 => :Covariate_Model_WT,
    :Value_2 => :Covariate_Model_WT_CRCL,
    :Value_3 => :Covariate_Model_WT_CRCL_SEX,
)
17×5 DataFrame
Row Metric Base_Model Covariate_Model_WT Covariate_Model_WT_CRCL Covariate_Model_WT_CRCL_SEX
String Any Any Any Any
1 Successful true true true true
2 Estimation Time 3.794 2.52 5.294 11.632
3 Subjects 30 30 30 30
4 Fixed Parameters 0 0 0 0
5 Optimized Parameters 8 8 10 13
6 DV Active Observations 540 540 540 540
7 DV Missing Observations 0 0 0 0
8 Total Active Observations 540 540 540 540
9 Total Missing Observations 0 0 0 0
10 Likelihood Approximation Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}} Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}} Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}} Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}
11 LogLikelihood (LL) -1901.82 -1896.46 -1892.41 -1888.36
12 -2LL 3803.63 3792.93 3784.82 3776.73
13 AIC 3819.63 3808.93 3804.82 3802.73
14 BIC 3853.96 3843.26 3847.73 3858.52
15 (η-shrinkage) η₁ -0.015 -0.014 -0.014 -0.013
16 (η-shrinkage) η₂ -0.013 -0.012 -0.012 -0.012
17 (ϵ-shrinkage) DV 0.056 0.056 0.056 0.056

We can also use specific functions to retrieve those. For example, in order to get a model’s AIC you can use the aic function:

aic(fit_base_model)
3819.629984952819
aic(fit_covariate_model_wt)
3808.9264607805967
aic(fit_covariate_model_wt_crcl)
3804.817947371705
aic(fit_covariate_model_wt_crcl_sex)
3802.7275243741956

We should favor lower values of AIC, hence, the covariate model with weight, creatinine clerance, and different sex effects on clearance should be preferred, i.e. covariate_model_wt_crcl_sex.

5.1 Goodness of Fit Plots

Additionally, we should inspect the goodness of fit of the model. This is done with the plotting function goodness_of_fit, which should be given a result from a inspect function. So, let’s first call inspect on our covariate_model_wt_crcl_sex candidate for best model:

inspect_covariate_model_wt_crcl_sex = inspect(fit_covariate_model_wt_crcl_sex)
[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating individual parameters.
[ Info: Evaluating dose control parameters.
[ Info: Done.
FittedPumasModelInspection

Fitting was successful: true
Likehood approximation used for weighted residuals: FOCE

And now we pass inspect_covariate_model_wt_crcl_sex to the goodness_of_fit plotting function:

goodness_of_fit(inspect_covariate_model_wt_crcl_sex)

The idea is that the population predictions (preds) capture the general tendency of the observations while the individual predictions (ipreds) should coincide much more closely with the observations. That is exactly what we are observing in the top row subplots in our goodness of fit plot.

Regarding the bottom row, on the left we have the weighted population residuals (wres) against time, and on the right we have the weighted individual residuals (iwres) against ipreds. Here we should not see any perceived pattern, which indicates that the error model in the model has a mean 0 and constant variance. Like before, this seems to be what we are observing in our goodness of fit plot.

Hence, our covariate model with allometric scaling and different sex creatinine clearance effectw on clearance is a good candidate for our final model.

6 Time-Varying Covariates

Pumas can handle time-varying covariates. This happens automatically if, when parsing a dataset, read_pumas detects that covariate values change over time.

6.1 painord Dataset

Here’s an example with an ordinal regression using the painord dataset from PharmaDatasets.jl. :painord is our observations measuring the perceived pain in a scale from 0 to 3, which we need to have the range shifted by 1 (1 to 4). Additionally, we’ll use the concentration in plasma, :conc as a covariate. Of course, :conc varies with time, thus, it is a time-varying covariate:

painord = dataset("pumas/pain_remed")
first(painord, 5)
5×8 DataFrame
Row id arm dose time conc painord dv remed
Int64 Int64 Int64 Float64 Float64 Int64 Int64 Int64
1 1 2 20 0.0 0.0 3 0 0
2 1 2 20 0.5 1.15578 1 1 0
3 1 2 20 1.0 1.37211 0 1 0
4 1 2 20 1.5 1.30058 0 1 0
5 1 2 20 2.0 1.19195 1 1 0 
@rtransform! painord :painord = :painord + 1;
describe(painord, :mean, :std, :first, :last, :eltype)
8×6 DataFrame
Row variable mean std first last eltype
Symbol Float64 Float64 Real Real DataType
1 id 80.5 46.1992 1 160 Int64
2 arm 1.5 1.11833 2 0 Int64
3 dose 26.25 31.9017 20 0 Int64
4 time 3.375 2.5183 0.0 8.0 Float64
5 conc 0.93018 1.49902 0.0 0.0 Float64
6 painord 2.50208 0.863839 4 4 Int64
7 dv 0.508333 0.500061 0 0 Int64
8 remed 0.059375 0.236387 0 0 Int64 
unique(painord.dose)
4-element Vector{Int64}:
 20
 80
  0
  5

As we can see we have 160 subjects were given either 0, 5, 20, or 80 units of a certain painkiller drug.

:conc is the drug concentration in plasma and :painord is the perceived pain in a scale from 1 to 4.

First, we’ll parse the painord dataset into a Population. Note that we’ll be using event_data=false since we do not have any dosing rows:

pop_ord =
    read_pumas(painord; observations = [:painord], covariates = [:conc], event_data = false)
Note

We won’t be going into the details of the ordinal regression model in this tutorial. We highly encourage you to take a look at the Ordinal Regression Pumas Tutorial for an in-depth explanation.

We’ll build an ordinal regression model declaring :conc as a covariate. In the @derived block we’ll state the the likelihood of :painord follows a Categorical distribution:

ordinal_model = @model begin
    @param begin
        b₁ ∈ RealDomain(; init = 0)
        b₂ ∈ RealDomain(; init = 1)
        b₃ ∈ RealDomain(; init = 1)
        slope ∈ RealDomain(; init = 0)
    end

    @covariates conc # time-varying

    @pre begin
        effect = slope * conc

        # Logit of cumulative probabilities
        lge₁ = b₁ + effect
        lge₂ = lge₁ - b₂
        lge₃ = lge₂ - b₃

        # Probabilities of >=1 and >=2 and >=3
        pge₁ = logistic(lge₁)
        pge₂ = logistic(lge₂)
        pge₃ = logistic(lge₃)

        # Probabilities of Y=1,2,3,4
        p₁ = 1.0 - pge₁
        p₂ = pge₁ - pge₂
        p₃ = pge₂ - pge₃
        p₄ = pge₃
    end

    @derived begin
        painord ~ @. Categorical(p₁, p₂, p₃, p₄)
    end
end
PumasModel
  Parameters: b₁, b₂, b₃, slope
  Random effects:
  Covariates: conc
  Dynamical system variables:
  Dynamical system type: No dynamical model
  Derived: painord
  Observed: painord

Finally we’ll fit our model using NaivePooled estimation method since it does not have any random-effects, i.e. no @random block:

ordinal_fit = fit(ordinal_model, pop_ord, init_params(ordinal_model), NaivePooled())
[ 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     3.103008e+03     7.031210e+02
 * time: 6.818771362304688e-5
     1     2.994747e+03     1.083462e+03
 * time: 0.005970001220703125
     2     2.406265e+03     1.884408e+02
 * time: 0.011843204498291016
     3     2.344175e+03     7.741610e+01
 * time: 0.01760101318359375
     4     2.323153e+03     2.907642e+01
 * time: 0.023398160934448242
     5     2.318222e+03     2.273295e+01
 * time: 0.029120206832885742
     6     2.316833e+03     1.390527e+01
 * time: 0.03482699394226074
     7     2.316425e+03     4.490883e+00
 * time: 0.04061007499694824
     8     2.316362e+03     9.374519e-01
 * time: 0.046227216720581055
     9     2.316356e+03     1.928785e-01
 * time: 0.05182600021362305
    10     2.316355e+03     3.119615e-02
 * time: 0.05752921104431152
    11     2.316355e+03     6.215513e-03
 * time: 0.06323504447937012
    12     2.316355e+03     8.313206e-04
 * time: 0.06901812553405762
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:              NaivePooled
Likelihood Optimizer:                         BFGS
Dynamical system type:          No dynamical model

Log-likelihood value:                   -2316.3554
Number of subjects:                            160
Number of parameters:         Fixed      Optimized
                                  0              4
Observation records:         Active        Missing
    painord:                   1920              0
    Total:                     1920              0

-------------------
          Estimate
-------------------
b₁         2.5112
b₂         2.1951
b₃         1.9643
slope     -0.38871
-------------------

As expected, the ordinal model fit estimates a negative effect of :conc on :painord measured by the slope parameter.

7 Missing Data in Covariates

The way how Pumas handles missing values inside covariates depends if the covariate is constant or time-varying. For both cases Pumas will interpolate the available values to fill in the missing values. If, for any subject, all of the covariate’s values are missing, Pumas will throw an error while parsing the data with read_pumas.

For both missing constant and time-varying covariates, Pumas, by default, does piece-wise constant interpolation with “next observation carried backward” (NOCB, NONMEM default). Of course for constant covariates the interpolated values over the missing values will be constant values. This can be adjusted with the covariates_direction keyword argument of read_pumas. The default value :right is NOCB and :left is “last observation carried forward” (LOCF, Monolix default).

Hence, for LOCF, you can use the following:

pop = read_pumas(pkdata; covariates_direction = :left)

along with any other required keyword arguments for column mapping.

Note

The same behavior for covariates_direction applies to time-varying covariates during the interpolation in the ODE solver. They will also be piece-wise constant interpolated following either NOCB or LOCF depending on the covariates_direction value.

8 Categorical Covariates

In some situations, you’ll find yourself with categorical covariates with multiple levels, instead of binary or continuous covariates. Categorical covariates are covariates that can take on a finite number of distinct values.

Pumas can easily address categorical covariates. In the @pre block you can use a nested if ... elseif ... else statement to handle the different categories.

For example:

@pre begin
    CL = if RACE == 1
        tvcl * (WT / 70)^dwtdcl * exp(η[1]) * drace1dcl
    elseif RACE == 2
        tvcl * (WT / 70)^dwtdcl * exp(η[1]) * drace2dcl
    elseif RACE == 3
        tvcl * (WT / 70)^dwtdcl * exp(η[1]) * drace3dcl
    end
end

Here we are conditioning the clearance (CL) on the RACE covariate by modulating which population-level parameter will be used for the clearance calculation: drace1dcl, drace2dcl, and drace3dcl.

There’s nothing wrong with the code above, but it can be a bit cumbersome to write and read. In order to make it more readable and maintainable, you can use the following example:

@pre begin
    raceoncl = race1cl^(race == 1) * race2cl^(race == 2) * race3cl^(race == 3)
    CL = tvcl * raceoncl
end

Here we are using the ^ operator to raise each race value to the power of the race1cl, race2cl, and race3cl values. If any of the race values is not equal to the race value, the result will be 1, otherwise it will be the respective race1cl, race2cl, or race3cl value.

9 Conclusion

This tutorial shows how to build covariate model in Pumas in a workflow approach. The main purpose was to inform how to:

  • parse covariate data into a Population
  • add covariate information into a model

We also went over what are the differences between constant and time-varying covariates and how does Pumas deal with missing data inside covariates.