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 ~ @. CombinedNormal(cp, σ_add, σ_prop)
    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, Ω = Diagonal([0.01, 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.03862500190734863
     1     3.110315e+03     9.706222e+02
 * time: 7.614857912063599
     2     2.831659e+03     7.817006e+02
 * time: 8.025831937789917
     3     2.405281e+03     2.923696e+02
 * time: 8.04499888420105
     4     2.370406e+03     3.032286e+02
 * time: 8.058904886245728
     5     2.313631e+03     3.126188e+02
 * time: 8.07377290725708
     6     2.263986e+03     2.946697e+02
 * time: 8.08785891532898
     7     2.201497e+03     2.519434e+02
 * time: 8.110696077346802
     8     2.151491e+03     1.868334e+02
 * time: 8.131824016571045
     9     2.104228e+03     1.156870e+02
 * time: 8.153064012527466
    10     2.086667e+03     8.564363e+01
 * time: 8.166585922241211
    11     2.077820e+03     8.122079e+01
 * time: 8.179959058761597
    12     2.073275e+03     7.699389e+01
 * time: 8.192631006240845
    13     2.066038e+03     6.982578e+01
 * time: 8.205197095870972
    14     2.051077e+03     9.768714e+01
 * time: 8.217514038085938
    15     2.026074e+03     1.119946e+02
 * time: 8.229926109313965
    16     1.999606e+03     8.197429e+01
 * time: 8.242930889129639
    17     1.982606e+03     5.004245e+01
 * time: 8.256038904190063
    18     1.973959e+03     6.738765e+01
 * time: 8.270637035369873
    19     1.969329e+03     7.915826e+01
 * time: 8.283421993255615
    20     1.963614e+03     8.304473e+01
 * time: 8.296058893203735
    21     1.952491e+03     5.562904e+01
 * time: 8.308881998062134
    22     1.934222e+03     3.768831e+01
 * time: 8.321918964385986
    23     1.926714e+03     5.596542e+01
 * time: 9.217858076095581
    24     1.923574e+03     4.898901e+01
 * time: 9.227950096130371
    25     1.914691e+03     1.764745e+01
 * time: 9.238723039627075
    26     1.907315e+03     1.300123e+01
 * time: 9.250267028808594
    27     1.905803e+03     1.502353e+01
 * time: 9.260088920593262
    28     1.904848e+03     1.092933e+01
 * time: 9.271151065826416
    29     1.904765e+03     1.121176e+01
 * time: 9.280560970306396
    30     1.904686e+03     1.122618e+01
 * time: 9.290011882781982
    31     1.904360e+03     1.139304e+01
 * time: 9.299659967422485
    32     1.903849e+03     1.540399e+01
 * time: 9.309117078781128
    33     1.903177e+03     1.433745e+01
 * time: 9.318711042404175
    34     1.902702e+03     1.076929e+01
 * time: 9.327713966369629
    35     1.902515e+03     1.064215e+01
 * time: 9.336878061294556
    36     1.902477e+03     9.724041e+00
 * time: 9.345686912536621
    37     1.902363e+03     6.745388e+00
 * time: 9.354907989501953
    38     1.902219e+03     6.082204e+00
 * time: 9.36432409286499
    39     1.902109e+03     3.617708e+00
 * time: 9.373845100402832
    40     1.902077e+03     1.440114e+00
 * time: 9.382976055145264
    41     1.902070e+03     1.431248e+00
 * time: 9.391899108886719
    42     1.902068e+03     1.484310e+00
 * time: 9.400404930114746
    43     1.902065e+03     1.501982e+00
 * time: 9.409008026123047
    44     1.902050e+03     1.474910e+00
 * time: 9.417813062667847
    45     1.902022e+03     1.928989e+00
 * time: 9.426722049713135
    46     1.901975e+03     2.573574e+00
 * time: 9.43538498878479
    47     1.901935e+03     1.843132e+00
 * time: 9.443902015686035
    48     1.901921e+03     5.848792e-01
 * time: 9.452617883682251
    49     1.901919e+03     5.505571e-01
 * time: 9.461019039154053
    50     1.901919e+03     5.543746e-01
 * time: 9.469356060028076
    51     1.901919e+03     5.625072e-01
 * time: 9.477771997451782
    52     1.901918e+03     5.753028e-01
 * time: 9.487117052078247
    53     1.901915e+03     5.920985e-01
 * time: 9.49755311012268
    54     1.901909e+03     9.386505e-01
 * time: 9.508377075195312
    55     1.901896e+03     1.456476e+00
 * time: 9.546391010284424
    56     1.901870e+03     1.923409e+00
 * time: 9.556219100952148
    57     1.901840e+03     1.799523e+00
 * time: 9.565995931625366
    58     1.901821e+03     8.979705e-01
 * time: 9.575829029083252
    59     1.901815e+03     1.548766e-01
 * time: 9.585295915603638
    60     1.901815e+03     4.957007e-02
 * time: 9.594965934753418
    61     1.901815e+03     7.184157e-03
 * time: 9.604183912277222
    62     1.901815e+03     6.778480e-04
 * time: 9.612416982650757

FittedPumasModel

Dynamical system type:                 Closed form

Number of subjects:                             30

Observation records:         Active        Missing
    DV:                         540              0
    Total:                      540              0

Number of parameters:      Constant      Optimized
                                  0              8

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                    -1901.815

------------------
         Estimate
------------------
tvcl     0.1542
tvvc     4.5856
tvq      0.042341
tvvp     3.7085
Ω₁,₁     0.26468
Ω₂,₂     0.10627
σ_add    0.032184
σ_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 ~ @. CombinedNormal(cp, σ_add, σ_prop)
    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: 2.09808349609375e-5
     1     2.846050e+03     1.128657e+03
 * time: 1.2706060409545898
     2     2.472982e+03     7.008264e+02
 * time: 1.2836251258850098
     3     2.180690e+03     2.312704e+02
 * time: 1.2963290214538574
     4     2.125801e+03     1.862929e+02
 * time: 1.3073680400848389
     5     2.103173e+03     1.320946e+02
 * time: 1.3179879188537598
     6     2.091136e+03     1.103035e+02
 * time: 1.3282830715179443
     7     2.081443e+03     1.091133e+02
 * time: 1.3378539085388184
     8     2.071793e+03     7.197675e+01
 * time: 1.347717046737671
     9     2.062706e+03     7.623310e+01
 * time: 1.35782790184021
    10     2.057515e+03     6.885476e+01
 * time: 1.3676199913024902
    11     2.051133e+03     6.368504e+01
 * time: 1.3771350383758545
    12     2.038626e+03     7.730243e+01
 * time: 1.3866279125213623
    13     2.019352e+03     1.136864e+02
 * time: 1.3962500095367432
    14     1.997136e+03     1.005748e+02
 * time: 1.4060399532318115
    15     1.983023e+03     6.831478e+01
 * time: 1.4161460399627686
    16     1.977700e+03     5.649783e+01
 * time: 1.4258348941802979
    17     1.974583e+03     6.357642e+01
 * time: 1.4353430271148682
    18     1.967292e+03     7.658974e+01
 * time: 1.445444107055664
    19     1.951045e+03     6.130573e+01
 * time: 1.4569919109344482
    20     1.935868e+03     4.845839e+01
 * time: 1.5835881233215332
    21     1.929356e+03     6.325782e+01
 * time: 1.5946400165557861
    22     1.925187e+03     3.142245e+01
 * time: 1.60475492477417
    23     1.923733e+03     4.623400e+01
 * time: 1.6147890090942383
    24     1.918498e+03     5.347738e+01
 * time: 1.6248259544372559
    25     1.912383e+03     5.849125e+01
 * time: 1.6359190940856934
    26     1.905510e+03     3.254038e+01
 * time: 1.6467959880828857
    27     1.903629e+03     2.905618e+01
 * time: 1.6564741134643555
    28     1.902833e+03     2.907696e+01
 * time: 1.6661810874938965
    29     1.902447e+03     2.746037e+01
 * time: 1.675076961517334
    30     1.899399e+03     1.930949e+01
 * time: 1.6849210262298584
    31     1.898705e+03     1.186361e+01
 * time: 1.6944739818572998
    32     1.898505e+03     1.050402e+01
 * time: 1.7038259506225586
    33     1.898446e+03     1.033938e+01
 * time: 1.7173280715942383
    34     1.897992e+03     1.241146e+01
 * time: 1.7261180877685547
    35     1.897462e+03     1.753723e+01
 * time: 1.7350420951843262
    36     1.896942e+03     1.449253e+01
 * time: 1.744157075881958
    37     1.896739e+03     5.424819e+00
 * time: 1.7532179355621338
    38     1.896705e+03     1.181392e+00
 * time: 1.7615330219268799
    39     1.896704e+03     1.217604e+00
 * time: 1.7700340747833252
    40     1.896703e+03     1.226705e+00
 * time: 1.7787559032440186
    41     1.896698e+03     1.236962e+00
 * time: 1.7872989177703857
    42     1.896687e+03     1.534546e+00
 * time: 1.84529709815979
    43     1.896664e+03     2.656196e+00
 * time: 1.8543689250946045
    44     1.896629e+03     3.132023e+00
 * time: 1.8634419441223145
    45     1.896597e+03     2.066098e+00
 * time: 1.872675895690918
    46     1.896586e+03     6.603511e-01
 * time: 1.881716012954712
    47     1.896585e+03     6.269104e-01
 * time: 1.8904991149902344
    48     1.896585e+03     5.980951e-01
 * time: 1.8989429473876953
    49     1.896584e+03     5.903886e-01
 * time: 1.9073660373687744
    50     1.896583e+03     6.188684e-01
 * time: 1.9158189296722412
    51     1.896578e+03     8.758629e-01
 * time: 1.924375057220459
    52     1.896568e+03     1.384124e+00
 * time: 1.9330179691314697
    53     1.896546e+03     1.986298e+00
 * time: 1.9415850639343262
    54     1.896510e+03     2.277187e+00
 * time: 1.9502530097961426
    55     1.896480e+03     1.614779e+00
 * time: 1.9588029384613037
    56     1.896465e+03     4.382984e-01
 * time: 1.9674720764160156
    57     1.896463e+03     1.014017e-01
 * time: 1.9759440422058105
    58     1.896463e+03     1.390678e-02
 * time: 1.9839179515838623
    59     1.896463e+03     2.102101e-03
 * time: 1.9913511276245117
    60     1.896463e+03     2.971897e-04
 * time: 1.998690128326416

FittedPumasModel

Dynamical system type:                 Closed form

Number of subjects:                             30

Observation records:         Active        Missing
    DV:                         540              0
    Total:                      540              0

Number of parameters:      Constant      Optimized
                                  0              8

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -1896.4632

------------------
         Estimate
------------------
tvcl     0.13915
tvvc     3.9754
tvq      0.041988
tvvp     3.5723
Ω₁,₁     0.23874
Ω₂,₂     0.081372
σ_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 ~ @. CombinedNormal(cp, σ_add, σ_prop)
    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, Ω = Diagonal([0.01, 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: 1.5020370483398438e-5
     1     3.572050e+03     1.302046e+03
 * time: 1.3685929775238037
     2     3.266947e+03     5.384244e+02
 * time: 1.3836669921875
     3     3.150635e+03     1.918079e+02
 * time: 1.396766185760498
     4     3.108122e+03     1.277799e+02
 * time: 1.4088850021362305
     5     3.091358e+03     8.838080e+01
 * time: 1.4203600883483887
     6     3.082997e+03     8.321689e+01
 * time: 1.4319860935211182
     7     3.076379e+03     8.167702e+01
 * time: 1.4436659812927246
     8     3.067422e+03     1.177822e+02
 * time: 1.4557640552520752
     9     3.048580e+03     2.526969e+02
 * time: 1.4694221019744873
    10     2.993096e+03     6.325396e+02
 * time: 1.4906821250915527
    11     2.971933e+03     7.263683e+02
 * time: 1.640273094177246
    12     2.913647e+03     8.338659e+02
 * time: 1.6823010444641113
    13     2.815755e+03     8.865146e+02
 * time: 1.7232441902160645
    14     2.675887e+03     7.849516e+02
 * time: 1.7576169967651367
    15     2.586444e+03     6.412063e+02
 * time: 1.7870020866394043
    16     2.407387e+03     6.174235e+02
 * time: 1.8126351833343506
    17     2.255853e+03     4.700647e+02
 * time: 1.8256750106811523
    18     2.187754e+03     2.525045e+02
 * time: 1.8387041091918945
    19     2.166014e+03     2.482732e+02
 * time: 1.9012269973754883
    20     2.091740e+03     2.102833e+02
 * time: 1.9135210514068604
    21     2.069405e+03     9.995015e+01
 * time: 1.9250860214233398
    22     2.061725e+03     9.755106e+01
 * time: 1.9360551834106445
    23     2.050858e+03     9.071420e+01
 * time: 1.9472072124481201
    24     2.038959e+03     1.248437e+02
 * time: 1.9581170082092285
    25     1.981945e+03     1.834009e+02
 * time: 1.9696919918060303
    26     1.960293e+03     1.861977e+02
 * time: 1.98746919631958
    27     1.937145e+03     1.143323e+02
 * time: 1.9985320568084717
    28     1.919599e+03     6.766229e+01
 * time: 2.00924015045166
    29     1.908338e+03     4.657224e+01
 * time: 2.0199360847473145
    30     1.902080e+03     3.527271e+01
 * time: 2.0306200981140137
    31     1.899564e+03     2.865779e+01
 * time: 2.0412189960479736
    32     1.898724e+03     2.563582e+01
 * time: 2.0515801906585693
    33     1.898161e+03     2.299606e+01
 * time: 2.083052158355713
    34     1.897586e+03     1.978071e+01
 * time: 2.093820095062256
    35     1.897090e+03     1.672129e+01
 * time: 2.1046500205993652
    36     1.896519e+03     1.866527e+01
 * time: 2.1156070232391357
    37     1.895771e+03     2.155412e+01
 * time: 2.126420021057129
    38     1.895040e+03     1.767141e+01
 * time: 2.137453079223633
    39     1.894666e+03     9.006873e+00
 * time: 2.148207187652588
    40     1.894576e+03     2.803611e+00
 * time: 2.158508062362671
    41     1.894559e+03     2.603244e+00
 * time: 2.168354034423828
    42     1.894553e+03     2.569313e+00
 * time: 2.178161144256592
    43     1.894545e+03     2.503113e+00
 * time: 2.187995195388794
    44     1.894527e+03     2.349539e+00
 * time: 2.197780132293701
    45     1.894484e+03     2.239653e+00
 * time: 2.2258059978485107
    46     1.894391e+03     4.277324e+00
 * time: 2.2360880374908447
    47     1.894224e+03     6.085670e+00
 * time: 2.2475881576538086
    48     1.894026e+03     5.656832e+00
 * time: 2.258169174194336
    49     1.893915e+03     3.551925e+00
 * time: 2.268712043762207
    50     1.893889e+03     3.696450e+00
 * time: 2.278994083404541
    51     1.893881e+03     3.565807e+00
 * time: 2.2889270782470703
    52     1.893873e+03     3.402550e+00
 * time: 2.298718214035034
    53     1.893849e+03     3.025868e+00
 * time: 2.3086371421813965
    54     1.893796e+03     3.616269e+00
 * time: 2.3186941146850586
    55     1.893678e+03     5.277021e+00
 * time: 2.358971118927002
    56     1.893500e+03     5.738210e+00
 * time: 2.3695151805877686
    57     1.893355e+03     3.514473e+00
 * time: 2.380251169204712
    58     1.893299e+03     2.176745e+00
 * time: 2.3910820484161377
    59     1.893289e+03     2.230959e+00
 * time: 2.4013171195983887
    60     1.893285e+03     2.282110e+00
 * time: 2.411375045776367
    61     1.893269e+03     2.395111e+00
 * time: 2.4214820861816406
    62     1.893235e+03     2.508353e+00
 * time: 2.4314889907836914
    63     1.893149e+03     3.181151e+00
 * time: 2.4414801597595215
    64     1.892978e+03     4.518295e+00
 * time: 2.4515881538391113
    65     1.892743e+03     4.891794e+00
 * time: 2.4791030883789062
    66     1.892560e+03     2.854770e+00
 * time: 2.489549160003662
    67     1.892492e+03     5.319881e-01
 * time: 2.5000250339508057
    68     1.892485e+03     3.170874e-01
 * time: 2.5102131366729736
    69     1.892485e+03     3.053453e-01
 * time: 2.5199050903320312
    70     1.892485e+03     2.989026e-01
 * time: 2.5293312072753906
    71     1.892485e+03     2.979731e-01
 * time: 2.5383780002593994
    72     1.892485e+03     2.963326e-01
 * time: 2.5476531982421875
    73     1.892484e+03     3.456800e-01
 * time: 2.5568630695343018
    74     1.892483e+03     6.755932e-01
 * time: 2.566376209259033
    75     1.892479e+03     1.153383e+00
 * time: 2.593351125717163
    76     1.892471e+03     1.831175e+00
 * time: 2.6035749912261963
    77     1.892456e+03     2.458147e+00
 * time: 2.613793134689331
    78     1.892438e+03     2.315779e+00
 * time: 2.6241180896759033
    79     1.892419e+03     1.130926e+00
 * time: 2.6341440677642822
    80     1.892413e+03     1.293234e-01
 * time: 2.6444411277770996
    81     1.892412e+03     1.470565e-01
 * time: 2.654313087463379
    82     1.892412e+03     1.267892e-01
 * time: 2.663856029510498
    83     1.892412e+03     1.278265e-01
 * time: 2.673272132873535
    84     1.892412e+03     1.268709e-01
 * time: 2.682443141937256
    85     1.892412e+03     1.261773e-01
 * time: 2.708437204360962
    86     1.892412e+03     1.244539e-01
 * time: 2.7179811000823975
    87     1.892412e+03     1.207614e-01
 * time: 2.7275149822235107
    88     1.892412e+03     1.243552e-01
 * time: 2.737248182296753
    89     1.892412e+03     2.023712e-01
 * time: 2.747756004333496
    90     1.892411e+03     2.967994e-01
 * time: 2.758338212966919
    91     1.892410e+03     3.389817e-01
 * time: 2.7685840129852295
    92     1.892409e+03     2.414201e-01
 * time: 2.7786710262298584
    93     1.892409e+03     7.754118e-02
 * time: 2.788763999938965
    94     1.892409e+03     5.508363e-03
 * time: 2.79887318611145
    95     1.892409e+03     2.812552e-03
 * time: 2.826968193054199
    96     1.892409e+03     6.800357e-04
 * time: 2.8367230892181396

FittedPumasModel

Dynamical system type:                 Closed form

Number of subjects:                             30

Observation records:         Active        Missing
    DV:                         540              0
    Total:                      540              0

Number of parameters:      Constant      Optimized
                                  0             10

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                    -1892.409

--------------------
           Estimate
--------------------
tvvc       3.9757
tvq        0.042176
tvvp       3.6431
tvcl_hep   0.058575
tvcl_ren   0.13369
Ω₁,₁       0.18298
Ω₂,₂       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 ~ @. CombinedNormal(cp, σ_add, σ_prop)
    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, Ω = Diagonal([0.01, 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: 1.7881393432617188e-5
     1     3.641387e+03     1.432080e+03
 * time: 1.4757299423217773
     2     3.290450e+03     5.274782e+02
 * time: 1.4923250675201416
     3     3.185512e+03     2.173676e+02
 * time: 1.506932020187378
     4     3.143009e+03     1.479653e+02
 * time: 1.5205919742584229
     5     3.128511e+03     8.980031e+01
 * time: 1.5340509414672852
     6     3.123188e+03     5.033125e+01
 * time: 1.7481350898742676
     7     3.120794e+03     4.279722e+01
 * time: 1.7605350017547607
     8     3.118627e+03     3.971051e+01
 * time: 1.773149013519287
     9     3.115300e+03     8.456587e+01
 * time: 1.7857918739318848
    10     3.109353e+03     1.350354e+02
 * time: 1.798780918121338
    11     3.095894e+03     1.998258e+02
 * time: 1.8123879432678223
    12     2.567465e+03     7.201745e+02
 * time: 1.9108178615570068
    13     2.529114e+03     7.046190e+02
 * time: 1.9851958751678467
    14     2.302037e+03     4.880204e+02
 * time: 2.0587759017944336
    15     2.236337e+03     2.222133e+02
 * time: 2.0751118659973145
    16     2.181176e+03     2.344816e+02
 * time: 2.0992319583892822
    17     2.118163e+03     6.376161e+02
 * time: 2.114422082901001
    18     2.003998e+03     9.598838e+01
 * time: 2.1295390129089355
    19     1.986094e+03     8.304443e+01
 * time: 2.143065929412842
    20     1.979945e+03     7.240571e+01
 * time: 2.1564159393310547
    21     1.976200e+03     6.750758e+01
 * time: 2.169174909591675
    22     1.945029e+03     7.750595e+01
 * time: 2.1828019618988037
    23     1.920404e+03     7.521072e+01
 * time: 2.1959969997406006
    24     1.903076e+03     4.974067e+01
 * time: 2.208997964859009
    25     1.898051e+03     2.838327e+01
 * time: 2.221829891204834
    26     1.897204e+03     2.016238e+01
 * time: 2.25811505317688
    27     1.896962e+03     1.015686e+01
 * time: 2.2713730335235596
    28     1.896789e+03     7.310338e+00
 * time: 2.2846059799194336
    29     1.896674e+03     7.545736e+00
 * time: 2.298079013824463
    30     1.896557e+03     7.875214e+00
 * time: 2.3113460540771484
    31     1.896411e+03     8.138772e+00
 * time: 2.3256449699401855
    32     1.896174e+03     1.087082e+01
 * time: 2.3386409282684326
    33     1.895763e+03     1.563926e+01
 * time: 2.351284980773926
    34     1.895106e+03     1.762106e+01
 * time: 2.364288091659546
    35     1.894377e+03     1.439104e+01
 * time: 2.377457857131958
    36     1.893916e+03     8.187571e+00
 * time: 2.3930740356445312
    37     1.893796e+03     5.171054e+00
 * time: 2.426542043685913
    38     1.893768e+03     4.974975e+00
 * time: 2.439379930496216
    39     1.893727e+03     7.233885e+00
 * time: 2.4523448944091797
    40     1.893648e+03     9.421195e+00
 * time: 2.465501070022583
    41     1.893504e+03     1.004714e+01
 * time: 2.4787259101867676
    42     1.893324e+03     6.975364e+00
 * time: 2.4918930530548096
    43     1.893189e+03     7.168677e+00
 * time: 2.504833936691284
    44     1.893132e+03     5.410856e+00
 * time: 2.5174078941345215
    45     1.893108e+03     4.228610e+00
 * time: 2.5301759243011475
    46     1.893086e+03     3.955321e+00
 * time: 2.5639760494232178
    47     1.893064e+03     3.508037e+00
 * time: 2.5775539875030518
    48     1.893036e+03     3.521010e+00
 * time: 2.5910089015960693
    49     1.892987e+03     3.440055e+00
 * time: 2.6043078899383545
    50     1.892888e+03     4.921734e+00
 * time: 2.6175289154052734
    51     1.892722e+03     7.806123e+00
 * time: 2.6308679580688477
    52     1.892549e+03     7.367283e+00
 * time: 2.644603967666626
    53     1.892467e+03     3.687521e+00
 * time: 2.6579740047454834
    54     1.892446e+03     3.383982e+00
 * time: 2.690937042236328
    55     1.892439e+03     3.384463e+00
 * time: 2.7027199268341064
    56     1.892419e+03     3.366585e+00
 * time: 2.7142210006713867
    57     1.892372e+03     3.712901e+00
 * time: 2.7263660430908203
    58     1.892258e+03     7.013905e+00
 * time: 2.7385849952697754
    59     1.892026e+03     1.056922e+01
 * time: 2.7508020401000977
    60     1.891729e+03     1.126894e+01
 * time: 2.763421058654785
    61     1.891558e+03     6.626390e+00
 * time: 2.7764689922332764
    62     1.891475e+03     3.785422e+00
 * time: 2.8090288639068604
    63     1.891461e+03     3.662750e+00
 * time: 2.822956085205078
    64     1.891446e+03     3.596873e+00
 * time: 2.8355650901794434
    65     1.891392e+03     3.413681e+00
 * time: 2.8483500480651855
    66     1.891271e+03     3.322480e+00
 * time: 2.8611950874328613
    67     1.890959e+03     3.991911e+00
 * time: 2.87431001663208
    68     1.890360e+03     6.504714e+00
 * time: 2.8875038623809814
    69     1.889705e+03     6.158568e+00
 * time: 2.9006710052490234
    70     1.889326e+03     2.571572e+00
 * time: 2.9135239124298096
    71     1.889264e+03     8.652629e-01
 * time: 2.946518898010254
    72     1.889245e+03     7.083531e-01
 * time: 2.9594008922576904
    73     1.889242e+03     7.399495e-01
 * time: 2.9721970558166504
    74     1.889242e+03     7.391629e-01
 * time: 2.984135866165161
    75     1.889241e+03     7.293367e-01
 * time: 2.9962899684906006
    76     1.889239e+03     7.130926e-01
 * time: 3.0085859298706055
    77     1.889233e+03     6.763961e-01
 * time: 3.0208969116210938
    78     1.889219e+03     6.853163e-01
 * time: 3.033747911453247
    79     1.889185e+03     9.625001e-01
 * time: 3.0667190551757812
    80     1.889118e+03     1.659399e+00
 * time: 3.0803709030151367
    81     1.889019e+03     2.321116e+00
 * time: 3.0940489768981934
    82     1.888850e+03     1.927098e+00
 * time: 3.1078438758850098
    83     1.888700e+03     4.913841e-01
 * time: 3.1213550567626953
    84     1.888636e+03     1.431532e+00
 * time: 3.1346499919891357
    85     1.888568e+03     7.084613e-01
 * time: 3.1481339931488037
    86     1.888513e+03     9.767219e-01
 * time: 3.160569906234741
    87     1.888496e+03     6.084234e-01
 * time: 3.207813024520874
    88     1.888487e+03     4.675043e-01
 * time: 3.2205419540405273
    89     1.888485e+03     4.145920e-01
 * time: 3.2326560020446777
    90     1.888484e+03     4.635758e-01
 * time: 3.2446069717407227
    91     1.888484e+03     4.872007e-01
 * time: 3.256016969680786
    92     1.888484e+03     5.034784e-01
 * time: 3.267333984375
    93     1.888484e+03     5.319932e-01
 * time: 3.2786879539489746
    94     1.888483e+03     5.729332e-01
 * time: 3.290195941925049
    95     1.888480e+03     6.334098e-01
 * time: 3.3012900352478027
    96     1.888474e+03     7.067096e-01
 * time: 3.330306053161621
    97     1.888460e+03     7.633173e-01
 * time: 3.341531991958618
    98     1.888433e+03     7.288444e-01
 * time: 3.3529820442199707
    99     1.888398e+03     6.721772e-01
 * time: 3.3645200729370117
   100     1.888385e+03     2.363196e-01
 * time: 3.376631021499634
   101     1.888382e+03     3.523493e-01
 * time: 3.3950090408325195
   102     1.888381e+03     2.371289e-01
 * time: 3.406406879425049
   103     1.888379e+03     9.721378e-02
 * time: 3.4173898696899414
   104     1.888377e+03     2.046632e-01
 * time: 3.4476449489593506
   105     1.888377e+03     1.619384e-01
 * time: 3.4582178592681885
   106     1.888376e+03     6.972147e-02
 * time: 3.468751907348633
   107     1.888376e+03     6.877586e-02
 * time: 3.478864908218384
   108     1.888376e+03     6.867761e-02
 * time: 3.488636016845703
   109     1.888376e+03     6.847873e-02
 * time: 3.4988489151000977
   110     1.888376e+03     6.807833e-02
 * time: 3.5092949867248535
   111     1.888376e+03     1.160184e-01
 * time: 3.519834041595459
   112     1.888376e+03     1.941013e-01
 * time: 3.548646926879883
   113     1.888376e+03     3.104284e-01
 * time: 3.5599100589752197
   114     1.888375e+03     4.477116e-01
 * time: 3.571368932723999
   115     1.888374e+03     5.279254e-01
 * time: 3.5828700065612793
   116     1.888373e+03     4.174944e-01
 * time: 3.5939290523529053
   117     1.888372e+03     1.669162e-01
 * time: 3.604851007461548
   118     1.888372e+03     2.661730e-02
 * time: 3.6155478954315186
   119     1.888372e+03     1.954690e-02
 * time: 3.6258840560913086
   120     1.888372e+03     1.969381e-02
 * time: 3.6358299255371094
   121     1.888372e+03     1.972538e-02
 * time: 3.6638309955596924
   122     1.888372e+03     1.975794e-02
 * time: 3.674304962158203
   123     1.888372e+03     1.979852e-02
 * time: 3.684870958328247
   124     1.888372e+03     1.983616e-02
 * time: 3.6953940391540527
   125     1.888372e+03     1.989178e-02
 * time: 3.705798864364624
   126     1.888372e+03     1.992609e-02
 * time: 3.715914011001587
   127     1.888372e+03     2.374528e-02
 * time: 3.726170063018799
   128     1.888372e+03     4.233811e-02
 * time: 3.736661911010742
   129     1.888371e+03     6.803769e-02
 * time: 3.7473058700561523
   130     1.888370e+03     9.993454e-02
 * time: 3.7759110927581787
   131     1.888369e+03     1.445654e-01
 * time: 3.7872278690338135
   132     1.888367e+03     1.764869e-01
 * time: 3.7980809211730957
   133     1.888366e+03     1.250137e-01
 * time: 3.8089818954467773
   134     1.888365e+03     8.191620e-02
 * time: 3.8196940422058105
   135     1.888364e+03     2.379575e-02
 * time: 3.8317480087280273
   136     1.888364e+03     5.500305e-02
 * time: 3.8426709175109863
   137     1.888364e+03     2.796876e-02
 * time: 3.853381872177124
   138     1.888364e+03     1.621753e-02
 * time: 3.883059024810791
   139     1.888364e+03     1.360564e-02
 * time: 3.894031047821045
   140     1.888364e+03     4.224404e-03
 * time: 3.9050889015197754
   141     1.888364e+03     1.416970e-03
 * time: 3.916106939315796
   142     1.888364e+03     1.720473e-03
 * time: 3.926668882369995
   143     1.888364e+03     7.975571e-04
 * time: 3.936660051345825

FittedPumasModel

Dynamical system type:                 Closed form

Number of subjects:                             30

Observation records:         Active        Missing
    DV:                         540              0
    Total:                      540              0

Number of parameters:      Constant      Optimized
                                  0             13

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -1888.3638

-----------------------
             Estimate
-----------------------
tvvc         3.976
tvq          0.04239
tvvp         3.725
tvcl_hep_M   4.5558e-7
tvcl_hep_F   0.13349
tvcl_ren_M   0.19378
tvcl_ren_F   0.042196
Ω₁,₁         0.14046
Ω₂,₂         0.081349
σ_add        0.032171
σ_prop       0.061007
dCRCL_M      0.94823
dCRCL_F      1.9407
-----------------------

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)
16×2 DataFrame
Row Metric Value
String Any
1 Estimation Time 9.612
2 Subjects 30
3 Fixed Parameters 0
4 Optimized Parameters 8
5 DV Active Observations 540
6 DV Missing Observations 0
7 Total Active Observations 540
8 Total Missing Observations 0
9 Likelihood Approximation FOCE
10 LogLikelihood (LL) -1901.82
11 -2LL 3803.63
12 AIC 3819.63
13 BIC 3853.96
14 (η-shrinkage) η₁ -0.015
15 (η-shrinkage) η₂ -0.013
16 (ϵ-shrinkage) DV 0.056 
metrics_table(fit_covariate_model_wt)
16×2 DataFrame
Row Metric Value
String Any
1 Estimation Time 1.999
2 Subjects 30
3 Fixed Parameters 0
4 Optimized Parameters 8
5 DV Active Observations 540
6 DV Missing Observations 0
7 Total Active Observations 540
8 Total Missing Observations 0
9 Likelihood Approximation FOCE
10 LogLikelihood (LL) -1896.46
11 -2LL 3792.93
12 AIC 3808.93
13 BIC 3843.26
14 (η-shrinkage) η₁ -0.014
15 (η-shrinkage) η₂ -0.012
16 (ϵ-shrinkage) DV 0.056 
metrics_table(fit_covariate_model_wt_crcl)
16×2 DataFrame
Row Metric Value
String Any
1 Estimation Time 2.837
2 Subjects 30
3 Fixed Parameters 0
4 Optimized Parameters 10
5 DV Active Observations 540
6 DV Missing Observations 0
7 Total Active Observations 540
8 Total Missing Observations 0
9 Likelihood Approximation FOCE
10 LogLikelihood (LL) -1892.41
11 -2LL 3784.82
12 AIC 3804.82
13 BIC 3847.73
14 (η-shrinkage) η₁ -0.014
15 (η-shrinkage) η₂ -0.012
16 (ϵ-shrinkage) DV 0.056 
metrics_table(fit_covariate_model_wt_crcl_sex)
16×2 DataFrame
Row Metric Value
String Any
1 Estimation Time 3.937
2 Subjects 30
3 Fixed Parameters 0
4 Optimized Parameters 13
5 DV Active Observations 540
6 DV Missing Observations 0
7 Total Active Observations 540
8 Total Missing Observations 0
9 Likelihood Approximation FOCE
10 LogLikelihood (LL) -1888.36
11 -2LL 3776.73
12 AIC 3802.73
13 BIC 3858.52
14 (η-shrinkage) η₁ -0.013
15 (η-shrinkage) η₂ -0.012
16 (ϵ-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,
)
16×5 DataFrame
Row Metric Base_Model Covariate_Model_WT Covariate_Model_WT_CRCL Covariate_Model_WT_CRCL_SEX
String Any Any Any Any
1 Estimation Time 9.612 1.999 2.837 3.937
2 Subjects 30 30 30 30
3 Fixed Parameters 0 0 0 0
4 Optimized Parameters 8 8 10 13
5 DV Active Observations 540 540 540 540
6 DV Missing Observations 0 0 0 0
7 Total Active Observations 540 540 540 540
8 Total Missing Observations 0 0 0 0
9 Likelihood Approximation FOCE FOCE FOCE FOCE
10 LogLikelihood (LL) -1901.82 -1896.46 -1892.41 -1888.36
11 -2LL 3803.63 3792.93 3784.82 3776.73
12 AIC 3819.63 3808.93 3804.82 3802.73
13 BIC 3853.96 3843.26 3847.73 3858.52
14 (η-shrinkage) η₁ -0.015 -0.014 -0.014 -0.013
15 (η-shrinkage) η₂ -0.013 -0.012 -0.012 -0.012
16 (ϵ-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.629985029065
aic(fit_covariate_model_wt)
3808.9264612693205
aic(fit_covariate_model_wt_crcl)
3804.817947427335
aic(fit_covariate_model_wt_crcl_sex)
3802.7275279302753

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 dose control parameters.
[ Info: Evaluating individual parameters.
[ Info: Done.

FittedPumasModelInspection

Likelihood 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: 5.602836608886719e-5
     1     2.994747e+03     1.083462e+03
 * time: 2.3060760498046875
     2     2.406265e+03     1.884408e+02
 * time: 2.311228036880493
     3     2.344175e+03     7.741610e+01
 * time: 2.3161098957061768
     4     2.323153e+03     2.907642e+01
 * time: 2.320791006088257
     5     2.318222e+03     2.273295e+01
 * time: 2.3261239528656006
     6     2.316833e+03     1.390527e+01
 * time: 2.3310060501098633
     7     2.316425e+03     4.490883e+00
 * time: 2.3357698917388916
     8     2.316362e+03     9.374519e-01
 * time: 2.3403818607330322
     9     2.316356e+03     1.928785e-01
 * time: 2.345396041870117
    10     2.316355e+03     3.119615e-02
 * time: 2.966567039489746
    11     2.316355e+03     6.215513e-03
 * time: 2.970592975616455
    12     2.316355e+03     8.313206e-04
 * time: 2.974536895751953

FittedPumasModel

Dynamical system type:          No dynamical model

Number of subjects:                            160

Observation records:         Active        Missing
    painord:                   1920              0
    Total:                     1920              0

Number of parameters:      Constant      Optimized
                                  0              4

Likelihood approximation:              NaivePooled
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -2316.3554

-----------------
        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

Pumas will impute missing covariate values by constant interpolation of the available covariate 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.

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

Hence, for NOCB, you can use the following:

pop = read_pumas(pkdata; covariates_direction = :right)

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 LOCF or NOCB 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.