using Dates
using Pumas
using PumasUtilities
using CairoMakie
using DataFramesMeta
using CSV
using PharmaDatasets
Covariate Models
Caution
Some functions in this tutorial are only available after you load the PumasUtilities package.
1 Covariate Model Building
In this tutorial we’ll cover a workflow around covariate model building. You’ll learn how to:
- 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.1 nlme_sample Dataset
For this tutorial we’ll use the nlme_sample dataset from PharmaDatasets.jl:
pkfile = dataset("nlme_sample", String)
pkdata = CSV.read(pkfile, DataFrame; missingstring = ["NA", ""])
first(pkdata, 5)5×15 DataFrame
| Row | ID | TIME | DV | AMT | EVID | CMT | RATE | WT | AGE | SEX | CRCL | GROUP | ROUTE | DURATION | OCC |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Int64 | Float64 | Float64? | Int64? | Int64 | Int64? | Int64 | Int64 | Int64 | String1 | Int64 | String7 | Float64 | Int64? | Int64 | |
| 1 | 1 | 0.0 | missing | 1000 | 1 | 1 | 500 | 90 | 47 | M | 75 | 1000 mg | Inf | 2 | 1 |
| 2 | 1 | 0.001 | 0.0667231 | missing | 0 | missing | 0 | 90 | 47 | M | 75 | 1000 mg | Inf | missing | 1 |
| 3 | 1 | 1.0 | 112.817 | missing | 0 | missing | 0 | 90 | 47 | M | 75 | 1000 mg | Inf | missing | 1 |
| 4 | 1 | 2.0 | 224.087 | missing | 0 | missing | 0 | 90 | 47 | M | 75 | 1000 mg | Inf | missing | 1 |
| 5 | 1 | 4.0 | 220.047 | missing | 0 | missing | 0 | 90 | 47 | M | 75 | 1000 mg | Inf | missing | 1 |
Note
The nlme_sample dataset has different missing values as the standard datasets in the PharmaDatasets.jl. That’s why we are first getting a String representation of the dataset as a CSV file as pkfile variable. Then, we use it to customize the missingstring keyword argument inside CSV.read function in order to have a working DataFrame for the nlme_sample dataset.
If you want to know more about data wrangling and how to read and write data in different formats, please check out the Data Wrangling Tutorials at tutorials.pumas.ai.
As you can see, the nlme_sample dataset has the standard PK dataset columns such as :ID, :TIME, :DV, :AMT, :EVID and :CMT. The dataset also contains the following list of covariates:
:WT: subject weight in kilograms:SEX: subject sex, either"F"or"M":CRCL: subject creatinine clearance:GROUP: subject dosing group, either"500 mg","750 mg", or"1000 mg"
And we’ll learn how to include them in our Pumas modeling workflows.
describe(pkdata, :mean, :std, :nunique, :first, :last, :eltype)15×7 DataFrame
| Row | variable | mean | std | nunique | first | last | eltype |
|---|---|---|---|---|---|---|---|
| Symbol | Union… | Union… | Union… | Any | Any | Type | |
| 1 | ID | 15.5 | 8.661 | 1 | 30 | Int64 | |
| 2 | TIME | 82.6527 | 63.2187 | 0.0 | 168.0 | Float64 | |
| 3 | DV | 157.315 | 110.393 | missing | missing | Union{Missing, Float64} | |
| 4 | AMT | 750.0 | 204.551 | 1000 | 500 | Union{Missing, Int64} | |
| 5 | EVID | 0.307692 | 0.461835 | 1 | 1 | Int64 | |
| 6 | CMT | 1.0 | 0.0 | 1 | 1 | Union{Missing, Int64} | |
| 7 | RATE | 115.385 | 182.218 | 500 | 250 | Int64 | |
| 8 | WT | 81.6 | 11.6051 | 90 | 96 | Int64 | |
| 9 | AGE | 40.0333 | 11.6479 | 47 | 56 | Int64 | |
| 10 | SEX | 2 | M | F | String1 | ||
| 11 | CRCL | 72.5667 | 26.6212 | 75 | 90 | Int64 | |
| 12 | GROUP | 3 | 1000 mg | 500 mg | String7 | ||
| 13 | ROUTE | Inf | NaN | Inf | Inf | Float64 | |
| 14 | DURATION | 2.0 | 0.0 | 2 | 2 | Union{Missing, Int64} | |
| 15 | OCC | 4.15385 | 2.62836 | 1 | 8 | Int64 |
Tip
As you can see, all these covariates are constant. That means, they do not vary over time. We’ll also cover time-varying covariates later in this tutorial.
1.2 Step 1 - Parse Data into a Population
The first step in our covariate model building workflow is to parse data into a Population. This is accomplished with the read_pumas function. Here we are to use the covariates keyword argument to pass a vector of column names to be parsed as covariates:
pop = read_pumas(
pkdata;
id = :ID,
time = :TIME,
amt = :AMT,
covariates = [:WT, :AGE, :SEX, :CRCL, :GROUP],
observations = [:DV],
cmt = :CMT,
evid = :EVID,
rate = :RATE,
)Population
Subjects: 30
Covariates: WT, AGE, SEX, CRCL, GROUP
Observations: DVDue to Pumas’ dynamic workflow capabilities, we only need to define our Population once. That is, we tell read_pumas to parse all the covariates available, even if we will be using none or a subset of those in our models.
This is a feature that highly increases workflow efficiency in developing and fitting models.
1.3 Step 2 - Base Model
The second step of our covariate model building workflow is to develop a base model, i.e., a model without any covariate effects on its parameters. This represents the null model against which covariate models can be tested after checking if covariate inclusion is helpful in our model.
Let’s create a combined-error simple 2-compartment base model:
base_model = @model begin
@param begin
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvq ∈ RealDomain(; lower = 0)
tvvp ∈ RealDomain(; lower = 0)
Ω ∈ PDiagDomain(2)
σ_add ∈ RealDomain(; lower = 0)
σ_prop ∈ RealDomain(; lower = 0)
end
@random begin
η ~ MvNormal(Ω)
end
@pre begin
CL = tvcl * exp(η[1])
Vc = tvvc * exp(η[2])
Q = tvq
Vp = tvvp
end
@dynamics Central1Periph1
@derived begin
cp := @. Central / Vc
DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
end
endPumasModel
Parameters: tvcl, tvvc, tvq, tvvp, Ω, σ_add, σ_prop
Random effects: η
Covariates:
Dynamical system variables: Central, Peripheral
Dynamical system type: Closed form
Derived: DV
Observed: DVAnd 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.023867130279541016
1 3.110315e+03 9.706222e+02
* time: 0.9255001544952393
2 2.831659e+03 7.817006e+02
* time: 1.0499961376190186
3 2.405281e+03 2.923696e+02
* time: 1.1079251766204834
4 2.370406e+03 3.032286e+02
* time: 1.1468040943145752
5 2.313631e+03 3.126188e+02
* time: 1.1874279975891113
6 2.263986e+03 2.946697e+02
* time: 1.2246451377868652
7 2.160182e+03 1.917599e+02
* time: 1.2819831371307373
8 2.112467e+03 1.588288e+02
* time: 1.3332149982452393
9 2.090339e+03 1.108334e+02
* time: 1.3695480823516846
10 2.078171e+03 8.108027e+01
* time: 1.4461190700531006
11 2.074517e+03 7.813928e+01
* time: 1.4788429737091064
12 2.066270e+03 7.044632e+01
* time: 1.5106379985809326
13 2.049660e+03 1.062584e+02
* time: 1.54243803024292
14 2.021965e+03 1.130570e+02
* time: 1.5754940509796143
15 1.994936e+03 7.825801e+01
* time: 1.609158992767334
16 1.979337e+03 5.318263e+01
* time: 1.6435000896453857
17 1.972141e+03 6.807046e+01
* time: 1.6773631572723389
18 1.967973e+03 7.896361e+01
* time: 1.710414171218872
19 1.962237e+03 8.343757e+01
* time: 1.7433831691741943
20 1.952791e+03 5.565304e+01
* time: 1.7826240062713623
21 1.935857e+03 3.923284e+01
* time: 1.8359081745147705
22 1.926254e+03 5.749643e+01
* time: 1.8722960948944092
23 1.922144e+03 4.306225e+01
* time: 1.9055631160736084
24 1.911566e+03 4.810496e+01
* time: 1.9419910907745361
25 1.906464e+03 4.324267e+01
* time: 1.9776661396026611
26 1.905339e+03 1.207954e+01
* time: 2.0093209743499756
27 1.905092e+03 1.093479e+01
* time: 2.0404410362243652
28 1.904957e+03 1.057034e+01
* time: 2.0720791816711426
29 1.904875e+03 1.060882e+01
* time: 2.1024529933929443
30 1.904459e+03 1.031525e+01
* time: 2.1344900131225586
31 1.903886e+03 1.041179e+01
* time: 2.190859079360962
32 1.903313e+03 1.135672e+01
* time: 2.225054979324341
33 1.903057e+03 1.075683e+01
* time: 2.256148099899292
34 1.902950e+03 1.091295e+01
* time: 2.28729510307312
35 1.902887e+03 1.042409e+01
* time: 2.3177881240844727
36 1.902640e+03 9.213043e+00
* time: 2.349245071411133
37 1.902364e+03 9.519321e+00
* time: 2.381190061569214
38 1.902156e+03 5.590984e+00
* time: 2.4127399921417236
39 1.902094e+03 1.811898e+00
* time: 2.443247079849243
40 1.902086e+03 1.644722e+00
* time: 2.4735281467437744
41 1.902084e+03 1.653520e+00
* time: 2.5027551651000977
42 1.902074e+03 1.720184e+00
* time: 2.535752058029175
43 1.902055e+03 2.619061e+00
* time: 2.580761194229126
44 1.902015e+03 3.519729e+00
* time: 2.6115171909332275
45 1.901962e+03 3.403372e+00
* time: 2.6417980194091797
46 1.901924e+03 1.945644e+00
* time: 2.672307014465332
47 1.901914e+03 6.273342e-01
* time: 2.7024340629577637
48 1.901913e+03 5.374557e-01
* time: 2.731743097305298
49 1.901913e+03 5.710007e-01
* time: 2.760017156600952
50 1.901913e+03 6.091390e-01
* time: 2.788767099380493
51 1.901912e+03 6.692417e-01
* time: 2.818356990814209
52 1.901909e+03 7.579620e-01
* time: 2.8478009700775146
53 1.901903e+03 8.798211e-01
* time: 2.87833309173584
54 1.901889e+03 1.002981e+00
* time: 2.930122137069702
55 1.901864e+03 1.495512e+00
* time: 2.960520029067993
56 1.901837e+03 1.664621e+00
* time: 2.989129066467285
57 1.901819e+03 8.601118e-01
* time: 3.018311023712158
58 1.901815e+03 4.525470e-02
* time: 3.04726505279541
59 1.901815e+03 1.294280e-02
* time: 3.074526071548462
60 1.901815e+03 2.876568e-03
* time: 3.1003291606903076
61 1.901815e+03 2.876568e-03
* time: 3.1742069721221924
62 1.901815e+03 2.876568e-03
* time: 3.251160144805908FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -1901.815
Number of subjects: 30
Number of parameters: Fixed Optimized
0 8
Observation records: Active Missing
DV: 540 0
Total: 540 0
---------------------
Estimate
---------------------
tvcl 0.1542
tvvc 4.5856
tvq 0.042341
tvvp 3.7082
Ω₁,₁ 0.26467
Ω₂,₂ 0.10627
σ_add 0.032183
σ_prop 0.061005
---------------------1.4 Step 3 - Covariate Model
The third step of our covariate model building workflow is to actually develop one or more covariate models. Let’s develop three covariate models:
- 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 ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
end
endPumasModel
Parameters: tvcl, tvvc, tvq, tvvp, Ω, σ_add, σ_prop
Random effects: η
Covariates: WT
Dynamical system variables: Central, Peripheral
Dynamical system type: Closed form
Derived: DV
Observed: DVOnce 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: 7.295608520507812e-5
1 2.846050e+03 1.128657e+03
* time: 0.07731103897094727
2 2.472982e+03 7.008264e+02
* time: 0.12812185287475586
3 2.180690e+03 2.312704e+02
* time: 0.17886686325073242
4 2.125801e+03 1.862929e+02
* time: 0.22588682174682617
5 2.103173e+03 1.320946e+02
* time: 0.2716488838195801
6 2.091136e+03 1.103035e+02
* time: 0.3157188892364502
7 2.081443e+03 1.091133e+02
* time: 0.3547179698944092
8 2.071793e+03 7.197675e+01
* time: 0.39705491065979004
9 2.062706e+03 7.623310e+01
* time: 0.48050594329833984
10 2.057515e+03 6.885476e+01
* time: 0.5197558403015137
11 2.051133e+03 6.368504e+01
* time: 0.5571529865264893
12 2.038626e+03 7.730243e+01
* time: 0.5969710350036621
13 2.019352e+03 1.136864e+02
* time: 0.6387348175048828
14 1.997136e+03 1.005748e+02
* time: 0.680793046951294
15 1.983023e+03 6.831478e+01
* time: 0.7253198623657227
16 1.977700e+03 5.649783e+01
* time: 0.7699840068817139
17 1.974583e+03 6.357642e+01
* time: 0.8132219314575195
18 1.967292e+03 7.658974e+01
* time: 0.8613448143005371
19 1.951045e+03 6.130573e+01
* time: 0.9947528839111328
20 1.935868e+03 4.845839e+01
* time: 1.073441982269287
21 1.929356e+03 6.325782e+01
* time: 1.115968942642212
22 1.925187e+03 3.142245e+01
* time: 1.154390811920166
23 1.923733e+03 4.623400e+01
* time: 1.1925559043884277
24 1.918498e+03 5.347738e+01
* time: 1.231017827987671
25 1.912383e+03 5.849125e+01
* time: 1.2760870456695557
26 1.905510e+03 3.254038e+01
* time: 1.319382905960083
27 1.903629e+03 2.905618e+01
* time: 1.356947898864746
28 1.902833e+03 2.907696e+01
* time: 1.3948018550872803
29 1.902447e+03 2.746037e+01
* time: 1.4332950115203857
30 1.899399e+03 1.930949e+01
* time: 1.498849868774414
31 1.898705e+03 1.186361e+01
* time: 1.5373878479003906
32 1.898505e+03 1.050402e+01
* time: 1.5747230052947998
33 1.898474e+03 1.042186e+01
* time: 1.6071648597717285
34 1.897992e+03 1.238729e+01
* time: 1.6406848430633545
35 1.897498e+03 1.729368e+01
* time: 1.6750450134277344
36 1.896954e+03 1.472554e+01
* time: 1.709561824798584
37 1.896744e+03 5.852825e+00
* time: 1.743783950805664
38 1.896705e+03 1.171353e+00
* time: 1.7751789093017578
39 1.896704e+03 1.216117e+00
* time: 1.807633876800537
40 1.896703e+03 1.230336e+00
* time: 1.8401670455932617
41 1.896698e+03 1.250715e+00
* time: 1.8932268619537354
42 1.896688e+03 1.449552e+00
* time: 1.9270319938659668
43 1.896666e+03 2.533300e+00
* time: 1.9605610370635986
44 1.896631e+03 3.075536e+00
* time: 1.995013952255249
45 1.896599e+03 2.139797e+00
* time: 2.0295979976654053
46 1.896587e+03 6.636031e-01
* time: 2.0644519329071045
47 1.896585e+03 6.303517e-01
* time: 2.098475933074951
48 1.896585e+03 5.995265e-01
* time: 2.1298019886016846
49 1.896584e+03 5.844159e-01
* time: 2.161263942718506
50 1.896583e+03 6.083858e-01
* time: 2.19382381439209
51 1.896579e+03 8.145326e-01
* time: 2.2262020111083984
52 1.896570e+03 1.293490e+00
* time: 2.2592320442199707
53 1.896549e+03 1.877705e+00
* time: 2.3152809143066406
54 1.896513e+03 2.217391e+00
* time: 2.348289966583252
55 1.896482e+03 1.658147e+00
* time: 2.38092303276062
56 1.896466e+03 5.207215e-01
* time: 2.4134879112243652
57 1.896463e+03 1.177468e-01
* time: 2.4453699588775635
58 1.896463e+03 1.678937e-02
* time: 2.475184917449951
59 1.896463e+03 2.666635e-03
* time: 2.503002882003784
60 1.896463e+03 2.666635e-03
* time: 2.5638539791107178FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -1896.4632
Number of subjects: 30
Number of parameters: Fixed Optimized
0 8
Observation records: Active Missing
DV: 540 0
Total: 540 0
---------------------
Estimate
---------------------
tvcl 0.13915
tvvc 3.9754
tvq 0.041988
tvvp 3.5722
Ω₁,₁ 0.23874
Ω₂,₂ 0.081371
σ_add 0.032174
σ_prop 0.061012
---------------------We can definitely see that, despite not increasing the parameters of the model, our loglikelihood is higher (implying lower objective function value). Furthermore, our additive and combined error values remain almost the same while our Ωs decreased for CL and Vc. This implies that the WT covariate is definitely assisting in a better model fit by capturing the variability in CL and Vc. We’ll compare models in the next step.
Let’s now try to incorporate into this model creatinine clearance (CRCL) effect on clearance (CL):
Tip
Like the tip above, you can derive covariates or perform any required transformations or scaling as a data wrangling step and pass the derived/transformed into read_pumas as a covariate. In this particular case, it is easy to create three more columns in the original data as:
@rtransform! pkdata :CRCL_CL = :CRCL / 100
@rtransform! pkdata :ASWT_CL = (:WT / 70)^0.75
@rtransform! pkdata :ASWT_Vc = (:WT / 70)Then, you bring these covariates into read_pumas and use them directly on renCL, CL and Vc. This will avoid computations of your transformations during the model fit.
covariate_model_wt_crcl = @model begin
@param begin
tvvc ∈ RealDomain(; lower = 0)
tvq ∈ RealDomain(; lower = 0)
tvvp ∈ RealDomain(; lower = 0)
tvcl_hep ∈ RealDomain(; lower = 0)
tvcl_ren ∈ RealDomain(; lower = 0)
Ω ∈ PDiagDomain(2)
σ_add ∈ RealDomain(; lower = 0)
σ_prop ∈ RealDomain(; lower = 0)
dCRCL ∈ RealDomain()
end
@random begin
η ~ MvNormal(Ω)
end
@covariates begin
WT
CRCL
end
@pre begin
hepCL = tvcl_hep * (WT / 70)^0.75
renCL = tvcl_ren * (CRCL / 100)^dCRCL
CL = (hepCL + renCL) * exp(η[1])
Vc = tvvc * (WT / 70) * exp(η[2])
Q = tvq
Vp = tvvp
end
@dynamics Central1Periph1
@derived begin
cp := @. Central / Vc
DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
end
endPumasModel
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: DVIn the covariate_model_wt_crcl model we are keeping our allometric scaling on WT from before. But we are also adding a new covariate creatinine clearance (CRCL), dividing clearance (CL) into hepatic (hepCL) and renal clearance (renCL), along with a new parameter dCRCL.
dCRCL is the exponent of the power function for the effect of creatinine clearance on renal clearance. In some models this parameter is fixed, however we’ll allow the model to estimate it.
This is a good example on how to add covariate coefficients such as dCRCL in any Pumas covariate model. Now, let’s fit this model. Note that we need a new initial parameters values’ list since the previous one we used doesn’t include dCRCL, tvcl_hep or tvcl_ren:
iparams_covariate_model_wt_crcl = (;
tvvc = 5,
tvcl_hep = 0.01,
tvcl_ren = 0.01,
tvq = 0.01,
tvvp = 10,
Ω = Diagonal([0.01, 0.01]),
σ_add = 0.1,
σ_prop = 0.1,
dCRCL = 0.9,
)(tvvc = 5,
tvcl_hep = 0.01,
tvcl_ren = 0.01,
tvq = 0.01,
tvvp = 10,
Ω = [0.01 0.0; 0.0 0.01],
σ_add = 0.1,
σ_prop = 0.1,
dCRCL = 0.9,)fit_covariate_model_wt_crcl =
fit(covariate_model_wt_crcl, pop, iparams_covariate_model_wt_crcl, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 8.554152e+03 5.650279e+03
* time: 7.605552673339844e-5
1 3.572050e+03 1.302046e+03
* time: 0.09332609176635742
2 3.266947e+03 5.384244e+02
* time: 0.15255498886108398
3 3.150635e+03 1.918079e+02
* time: 0.20343399047851562
4 3.108122e+03 1.277799e+02
* time: 0.25118112564086914
5 3.091358e+03 8.838080e+01
* time: 0.29942798614501953
6 3.082997e+03 8.321689e+01
* time: 0.4316599369049072
7 3.076379e+03 8.167702e+01
* time: 0.4797680377960205
8 3.067422e+03 1.177822e+02
* time: 0.5289099216461182
9 3.048580e+03 2.526969e+02
* time: 0.5879161357879639
10 2.993096e+03 6.325396e+02
* time: 0.6720559597015381
11 2.965744e+03 7.634718e+02
* time: 0.8442749977111816
12 2.921212e+03 9.704020e+02
* time: 0.951228141784668
13 2.553649e+03 6.642510e+02
* time: 1.084792137145996
14 2.319495e+03 3.204711e+02
* time: 1.2024660110473633
15 2.183040e+03 2.174905e+02
* time: 1.293044090270996
16 2.157621e+03 3.150983e+02
* time: 1.3458540439605713
17 2.132395e+03 2.847614e+02
* time: 1.3980381488800049
18 2.084799e+03 1.563370e+02
* time: 1.4487149715423584
19 2.071497e+03 1.006429e+02
* time: 1.498534917831421
20 2.064983e+03 9.753313e+01
* time: 1.5517489910125732
21 2.048289e+03 9.230405e+01
* time: 1.6274900436401367
22 2.020938e+03 1.292359e+02
* time: 1.6767370700836182
23 1.983888e+03 2.284990e+02
* time: 1.725870132446289
24 1.962132e+03 1.220188e+02
* time: 1.7720561027526855
25 1.945947e+03 1.035894e+02
* time: 1.815904140472412
26 1.917782e+03 8.246698e+01
* time: 1.8625199794769287
27 1.905967e+03 5.408054e+01
* time: 1.9142839908599854
28 1.898569e+03 2.172222e+01
* time: 1.966433048248291
29 1.897473e+03 1.689350e+01
* time: 2.0156919956207275
30 1.897019e+03 1.666689e+01
* time: 2.068406105041504
31 1.896796e+03 1.699751e+01
* time: 2.149477958679199
32 1.896559e+03 1.645900e+01
* time: 2.1920220851898193
33 1.896223e+03 1.415504e+01
* time: 2.234653949737549
34 1.895773e+03 1.630081e+01
* time: 2.27699613571167
35 1.895309e+03 1.723930e+01
* time: 2.3201351165771484
36 1.895004e+03 1.229983e+01
* time: 2.3638460636138916
37 1.894871e+03 5.385102e+00
* time: 2.408555030822754
38 1.894827e+03 3.465463e+00
* time: 2.45468807220459
39 1.894816e+03 3.387474e+00
* time: 2.4993269443511963
40 1.894807e+03 3.295388e+00
* time: 2.543868064880371
41 1.894786e+03 3.089194e+00
* time: 2.589340925216675
42 1.894737e+03 2.928080e+00
* time: 2.6603620052337646
43 1.894624e+03 3.088723e+00
* time: 2.7032151222229004
44 1.894413e+03 3.493791e+00
* time: 2.7456071376800537
45 1.894127e+03 3.142865e+00
* time: 2.7881081104278564
46 1.893933e+03 2.145253e+00
* time: 2.8310229778289795
47 1.893888e+03 2.172800e+00
* time: 2.8732099533081055
48 1.893880e+03 2.180509e+00
* time: 2.9155540466308594
49 1.893876e+03 2.134257e+00
* time: 2.960113048553467
50 1.893868e+03 2.032354e+00
* time: 3.005000114440918
51 1.893846e+03 1.760874e+00
* time: 3.050462007522583
52 1.893796e+03 1.779016e+00
* time: 3.0964860916137695
53 1.893694e+03 2.018956e+00
* time: 3.1740500926971436
54 1.893559e+03 2.366854e+00
* time: 3.2168729305267334
55 1.893474e+03 3.690142e+00
* time: 3.2593400478363037
56 1.893446e+03 3.675109e+00
* time: 3.3012919425964355
57 1.893439e+03 3.426419e+00
* time: 3.341876983642578
58 1.893429e+03 3.183164e+00
* time: 3.381722927093506
59 1.893398e+03 2.695171e+00
* time: 3.4225690364837646
60 1.893328e+03 2.753548e+00
* time: 3.4681220054626465
61 1.893169e+03 3.589748e+00
* time: 3.5147390365600586
62 1.892920e+03 3.680718e+00
* time: 3.561738967895508
63 1.892667e+03 2.568107e+00
* time: 3.608362913131714
64 1.892514e+03 1.087910e+00
* time: 3.6783320903778076
65 1.892493e+03 3.287296e-01
* time: 3.720505952835083
66 1.892492e+03 2.967465e-01
* time: 3.7618961334228516
67 1.892492e+03 3.020682e-01
* time: 3.8059351444244385
68 1.892491e+03 3.034704e-01
* time: 3.843632936477661
69 1.892491e+03 3.091846e-01
* time: 3.8818371295928955
70 1.892491e+03 3.224170e-01
* time: 3.9196879863739014
71 1.892490e+03 6.494197e-01
* time: 3.9619200229644775
72 1.892488e+03 1.115188e+00
* time: 4.00527811050415
73 1.892483e+03 1.838833e+00
* time: 4.048789024353027
74 1.892472e+03 2.765371e+00
* time: 4.093003034591675
75 1.892452e+03 3.463807e+00
* time: 4.140151023864746
76 1.892431e+03 2.805270e+00
* time: 4.215433120727539
77 1.892411e+03 5.758916e-01
* time: 4.258139133453369
78 1.892410e+03 1.434041e-01
* time: 4.2993340492248535
79 1.892409e+03 1.639246e-01
* time: 4.339478969573975
80 1.892409e+03 1.145856e-01
* time: 4.378265142440796
81 1.892409e+03 3.966861e-02
* time: 4.417387962341309
82 1.892409e+03 3.550808e-02
* time: 4.457374095916748
83 1.892409e+03 3.456241e-02
* time: 4.497880935668945
84 1.892409e+03 3.114018e-02
* time: 4.539130926132202
85 1.892409e+03 4.080806e-02
* time: 4.580549001693726
86 1.892409e+03 6.722726e-02
* time: 4.623697996139526
87 1.892409e+03 1.006791e-01
* time: 4.690177917480469
88 1.892409e+03 1.303988e-01
* time: 4.730164051055908
89 1.892409e+03 1.228919e-01
* time: 4.7690651416778564
90 1.892409e+03 6.433813e-02
* time: 4.8077380657196045
91 1.892409e+03 1.314164e-02
* time: 4.846837043762207
92 1.892409e+03 4.929931e-04
* time: 4.885097980499268FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -1892.409
Number of subjects: 30
Number of parameters: Fixed Optimized
0 10
Observation records: Active Missing
DV: 540 0
Total: 540 0
-----------------------
Estimate
-----------------------
tvvc 3.9757
tvq 0.042177
tvvp 3.6434
tvcl_hep 0.058572
tvcl_ren 0.1337
Ω₁,₁ 0.18299
Ω₂,₂ 0.081353
σ_add 0.032174
σ_prop 0.06101
dCRCL 1.1331
-----------------------As before, our loglikelihood is higher (implying lower objective function value). Furthermore, our additive and combined error values remain almost the same while our Ω on CL, Ω₁,₁, decreased. This implies that the CRCL covariate with an estimated exponent dCRCL is definitely assisting in a better model fit.
Finally let’s include a separated CL model based on sex as a third covariate model:
covariate_model_wt_crcl_sex = @model begin
@param begin
tvvc ∈ RealDomain(; lower = 0)
tvq ∈ RealDomain(; lower = 0)
tvvp ∈ RealDomain(; lower = 0)
tvcl_hep_M ∈ RealDomain(; lower = 0)
tvcl_hep_F ∈ RealDomain(; lower = 0)
tvcl_ren_M ∈ RealDomain(; lower = 0)
tvcl_ren_F ∈ RealDomain(; lower = 0)
Ω ∈ PDiagDomain(2)
σ_add ∈ RealDomain(; lower = 0)
σ_prop ∈ RealDomain(; lower = 0)
dCRCL_M ∈ RealDomain()
dCRCL_F ∈ RealDomain()
end
@random begin
η ~ MvNormal(Ω)
end
@covariates begin
WT
CRCL
SEX
end
@pre begin
tvcl_hep = ifelse(SEX == "M", tvcl_hep_M, tvcl_hep_F)
tvcl_ren = ifelse(SEX == "M", tvcl_ren_M, tvcl_ren_F)
dCRCL = ifelse(SEX == "M", dCRCL_M, dCRCL_F)
hepCL = tvcl_hep * (WT / 70)^0.75
renCL = tvcl_ren * (CRCL / 100)^dCRCL
CL = (hepCL + renCL) * exp(η[1])
Vc = tvvc * (WT / 70) * exp(η[2])
Q = tvq
Vp = tvvp
end
@dynamics Central1Periph1
@derived begin
cp := @. Central / Vc
DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
end
endPumasModel
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: DVIn 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: 6.794929504394531e-5
1 3.641387e+03 1.432080e+03
* time: 0.09427189826965332
2 3.290450e+03 5.274782e+02
* time: 0.16293883323669434
3 3.185512e+03 2.173676e+02
* time: 0.22732996940612793
4 3.143009e+03 1.479653e+02
* time: 0.36737895011901855
5 3.128511e+03 8.980031e+01
* time: 0.41657590866088867
6 3.123188e+03 5.033125e+01
* time: 0.4649159908294678
7 3.120794e+03 4.279722e+01
* time: 0.5125699043273926
8 3.118627e+03 3.971051e+01
* time: 0.5610959529876709
9 3.115300e+03 8.456587e+01
* time: 0.6102759838104248
10 3.109353e+03 1.350354e+02
* time: 0.6609559059143066
11 3.095894e+03 1.998258e+02
* time: 0.7239890098571777
12 2.988214e+03 4.366433e+02
* time: 0.823145866394043
13 2.896081e+03 5.505943e+02
* time: 1.0892858505249023
14 2.652467e+03 7.300323e+02
* time: 2.380638837814331
15 2.560937e+03 6.973661e+02
* time: 2.7241969108581543
16 2.254941e+03 2.740033e+02
* time: 2.785456895828247
17 2.222509e+03 2.034303e+02
* time: 2.838122844696045
18 2.171255e+03 2.449580e+02
* time: 2.8910439014434814
19 2.024532e+03 1.121511e+02
* time: 2.951462984085083
20 1.993723e+03 1.042814e+02
* time: 3.010756015777588
21 1.985113e+03 8.079014e+01
* time: 3.0688118934631348
22 1.976757e+03 7.054196e+01
* time: 3.1269760131835938
23 1.969970e+03 6.070322e+01
* time: 3.1833748817443848
24 1.961095e+03 6.810782e+01
* time: 3.265570878982544
25 1.947983e+03 8.116920e+01
* time: 3.3122010231018066
26 1.930371e+03 8.530051e+01
* time: 3.358163833618164
27 1.910209e+03 6.993170e+01
* time: 3.4044368267059326
28 1.899107e+03 3.362640e+01
* time: 3.450856924057007
29 1.898022e+03 2.642220e+01
* time: 3.494062900543213
30 1.897055e+03 1.213144e+01
* time: 3.5376880168914795
31 1.896596e+03 7.773239e+00
* time: 3.584251880645752
32 1.896538e+03 7.997039e+00
* time: 3.6305599212646484
33 1.896451e+03 8.160909e+00
* time: 3.6791770458221436
34 1.896283e+03 8.237721e+00
* time: 3.758884906768799
35 1.895903e+03 1.520219e+01
* time: 3.8037729263305664
36 1.895272e+03 2.358916e+01
* time: 3.849540948867798
37 1.894536e+03 2.461296e+01
* time: 3.8944108486175537
38 1.893995e+03 1.546128e+01
* time: 3.938717842102051
39 1.893858e+03 6.976137e+00
* time: 3.981987953186035
40 1.893833e+03 6.019466e+00
* time: 4.024235010147095
41 1.893786e+03 3.827201e+00
* time: 4.067878007888794
42 1.893714e+03 3.323412e+00
* time: 4.114794969558716
43 1.893592e+03 3.215150e+00
* time: 4.163991928100586
44 1.893435e+03 6.534965e+00
* time: 4.212926864624023
45 1.893286e+03 7.424154e+00
* time: 4.290812969207764
46 1.893190e+03 5.552627e+00
* time: 4.335752964019775
47 1.893139e+03 3.222316e+00
* time: 4.379441022872925
48 1.893120e+03 3.015339e+00
* time: 4.421304941177368
49 1.893107e+03 3.244809e+00
* time: 4.462121963500977
50 1.893080e+03 6.163100e+00
* time: 4.5031418800354
51 1.893027e+03 9.824713e+00
* time: 4.544464826583862
52 1.892912e+03 1.390100e+01
* time: 4.586113929748535
53 1.892734e+03 1.510937e+01
* time: 4.63167405128479
54 1.892561e+03 1.008563e+01
* time: 4.678915977478027
55 1.892485e+03 3.730668e+00
* time: 4.727118015289307
56 1.892471e+03 3.380261e+00
* time: 4.7983479499816895
57 1.892463e+03 3.167904e+00
* time: 4.84109091758728
58 1.892441e+03 4.152065e+00
* time: 4.882704973220825
59 1.892391e+03 7.355996e+00
* time: 4.923292875289917
60 1.892268e+03 1.195397e+01
* time: 4.96540904045105
61 1.892026e+03 1.640783e+01
* time: 5.008277893066406
62 1.891735e+03 1.593576e+01
* time: 5.0517120361328125
63 1.891569e+03 8.316423e+00
* time: 5.094996929168701
64 1.891494e+03 3.948212e+00
* time: 5.1398069858551025
65 1.891481e+03 3.911593e+00
* time: 5.186386823654175
66 1.891457e+03 3.875559e+00
* time: 5.23212194442749
67 1.891405e+03 3.811247e+00
* time: 5.306467056274414
68 1.891262e+03 3.657045e+00
* time: 5.349447965621948
69 1.890930e+03 4.957405e+00
* time: 5.39349102973938
70 1.890317e+03 6.657726e+00
* time: 5.436975002288818
71 1.889660e+03 6.086302e+00
* time: 5.481097936630249
72 1.889303e+03 2.270929e+00
* time: 5.5242509841918945
73 1.889253e+03 7.695301e-01
* time: 5.566597938537598
74 1.889252e+03 7.382144e-01
* time: 5.608165979385376
75 1.889251e+03 7.187898e-01
* time: 5.652859926223755
76 1.889251e+03 7.215047e-01
* time: 5.7000508308410645
77 1.889250e+03 7.235155e-01
* time: 5.744102954864502
78 1.889249e+03 7.246818e-01
* time: 5.813565969467163
79 1.889244e+03 7.257796e-01
* time: 5.8559348583221436
80 1.889233e+03 7.198190e-01
* time: 5.898023843765259
81 1.889204e+03 1.089029e+00
* time: 5.939227819442749
82 1.889142e+03 1.801601e+00
* time: 5.981328010559082
83 1.889043e+03 2.967917e+00
* time: 6.023749828338623
84 1.888889e+03 2.965856e+00
* time: 6.065776824951172
85 1.888705e+03 5.933557e-01
* time: 6.108080863952637
86 1.888655e+03 9.577696e-01
* time: 6.152273893356323
87 1.888582e+03 1.498494e+00
* time: 6.1989428997039795
88 1.888533e+03 1.502753e+00
* time: 6.245723009109497
89 1.888490e+03 1.184664e+00
* time: 6.3207080364227295
90 1.888480e+03 6.684517e-01
* time: 6.363349914550781
91 1.888476e+03 3.680034e-01
* time: 6.405491828918457
92 1.888476e+03 4.720040e-01
* time: 6.445335865020752
93 1.888476e+03 4.768646e-01
* time: 6.484545946121216
94 1.888475e+03 4.736673e-01
* time: 6.523632049560547
95 1.888475e+03 4.552764e-01
* time: 6.5625059604644775
96 1.888474e+03 5.193733e-01
* time: 6.601455926895142
97 1.888473e+03 8.850112e-01
* time: 6.64172101020813
98 1.888468e+03 1.461600e+00
* time: 6.685330867767334
99 1.888458e+03 2.209127e+00
* time: 6.72878885269165
100 1.888437e+03 2.961239e+00
* time: 6.774442911148071
101 1.888407e+03 2.978463e+00
* time: 6.845628023147583
102 1.888384e+03 1.707201e+00
* time: 6.891321897506714
103 1.888381e+03 6.199354e-01
* time: 6.9339599609375
104 1.888380e+03 5.170909e-01
* time: 6.974484920501709
105 1.888378e+03 1.037408e-01
* time: 7.0147459506988525
106 1.888378e+03 8.473247e-02
* time: 7.052626848220825
107 1.888378e+03 8.364978e-02
* time: 7.090404033660889
108 1.888378e+03 8.080446e-02
* time: 7.128780841827393
109 1.888378e+03 7.873905e-02
* time: 7.170028924942017
110 1.888378e+03 7.798398e-02
* time: 7.212401866912842
111 1.888378e+03 7.788170e-02
* time: 7.252710819244385
112 1.888378e+03 7.776461e-02
* time: 7.32112193107605
113 1.888378e+03 9.023662e-02
* time: 7.360683917999268
114 1.888378e+03 1.631390e-01
* time: 7.400522947311401
115 1.888378e+03 2.768731e-01
* time: 7.439647912979126
116 1.888377e+03 4.462386e-01
* time: 7.479233980178833
117 1.888377e+03 6.643297e-01
* time: 7.519218921661377
118 1.888375e+03 8.433397e-01
* time: 7.559473037719727
119 1.888374e+03 7.596814e-01
* time: 7.5995728969573975
120 1.888373e+03 3.638119e-01
* time: 7.640797853469849
121 1.888372e+03 8.306034e-02
* time: 7.686378002166748
122 1.888372e+03 2.084513e-02
* time: 7.7286388874053955
123 1.888372e+03 2.056419e-02
* time: 7.796108961105347
124 1.888372e+03 2.044080e-02
* time: 7.834166049957275
125 1.888372e+03 2.035196e-02
* time: 7.870481014251709
126 1.888372e+03 2.021264e-02
* time: 7.906515836715698
127 1.888372e+03 1.998163e-02
* time: 7.942703008651733
128 1.888372e+03 3.161638e-02
* time: 7.979389905929565
129 1.888372e+03 5.509218e-02
* time: 8.01636004447937
130 1.888372e+03 9.275848e-02
* time: 8.053524017333984
131 1.888372e+03 1.528749e-01
* time: 8.091158866882324
132 1.888372e+03 2.461766e-01
* time: 8.130915880203247
133 1.888372e+03 3.799362e-01
* time: 8.173759937286377
134 1.888371e+03 5.311665e-01
* time: 8.219302892684937
135 1.888369e+03 6.019039e-01
* time: 8.293083906173706
136 1.888367e+03 4.664896e-01
* time: 8.33390188217163
137 1.888366e+03 1.404934e-01
* time: 8.374315977096558
138 1.888365e+03 8.513331e-02
* time: 8.414265871047974
139 1.888364e+03 1.244077e-01
* time: 8.452828884124756
140 1.888364e+03 1.028085e-01
* time: 8.4911949634552
141 1.888364e+03 5.162231e-02
* time: 8.530449867248535
142 1.888364e+03 5.149630e-02
* time: 8.569143056869507
143 1.888364e+03 3.147284e-02
* time: 8.60775899887085
144 1.888364e+03 2.104766e-02
* time: 8.64875602722168
145 1.888364e+03 6.539151e-03
* time: 8.688612937927246
146 1.888364e+03 2.540196e-03
* time: 8.728418827056885
147 1.888364e+03 4.362661e-03
* time: 8.790963888168335
148 1.888364e+03 3.034416e-03
* time: 8.827704906463623
149 1.888364e+03 5.953892e-04
* time: 8.863325834274292FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -1888.3638
Number of subjects: 30
Number of parameters: Fixed Optimized
0 13
Observation records: Active Missing
DV: 540 0
Total: 540 0
--------------------------
Estimate
--------------------------
tvvc 3.976
tvq 0.04239
tvvp 3.7249
tvcl_hep_M 1.7175e-7
tvcl_hep_F 0.13348
tvcl_ren_M 0.19378
tvcl_ren_F 0.042211
Ω₁,₁ 0.14046
Ω₂,₂ 0.081349
σ_add 0.032171
σ_prop 0.061007
dCRCL_M 0.94821
dCRCL_F 1.9405
--------------------------As before, our loglikelihood is higher (implying lower objective function value). This is expected since we also added six new parameters to the model.
1.5 Step 4 - Model Comparison
Now that we’ve fitted all of our models we need to compare them and choose one for our final model.
We begin by analyzing the model metrics. This can be done with the metrics_table function:
metrics_table(fit_base_model)17×2 DataFrame
| Row | Metric | Value |
|---|---|---|
| String | Any | |
| 1 | Successful | true |
| 2 | Estimation Time | 3.251 |
| 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.564 |
| 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.885 |
| 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 | 8.863 |
| 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.251 | 2.564 | 4.885 | 8.863 |
| 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.629984952819aic(fit_covariate_model_wt)3808.9264607805967aic(fit_covariate_model_wt_crcl)3804.817947371705aic(fit_covariate_model_wt_crcl_sex)3802.7275243741956We should favor lower values of AIC, hence, the covariate model with weight, creatinine clerance, and different sex effects on clearance should be preferred, i.e. covariate_model_wt_crcl_sex.
1.5.1 Goodness of Fit Plots
Additionally, we should inspect the goodness of fit of the model. This is done with the plotting function goodness_of_fit, which should be given a result from a inspect function. So, let’s first call inspect on our covariate_model_wt_crcl_sex candidate for best model:
inspect_covariate_model_wt_crcl_sex = inspect(fit_covariate_model_wt_crcl_sex)[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating individual parameters.
[ Info: Evaluating dose control parameters.
[ Info: Done.FittedPumasModelInspection
Fitting was successful: true
Likehood approximation used for weighted residuals : Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}(Optim.NewtonTrustRegion{Float64}(1.0, 100.0, 1.4901161193847656e-8, 0.1, 0.25, 0.75, false), Optim.Options(x_abstol = 0.0, x_reltol = 0.0, f_abstol = 0.0, f_reltol = 0.0, g_abstol = 1.0e-5, g_reltol = 1.0e-8, outer_x_abstol = 0.0, outer_x_reltol = 0.0, outer_f_abstol = 0.0, outer_f_reltol = 0.0, outer_g_abstol = 1.0e-8, outer_g_reltol = 1.0e-8, f_calls_limit = 0, g_calls_limit = 0, h_calls_limit = 0, allow_f_increases = false, allow_outer_f_increases = true, successive_f_tol = 1, iterations = 1000, outer_iterations = 1000, store_trace = false, trace_simplex = false, show_trace = false, extended_trace = false, show_warnings = true, show_every = 1, time_limit = NaN, )
)And now we pass inspect_covariate_model_wt_crcl_sex to the goodness_of_fit plotting function:
goodness_of_fit(inspect_covariate_model_wt_crcl_sex)
The idea is that the population predictions (preds) capture the general tendency of the observations while the individual predictions (ipreds) should coincide much more closely with the observations. That is exactly what we are observing in the top row subplots in our goodness of fit plot.
Regarding the bottom row, on the left we have the weighted population residuals (wres) against time, and on the right we have the weighted individual residuals (iwres) against ipreds. Here we should not see any perceived pattern, which indicates that the error model in the model has a mean 0 and constant variance. Like before, this seems to be what we are observing in our goodness of fit plot.
Hence, our covariate model with allometric scaling and different sex creatinine clearance effectw on clearance is a good candidate for our final model.
1.6 Time-Varying Covariates
Pumas can handle time-varying covariates. This happens automatically if, when parsing a dataset, read_pumas detects that covariate values change over time.
1.6.1 painord Dataset
Here’s an example with an ordinal regression using the painord dataset from PharmaDatasets.jl. :painord is our observations measuring the perceived pain in a scale from 0 to 3, which we need to have the range shifted by 1 (1 to 4). Additionally, we’ll use the concentration in plasma, :conc as a covariate. Of course, :conc varies with time, thus, it is a time-varying covariate:
painord = dataset("pumas/pain_remed")
first(painord, 5)5×8 DataFrame
| Row | id | arm | dose | time | conc | painord | dv | remed |
|---|---|---|---|---|---|---|---|---|
| Int64 | Int64 | Int64 | Float64 | Float64 | Int64 | Int64 | Int64 | |
| 1 | 1 | 2 | 20 | 0.0 | 0.0 | 3 | 0 | 0 |
| 2 | 1 | 2 | 20 | 0.5 | 1.15578 | 1 | 1 | 0 |
| 3 | 1 | 2 | 20 | 1.0 | 1.37211 | 0 | 1 | 0 |
| 4 | 1 | 2 | 20 | 1.5 | 1.30058 | 0 | 1 | 0 |
| 5 | 1 | 2 | 20 | 2.0 | 1.19195 | 1 | 1 | 0 |
@rtransform! painord :painord = :painord + 1;
describe(painord, :mean, :std, :first, :last, :eltype)8×6 DataFrame
| Row | variable | mean | std | first | last | eltype |
|---|---|---|---|---|---|---|
| Symbol | Float64 | Float64 | Real | Real | DataType | |
| 1 | id | 80.5 | 46.1992 | 1 | 160 | Int64 |
| 2 | arm | 1.5 | 1.11833 | 2 | 0 | Int64 |
| 3 | dose | 26.25 | 31.9017 | 20 | 0 | Int64 |
| 4 | time | 3.375 | 2.5183 | 0.0 | 8.0 | Float64 |
| 5 | conc | 0.93018 | 1.49902 | 0.0 | 0.0 | Float64 |
| 6 | painord | 2.50208 | 0.863839 | 4 | 4 | Int64 |
| 7 | dv | 0.508333 | 0.500061 | 0 | 0 | Int64 |
| 8 | remed | 0.059375 | 0.236387 | 0 | 0 | Int64 |
unique(painord.dose)4-element Vector{Int64}:
20
80
0
5As 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
endPumasModel
Parameters: b₁, b₂, b₃, slope
Random effects:
Covariates: conc
Dynamical system variables:
Dynamical system type: No dynamical model
Derived: painord
Observed: painordFinally 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: 7.295608520507812e-5
1 2.994747e+03 1.083462e+03
* time: 0.004210948944091797
2 2.406265e+03 1.884408e+02
* time: 0.008188009262084961
3 2.344175e+03 7.741610e+01
* time: 0.012206792831420898
4 2.323153e+03 2.907642e+01
* time: 0.016251802444458008
5 2.318222e+03 2.273295e+01
* time: 0.020190954208374023
6 2.316833e+03 1.390527e+01
* time: 0.024205923080444336
7 2.316425e+03 4.490883e+00
* time: 0.028472900390625
8 2.316362e+03 9.374519e-01
* time: 0.03273797035217285
9 2.316356e+03 1.928785e-01
* time: 0.03696393966674805
10 2.316355e+03 3.119615e-02
* time: 0.04159998893737793
11 2.316355e+03 6.215513e-03
* time: 0.04570198059082031
12 2.316355e+03 8.313206e-04
* time: 0.05024600028991699FittedPumasModel
Successful minimization: true
Likelihood approximation: NaivePooled
Likelihood Optimizer: BFGS
Dynamical system type: No dynamical model
Log-likelihood value: -2316.3554
Number of subjects: 160
Number of parameters: Fixed Optimized
0 4
Observation records: Active Missing
painord: 1920 0
Total: 1920 0
-------------------
Estimate
-------------------
b₁ 2.5112
b₂ 2.1951
b₃ 1.9643
slope -0.38871
-------------------As expected, the ordinal model fit estimates a negative effect of :conc on :painord measured by the slope parameter.
1.7 Missing Data in Covariates
The way how Pumas handles missing values inside covariates depends if the covariate is constant or time-varying. For both cases Pumas will interpolate the available values to fill in the missing values. If, for any subject, all of the covariate’s values are missing, Pumas will throw an error while parsing the data with read_pumas.
For both missing constant and time-varying covariates, Pumas, by default, does piece-wise constant interpolation with “next observation carried backward” (NOCB, NONMEM default). Of course for constant covariates the interpolated values over the missing values will be constant values. This can be adjusted with the covariates_direction keyword argument of read_pumas. The default value :right is NOCB and :left is “last observation carried forward” (LOCF, Monolix default).
Hence, for LOCF, you can use the following:
pop = read_pumas(pkdata; covariates_direction = :left)along with any other required keyword arguments for column mapping.
Note
The same behavior for covariates_direction applies to time-varying covariates during the interpolation in the ODE solver. They will also be piece-wise constant interpolated following either NOCB or LOCF depending on the covariates_direction value.
1.8 Categorical Covariates
In some situations, you’ll find yourself with categorical covariates with multiple levels, instead of binary or continuous covariates. Categorical covariates are covariates that can take on a finite number of distinct values.
Pumas can easily address categorical covariates. In the @pre block you can use a nested if ... elseif ... else statement to handle the different categories.
For example:
@pre begin
CL = if RACE == 1
tvcl * (WT / 70)^dwtdcl * exp(η[1]) * drace1dcl
elseif RACE == 2
tvcl * (WT / 70)^dwtdcl * exp(η[1]) * drace2dcl
elseif RACE == 3
tvcl * (WT / 70)^dwtdcl * exp(η[1]) * drace3dcl
end
endHere 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
endHere we are using the ^ operator to raise each race value to the power of the race1cl, race2cl, and race3cl values. If any of the race values is not equal to the race value, the result will be 1, otherwise it will be the respective race1cl, race2cl, or race3cl value.
1.9 Conclusion
This tutorial shows how to build covariate model in Pumas in a workflow approach. The main purpose was to inform how to:
- parse covariate data into a
Population - add covariate information into a model
We also went over what are the differences between constant and time-varying covariates and how does Pumas deal with missing data inside covariates.