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:

include covariates in your model
parse data with covariates
understand the difference between constant and time-varying covariates
handle continuous and categorical covariates
deal with missing data in your covariates
deal with categorical covariates

1 `nlme_sample` Dataset

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

pkfile = dataset("nlme_sample", String)
pkdata = CSV.read(pkfile, DataFrame; missingstring = ["NA", ""])
first(pkdata, 5)

5×15 DataFrame

Row	ID	TIME	DV	AMT	EVID	CMT	RATE	WT	AGE	SEX	CRCL	GROUP	ROUTE	DURATION	OCC
	Int64	Float64	Float64?	Int64?	Int64	Int64?	Int64	Int64	Int64	String1	Int64	String7	Float64	Int64?	Int64
1	1	0.0	missing	1000	1	1	500	90	47	M	75	1000 mg	Inf	2	1
2	1	0.001	0.0667231	missing	0	missing	0	90	47	M	75	1000 mg	Inf	missing	1
3	1	1.0	112.817	missing	0	missing	0	90	47	M	75	1000 mg	Inf	missing	1
4	1	2.0	224.087	missing	0	missing	0	90	47	M	75	1000 mg	Inf	missing	1
5	1	4.0	220.047	missing	0	missing	0	90	47	M	75	1000 mg	Inf	missing	1

Note

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

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

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

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

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

describe(pkdata, :mean, :std, :nunique, :first, :last, :eltype)

15×7 DataFrame

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

Tip

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

2 Step 1 - Parse Data into a `Population`

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

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

Population
  Subjects: 30
  Covariates: WT, AGE, SEX, CRCL, GROUP
  Observations: DV

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

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

3 Step 2 - Base Model

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

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

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

    @random begin
        η ~ MvNormal(Ω)
    end

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

    @dynamics Central1Periph1

    @derived begin
        cp := @. Central / Vc
        DV ~ @. CombinedNormal(cp, σ_add, σ_prop)
    end
end

PumasModel
  Parameters: tvcl, tvvc, tvq, tvvp, Ω, σ_add, σ_prop
  Random effects: η
  Covariates:
  Dynamical system variables: Central, Peripheral
  Dynamical system type: Closed form
  Derived: DV
  Observed: DV

And fit it accordingly:

iparams_base_model = (;
    tvvc = 5,
    tvcl = 0.02,
    tvq = 0.01,
    tvvp = 10,
    Ω = Diagonal([0.01, 0.01]),
    σ_add = 0.1,
    σ_prop = 0.1,
)

(tvvc = 5, tvcl = 0.02, tvq = 0.01, tvvp = 10, Ω = [0.01 0.0; 0.0 0.01], σ_add = 0.1, σ_prop = 0.1)

fit_base_model = fit(base_model, pop, iparams_base_model, FOCE())

