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.02935004234313965
1 3.110315e+03 9.706222e+02
* time: 1.1769850254058838
2 2.831659e+03 7.817006e+02
* time: 1.2634670734405518
3 2.405281e+03 2.923696e+02
* time: 1.3285319805145264
4 2.370406e+03 3.032286e+02
* time: 1.3738369941711426
5 2.313631e+03 3.126188e+02
* time: 1.4200849533081055
6 2.263986e+03 2.946697e+02
* time: 1.4630460739135742
7 2.160182e+03 1.917599e+02
* time: 1.533384084701538
8 2.112467e+03 1.588288e+02
* time: 1.6376910209655762
9 2.090339e+03 1.108334e+02
* time: 1.6771059036254883
10 2.078171e+03 8.108027e+01
* time: 1.7159090042114258
11 2.074517e+03 7.813928e+01
* time: 1.7511959075927734
12 2.066270e+03 7.044632e+01
* time: 1.7862639427185059
13 2.049660e+03 1.062584e+02
* time: 1.821418046951294
14 2.021965e+03 1.130570e+02
* time: 1.858064889907837
15 1.994936e+03 7.825801e+01
* time: 1.8950130939483643
16 1.979337e+03 5.318263e+01
* time: 1.935333013534546
17 1.972141e+03 6.807046e+01
* time: 1.9796230792999268
18 1.967973e+03 7.896361e+01
* time: 2.0380210876464844
19 1.962237e+03 8.343757e+01
* time: 2.0758769512176514
20 1.952791e+03 5.565304e+01
* time: 2.114711046218872
21 1.935857e+03 3.923284e+01
* time: 2.153898000717163
22 1.926254e+03 5.749643e+01
* time: 2.192288875579834
23 1.922144e+03 4.306225e+01
* time: 2.2283780574798584
24 1.911566e+03 4.810496e+01
* time: 2.2677738666534424
25 1.906464e+03 4.324267e+01
* time: 2.3067610263824463
26 1.905339e+03 1.207954e+01
* time: 2.3413848876953125
27 1.905092e+03 1.093479e+01
* time: 2.3747599124908447
28 1.904957e+03 1.057034e+01
* time: 2.409632921218872
29 1.904875e+03 1.060882e+01
* time: 2.4610989093780518
30 1.904459e+03 1.031525e+01
* time: 2.4972469806671143
31 1.903886e+03 1.041179e+01
* time: 2.5325279235839844
32 1.903313e+03 1.135672e+01
* time: 2.567986011505127
33 1.903057e+03 1.075683e+01
* time: 2.6014559268951416
34 1.902950e+03 1.091295e+01
* time: 2.635232925415039
35 1.902887e+03 1.042409e+01
* time: 2.667941093444824
36 1.902640e+03 9.213043e+00
* time: 2.7020180225372314
37 1.902364e+03 9.519321e+00
* time: 2.736999034881592
38 1.902156e+03 5.590984e+00
* time: 2.7715940475463867
39 1.902094e+03 1.811898e+00
* time: 2.806452989578247
40 1.902086e+03 1.644722e+00
* time: 2.857564926147461
41 1.902084e+03 1.653520e+00
* time: 2.8898918628692627
42 1.902074e+03 1.720184e+00
* time: 2.922653913497925
43 1.902055e+03 2.619061e+00
* time: 2.9552340507507324
44 1.902015e+03 3.519729e+00
* time: 2.9891300201416016
45 1.901962e+03 3.403372e+00
* time: 3.021714925765991
46 1.901924e+03 1.945644e+00
* time: 3.0543620586395264
47 1.901914e+03 6.273342e-01
* time: 3.0869898796081543
48 1.901913e+03 5.374557e-01
* time: 3.1189348697662354
49 1.901913e+03 5.710007e-01
* time: 3.1495790481567383
50 1.901913e+03 6.091390e-01
* time: 3.1804840564727783
51 1.901912e+03 6.692417e-01
* time: 3.2299580574035645
52 1.901909e+03 7.579620e-01
* time: 3.2625069618225098
53 1.901903e+03 8.798211e-01
* time: 3.2947909832000732
54 1.901889e+03 1.002981e+00
* time: 3.3270270824432373
55 1.901864e+03 1.495512e+00
* time: 3.359454870223999
56 1.901837e+03 1.664621e+00
* time: 3.3905060291290283
57 1.901819e+03 8.601118e-01
* time: 3.4224190711975098
58 1.901815e+03 4.525470e-02
* time: 3.45468807220459
59 1.901815e+03 1.294280e-02
* time: 3.4860939979553223
60 1.901815e+03 2.876568e-03
* time: 3.5148510932922363
61 1.901815e+03 2.876568e-03
* time: 3.5974059104919434
62 1.901815e+03 2.876568e-03
* time: 3.6968939304351807FittedPumasModel
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: 6.103515625e-5
1 2.846050e+03 1.128657e+03
* time: 0.07121706008911133
2 2.472982e+03 7.008264e+02
* time: 0.11736798286437988
3 2.180690e+03 2.312704e+02
* time: 0.16271114349365234
4 2.125801e+03 1.862929e+02
* time: 0.20327496528625488
5 2.103173e+03 1.320946e+02
* time: 0.2435779571533203
6 2.091136e+03 1.103035e+02
* time: 0.32127809524536133
7 2.081443e+03 1.091133e+02
* time: 0.35655808448791504
8 2.071793e+03 7.197675e+01
* time: 0.392578125
9 2.062706e+03 7.623310e+01
* time: 0.43067193031311035
10 2.057515e+03 6.885476e+01
* time: 0.4664490222930908
11 2.051133e+03 6.368504e+01
* time: 0.5010969638824463
12 2.038626e+03 7.730243e+01
* time: 0.535750150680542
13 2.019352e+03 1.136864e+02
* time: 0.5711359977722168
14 1.997136e+03 1.005748e+02
* time: 0.6070339679718018
15 1.983023e+03 6.831478e+01
* time: 0.6450009346008301
16 1.977700e+03 5.649783e+01
* time: 0.6816380023956299
17 1.974583e+03 6.357642e+01
* time: 0.7499029636383057
18 1.967292e+03 7.658974e+01
* time: 0.7880001068115234
19 1.951045e+03 6.130573e+01
* time: 0.8325300216674805
20 1.935868e+03 4.845839e+01
* time: 0.871906042098999
21 1.929356e+03 6.325782e+01
* time: 0.9122159481048584
22 1.925187e+03 3.142245e+01
* time: 0.950484037399292
23 1.923733e+03 4.623400e+01
* time: 0.9869439601898193
24 1.918498e+03 5.347738e+01
* time: 1.0244760513305664
25 1.912383e+03 5.849125e+01
* time: 1.0669541358947754
26 1.905510e+03 3.254038e+01
* time: 1.111915111541748
27 1.903629e+03 2.905618e+01
* time: 1.1695449352264404
28 1.902833e+03 2.907696e+01
* time: 1.2066810131072998
29 1.902447e+03 2.746037e+01
* time: 1.2401659488677979
30 1.899399e+03 1.930949e+01
* time: 1.2772619724273682
31 1.898705e+03 1.186361e+01
* time: 1.313758134841919
32 1.898505e+03 1.050402e+01
* time: 1.3499641418457031
33 1.898474e+03 1.042186e+01
* time: 1.3812201023101807
34 1.897992e+03 1.238729e+01
* time: 1.4143741130828857
35 1.897498e+03 1.729368e+01
* time: 1.4484729766845703
36 1.896954e+03 1.472554e+01
* time: 1.4816229343414307
37 1.896744e+03 5.852825e+00
* time: 1.517406940460205
38 1.896705e+03 1.171353e+00
* time: 1.5659239292144775
39 1.896704e+03 1.216117e+00
* time: 1.5972459316253662
40 1.896703e+03 1.230336e+00
* time: 1.6280810832977295
41 1.896698e+03 1.250715e+00
* time: 1.6588480472564697
42 1.896688e+03 1.449552e+00
* time: 1.68961501121521
43 1.896666e+03 2.533300e+00
* time: 1.7208261489868164
44 1.896631e+03 3.075536e+00
* time: 1.7521800994873047
45 1.896599e+03 2.139797e+00
* time: 1.7832670211791992
46 1.896587e+03 6.636031e-01
* time: 1.8157551288604736
47 1.896585e+03 6.303517e-01
* time: 1.8475430011749268
48 1.896585e+03 5.995265e-01
* time: 1.876878023147583
49 1.896584e+03 5.844159e-01
* time: 1.9260039329528809
50 1.896583e+03 6.083858e-01
* time: 1.9573280811309814
51 1.896579e+03 8.145326e-01
* time: 1.9879429340362549
52 1.896570e+03 1.293490e+00
* time: 2.019066095352173
53 1.896549e+03 1.877705e+00
* time: 2.05019211769104
54 1.896513e+03 2.217391e+00
* time: 2.081171989440918
55 1.896482e+03 1.658147e+00
* time: 2.111686944961548
56 1.896466e+03 5.207215e-01
* time: 2.1434619426727295
57 1.896463e+03 1.177468e-01
* time: 2.175405979156494
58 1.896463e+03 1.678937e-02
* time: 2.2050631046295166
59 1.896463e+03 2.666635e-03
* time: 2.232412099838257
60 1.896463e+03 2.666635e-03
* time: 2.2954370975494385FittedPumasModel
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: 5.602836608886719e-5
1 3.572050e+03 1.302046e+03
* time: 0.10300612449645996
2 3.266947e+03 5.384244e+02
* time: 0.16959404945373535
3 3.150635e+03 1.918079e+02
* time: 0.22528409957885742
4 3.108122e+03 1.277799e+02
* time: 0.34724998474121094
5 3.091358e+03 8.838080e+01
* time: 0.39362001419067383
6 3.082997e+03 8.321689e+01
* time: 0.43983006477355957
7 3.076379e+03 8.167702e+01
* time: 0.4876370429992676
8 3.067422e+03 1.177822e+02
* time: 0.5364120006561279
9 3.048580e+03 2.526969e+02
* time: 0.5916790962219238
10 2.993096e+03 6.325396e+02
* time: 0.6786770820617676
11 2.965744e+03 7.634718e+02
* time: 0.9084429740905762
12 2.921212e+03 9.704020e+02
* time: 1.0079030990600586
13 2.553649e+03 6.642510e+02
* time: 1.0840721130371094
14 2.319495e+03 3.204711e+02
* time: 1.2045221328735352
15 2.183040e+03 2.174905e+02
* time: 1.3020451068878174
16 2.157621e+03 3.150983e+02
* time: 1.356050968170166
17 2.132395e+03 2.847614e+02
* time: 1.4076321125030518
18 2.084799e+03 1.563370e+02
* time: 1.4855740070343018
19 2.071497e+03 1.006429e+02
* time: 1.5309531688690186
20 2.064983e+03 9.753313e+01
* time: 1.578489065170288
21 2.048289e+03 9.230405e+01
* time: 1.6272411346435547
22 2.020938e+03 1.292359e+02
* time: 1.6828770637512207
23 1.983888e+03 2.284990e+02
* time: 1.7361149787902832
24 1.962132e+03 1.220188e+02
* time: 1.7900409698486328
25 1.945947e+03 1.035894e+02
* time: 1.844547986984253
26 1.917782e+03 8.246698e+01
* time: 1.9059619903564453
27 1.905967e+03 5.408054e+01
* time: 2.011463165283203
28 1.898569e+03 2.172222e+01
* time: 2.082237958908081
29 1.897473e+03 1.689350e+01
* time: 2.165977954864502
30 1.897019e+03 1.666689e+01
* time: 2.2126760482788086
31 1.896796e+03 1.699751e+01
* time: 2.2585320472717285
32 1.896559e+03 1.645900e+01
* time: 2.304208993911743
33 1.896223e+03 1.415504e+01
* time: 2.3487319946289062
34 1.895773e+03 1.630081e+01
* time: 2.3929941654205322
35 1.895309e+03 1.723930e+01
* time: 2.441127061843872
36 1.895004e+03 1.229983e+01
* time: 2.490797996520996
37 1.894871e+03 5.385102e+00
* time: 2.5393850803375244
38 1.894827e+03 3.465463e+00
* time: 2.5885801315307617
39 1.894816e+03 3.387474e+00
* time: 2.662961006164551
40 1.894807e+03 3.295388e+00
* time: 2.7045440673828125
41 1.894786e+03 3.089194e+00
* time: 2.7468440532684326
42 1.894737e+03 2.928080e+00
* time: 2.789769172668457
43 1.894624e+03 3.088723e+00
* time: 2.8348071575164795
44 1.894413e+03 3.493791e+00
* time: 2.87839412689209
45 1.894127e+03 3.142865e+00
* time: 2.9214179515838623
46 1.893933e+03 2.145253e+00
* time: 2.967803955078125
47 1.893888e+03 2.172800e+00
* time: 3.0156800746917725
48 1.893880e+03 2.180509e+00
* time: 3.061156988143921
49 1.893876e+03 2.134257e+00
* time: 3.106034994125366
50 1.893868e+03 2.032354e+00
* time: 3.1510660648345947
51 1.893846e+03 1.760874e+00
* time: 3.226155996322632
52 1.893796e+03 1.779016e+00
* time: 3.267941951751709
53 1.893694e+03 2.018956e+00
* time: 3.3100521564483643
54 1.893559e+03 2.366854e+00
* time: 3.3517560958862305
55 1.893474e+03 3.690142e+00
* time: 3.3935561180114746
56 1.893446e+03 3.675109e+00
* time: 3.4343600273132324
57 1.893439e+03 3.426419e+00
* time: 3.4769749641418457
58 1.893429e+03 3.183164e+00
* time: 3.522099018096924
59 1.893398e+03 2.695171e+00
* time: 3.566749095916748
60 1.893328e+03 2.753548e+00
* time: 3.6116859912872314
61 1.893169e+03 3.589748e+00
* time: 3.658092975616455
62 1.892920e+03 3.680718e+00
* time: 3.728778123855591
63 1.892667e+03 2.568107e+00
* time: 3.7706100940704346
64 1.892514e+03 1.087910e+00
* time: 3.8120369911193848
65 1.892493e+03 3.287296e-01
* time: 3.8526649475097656
66 1.892492e+03 2.967465e-01
* time: 3.8928849697113037
67 1.892492e+03 3.020682e-01
* time: 3.931666135787964
68 1.892491e+03 3.034704e-01
* time: 3.970111131668091
69 1.892491e+03 3.091846e-01
* time: 4.0128560066223145
70 1.892491e+03 3.224170e-01
* time: 4.0545690059661865
71 1.892490e+03 6.494197e-01
* time: 4.097576141357422
72 1.892488e+03 1.115188e+00
* time: 4.141026020050049
73 1.892483e+03 1.838833e+00
* time: 4.214269161224365
74 1.892472e+03 2.765371e+00
* time: 4.254464149475098
75 1.892452e+03 3.463807e+00
* time: 4.294823169708252
76 1.892431e+03 2.805270e+00
* time: 4.334831953048706
77 1.892411e+03 5.758916e-01
* time: 4.37568998336792
78 1.892410e+03 1.434041e-01
* time: 4.414981126785278
79 1.892409e+03 1.639246e-01
* time: 4.4541871547698975
80 1.892409e+03 1.145856e-01
* time: 4.497282981872559
81 1.892409e+03 3.966861e-02
* time: 4.538539171218872
82 1.892409e+03 3.550808e-02
* time: 4.578602075576782
83 1.892409e+03 3.456241e-02
* time: 4.616714000701904
84 1.892409e+03 3.114018e-02
* time: 4.6818201541900635
85 1.892409e+03 4.080806e-02
* time: 4.718448162078857
86 1.892409e+03 6.722726e-02
* time: 4.755295038223267
87 1.892409e+03 1.006791e-01
* time: 4.791882038116455
88 1.892409e+03 1.303988e-01
* time: 4.8290581703186035
89 1.892409e+03 1.228919e-01
* time: 4.865770101547241
90 1.892409e+03 6.433813e-02
* time: 4.90289306640625
91 1.892409e+03 1.314164e-02
* time: 4.94144606590271
92 1.892409e+03 4.929931e-04
* time: 4.98210597038269FittedPumasModel
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: 7.510185241699219e-5
1 3.641387e+03 1.432080e+03
* time: 0.1184539794921875
2 3.290450e+03 5.274782e+02
* time: 0.33504700660705566
3 3.185512e+03 2.173676e+02
* time: 0.3999791145324707
4 3.143009e+03 1.479653e+02
* time: 0.46003198623657227
5 3.128511e+03 8.980031e+01
* time: 0.5174679756164551
6 3.123188e+03 5.033125e+01
* time: 0.5743131637573242
7 3.120794e+03 4.279722e+01
* time: 0.6425201892852783
8 3.118627e+03 3.971051e+01
* time: 0.7150530815124512
9 3.115300e+03 8.456587e+01
* time: 0.7919900417327881
10 3.109353e+03 1.350354e+02
* time: 0.8742451667785645
11 3.095894e+03 1.998258e+02
* time: 0.9604899883270264
12 2.988214e+03 4.366433e+02
* time: 1.1381890773773193
13 2.896081e+03 5.505943e+02
* time: 1.3841910362243652
14 2.652467e+03 7.300323e+02
* time: 3.040929079055786
15 2.560937e+03 6.973661e+02
* time: 3.3954989910125732
16 2.254941e+03 2.740033e+02
* time: 3.4765830039978027
17 2.222509e+03 2.034303e+02
* time: 3.549384117126465
18 2.171255e+03 2.449580e+02
* time: 3.6287879943847656
19 2.024532e+03 1.121511e+02
* time: 3.7054200172424316
20 1.993723e+03 1.042814e+02
* time: 3.7758371829986572
21 1.985113e+03 8.079014e+01
* time: 3.887111186981201
22 1.976757e+03 7.054196e+01
* time: 3.9451889991760254
23 1.969970e+03 6.070322e+01
* time: 4.0024189949035645
24 1.961095e+03 6.810782e+01
* time: 4.0570220947265625
25 1.947983e+03 8.116920e+01
* time: 4.112661123275757
26 1.930371e+03 8.530051e+01
* time: 4.168943166732788
27 1.910209e+03 6.993170e+01
* time: 4.225691080093384
28 1.899107e+03 3.362640e+01
* time: 4.2855870723724365
29 1.898022e+03 2.642220e+01
* time: 4.351611137390137
30 1.897055e+03 1.213144e+01
* time: 4.418753147125244
31 1.896596e+03 7.773239e+00
* time: 4.486931085586548
32 1.896538e+03 7.997039e+00
* time: 4.576246023178101
33 1.896451e+03 8.160909e+00
* time: 4.62873911857605
34 1.896283e+03 8.237721e+00
* time: 4.684401988983154
35 1.895903e+03 1.520219e+01
* time: 4.742570161819458
36 1.895272e+03 2.358916e+01
* time: 4.801696062088013
37 1.894536e+03 2.461296e+01
* time: 4.860040187835693
38 1.893995e+03 1.546128e+01
* time: 4.916942119598389
39 1.893858e+03 6.976137e+00
* time: 4.976181983947754
40 1.893833e+03 6.019466e+00
* time: 5.03357720375061
41 1.893786e+03 3.827201e+00
* time: 5.093119144439697
42 1.893714e+03 3.323412e+00
* time: 5.18989896774292
43 1.893592e+03 3.215150e+00
* time: 5.247416973114014
44 1.893435e+03 6.534965e+00
* time: 5.305582046508789
45 1.893286e+03 7.424154e+00
* time: 5.3625781536102295
46 1.893190e+03 5.552627e+00
* time: 5.409389972686768
47 1.893139e+03 3.222316e+00
* time: 5.455549001693726
48 1.893120e+03 3.015339e+00
* time: 5.5014941692352295
49 1.893107e+03 3.244809e+00
* time: 5.546709060668945
50 1.893080e+03 6.163100e+00
* time: 5.60087513923645
51 1.893027e+03 9.824713e+00
* time: 5.6551690101623535
52 1.892912e+03 1.390100e+01
* time: 5.710464000701904
53 1.892734e+03 1.510937e+01
* time: 5.798603057861328
54 1.892561e+03 1.008563e+01
* time: 5.853283166885376
55 1.892485e+03 3.730668e+00
* time: 5.9068989753723145
56 1.892471e+03 3.380261e+00
* time: 5.9567461013793945
57 1.892463e+03 3.167904e+00
* time: 6.003578186035156
58 1.892441e+03 4.152065e+00
* time: 6.050342082977295
59 1.892391e+03 7.355996e+00
* time: 6.098505020141602
60 1.892268e+03 1.195397e+01
* time: 6.146061182022095
61 1.892026e+03 1.640783e+01
* time: 6.197158098220825
62 1.891735e+03 1.593576e+01
* time: 6.249678134918213
63 1.891569e+03 8.316423e+00
* time: 6.302676200866699
64 1.891494e+03 3.948212e+00
* time: 6.387164115905762
65 1.891481e+03 3.911593e+00
* time: 6.435155153274536
66 1.891457e+03 3.875559e+00
* time: 6.483137130737305
67 1.891405e+03 3.811247e+00
* time: 6.532414197921753
68 1.891262e+03 3.657045e+00
* time: 6.580892086029053
69 1.890930e+03 4.957405e+00
* time: 6.631775140762329
70 1.890317e+03 6.657726e+00
* time: 6.683815002441406
71 1.889660e+03 6.086302e+00
* time: 6.73779296875
72 1.889303e+03 2.270929e+00
* time: 6.792219161987305
73 1.889253e+03 7.695301e-01
* time: 6.851237058639526
74 1.889252e+03 7.382144e-01
* time: 6.907608985900879
75 1.889251e+03 7.187898e-01
* time: 6.987549066543579
76 1.889251e+03 7.215047e-01
* time: 7.036225080490112
77 1.889250e+03 7.235155e-01
* time: 7.086263179779053
78 1.889249e+03 7.246818e-01
* time: 7.13473105430603
79 1.889244e+03 7.257796e-01
* time: 7.183533191680908
80 1.889233e+03 7.198190e-01
* time: 7.232514142990112
81 1.889204e+03 1.089029e+00
* time: 7.282849073410034
82 1.889142e+03 1.801601e+00
* time: 7.335432052612305
83 1.889043e+03 2.967917e+00
* time: 7.391189098358154
84 1.888889e+03 2.965856e+00
* time: 7.449195146560669
85 1.888705e+03 5.933557e-01
* time: 7.507081031799316
86 1.888655e+03 9.577696e-01
* time: 7.595210075378418
87 1.888582e+03 1.498494e+00
* time: 7.646392107009888
88 1.888533e+03 1.502753e+00
* time: 7.696767091751099
89 1.888490e+03 1.184664e+00
* time: 7.745906114578247
90 1.888480e+03 6.684517e-01
* time: 7.793932199478149
91 1.888476e+03 3.680034e-01
* time: 7.84133505821228
92 1.888476e+03 4.720040e-01
* time: 7.887675046920776
93 1.888476e+03 4.768646e-01
* time: 7.933869123458862
94 1.888475e+03 4.736673e-01
* time: 7.982248067855835
95 1.888475e+03 4.552764e-01
* time: 8.03174614906311
96 1.888474e+03 5.193733e-01
* time: 8.080312967300415
97 1.888473e+03 8.850112e-01
* time: 8.13499116897583
98 1.888468e+03 1.461600e+00
* time: 8.217288970947266
99 1.888458e+03 2.209127e+00
* time: 8.268439054489136
100 1.888437e+03 2.961239e+00
* time: 8.32102108001709
101 1.888407e+03 2.978463e+00
* time: 8.373000144958496
102 1.888384e+03 1.707201e+00
* time: 8.424444198608398
103 1.888381e+03 6.199354e-01
* time: 8.476975202560425
104 1.888380e+03 5.170909e-01
* time: 8.529445171356201
105 1.888378e+03 1.037408e-01
* time: 8.584989070892334
106 1.888378e+03 8.473247e-02
* time: 8.640403032302856
107 1.888378e+03 8.364978e-02
* time: 8.69565200805664
108 1.888378e+03 8.080446e-02
* time: 8.751947164535522
109 1.888378e+03 7.873905e-02
* time: 8.839109182357788
110 1.888378e+03 7.798398e-02
* time: 8.885944128036499
111 1.888378e+03 7.788170e-02
* time: 8.937448978424072
112 1.888378e+03 7.776461e-02
* time: 8.988469123840332
113 1.888378e+03 9.023662e-02
* time: 9.035887002944946
114 1.888378e+03 1.631390e-01
* time: 9.08444619178772
115 1.888378e+03 2.768731e-01
* time: 9.134692192077637
116 1.888377e+03 4.462386e-01
* time: 9.183760166168213
117 1.888377e+03 6.643297e-01
* time: 9.238860130310059
118 1.888375e+03 8.433397e-01
* time: 9.294308185577393
119 1.888374e+03 7.596814e-01
* time: 9.349895000457764
120 1.888373e+03 3.638119e-01
* time: 9.43828797340393
121 1.888372e+03 8.306034e-02
* time: 9.488901138305664
122 1.888372e+03 2.084513e-02
* time: 9.534430027008057
123 1.888372e+03 2.056419e-02
* time: 9.579140186309814
124 1.888372e+03 2.044080e-02
* time: 9.623831033706665
125 1.888372e+03 2.035196e-02
* time: 9.665657043457031
126 1.888372e+03 2.021264e-02
* time: 9.708340167999268
127 1.888372e+03 1.998163e-02
* time: 9.75341010093689
128 1.888372e+03 3.161638e-02
* time: 9.803461074829102
129 1.888372e+03 5.509218e-02
* time: 9.854643106460571
130 1.888372e+03 9.275848e-02
* time: 9.905915021896362
131 1.888372e+03 1.528749e-01
* time: 9.957207202911377
132 1.888372e+03 2.461766e-01
* time: 10.039924144744873
133 1.888372e+03 3.799362e-01
* time: 10.085792064666748
134 1.888371e+03 5.311665e-01
* time: 10.132667064666748
135 1.888369e+03 6.019039e-01
* time: 10.178357124328613
136 1.888367e+03 4.664896e-01
* time: 10.22363018989563
137 1.888366e+03 1.404934e-01
* time: 10.268396139144897
138 1.888365e+03 8.513331e-02
* time: 10.31547999382019
139 1.888364e+03 1.244077e-01
* time: 10.36046814918518
140 1.888364e+03 1.028085e-01
* time: 10.405840158462524
141 1.888364e+03 5.162231e-02
* time: 10.45581603050232
142 1.888364e+03 5.149630e-02
* time: 10.505033016204834
143 1.888364e+03 3.147284e-02
* time: 10.584568977355957
144 1.888364e+03 2.104766e-02
* time: 10.631473064422607
145 1.888364e+03 6.539151e-03
* time: 10.67688512802124
146 1.888364e+03 2.540196e-03
* time: 10.721837997436523
147 1.888364e+03 4.362661e-03
* time: 10.764794111251831
148 1.888364e+03 3.034416e-03
* time: 10.809429168701172
149 1.888364e+03 5.953892e-04
* time: 10.854665040969849FittedPumasModel
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.697 |
| 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.296 |
| 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.982 |
| 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 | 10.855 |
| 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.697 | 2.296 | 4.982 | 10.855 |
| 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: 8.392333984375e-5
1 2.994747e+03 1.083462e+03
* time: 0.0056498050689697266
2 2.406265e+03 1.884408e+02
* time: 0.01116180419921875
3 2.344175e+03 7.741610e+01
* time: 0.01638197898864746
4 2.323153e+03 2.907642e+01
* time: 0.021620988845825195
5 2.318222e+03 2.273295e+01
* time: 0.026874780654907227
6 2.316833e+03 1.390527e+01
* time: 0.032158851623535156
7 2.316425e+03 4.490883e+00
* time: 0.0375058650970459
8 2.316362e+03 9.374519e-01
* time: 0.042964935302734375
9 2.316356e+03 1.928785e-01
* time: 0.04831385612487793
10 2.316355e+03 3.119615e-02
* time: 0.05378580093383789
11 2.316355e+03 6.215513e-03
* time: 0.059282779693603516
12 2.316355e+03 8.313206e-04
* time: 0.06468796730041504FittedPumasModel
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.