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.022853851318359375
1 3.110315e+03 9.706222e+02
* time: 1.0462698936462402
2 2.831659e+03 7.817006e+02
* time: 1.1419320106506348
3 2.405281e+03 2.923696e+02
* time: 1.2141389846801758
4 2.370406e+03 3.032286e+02
* time: 1.2662639617919922
5 2.313631e+03 3.126188e+02
* time: 1.3207719326019287
6 2.263986e+03 2.946697e+02
* time: 1.4267120361328125
7 2.160182e+03 1.917599e+02
* time: 1.5041699409484863
8 2.112467e+03 1.588288e+02
* time: 1.5744109153747559
9 2.090339e+03 1.108334e+02
* time: 1.6208808422088623
10 2.078171e+03 8.108027e+01
* time: 1.667639970779419
11 2.074517e+03 7.813928e+01
* time: 1.708946943283081
12 2.066270e+03 7.044632e+01
* time: 1.7506470680236816
13 2.049660e+03 1.062584e+02
* time: 1.7920329570770264
14 2.021965e+03 1.130570e+02
* time: 1.8333568572998047
15 1.994936e+03 7.825801e+01
* time: 1.9092810153961182
16 1.979337e+03 5.318263e+01
* time: 1.944929838180542
17 1.972141e+03 6.807046e+01
* time: 1.9796228408813477
18 1.967973e+03 7.896361e+01
* time: 2.014770030975342
19 1.962237e+03 8.343757e+01
* time: 2.051054000854492
20 1.952791e+03 5.565304e+01
* time: 2.086704969406128
21 1.935857e+03 3.923284e+01
* time: 2.122864007949829
22 1.926254e+03 5.749643e+01
* time: 2.1586759090423584
23 1.922144e+03 4.306225e+01
* time: 2.1918258666992188
24 1.911566e+03 4.810496e+01
* time: 2.2284798622131348
25 1.906464e+03 4.324267e+01
* time: 2.2872068881988525
26 1.905339e+03 1.207954e+01
* time: 2.3207149505615234
27 1.905092e+03 1.093479e+01
* time: 2.3525798320770264
28 1.904957e+03 1.057034e+01
* time: 2.3846189975738525
29 1.904875e+03 1.060882e+01
* time: 2.4150240421295166
30 1.904459e+03 1.031525e+01
* time: 2.4471049308776855
31 1.903886e+03 1.041179e+01
* time: 2.479979991912842
32 1.903313e+03 1.135672e+01
* time: 2.5128300189971924
33 1.903057e+03 1.075683e+01
* time: 2.544301986694336
34 1.902950e+03 1.091295e+01
* time: 2.5775818824768066
35 1.902887e+03 1.042409e+01
* time: 2.6091790199279785
36 1.902640e+03 9.213043e+00
* time: 2.6653878688812256
37 1.902364e+03 9.519321e+00
* time: 2.6986310482025146
38 1.902156e+03 5.590984e+00
* time: 2.73110294342041
39 1.902094e+03 1.811898e+00
* time: 2.7624669075012207
40 1.902086e+03 1.644722e+00
* time: 2.7929978370666504
41 1.902084e+03 1.653520e+00
* time: 2.822460889816284
42 1.902074e+03 1.720184e+00
* time: 2.8532350063323975
43 1.902055e+03 2.619061e+00
* time: 2.883758068084717
44 1.902015e+03 3.519729e+00
* time: 2.9161360263824463
45 1.901962e+03 3.403372e+00
* time: 2.947866916656494
46 1.901924e+03 1.945644e+00
* time: 2.979714870452881
47 1.901914e+03 6.273342e-01
* time: 3.0127298831939697
48 1.901913e+03 5.374557e-01
* time: 3.0613958835601807
49 1.901913e+03 5.710007e-01
* time: 3.0928280353546143
50 1.901913e+03 6.091390e-01
* time: 3.124202013015747
51 1.901912e+03 6.692417e-01
* time: 3.155704975128174
52 1.901909e+03 7.579620e-01
* time: 3.187429904937744
53 1.901903e+03 8.798211e-01
* time: 3.219697952270508
54 1.901889e+03 1.002981e+00
* time: 3.252347946166992
55 1.901864e+03 1.495512e+00
* time: 3.2847700119018555
56 1.901837e+03 1.664621e+00
* time: 3.316210985183716
57 1.901819e+03 8.601118e-01
* time: 3.348296880722046
58 1.901815e+03 4.525470e-02
* time: 3.380599021911621
59 1.901815e+03 1.294280e-02
* time: 3.432102918624878
60 1.901815e+03 2.876568e-03
* time: 3.461632013320923
61 1.901815e+03 2.876568e-03
* time: 3.54441499710083
62 1.901815e+03 2.876568e-03
* time: 3.628018856048584FittedPumasModel
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.008148193359375e-5
1 2.846050e+03 1.128657e+03
* time: 0.08632898330688477
2 2.472982e+03 7.008264e+02
* time: 0.14580106735229492
3 2.180690e+03 2.312704e+02
* time: 0.24759197235107422
4 2.125801e+03 1.862929e+02
* time: 0.2957630157470703
5 2.103173e+03 1.320946e+02
* time: 0.3419220447540283
6 2.091136e+03 1.103035e+02
* time: 0.38486695289611816
7 2.081443e+03 1.091133e+02
* time: 0.4247100353240967
8 2.071793e+03 7.197675e+01
* time: 0.46649694442749023
9 2.062706e+03 7.623310e+01
* time: 0.5104320049285889
10 2.057515e+03 6.885476e+01
* time: 0.5546009540557861
11 2.051133e+03 6.368504e+01
* time: 0.5950629711151123
12 2.038626e+03 7.730243e+01
* time: 0.6355991363525391
13 2.019352e+03 1.136864e+02
* time: 0.7162010669708252
14 1.997136e+03 1.005748e+02
* time: 0.754997968673706
15 1.983023e+03 6.831478e+01
* time: 0.7964339256286621
16 1.977700e+03 5.649783e+01
* time: 0.8347539901733398
17 1.974583e+03 6.357642e+01
* time: 0.8718240261077881
18 1.967292e+03 7.658974e+01
* time: 0.9117641448974609
19 1.951045e+03 6.130573e+01
* time: 0.9575850963592529
20 1.935868e+03 4.845839e+01
* time: 0.99788498878479
21 1.929356e+03 6.325782e+01
* time: 1.0387070178985596
22 1.925187e+03 3.142245e+01
* time: 1.0762081146240234
23 1.923733e+03 4.623400e+01
* time: 1.1147370338439941
24 1.918498e+03 5.347738e+01
* time: 1.1704990863800049
25 1.912383e+03 5.849125e+01
* time: 1.2135579586029053
26 1.905510e+03 3.254038e+01
* time: 1.2541630268096924
27 1.903629e+03 2.905618e+01
* time: 1.2914140224456787
28 1.902833e+03 2.907696e+01
* time: 1.3274199962615967
29 1.902447e+03 2.746037e+01
* time: 1.3603029251098633
30 1.899399e+03 1.930949e+01
* time: 1.3964619636535645
31 1.898705e+03 1.186361e+01
* time: 1.431649923324585
32 1.898505e+03 1.050402e+01
* time: 1.4660649299621582
33 1.898474e+03 1.042186e+01
* time: 1.4965879917144775
34 1.897992e+03 1.238729e+01
* time: 1.5510480403900146
35 1.897498e+03 1.729368e+01
* time: 1.5843441486358643
36 1.896954e+03 1.472554e+01
* time: 1.617461919784546
37 1.896744e+03 5.852825e+00
* time: 1.649656057357788
38 1.896705e+03 1.171353e+00
* time: 1.679095983505249
39 1.896704e+03 1.216117e+00
* time: 1.709317922592163
40 1.896703e+03 1.230336e+00
* time: 1.7389841079711914
41 1.896698e+03 1.250715e+00
* time: 1.7696609497070312
42 1.896688e+03 1.449552e+00
* time: 1.8012731075286865
43 1.896666e+03 2.533300e+00
* time: 1.8319931030273438
44 1.896631e+03 3.075536e+00
* time: 1.8627030849456787
45 1.896599e+03 2.139797e+00
* time: 1.8956079483032227
46 1.896587e+03 6.636031e-01
* time: 1.9433300495147705
47 1.896585e+03 6.303517e-01
* time: 1.97493314743042
48 1.896585e+03 5.995265e-01
* time: 2.004667043685913
49 1.896584e+03 5.844159e-01
* time: 2.0348191261291504
50 1.896583e+03 6.083858e-01
* time: 2.065570116043091
51 1.896579e+03 8.145326e-01
* time: 2.096142053604126
52 1.896570e+03 1.293490e+00
* time: 2.1271090507507324
53 1.896549e+03 1.877705e+00
* time: 2.15854811668396
54 1.896513e+03 2.217391e+00
* time: 2.1897449493408203
55 1.896482e+03 1.658147e+00
* time: 2.2208950519561768
56 1.896466e+03 5.207215e-01
* time: 2.252514123916626
57 1.896463e+03 1.177468e-01
* time: 2.3054800033569336
58 1.896463e+03 1.678937e-02
* time: 2.3351259231567383
59 1.896463e+03 2.666635e-03
* time: 2.3621439933776855
60 1.896463e+03 2.666635e-03
* time: 2.4204909801483154FittedPumasModel
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: 6.794929504394531e-5
1 3.572050e+03 1.302046e+03
* time: 0.17165684700012207
2 3.266947e+03 5.384244e+02
* time: 0.23880887031555176
3 3.150635e+03 1.918079e+02
* time: 0.297792911529541
4 3.108122e+03 1.277799e+02
* time: 0.3911118507385254
5 3.091358e+03 8.838080e+01
* time: 0.4435398578643799
6 3.082997e+03 8.321689e+01
* time: 0.49619603157043457
7 3.076379e+03 8.167702e+01
* time: 0.548677921295166
8 3.067422e+03 1.177822e+02
* time: 0.6121079921722412
9 3.048580e+03 2.526969e+02
* time: 0.6740949153900146
10 2.993096e+03 6.325396e+02
* time: 0.7691209316253662
11 2.965744e+03 7.634718e+02
* time: 0.9668519496917725
12 2.921212e+03 9.704020e+02
* time: 1.0756778717041016
13 2.553649e+03 6.642510e+02
* time: 1.1699910163879395
14 2.319495e+03 3.204711e+02
* time: 1.302757978439331
15 2.183040e+03 2.174905e+02
* time: 1.4131379127502441
16 2.157621e+03 3.150983e+02
* time: 1.4680969715118408
17 2.132395e+03 2.847614e+02
* time: 1.5204799175262451
18 2.084799e+03 1.563370e+02
* time: 1.571202039718628
19 2.071497e+03 1.006429e+02
* time: 1.6296749114990234
20 2.064983e+03 9.753313e+01
* time: 1.6803638935089111
21 2.048289e+03 9.230405e+01
* time: 1.7326679229736328
22 2.020938e+03 1.292359e+02
* time: 1.7836878299713135
23 1.983888e+03 2.284990e+02
* time: 1.844878911972046
24 1.962132e+03 1.220188e+02
* time: 1.8890049457550049
25 1.945947e+03 1.035894e+02
* time: 1.9319298267364502
26 1.917782e+03 8.246698e+01
* time: 1.9774789810180664
27 1.905967e+03 5.408054e+01
* time: 2.031255006790161
28 1.898569e+03 2.172222e+01
* time: 2.0751659870147705
29 1.897473e+03 1.689350e+01
* time: 2.116830825805664
30 1.897019e+03 1.666689e+01
* time: 2.157611846923828
31 1.896796e+03 1.699751e+01
* time: 2.198110818862915
32 1.896559e+03 1.645900e+01
* time: 2.2505300045013428
33 1.896223e+03 1.415504e+01
* time: 2.292948007583618
34 1.895773e+03 1.630081e+01
* time: 2.3363430500030518
35 1.895309e+03 1.723930e+01
* time: 2.3798139095306396
36 1.895004e+03 1.229983e+01
* time: 2.431720018386841
37 1.894871e+03 5.385102e+00
* time: 2.4741618633270264
38 1.894827e+03 3.465463e+00
* time: 2.5163609981536865
39 1.894816e+03 3.387474e+00
* time: 2.556527853012085
40 1.894807e+03 3.295388e+00
* time: 2.596411943435669
41 1.894786e+03 3.089194e+00
* time: 2.6451399326324463
42 1.894737e+03 2.928080e+00
* time: 2.6863210201263428
43 1.894624e+03 3.088723e+00
* time: 2.7281298637390137
44 1.894413e+03 3.493791e+00
* time: 2.7698628902435303
45 1.894127e+03 3.142865e+00
* time: 2.8195278644561768
46 1.893933e+03 2.145253e+00
* time: 2.862872838973999
47 1.893888e+03 2.172800e+00
* time: 2.9045238494873047
48 1.893880e+03 2.180509e+00
* time: 2.94720196723938
49 1.893876e+03 2.134257e+00
* time: 2.987074851989746
50 1.893868e+03 2.032354e+00
* time: 3.0392038822174072
51 1.893846e+03 1.760874e+00
* time: 3.0790200233459473
52 1.893796e+03 1.779016e+00
* time: 3.119339942932129
53 1.893694e+03 2.018956e+00
* time: 3.1597280502319336
54 1.893559e+03 2.366854e+00
* time: 3.207430839538574
55 1.893474e+03 3.690142e+00
* time: 3.2476818561553955
56 1.893446e+03 3.675109e+00
* time: 3.286834955215454
57 1.893439e+03 3.426419e+00
* time: 3.326249837875366
58 1.893429e+03 3.183164e+00
* time: 3.3639039993286133
59 1.893398e+03 2.695171e+00
* time: 3.410275936126709
60 1.893328e+03 2.753548e+00
* time: 3.4492859840393066
61 1.893169e+03 3.589748e+00
* time: 3.4890220165252686
62 1.892920e+03 3.680718e+00
* time: 3.5291659832000732
63 1.892667e+03 2.568107e+00
* time: 3.576340913772583
64 1.892514e+03 1.087910e+00
* time: 3.6155779361724854
65 1.892493e+03 3.287296e-01
* time: 3.6542069911956787
66 1.892492e+03 2.967465e-01
* time: 3.6918649673461914
67 1.892492e+03 3.020682e-01
* time: 3.735811948776245
68 1.892491e+03 3.034704e-01
* time: 3.770577907562256
69 1.892491e+03 3.091846e-01
* time: 3.80651593208313
70 1.892491e+03 3.224170e-01
* time: 3.842834949493408
71 1.892490e+03 6.494197e-01
* time: 3.878938913345337
72 1.892488e+03 1.115188e+00
* time: 3.92323899269104
73 1.892483e+03 1.838833e+00
* time: 3.9602699279785156
74 1.892472e+03 2.765371e+00
* time: 3.9979000091552734
75 1.892452e+03 3.463807e+00
* time: 4.03542685508728
76 1.892431e+03 2.805270e+00
* time: 4.072849988937378
77 1.892411e+03 5.758916e-01
* time: 4.119343996047974
78 1.892410e+03 1.434041e-01
* time: 4.1562819480896
79 1.892409e+03 1.639246e-01
* time: 4.1926429271698
80 1.892409e+03 1.145856e-01
* time: 4.2280988693237305
81 1.892409e+03 3.966861e-02
* time: 4.270687818527222
82 1.892409e+03 3.550808e-02
* time: 4.304358005523682
83 1.892409e+03 3.456241e-02
* time: 4.337874889373779
84 1.892409e+03 3.114018e-02
* time: 4.370725870132446
85 1.892409e+03 4.080806e-02
* time: 4.403383016586304
86 1.892409e+03 6.722726e-02
* time: 4.444466829299927
87 1.892409e+03 1.006791e-01
* time: 4.477787971496582
88 1.892409e+03 1.303988e-01
* time: 4.511816024780273
89 1.892409e+03 1.228919e-01
* time: 4.5449018478393555
90 1.892409e+03 6.433813e-02
* time: 4.578007936477661
91 1.892409e+03 1.314164e-02
* time: 4.619120836257935
92 1.892409e+03 4.929931e-04
* time: 4.651192903518677FittedPumasModel
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: 5.698204040527344e-5
1 3.641387e+03 1.432080e+03
* time: 0.08399701118469238
2 3.290450e+03 5.274782e+02
* time: 0.1405479907989502
3 3.185512e+03 2.173676e+02
* time: 0.1915149688720703
4 3.143009e+03 1.479653e+02
* time: 0.27722787857055664
5 3.128511e+03 8.980031e+01
* time: 0.3225889205932617
6 3.123188e+03 5.033125e+01
* time: 0.36745786666870117
7 3.120794e+03 4.279722e+01
* time: 0.41774582862854004
8 3.118627e+03 3.971051e+01
* time: 0.4665567874908447
9 3.115300e+03 8.456587e+01
* time: 0.5372278690338135
10 3.109353e+03 1.350354e+02
* time: 0.5835208892822266
11 3.095894e+03 1.998258e+02
* time: 0.632476806640625
12 2.988214e+03 4.366433e+02
* time: 0.7082998752593994
13 2.896081e+03 5.505943e+02
* time: 0.906702995300293
14 2.652467e+03 7.300323e+02
* time: 2.092374801635742
15 2.560937e+03 6.973661e+02
* time: 2.3706278800964355
16 2.254941e+03 2.740033e+02
* time: 2.428736925125122
17 2.222509e+03 2.034303e+02
* time: 2.478142023086548
18 2.171255e+03 2.449580e+02
* time: 2.527132987976074
19 2.024532e+03 1.121511e+02
* time: 2.588225841522217
20 1.993723e+03 1.042814e+02
* time: 2.6339569091796875
21 1.985113e+03 8.079014e+01
* time: 2.679111957550049
22 1.976757e+03 7.054196e+01
* time: 2.7230989933013916
23 1.969970e+03 6.070322e+01
* time: 2.7658679485321045
24 1.961095e+03 6.810782e+01
* time: 2.808180809020996
25 1.947983e+03 8.116920e+01
* time: 2.8641319274902344
26 1.930371e+03 8.530051e+01
* time: 2.9077839851379395
27 1.910209e+03 6.993170e+01
* time: 2.951231002807617
28 1.899107e+03 3.362640e+01
* time: 2.994947910308838
29 1.898022e+03 2.642220e+01
* time: 3.035061836242676
30 1.897055e+03 1.213144e+01
* time: 3.0899558067321777
31 1.896596e+03 7.773239e+00
* time: 3.1325888633728027
32 1.896538e+03 7.997039e+00
* time: 3.1738059520721436
33 1.896451e+03 8.160909e+00
* time: 3.212535858154297
34 1.896283e+03 8.237721e+00
* time: 3.2520408630371094
35 1.895903e+03 1.520219e+01
* time: 3.2920379638671875
36 1.895272e+03 2.358916e+01
* time: 3.3471109867095947
37 1.894536e+03 2.461296e+01
* time: 3.389657974243164
38 1.893995e+03 1.546128e+01
* time: 3.431159019470215
39 1.893858e+03 6.976137e+00
* time: 3.471750020980835
40 1.893833e+03 6.019466e+00
* time: 3.5104639530181885
41 1.893786e+03 3.827201e+00
* time: 3.5622618198394775
42 1.893714e+03 3.323412e+00
* time: 3.6038658618927
43 1.893592e+03 3.215150e+00
* time: 3.645493984222412
44 1.893435e+03 6.534965e+00
* time: 3.687417984008789
45 1.893286e+03 7.424154e+00
* time: 3.7291979789733887
46 1.893190e+03 5.552627e+00
* time: 3.769570827484131
47 1.893139e+03 3.222316e+00
* time: 3.823171854019165
48 1.893120e+03 3.015339e+00
* time: 3.863693952560425
49 1.893107e+03 3.244809e+00
* time: 3.9023988246917725
50 1.893080e+03 6.163100e+00
* time: 3.940603017807007
51 1.893027e+03 9.824713e+00
* time: 3.979189872741699
52 1.892912e+03 1.390100e+01
* time: 4.018296003341675
53 1.892734e+03 1.510937e+01
* time: 4.0729079246521
54 1.892561e+03 1.008563e+01
* time: 4.112792015075684
55 1.892485e+03 3.730668e+00
* time: 4.152537822723389
56 1.892471e+03 3.380261e+00
* time: 4.191667795181274
57 1.892463e+03 3.167904e+00
* time: 4.229518890380859
58 1.892441e+03 4.152065e+00
* time: 4.279533863067627
59 1.892391e+03 7.355996e+00
* time: 4.319795846939087
60 1.892268e+03 1.195397e+01
* time: 4.359622001647949
61 1.892026e+03 1.640783e+01
* time: 4.399518013000488
62 1.891735e+03 1.593576e+01
* time: 4.439692974090576
63 1.891569e+03 8.316423e+00
* time: 4.480232000350952
64 1.891494e+03 3.948212e+00
* time: 4.532528877258301
65 1.891481e+03 3.911593e+00
* time: 4.571582794189453
66 1.891457e+03 3.875559e+00
* time: 4.610850811004639
67 1.891405e+03 3.811247e+00
* time: 4.650377988815308
68 1.891262e+03 3.657045e+00
* time: 4.690889835357666
69 1.890930e+03 4.957405e+00
* time: 4.730949878692627
70 1.890317e+03 6.657726e+00
* time: 4.78519082069397
71 1.889660e+03 6.086302e+00
* time: 4.827096939086914
72 1.889303e+03 2.270929e+00
* time: 4.867990970611572
73 1.889253e+03 7.695301e-01
* time: 4.907615900039673
74 1.889252e+03 7.382144e-01
* time: 4.945988893508911
75 1.889251e+03 7.187898e-01
* time: 4.997506856918335
76 1.889251e+03 7.215047e-01
* time: 5.035559892654419
77 1.889250e+03 7.235155e-01
* time: 5.073261022567749
78 1.889249e+03 7.246818e-01
* time: 5.110736846923828
79 1.889244e+03 7.257796e-01
* time: 5.1492719650268555
80 1.889233e+03 7.198190e-01
* time: 5.188949823379517
81 1.889204e+03 1.089029e+00
* time: 5.240418910980225
82 1.889142e+03 1.801601e+00
* time: 5.282528877258301
83 1.889043e+03 2.967917e+00
* time: 5.324095010757446
84 1.888889e+03 2.965856e+00
* time: 5.364868879318237
85 1.888705e+03 5.933557e-01
* time: 5.406105995178223
86 1.888655e+03 9.577696e-01
* time: 5.446254014968872
87 1.888582e+03 1.498494e+00
* time: 5.50121283531189
88 1.888533e+03 1.502753e+00
* time: 5.54264497756958
89 1.888490e+03 1.184664e+00
* time: 5.583887815475464
90 1.888480e+03 6.684517e-01
* time: 5.624224901199341
91 1.888476e+03 3.680034e-01
* time: 5.66509485244751
92 1.888476e+03 4.720040e-01
* time: 5.716404914855957
93 1.888476e+03 4.768646e-01
* time: 5.756554841995239
94 1.888475e+03 4.736673e-01
* time: 5.795373916625977
95 1.888475e+03 4.552764e-01
* time: 5.834204912185669
96 1.888474e+03 5.193733e-01
* time: 5.872846841812134
97 1.888473e+03 8.850112e-01
* time: 5.91143798828125
98 1.888468e+03 1.461600e+00
* time: 5.962602853775024
99 1.888458e+03 2.209127e+00
* time: 6.0014328956604
100 1.888437e+03 2.961239e+00
* time: 6.040201902389526
101 1.888407e+03 2.978463e+00
* time: 6.0788328647613525
102 1.888384e+03 1.707201e+00
* time: 6.118496894836426
103 1.888381e+03 6.199354e-01
* time: 6.159245014190674
104 1.888380e+03 5.170909e-01
* time: 6.213109970092773
105 1.888378e+03 1.037408e-01
* time: 6.2544519901275635
106 1.888378e+03 8.473247e-02
* time: 6.292920827865601
107 1.888378e+03 8.364978e-02
* time: 6.331905841827393
108 1.888378e+03 8.080446e-02
* time: 6.370759963989258
109 1.888378e+03 7.873905e-02
* time: 6.407993793487549
110 1.888378e+03 7.798398e-02
* time: 6.458576917648315
111 1.888378e+03 7.788170e-02
* time: 6.4956488609313965
112 1.888378e+03 7.776461e-02
* time: 6.532890796661377
113 1.888378e+03 9.023662e-02
* time: 6.570108890533447
114 1.888378e+03 1.631390e-01
* time: 6.607915878295898
115 1.888378e+03 2.768731e-01
* time: 6.644001007080078
116 1.888377e+03 4.462386e-01
* time: 6.694090843200684
117 1.888377e+03 6.643297e-01
* time: 6.732261896133423
118 1.888375e+03 8.433397e-01
* time: 6.7696449756622314
119 1.888374e+03 7.596814e-01
* time: 6.807506799697876
120 1.888373e+03 3.638119e-01
* time: 6.844988822937012
121 1.888372e+03 8.306034e-02
* time: 6.881482839584351
122 1.888372e+03 2.084513e-02
* time: 6.9306960105896
123 1.888372e+03 2.056419e-02
* time: 6.966517925262451
124 1.888372e+03 2.044080e-02
* time: 7.000717878341675
125 1.888372e+03 2.035196e-02
* time: 7.0351879596710205
126 1.888372e+03 2.021264e-02
* time: 7.071230888366699
127 1.888372e+03 1.998163e-02
* time: 7.107580900192261
128 1.888372e+03 3.161638e-02
* time: 7.1580469608306885
129 1.888372e+03 5.509218e-02
* time: 7.198298931121826
130 1.888372e+03 9.275848e-02
* time: 7.24151086807251
131 1.888372e+03 1.528749e-01
* time: 7.281561851501465
132 1.888372e+03 2.461766e-01
* time: 7.321044921875
133 1.888372e+03 3.799362e-01
* time: 7.361806869506836
134 1.888371e+03 5.311665e-01
* time: 7.418967962265015
135 1.888369e+03 6.019039e-01
* time: 7.461224794387817
136 1.888367e+03 4.664896e-01
* time: 7.502794981002808
137 1.888366e+03 1.404934e-01
* time: 7.54407000541687
138 1.888365e+03 8.513331e-02
* time: 7.584707975387573
139 1.888364e+03 1.244077e-01
* time: 7.624077796936035
140 1.888364e+03 1.028085e-01
* time: 7.677246809005737
141 1.888364e+03 5.162231e-02
* time: 7.7190399169921875
142 1.888364e+03 5.149630e-02
* time: 7.759603977203369
143 1.888364e+03 3.147284e-02
* time: 7.799970865249634
144 1.888364e+03 2.104766e-02
* time: 7.837251901626587
145 1.888364e+03 6.539151e-03
* time: 7.872690916061401
146 1.888364e+03 2.540196e-03
* time: 7.924468994140625
147 1.888364e+03 4.362661e-03
* time: 7.962589979171753
148 1.888364e+03 3.034416e-03
* time: 7.9999167919158936
149 1.888364e+03 5.953892e-04
* time: 8.036953926086426FittedPumasModel
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.628 |
| 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.421 |
| 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.651 |
| 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.037 |
| 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.628 | 2.421 | 4.651 | 8.037 |
| 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: 5.984306335449219e-5
1 2.994747e+03 1.083462e+03
* time: 0.005609989166259766
2 2.406265e+03 1.884408e+02
* time: 0.01100301742553711
3 2.344175e+03 7.741610e+01
* time: 0.016587018966674805
4 2.323153e+03 2.907642e+01
* time: 0.02209782600402832
5 2.318222e+03 2.273295e+01
* time: 0.027559995651245117
6 2.316833e+03 1.390527e+01
* time: 0.033174991607666016
7 2.316425e+03 4.490883e+00
* time: 0.03866982460021973
8 2.316362e+03 9.374519e-01
* time: 0.04416084289550781
9 2.316356e+03 1.928785e-01
* time: 0.04961395263671875
10 2.316355e+03 3.119615e-02
* time: 0.05507183074951172
11 2.316355e+03 6.215513e-03
* time: 0.06056499481201172
12 2.316355e+03 8.313206e-04
* time: 0.06600403785705566FittedPumasModel
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.