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.

1 Covariate Model Building

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

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

1.3 Step 2 - Base Model

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

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

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

    @random begin
        η ~ MvNormal(Ω)
    end

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

    @dynamics Central1Periph1

    @derived begin
        cp := @. Central / Vc
        DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
    end
end
PumasModel
  Parameters: tvcl, tvvc, tvq, tvvp, Ω, σ_add, σ_prop
  Random effects: η
  Covariates: 
  Dynamical variables: Central, Peripheral
  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.020631790161132812
     1     3.110315e+03     9.706222e+02
 * time: 0.4215998649597168
     2     2.831659e+03     7.817006e+02
 * time: 0.43834400177001953
     3     2.405281e+03     2.923696e+02
 * time: 0.47801780700683594
     4     2.370406e+03     3.032286e+02
 * time: 0.4879429340362549
     5     2.313631e+03     3.126188e+02
 * time: 0.49889183044433594
     6     2.263986e+03     2.946697e+02
 * time: 0.5264909267425537
     7     2.160182e+03     1.917599e+02
 * time: 0.5395100116729736
     8     2.112467e+03     1.588288e+02
 * time: 0.5506808757781982
     9     2.090339e+03     1.108334e+02
 * time: 0.5589849948883057
    10     2.078171e+03     8.108027e+01
 * time: 0.5751988887786865
    11     2.074517e+03     7.813928e+01
 * time: 0.5833818912506104
    12     2.066270e+03     7.044632e+01
 * time: 0.591113805770874
    13     2.049660e+03     1.062584e+02
 * time: 0.5992169380187988
    14     2.021965e+03     1.130570e+02
 * time: 0.6150429248809814
    15     1.994936e+03     7.825801e+01
 * time: 0.6235218048095703
    16     1.979337e+03     5.318263e+01
 * time: 0.631633996963501
    17     1.972141e+03     6.807046e+01
 * time: 0.6396708488464355
    18     1.967973e+03     7.896361e+01
 * time: 0.6551549434661865
    19     1.962237e+03     8.343757e+01
 * time: 0.6635808944702148
    20     1.952791e+03     5.565304e+01
 * time: 0.6719658374786377
    21     1.935857e+03     3.923284e+01
 * time: 0.6805219650268555
    22     1.926254e+03     5.749643e+01
 * time: 0.6964828968048096
    23     1.922144e+03     4.306225e+01
 * time: 0.7047228813171387
    24     1.911566e+03     4.810496e+01
 * time: 0.712993860244751
    25     1.906464e+03     4.324267e+01
 * time: 0.7210898399353027
    26     1.905339e+03     1.207954e+01
 * time: 0.7364490032196045
    27     1.905092e+03     1.093479e+01
 * time: 0.7442679405212402
    28     1.904957e+03     1.057034e+01
 * time: 0.7518928050994873
    29     1.904875e+03     1.060882e+01
 * time: 0.7591288089752197
    30     1.904459e+03     1.031525e+01
 * time: 0.7741978168487549
    31     1.903886e+03     1.041179e+01
 * time: 0.7825889587402344
    32     1.903313e+03     1.135672e+01
 * time: 0.7905318737030029
    33     1.903057e+03     1.075683e+01
 * time: 0.7981619834899902
    34     1.902950e+03     1.091295e+01
 * time: 0.8131418228149414
    35     1.902887e+03     1.042409e+01
 * time: 0.8214499950408936
    36     1.902640e+03     9.213043e+00
 * time: 0.8291938304901123
    37     1.902364e+03     9.519321e+00
 * time: 0.8367469310760498
    38     1.902156e+03     5.590984e+00
 * time: 0.8513679504394531
    39     1.902094e+03     1.811898e+00
 * time: 0.8595700263977051
    40     1.902086e+03     1.644722e+00
 * time: 0.8672277927398682
    41     1.902084e+03     1.653520e+00
 * time: 0.8764729499816895
    42     1.902074e+03     1.720184e+00
 * time: 0.8838658332824707
    43     1.902055e+03     2.619061e+00
 * time: 0.9036688804626465
    44     1.902015e+03     3.519729e+00
 * time: 0.9114398956298828
    45     1.901962e+03     3.403372e+00
 * time: 0.9188809394836426
    46     1.901924e+03     1.945644e+00
 * time: 0.926196813583374
    47     1.901914e+03     6.273342e-01
 * time: 0.9405670166015625
    48     1.901913e+03     5.374557e-01
 * time: 0.9485499858856201
    49     1.901913e+03     5.710007e-01
 * time: 0.9557828903198242
    50     1.901913e+03     6.091390e-01
 * time: 0.9628889560699463
    51     1.901912e+03     6.692417e-01
 * time: 0.9774580001831055
    52     1.901909e+03     7.579620e-01
 * time: 0.9854249954223633
    53     1.901903e+03     8.798211e-01
 * time: 0.9930980205535889
    54     1.901889e+03     1.002981e+00
 * time: 1.000459909439087
    55     1.901864e+03     1.495512e+00
 * time: 1.0077078342437744
    56     1.901837e+03     1.664621e+00
 * time: 1.0220239162445068
    57     1.901819e+03     8.601119e-01
 * time: 1.0298428535461426
    58     1.901815e+03     4.525470e-02
 * time: 1.037226915359497
    59     1.901815e+03     1.294280e-02
 * time: 1.0441668033599854
    60     1.901815e+03     2.876567e-03
 * time: 1.057811975479126
    61     1.901815e+03     2.876567e-03
 * time: 1.0733819007873535
    62     1.901815e+03     2.876567e-03
 * time: 1.0900278091430664
