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, Ω = [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.04379701614379883
     1     3.110315e+03     9.706222e+02
 * time: 2.837930917739868
     2     2.831659e+03     7.817006e+02
 * time: 2.886122941970825
     3     2.405281e+03     2.923696e+02
 * time: 2.9242238998413086
     4     2.370406e+03     3.032286e+02
 * time: 2.9528019428253174
     5     2.313631e+03     3.126188e+02
 * time: 4.58235502243042
     6     2.263986e+03     2.946697e+02
 * time: 4.612850904464722
     7     2.160182e+03     1.917599e+02
 * time: 4.651752948760986
     8     2.112467e+03     1.588288e+02
 * time: 4.68664288520813
     9     2.090339e+03     1.108334e+02
 * time: 4.712116003036499
    10     2.078171e+03     8.108027e+01
 * time: 4.737342834472656
    11     2.074517e+03     7.813928e+01
 * time: 4.758900880813599
    12     2.066270e+03     7.044632e+01
 * time: 4.7811338901519775
    13     2.049660e+03     1.062584e+02
 * time: 4.806419849395752
    14     2.021965e+03     1.130570e+02
 * time: 4.8346848487854
    15     1.994936e+03     7.825801e+01
 * time: 4.862445831298828
    16     1.979337e+03     5.318263e+01
 * time: 4.88920783996582
    17     1.972141e+03     6.807046e+01
 * time: 4.915241003036499
    18     1.967973e+03     7.896361e+01
 * time: 4.940997838973999
    19     1.962237e+03     8.343757e+01
 * time: 4.966983795166016
    20     1.952791e+03     5.565304e+01
 * time: 4.9939048290252686
    21     1.935857e+03     3.923284e+01
 * time: 5.084208965301514
    22     1.926254e+03     5.749643e+01
 * time: 5.109196901321411
    23     1.922144e+03     4.306225e+01
 * time: 5.132999897003174
    24     1.911566e+03     4.810496e+01
 * time: 5.157912969589233
    25     1.906464e+03     4.324267e+01
 * time: 5.183523893356323
    26     1.905339e+03     1.207954e+01
 * time: 5.2057578563690186
    27     1.905092e+03     1.093479e+01
 * time: 5.227697849273682
    28     1.904957e+03     1.057034e+01
 * time: 5.249720811843872
    29     1.904875e+03     1.060882e+01
 * time: 5.270972013473511
    30     1.904459e+03     1.031525e+01
 * time: 5.292738914489746
    31     1.903886e+03     1.041179e+01
 * time: 5.3155457973480225
    32     1.903313e+03     1.135672e+01
 * time: 5.3388848304748535
    33     1.903057e+03     1.075683e+01
 * time: 5.382948875427246
    34     1.902950e+03     1.091295e+01
 * time: 5.4069390296936035
    35     1.902887e+03     1.042409e+01
 * time: 5.430078029632568
    36     1.902640e+03     9.213043e+00
 * time: 5.453382968902588
    37     1.902364e+03     9.519321e+00
 * time: 5.4765849113464355
    38     1.902156e+03     5.590984e+00
 * time: 5.499449014663696
    39     1.902094e+03     1.811898e+00
 * time: 5.521169900894165
    40     1.902086e+03     1.644722e+00
 * time: 5.542716979980469
    41     1.902084e+03     1.653520e+00
 * time: 5.5638158321380615
    42     1.902074e+03     1.720184e+00
 * time: 5.585951805114746
    43     1.902055e+03     2.619061e+00
 * time: 5.6267218589782715
    44     1.902015e+03     3.519729e+00
 * time: 5.649840831756592
    45     1.901962e+03     3.403372e+00
 * time: 5.672708988189697
    46     1.901924e+03     1.945644e+00
 * time: 5.696163892745972
    47     1.901914e+03     6.273342e-01
 * time: 5.717460870742798
    48     1.901913e+03     5.374557e-01
 * time: 5.738102912902832
    49     1.901913e+03     5.710007e-01
 * time: 5.758168935775757
    50     1.901913e+03     6.091390e-01
 * time: 5.794642925262451
    51     1.901912e+03     6.692417e-01
 * time: 5.816102027893066
    52     1.901909e+03     7.579620e-01
 * time: 5.8375959396362305
    53     1.901903e+03     8.798211e-01
 * time: 5.858814001083374
    54     1.901889e+03     1.002981e+00
 * time: 5.879491806030273
    55     1.901864e+03     1.495512e+00
 * time: 5.898345947265625
    56     1.901837e+03     1.664621e+00
 * time: 5.917364835739136
    57     1.901819e+03     8.601119e-01
 * time: 5.951814889907837
    58     1.901815e+03     4.525470e-02
 * time: 5.973102807998657
    59     1.901815e+03     1.294280e-02
 * time: 5.994143962860107
    60     1.901815e+03     2.876567e-03
 * time: 6.013682842254639
    61     1.901815e+03     2.876567e-03
 * time: 6.067133903503418
    62     1.901815e+03     2.876513e-03
 * time: 6.099836826324463
    63     1.901815e+03     2.876513e-03
 * time: 6.142007827758789
    64     1.901815e+03     2.876513e-03
 * time: 6.191589832305908

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:                      NoXChange
Log-likelihood value:                    -1901.815

------------------
         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 ~ @. 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: 4.291534423828125e-5
     1     2.846050e+03     1.128657e+03
 * time: 1.226660966873169
     2     2.472982e+03     7.008264e+02
 * time: 1.2602200508117676
     3     2.180690e+03     2.312704e+02
 * time: 1.2930009365081787
     4     2.125801e+03     1.862929e+02
 * time: 1.5326738357543945
     5     2.103173e+03     1.320946e+02
 * time: 1.5624480247497559
     6     2.091136e+03     1.103035e+02
 * time: 1.5905818939208984
     7     2.081443e+03     1.091133e+02
 * time: 1.6175260543823242
     8     2.071793e+03     7.197675e+01
 * time: 1.644704818725586
     9     2.062706e+03     7.623310e+01
 * time: 1.671888828277588
    10     2.057515e+03     6.885476e+01
 * time: 1.6989648342132568
    11     2.051133e+03     6.368504e+01
 * time: 1.7252819538116455
    12     2.038626e+03     7.730243e+01
 * time: 1.7524080276489258
    13     2.019352e+03     1.136864e+02
 * time: 1.8281619548797607
    14     1.997136e+03     1.005748e+02
 * time: 1.8575918674468994
    15     1.983023e+03     6.831478e+01
 * time: 1.8864669799804688
    16     1.977700e+03     5.649783e+01
 * time: 1.914154052734375
    17     1.974583e+03     6.357642e+01
 * time: 1.9412689208984375
    18     1.967292e+03     7.658974e+01
 * time: 1.9698028564453125
    19     1.951045e+03     6.130573e+01
 * time: 2.0018138885498047
    20     1.935868e+03     4.845839e+01
 * time: 2.048818826675415
    21     1.929356e+03     6.325782e+01
 * time: 2.0791728496551514
    22     1.925187e+03     3.142245e+01
 * time: 2.1074140071868896
    23     1.923733e+03     4.623400e+01
 * time: 2.134768009185791
    24     1.918498e+03     5.347738e+01
 * time: 2.162644863128662
    25     1.912383e+03     5.849125e+01
 * time: 2.2121808528900146
    26     1.905510e+03     3.254038e+01
 * time: 2.2437548637390137
    27     1.903629e+03     2.905618e+01
 * time: 2.271888017654419
    28     1.902833e+03     2.907696e+01
 * time: 2.298275947570801
    29     1.902447e+03     2.746037e+01
 * time: 2.3415050506591797
    30     1.899399e+03     1.930949e+01
 * time: 2.370429039001465
    31     1.898705e+03     1.186361e+01
 * time: 2.3977789878845215
    32     1.898505e+03     1.050402e+01
 * time: 2.4245200157165527
    33     1.898474e+03     1.042186e+01
 * time: 2.465268850326538
    34     1.897992e+03     1.238729e+01
 * time: 2.491096019744873
    35     1.897498e+03     1.729368e+01
 * time: 2.5165200233459473
    36     1.896954e+03     1.472554e+01
 * time: 2.54145884513855
    37     1.896744e+03     5.852825e+00
 * time: 2.5800299644470215
    38     1.896705e+03     1.171353e+00
 * time: 2.603317975997925
    39     1.896704e+03     1.216117e+00
 * time: 2.6269028186798096
    40     1.896703e+03     1.230336e+00
 * time: 2.6651790142059326
    41     1.896698e+03     1.250715e+00
 * time: 2.691074848175049
    42     1.896688e+03     1.449552e+00
 * time: 2.717634916305542
    43     1.896666e+03     2.533300e+00
 * time: 2.7437899112701416
    44     1.896631e+03     3.075537e+00
 * time: 2.7846100330352783
    45     1.896599e+03     2.139797e+00
 * time: 2.811224937438965
    46     1.896587e+03     6.636030e-01
 * time: 2.837601900100708
    47     1.896585e+03     6.303517e-01
 * time: 2.8645098209381104
    48     1.896585e+03     5.995265e-01
 * time: 2.9051530361175537
    49     1.896584e+03     5.844159e-01
 * time: 2.9314498901367188
    50     1.896583e+03     6.083858e-01
 * time: 2.9575438499450684
    51     1.896579e+03     8.145327e-01
 * time: 2.9828479290008545
    52     1.896570e+03     1.293490e+00
 * time: 3.0232059955596924
    53     1.896549e+03     1.877705e+00
 * time: 3.0490479469299316
    54     1.896513e+03     2.217392e+00
 * time: 3.0738158226013184
    55     1.896482e+03     1.658148e+00
 * time: 3.0988218784332275
    56     1.896466e+03     5.207218e-01
 * time: 3.1394619941711426
    57     1.896463e+03     1.177468e-01
 * time: 3.1660079956054688
    58     1.896463e+03     1.678937e-02
 * time: 3.1908328533172607
    59     1.896463e+03     2.666636e-03
 * time: 3.2139148712158203
    60     1.896463e+03     2.666636e-03
 * time: 3.280349016189575
    61     1.896463e+03     2.666636e-03
 * time: 3.3151779174804688

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:                      NoXChange
Log-likelihood value:                   -1896.4632

------------------
         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 ~ @. 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, Ω = [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: 2.6941299438476562e-5
     1     3.572050e+03     1.302046e+03
 * time: 1.2338230609893799
     2     3.266947e+03     5.384244e+02
 * time: 1.5520520210266113
     3     3.150635e+03     1.918079e+02
 * time: 1.585827112197876
     4     3.108122e+03     1.277799e+02
 * time: 1.618001937866211
     5     3.091358e+03     8.838080e+01
 * time: 1.649873971939087
     6     3.082997e+03     8.321689e+01
 * time: 1.68184494972229
     7     3.076379e+03     8.167702e+01
 * time: 1.715338945388794
     8     3.067422e+03     1.177822e+02
 * time: 1.7501099109649658
     9     3.048580e+03     2.526969e+02
 * time: 1.7883410453796387
    10     2.993096e+03     6.325396e+02
 * time: 1.9044280052185059
    11     2.965744e+03     7.634718e+02
 * time: 1.9996769428253174
    12     2.921212e+03     9.704020e+02
 * time: 2.0596699714660645
    13     2.553649e+03     6.642510e+02
 * time: 2.10904598236084
    14     2.319495e+03     3.204711e+02
 * time: 2.2031800746917725
    15     2.183040e+03     2.174905e+02
 * time: 2.2564148902893066
    16     2.157621e+03     3.150983e+02
 * time: 2.288645029067993
    17     2.132395e+03     2.847614e+02
 * time: 2.342194080352783
    18     2.084799e+03     1.563370e+02
 * time: 2.37426495552063
    19     2.071497e+03     1.006429e+02
 * time: 2.4055819511413574
    20     2.064983e+03     9.753313e+01
 * time: 2.4369800090789795
    21     2.048289e+03     9.230405e+01
 * time: 2.4680569171905518
    22     2.020938e+03     1.292359e+02
 * time: 2.5160109996795654
    23     1.983888e+03     2.284990e+02
 * time: 2.5489389896392822
    24     1.962132e+03     1.220188e+02
 * time: 2.5794100761413574
    25     1.945947e+03     1.035894e+02
 * time: 2.6231400966644287
    26     1.917782e+03     8.246698e+01
 * time: 2.654670000076294
    27     1.905967e+03     5.408054e+01
 * time: 2.685408115386963
    28     1.898569e+03     2.172222e+01
 * time: 2.7293810844421387
    29     1.897473e+03     1.689350e+01
 * time: 2.7600479125976562
    30     1.897019e+03     1.666689e+01
 * time: 2.7889161109924316
    31     1.896796e+03     1.699751e+01
 * time: 2.8313279151916504
    32     1.896559e+03     1.645900e+01
 * time: 2.860243082046509
    33     1.896223e+03     1.415504e+01
 * time: 2.888284921646118
    34     1.895773e+03     1.630081e+01
 * time: 2.929616928100586
    35     1.895309e+03     1.723930e+01
 * time: 2.9589450359344482
    36     1.895004e+03     1.229983e+01
 * time: 2.987394094467163
    37     1.894871e+03     5.385102e+00
 * time: 3.0283870697021484
    38     1.894827e+03     3.465463e+00
 * time: 3.0567281246185303
    39     1.894816e+03     3.387474e+00
 * time: 3.0835471153259277
    40     1.894807e+03     3.295388e+00
 * time: 3.123141050338745
    41     1.894786e+03     3.089194e+00
 * time: 3.1508140563964844
    42     1.894737e+03     2.928080e+00
 * time: 3.1782119274139404
    43     1.894624e+03     3.088723e+00
 * time: 3.2188289165496826
    44     1.894413e+03     3.493791e+00
 * time: 3.247464895248413
    45     1.894127e+03     3.142865e+00
 * time: 3.275094985961914
    46     1.893933e+03     2.145253e+00
 * time: 3.3172190189361572
    47     1.893888e+03     2.172800e+00
 * time: 3.3454079627990723
    48     1.893880e+03     2.180509e+00
 * time: 3.372529983520508
    49     1.893876e+03     2.134257e+00
 * time: 3.41145396232605
    50     1.893868e+03     2.032354e+00
 * time: 3.439038038253784
    51     1.893846e+03     1.760874e+00
 * time: 3.4662420749664307
    52     1.893796e+03     1.779016e+00
 * time: 3.5056540966033936
    53     1.893694e+03     2.018956e+00
 * time: 3.5337259769439697
    54     1.893559e+03     2.366854e+00
 * time: 3.5613460540771484
    55     1.893474e+03     3.690142e+00
 * time: 3.588714122772217
    56     1.893446e+03     3.675109e+00
 * time: 3.6294400691986084
    57     1.893439e+03     3.426419e+00
 * time: 3.6565299034118652
    58     1.893429e+03     3.183164e+00
 * time: 3.68306303024292
    59     1.893398e+03     2.695171e+00
 * time: 3.7232260704040527
    60     1.893328e+03     2.753548e+00
 * time: 3.750899076461792
    61     1.893169e+03     3.589748e+00
 * time: 3.778485059738159
    62     1.892920e+03     3.680718e+00
 * time: 3.8202719688415527
    63     1.892667e+03     2.568107e+00
 * time: 3.8484549522399902
    64     1.892514e+03     1.087910e+00
 * time: 3.87595796585083
    65     1.892493e+03     3.287296e-01
 * time: 3.916156053543091
    66     1.892492e+03     2.967465e-01
 * time: 3.943542003631592
    67     1.892492e+03     3.020682e-01
 * time: 3.9697799682617188
    68     1.892491e+03     3.034704e-01
 * time: 4.007508039474487
    69     1.892491e+03     3.091846e-01
 * time: 4.033973932266235
    70     1.892491e+03     3.224170e-01
 * time: 4.059741973876953
    71     1.892490e+03     6.494197e-01
 * time: 4.098314046859741
    72     1.892488e+03     1.115188e+00
 * time: 4.125349044799805
    73     1.892483e+03     1.838833e+00
 * time: 4.151883125305176
    74     1.892472e+03     2.765371e+00
 * time: 4.191416025161743
    75     1.892452e+03     3.463807e+00
 * time: 4.219326972961426
    76     1.892431e+03     2.805270e+00
 * time: 4.246336936950684
    77     1.892411e+03     5.758916e-01
 * time: 4.286294937133789
    78     1.892410e+03     1.434041e-01
 * time: 4.315022945404053
    79     1.892409e+03     1.639246e-01
 * time: 4.341756105422974
    80     1.892409e+03     1.145856e-01
 * time: 4.380342960357666
    81     1.892409e+03     3.966861e-02
 * time: 4.407128095626831
    82     1.892409e+03     3.550808e-02
 * time: 4.432642936706543
    83     1.892409e+03     3.456241e-02
 * time: 4.457125902175903
    84     1.892409e+03     3.114018e-02
 * time: 4.495525121688843
    85     1.892409e+03     4.080806e-02
 * time: 4.5211029052734375
    86     1.892409e+03     6.722726e-02
 * time: 4.546617031097412
    87     1.892409e+03     1.006791e-01
 * time: 4.587049961090088
    88     1.892409e+03     1.303988e-01
 * time: 4.613785028457642
    89     1.892409e+03     1.228919e-01
 * time: 4.639107942581177
    90     1.892409e+03     6.433813e-02
 * time: 4.678416013717651
    91     1.892409e+03     1.314164e-02
 * time: 4.704900026321411
    92     1.892409e+03     4.929931e-04
 * time: 4.730042934417725

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.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 ~ @. 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, Ω = [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: 3.123283386230469e-5
     1     3.641387e+03     1.432080e+03
 * time: 1.0070462226867676
     2     3.290450e+03     5.274782e+02
 * time: 1.0445802211761475
     3     3.185512e+03     2.173676e+02
 * time: 1.2887182235717773
     4     3.143009e+03     1.479653e+02
 * time: 1.3216822147369385
     5     3.128511e+03     8.980031e+01
 * time: 1.3522751331329346
     6     3.123188e+03     5.033125e+01
 * time: 1.3820011615753174
     7     3.120794e+03     4.279722e+01
 * time: 1.4113051891326904
     8     3.118627e+03     3.971051e+01
 * time: 1.4413340091705322
     9     3.115300e+03     8.456587e+01
 * time: 1.473046064376831
    10     3.109353e+03     1.350354e+02
 * time: 1.57078218460083
    11     3.095894e+03     1.998258e+02
 * time: 1.6036090850830078
    12     2.988214e+03     4.366433e+02
 * time: 1.6469831466674805
    13     2.896081e+03     5.505943e+02
 * time: 1.7400932312011719
    14     2.652467e+03     7.300323e+02
 * time: 2.3210952281951904
    15     2.560937e+03     6.973661e+02
 * time: 2.507571220397949
    16     2.254941e+03     2.740033e+02
 * time: 2.542879104614258
    17     2.222509e+03     2.034303e+02
 * time: 2.586560010910034
    18     2.171255e+03     2.449580e+02
 * time: 2.618098020553589
    19     2.024532e+03     1.121511e+02
 * time: 2.64872407913208
    20     1.993723e+03     1.042814e+02
 * time: 2.690258026123047
    21     1.985113e+03     8.079014e+01
 * time: 2.719144105911255
    22     1.976757e+03     7.054196e+01
 * time: 2.759411096572876
    23     1.969970e+03     6.070322e+01
 * time: 2.7877161502838135
    24     1.961095e+03     6.810782e+01
 * time: 2.8160111904144287
    25     1.947983e+03     8.116920e+01
 * time: 2.8578200340270996
    26     1.930371e+03     8.530051e+01
 * time: 2.886719226837158
    27     1.910209e+03     6.993170e+01
 * time: 2.9276671409606934
    28     1.899107e+03     3.362640e+01
 * time: 2.9571280479431152
    29     1.898022e+03     2.642220e+01
 * time: 2.9850242137908936
    30     1.897055e+03     1.213144e+01
 * time: 3.025090217590332
    31     1.896596e+03     7.773239e+00
 * time: 3.053441047668457
    32     1.896538e+03     7.997039e+00
 * time: 3.080660104751587
    33     1.896451e+03     8.160909e+00
 * time: 3.1197760105133057
    34     1.896283e+03     8.237721e+00
 * time: 3.1473770141601562
    35     1.895903e+03     1.520219e+01
 * time: 3.1866660118103027
    36     1.895272e+03     2.358916e+01
 * time: 3.218550205230713
    37     1.894536e+03     2.461296e+01
 * time: 3.2634291648864746
    38     1.893995e+03     1.546128e+01
 * time: 3.322864055633545
    39     1.893858e+03     6.976137e+00
 * time: 3.365976095199585
    40     1.893833e+03     6.019466e+00
 * time: 3.4074831008911133
    41     1.893786e+03     3.827201e+00
 * time: 3.464370012283325
    42     1.893714e+03     3.323412e+00
 * time: 3.5027430057525635
    43     1.893592e+03     3.215150e+00
 * time: 3.542478084564209
    44     1.893435e+03     6.534965e+00
 * time: 3.570855140686035
    45     1.893286e+03     7.424154e+00
 * time: 3.599133014678955
    46     1.893190e+03     5.552627e+00
 * time: 3.6405441761016846
    47     1.893139e+03     3.222316e+00
 * time: 3.6689670085906982
    48     1.893120e+03     3.015339e+00
 * time: 3.707782030105591
    49     1.893107e+03     3.244809e+00
 * time: 3.7349741458892822
    50     1.893080e+03     6.163100e+00
 * time: 3.7619781494140625
    51     1.893027e+03     9.824713e+00
 * time: 3.801140069961548
    52     1.892912e+03     1.390100e+01
 * time: 3.8293511867523193
    53     1.892734e+03     1.510937e+01
 * time: 3.8574681282043457
    54     1.892561e+03     1.008563e+01
 * time: 3.8966221809387207
    55     1.892485e+03     3.730668e+00
 * time: 3.9246761798858643
    56     1.892471e+03     3.380261e+00
 * time: 3.9639980792999268
    57     1.892463e+03     3.167904e+00
 * time: 3.9913392066955566
    58     1.892441e+03     4.152065e+00
 * time: 4.018607139587402
    59     1.892391e+03     7.355996e+00
 * time: 4.069102048873901
    60     1.892268e+03     1.195397e+01
 * time: 4.097748041152954
    61     1.892026e+03     1.640783e+01
 * time: 4.126658201217651
    62     1.891735e+03     1.593576e+01
 * time: 4.169431209564209
    63     1.891569e+03     8.316423e+00
 * time: 4.198472023010254
    64     1.891494e+03     3.948212e+00
 * time: 4.2385852336883545
    65     1.891481e+03     3.911593e+00
 * time: 4.267411231994629
    66     1.891457e+03     3.875559e+00
 * time: 4.296016216278076
    67     1.891405e+03     3.811247e+00
 * time: 4.338113069534302
    68     1.891262e+03     3.657045e+00
 * time: 4.366690158843994
    69     1.890930e+03     4.957405e+00
 * time: 4.395815134048462
    70     1.890317e+03     6.657726e+00
 * time: 4.438716173171997
    71     1.889660e+03     6.086302e+00
 * time: 4.469796180725098
    72     1.889303e+03     2.270929e+00
 * time: 4.512453079223633
    73     1.889253e+03     7.695301e-01
 * time: 4.541966199874878
    74     1.889252e+03     7.382144e-01
 * time: 4.571489095687866
    75     1.889251e+03     7.187898e-01
 * time: 4.613083124160767
    76     1.889251e+03     7.215047e-01
 * time: 4.641733169555664
    77     1.889250e+03     7.235155e-01
 * time: 4.670827150344849
    78     1.889249e+03     7.246818e-01
 * time: 4.712966203689575
    79     1.889244e+03     7.257796e-01
 * time: 4.742214202880859
    80     1.889233e+03     7.198189e-01
 * time: 4.771003007888794
    81     1.889204e+03     1.089029e+00
 * time: 4.8128721714019775
    82     1.889142e+03     1.801601e+00
 * time: 4.844545125961304
    83     1.889043e+03     2.967917e+00
 * time: 4.889618158340454
    84     1.888889e+03     2.965856e+00
 * time: 4.920137166976929
    85     1.888705e+03     5.933554e-01
 * time: 4.95054817199707
    86     1.888655e+03     9.577700e-01
 * time: 4.993028163909912
    87     1.888582e+03     1.498495e+00
 * time: 5.023543119430542
    88     1.888533e+03     1.502750e+00
 * time: 5.053840160369873
    89     1.888490e+03     1.184664e+00
 * time: 5.098099231719971
    90     1.888480e+03     6.684512e-01
 * time: 5.129159212112427
    91     1.888476e+03     3.680030e-01
 * time: 5.172979116439819
    92     1.888476e+03     4.720039e-01
 * time: 5.203711032867432
    93     1.888476e+03     4.768646e-01
 * time: 5.23405909538269
    94     1.888475e+03     4.736674e-01
 * time: 5.277422189712524
    95     1.888475e+03     4.552766e-01
 * time: 5.307950019836426
    96     1.888474e+03     5.193717e-01
 * time: 5.339208126068115
    97     1.888473e+03     8.850084e-01
 * time: 5.384765148162842
    98     1.888468e+03     1.461596e+00
 * time: 5.416812181472778
    99     1.888458e+03     2.209122e+00
 * time: 5.458933115005493
   100     1.888437e+03     2.961233e+00
 * time: 5.488408088684082
   101     1.888407e+03     2.978461e+00
 * time: 5.517944097518921
   102     1.888384e+03     1.707197e+00
 * time: 5.563612222671509
   103     1.888381e+03     6.198617e-01
 * time: 5.596243143081665
   104     1.888380e+03     5.171253e-01
 * time: 5.627691030502319
   105     1.888378e+03     1.037235e-01
 * time: 5.67261004447937
   106     1.888378e+03     8.473258e-02
 * time: 5.702542066574097
   107     1.888378e+03     8.364952e-02
 * time: 5.745103120803833
   108     1.888378e+03     8.080437e-02
 * time: 5.775638103485107
   109     1.888378e+03     7.873895e-02
 * time: 5.805823087692261
   110     1.888378e+03     7.798398e-02
 * time: 5.849322080612183
   111     1.888378e+03     7.788171e-02
 * time: 5.87748122215271
   112     1.888378e+03     7.776461e-02
 * time: 5.905781030654907
   113     1.888378e+03     9.023510e-02
 * time: 5.947090148925781
   114     1.888378e+03     1.631350e-01
 * time: 5.976325035095215
   115     1.888378e+03     2.768651e-01
 * time: 6.005245208740234
   116     1.888377e+03     4.462240e-01
 * time: 6.047398090362549
   117     1.888377e+03     6.643039e-01
 * time: 6.076444149017334
   118     1.888375e+03     8.432956e-01
 * time: 6.117724180221558
   119     1.888374e+03     7.596136e-01
 * time: 6.147057056427002
   120     1.888373e+03     3.637586e-01
 * time: 6.176208019256592
   121     1.888372e+03     8.304420e-02
 * time: 6.217543125152588
   122     1.888372e+03     2.084519e-02
 * time: 6.246016025543213
   123     1.888372e+03     2.056413e-02
 * time: 6.2736241817474365
   124     1.888372e+03     2.044077e-02
 * time: 6.313430070877075
   125     1.888372e+03     2.035198e-02
 * time: 6.341545104980469
   126     1.888372e+03     2.021268e-02
 * time: 6.381649017333984
   127     1.888372e+03     1.998173e-02
 * time: 6.408853054046631
   128     1.888372e+03     3.162545e-02
 * time: 6.436462163925171
   129     1.888372e+03     5.510789e-02
 * time: 6.477373123168945
   130     1.888372e+03     9.278492e-02
 * time: 6.505300998687744
   131     1.888372e+03     1.529182e-01
 * time: 6.5339391231536865
   132     1.888372e+03     2.462454e-01
 * time: 6.57564902305603
   133     1.888372e+03     3.800394e-01
 * time: 6.60486102104187
   134     1.888371e+03     5.313041e-01
 * time: 6.633714199066162
   135     1.888369e+03     6.020486e-01
 * time: 6.675993204116821
   136     1.888367e+03     4.665794e-01
 * time: 6.705169200897217
   137     1.888366e+03     1.404900e-01
 * time: 6.746726036071777
   138     1.888365e+03     8.513229e-02
 * time: 6.775470018386841
   139     1.888364e+03     1.244490e-01
 * time: 6.80377721786499
   140     1.888364e+03     1.028375e-01
 * time: 6.845769166946411
   141     1.888364e+03     5.164409e-02
 * time: 6.874897003173828
   142     1.888364e+03     5.147608e-02
 * time: 6.903462171554565
   143     1.888364e+03     3.147211e-02
 * time: 6.964107036590576
   144     1.888364e+03     2.104429e-02
 * time: 6.992862224578857
   145     1.888364e+03     6.544008e-03
 * time: 7.020952224731445
   146     1.888364e+03     2.536819e-03
 * time: 7.062792062759399
   147     1.888364e+03     4.361067e-03
 * time: 7.090376138687134
   148     1.888364e+03     3.035268e-03
 * time: 7.130414009094238
   149     1.888364e+03     5.968930e-04
 * time: 7.1570351123809814

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.7249
tvcl_hep_M   1.7174e-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 6.192
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 3.315
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 4.73
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 7.157
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 6.192 3.315 4.73 7.157
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.6299849303396
aic(fit_covariate_model_wt)
3808.9264607805962
aic(fit_covariate_model_wt_crcl)
3804.817947371706
aic(fit_covariate_model_wt_crcl_sex)
3802.7275243739386

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: 2.5987625122070312e-5
     1     2.994747e+03     1.083462e+03
 * time: 2.5820400714874268
     2     2.406265e+03     1.884408e+02
 * time: 2.5878090858459473
     3     2.344175e+03     7.741610e+01
 * time: 2.5927979946136475
     4     2.323153e+03     2.907642e+01
 * time: 2.5978829860687256
     5     2.318222e+03     2.273295e+01
 * time: 2.603188991546631
     6     2.316833e+03     1.390527e+01
 * time: 2.608703136444092
     7     2.316425e+03     4.490883e+00
 * time: 2.614586114883423
     8     2.316362e+03     9.374519e-01
 * time: 2.6201770305633545
     9     2.316356e+03     1.928785e-01
 * time: 2.6257810592651367
    10     2.316355e+03     3.119615e-02
 * time: 2.6314029693603516
    11     2.316355e+03     6.215513e-03
 * time: 2.6370720863342285
    12     2.316355e+03     8.313206e-04
 * time: 2.643019914627075

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.