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.037613868713378906
1 3.110315e+03 9.706222e+02
* time: 1.6927227973937988
2 2.831659e+03 7.817006e+02
* time: 1.8174989223480225
3 2.405281e+03 2.923696e+02
* time: 1.9134018421173096
4 2.370406e+03 3.032286e+02
* time: 1.9823029041290283
5 2.313631e+03 3.126188e+02
* time: 2.053475856781006
6 2.263986e+03 2.946697e+02
* time: 2.1186938285827637
7 2.160182e+03 1.917599e+02
* time: 2.4834418296813965
8 2.112467e+03 1.588288e+02
* time: 2.5671849250793457
9 2.090339e+03 1.108334e+02
* time: 2.6232829093933105
10 2.078171e+03 8.108027e+01
* time: 2.6788039207458496
11 2.074517e+03 7.813928e+01
* time: 2.7246339321136475
12 2.066270e+03 7.044632e+01
* time: 2.7693469524383545
13 2.049660e+03 1.062584e+02
* time: 2.8156960010528564
14 2.021965e+03 1.130570e+02
* time: 2.8640048503875732
15 1.994936e+03 7.825801e+01
* time: 2.923285961151123
16 1.979337e+03 5.318263e+01
* time: 2.9819259643554688
17 1.972141e+03 6.807046e+01
* time: 3.041278839111328
18 1.967973e+03 7.896361e+01
* time: 3.099066972732544
19 1.962237e+03 8.343757e+01
* time: 3.301644802093506
20 1.952791e+03 5.565304e+01
* time: 3.365363836288452
21 1.935857e+03 3.923284e+01
* time: 3.4615519046783447
22 1.926254e+03 5.749643e+01
* time: 3.5150527954101562
23 1.922144e+03 4.306225e+01
* time: 3.5657858848571777
24 1.911566e+03 4.810496e+01
* time: 3.618868827819824
25 1.906464e+03 4.324267e+01
* time: 3.672146797180176
26 1.905339e+03 1.207954e+01
* time: 3.7202658653259277
27 1.905092e+03 1.093479e+01
* time: 3.767251968383789
28 1.904957e+03 1.057034e+01
* time: 3.8143060207366943
29 1.904875e+03 1.060882e+01
* time: 3.860481023788452
30 1.904459e+03 1.031525e+01
* time: 3.9099488258361816
31 1.903886e+03 1.041179e+01
* time: 3.9600419998168945
32 1.903313e+03 1.135672e+01
* time: 4.080448865890503
33 1.903057e+03 1.075683e+01
* time: 4.124199867248535
34 1.902950e+03 1.091295e+01
* time: 4.1670379638671875
35 1.902887e+03 1.042409e+01
* time: 4.208613872528076
36 1.902640e+03 9.213043e+00
* time: 4.253049850463867
37 1.902364e+03 9.519321e+00
* time: 4.296263933181763
38 1.902156e+03 5.590984e+00
* time: 4.339586973190308
39 1.902094e+03 1.811898e+00
* time: 4.381480932235718
40 1.902086e+03 1.644722e+00
* time: 4.423092842102051
41 1.902084e+03 1.653520e+00
* time: 4.462827920913696
42 1.902074e+03 1.720184e+00
* time: 4.505136966705322
43 1.902055e+03 2.619061e+00
* time: 4.54554295539856
44 1.902015e+03 3.519729e+00
* time: 4.586732864379883
45 1.901962e+03 3.403372e+00
* time: 4.631021022796631
46 1.901924e+03 1.945644e+00
* time: 4.698637008666992
47 1.901914e+03 6.273342e-01
* time: 4.740691900253296
48 1.901913e+03 5.374557e-01
* time: 4.779408931732178
49 1.901913e+03 5.710007e-01
* time: 4.816184997558594
50 1.901913e+03 6.091390e-01
* time: 4.853482961654663
51 1.901912e+03 6.692417e-01
* time: 4.891979932785034
52 1.901909e+03 7.579620e-01
* time: 4.931056022644043
53 1.901903e+03 8.798211e-01
* time: 4.9715728759765625
54 1.901889e+03 1.002981e+00
* time: 5.010846853256226
55 1.901864e+03 1.495512e+00
* time: 5.050752878189087
56 1.901837e+03 1.664621e+00
* time: 5.088692903518677
57 1.901819e+03 8.601118e-01
* time: 5.127067804336548
58 1.901815e+03 4.525470e-02
* time: 5.166136980056763
59 1.901815e+03 1.294280e-02
* time: 5.209069013595581
60 1.901815e+03 2.876568e-03
* time: 5.2988879680633545
61 1.901815e+03 2.876568e-03
* time: 5.399693012237549
62 1.901815e+03 2.876568e-03
* time: 5.503976821899414FittedPumasModel
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: 0.00010704994201660156
1 2.846050e+03 1.128657e+03
* time: 0.0870969295501709
2 2.472982e+03 7.008264e+02
* time: 0.14429807662963867
3 2.180690e+03 2.312704e+02
* time: 0.20038104057312012
4 2.125801e+03 1.862929e+02
* time: 0.24765992164611816
5 2.103173e+03 1.320946e+02
* time: 0.29416799545288086
6 2.091136e+03 1.103035e+02
* time: 0.3423731327056885
7 2.081443e+03 1.091133e+02
* time: 0.387653112411499
8 2.071793e+03 7.197675e+01
* time: 0.43422913551330566
9 2.062706e+03 7.623310e+01
* time: 0.4820399284362793
10 2.057515e+03 6.885476e+01
* time: 0.5353949069976807
11 2.051133e+03 6.368504e+01
* time: 0.6937949657440186
12 2.038626e+03 7.730243e+01
* time: 0.7380599975585938
13 2.019352e+03 1.136864e+02
* time: 0.7822999954223633
14 1.997136e+03 1.005748e+02
* time: 0.8268520832061768
15 1.983023e+03 6.831478e+01
* time: 0.8729031085968018
16 1.977700e+03 5.649783e+01
* time: 0.9174380302429199
17 1.974583e+03 6.357642e+01
* time: 0.9641721248626709
18 1.967292e+03 7.658974e+01
* time: 1.0170819759368896
19 1.951045e+03 6.130573e+01
* time: 1.0707659721374512
20 1.935868e+03 4.845839e+01
* time: 1.1178460121154785
21 1.929356e+03 6.325782e+01
* time: 1.1662590503692627
22 1.925187e+03 3.142245e+01
* time: 1.211864948272705
23 1.923733e+03 4.623400e+01
* time: 1.2562329769134521
24 1.918498e+03 5.347738e+01
* time: 1.3574779033660889
25 1.912383e+03 5.849125e+01
* time: 1.4085490703582764
26 1.905510e+03 3.254038e+01
* time: 1.4632880687713623
27 1.903629e+03 2.905618e+01
* time: 1.508781909942627
28 1.902833e+03 2.907696e+01
* time: 1.5532031059265137
29 1.902447e+03 2.746037e+01
* time: 1.593545913696289
30 1.899399e+03 1.930949e+01
* time: 1.6377909183502197
31 1.898705e+03 1.186361e+01
* time: 1.6822879314422607
32 1.898505e+03 1.050402e+01
* time: 1.7258970737457275
33 1.898474e+03 1.042186e+01
* time: 1.7653040885925293
34 1.897992e+03 1.238729e+01
* time: 1.8097329139709473
35 1.897498e+03 1.729368e+01
* time: 1.8515160083770752
36 1.896954e+03 1.472554e+01
* time: 1.8983159065246582
37 1.896744e+03 5.852825e+00
* time: 1.9947080612182617
38 1.896705e+03 1.171353e+00
* time: 2.029391050338745
39 1.896704e+03 1.216117e+00
* time: 2.065200090408325
40 1.896703e+03 1.230336e+00
* time: 2.100904941558838
41 1.896698e+03 1.250715e+00
* time: 2.1370201110839844
42 1.896688e+03 1.449552e+00
* time: 2.1727540493011475
43 1.896666e+03 2.533300e+00
* time: 2.2095561027526855
44 1.896631e+03 3.075536e+00
* time: 2.2464699745178223
45 1.896599e+03 2.139797e+00
* time: 2.2826199531555176
46 1.896587e+03 6.636031e-01
* time: 2.319087028503418
47 1.896585e+03 6.303517e-01
* time: 2.3556110858917236
48 1.896585e+03 5.995265e-01
* time: 2.3899149894714355
49 1.896584e+03 5.844159e-01
* time: 2.423732042312622
50 1.896583e+03 6.083858e-01
* time: 2.4617679119110107
51 1.896579e+03 8.145326e-01
* time: 2.529026985168457
52 1.896570e+03 1.293490e+00
* time: 2.5647029876708984
53 1.896549e+03 1.877705e+00
* time: 2.6010890007019043
54 1.896513e+03 2.217391e+00
* time: 2.6366829872131348
55 1.896482e+03 1.658147e+00
* time: 2.6729581356048584
56 1.896466e+03 5.207215e-01
* time: 2.7091219425201416
57 1.896463e+03 1.177468e-01
* time: 2.7436161041259766
58 1.896463e+03 1.678937e-02
* time: 2.7815821170806885
59 1.896463e+03 2.666635e-03
* time: 2.8124070167541504
60 1.896463e+03 2.666635e-03
* time: 2.8821511268615723FittedPumasModel
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: 8.296966552734375e-5
1 3.572050e+03 1.302046e+03
* time: 0.1417100429534912
2 3.266947e+03 5.384244e+02
* time: 0.22774195671081543
3 3.150635e+03 1.918079e+02
* time: 0.3013448715209961
4 3.108122e+03 1.277799e+02
* time: 0.37005186080932617
5 3.091358e+03 8.838080e+01
* time: 0.4475870132446289
6 3.082997e+03 8.321689e+01
* time: 0.6433870792388916
7 3.076379e+03 8.167702e+01
* time: 0.7018990516662598
8 3.067422e+03 1.177822e+02
* time: 0.7698838710784912
9 3.048580e+03 2.526969e+02
* time: 0.8349480628967285
10 2.993096e+03 6.325396e+02
* time: 0.937824010848999
11 2.965744e+03 7.634718e+02
* time: 1.16143798828125
12 2.921212e+03 9.704020e+02
* time: 1.2975618839263916
13 2.553649e+03 6.642510e+02
* time: 1.4014148712158203
14 2.319495e+03 3.204711e+02
* time: 1.628066062927246
15 2.183040e+03 2.174905e+02
* time: 1.7413899898529053
16 2.157621e+03 3.150983e+02
* time: 1.7996799945831299
17 2.132395e+03 2.847614e+02
* time: 1.8583929538726807
18 2.084799e+03 1.563370e+02
* time: 1.9126009941101074
19 2.071497e+03 1.006429e+02
* time: 1.9681079387664795
20 2.064983e+03 9.753313e+01
* time: 2.0285098552703857
21 2.048289e+03 9.230405e+01
* time: 2.097430944442749
22 2.020938e+03 1.292359e+02
* time: 2.1566848754882812
23 1.983888e+03 2.284990e+02
* time: 2.223214864730835
24 1.962132e+03 1.220188e+02
* time: 2.2848639488220215
25 1.945947e+03 1.035894e+02
* time: 2.3785290718078613
26 1.917782e+03 8.246698e+01
* time: 2.431752920150757
27 1.905967e+03 5.408054e+01
* time: 2.4941189289093018
28 1.898569e+03 2.172222e+01
* time: 2.5461950302124023
29 1.897473e+03 1.689350e+01
* time: 2.5955779552459717
30 1.897019e+03 1.666689e+01
* time: 2.6452879905700684
31 1.896796e+03 1.699751e+01
* time: 2.6950109004974365
32 1.896559e+03 1.645900e+01
* time: 2.749758005142212
33 1.896223e+03 1.415504e+01
* time: 2.8059210777282715
34 1.895773e+03 1.630081e+01
* time: 2.859092950820923
35 1.895309e+03 1.723930e+01
* time: 2.9123990535736084
36 1.895004e+03 1.229983e+01
* time: 2.9631779193878174
37 1.894871e+03 5.385102e+00
* time: 3.011909008026123
38 1.894827e+03 3.465463e+00
* time: 3.100877046585083
39 1.894816e+03 3.387474e+00
* time: 3.1442930698394775
40 1.894807e+03 3.295388e+00
* time: 3.1884028911590576
41 1.894786e+03 3.089194e+00
* time: 3.233250856399536
42 1.894737e+03 2.928080e+00
* time: 3.2783100605010986
43 1.894624e+03 3.088723e+00
* time: 3.3236958980560303
44 1.894413e+03 3.493791e+00
* time: 3.370123863220215
45 1.894127e+03 3.142865e+00
* time: 3.4217710494995117
46 1.893933e+03 2.145253e+00
* time: 3.472043037414551
47 1.893888e+03 2.172800e+00
* time: 3.522486925125122
48 1.893880e+03 2.180509e+00
* time: 3.5744850635528564
49 1.893876e+03 2.134257e+00
* time: 3.623337984085083
50 1.893868e+03 2.032354e+00
* time: 3.669792890548706
51 1.893846e+03 1.760874e+00
* time: 3.752027988433838
52 1.893796e+03 1.779016e+00
* time: 3.7953689098358154
53 1.893694e+03 2.018956e+00
* time: 3.8391449451446533
54 1.893559e+03 2.366854e+00
* time: 3.881906032562256
55 1.893474e+03 3.690142e+00
* time: 3.928046941757202
56 1.893446e+03 3.675109e+00
* time: 3.9780049324035645
57 1.893439e+03 3.426419e+00
* time: 4.02300500869751
58 1.893429e+03 3.183164e+00
* time: 4.073168992996216
59 1.893398e+03 2.695171e+00
* time: 4.125269889831543
60 1.893328e+03 2.753548e+00
* time: 4.181273937225342
61 1.893169e+03 3.589748e+00
* time: 4.235162973403931
62 1.892920e+03 3.680718e+00
* time: 4.289475917816162
63 1.892667e+03 2.568107e+00
* time: 4.344938039779663
64 1.892514e+03 1.087910e+00
* time: 4.452255964279175
65 1.892493e+03 3.287296e-01
* time: 4.497959852218628
66 1.892492e+03 2.967465e-01
* time: 4.542850971221924
67 1.892492e+03 3.020682e-01
* time: 4.586003065109253
68 1.892491e+03 3.034704e-01
* time: 4.6264259815216064
69 1.892491e+03 3.091846e-01
* time: 4.668262958526611
70 1.892491e+03 3.224170e-01
* time: 4.7104880809783936
71 1.892490e+03 6.494197e-01
* time: 4.755203008651733
72 1.892488e+03 1.115188e+00
* time: 4.802669048309326
73 1.892483e+03 1.838833e+00
* time: 4.849531888961792
74 1.892472e+03 2.765371e+00
* time: 4.896260976791382
75 1.892452e+03 3.463807e+00
* time: 4.944777011871338
76 1.892431e+03 2.805270e+00
* time: 4.99278998374939
77 1.892411e+03 5.758916e-01
* time: 5.057107925415039
78 1.892410e+03 1.434041e-01
* time: 5.155796051025391
79 1.892409e+03 1.639246e-01
* time: 5.211509943008423
80 1.892409e+03 1.145856e-01
* time: 5.260776996612549
81 1.892409e+03 3.966861e-02
* time: 5.305372953414917
82 1.892409e+03 3.550808e-02
* time: 5.350414037704468
83 1.892409e+03 3.456241e-02
* time: 5.398382902145386
84 1.892409e+03 3.114018e-02
* time: 5.440432071685791
85 1.892409e+03 4.080806e-02
* time: 5.486232042312622
86 1.892409e+03 6.722726e-02
* time: 5.5426318645477295
87 1.892409e+03 1.006791e-01
* time: 5.594470024108887
88 1.892409e+03 1.303988e-01
* time: 5.647830009460449
89 1.892409e+03 1.228919e-01
* time: 5.691498041152954
90 1.892409e+03 6.433813e-02
* time: 5.737854957580566
91 1.892409e+03 1.314164e-02
* time: 5.789722919464111
92 1.892409e+03 4.929931e-04
* time: 5.881663084030151FittedPumasModel
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: 9.679794311523438e-5
1 3.641387e+03 1.432080e+03
* time: 0.13288593292236328
2 3.290450e+03 5.274782e+02
* time: 0.2219858169555664
3 3.185512e+03 2.173676e+02
* time: 0.3049478530883789
4 3.143009e+03 1.479653e+02
* time: 0.3828868865966797
5 3.128511e+03 8.980031e+01
* time: 0.6217238903045654
6 3.123188e+03 5.033125e+01
* time: 0.6880419254302979
7 3.120794e+03 4.279722e+01
* time: 0.7538299560546875
8 3.118627e+03 3.971051e+01
* time: 0.8227658271789551
9 3.115300e+03 8.456587e+01
* time: 0.8885397911071777
10 3.109353e+03 1.350354e+02
* time: 0.9711599349975586
11 3.095894e+03 1.998258e+02
* time: 1.1093218326568604
12 2.988214e+03 4.366433e+02
* time: 1.2924718856811523
13 2.896081e+03 5.505943e+02
* time: 1.815337896347046
14 2.652467e+03 7.300323e+02
* time: 5.8744797706604
15 2.560937e+03 6.973661e+02
* time: 6.554111957550049
16 2.254941e+03 2.740033e+02
* time: 6.674770832061768
17 2.222509e+03 2.034303e+02
* time: 6.7637619972229
18 2.171255e+03 2.449580e+02
* time: 6.894490003585815
19 2.024532e+03 1.121511e+02
* time: 7.167182922363281
20 1.993723e+03 1.042814e+02
* time: 7.303891897201538
21 1.985113e+03 8.079014e+01
* time: 7.4368767738342285
22 1.976757e+03 7.054196e+01
* time: 7.5667197704315186
23 1.969970e+03 6.070322e+01
* time: 7.692018985748291
24 1.961095e+03 6.810782e+01
* time: 7.813476800918579
25 1.947983e+03 8.116920e+01
* time: 7.935450792312622
26 1.930371e+03 8.530051e+01
* time: 8.055539846420288
27 1.910209e+03 6.993170e+01
* time: 8.183430910110474
28 1.899107e+03 3.362640e+01
* time: 8.343338966369629
29 1.898022e+03 2.642220e+01
* time: 8.501652002334595
30 1.897055e+03 1.213144e+01
* time: 8.673544883728027
31 1.896596e+03 7.773239e+00
* time: 8.91552186012268
32 1.896538e+03 7.997039e+00
* time: 9.066603899002075
33 1.896451e+03 8.160909e+00
* time: 9.211074829101562
34 1.896283e+03 8.237721e+00
* time: 9.357908964157104
35 1.895903e+03 1.520219e+01
* time: 9.503371000289917
36 1.895272e+03 2.358916e+01
* time: 9.65067195892334
37 1.894536e+03 2.461296e+01
* time: 9.800596952438354
38 1.893995e+03 1.546128e+01
* time: 9.94232988357544
39 1.893858e+03 6.976137e+00
* time: 10.081135988235474
40 1.893833e+03 6.019466e+00
* time: 10.220865964889526
41 1.893786e+03 3.827201e+00
* time: 10.330291986465454
42 1.893714e+03 3.323412e+00
* time: 10.418009996414185
43 1.893592e+03 3.215150e+00
* time: 10.49972677230835
44 1.893435e+03 6.534965e+00
* time: 10.646251916885376
45 1.893286e+03 7.424154e+00
* time: 10.735177993774414
46 1.893190e+03 5.552627e+00
* time: 10.852666854858398
47 1.893139e+03 3.222316e+00
* time: 10.980474948883057
48 1.893120e+03 3.015339e+00
* time: 11.104321956634521
49 1.893107e+03 3.244809e+00
* time: 11.178218841552734
50 1.893080e+03 6.163100e+00
* time: 11.253757953643799
51 1.893027e+03 9.824713e+00
* time: 11.344467878341675
52 1.892912e+03 1.390100e+01
* time: 11.417996883392334
53 1.892734e+03 1.510937e+01
* time: 11.526832818984985
54 1.892561e+03 1.008563e+01
* time: 11.60513687133789
55 1.892485e+03 3.730668e+00
* time: 11.69468879699707
56 1.892471e+03 3.380261e+00
* time: 11.842456817626953
57 1.892463e+03 3.167904e+00
* time: 11.92451786994934
58 1.892441e+03 4.152065e+00
* time: 11.996886968612671
59 1.892391e+03 7.355996e+00
* time: 12.100000858306885
60 1.892268e+03 1.195397e+01
* time: 12.224378824234009
61 1.892026e+03 1.640783e+01
* time: 12.354666948318481
62 1.891735e+03 1.593576e+01
* time: 12.474936962127686
63 1.891569e+03 8.316423e+00
* time: 12.616017818450928
64 1.891494e+03 3.948212e+00
* time: 12.748363971710205
65 1.891481e+03 3.911593e+00
* time: 12.883448839187622
66 1.891457e+03 3.875559e+00
* time: 13.024351835250854
67 1.891405e+03 3.811247e+00
* time: 13.17592477798462
68 1.891262e+03 3.657045e+00
* time: 13.313347816467285
69 1.890930e+03 4.957405e+00
* time: 13.535538911819458
70 1.890317e+03 6.657726e+00
* time: 13.665281772613525
71 1.889660e+03 6.086302e+00
* time: 13.786556959152222
72 1.889303e+03 2.270929e+00
* time: 13.903620958328247
73 1.889253e+03 7.695301e-01
* time: 14.016101837158203
74 1.889252e+03 7.382144e-01
* time: 14.126113891601562
75 1.889251e+03 7.187898e-01
* time: 14.23475694656372
76 1.889251e+03 7.215047e-01
* time: 14.302454948425293
77 1.889250e+03 7.235155e-01
* time: 14.37543797492981
78 1.889249e+03 7.246818e-01
* time: 14.445417881011963
79 1.889244e+03 7.257796e-01
* time: 14.529278993606567
80 1.889233e+03 7.198190e-01
* time: 14.62271785736084
81 1.889204e+03 1.089029e+00
* time: 14.699723958969116
82 1.889142e+03 1.801601e+00
* time: 14.840205907821655
83 1.889043e+03 2.967917e+00
* time: 14.917443990707397
84 1.888889e+03 2.965856e+00
* time: 15.037205934524536
85 1.888705e+03 5.933557e-01
* time: 15.168643951416016
86 1.888655e+03 9.577696e-01
* time: 15.285396814346313
87 1.888582e+03 1.498494e+00
* time: 15.404414892196655
88 1.888533e+03 1.502753e+00
* time: 15.517448902130127
89 1.888490e+03 1.184664e+00
* time: 15.632523775100708
90 1.888480e+03 6.684517e-01
* time: 15.741469860076904
91 1.888476e+03 3.680034e-01
* time: 15.84157681465149
92 1.888476e+03 4.720040e-01
* time: 15.939613819122314
93 1.888476e+03 4.768646e-01
* time: 16.022581815719604
94 1.888475e+03 4.736673e-01
* time: 16.116747856140137
95 1.888475e+03 4.552764e-01
* time: 16.274987936019897
96 1.888474e+03 5.193733e-01
* time: 16.394551992416382
97 1.888473e+03 8.850112e-01
* time: 16.515661001205444
98 1.888468e+03 1.461600e+00
* time: 16.637179851531982
99 1.888458e+03 2.209127e+00
* time: 16.745754957199097
100 1.888437e+03 2.961239e+00
* time: 16.87685489654541
101 1.888407e+03 2.978463e+00
* time: 17.000081777572632
102 1.888384e+03 1.707201e+00
* time: 17.1226749420166
103 1.888381e+03 6.199354e-01
* time: 17.244938850402832
104 1.888380e+03 5.170909e-01
* time: 17.369913816452026
105 1.888378e+03 1.037408e-01
* time: 17.501852989196777
106 1.888378e+03 8.473247e-02
* time: 17.623881816864014
107 1.888378e+03 8.364978e-02
* time: 17.74494481086731
108 1.888378e+03 8.080446e-02
* time: 17.937461853027344
109 1.888378e+03 7.873905e-02
* time: 18.055684804916382
110 1.888378e+03 7.798398e-02
* time: 18.16845393180847
111 1.888378e+03 7.788170e-02
* time: 18.254133939743042
112 1.888378e+03 7.776461e-02
* time: 18.35130000114441
113 1.888378e+03 9.023662e-02
* time: 18.42649483680725
114 1.888378e+03 1.631390e-01
* time: 18.514619827270508
115 1.888378e+03 2.768731e-01
* time: 18.616392850875854
116 1.888377e+03 4.462386e-01
* time: 18.709674835205078
117 1.888377e+03 6.643297e-01
* time: 18.804588794708252
118 1.888375e+03 8.433397e-01
* time: 18.886454820632935
119 1.888374e+03 7.596814e-01
* time: 18.983346939086914
120 1.888373e+03 3.638119e-01
* time: 19.093159914016724
121 1.888372e+03 8.306034e-02
* time: 19.206685781478882
122 1.888372e+03 2.084513e-02
* time: 19.454983949661255
123 1.888372e+03 2.056419e-02
* time: 19.57160782814026
124 1.888372e+03 2.044080e-02
* time: 19.686099767684937
125 1.888372e+03 2.035196e-02
* time: 19.800028800964355
126 1.888372e+03 2.021264e-02
* time: 19.91709589958191
127 1.888372e+03 1.998163e-02
* time: 20.004896879196167
128 1.888372e+03 3.161638e-02
* time: 20.118181943893433
129 1.888372e+03 5.509218e-02
* time: 20.212055921554565
130 1.888372e+03 9.275848e-02
* time: 20.33508276939392
131 1.888372e+03 1.528749e-01
* time: 20.46320676803589
132 1.888372e+03 2.461766e-01
* time: 20.603075981140137
133 1.888372e+03 3.799362e-01
* time: 20.74925398826599
134 1.888371e+03 5.311665e-01
* time: 20.895827770233154
135 1.888369e+03 6.019039e-01
* time: 21.105527877807617
136 1.888367e+03 4.664896e-01
* time: 21.24398183822632
137 1.888366e+03 1.404934e-01
* time: 21.36772894859314
138 1.888365e+03 8.513331e-02
* time: 21.486042976379395
139 1.888364e+03 1.244077e-01
* time: 21.61138081550598
140 1.888364e+03 1.028085e-01
* time: 21.74813485145569
141 1.888364e+03 5.162231e-02
* time: 21.88889193534851
142 1.888364e+03 5.149630e-02
* time: 22.01035189628601
143 1.888364e+03 3.147284e-02
* time: 22.1330828666687
144 1.888364e+03 2.104766e-02
* time: 22.253137826919556
145 1.888364e+03 6.539151e-03
* time: 22.372392892837524
146 1.888364e+03 2.540196e-03
* time: 22.4905948638916
147 1.888364e+03 4.362661e-03
* time: 22.604761838912964
148 1.888364e+03 3.034416e-03
* time: 22.70758295059204
149 1.888364e+03 5.953892e-04
* time: 22.885556936264038FittedPumasModel
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 | 5.504 |
| 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.882 |
| 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 | 5.882 |
| 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 | 22.886 |
| 3 | Subjects | 30 |
| 4 | Fixed Parameters | 0 |
| 5 | Optimized Parameters | 13 |
| 6 | DV Active Observations | 540 |
| 7 | DV Missing Observations | 0 |
| 8 | Total Active Observations | 540 |
| 9 | Total Missing Observations | 0 |
| 10 | Likelihood Approximation | Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}} |
| 11 | LogLikelihood (LL) | -1888.36 |
| 12 | -2LL | 3776.73 |
| 13 | AIC | 3802.73 |
| 14 | BIC | 3858.52 |
| 15 | (η-shrinkage) η₁ | -0.013 |
| 16 | (η-shrinkage) η₂ | -0.012 |
| 17 | (ϵ-shrinkage) DV | 0.056 |
metrics_table outputs all of the model metrics we might be interested with respect to a certain model. That includes metadata such as estimation time, number of subjects, how many parameters were optimized and fixed, and number of observations. It also includes common model metrics like AIC, BIC, objective function value with constant (-2 loglikelihood), and so on.
We can also do an innerjoin (check our Data Wrangling Tutorials) to get all metrics into a single DataFrame:
all_metrics = innerjoin(
metrics_table(fit_base_model),
metrics_table(fit_covariate_model_wt),
metrics_table(fit_covariate_model_wt_crcl),
metrics_table(fit_covariate_model_wt_crcl_sex);
on = :Metric,
makeunique = true,
);
rename!(
all_metrics,
:Value => :Base_Model,
:Value_1 => :Covariate_Model_WT,
:Value_2 => :Covariate_Model_WT_CRCL,
:Value_3 => :Covariate_Model_WT_CRCL_SEX,
)17×5 DataFrame
| Row | Metric | Base_Model | Covariate_Model_WT | Covariate_Model_WT_CRCL | Covariate_Model_WT_CRCL_SEX |
|---|---|---|---|---|---|
| String | Any | Any | Any | Any | |
| 1 | Successful | true | true | true | true |
| 2 | Estimation Time | 5.504 | 2.882 | 5.882 | 22.886 |
| 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: 0.00011396408081054688
1 2.994747e+03 1.083462e+03
* time: 0.007758140563964844
2 2.406265e+03 1.884408e+02
* time: 0.016417980194091797
3 2.344175e+03 7.741610e+01
* time: 0.02510213851928711
4 2.323153e+03 2.907642e+01
* time: 0.03370499610900879
5 2.318222e+03 2.273295e+01
* time: 0.042832136154174805
6 2.316833e+03 1.390527e+01
* time: 0.05167698860168457
7 2.316425e+03 4.490883e+00
* time: 0.06068015098571777
8 2.316362e+03 9.374519e-01
* time: 0.06995916366577148
9 2.316356e+03 1.928785e-01
* time: 0.07933211326599121
10 2.316355e+03 3.119615e-02
* time: 0.0887911319732666
11 2.316355e+03 6.215513e-03
* time: 0.10011816024780273
12 2.316355e+03 8.313206e-04
* time: 0.11786103248596191FittedPumasModel
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.