FittedPumasModel

Successful minimization:                      true

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

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

1.4 Step 3 - Covariate Model

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

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

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

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

Tip

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

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

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

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

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates begin
        WT
    end

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

    @dynamics Central1Periph1

    @derived begin
        cp := @. Central / Vc
        DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
    end
end
PumasModel
  Parameters: tvcl, tvvc, tvq, tvvp, Ω, σ_add, σ_prop
  Random effects: η
  Covariates: WT
  Dynamical variables: Central, Peripheral
  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: 9.107589721679688e-5
     1     2.846050e+03     1.128657e+03
 * time: 0.040435075759887695
     2     2.472982e+03     7.008264e+02
 * time: 0.050436973571777344
     3     2.180690e+03     2.312704e+02
 * time: 0.05978798866271973
     4     2.125801e+03     1.862929e+02
 * time: 0.0868680477142334
     5     2.103173e+03     1.320946e+02
 * time: 0.09629297256469727
     6     2.091136e+03     1.103035e+02
 * time: 0.10457992553710938
     7     2.081443e+03     1.091133e+02
 * time: 0.11224794387817383
     8     2.071793e+03     7.197675e+01
 * time: 0.13861513137817383
     9     2.062706e+03     7.623310e+01
 * time: 0.14760994911193848
    10     2.057515e+03     6.885476e+01
 * time: 0.15587711334228516
    11     2.051133e+03     6.368504e+01
 * time: 0.16358304023742676
    12     2.038626e+03     7.730243e+01
 * time: 0.1787581443786621
    13     2.019352e+03     1.136864e+02
 * time: 0.18777799606323242
    14     1.997136e+03     1.005748e+02
 * time: 0.195997953414917
    15     1.983023e+03     6.831478e+01
 * time: 0.20414400100708008
    16     1.977700e+03     5.649783e+01
 * time: 0.2192249298095703
    17     1.974583e+03     6.357642e+01
 * time: 0.22799301147460938
    18     1.967292e+03     7.658974e+01
 * time: 0.23650097846984863
    19     1.951045e+03     6.130573e+01
 * time: 0.24562907218933105
    20     1.935868e+03     4.845839e+01
 * time: 0.2613980770111084
    21     1.929356e+03     6.325782e+01
 * time: 0.27101802825927734
    22     1.925187e+03     3.142245e+01
 * time: 0.2792649269104004
    23     1.923733e+03     4.623400e+01
 * time: 0.2869880199432373
    24     1.918498e+03     5.347738e+01
 * time: 0.3021869659423828
    25     1.912383e+03     5.849125e+01
 * time: 0.3120839595794678
    26     1.905510e+03     3.254038e+01
 * time: 0.3208010196685791
    27     1.903629e+03     2.905618e+01
 * time: 0.3285551071166992
    28     1.902833e+03     2.907696e+01
 * time: 0.3439500331878662
    29     1.902447e+03     2.746037e+01
 * time: 0.35215306282043457
    30     1.899399e+03     1.930949e+01
 * time: 0.3603339195251465
    31     1.898705e+03     1.186361e+01
 * time: 0.3680839538574219
    32     1.898505e+03     1.050402e+01
 * time: 0.3830990791320801
    33     1.898474e+03     1.042186e+01
 * time: 0.391085147857666
    34     1.897992e+03     1.238729e+01
 * time: 0.398684024810791
    35     1.897498e+03     1.729368e+01
 * time: 0.4061450958251953
    36     1.896954e+03     1.472554e+01
 * time: 0.4139590263366699
    37     1.896744e+03     5.852825e+00
 * time: 0.4296259880065918
    38     1.896705e+03     1.171353e+00
 * time: 0.4370110034942627
    39     1.896704e+03     1.216117e+00
 * time: 0.44471192359924316
    40     1.896703e+03     1.230336e+00
 * time: 0.4517550468444824
    41     1.896698e+03     1.250715e+00
 * time: 0.4663350582122803
    42     1.896688e+03     1.449552e+00
 * time: 0.47430896759033203
    43     1.896666e+03     2.533300e+00
 * time: 0.48177313804626465
    44     1.896631e+03     3.075537e+00
 * time: 0.48914504051208496
    45     1.896599e+03     2.139797e+00
 * time: 0.5038719177246094
    46     1.896587e+03     6.636030e-01
 * time: 0.5123310089111328
    47     1.896585e+03     6.303517e-01
 * time: 0.5198140144348145
    48     1.896585e+03     5.995265e-01
 * time: 0.526871919631958
    49     1.896584e+03     5.844159e-01
 * time: 0.5338249206542969
    50     1.896583e+03     6.083858e-01
 * time: 0.5487279891967773
    51     1.896579e+03     8.145327e-01
 * time: 0.5564429759979248
    52     1.896570e+03     1.293490e+00
 * time: 0.5639610290527344
    53     1.896549e+03     1.877705e+00
 * time: 0.5711479187011719
    54     1.896513e+03     2.217392e+00
 * time: 0.5856101512908936
    55     1.896482e+03     1.658148e+00
 * time: 0.5938820838928223
    56     1.896466e+03     5.207218e-01
 * time: 0.6014001369476318
    57     1.896463e+03     1.177468e-01
 * time: 0.6086559295654297
    58     1.896463e+03     1.678937e-02
 * time: 0.6154670715332031
    59     1.896463e+03     2.666636e-03
 * time: 0.6302881240844727
    60     1.896463e+03     2.666636e-03
 * time: 0.6438541412353516
    61     1.896463e+03     2.666636e-03
 * time: 0.6648819446563721