[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter     Function value   Gradient norm 
     0     8.300164e+03     4.360770e+03
 * time: 0.04153299331665039
     1     3.110315e+03     9.706222e+02
 * time: 3.5447771549224854
     2     2.831659e+03     7.817006e+02
 * time: 3.5882480144500732
     3     2.405281e+03     2.923696e+02
 * time: 3.6240391731262207
     4     2.370406e+03     3.032286e+02
 * time: 3.6508941650390625
     5     2.313631e+03     3.126188e+02
 * time: 3.6782450675964355
     6     2.263986e+03     2.946697e+02
 * time: 3.7031450271606445
     7     2.160182e+03     1.917599e+02
 * time: 3.851119041442871
     8     2.112467e+03     1.588288e+02
 * time: 3.8982291221618652
     9     2.090339e+03     1.108334e+02
 * time: 3.921316146850586
    10     2.078171e+03     8.108027e+01
 * time: 3.9444730281829834
    11     2.074517e+03     7.813928e+01
 * time: 3.9658491611480713
    12     2.066270e+03     7.044632e+01
 * time: 3.9873311519622803
    13     2.049660e+03     1.062584e+02
 * time: 4.008433103561401
    14     2.021965e+03     1.130570e+02
 * time: 4.030004978179932
    15     1.994936e+03     7.825801e+01
 * time: 4.102138042449951
    16     1.979337e+03     5.318263e+01
 * time: 4.125488042831421
    17     1.972141e+03     6.807046e+01
 * time: 4.148187160491943
    18     1.967973e+03     7.896361e+01
 * time: 4.1705591678619385
    19     1.962237e+03     8.343757e+01
 * time: 4.192453145980835
    20     1.952791e+03     5.565304e+01
 * time: 4.21525502204895
    21     1.935857e+03     3.923284e+01
 * time: 4.238150119781494
    22     1.926254e+03     5.749643e+01
 * time: 4.275330066680908
    23     1.922144e+03     4.306225e+01
 * time: 4.297933101654053
    24     1.911566e+03     4.810496e+01
 * time: 4.321612119674683
    25     1.906464e+03     4.324267e+01
 * time: 4.344896078109741
    26     1.905339e+03     1.207954e+01
 * time: 4.3661229610443115
    27     1.905092e+03     1.093479e+01
 * time: 4.38704514503479
    28     1.904957e+03     1.057034e+01
 * time: 4.408066034317017
    29     1.904875e+03     1.060882e+01
 * time: 4.440562963485718
    30     1.904459e+03     1.031525e+01
 * time: 4.462398052215576
    31     1.903886e+03     1.041179e+01
 * time: 4.484292030334473
    32     1.903313e+03     1.135672e+01
 * time: 4.505895137786865
    33     1.903057e+03     1.075683e+01
 * time: 4.526810169219971
    34     1.902950e+03     1.091295e+01
 * time: 4.558855056762695
    35     1.902887e+03     1.042409e+01
 * time: 4.580183029174805
    36     1.902640e+03     9.213043e+00
 * time: 4.601253032684326
    37     1.902364e+03     9.519321e+00
 * time: 4.622502088546753
    38     1.902156e+03     5.590984e+00
 * time: 4.64333701133728
    39     1.902094e+03     1.811898e+00
 * time: 4.674285173416138
    40     1.902086e+03     1.644722e+00
 * time: 4.695549011230469
    41     1.902084e+03     1.653520e+00
 * time: 4.715636968612671
    42     1.902074e+03     1.720184e+00
 * time: 4.746025085449219
    43     1.902055e+03     2.619061e+00
 * time: 4.766674995422363
    44     1.902015e+03     3.519729e+00
 * time: 4.7875800132751465
    45     1.901962e+03     3.403372e+00
 * time: 4.808030128479004
    46     1.901924e+03     1.945644e+00
 * time: 4.838848114013672
    47     1.901914e+03     6.273342e-01
 * time: 4.859898090362549
    48     1.901913e+03     5.374557e-01
 * time: 4.880105972290039
    49     1.901913e+03     5.710007e-01
 * time: 4.899749994277954
    50     1.901913e+03     6.091390e-01
 * time: 4.92940616607666
    51     1.901912e+03     6.692417e-01
 * time: 4.949562072753906
    52     1.901909e+03     7.579620e-01
 * time: 4.969468116760254
    53     1.901903e+03     8.798211e-01
 * time: 4.998855113983154
    54     1.901889e+03     1.002981e+00
 * time: 5.019103050231934
    55     1.901864e+03     1.495512e+00
 * time: 5.039591073989868
    56     1.901837e+03     1.664621e+00
 * time: 5.059282064437866
    57     1.901819e+03     8.601119e-01
 * time: 5.088977098464966
    58     1.901815e+03     4.525470e-02
 * time: 5.109663009643555
    59     1.901815e+03     1.294280e-02
 * time: 5.129894018173218
    60     1.901815e+03     2.876567e-03
 * time: 5.148813962936401
    61     1.901815e+03     2.876567e-03
 * time: 5.2155821323394775
    62     1.901815e+03     2.876523e-03
 * time: 5.271620988845825
    63     1.901815e+03     2.876512e-03
 * time: 5.3117899894714355
    64     1.901815e+03     2.876512e-03
 * time: 5.370476961135864
    65     1.901815e+03     2.876512e-03
 * time: 5.41736912727356

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:

allometric scaling based on weight
clearance effect based on creatinine clearance
separated error model based on sex

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

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

Tip

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

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

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

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

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates begin
        WT
    end

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

    @dynamics Central1Periph1

    @derived begin
        cp := @. Central / Vc
        DV ~ @. CombinedNormal(cp, σ_add, σ_prop)
    end
end

PumasModel
  Parameters: tvcl, tvvc, tvq, tvvp, Ω, σ_add, σ_prop
  Random effects: η
  Covariates: WT
  Dynamical system variables: Central, Peripheral
  Dynamical system type: Closed form
  Derived: DV
  Observed: DV

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

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

fit_covariate_model_wt = fit(covariate_model_wt, pop, iparams_base_model, FOCE())

[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter     Function value   Gradient norm 
     0     7.695401e+03     4.898919e+03
 * time: 1.4066696166992188e-5
     1     2.846050e+03     1.128657e+03
 * time: 1.1916320323944092
     2     2.472982e+03     7.008264e+02
 * time: 1.229189157485962
     3     2.180690e+03     2.312704e+02
 * time: 1.3424561023712158
     4     2.125801e+03     1.862929e+02
 * time: 1.3666880130767822
     5     2.103173e+03     1.320946e+02
 * time: 1.3896420001983643
     6     2.091136e+03     1.103035e+02
 * time: 1.4117610454559326
     7     2.081443e+03     1.091133e+02
 * time: 1.433312177658081
     8     2.071793e+03     7.197675e+01
 * time: 1.4557201862335205
     9     2.062706e+03     7.623310e+01
 * time: 1.4792242050170898
    10     2.057515e+03     6.885476e+01
 * time: 1.5395112037658691
    11     2.051133e+03     6.368504e+01
 * time: 1.5617871284484863
    12     2.038626e+03     7.730243e+01
 * time: 1.5835020542144775
    13     2.019352e+03     1.136864e+02
 * time: 1.6054801940917969
    14     1.997136e+03     1.005748e+02
 * time: 1.627824068069458
    15     1.983023e+03     6.831478e+01
 * time: 1.6508500576019287
    16     1.977700e+03     5.649783e+01
 * time: 1.6863172054290771
    17     1.974583e+03     6.357642e+01
 * time: 1.709043025970459
    18     1.967292e+03     7.658974e+01
 * time: 1.7325310707092285
    19     1.951045e+03     6.130573e+01
 * time: 1.7578961849212646
    20     1.935868e+03     4.845839e+01
 * time: 1.7810051441192627
    21     1.929356e+03     6.325782e+01
 * time: 1.8164491653442383
    22     1.925187e+03     3.142245e+01
 * time: 1.8393731117248535
    23     1.923733e+03     4.623400e+01
 * time: 1.861497163772583
    24     1.918498e+03     5.347738e+01
 * time: 1.894430160522461
    25     1.912383e+03     5.849125e+01
 * time: 1.9187970161437988
    26     1.905510e+03     3.254038e+01
 * time: 1.9427549839019775
    27     1.903629e+03     2.905618e+01
 * time: 1.9646670818328857
    28     1.902833e+03     2.907696e+01
 * time: 1.9973640441894531
    29     1.902447e+03     2.746037e+01
 * time: 2.018665075302124
    30     1.899399e+03     1.930949e+01
 * time: 2.0408260822296143
    31     1.898705e+03     1.186361e+01
 * time: 2.072859048843384
    32     1.898505e+03     1.050402e+01
 * time: 2.094892978668213
    33     1.898474e+03     1.042186e+01
 * time: 2.118039131164551
    34     1.897992e+03     1.238729e+01
 * time: 2.1403660774230957
    35     1.897498e+03     1.729368e+01
 * time: 2.173579216003418
    36     1.896954e+03     1.472554e+01
 * time: 2.1953961849212646
    37     1.896744e+03     5.852825e+00
 * time: 2.2164862155914307
    38     1.896705e+03     1.171353e+00
 * time: 2.235924005508423
    39     1.896704e+03     1.216117e+00
 * time: 2.266849994659424
    40     1.896703e+03     1.230336e+00
 * time: 2.287353038787842
    41     1.896698e+03     1.250715e+00
 * time: 2.3075220584869385
    42     1.896688e+03     1.449552e+00
 * time: 2.3376951217651367
    43     1.896666e+03     2.533300e+00
 * time: 2.358933210372925
    44     1.896631e+03     3.075537e+00
 * time: 2.3799471855163574
    45     1.896599e+03     2.139797e+00
 * time: 2.4002161026000977
    46     1.896587e+03     6.636030e-01
 * time: 2.43083119392395
    47     1.896585e+03     6.303517e-01
 * time: 2.451341152191162
    48     1.896585e+03     5.995265e-01
 * time: 2.470959186553955
    49     1.896584e+03     5.844159e-01
 * time: 2.4908339977264404
    50     1.896583e+03     6.083858e-01
 * time: 2.5206010341644287
    51     1.896579e+03     8.145327e-01
 * time: 2.540472984313965
    52     1.896570e+03     1.293490e+00
 * time: 2.560474157333374
    53     1.896549e+03     1.877705e+00
 * time: 2.5803580284118652
    54     1.896513e+03     2.217392e+00
 * time: 2.6100080013275146
    55     1.896482e+03     1.658148e+00
 * time: 2.6306521892547607
    56     1.896466e+03     5.207218e-01
 * time: 2.650658130645752
    57     1.896463e+03     1.177468e-01
 * time: 2.6805331707000732
    58     1.896463e+03     1.678937e-02
 * time: 2.699949026107788
    59     1.896463e+03     2.666636e-03
 * time: 2.718410015106201
    60     1.896463e+03     2.666636e-03
 * time: 2.777644157409668
    61     1.896463e+03     2.666636e-03
 * time: 2.8227171897888184

FittedPumasModel

Dynamical system type:                 Closed form

Number of subjects:                             30

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

Number of parameters:      Constant      Optimized
                                  0              8

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                      NoXChange
Log-likelihood value:                   -1896.4632

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

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

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

Tip

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

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

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

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

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates begin
        WT
        CRCL
    end

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

    @dynamics Central1Periph1

    @derived begin
        cp := @. Central / Vc
        DV ~ @. CombinedNormal(cp, σ_add, σ_prop)
    end
end

PumasModel
  Parameters: tvvc, tvq, tvvp, tvcl_hep, tvcl_ren, Ω, σ_add, σ_prop, dCRCL
  Random effects: η
  Covariates: WT, CRCL
  Dynamical system variables: Central, Peripheral
  Dynamical system type: Closed form
  Derived: DV
  Observed: DV

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

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

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

iparams_covariate_model_wt_crcl = (;
    tvvc = 5,
    tvcl_hep = 0.01,
    tvcl_ren = 0.01,
    tvq = 0.01,
    tvvp = 10,
    Ω = Diagonal([0.01, 0.01]),
    σ_add = 0.1,
    σ_prop = 0.1,
    dCRCL = 0.9,
)

(tvvc = 5, tvcl_hep = 0.01, tvcl_ren = 0.01, tvq = 0.01, tvvp = 10, Ω = [0.01 0.0; 0.0 0.01], σ_add = 0.1, σ_prop = 0.1, dCRCL = 0.9)

fit_covariate_model_wt_crcl =
    fit(covariate_model_wt_crcl, pop, iparams_covariate_model_wt_crcl, FOCE())

[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter     Function value   Gradient norm 
     0     8.554152e+03     5.650279e+03
 * time: 1.4066696166992188e-5
     1     3.572050e+03     1.302046e+03
 * time: 1.2454960346221924
     2     3.266947e+03     5.384244e+02
 * time: 1.2784450054168701
     3     3.150635e+03     1.918079e+02
 * time: 1.4240050315856934
     4     3.108122e+03     1.277799e+02
 * time: 1.4512109756469727
     5     3.091358e+03     8.838080e+01
 * time: 1.477060079574585
     6     3.082997e+03     8.321689e+01
 * time: 1.5026259422302246
     7     3.076379e+03     8.167702e+01
 * time: 1.5278339385986328
     8     3.067422e+03     1.177822e+02
 * time: 1.5547630786895752
     9     3.048580e+03     2.526969e+02
 * time: 1.6309540271759033
    10     2.993096e+03     6.325396e+02
 * time: 1.6724660396575928
    11     2.965744e+03     7.634718e+02
 * time: 1.7468268871307373
    12     2.921212e+03     9.704020e+02
 * time: 1.8082950115203857
    13     2.553649e+03     6.642510e+02
 * time: 1.8465840816497803
    14     2.319495e+03     3.204711e+02
 * time: 1.9021408557891846
    15     2.183040e+03     2.174905e+02
 * time: 1.9585418701171875
    16     2.157621e+03     3.150983e+02
 * time: 1.9865999221801758
    17     2.132395e+03     2.847614e+02
 * time: 2.0131380558013916
    18     2.084799e+03     1.563370e+02
 * time: 2.039280891418457
    19     2.071497e+03     1.006429e+02
 * time: 2.0776400566101074
    20     2.064983e+03     9.753313e+01
 * time: 2.1042048931121826
    21     2.048289e+03     9.230405e+01
 * time: 2.1305770874023438
    22     2.020938e+03     1.292359e+02
 * time: 2.1679649353027344
    23     1.983888e+03     2.284990e+02
 * time: 2.195605993270874
    24     1.962132e+03     1.220188e+02
 * time: 2.2216129302978516
    25     1.945947e+03     1.035894e+02
 * time: 2.2571699619293213
    26     1.917782e+03     8.246698e+01
 * time: 2.283421039581299
    27     1.905967e+03     5.408054e+01
 * time: 2.3098669052124023
    28     1.898569e+03     2.172222e+01
 * time: 2.3475189208984375
    29     1.897473e+03     1.689350e+01
 * time: 2.3730289936065674
    30     1.897019e+03     1.666689e+01
 * time: 2.408107042312622
    31     1.896796e+03     1.699751e+01
 * time: 2.433187961578369
    32     1.896559e+03     1.645900e+01
 * time: 2.4581730365753174
    33     1.896223e+03     1.415504e+01
 * time: 2.493090867996216
    34     1.895773e+03     1.630081e+01
 * time: 2.518735885620117
    35     1.895309e+03     1.723930e+01
 * time: 2.544666051864624
    36     1.895004e+03     1.229983e+01
 * time: 2.580965995788574
    37     1.894871e+03     5.385102e+00
 * time: 2.6072728633880615
    38     1.894827e+03     3.465463e+00
 * time: 2.6323418617248535
    39     1.894816e+03     3.387474e+00
 * time: 2.6564719676971436
    40     1.894807e+03     3.295388e+00
 * time: 2.692502021789551
    41     1.894786e+03     3.089194e+00
 * time: 2.7172179222106934
    42     1.894737e+03     2.928080e+00
 * time: 2.7415759563446045
    43     1.894624e+03     3.088723e+00
 * time: 2.777448892593384
    44     1.894413e+03     3.493791e+00
 * time: 2.802393913269043
    45     1.894127e+03     3.142865e+00
 * time: 2.8269240856170654
    46     1.893933e+03     2.145253e+00
 * time: 2.862438917160034
    47     1.893888e+03     2.172800e+00
 * time: 2.8876190185546875
    48     1.893880e+03     2.180509e+00
 * time: 2.912061929702759
    49     1.893876e+03     2.134257e+00
 * time: 2.9466419219970703
    50     1.893868e+03     2.032354e+00
 * time: 2.9708058834075928
    51     1.893846e+03     1.760874e+00
 * time: 2.9948790073394775
    52     1.893796e+03     1.779016e+00
 * time: 3.0293359756469727
    53     1.893694e+03     2.018956e+00
 * time: 3.0543808937072754
    54     1.893559e+03     2.366854e+00
 * time: 3.078300952911377
    55     1.893474e+03     3.690142e+00
 * time: 3.1137728691101074
    56     1.893446e+03     3.675109e+00
 * time: 3.1386849880218506
    57     1.893439e+03     3.426419e+00
 * time: 3.1625258922576904
    58     1.893429e+03     3.183164e+00
 * time: 3.196418046951294
    59     1.893398e+03     2.695171e+00
 * time: 3.2209630012512207
    60     1.893328e+03     2.753548e+00
 * time: 3.2451930046081543
    61     1.893169e+03     3.589748e+00
 * time: 3.280395984649658
    62     1.892920e+03     3.680718e+00
 * time: 3.3065290451049805
    63     1.892667e+03     2.568107e+00
 * time: 3.3311538696289062
    64     1.892514e+03     1.087910e+00
 * time: 3.365957021713257
    65     1.892493e+03     3.287296e-01
 * time: 3.390516996383667
    66     1.892492e+03     2.967465e-01
 * time: 3.4148240089416504
    67     1.892492e+03     3.020682e-01
 * time: 3.4489879608154297
    68     1.892491e+03     3.034704e-01
 * time: 3.4721899032592773
    69     1.892491e+03     3.091846e-01
 * time: 3.4953200817108154
    70     1.892491e+03     3.224170e-01
 * time: 3.528428077697754
    71     1.892490e+03     6.494197e-01
 * time: 3.5523200035095215
    72     1.892488e+03     1.115188e+00
 * time: 3.57588791847229
    73     1.892483e+03     1.838833e+00
 * time: 3.609438896179199
    74     1.892472e+03     2.765371e+00
 * time: 3.6339519023895264
    75     1.892452e+03     3.463807e+00
 * time: 3.6583268642425537
    76     1.892431e+03     2.805270e+00
 * time: 3.6923630237579346
    77     1.892411e+03     5.758916e-01
 * time: 3.7176458835601807
    78     1.892410e+03     1.434041e-01
 * time: 3.7419240474700928
    79     1.892409e+03     1.639246e-01
 * time: 3.775852918624878
    80     1.892409e+03     1.145856e-01
 * time: 3.799983024597168
    81     1.892409e+03     3.966861e-02
 * time: 3.823279857635498
    82     1.892409e+03     3.550808e-02
 * time: 3.857668876647949
    83     1.892409e+03     3.456241e-02
 * time: 3.880631923675537
    84     1.892409e+03     3.114018e-02
 * time: 3.903493881225586
    85     1.892409e+03     4.080806e-02
 * time: 3.936255931854248
    86     1.892409e+03     6.722726e-02
 * time: 3.9599790573120117
    87     1.892409e+03     1.006791e-01
 * time: 3.9833359718322754
    88     1.892409e+03     1.303988e-01
 * time: 4.006610870361328
    89     1.892409e+03     1.228919e-01
 * time: 4.04027795791626
    90     1.892409e+03     6.433813e-02
 * time: 4.063938856124878
    91     1.892409e+03     1.314164e-02
 * time: 4.087267875671387
    92     1.892409e+03     4.929931e-04
 * time: 4.120559930801392

FittedPumasModel

Dynamical system type:                 Closed form

Number of subjects:                             30

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

Number of parameters:      Constant      Optimized
                                  0             10

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                    -1892.409

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

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

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

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

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates begin
        WT
        CRCL
        SEX
    end

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

    @dynamics Central1Periph1

    @derived begin
        cp := @. Central / Vc
        DV ~ @. CombinedNormal(cp, σ_add, σ_prop)
    end
end

PumasModel
  Parameters: tvvc, tvq, tvvp, tvcl_hep_M, tvcl_hep_F, tvcl_ren_M, tvcl_ren_F, Ω, σ_add, σ_prop, dCRCL_M, dCRCL_F
  Random effects: η
  Covariates: WT, CRCL, SEX
  Dynamical system variables: Central, Peripheral
  Dynamical system type: Closed form
  Derived: DV
  Observed: DV

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

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

iparams_covariate_model_wt_crcl_sex = (;
    tvvc = 5,
    tvcl_hep_M = 0.01,
    tvcl_hep_F = 0.01,
    tvcl_ren_M = 0.01,
    tvcl_ren_F = 0.01,
    tvq = 0.01,
    tvvp = 10,
    Ω = Diagonal([0.01, 0.01]),
    σ_add = 0.1,
    σ_prop = 0.1,
    dCRCL_M = 0.9,
    dCRCL_F = 0.9,
)

(tvvc = 5, tvcl_hep_M = 0.01, tvcl_hep_F = 0.01, tvcl_ren_M = 0.01, tvcl_ren_F = 0.01, tvq = 0.01, tvvp = 10, Ω = [0.01 0.0; 0.0 0.01], σ_add = 0.1, σ_prop = 0.1, dCRCL_M = 0.9, dCRCL_F = 0.9)

fit_covariate_model_wt_crcl_sex =
    fit(covariate_model_wt_crcl_sex, pop, iparams_covariate_model_wt_crcl_sex, FOCE())

[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter     Function value   Gradient norm 
     0     8.554152e+03     5.650279e+03
 * time: 2.09808349609375e-5
     1     3.641387e+03     1.432080e+03
 * time: 1.5839121341705322
     2     3.290450e+03     5.274782e+02
 * time: 1.6208200454711914
     3     3.185512e+03     2.173676e+02
 * time: 1.6551411151885986
     4     3.143009e+03     1.479653e+02
 * time: 1.6878199577331543
     5     3.128511e+03     8.980031e+01
 * time: 1.7197270393371582
     6     3.123188e+03     5.033125e+01
 * time: 1.7514209747314453
     7     3.120794e+03     4.279722e+01
 * time: 1.8529109954833984
     8     3.118627e+03     3.971051e+01
 * time: 1.8856861591339111
     9     3.115300e+03     8.456587e+01
 * time: 1.9188530445098877
    10     3.109353e+03     1.350354e+02
 * time: 1.951909065246582
    11     3.095894e+03     1.998258e+02
 * time: 1.9859609603881836
    12     2.988214e+03     4.366433e+02
 * time: 2.0309650897979736
    13     2.896081e+03     5.505943e+02
 * time: 2.1456849575042725
    14     2.652467e+03     7.300323e+02
 * time: 2.781512975692749
    15     2.560937e+03     6.973661e+02
 * time: 2.9837210178375244
    16     2.254941e+03     2.740033e+02
 * time: 3.022353172302246
    17     2.222509e+03     2.034303e+02
 * time: 3.070847988128662
    18     2.171255e+03     2.449580e+02
 * time: 3.1061251163482666
    19     2.024532e+03     1.121511e+02
 * time: 3.152337074279785
    20     1.993723e+03     1.042814e+02
 * time: 3.1856751441955566
    21     1.985113e+03     8.079014e+01
 * time: 3.218414068222046
    22     1.976757e+03     7.054196e+01
 * time: 3.263604164123535
    23     1.969970e+03     6.070322e+01
 * time: 3.2961840629577637
    24     1.961095e+03     6.810782e+01
 * time: 3.3399341106414795
    25     1.947983e+03     8.116920e+01
 * time: 3.3723151683807373
    26     1.930371e+03     8.530051e+01
 * time: 3.404478073120117
    27     1.910209e+03     6.993170e+01
 * time: 3.4489941596984863
    28     1.899107e+03     3.362640e+01
 * time: 3.481872081756592
    29     1.898022e+03     2.642220e+01
 * time: 3.512892007827759
    30     1.897055e+03     1.213144e+01
 * time: 3.556107997894287
    31     1.896596e+03     7.773239e+00
 * time: 3.5874831676483154
    32     1.896538e+03     7.997039e+00
 * time: 3.629601001739502
    33     1.896451e+03     8.160909e+00
 * time: 3.6602940559387207
    34     1.896283e+03     8.237721e+00
 * time: 3.690864086151123
    35     1.895903e+03     1.520219e+01
 * time: 3.733293056488037
    36     1.895272e+03     2.358916e+01
 * time: 3.7649471759796143
    37     1.894536e+03     2.461296e+01
 * time: 3.7971110343933105
    38     1.893995e+03     1.546128e+01
 * time: 3.840294122695923
    39     1.893858e+03     6.976137e+00
 * time: 3.871062994003296
    40     1.893833e+03     6.019466e+00
 * time: 3.913194179534912
    41     1.893786e+03     3.827201e+00
 * time: 3.944094181060791
    42     1.893714e+03     3.323412e+00
 * time: 3.975314140319824
    43     1.893592e+03     3.215150e+00
 * time: 4.017865180969238
    44     1.893435e+03     6.534965e+00
 * time: 4.049548149108887
    45     1.893286e+03     7.424154e+00
 * time: 4.08037805557251
    46     1.893190e+03     5.552627e+00
 * time: 4.122640132904053
    47     1.893139e+03     3.222316e+00
 * time: 4.153318166732788
    48     1.893120e+03     3.015339e+00
 * time: 4.194887161254883
    49     1.893107e+03     3.244809e+00
 * time: 4.224240064620972
    50     1.893080e+03     6.163100e+00
 * time: 4.253365993499756
    51     1.893027e+03     9.824713e+00
 * time: 4.296349048614502
    52     1.892912e+03     1.390100e+01
 * time: 4.327911138534546
    53     1.892734e+03     1.510937e+01
 * time: 4.359250068664551
    54     1.892561e+03     1.008563e+01
 * time: 4.403414011001587
    55     1.892485e+03     3.730668e+00
 * time: 4.43428111076355
    56     1.892471e+03     3.380261e+00
 * time: 4.4768900871276855
    57     1.892463e+03     3.167904e+00
 * time: 4.507188081741333
    58     1.892441e+03     4.152065e+00
 * time: 4.537341117858887
    59     1.892391e+03     7.355996e+00
 * time: 4.580203056335449
    60     1.892268e+03     1.195397e+01
 * time: 4.611562967300415
    61     1.892026e+03     1.640783e+01
 * time: 4.642657041549683
    62     1.891735e+03     1.593576e+01
 * time: 4.68657112121582
    63     1.891569e+03     8.316423e+00
 * time: 4.718667984008789
    64     1.891494e+03     3.948212e+00
 * time: 4.762048959732056
    65     1.891481e+03     3.911593e+00
 * time: 4.79320502281189
    66     1.891457e+03     3.875559e+00
 * time: 4.824079990386963
    67     1.891405e+03     3.811247e+00
 * time: 4.86665415763855
    68     1.891262e+03     3.657045e+00
 * time: 4.897396087646484
    69     1.890930e+03     4.957405e+00
 * time: 4.928711175918579
    70     1.890317e+03     6.657726e+00
 * time: 4.973222017288208
    71     1.889660e+03     6.086302e+00
 * time: 5.005815029144287
    72     1.889303e+03     2.270929e+00
 * time: 5.050161123275757
    73     1.889253e+03     7.695301e-01
 * time: 5.08139705657959
    74     1.889252e+03     7.382144e-01
 * time: 5.113001108169556
    75     1.889251e+03     7.187898e-01
 * time: 5.155567169189453
    76     1.889251e+03     7.215047e-01
 * time: 5.185218095779419
    77     1.889250e+03     7.235155e-01
 * time: 5.214817047119141
    78     1.889249e+03     7.246818e-01
 * time: 5.256647109985352
    79     1.889244e+03     7.257796e-01
 * time: 5.290833950042725
    80     1.889233e+03     7.198190e-01
 * time: 5.321987152099609
    81     1.889204e+03     1.089029e+00
 * time: 5.367544174194336
    82     1.889142e+03     1.801601e+00
 * time: 5.3989479541778564
    83     1.889043e+03     2.967917e+00
 * time: 5.4421820640563965
    84     1.888889e+03     2.965856e+00
 * time: 5.47314715385437
    85     1.888705e+03     5.933555e-01
 * time: 5.504568099975586
    86     1.888655e+03     9.577698e-01
 * time: 5.546908140182495
    87     1.888582e+03     1.498494e+00
 * time: 5.577831983566284
    88     1.888533e+03     1.502751e+00
 * time: 5.608208179473877
    89     1.888490e+03     1.184664e+00
 * time: 5.650274991989136
    90     1.888480e+03     6.684515e-01
 * time: 5.681398153305054
    91     1.888476e+03     3.680032e-01
 * time: 5.723778963088989
    92     1.888476e+03     4.720040e-01
 * time: 5.754248142242432
    93     1.888476e+03     4.768646e-01
 * time: 5.784105062484741
    94     1.888475e+03     4.736674e-01
 * time: 5.825283050537109
    95     1.888475e+03     4.552765e-01
 * time: 5.855190992355347
    96     1.888474e+03     5.193725e-01
 * time: 5.885048151016235
    97     1.888473e+03     8.850098e-01
 * time: 5.926819086074829
    98     1.888468e+03     1.461598e+00
 * time: 5.956557035446167
    99     1.888458e+03     2.209124e+00
 * time: 5.997897148132324
   100     1.888437e+03     2.961236e+00
 * time: 6.028614044189453
   101     1.888407e+03     2.978462e+00
 * time: 6.059373140335083
   102     1.888384e+03     1.707199e+00
 * time: 6.10138201713562
   103     1.888381e+03     6.198983e-01
 * time: 6.132280111312866
   104     1.888380e+03     5.171082e-01
 * time: 6.162866115570068
   105     1.888378e+03     1.037321e-01
 * time: 6.203780174255371
   106     1.888378e+03     8.473253e-02
 * time: 6.233283042907715
   107     1.888378e+03     8.364965e-02
 * time: 6.262382984161377
   108     1.888378e+03     8.080442e-02
 * time: 6.303958177566528
   109     1.888378e+03     7.873900e-02
 * time: 6.33322811126709
   110     1.888378e+03     7.798398e-02
 * time: 6.373461008071899
   111     1.888378e+03     7.788170e-02
 * time: 6.401585102081299
   112     1.888378e+03     7.776461e-02
 * time: 6.430546045303345
   113     1.888378e+03     9.023586e-02
 * time: 6.471481084823608
   114     1.888378e+03     1.631370e-01
 * time: 6.500964164733887
   115     1.888378e+03     2.768691e-01
 * time: 6.529466152191162
   116     1.888377e+03     4.462313e-01
 * time: 6.5701329708099365
   117     1.888377e+03     6.643167e-01
 * time: 6.599853992462158
   118     1.888375e+03     8.433175e-01
 * time: 6.629081964492798
   119     1.888374e+03     7.596472e-01
 * time: 6.671010971069336
   120     1.888373e+03     3.637851e-01
 * time: 6.701253175735474
   121     1.888372e+03     8.305222e-02
 * time: 6.741498947143555
   122     1.888372e+03     2.084516e-02
 * time: 6.770802021026611
   123     1.888372e+03     2.056416e-02
 * time: 6.799774169921875
   124     1.888372e+03     2.044079e-02
 * time: 6.839878082275391
   125     1.888372e+03     2.035197e-02
 * time: 6.867856025695801
   126     1.888372e+03     2.021266e-02
 * time: 6.89572811126709
   127     1.888372e+03     1.998168e-02
 * time: 6.935342073440552
   128     1.888372e+03     3.162095e-02
 * time: 6.964091062545776
   129     1.888372e+03     5.510010e-02
 * time: 6.992335081100464
   130     1.888372e+03     9.277181e-02
 * time: 7.032448053359985
   131     1.888372e+03     1.528968e-01
 * time: 7.061820030212402
   132     1.888372e+03     2.462113e-01
 * time: 7.1024861335754395
   133     1.888372e+03     3.799882e-01
 * time: 7.132159948348999
   134     1.888371e+03     5.312359e-01
 * time: 7.1618499755859375
   135     1.888369e+03     6.019768e-01
 * time: 7.204041957855225
   136     1.888367e+03     4.665349e-01
 * time: 7.233881950378418
   137     1.888366e+03     1.404917e-01
 * time: 7.26344895362854
   138     1.888365e+03     8.513280e-02
 * time: 7.304628133773804
   139     1.888364e+03     1.244286e-01
 * time: 7.334166049957275
   140     1.888364e+03     1.028231e-01
 * time: 7.363621950149536
   141     1.888364e+03     5.163329e-02
 * time: 7.405181169509888
   142     1.888364e+03     5.148613e-02
 * time: 7.434988021850586
   143     1.888364e+03     3.147247e-02
 * time: 7.475337028503418
   144     1.888364e+03     2.104597e-02
 * time: 7.504127025604248
   145     1.888364e+03     6.541603e-03
 * time: 7.532711029052734
   146     1.888364e+03     2.538497e-03
 * time: 7.5729920864105225
   147     1.888364e+03     4.361855e-03
 * time: 7.601240158081055
   148     1.888364e+03     3.034846e-03
 * time: 7.629408121109009
   149     1.888364e+03     5.961481e-04
 * time: 7.6694371700286865

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

metrics_table(fit_covariate_model_wt)

17×2 DataFrame

Row	Metric	Value
	String	Any
1	Successful	true
2	Estimation Time	2.823
3	Subjects	30
4	Fixed Parameters	0
5	Optimized Parameters	8
6	DV Active Observations	540
7	DV Missing Observations	0
8	Total Active Observations	540
9	Total Missing Observations	0
10	Likelihood Approximation	Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}
11	LogLikelihood (LL)	-1896.46
12	-2LL	3792.93
13	AIC	3808.93
14	BIC	3843.26
15	(η-shrinkage) η₁	-0.014
16	(η-shrinkage) η₂	-0.012
17	(ϵ-shrinkage) DV	0.056

metrics_table(fit_covariate_model_wt_crcl)

17×2 DataFrame

Row	Metric	Value
	String	Any
1	Successful	true
2	Estimation Time	4.121
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.67
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	5.417	2.823	4.121	7.67
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.6299849518064

aic(fit_covariate_model_wt)

3808.9264607805962

aic(fit_covariate_model_wt_crcl)

3804.817947371705

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: 1.811981201171875e-5
     1     2.994747e+03     1.083462e+03
 * time: 2.09023118019104
     2     2.406265e+03     1.884408e+02
 * time: 2.097912073135376
     3     2.344175e+03     7.741610e+01
 * time: 2.1054000854492188
     4     2.323153e+03     2.907642e+01
 * time: 2.1128220558166504
     5     2.318222e+03     2.273295e+01
 * time: 2.12027907371521
     6     2.316833e+03     1.390527e+01
 * time: 2.1276471614837646
     7     2.316425e+03     4.490883e+00
 * time: 2.135193109512329
     8     2.316362e+03     9.374519e-01
 * time: 2.14237117767334
     9     2.316356e+03     1.928785e-01
 * time: 2.149466037750244
    10     2.316355e+03     3.119615e-02
 * time: 2.1564791202545166
    11     2.316355e+03     6.215513e-03
 * time: 2.163626194000244
    12     2.316355e+03     8.313206e-04
 * time: 2.170772075653076

FittedPumasModel

Dynamical system type:          No dynamical model

Number of subjects:                            160

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

Number of parameters:      Constant      Optimized
                                  0              4

Likelihood approximation:              NaivePooled
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -2316.3554

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

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

7 Missing Data in Covariates

Pumas will impute missing covariate values by constant interpolation of the available covariate values. If, for any subject, all of the covariate’s values are missing, Pumas will throw an error while parsing the data with read_pumas.

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

Hence, for NOCB, you can use the following:

pop = read_pumas(pkdata; covariates_direction = :right)

along with any other required keyword arguments for column mapping.

Note

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

8 Categorical Covariates

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

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

For example:

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

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

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

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

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

9 Conclusion

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

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

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

Reuse

CC BY-SA 4.0

1 nlme_sample Dataset

2 Step 1 - Parse Data into a Population

3 Step 2 - Base Model

4 Step 3 - Covariate Model

5 Step 4 - Model Comparison

5.1 Goodness of Fit Plots

6 Time-Varying Covariates

6.1 painord Dataset

7 Missing Data in Covariates

8 Categorical Covariates

9 Conclusion

Reuse

1 `nlme_sample` Dataset

2 Step 1 - Parse Data into a `Population`

6.1 `painord` Dataset