Covariate Models

Authors

Jose Storopoli

Joel Owen 

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

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

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

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

1 nlme_sample Dataset

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

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

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

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

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

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

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

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

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

2 Step 1 - Parse Data into a Population

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

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

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

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

3 Step 2 - Base Model

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

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

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

    @random begin
        η ~ MvNormal(Ω)
    end

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

    @dynamics Central1Periph1

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

And fit it accordingly:

iparams_base_model = (;
    tvvc = 5,
    tvcl = 0.02,
    tvq = 0.01,
    tvvp = 10,
    Ω = Diagonal([0.01, 0.01]),
    σ_add = 0.1,
    σ_prop = 0.1,
)
(tvvc = 5,
 tvcl = 0.02,
 tvq = 0.01,
 tvvp = 10,
 Ω = [0.01 0.0; 0.0 0.01],
 σ_add = 0.1,
 σ_prop = 0.1,)
fit_base_model = fit(base_model, pop, iparams_base_model, FOCE())
[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter     Function value   Gradient norm 
     0     8.300164e+03     4.360770e+03
 * time: 0.04593396186828613
     1     3.110315e+03     9.706222e+02
 * time: 1.9061741828918457
     2     2.831659e+03     7.817006e+02
 * time: 1.9460821151733398
     3     2.405281e+03     2.923696e+02
 * time: 1.9768011569976807
     4     2.370406e+03     3.032286e+02
 * time: 1.9993419647216797
     5     2.313631e+03     3.126188e+02
 * time: 2.022578001022339
     6     2.263986e+03     2.946697e+02
 * time: 2.0437850952148438
     7     2.160182e+03     1.917599e+02
 * time: 2.0731520652770996
     8     2.112467e+03     1.588288e+02
 * time: 2.1554691791534424
     9     2.090339e+03     1.108334e+02
 * time: 2.17637300491333
    10     2.078171e+03     8.108027e+01
 * time: 2.198012113571167
    11     2.074517e+03     7.813928e+01
 * time: 2.2172281742095947
    12     2.066270e+03     7.044632e+01
 * time: 2.236515998840332
    13     2.049660e+03     1.062584e+02
 * time: 2.256103038787842
    14     2.021965e+03     1.130570e+02
 * time: 2.2757561206817627
    15     1.994936e+03     7.825801e+01
 * time: 2.2952780723571777
    16     1.979337e+03     5.318263e+01
 * time: 2.315155029296875
    17     1.972141e+03     6.807046e+01
 * time: 2.368110179901123
    18     1.967973e+03     7.896361e+01
 * time: 2.388524055480957
    19     1.962237e+03     8.343757e+01
 * time: 2.408830165863037
    20     1.952791e+03     5.565304e+01
 * time: 2.4297521114349365
    21     1.935857e+03     3.923284e+01
 * time: 2.450730085372925
    22     1.926254e+03     5.749643e+01
 * time: 2.4711790084838867
    23     1.922144e+03     4.306225e+01
 * time: 2.490638017654419
    24     1.911566e+03     4.810496e+01
 * time: 2.523993968963623
    25     1.906464e+03     4.324267e+01
 * time: 2.5450830459594727
    26     1.905339e+03     1.207954e+01
 * time: 2.5643789768218994
    27     1.905092e+03     1.093479e+01
 * time: 2.5825321674346924
    28     1.904957e+03     1.057034e+01
 * time: 2.601018190383911
    29     1.904875e+03     1.060882e+01
 * time: 2.6301000118255615
    30     1.904459e+03     1.031525e+01
 * time: 2.648987054824829
    31     1.903886e+03     1.041179e+01
 * time: 2.6681690216064453
    32     1.903313e+03     1.135672e+01
 * time: 2.6874771118164062
    33     1.903057e+03     1.075683e+01
 * time: 2.7059450149536133
    34     1.902950e+03     1.091295e+01
 * time: 2.7354249954223633
    35     1.902887e+03     1.042409e+01
 * time: 2.75412917137146
    36     1.902640e+03     9.213043e+00
 * time: 2.7729949951171875
    37     1.902364e+03     9.519321e+00
 * time: 2.79206919670105
    38     1.902156e+03     5.590984e+00
 * time: 2.820669174194336
    39     1.902094e+03     1.811898e+00
 * time: 2.8395111560821533
    40     1.902086e+03     1.644722e+00
 * time: 2.858236074447632
    41     1.902084e+03     1.653520e+00
 * time: 2.876270055770874
    42     1.902074e+03     1.720184e+00
 * time: 2.9045569896698
    43     1.902055e+03     2.619061e+00
 * time: 2.9231741428375244
    44     1.902015e+03     3.519729e+00
 * time: 2.941899061203003
    45     1.901962e+03     3.403372e+00
 * time: 2.9605581760406494
    46     1.901924e+03     1.945644e+00
 * time: 2.9891810417175293
    47     1.901914e+03     6.273342e-01
 * time: 3.007956027984619
    48     1.901913e+03     5.374557e-01
 * time: 3.026482105255127
    49     1.901913e+03     5.710007e-01
 * time: 3.0439460277557373
    50     1.901913e+03     6.091390e-01
 * time: 3.0710220336914062
    51     1.901912e+03     6.692417e-01
 * time: 3.0890121459960938
    52     1.901909e+03     7.579620e-01
 * time: 3.1068010330200195
    53     1.901903e+03     8.798211e-01
 * time: 3.124962091445923
    54     1.901889e+03     1.002981e+00
 * time: 3.153264045715332
    55     1.901864e+03     1.495512e+00
 * time: 3.1718859672546387
    56     1.901837e+03     1.664621e+00
 * time: 3.1899120807647705
    57     1.901819e+03     8.601119e-01
 * time: 3.207961082458496
    58     1.901815e+03     4.525470e-02
 * time: 3.235589027404785
    59     1.901815e+03     1.294280e-02
 * time: 3.2533321380615234
    60     1.901815e+03     2.876567e-03
 * time: 3.2700130939483643
    61     1.901815e+03     2.876567e-03
 * time: 3.3175711631774902
    62     1.901815e+03     2.876567e-03
 * time: 3.3527050018310547

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:              NoObjectiveChange
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 ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
    end
end
PumasModel
  Parameters: tvcl, tvvc, tvq, tvvp, Ω, σ_add, σ_prop
  Random effects: η
  Covariates: WT
  Dynamical system variables: Central, Peripheral
  Dynamical system type: Closed form
  Derived: DV
  Observed: DV

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

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

fit_covariate_model_wt = fit(covariate_model_wt, pop, iparams_base_model, FOCE())
[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter     Function value   Gradient norm 
     0     7.695401e+03     4.898919e+03
 * time: 2.384185791015625e-5
     1     2.846050e+03     1.128657e+03
 * time: 0.5385689735412598
     2     2.472982e+03     7.008264e+02
 * time: 0.6559848785400391
     3     2.180690e+03     2.312704e+02
 * time: 0.6808948516845703
     4     2.125801e+03     1.862929e+02
 * time: 0.7026798725128174
     5     2.103173e+03     1.320946e+02
 * time: 0.7238888740539551
     6     2.091136e+03     1.103035e+02
 * time: 0.7441959381103516
     7     2.081443e+03     1.091133e+02
 * time: 0.7636978626251221
     8     2.071793e+03     7.197675e+01
 * time: 0.7837240695953369
     9     2.062706e+03     7.623310e+01
 * time: 0.8381810188293457
    10     2.057515e+03     6.885476e+01
 * time: 0.8587050437927246
    11     2.051133e+03     6.368504e+01
 * time: 0.8782858848571777
    12     2.038626e+03     7.730243e+01
 * time: 0.8979129791259766
    13     2.019352e+03     1.136864e+02
 * time: 0.9178669452667236
    14     1.997136e+03     1.005748e+02
 * time: 0.9378490447998047
    15     1.983023e+03     6.831478e+01
 * time: 0.9709649085998535
    16     1.977700e+03     5.649783e+01
 * time: 0.9914679527282715
    17     1.974583e+03     6.357642e+01
 * time: 1.0115289688110352
    18     1.967292e+03     7.658974e+01
 * time: 1.032456874847412
    19     1.951045e+03     6.130573e+01
 * time: 1.0556640625
    20     1.935868e+03     4.845839e+01
 * time: 1.0878279209136963
    21     1.929356e+03     6.325782e+01
 * time: 1.1098530292510986
    22     1.925187e+03     3.142245e+01
 * time: 1.1301639080047607
    23     1.923733e+03     4.623400e+01
 * time: 1.1500270366668701
    24     1.918498e+03     5.347738e+01
 * time: 1.1806020736694336
    25     1.912383e+03     5.849125e+01
 * time: 1.2030160427093506
    26     1.905510e+03     3.254038e+01
 * time: 1.2247228622436523
    27     1.903629e+03     2.905618e+01
 * time: 1.2547860145568848
    28     1.902833e+03     2.907696e+01
 * time: 1.2750270366668701
    29     1.902447e+03     2.746037e+01
 * time: 1.2941889762878418
    30     1.899399e+03     1.930949e+01
 * time: 1.3145909309387207
    31     1.898705e+03     1.186361e+01
 * time: 1.3447389602661133
    32     1.898505e+03     1.050402e+01
 * time: 1.3646018505096436
    33     1.898474e+03     1.042186e+01
 * time: 1.382849931716919
    34     1.897992e+03     1.238729e+01
 * time: 1.4015629291534424
    35     1.897498e+03     1.729368e+01
 * time: 1.430527925491333
    36     1.896954e+03     1.472554e+01
 * time: 1.4500210285186768
    37     1.896744e+03     5.852824e+00
 * time: 1.4691429138183594
    38     1.896705e+03     1.171353e+00
 * time: 1.4868528842926025
    39     1.896704e+03     1.216117e+00
 * time: 1.515204906463623
    40     1.896703e+03     1.230336e+00
 * time: 1.5335350036621094
    41     1.896698e+03     1.250715e+00
 * time: 1.5519840717315674
    42     1.896688e+03     1.449552e+00
 * time: 1.5702450275421143
    43     1.896666e+03     2.533300e+00
 * time: 1.5988290309906006
    44     1.896631e+03     3.075536e+00
 * time: 1.6176040172576904
    45     1.896599e+03     2.139797e+00
 * time: 1.6363379955291748
    46     1.896587e+03     6.636030e-01
 * time: 1.6550359725952148
    47     1.896585e+03     6.303517e-01
 * time: 1.68361496925354
    48     1.896585e+03     5.995265e-01
 * time: 1.7016799449920654
    49     1.896584e+03     5.844159e-01
 * time: 1.7196369171142578
    50     1.896583e+03     6.083858e-01
 * time: 1.737612009048462
    51     1.896579e+03     8.145327e-01
 * time: 1.7656750679016113
    52     1.896570e+03     1.293490e+00
 * time: 1.7845509052276611
    53     1.896549e+03     1.877705e+00
 * time: 1.8031020164489746
    54     1.896513e+03     2.217392e+00
 * time: 1.8215129375457764
    55     1.896482e+03     1.658148e+00
 * time: 1.8497788906097412
    56     1.896466e+03     5.207217e-01
 * time: 1.8683269023895264
    57     1.896463e+03     1.177468e-01
 * time: 1.8866899013519287
    58     1.896463e+03     1.678937e-02
 * time: 1.9039969444274902
    59     1.896463e+03     2.666635e-03
 * time: 1.930351972579956
    60     1.896463e+03     2.666635e-03
 * time: 1.9652109146118164
    61     1.896463e+03     2.666635e-03
 * time: 1.9920899868011475

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 ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
    end
end
PumasModel
  Parameters: tvvc, tvq, tvvp, tvcl_hep, tvcl_ren, Ω, σ_add, σ_prop, dCRCL
  Random effects: η
  Covariates: WT, CRCL
  Dynamical system variables: Central, Peripheral
  Dynamical system type: Closed form
  Derived: DV
  Observed: DV

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

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

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

iparams_covariate_model_wt_crcl = (;
    tvvc = 5,
    tvcl_hep = 0.01,
    tvcl_ren = 0.01,
    tvq = 0.01,
    tvvp = 10,
    Ω = Diagonal([0.01, 0.01]),
    σ_add = 0.1,
    σ_prop = 0.1,
    dCRCL = 0.9,
)
(tvvc = 5,
 tvcl_hep = 0.01,
 tvcl_ren = 0.01,
 tvq = 0.01,
 tvvp = 10,
 Ω = [0.01 0.0; 0.0 0.01],
 σ_add = 0.1,
 σ_prop = 0.1,
 dCRCL = 0.9,)
fit_covariate_model_wt_crcl =
    fit(covariate_model_wt_crcl, pop, iparams_covariate_model_wt_crcl, FOCE())
[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter     Function value   Gradient norm 
     0     8.554152e+03     5.650279e+03
 * time: 3.695487976074219e-5
     1     3.572050e+03     1.302046e+03
 * time: 0.7907719612121582
     2     3.266947e+03     5.384244e+02
 * time: 1.0262470245361328
     3     3.150635e+03     1.918079e+02
 * time: 1.055945873260498
     4     3.108122e+03     1.277799e+02
 * time: 1.083669900894165
     5     3.091358e+03     8.838080e+01
 * time: 1.1098239421844482
     6     3.082997e+03     8.321689e+01
 * time: 1.1364469528198242
     7     3.076379e+03     8.167702e+01
 * time: 1.1626908779144287
     8     3.067422e+03     1.177822e+02
 * time: 1.1903748512268066
     9     3.048580e+03     2.526969e+02
 * time: 1.2786378860473633
    10     2.993096e+03     6.325396e+02
 * time: 1.3217649459838867
    11     2.965744e+03     7.634718e+02
 * time: 1.3999059200286865
    12     2.921212e+03     9.704020e+02
 * time: 1.4477999210357666
    13     2.553649e+03     6.642510e+02
 * time: 1.504417896270752
    14     2.319495e+03     3.204711e+02
 * time: 1.5618150234222412
    15     2.183040e+03     2.174905e+02
 * time: 1.606086015701294
    16     2.157621e+03     3.150983e+02
 * time: 1.633354902267456
    17     2.132395e+03     2.847614e+02
 * time: 1.67510986328125
    18     2.084799e+03     1.563370e+02
 * time: 1.7030329704284668
    19     2.071497e+03     1.006429e+02
 * time: 1.7289249897003174
    20     2.064983e+03     9.753313e+01
 * time: 1.754338026046753
    21     2.048289e+03     9.230405e+01
 * time: 1.7929918766021729
    22     2.020938e+03     1.292359e+02
 * time: 1.8195610046386719
    23     1.983888e+03     2.284990e+02
 * time: 1.8468818664550781
    24     1.962132e+03     1.220188e+02
 * time: 1.87260103225708
    25     1.945947e+03     1.035894e+02
 * time: 1.9090468883514404
    26     1.917782e+03     8.246698e+01
 * time: 1.9356598854064941
    27     1.905967e+03     5.408054e+01
 * time: 1.9621219635009766
    28     1.898569e+03     2.172222e+01
 * time: 1.9993810653686523
    29     1.897473e+03     1.689350e+01
 * time: 2.0248069763183594
    30     1.897019e+03     1.666689e+01
 * time: 2.049522876739502
    31     1.896796e+03     1.699751e+01
 * time: 2.0847270488739014
    32     1.896559e+03     1.645900e+01
 * time: 2.110153913497925
    33     1.896223e+03     1.415504e+01
 * time: 2.134842872619629
    34     1.895773e+03     1.630081e+01
 * time: 2.16812801361084
    35     1.895309e+03     1.723930e+01
 * time: 2.191890001296997
    36     1.895004e+03     1.229983e+01
 * time: 2.2155280113220215
    37     1.894871e+03     5.385102e+00
 * time: 2.2386488914489746
    38     1.894827e+03     3.465463e+00
 * time: 2.272165060043335
    39     1.894816e+03     3.387474e+00
 * time: 2.2950658798217773
    40     1.894807e+03     3.295388e+00
 * time: 2.3171470165252686
    41     1.894786e+03     3.089194e+00
 * time: 2.3496038913726807
    42     1.894737e+03     2.928080e+00
 * time: 2.371971845626831
    43     1.894624e+03     3.088723e+00
 * time: 2.3944320678710938
    44     1.894413e+03     3.493791e+00
 * time: 2.4268798828125
    45     1.894127e+03     3.142865e+00
 * time: 2.4501750469207764
    46     1.893933e+03     2.145253e+00
 * time: 2.4729180335998535
    47     1.893888e+03     2.172800e+00
 * time: 2.495419979095459
    48     1.893880e+03     2.180509e+00
 * time: 2.528125047683716
    49     1.893876e+03     2.134257e+00
 * time: 2.5501489639282227
    50     1.893868e+03     2.032354e+00
 * time: 2.572143077850342
    51     1.893846e+03     1.760874e+00
 * time: 2.604569911956787
    52     1.893796e+03     1.779016e+00
 * time: 2.6273109912872314
    53     1.893694e+03     2.018956e+00
 * time: 2.6500909328460693
    54     1.893559e+03     2.366854e+00
 * time: 2.6823549270629883
    55     1.893474e+03     3.690142e+00
 * time: 2.7051680088043213
    56     1.893446e+03     3.675109e+00
 * time: 2.727540969848633
    57     1.893439e+03     3.426419e+00
 * time: 2.7487540245056152
    58     1.893429e+03     3.183164e+00
 * time: 2.7806448936462402
    59     1.893398e+03     2.695171e+00
 * time: 2.8034799098968506
    60     1.893328e+03     2.753548e+00
 * time: 2.826385974884033
    61     1.893169e+03     3.589748e+00
 * time: 2.8604068756103516
    62     1.892920e+03     3.680718e+00
 * time: 2.883957862854004
    63     1.892667e+03     2.568107e+00
 * time: 2.906899929046631
    64     1.892514e+03     1.087910e+00
 * time: 2.9419338703155518
    65     1.892493e+03     3.287296e-01
 * time: 2.966784954071045
    66     1.892492e+03     2.967465e-01
 * time: 2.989466905593872
    67     1.892492e+03     3.020682e-01
 * time: 3.011033058166504
    68     1.892491e+03     3.034704e-01
 * time: 3.0428829193115234
    69     1.892491e+03     3.091846e-01
 * time: 3.0647499561309814
    70     1.892491e+03     3.224170e-01
 * time: 3.0860040187835693
    71     1.892490e+03     6.494197e-01
 * time: 3.1179280281066895
    72     1.892488e+03     1.115188e+00
 * time: 3.1401939392089844
    73     1.892483e+03     1.838833e+00
 * time: 3.162276029586792
    74     1.892472e+03     2.765371e+00
 * time: 3.184556007385254
    75     1.892452e+03     3.463807e+00
 * time: 3.217589855194092
    76     1.892431e+03     2.805270e+00
 * time: 3.240056037902832
    77     1.892411e+03     5.758916e-01
 * time: 3.262515068054199
    78     1.892410e+03     1.434041e-01
 * time: 3.295076847076416
    79     1.892409e+03     1.639246e-01
 * time: 3.3179469108581543
    80     1.892409e+03     1.145856e-01
 * time: 3.3397669792175293
    81     1.892409e+03     3.966861e-02
 * time: 3.370898962020874
    82     1.892409e+03     3.550808e-02
 * time: 3.3920600414276123
    83     1.892409e+03     3.456241e-02
 * time: 3.412559986114502
    84     1.892409e+03     3.114018e-02
 * time: 3.433134078979492
    85     1.892409e+03     4.080806e-02
 * time: 3.464169979095459
    86     1.892409e+03     6.722726e-02
 * time: 3.485511064529419
    87     1.892409e+03     1.006791e-01
 * time: 3.506422996520996
    88     1.892409e+03     1.303988e-01
 * time: 3.5376529693603516
    89     1.892409e+03     1.228919e-01
 * time: 3.5587570667266846
    90     1.892409e+03     6.433813e-02
 * time: 3.57991099357605
    91     1.892409e+03     1.314164e-02
 * time: 3.6008970737457275
    92     1.892409e+03     4.929931e-04
 * time: 3.631732940673828

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 ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
    end
end
PumasModel
  Parameters: tvvc, tvq, tvvp, tvcl_hep_M, tvcl_hep_F, tvcl_ren_M, tvcl_ren_F, Ω, σ_add, σ_prop, dCRCL_M, dCRCL_F
  Random effects: η
  Covariates: WT, CRCL, SEX
  Dynamical system variables: Central, Peripheral
  Dynamical system type: Closed form
  Derived: DV
  Observed: DV

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

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

iparams_covariate_model_wt_crcl_sex = (;
    tvvc = 5,
    tvcl_hep_M = 0.01,
    tvcl_hep_F = 0.01,
    tvcl_ren_M = 0.01,
    tvcl_ren_F = 0.01,
    tvq = 0.01,
    tvvp = 10,
    Ω = Diagonal([0.01, 0.01]),
    σ_add = 0.1,
    σ_prop = 0.1,
    dCRCL_M = 0.9,
    dCRCL_F = 0.9,
)
(tvvc = 5,
 tvcl_hep_M = 0.01,
 tvcl_hep_F = 0.01,
 tvcl_ren_M = 0.01,
 tvcl_ren_F = 0.01,
 tvq = 0.01,
 tvvp = 10,
 Ω = [0.01 0.0; 0.0 0.01],
 σ_add = 0.1,
 σ_prop = 0.1,
 dCRCL_M = 0.9,
 dCRCL_F = 0.9,)
fit_covariate_model_wt_crcl_sex =
    fit(covariate_model_wt_crcl_sex, pop, iparams_covariate_model_wt_crcl_sex, FOCE())
[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter     Function value   Gradient norm 
     0     8.554152e+03     5.650279e+03
 * time: 1.5020370483398438e-5
     1     3.641387e+03     1.432080e+03
 * time: 0.5748898983001709
     2     3.290450e+03     5.274782e+02
 * time: 0.6103770732879639
     3     3.185512e+03     2.173676e+02
 * time: 0.6424679756164551
     4     3.143009e+03     1.479653e+02
 * time: 0.8500308990478516
     5     3.128511e+03     8.980031e+01
 * time: 0.8791029453277588
     6     3.123188e+03     5.033125e+01
 * time: 0.9088430404663086
     7     3.120794e+03     4.279722e+01
 * time: 0.9379138946533203
     8     3.118627e+03     3.971051e+01
 * time: 0.9672050476074219
     9     3.115300e+03     8.456587e+01
 * time: 0.9967479705810547
    10     3.109353e+03     1.350354e+02
 * time: 1.0273690223693848
    11     3.095894e+03     1.998258e+02
 * time: 1.4164199829101562
    12     2.988214e+03     4.366433e+02
 * time: 1.4634170532226562
    13     2.896081e+03     5.505943e+02
 * time: 1.568173885345459
    14     2.652467e+03     7.300323e+02
 * time: 3.3959848880767822
    15     2.560937e+03     6.973661e+02
 * time: 3.538709878921509
    16     2.254941e+03     2.740033e+02
 * time: 3.620448112487793
    17     2.222509e+03     2.034303e+02
 * time: 3.6524150371551514
    18     2.171255e+03     2.449580e+02
 * time: 3.6836330890655518
    19     2.024532e+03     1.121511e+02
 * time: 3.713408946990967
    20     1.993723e+03     1.042814e+02
 * time: 3.741520881652832
    21     1.985113e+03     8.079014e+01
 * time: 3.7687458992004395
    22     1.976757e+03     7.054196e+01
 * time: 3.7956950664520264
    23     1.969970e+03     6.070322e+01
 * time: 3.8226358890533447
    24     1.961095e+03     6.810782e+01
 * time: 3.8528048992156982
    25     1.947983e+03     8.116920e+01
 * time: 3.8855741024017334
    26     1.930371e+03     8.530051e+01
 * time: 3.9367198944091797
    27     1.910209e+03     6.993170e+01
 * time: 3.965648889541626
    28     1.899107e+03     3.362640e+01
 * time: 3.9943079948425293
    29     1.898022e+03     2.642220e+01
 * time: 4.021176099777222
    30     1.897055e+03     1.213144e+01
 * time: 4.048240900039673
    31     1.896596e+03     7.773239e+00
 * time: 4.075054883956909
    32     1.896538e+03     7.997039e+00
 * time: 4.101284027099609
    33     1.896451e+03     8.160909e+00
 * time: 4.12702488899231
    34     1.896283e+03     8.237721e+00
 * time: 4.152663946151733
    35     1.895903e+03     1.520219e+01
 * time: 4.195198059082031
    36     1.895272e+03     2.358916e+01
 * time: 4.222898960113525
    37     1.894536e+03     2.461296e+01
 * time: 4.251069068908691
    38     1.893995e+03     1.546128e+01
 * time: 4.278095006942749
    39     1.893858e+03     6.976137e+00
 * time: 4.304495096206665
    40     1.893833e+03     6.019466e+00
 * time: 4.329987049102783
    41     1.893786e+03     3.827201e+00
 * time: 4.355616092681885
    42     1.893714e+03     3.323412e+00
 * time: 4.397566080093384
    43     1.893592e+03     3.215150e+00
 * time: 4.424911975860596
    44     1.893435e+03     6.534965e+00
 * time: 4.451607942581177
    45     1.893286e+03     7.424154e+00
 * time: 4.478426933288574
    46     1.893190e+03     5.552627e+00
 * time: 4.504862070083618
    47     1.893139e+03     3.222316e+00
 * time: 4.543848991394043
    48     1.893120e+03     3.015339e+00
 * time: 4.569974899291992
    49     1.893107e+03     3.244809e+00
 * time: 4.596102952957153
    50     1.893080e+03     6.163100e+00
 * time: 4.621737957000732
    51     1.893027e+03     9.824713e+00
 * time: 4.648016929626465
    52     1.892912e+03     1.390100e+01
 * time: 4.68916392326355
    53     1.892734e+03     1.510937e+01
 * time: 4.717516899108887
    54     1.892561e+03     1.008563e+01
 * time: 4.744771957397461
    55     1.892485e+03     3.730668e+00
 * time: 4.771629095077515
    56     1.892471e+03     3.380261e+00
 * time: 4.808803081512451
    57     1.892463e+03     3.167904e+00
 * time: 4.834775924682617
    58     1.892441e+03     4.152065e+00
 * time: 4.860971927642822
    59     1.892391e+03     7.355996e+00
 * time: 4.903656005859375
    60     1.892268e+03     1.195397e+01
 * time: 4.945604085922241
    61     1.892026e+03     1.640783e+01
 * time: 4.999799013137817
    62     1.891735e+03     1.593576e+01
 * time: 5.032012939453125
    63     1.891569e+03     8.316423e+00
 * time: 5.057825088500977
    64     1.891494e+03     3.948212e+00
 * time: 5.094204902648926
    65     1.891481e+03     3.911593e+00
 * time: 5.121423959732056
    66     1.891457e+03     3.875559e+00
 * time: 5.148003101348877
    67     1.891405e+03     3.811247e+00
 * time: 5.184107065200806
    68     1.891262e+03     3.657045e+00
 * time: 5.210493087768555
    69     1.890930e+03     4.957405e+00
 * time: 5.246721982955933
    70     1.890317e+03     6.657726e+00
 * time: 5.273772954940796
    71     1.889660e+03     6.086302e+00
 * time: 5.3006439208984375
    72     1.889303e+03     2.270929e+00
 * time: 5.337035894393921
    73     1.889253e+03     7.695301e-01
 * time: 5.363384962081909
    74     1.889252e+03     7.382144e-01
 * time: 5.388885021209717
    75     1.889251e+03     7.187898e-01
 * time: 5.4250168800354
    76     1.889251e+03     7.215047e-01
 * time: 5.449865102767944
    77     1.889250e+03     7.235155e-01
 * time: 5.4845709800720215
    78     1.889249e+03     7.246818e-01
 * time: 5.510447978973389
    79     1.889244e+03     7.257796e-01
 * time: 5.536461114883423
    80     1.889233e+03     7.198190e-01
 * time: 5.5738489627838135
    81     1.889204e+03     1.089029e+00
 * time: 5.6005189418792725
    82     1.889142e+03     1.801601e+00
 * time: 5.636109113693237
    83     1.889043e+03     2.967917e+00
 * time: 5.662697076797485
    84     1.888889e+03     2.965856e+00
 * time: 5.689665079116821
    85     1.888705e+03     5.933555e-01
 * time: 5.726285934448242
    86     1.888655e+03     9.577698e-01
 * time: 5.752593040466309
    87     1.888582e+03     1.498494e+00
 * time: 5.778956890106201
    88     1.888533e+03     1.502751e+00
 * time: 5.815237045288086
    89     1.888490e+03     1.184664e+00
 * time: 5.842037916183472
    90     1.888480e+03     6.684515e-01
 * time: 5.8782970905303955
    91     1.888476e+03     3.680032e-01
 * time: 5.904489040374756
    92     1.888476e+03     4.720040e-01
 * time: 5.930001974105835
    93     1.888476e+03     4.768646e-01
 * time: 5.96504807472229
    94     1.888475e+03     4.736674e-01
 * time: 5.990714073181152
    95     1.888475e+03     4.552765e-01
 * time: 6.016721963882446
    96     1.888474e+03     5.193725e-01
 * time: 6.052011013031006
    97     1.888473e+03     8.850098e-01
 * time: 6.077670097351074
    98     1.888468e+03     1.461598e+00
 * time: 6.112586975097656
    99     1.888458e+03     2.209124e+00
 * time: 6.138339996337891
   100     1.888437e+03     2.961236e+00
 * time: 6.164604902267456
   101     1.888407e+03     2.978462e+00
 * time: 6.200401067733765
   102     1.888384e+03     1.707199e+00
 * time: 6.226569890975952
   103     1.888381e+03     6.198984e-01
 * time: 6.252886056900024
   104     1.888380e+03     5.171082e-01
 * time: 6.288728952407837
   105     1.888378e+03     1.037321e-01
 * time: 6.3151938915252686
   106     1.888378e+03     8.473253e-02
 * time: 6.3506550788879395
   107     1.888378e+03     8.364965e-02
 * time: 6.376274108886719
   108     1.888378e+03     8.080442e-02
 * time: 6.402374029159546
   109     1.888378e+03     7.873900e-02
 * time: 6.437762975692749
   110     1.888378e+03     7.798398e-02
 * time: 6.463335990905762
   111     1.888378e+03     7.788170e-02
 * time: 6.488385915756226
   112     1.888378e+03     7.776461e-02
 * time: 6.5234410762786865
   113     1.888378e+03     9.023586e-02
 * time: 6.548877954483032
   114     1.888378e+03     1.631370e-01
 * time: 6.584585905075073
   115     1.888378e+03     2.768691e-01
 * time: 6.610511064529419
   116     1.888377e+03     4.462313e-01
 * time: 6.636148929595947
   117     1.888377e+03     6.643167e-01
 * time: 6.671735048294067
   118     1.888375e+03     8.433175e-01
 * time: 6.697432994842529
   119     1.888374e+03     7.596473e-01
 * time: 6.7236809730529785
   120     1.888373e+03     3.637851e-01
 * time: 6.760334014892578
   121     1.888372e+03     8.305224e-02
 * time: 6.786334991455078
   122     1.888372e+03     2.084516e-02
 * time: 6.821507930755615
   123     1.888372e+03     2.056416e-02
 * time: 6.846122980117798
   124     1.888372e+03     2.044079e-02
 * time: 6.869066953659058
   125     1.888372e+03     2.035197e-02
 * time: 6.900161027908325
   126     1.888372e+03     2.021266e-02
 * time: 6.92298698425293
   127     1.888372e+03     1.998168e-02
 * time: 6.945894956588745
   128     1.888372e+03     3.162094e-02
 * time: 6.9776930809021
   129     1.888372e+03     5.510007e-02
 * time: 7.001025915145874
   130     1.888372e+03     9.277177e-02
 * time: 7.024463891983032
   131     1.888372e+03     1.528967e-01
 * time: 7.057027101516724
   132     1.888372e+03     2.462112e-01
 * time: 7.080802917480469
   133     1.888372e+03     3.799880e-01
 * time: 7.113508939743042
   134     1.888371e+03     5.312357e-01
 * time: 7.137594938278198
   135     1.888369e+03     6.019766e-01
 * time: 7.1617701053619385
   136     1.888367e+03     4.665348e-01
 * time: 7.194713115692139
   137     1.888366e+03     1.404917e-01
 * time: 7.218880891799927
   138     1.888365e+03     8.513280e-02
 * time: 7.24288010597229
   139     1.888364e+03     1.244285e-01
 * time: 7.275374889373779
   140     1.888364e+03     1.028231e-01
 * time: 7.298995018005371
   141     1.888364e+03     5.163326e-02
 * time: 7.331707954406738
   142     1.888364e+03     5.148616e-02
 * time: 7.3558509349823
   143     1.888364e+03     3.147247e-02
 * time: 7.3795740604400635
   144     1.888364e+03     2.104597e-02
 * time: 7.412302017211914
   145     1.888364e+03     6.541596e-03
 * time: 7.435370922088623
   146     1.888364e+03     2.538502e-03
 * time: 7.460169076919556
   147     1.888364e+03     4.361857e-03
 * time: 7.491796970367432
   148     1.888364e+03     3.034845e-03
 * time: 7.514419078826904
   149     1.888364e+03     5.961460e-04
 * time: 7.545598030090332

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 3.353
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 1.992
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 3.632
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.546
3 Subjects 30
4 Fixed Parameters 0
5 Optimized Parameters 13
6 DV Active Observations 540
7 DV Missing Observations 0
8 Total Active Observations 540
9 Total Missing Observations 0
10 Likelihood Approximation Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}
11 LogLikelihood (LL) -1888.36
12 -2LL 3776.73
13 AIC 3802.73
14 BIC 3858.52
15 (η-shrinkage) η₁ -0.013
16 (η-shrinkage) η₂ -0.012
17 (ϵ-shrinkage) DV 0.056

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

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

all_metrics = innerjoin(
    metrics_table(fit_base_model),
    metrics_table(fit_covariate_model_wt),
    metrics_table(fit_covariate_model_wt_crcl),
    metrics_table(fit_covariate_model_wt_crcl_sex);
    on = :Metric,
    makeunique = true,
);
rename!(
    all_metrics,
    :Value => :Base_Model,
    :Value_1 => :Covariate_Model_WT,
    :Value_2 => :Covariate_Model_WT_CRCL,
    :Value_3 => :Covariate_Model_WT_CRCL_SEX,
)
17×5 DataFrame
Row Metric Base_Model Covariate_Model_WT Covariate_Model_WT_CRCL Covariate_Model_WT_CRCL_SEX
String Any Any Any Any
1 Successful true true true true
2 Estimation Time 3.353 1.992 3.632 7.546
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.6299849528186
aic(fit_covariate_model_wt)
3808.9264607805967
aic(fit_covariate_model_wt_crcl)
3804.8179473717055
aic(fit_covariate_model_wt_crcl_sex)
3802.7275243740673

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.002716064453125e-5
     1     2.994747e+03     1.083462e+03
 * time: 0.533405065536499
     2     2.406265e+03     1.884408e+02
 * time: 0.5377020835876465
     3     2.344175e+03     7.741610e+01
 * time: 0.5419151782989502
     4     2.323153e+03     2.907642e+01
 * time: 0.5461940765380859
     5     2.318222e+03     2.273295e+01
 * time: 0.5503830909729004
     6     2.316833e+03     1.390527e+01
 * time: 0.5543379783630371
     7     2.316425e+03     4.490883e+00
 * time: 0.5582060813903809
     8     2.316362e+03     9.374519e-01
 * time: 0.5620231628417969
     9     2.316356e+03     1.928785e-01
 * time: 0.5658309459686279
    10     2.316355e+03     3.119615e-02
 * time: 0.5696070194244385
    11     2.316355e+03     6.215513e-03
 * time: 0.5733990669250488
    12     2.316355e+03     8.313206e-04
 * time: 0.5773200988769531

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

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

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

Hence, for LOCF, you can use the following:

pop = read_pumas(pkdata; covariates_direction = :left)

along with any other required keyword arguments for column mapping.

Note

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

8 Categorical Covariates

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

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

For example:

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

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

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

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

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

9 Conclusion

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

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

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