FittedPumasModel

Successful minimization:                      true

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

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

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

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

Tip

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

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

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

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

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates begin
        WT
        CRCL
    end

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

    @dynamics Central1Periph1

    @derived begin
        cp := @. Central / Vc
        DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
    end
end
PumasModel
  Parameters: tvvc, tvq, tvvp, tvcl_hep, tvcl_ren, Ω, σ_add, σ_prop, dCRCL
  Random effects: η
  Covariates: WT, CRCL
  Dynamical variables: Central, Peripheral
  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: 7.700920104980469e-5
     1     3.572050e+03     1.302046e+03
 * time: 0.04581093788146973
     2     3.266947e+03     5.384244e+02
 * time: 0.06013202667236328
     3     3.150635e+03     1.918079e+02
 * time: 0.09422183036804199
     4     3.108122e+03     1.277799e+02
 * time: 0.10559296607971191
     5     3.091358e+03     8.838080e+01
 * time: 0.11717796325683594
     6     3.082997e+03     8.321689e+01
 * time: 0.14975285530090332
     7     3.076379e+03     8.167702e+01
 * time: 0.16101789474487305
     8     3.067422e+03     1.177822e+02
 * time: 0.1725449562072754
     9     3.048580e+03     2.526969e+02
 * time: 0.195969820022583
    10     2.993096e+03     6.325396e+02
 * time: 0.21298003196716309
    11     2.965744e+03     7.634718e+02
 * time: 0.26103782653808594
    12     2.921212e+03     9.704020e+02
 * time: 0.2821059226989746
    13     2.553649e+03     6.642510e+02
 * time: 0.309643030166626
    14     2.319495e+03     3.204711e+02
 * time: 0.3347358703613281
    15     2.183040e+03     2.174905e+02
 * time: 0.3637678623199463
    16     2.157621e+03     3.150983e+02
 * time: 0.375777006149292
    17     2.132395e+03     2.847614e+02
 * time: 0.3979058265686035
    18     2.084799e+03     1.563370e+02
 * time: 0.4089479446411133
    19     2.071497e+03     1.006429e+02
 * time: 0.41950201988220215
    20     2.064983e+03     9.753313e+01
 * time: 0.43219900131225586
    21     2.048289e+03     9.230405e+01
 * time: 0.4532029628753662
    22     2.020938e+03     1.292359e+02
 * time: 0.46386194229125977
    23     1.983888e+03     2.284990e+02
 * time: 0.4767940044403076
    24     1.962132e+03     1.220188e+02
 * time: 0.49767589569091797
    25     1.945947e+03     1.035894e+02
 * time: 0.5079150199890137
    26     1.917782e+03     8.246698e+01
 * time: 0.5200209617614746
    27     1.905967e+03     5.408054e+01
 * time: 0.5409879684448242
    28     1.898569e+03     2.172222e+01
 * time: 0.5519859790802002
    29     1.897473e+03     1.689350e+01
 * time: 0.5630238056182861
    30     1.897019e+03     1.666689e+01
 * time: 0.5839059352874756
    31     1.896796e+03     1.699751e+01
 * time: 0.5944638252258301
    32     1.896559e+03     1.645900e+01
 * time: 0.6045289039611816
    33     1.896223e+03     1.415504e+01
 * time: 0.6166138648986816
    34     1.895773e+03     1.630081e+01
 * time: 0.6369359493255615
    35     1.895309e+03     1.723930e+01
 * time: 0.6473908424377441
    36     1.895004e+03     1.229983e+01
 * time: 0.6593899726867676
    37     1.894871e+03     5.385102e+00
 * time: 0.6798548698425293
    38     1.894827e+03     3.465463e+00
 * time: 0.6902070045471191
    39     1.894816e+03     3.387474e+00
 * time: 0.7005548477172852
    40     1.894807e+03     3.295388e+00
 * time: 0.7211658954620361
    41     1.894786e+03     3.089194e+00
 * time: 0.7316398620605469
    42     1.894737e+03     2.928080e+00
 * time: 0.7415218353271484
    43     1.894624e+03     3.088723e+00
 * time: 0.7532608509063721
    44     1.894413e+03     3.493791e+00
 * time: 0.7736618518829346
    45     1.894127e+03     3.142865e+00
 * time: 0.7839508056640625
    46     1.893933e+03     2.145253e+00
 * time: 0.7951059341430664
    47     1.893888e+03     2.172800e+00
 * time: 0.8157448768615723
    48     1.893880e+03     2.180509e+00
 * time: 0.8262758255004883
    49     1.893876e+03     2.134257e+00
 * time: 0.8358478546142578
    50     1.893868e+03     2.032354e+00
 * time: 0.8474688529968262
    51     1.893846e+03     1.760874e+00
 * time: 0.8670518398284912
    52     1.893796e+03     1.779016e+00
 * time: 0.8770928382873535
    53     1.893694e+03     2.018956e+00
 * time: 0.8881638050079346
    54     1.893559e+03     2.366854e+00
 * time: 0.9086039066314697
    55     1.893474e+03     3.690142e+00
 * time: 0.9190428256988525
    56     1.893446e+03     3.675109e+00
 * time: 0.9290549755096436
    57     1.893439e+03     3.426419e+00
 * time: 0.9426689147949219
    58     1.893429e+03     3.183164e+00
 * time: 0.9800939559936523
    59     1.893398e+03     2.695171e+00
 * time: 0.9899828433990479
    60     1.893328e+03     2.753548e+00
 * time: 1.001209020614624
    61     1.893169e+03     3.589748e+00
 * time: 1.0216999053955078
    62     1.892920e+03     3.680718e+00
 * time: 1.0322628021240234
    63     1.892667e+03     2.568107e+00
 * time: 1.042766809463501
    64     1.892514e+03     1.087910e+00
 * time: 1.0549290180206299
    65     1.892493e+03     3.287296e-01
 * time: 1.0744478702545166
    66     1.892492e+03     2.967465e-01
 * time: 1.0841658115386963
    67     1.892492e+03     3.020682e-01
 * time: 1.0953810214996338
    68     1.892491e+03     3.034704e-01
 * time: 1.1148138046264648
    69     1.892491e+03     3.091846e-01
 * time: 1.1246488094329834
    70     1.892491e+03     3.224170e-01
 * time: 1.1346499919891357
    71     1.892490e+03     6.494197e-01
 * time: 1.1462018489837646
    72     1.892488e+03     1.115188e+00
 * time: 1.1654298305511475
    73     1.892483e+03     1.838833e+00
 * time: 1.175055980682373
    74     1.892472e+03     2.765371e+00
 * time: 1.186194896697998
    75     1.892452e+03     3.463807e+00
 * time: 1.2063708305358887
    76     1.892431e+03     2.805270e+00
 * time: 1.2165210247039795
    77     1.892411e+03     5.758916e-01
 * time: 1.226820945739746
    78     1.892410e+03     1.434041e-01
 * time: 1.2385399341583252
    79     1.892409e+03     1.639246e-01
 * time: 1.257807970046997
    80     1.892409e+03     1.145856e-01
 * time: 1.267244815826416
    81     1.892409e+03     3.966861e-02
 * time: 1.2780449390411377
    82     1.892409e+03     3.550808e-02
 * time: 1.2980358600616455
    83     1.892409e+03     3.456241e-02
 * time: 1.3074498176574707
    84     1.892409e+03     3.114018e-02
 * time: 1.3166329860687256
    85     1.892409e+03     4.080806e-02
 * time: 1.3277568817138672
    86     1.892409e+03     6.722726e-02
 * time: 1.346925973892212
    87     1.892409e+03     1.006791e-01
 * time: 1.3565759658813477
    88     1.892409e+03     1.303988e-01
 * time: 1.3668968677520752
    89     1.892409e+03     1.228919e-01
 * time: 1.3782248497009277
    90     1.892409e+03     6.433813e-02
 * time: 1.3972728252410889
    91     1.892409e+03     1.314164e-02
 * time: 1.4066219329833984
    92     1.892409e+03     4.929931e-04
 * time: 1.4171359539031982
FittedPumasModel

Successful minimization:                      true

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

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

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

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

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

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates begin
        WT
        CRCL
        SEX
    end

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

    @dynamics Central1Periph1

    @derived begin
        cp := @. Central / Vc
        DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
    end
end
PumasModel
  Parameters: tvvc, tvq, tvvp, tvcl_hep_M, tvcl_hep_F, tvcl_ren_M, tvcl_ren_F, Ω, σ_add, σ_prop, dCRCL_M, dCRCL_F
  Random effects: η
  Covariates: WT, CRCL, SEX
  Dynamical variables: Central, Peripheral
  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: 9.107589721679688e-5
     1     3.641387e+03     1.432080e+03
 * time: 0.021095991134643555
     2     3.290450e+03     5.274782e+02
 * time: 0.061689138412475586
     3     3.185512e+03     2.173676e+02
 * time: 0.07662701606750488
     4     3.143009e+03     1.479653e+02
 * time: 0.09034299850463867
     5     3.128511e+03     8.980031e+01
 * time: 0.12900304794311523
     6     3.123188e+03     5.033125e+01
 * time: 0.1422100067138672
     7     3.120794e+03     4.279722e+01
 * time: 0.16795015335083008
     8     3.118627e+03     3.971051e+01
 * time: 0.18190908432006836
     9     3.115300e+03     8.456587e+01
 * time: 0.19435405731201172
    10     3.109353e+03     1.350354e+02
 * time: 0.21968603134155273
    11     3.095894e+03     1.998258e+02
 * time: 0.23378801345825195
    12     2.988214e+03     4.366433e+02
 * time: 0.2646331787109375
    13     2.896081e+03     5.505943e+02
 * time: 0.3252739906311035
    14     2.652467e+03     7.300323e+02
 * time: 0.6229751110076904
    15     2.560937e+03     6.973661e+02
 * time: 0.6995081901550293
    16     2.254941e+03     2.740033e+02
 * time: 0.7147011756896973
    17     2.222509e+03     2.034303e+02
 * time: 0.7361791133880615
    18     2.171255e+03     2.449580e+02
 * time: 0.7499501705169678
    19     2.024532e+03     1.121511e+02
 * time: 0.7626512050628662
    20     1.993723e+03     1.042814e+02
 * time: 0.7837581634521484
    21     1.985113e+03     8.079014e+01
 * time: 0.796018123626709
    22     1.976757e+03     7.054196e+01
 * time: 0.8158860206604004
    23     1.969970e+03     6.070322e+01
 * time: 0.8304719924926758
    24     1.961095e+03     6.810782e+01
 * time: 0.8420100212097168
    25     1.947983e+03     8.116920e+01
 * time: 0.8622560501098633
    26     1.930371e+03     8.530051e+01
 * time: 0.8745801448822021
    27     1.910209e+03     6.993170e+01
 * time: 0.8862481117248535
    28     1.899107e+03     3.362640e+01
 * time: 0.9074411392211914
    29     1.898022e+03     2.642220e+01
 * time: 0.9187901020050049
    30     1.897055e+03     1.213144e+01
 * time: 0.938621997833252
    31     1.896596e+03     7.773239e+00
 * time: 0.9510290622711182
    32     1.896538e+03     7.997039e+00
 * time: 0.9626240730285645
    33     1.896451e+03     8.160909e+00
 * time: 0.9858531951904297
    34     1.896283e+03     8.237721e+00
 * time: 0.9979050159454346
    35     1.895903e+03     1.520219e+01
 * time: 1.009296178817749
    36     1.895272e+03     2.358916e+01
 * time: 1.032461166381836
    37     1.894536e+03     2.461296e+01
 * time: 1.044175148010254
    38     1.893995e+03     1.546128e+01
 * time: 1.0639710426330566
    39     1.893858e+03     6.976137e+00
 * time: 1.0760631561279297
    40     1.893833e+03     6.019466e+00
 * time: 1.0871641635894775
    41     1.893786e+03     3.827201e+00
 * time: 1.1068670749664307
    42     1.893714e+03     3.323412e+00
 * time: 1.1205620765686035
    43     1.893592e+03     3.215150e+00
 * time: 1.1323370933532715
    44     1.893435e+03     6.534965e+00
 * time: 1.154712200164795
    45     1.893286e+03     7.424154e+00
 * time: 1.166260004043579
    46     1.893190e+03     5.552627e+00
 * time: 1.1857709884643555
    47     1.893139e+03     3.222316e+00
 * time: 1.1978261470794678
    48     1.893120e+03     3.015339e+00
 * time: 1.2090981006622314
    49     1.893107e+03     3.244809e+00
 * time: 1.2285301685333252
    50     1.893080e+03     6.163100e+00
 * time: 1.2400519847869873
    51     1.893027e+03     9.824713e+00
 * time: 1.252007007598877
    52     1.892912e+03     1.390100e+01
 * time: 1.2746450901031494
    53     1.892734e+03     1.510937e+01
 * time: 1.286067008972168
    54     1.892561e+03     1.008563e+01
 * time: 1.305774211883545
    55     1.892485e+03     3.730668e+00
 * time: 1.3177900314331055
    56     1.892471e+03     3.380261e+00
 * time: 1.3287560939788818
    57     1.892463e+03     3.167904e+00
 * time: 1.3483970165252686
    58     1.892441e+03     4.152065e+00
 * time: 1.360140085220337
    59     1.892391e+03     7.355996e+00
 * time: 1.371016025543213
    60     1.892268e+03     1.195397e+01
 * time: 1.3910419940948486
    61     1.892026e+03     1.640783e+01
 * time: 1.4029350280761719
    62     1.891735e+03     1.593576e+01
 * time: 1.4161651134490967
    63     1.891569e+03     8.316423e+00
 * time: 1.4548251628875732
    64     1.891494e+03     3.948212e+00
 * time: 1.4659650325775146
    65     1.891481e+03     3.911593e+00
 * time: 1.48537015914917
    66     1.891457e+03     3.875559e+00
 * time: 1.497251033782959
    67     1.891405e+03     3.811247e+00
 * time: 1.508350133895874
    68     1.891262e+03     3.657045e+00
 * time: 1.5281062126159668
    69     1.890930e+03     4.957405e+00
 * time: 1.5400161743164062
    70     1.890317e+03     6.657726e+00
 * time: 1.5512981414794922
    71     1.889660e+03     6.086302e+00
 * time: 1.5723001956939697
    72     1.889303e+03     2.270929e+00
 * time: 1.583921194076538
    73     1.889253e+03     7.695301e-01
 * time: 1.6034882068634033
    74     1.889252e+03     7.382144e-01
 * time: 1.6248860359191895
    75     1.889251e+03     7.187898e-01
 * time: 1.6442670822143555
    76     1.889251e+03     7.215047e-01
 * time: 1.678934097290039
    77     1.889250e+03     7.235155e-01
 * time: 1.6906661987304688
    78     1.889249e+03     7.246818e-01
 * time: 1.701167106628418
    79     1.889244e+03     7.257796e-01
 * time: 1.7204670906066895
    80     1.889233e+03     7.198189e-01
 * time: 1.731783151626587
    81     1.889204e+03     1.089029e+00
 * time: 1.7426750659942627
    82     1.889142e+03     1.801602e+00
 * time: 1.7635581493377686
    83     1.889043e+03     2.967917e+00
 * time: 1.7747540473937988
    84     1.888889e+03     2.965856e+00
 * time: 1.7945561408996582
    85     1.888705e+03     5.933552e-01
 * time: 1.807574987411499
    86     1.888655e+03     9.577703e-01
 * time: 1.8187310695648193
    87     1.888582e+03     1.498495e+00
 * time: 1.8518280982971191
    88     1.888533e+03     1.502748e+00
 * time: 1.8635461330413818
    89     1.888490e+03     1.184664e+00
 * time: 1.8744430541992188
    90     1.888480e+03     6.684509e-01
 * time: 1.8945322036743164
    91     1.888476e+03     3.680027e-01
 * time: 1.9059171676635742
    92     1.888476e+03     4.720038e-01
 * time: 1.9164481163024902
    93     1.888476e+03     4.768645e-01
 * time: 1.9362990856170654
    94     1.888475e+03     4.736675e-01
 * time: 1.9469380378723145
    95     1.888475e+03     4.552767e-01
 * time: 1.9748170375823975
    96     1.888474e+03     5.193706e-01
 * time: 1.9863600730895996
    97     1.888473e+03     8.850065e-01
 * time: 1.996847152709961
    98     1.888468e+03     1.461593e+00
 * time: 2.023090124130249
    99     1.888458e+03     2.209118e+00
 * time: 2.0345849990844727
   100     1.888437e+03     2.961228e+00
 * time: 2.045243978500366
   101     1.888407e+03     2.978460e+00
 * time: 2.064877986907959
   102     1.888384e+03     1.707193e+00
 * time: 2.0762741565704346
   103     1.888381e+03     6.198122e-01
 * time: 2.0869832038879395
   104     1.888380e+03     5.171484e-01
 * time: 2.1065421104431152
   105     1.888378e+03     1.037119e-01
 * time: 2.1174800395965576
   106     1.888378e+03     8.473266e-02
 * time: 2.128180980682373
   107     1.888378e+03     8.364935e-02
 * time: 2.1476972103118896
   108     1.888378e+03     8.080431e-02
 * time: 2.160701036453247
   109     1.888378e+03     7.873888e-02
 * time: 2.179978132247925
   110     1.888378e+03     7.798398e-02
 * time: 2.1910030841827393
   111     1.888378e+03     7.788171e-02
 * time: 2.201259136199951
   112     1.888378e+03     7.776462e-02
 * time: 2.220189094543457
   113     1.888378e+03     9.023408e-02
 * time: 2.2313292026519775
   114     1.888378e+03     1.631324e-01
 * time: 2.2416141033172607
   115     1.888378e+03     2.768598e-01
 * time: 2.2603721618652344
   116     1.888377e+03     4.462142e-01
 * time: 2.2715260982513428
   117     1.888377e+03     6.642866e-01
 * time: 2.2820191383361816
   118     1.888375e+03     8.432659e-01
 * time: 2.301259994506836
   119     1.888374e+03     7.595679e-01
 * time: 2.3123509883880615
   120     1.888373e+03     3.637227e-01
 * time: 2.323167085647583
   121     1.888372e+03     8.303332e-02
 * time: 2.3426430225372314
   122     1.888372e+03     2.084522e-02
 * time: 2.353464126586914
   123     1.888372e+03     2.056408e-02
 * time: 2.3723390102386475
   124     1.888372e+03     2.044076e-02
 * time: 2.3855230808258057
   125     1.888372e+03     2.035198e-02
 * time: 2.3956310749053955
   126     1.888372e+03     2.021271e-02
 * time: 2.4180331230163574
   127     1.888372e+03     1.998180e-02
 * time: 2.4289710521698
   128     1.888372e+03     3.163155e-02
 * time: 2.439171075820923
   129     1.888372e+03     5.511846e-02
 * time: 2.4575741291046143
   130     1.888372e+03     9.280271e-02
 * time: 2.4688570499420166
   131     1.888372e+03     1.529474e-01
 * time: 2.479074001312256
   132     1.888372e+03     2.462917e-01
 * time: 2.4976611137390137
   133     1.888372e+03     3.801088e-01
 * time: 2.5088260173797607
   134     1.888371e+03     5.313967e-01
 * time: 2.5192811489105225
   135     1.888369e+03     6.021459e-01
 * time: 2.538605213165283
   136     1.888367e+03     4.666397e-01
 * time: 2.5515010356903076
   137     1.888366e+03     1.404876e-01
 * time: 2.5621700286865234
   138     1.888365e+03     8.513159e-02
 * time: 2.5817651748657227
   139     1.888364e+03     1.244768e-01
 * time: 2.592377185821533
   140     1.888364e+03     1.028569e-01
 * time: 2.603041172027588
   141     1.888364e+03     5.165873e-02
 * time: 2.6254870891571045
   142     1.888364e+03     5.146239e-02
 * time: 2.63581919670105
   143     1.888364e+03     3.147161e-02
 * time: 2.654792070388794
   144     1.888364e+03     2.104201e-02
 * time: 2.6658191680908203
   145     1.888364e+03     6.547270e-03
 * time: 2.6758601665496826
   146     1.888364e+03     2.534557e-03
 * time: 2.69449520111084
   147     1.888364e+03     4.359993e-03
 * time: 2.7052910327911377
   148     1.888364e+03     3.035840e-03
 * time: 2.715224027633667
   149     1.888364e+03     5.979031e-04
 * time: 2.7333672046661377
FittedPumasModel

Successful minimization:                      true

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

--------------------------
               Estimate
--------------------------
tvvc            3.976
tvq             0.04239
tvvp            3.7249
tvcl_hep_M      1.7173e-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.

1.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)
17×2 DataFrame
Row Metric Value
String Any
1 Successful true
2 Estimation Time 1.09
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
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 0.665
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
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 1.417
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
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 2.733
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
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 1.09 0.665 1.417 2.733
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 Pumas.FOCE Pumas.FOCE Pumas.FOCE
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.6299849528186
aic(fit_covariate_model_wt)
3808.926460780597
aic(fit_covariate_model_wt_crcl)
3804.817947371705
aic(fit_covariate_model_wt_crcl_sex)
3802.7275243737645

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.

1.5.1 Goodness of Fit Plots

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

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

Fitting was successful: true
Likehood approximation used for weighted residuals : Pumas.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.

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

1.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 variables: 
  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: 9.107589721679688e-5
     1     2.994747e+03     1.083462e+03
 * time: 0.005037069320678711
     2     2.406265e+03     1.884408e+02
 * time: 0.010031938552856445
     3     2.344175e+03     7.741610e+01
 * time: 0.06113100051879883
     4     2.323153e+03     2.907642e+01
 * time: 0.06500005722045898
     5     2.318222e+03     2.273295e+01
 * time: 0.0687859058380127
     6     2.316833e+03     1.390527e+01
 * time: 0.07235097885131836
     7     2.316425e+03     4.490883e+00
 * time: 0.07584404945373535
     8     2.316362e+03     9.374519e-01
 * time: 0.08009195327758789
     9     2.316356e+03     1.928785e-01
 * time: 0.08470296859741211
    10     2.316355e+03     3.119615e-02
 * time: 0.08965396881103516
    11     2.316355e+03     6.215513e-03
 * time: 0.09449005126953125
    12     2.316355e+03     8.313206e-04
 * time: 0.1421360969543457
FittedPumasModel

Successful minimization:                      true

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

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

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

1.7 Missing Data in Covariates

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

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

Hence, for LOCF, you can use the following:

pop = read_pumas(pkdata; covariates_direction = :left)

along with any other required keyword arguments for column mapping.

Note

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

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

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