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.0279238224029541
1 3.110315e+03 9.706222e+02
* time: 0.6026949882507324
2 2.831659e+03 7.817006e+02
* time: 0.6305928230285645
3 2.405281e+03 2.923696e+02
* time: 0.6802480220794678
4 2.370406e+03 3.032286e+02
* time: 0.6976919174194336
5 2.313631e+03 3.126188e+02
* time: 0.7147059440612793
6 2.263986e+03 2.946697e+02
* time: 0.7301239967346191
7 2.160182e+03 1.917599e+02
* time: 0.7728569507598877
8 2.112467e+03 1.588288e+02
* time: 0.7934348583221436
9 2.090339e+03 1.108334e+02
* time: 0.8083679676055908
10 2.078171e+03 8.108027e+01
* time: 0.8231790065765381
11 2.074517e+03 7.813928e+01
* time: 0.846707820892334
12 2.066270e+03 7.044632e+01
* time: 0.8611040115356445
13 2.049660e+03 1.062584e+02
* time: 0.8750689029693604
14 2.021965e+03 1.130570e+02
* time: 0.8894739151000977
15 1.994936e+03 7.825801e+01
* time: 0.9128170013427734
16 1.979337e+03 5.318263e+01
* time: 0.927901029586792
17 1.972141e+03 6.807046e+01
* time: 0.9424159526824951
18 1.967973e+03 7.896361e+01
* time: 0.9568078517913818
19 1.962237e+03 8.343757e+01
* time: 0.9799609184265137
20 1.952791e+03 5.565304e+01
* time: 0.995589017868042
21 1.935857e+03 3.923284e+01
* time: 1.0107738971710205
22 1.926254e+03 5.749643e+01
* time: 1.0256268978118896
23 1.922144e+03 4.306225e+01
* time: 1.0482258796691895
24 1.911566e+03 4.810496e+01
* time: 1.0639898777008057
25 1.906464e+03 4.324267e+01
* time: 1.0788748264312744
26 1.905339e+03 1.207954e+01
* time: 1.0924458503723145
27 1.905092e+03 1.093479e+01
* time: 1.106562852859497
28 1.904957e+03 1.057034e+01
* time: 1.1294119358062744
29 1.904875e+03 1.060882e+01
* time: 1.142657995223999
30 1.904459e+03 1.031525e+01
* time: 1.1562118530273438
31 1.903886e+03 1.041179e+01
* time: 1.169867992401123
32 1.903313e+03 1.135672e+01
* time: 1.1927168369293213
33 1.903057e+03 1.075683e+01
* time: 1.2063839435577393
34 1.902950e+03 1.091295e+01
* time: 1.2195990085601807
35 1.902887e+03 1.042409e+01
* time: 1.2324159145355225
36 1.902640e+03 9.213043e+00
* time: 1.255086898803711
37 1.902364e+03 9.519321e+00
* time: 1.2693078517913818
38 1.902156e+03 5.590984e+00
* time: 1.2828140258789062
39 1.902094e+03 1.811898e+00
* time: 1.2959198951721191
40 1.902086e+03 1.644722e+00
* time: 1.3088529109954834
41 1.902084e+03 1.653520e+00
* time: 1.330941915512085
42 1.902074e+03 1.720184e+00
* time: 1.3446288108825684
43 1.902055e+03 2.619061e+00
* time: 1.3576648235321045
44 1.902015e+03 3.519729e+00
* time: 1.3706889152526855
45 1.901962e+03 3.403372e+00
* time: 1.3926680088043213
46 1.901924e+03 1.945644e+00
* time: 1.4064009189605713
47 1.901914e+03 6.273342e-01
* time: 1.419510841369629
48 1.901913e+03 5.374557e-01
* time: 1.4322710037231445
49 1.901913e+03 5.710007e-01
* time: 1.4447009563446045
50 1.901913e+03 6.091390e-01
* time: 1.4661178588867188
51 1.901912e+03 6.692417e-01
* time: 1.47940993309021
52 1.901909e+03 7.579620e-01
* time: 1.4922230243682861
53 1.901903e+03 8.798211e-01
* time: 1.5050609111785889
54 1.901889e+03 1.002981e+00
* time: 1.518035888671875
55 1.901864e+03 1.495512e+00
* time: 1.5404648780822754
56 1.901837e+03 1.664621e+00
* time: 1.5532889366149902
57 1.901819e+03 8.601119e-01
* time: 1.5661308765411377
58 1.901815e+03 4.525470e-02
* time: 1.5790019035339355
59 1.901815e+03 1.294280e-02
* time: 1.601274013519287
60 1.901815e+03 2.876567e-03
* time: 1.6138548851013184
61 1.901815e+03 2.876567e-03
* time: 1.640523910522461
62 1.901815e+03 2.876567e-03
* time: 1.6597769260406494FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -1901.815
Number of subjects: 30
Number of parameters: Fixed Optimized
0 8
Observation records: Active Missing
DV: 540 0
Total: 540 0
---------------------
Estimate
---------------------
tvcl 0.1542
tvvc 4.5856
tvq 0.042341
tvvp 3.7082
Ω₁,₁ 0.26467
Ω₂,₂ 0.10627
σ_add 0.032183
σ_prop 0.061005
---------------------1.4 Step 3 - Covariate Model
The third step of our covariate model building workflow is to actually develop one or more covariate models. Let’s develop three covariate models:
- allometric scaling based on weight
- clearance effect based on creatinine clearance
- separated error model based on sex
To include covariates in a Pumas model we need to first include them in the @covariates block. Then, we are free to use them inside the @pre block
Here’s our covariate model with allometric scaling based on weight:
Tip
When building covariate models, unlike in NONMEM, it is highly recommended to derive covariates or perform any required transformations or scaling as a data wrangling step and pass the derived/transformed into read_pumas as a covariate. In this particular case, it is easy to create two columns in the original data as:
@rtransform! pkdata :ASWT_CL = (:WT / 70)^0.75
@rtransform! pkdata :ASWT_Vc = (:WT / 70)Then, you bring these covariates into read_pumas and use them directly on CL and Vc. This will avoid computations of your transformations during the model fit.
covariate_model_wt = @model begin
@param begin
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvq ∈ RealDomain(; lower = 0)
tvvp ∈ RealDomain(; lower = 0)
Ω ∈ PDiagDomain(2)
σ_add ∈ RealDomain(; lower = 0)
σ_prop ∈ RealDomain(; lower = 0)
end
@random begin
η ~ MvNormal(Ω)
end
@covariates begin
WT
end
@pre begin
CL = tvcl * (WT / 70)^0.75 * exp(η[1])
Vc = tvvc * (WT / 70) * exp(η[2])
Q = tvq
Vp = tvvp
end
@dynamics Central1Periph1
@derived begin
cp := @. Central / Vc
DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
end
endPumasModel
Parameters: tvcl, tvvc, tvq, tvvp, Ω, σ_add, σ_prop
Random effects: η
Covariates: WT
Dynamical system variables: Central, Peripheral
Dynamical system type: Closed form
Derived: DV
Observed: DVOnce we included the WT covariate in the @covariates block we can use the WT values inside the @pre block. For both clearance (CL) and volume of the central compartment (Vc), we are allometric scaling by the WT value by the mean weight 70 and, in the case of CL using an allometric exponent with value 0.75.
Let’s fit our covariate_model_wt. Notice that we have not added any new parameters inside the @param block, thus, we can use the same iparams_base_model initial parameters values’ list:
fit_covariate_model_wt = fit(covariate_model_wt, pop, iparams_base_model, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 7.695401e+03 4.898919e+03
* time: 7.081031799316406e-5
1 2.846050e+03 1.128657e+03
* time: 0.025653839111328125
2 2.472982e+03 7.008264e+02
* time: 0.043434858322143555
3 2.180690e+03 2.312704e+02
* time: 0.061018943786621094
4 2.125801e+03 1.862929e+02
* time: 0.09877300262451172
5 2.103173e+03 1.320946e+02
* time: 0.11426782608032227
6 2.091136e+03 1.103035e+02
* time: 0.12880992889404297
7 2.081443e+03 1.091133e+02
* time: 0.14424395561218262
8 2.071793e+03 7.197675e+01
* time: 0.17983794212341309
9 2.062706e+03 7.623310e+01
* time: 0.19489192962646484
10 2.057515e+03 6.885476e+01
* time: 0.20925188064575195
11 2.051133e+03 6.368504e+01
* time: 0.22362089157104492
12 2.038626e+03 7.730243e+01
* time: 0.24757790565490723
13 2.019352e+03 1.136864e+02
* time: 0.2624669075012207
14 1.997136e+03 1.005748e+02
* time: 0.27692294120788574
15 1.983023e+03 6.831478e+01
* time: 0.29218196868896484
16 1.977700e+03 5.649783e+01
* time: 0.3159308433532715
17 1.974583e+03 6.357642e+01
* time: 0.33077001571655273
18 1.967292e+03 7.658974e+01
* time: 0.3458068370819092
19 1.951045e+03 6.130573e+01
* time: 0.36245083808898926
20 1.935868e+03 4.845839e+01
* time: 0.3866608142852783
21 1.929356e+03 6.325782e+01
* time: 0.40264892578125
22 1.925187e+03 3.142245e+01
* time: 0.4171409606933594
23 1.923733e+03 4.623400e+01
* time: 0.43126893043518066
24 1.918498e+03 5.347738e+01
* time: 0.4456198215484619
25 1.912383e+03 5.849125e+01
* time: 0.47100090980529785
26 1.905510e+03 3.254038e+01
* time: 0.4879150390625
27 1.903629e+03 2.905618e+01
* time: 0.5019838809967041
28 1.902833e+03 2.907696e+01
* time: 0.5166079998016357
29 1.902447e+03 2.746037e+01
* time: 0.5396809577941895
30 1.899399e+03 1.930949e+01
* time: 0.5543379783630371
31 1.898705e+03 1.186361e+01
* time: 0.5684778690338135
32 1.898505e+03 1.050402e+01
* time: 0.5826258659362793
33 1.898474e+03 1.042186e+01
* time: 0.6046469211578369
34 1.897992e+03 1.238729e+01
* time: 0.6183979511260986
35 1.897498e+03 1.729368e+01
* time: 0.6319069862365723
36 1.896954e+03 1.472554e+01
* time: 0.6454160213470459
37 1.896744e+03 5.852825e+00
* time: 0.6681599617004395
38 1.896705e+03 1.171353e+00
* time: 0.6817889213562012
39 1.896704e+03 1.216117e+00
* time: 0.6949539184570312
40 1.896703e+03 1.230336e+00
* time: 0.7077620029449463
41 1.896698e+03 1.250715e+00
* time: 0.7208259105682373
42 1.896688e+03 1.449552e+00
* time: 0.7429599761962891
43 1.896666e+03 2.533300e+00
* time: 0.7568008899688721
44 1.896631e+03 3.075537e+00
* time: 0.7700760364532471
45 1.896599e+03 2.139797e+00
* time: 0.7833578586578369
46 1.896587e+03 6.636030e-01
* time: 0.7967979907989502
47 1.896585e+03 6.303517e-01
* time: 0.8194818496704102
48 1.896585e+03 5.995265e-01
* time: 0.8323960304260254
49 1.896584e+03 5.844159e-01
* time: 0.84499192237854
50 1.896583e+03 6.083858e-01
* time: 0.8578739166259766
51 1.896579e+03 8.145327e-01
* time: 0.8800079822540283
52 1.896570e+03 1.293490e+00
* time: 0.8937959671020508
53 1.896549e+03 1.877705e+00
* time: 0.9069709777832031
54 1.896513e+03 2.217392e+00
* time: 0.9200129508972168
55 1.896482e+03 1.658148e+00
* time: 0.9330549240112305
56 1.896466e+03 5.207218e-01
* time: 0.9555950164794922
57 1.896463e+03 1.177468e-01
* time: 0.969005823135376
58 1.896463e+03 1.678937e-02
* time: 0.9825289249420166
59 1.896463e+03 2.666636e-03
* time: 0.9942989349365234
60 1.896463e+03 2.666636e-03
* time: 1.0295848846435547
61 1.896463e+03 2.666636e-03
* time: 1.0521998405456543FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -1896.4632
Number of subjects: 30
Number of parameters: Fixed Optimized
0 8
Observation records: Active Missing
DV: 540 0
Total: 540 0
---------------------
Estimate
---------------------
tvcl 0.13915
tvvc 3.9754
tvq 0.041988
tvvp 3.5722
Ω₁,₁ 0.23874
Ω₂,₂ 0.081371
σ_add 0.032174
σ_prop 0.061012
---------------------We can definitely see that, despite not increasing the parameters of the model, our loglikelihood is higher (implying lower objective function value). Furthermore, our additive and combined error values remain almost the same while our Ωs decreased for CL and Vc. This implies that the WT covariate is definitely assisting in a better model fit by capturing the variability in CL and Vc. We’ll compare models in the next step.
Let’s now try to incorporate into this model creatinine clearance (CRCL) effect on clearance (CL):
Tip
Like the tip above, you can derive covariates or perform any required transformations or scaling as a data wrangling step and pass the derived/transformed into read_pumas as a covariate. In this particular case, it is easy to create three more columns in the original data as:
@rtransform! pkdata :CRCL_CL = :CRCL / 100
@rtransform! pkdata :ASWT_CL = (:WT / 70)^0.75
@rtransform! pkdata :ASWT_Vc = (:WT / 70)Then, you bring these covariates into read_pumas and use them directly on renCL, CL and Vc. This will avoid computations of your transformations during the model fit.
covariate_model_wt_crcl = @model begin
@param begin
tvvc ∈ RealDomain(; lower = 0)
tvq ∈ RealDomain(; lower = 0)
tvvp ∈ RealDomain(; lower = 0)
tvcl_hep ∈ RealDomain(; lower = 0)
tvcl_ren ∈ RealDomain(; lower = 0)
Ω ∈ PDiagDomain(2)
σ_add ∈ RealDomain(; lower = 0)
σ_prop ∈ RealDomain(; lower = 0)
dCRCL ∈ RealDomain()
end
@random begin
η ~ MvNormal(Ω)
end
@covariates begin
WT
CRCL
end
@pre begin
hepCL = tvcl_hep * (WT / 70)^0.75
renCL = tvcl_ren * (CRCL / 100)^dCRCL
CL = (hepCL + renCL) * exp(η[1])
Vc = tvvc * (WT / 70) * exp(η[2])
Q = tvq
Vp = tvvp
end
@dynamics Central1Periph1
@derived begin
cp := @. Central / Vc
DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
end
endPumasModel
Parameters: tvvc, tvq, tvvp, tvcl_hep, tvcl_ren, Ω, σ_add, σ_prop, dCRCL
Random effects: η
Covariates: WT, CRCL
Dynamical system variables: Central, Peripheral
Dynamical system type: Closed form
Derived: DV
Observed: DVIn the covariate_model_wt_crcl model we are keeping our allometric scaling on WT from before. But we are also adding a new covariate creatinine clearance (CRCL), dividing clearance (CL) into hepatic (hepCL) and renal clearance (renCL), along with a new parameter dCRCL.
dCRCL is the exponent of the power function for the effect of creatinine clearance on renal clearance. In some models this parameter is fixed, however we’ll allow the model to estimate it.
This is a good example on how to add covariate coefficients such as dCRCL in any Pumas covariate model. Now, let’s fit this model. Note that we need a new initial parameters values’ list since the previous one we used doesn’t include dCRCL, tvcl_hep or tvcl_ren:
iparams_covariate_model_wt_crcl = (;
tvvc = 5,
tvcl_hep = 0.01,
tvcl_ren = 0.01,
tvq = 0.01,
tvvp = 10,
Ω = Diagonal([0.01, 0.01]),
σ_add = 0.1,
σ_prop = 0.1,
dCRCL = 0.9,
)(tvvc = 5,
tvcl_hep = 0.01,
tvcl_ren = 0.01,
tvq = 0.01,
tvvp = 10,
Ω = [0.01 0.0; 0.0 0.01],
σ_add = 0.1,
σ_prop = 0.1,
dCRCL = 0.9,)fit_covariate_model_wt_crcl =
fit(covariate_model_wt_crcl, pop, iparams_covariate_model_wt_crcl, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 8.554152e+03 5.650279e+03
* time: 7.700920104980469e-5
1 3.572050e+03 1.302046e+03
* time: 0.032488107681274414
2 3.266947e+03 5.384244e+02
* time: 0.09417104721069336
3 3.150635e+03 1.918079e+02
* time: 0.11456608772277832
4 3.108122e+03 1.277799e+02
* time: 0.1335740089416504
5 3.091358e+03 8.838080e+01
* time: 0.176192045211792
6 3.082997e+03 8.321689e+01
* time: 0.1950209140777588
7 3.076379e+03 8.167702e+01
* time: 0.2130570411682129
8 3.067422e+03 1.177822e+02
* time: 0.23166799545288086
9 3.048580e+03 2.526969e+02
* time: 0.2726891040802002
10 2.993096e+03 6.325396e+02
* time: 0.3014490604400635
11 2.965744e+03 7.634718e+02
* time: 0.36688899993896484
12 2.921212e+03 9.704020e+02
* time: 0.4105050563812256
13 2.553649e+03 6.642510e+02
* time: 0.4382600784301758
14 2.319495e+03 3.204711e+02
* time: 0.477100133895874
15 2.183040e+03 2.174905e+02
* time: 0.519503116607666
16 2.157621e+03 3.150983e+02
* time: 0.5383810997009277
17 2.132395e+03 2.847614e+02
* time: 0.5667569637298584
18 2.084799e+03 1.563370e+02
* time: 0.5855779647827148
19 2.071497e+03 1.006429e+02
* time: 0.6030600070953369
20 2.064983e+03 9.753313e+01
* time: 0.620703935623169
21 2.048289e+03 9.230405e+01
* time: 0.6491689682006836
22 2.020938e+03 1.292359e+02
* time: 0.6684269905090332
23 1.983888e+03 2.284990e+02
* time: 0.6870651245117188
24 1.962132e+03 1.220188e+02
* time: 0.714818000793457
25 1.945947e+03 1.035894e+02
* time: 0.7324299812316895
26 1.917782e+03 8.246698e+01
* time: 0.7502589225769043
27 1.905967e+03 5.408054e+01
* time: 0.7680189609527588
28 1.898569e+03 2.172222e+01
* time: 0.7961099147796631
29 1.897473e+03 1.689350e+01
* time: 0.8132381439208984
30 1.897019e+03 1.666689e+01
* time: 0.8301820755004883
31 1.896796e+03 1.699751e+01
* time: 0.8570709228515625
32 1.896559e+03 1.645900e+01
* time: 0.8741130828857422
33 1.896223e+03 1.415504e+01
* time: 0.8909780979156494
34 1.895773e+03 1.630081e+01
* time: 0.9078121185302734
35 1.895309e+03 1.723930e+01
* time: 0.9354250431060791
36 1.895004e+03 1.229983e+01
* time: 0.9527490139007568
37 1.894871e+03 5.385102e+00
* time: 0.9697849750518799
38 1.894827e+03 3.465463e+00
* time: 0.9966869354248047
39 1.894816e+03 3.387474e+00
* time: 1.013477087020874
40 1.894807e+03 3.295388e+00
* time: 1.0298500061035156
41 1.894786e+03 3.089194e+00
* time: 1.046112060546875
42 1.894737e+03 2.928080e+00
* time: 1.0729360580444336
43 1.894624e+03 3.088723e+00
* time: 1.0898890495300293
44 1.894413e+03 3.493791e+00
* time: 1.106734037399292
45 1.894127e+03 3.142865e+00
* time: 1.133085012435913
46 1.893933e+03 2.145253e+00
* time: 1.1505990028381348
47 1.893888e+03 2.172800e+00
* time: 1.1684250831604004
48 1.893880e+03 2.180509e+00
* time: 1.1845860481262207
49 1.893876e+03 2.134257e+00
* time: 1.2109260559082031
50 1.893868e+03 2.032354e+00
* time: 1.2275280952453613
51 1.893846e+03 1.760874e+00
* time: 1.2437870502471924
52 1.893796e+03 1.779016e+00
* time: 1.2696559429168701
53 1.893694e+03 2.018956e+00
* time: 1.2870440483093262
54 1.893559e+03 2.366854e+00
* time: 1.303666114807129
55 1.893474e+03 3.690142e+00
* time: 1.3203320503234863
56 1.893446e+03 3.675109e+00
* time: 1.3467950820922852
57 1.893439e+03 3.426419e+00
* time: 1.3634729385375977
58 1.893429e+03 3.183164e+00
* time: 1.3797481060028076
59 1.893398e+03 2.695171e+00
* time: 1.3960261344909668
60 1.893328e+03 2.753548e+00
* time: 1.4226670265197754
61 1.893169e+03 3.589748e+00
* time: 1.4394960403442383
62 1.892920e+03 3.680718e+00
* time: 1.4563019275665283
63 1.892667e+03 2.568107e+00
* time: 1.48288893699646
64 1.892514e+03 1.087910e+00
* time: 1.5002789497375488
65 1.892493e+03 3.287296e-01
* time: 1.5166590213775635
66 1.892492e+03 2.967465e-01
* time: 1.5329210758209229
67 1.892492e+03 3.020682e-01
* time: 1.5590710639953613
68 1.892491e+03 3.034704e-01
* time: 1.5745790004730225
69 1.892491e+03 3.091846e-01
* time: 1.5901849269866943
70 1.892491e+03 3.224170e-01
* time: 1.6152980327606201
71 1.892490e+03 6.494197e-01
* time: 1.6320760250091553
72 1.892488e+03 1.115188e+00
* time: 1.648164987564087
73 1.892483e+03 1.838833e+00
* time: 1.665349006652832
74 1.892472e+03 2.765371e+00
* time: 1.6913931369781494
75 1.892452e+03 3.463807e+00
* time: 1.7084369659423828
76 1.892431e+03 2.805270e+00
* time: 1.7247021198272705
77 1.892411e+03 5.758916e-01
* time: 1.7411830425262451
78 1.892410e+03 1.434041e-01
* time: 1.7676920890808105
79 1.892409e+03 1.639246e-01
* time: 1.783797025680542
80 1.892409e+03 1.145856e-01
* time: 1.7992839813232422
81 1.892409e+03 3.966861e-02
* time: 1.824181079864502
82 1.892409e+03 3.550808e-02
* time: 1.8402209281921387
83 1.892409e+03 3.456241e-02
* time: 1.8550729751586914
84 1.892409e+03 3.114018e-02
* time: 1.8698971271514893
85 1.892409e+03 4.080806e-02
* time: 1.8945250511169434
86 1.892409e+03 6.722726e-02
* time: 1.9106299877166748
87 1.892409e+03 1.006791e-01
* time: 1.9259920120239258
88 1.892409e+03 1.303988e-01
* time: 1.941357135772705
89 1.892409e+03 1.228919e-01
* time: 1.9663100242614746
90 1.892409e+03 6.433813e-02
* time: 1.9821829795837402
91 1.892409e+03 1.314164e-02
* time: 1.9976129531860352
92 1.892409e+03 4.929931e-04
* time: 2.0126209259033203FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -1892.409
Number of subjects: 30
Number of parameters: Fixed Optimized
0 10
Observation records: Active Missing
DV: 540 0
Total: 540 0
-----------------------
Estimate
-----------------------
tvvc 3.9757
tvq 0.042177
tvvp 3.6434
tvcl_hep 0.058572
tvcl_ren 0.1337
Ω₁,₁ 0.18299
Ω₂,₂ 0.081353
σ_add 0.032174
σ_prop 0.06101
dCRCL 1.1331
-----------------------As before, our loglikelihood is higher (implying lower objective function value). Furthermore, our additive and combined error values remain almost the same while our Ω on CL, Ω₁,₁, decreased. This implies that the CRCL covariate with an estimated exponent dCRCL is definitely assisting in a better model fit.
Finally let’s include a separated CL model based on sex as a third covariate model:
covariate_model_wt_crcl_sex = @model begin
@param begin
tvvc ∈ RealDomain(; lower = 0)
tvq ∈ RealDomain(; lower = 0)
tvvp ∈ RealDomain(; lower = 0)
tvcl_hep_M ∈ RealDomain(; lower = 0)
tvcl_hep_F ∈ RealDomain(; lower = 0)
tvcl_ren_M ∈ RealDomain(; lower = 0)
tvcl_ren_F ∈ RealDomain(; lower = 0)
Ω ∈ PDiagDomain(2)
σ_add ∈ RealDomain(; lower = 0)
σ_prop ∈ RealDomain(; lower = 0)
dCRCL_M ∈ RealDomain()
dCRCL_F ∈ RealDomain()
end
@random begin
η ~ MvNormal(Ω)
end
@covariates begin
WT
CRCL
SEX
end
@pre begin
tvcl_hep = ifelse(SEX == "M", tvcl_hep_M, tvcl_hep_F)
tvcl_ren = ifelse(SEX == "M", tvcl_ren_M, tvcl_ren_F)
dCRCL = ifelse(SEX == "M", dCRCL_M, dCRCL_F)
hepCL = tvcl_hep * (WT / 70)^0.75
renCL = tvcl_ren * (CRCL / 100)^dCRCL
CL = (hepCL + renCL) * exp(η[1])
Vc = tvvc * (WT / 70) * exp(η[2])
Q = tvq
Vp = tvvp
end
@dynamics Central1Periph1
@derived begin
cp := @. Central / Vc
DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
end
endPumasModel
Parameters: tvvc, tvq, tvvp, tvcl_hep_M, tvcl_hep_F, tvcl_ren_M, tvcl_ren_F, Ω, σ_add, σ_prop, dCRCL_M, dCRCL_F
Random effects: η
Covariates: WT, CRCL, SEX
Dynamical system variables: Central, Peripheral
Dynamical system type: Closed form
Derived: DV
Observed: DVIn the covariate_model_wt_crcl_sex model we are keeping our allometric scaling on WT, the CRCL effect on renCL, and the breakdown of CL into hepCL and renCL from before. However we are separating them with different values by sex. Hence, we have a new covariate SEX and six new parameters in the @param block by expanding tvcl_hep, tvcl_ren, and dCRCL into male (suffix M) and female (suffix F).
This is a good example on how to add create binary values based on covariate values such as SEX inside the @pre block with the ifelse function. Now, let’s fit this model. Note that we need a new initial parameters values’ list since the previous one we used had a single tvcl_hep, tvcl_ren, and dCRCL:
iparams_covariate_model_wt_crcl_sex = (;
tvvc = 5,
tvcl_hep_M = 0.01,
tvcl_hep_F = 0.01,
tvcl_ren_M = 0.01,
tvcl_ren_F = 0.01,
tvq = 0.01,
tvvp = 10,
Ω = Diagonal([0.01, 0.01]),
σ_add = 0.1,
σ_prop = 0.1,
dCRCL_M = 0.9,
dCRCL_F = 0.9,
)(tvvc = 5,
tvcl_hep_M = 0.01,
tvcl_hep_F = 0.01,
tvcl_ren_M = 0.01,
tvcl_ren_F = 0.01,
tvq = 0.01,
tvvp = 10,
Ω = [0.01 0.0; 0.0 0.01],
σ_add = 0.1,
σ_prop = 0.1,
dCRCL_M = 0.9,
dCRCL_F = 0.9,)fit_covariate_model_wt_crcl_sex =
fit(covariate_model_wt_crcl_sex, pop, iparams_covariate_model_wt_crcl_sex, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 8.554152e+03 5.650279e+03
* time: 6.818771362304688e-5
1 3.641387e+03 1.432080e+03
* time: 0.033859968185424805
2 3.290450e+03 5.274782e+02
* time: 0.057485103607177734
3 3.185512e+03 2.173676e+02
* time: 0.10406804084777832
4 3.143009e+03 1.479653e+02
* time: 0.12502408027648926
5 3.128511e+03 8.980031e+01
* time: 0.1683039665222168
6 3.123188e+03 5.033125e+01
* time: 0.1888871192932129
7 3.120794e+03 4.279722e+01
* time: 0.20848703384399414
8 3.118627e+03 3.971051e+01
* time: 0.23975801467895508
9 3.115300e+03 8.456587e+01
* time: 0.26175808906555176
10 3.109353e+03 1.350354e+02
* time: 0.2932701110839844
11 3.095894e+03 1.998258e+02
* time: 0.3161160945892334
12 2.988214e+03 4.366433e+02
* time: 0.34598708152770996
13 2.896081e+03 5.505943e+02
* time: 0.4241831302642822
14 2.652467e+03 7.300323e+02
* time: 0.8599610328674316
15 2.560937e+03 6.973661e+02
* time: 0.9765291213989258
16 2.254941e+03 2.740033e+02
* time: 1.001568078994751
17 2.222509e+03 2.034303e+02
* time: 1.034703016281128
18 2.171255e+03 2.449580e+02
* time: 1.0576250553131104
19 2.024532e+03 1.121511e+02
* time: 1.0790941715240479
20 1.993723e+03 1.042814e+02
* time: 1.1101691722869873
21 1.985113e+03 8.079014e+01
* time: 1.1306531429290771
22 1.976757e+03 7.054196e+01
* time: 1.1506831645965576
23 1.969970e+03 6.070322e+01
* time: 1.181730031967163
24 1.961095e+03 6.810782e+01
* time: 1.2013921737670898
25 1.947983e+03 8.116920e+01
* time: 1.2312390804290771
26 1.930371e+03 8.530051e+01
* time: 1.2533211708068848
27 1.910209e+03 6.993170e+01
* time: 1.2731549739837646
28 1.899107e+03 3.362640e+01
* time: 1.3043529987335205
29 1.898022e+03 2.642220e+01
* time: 1.3238871097564697
30 1.897055e+03 1.213144e+01
* time: 1.3428239822387695
31 1.896596e+03 7.773239e+00
* time: 1.3733429908752441
32 1.896538e+03 7.997039e+00
* time: 1.3923230171203613
33 1.896451e+03 8.160909e+00
* time: 1.421138048171997
34 1.896283e+03 8.237721e+00
* time: 1.4410719871520996
35 1.895903e+03 1.520219e+01
* time: 1.4602181911468506
36 1.895272e+03 2.358916e+01
* time: 1.4899241924285889
37 1.894536e+03 2.461296e+01
* time: 1.5104851722717285
38 1.893995e+03 1.546128e+01
* time: 1.5296471118927002
39 1.893858e+03 6.976137e+00
* time: 1.5593390464782715
40 1.893833e+03 6.019466e+00
* time: 1.578571081161499
41 1.893786e+03 3.827201e+00
* time: 1.5973241329193115
42 1.893714e+03 3.323412e+00
* time: 1.6274011135101318
43 1.893592e+03 3.215150e+00
* time: 1.6469380855560303
44 1.893435e+03 6.534965e+00
* time: 1.6762821674346924
45 1.893286e+03 7.424154e+00
* time: 1.696608066558838
46 1.893190e+03 5.552627e+00
* time: 1.715703010559082
47 1.893139e+03 3.222316e+00
* time: 1.7450101375579834
48 1.893120e+03 3.015339e+00
* time: 1.7658531665802002
49 1.893107e+03 3.244809e+00
* time: 1.7840640544891357
50 1.893080e+03 6.163100e+00
* time: 1.8127729892730713
51 1.893027e+03 9.824713e+00
* time: 1.8320059776306152
52 1.892912e+03 1.390100e+01
* time: 1.8505380153656006
53 1.892734e+03 1.510937e+01
* time: 1.8801581859588623
54 1.892561e+03 1.008563e+01
* time: 1.8993110656738281
55 1.892485e+03 3.730668e+00
* time: 1.9179580211639404
56 1.892471e+03 3.380261e+00
* time: 1.947432041168213
57 1.892463e+03 3.167904e+00
* time: 1.9659621715545654
58 1.892441e+03 4.152065e+00
* time: 1.9946630001068115
59 1.892391e+03 7.355996e+00
* time: 2.0144331455230713
60 1.892268e+03 1.195397e+01
* time: 2.033198118209839
61 1.892026e+03 1.640783e+01
* time: 2.0625619888305664
62 1.891735e+03 1.593576e+01
* time: 2.0827159881591797
63 1.891569e+03 8.316423e+00
* time: 2.101708173751831
64 1.891494e+03 3.948212e+00
* time: 2.1309990882873535
65 1.891481e+03 3.911593e+00
* time: 2.150083065032959
66 1.891457e+03 3.875559e+00
* time: 2.1685101985931396
67 1.891405e+03 3.811247e+00
* time: 2.1978611946105957
68 1.891262e+03 3.657045e+00
* time: 2.2169201374053955
69 1.890930e+03 4.957405e+00
* time: 2.2358691692352295
70 1.890317e+03 6.657726e+00
* time: 2.267695188522339
71 1.889660e+03 6.086302e+00
* time: 2.2872371673583984
72 1.889303e+03 2.270929e+00
* time: 2.316493034362793
73 1.889253e+03 7.695301e-01
* time: 2.3361010551452637
74 1.889252e+03 7.382144e-01
* time: 2.3545501232147217
75 1.889251e+03 7.187898e-01
* time: 2.383234977722168
76 1.889251e+03 7.215047e-01
* time: 2.4019830226898193
77 1.889250e+03 7.235155e-01
* time: 2.4199981689453125
78 1.889249e+03 7.246818e-01
* time: 2.4486000537872314
79 1.889244e+03 7.257796e-01
* time: 2.467674970626831
80 1.889233e+03 7.198190e-01
* time: 2.4858431816101074
81 1.889204e+03 1.089029e+00
* time: 2.5151519775390625
82 1.889142e+03 1.801601e+00
* time: 2.5346031188964844
83 1.889043e+03 2.967917e+00
* time: 2.5534210205078125
84 1.888889e+03 2.965856e+00
* time: 2.5834860801696777
85 1.888705e+03 5.933554e-01
* time: 2.6027331352233887
86 1.888655e+03 9.577699e-01
* time: 2.6215851306915283
87 1.888582e+03 1.498494e+00
* time: 2.65134596824646
88 1.888533e+03 1.502750e+00
* time: 2.6701090335845947
89 1.888490e+03 1.184664e+00
* time: 2.699054002761841
90 1.888480e+03 6.684513e-01
* time: 2.7184300422668457
91 1.888476e+03 3.680030e-01
* time: 2.7368640899658203
92 1.888476e+03 4.720039e-01
* time: 2.7667391300201416
93 1.888476e+03 4.768646e-01
* time: 2.7856791019439697
94 1.888475e+03 4.736674e-01
* time: 2.8038671016693115
95 1.888475e+03 4.552766e-01
* time: 2.8329291343688965
96 1.888474e+03 5.193719e-01
* time: 2.8519439697265625
97 1.888473e+03 8.850088e-01
* time: 2.870223045349121
98 1.888468e+03 1.461597e+00
* time: 2.899415969848633
99 1.888458e+03 2.209123e+00
* time: 2.918661117553711
100 1.888437e+03 2.961234e+00
* time: 2.9378151893615723
101 1.888407e+03 2.978462e+00
* time: 2.9676690101623535
102 1.888384e+03 1.707197e+00
* time: 2.9866549968719482
103 1.888381e+03 6.198730e-01
* time: 3.0054140090942383
104 1.888380e+03 5.171201e-01
* time: 3.035914182662964
105 1.888378e+03 1.037261e-01
* time: 3.0552830696105957
106 1.888378e+03 8.473257e-02
* time: 3.0837180614471436
107 1.888378e+03 8.364956e-02
* time: 3.1029460430145264
108 1.888378e+03 8.080438e-02
* time: 3.1213481426239014
109 1.888378e+03 7.873896e-02
* time: 3.149738073348999
110 1.888378e+03 7.798398e-02
* time: 3.168653964996338
111 1.888378e+03 7.788171e-02
* time: 3.1863410472869873
112 1.888378e+03 7.776461e-02
* time: 3.214738130569458
113 1.888378e+03 9.023533e-02
* time: 3.233631134033203
114 1.888378e+03 1.631356e-01
* time: 3.253178119659424
115 1.888378e+03 2.768664e-01
* time: 3.2821950912475586
116 1.888377e+03 4.462262e-01
* time: 3.3021559715270996
117 1.888377e+03 6.643078e-01
* time: 3.3210999965667725
118 1.888375e+03 8.433023e-01
* time: 3.350633144378662
119 1.888374e+03 7.596239e-01
* time: 3.3707051277160645
120 1.888373e+03 3.637667e-01
* time: 3.3897640705108643
121 1.888372e+03 8.304667e-02
* time: 3.4195220470428467
122 1.888372e+03 2.084518e-02
* time: 3.438681125640869
123 1.888372e+03 2.056414e-02
* time: 3.4570000171661377
124 1.888372e+03 2.044078e-02
* time: 3.4860470294952393
125 1.888372e+03 2.035197e-02
* time: 3.5039350986480713
126 1.888372e+03 2.021268e-02
* time: 3.5218050479888916
127 1.888372e+03 1.998172e-02
* time: 3.5513060092926025
128 1.888372e+03 3.162406e-02
* time: 3.5704381465911865
129 1.888372e+03 5.510549e-02
* time: 3.5888020992279053
130 1.888372e+03 9.278088e-02
* time: 3.6186230182647705
131 1.888372e+03 1.529116e-01
* time: 3.6376161575317383
132 1.888372e+03 2.462349e-01
* time: 3.656409978866577
133 1.888372e+03 3.800236e-01
* time: 3.686540126800537
134 1.888371e+03 5.312831e-01
* time: 3.7055461406707764
135 1.888369e+03 6.020265e-01
* time: 3.7248151302337646
136 1.888367e+03 4.665657e-01
* time: 3.756834030151367
137 1.888366e+03 1.404905e-01
* time: 3.7757301330566406
138 1.888365e+03 8.513244e-02
* time: 3.805065155029297
139 1.888364e+03 1.244427e-01
* time: 3.824846029281616
140 1.888364e+03 1.028331e-01
* time: 3.8435370922088623
141 1.888364e+03 5.164076e-02
* time: 3.8726611137390137
142 1.888364e+03 5.147918e-02
* time: 3.8913819789886475
143 1.888364e+03 3.147222e-02
* time: 3.9089601039886475
144 1.888364e+03 2.104481e-02
* time: 3.936915159225464
145 1.888364e+03 6.543267e-03
* time: 3.954979181289673
146 1.888364e+03 2.537332e-03
* time: 3.972551107406616
147 1.888364e+03 4.361311e-03
* time: 4.000972032546997
148 1.888364e+03 3.035139e-03
* time: 4.0191919803619385
149 1.888364e+03 5.966636e-04
* time: 4.0361950397491455FittedPumasModel
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.7174e-7
tvcl_hep_F 0.13348
tvcl_ren_M 0.19378
tvcl_ren_F 0.042211
Ω₁,₁ 0.14046
Ω₂,₂ 0.081349
σ_add 0.032171
σ_prop 0.061007
dCRCL_M 0.94821
dCRCL_F 1.9405
--------------------------As before, our loglikelihood is higher (implying lower objective function value). This is expected since we also added six new parameters to the model.
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 | 1.66 |
| 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 | 1.052 |
| 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 | 2.013 |
| 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 | 4.036 |
| 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 | 1.66 | 1.052 | 2.013 | 4.036 |
| 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.8179473717055aic(fit_covariate_model_wt_crcl_sex)3802.7275243739778We 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_every = 1, time_limit = NaN, )
)And now we pass inspect_covariate_model_wt_crcl_sex to the goodness_of_fit plotting function:
goodness_of_fit(inspect_covariate_model_wt_crcl_sex)
The idea is that the population predictions (preds) capture the general tendency of the observations while the individual predictions (ipreds) should coincide much more closely with the observations. That is exactly what we are observing in the top row subplots in our goodness of fit plot.
Regarding the bottom row, on the left we have the weighted population residuals (wres) against time, and on the right we have the weighted individual residuals (iwres) against ipreds. Here we should not see any perceived pattern, which indicates that the error model in the model has a mean 0 and constant variance. Like before, this seems to be what we are observing in our goodness of fit plot.
Hence, our covariate model with allometric scaling and different sex creatinine clearance effectw on clearance is a good candidate for our final model.
1.6 Time-Varying Covariates
Pumas can handle time-varying covariates. This happens automatically if, when parsing a dataset, read_pumas detects that covariate values change over time.
1.6.1 painord Dataset
Here’s an example with an ordinal regression using the painord dataset from PharmaDatasets.jl. :painord is our observations measuring the perceived pain in a scale from 0 to 3, which we need to have the range shifted by 1 (1 to 4). Additionally, we’ll use the concentration in plasma, :conc as a covariate. Of course, :conc varies with time, thus, it is a time-varying covariate:
painord = dataset("pumas/pain_remed")
first(painord, 5)5×8 DataFrame
| Row | id | arm | dose | time | conc | painord | dv | remed |
|---|---|---|---|---|---|---|---|---|
| Int64 | Int64 | Int64 | Float64 | Float64 | Int64 | Int64 | Int64 | |
| 1 | 1 | 2 | 20 | 0.0 | 0.0 | 3 | 0 | 0 |
| 2 | 1 | 2 | 20 | 0.5 | 1.15578 | 1 | 1 | 0 |
| 3 | 1 | 2 | 20 | 1.0 | 1.37211 | 0 | 1 | 0 |
| 4 | 1 | 2 | 20 | 1.5 | 1.30058 | 0 | 1 | 0 |
| 5 | 1 | 2 | 20 | 2.0 | 1.19195 | 1 | 1 | 0 |
@rtransform! painord :painord = :painord + 1;
describe(painord, :mean, :std, :first, :last, :eltype)8×6 DataFrame
| Row | variable | mean | std | first | last | eltype |
|---|---|---|---|---|---|---|
| Symbol | Float64 | Float64 | Real | Real | DataType | |
| 1 | id | 80.5 | 46.1992 | 1 | 160 | Int64 |
| 2 | arm | 1.5 | 1.11833 | 2 | 0 | Int64 |
| 3 | dose | 26.25 | 31.9017 | 20 | 0 | Int64 |
| 4 | time | 3.375 | 2.5183 | 0.0 | 8.0 | Float64 |
| 5 | conc | 0.93018 | 1.49902 | 0.0 | 0.0 | Float64 |
| 6 | painord | 2.50208 | 0.863839 | 4 | 4 | Int64 |
| 7 | dv | 0.508333 | 0.500061 | 0 | 0 | Int64 |
| 8 | remed | 0.059375 | 0.236387 | 0 | 0 | Int64 |
unique(painord.dose)4-element Vector{Int64}:
20
80
0
5As we can see we have 160 subjects were given either 0, 5, 20, or 80 units of a certain painkiller drug.
:conc is the drug concentration in plasma and :painord is the perceived pain in a scale from 1 to 4.
First, we’ll parse the painord dataset into a Population. Note that we’ll be using event_data=false since we do not have any dosing rows:
pop_ord =
read_pumas(painord; observations = [:painord], covariates = [:conc], event_data = false)
Note
We won’t be going into the details of the ordinal regression model in this tutorial. We highly encourage you to take a look at the Ordinal Regression Pumas Tutorial for an in-depth explanation.
We’ll build an ordinal regression model declaring :conc as a covariate. In the @derived block we’ll state the the likelihood of :painord follows a Categorical distribution:
ordinal_model = @model begin
@param begin
b₁ ∈ RealDomain(; init = 0)
b₂ ∈ RealDomain(; init = 1)
b₃ ∈ RealDomain(; init = 1)
slope ∈ RealDomain(; init = 0)
end
@covariates conc # time-varying
@pre begin
effect = slope * conc
# Logit of cumulative probabilities
lge₁ = b₁ + effect
lge₂ = lge₁ - b₂
lge₃ = lge₂ - b₃
# Probabilities of >=1 and >=2 and >=3
pge₁ = logistic(lge₁)
pge₂ = logistic(lge₂)
pge₃ = logistic(lge₃)
# Probabilities of Y=1,2,3,4
p₁ = 1.0 - pge₁
p₂ = pge₁ - pge₂
p₃ = pge₂ - pge₃
p₄ = pge₃
end
@derived begin
painord ~ @. Categorical(p₁, p₂, p₃, p₄)
end
endPumasModel
Parameters: b₁, b₂, b₃, slope
Random effects:
Covariates: conc
Dynamical system variables:
Dynamical system type: No dynamical model
Derived: painord
Observed: painordFinally we’ll fit our model using NaivePooled estimation method since it does not have any random-effects, i.e. no @random block:
ordinal_fit = fit(ordinal_model, pop_ord, init_params(ordinal_model), NaivePooled())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 3.103008e+03 7.031210e+02
* time: 7.390975952148438e-5
1 2.994747e+03 1.083462e+03
* time: 0.005884885787963867
2 2.406265e+03 1.884408e+02
* time: 0.010658025741577148
3 2.344175e+03 7.741610e+01
* time: 0.07176780700683594
4 2.323153e+03 2.907642e+01
* time: 0.07635378837585449
5 2.318222e+03 2.273295e+01
* time: 0.08074188232421875
6 2.316833e+03 1.390527e+01
* time: 0.08538699150085449
7 2.316425e+03 4.490883e+00
* time: 0.08992791175842285
8 2.316362e+03 9.374519e-01
* time: 0.09516787528991699
9 2.316356e+03 1.928785e-01
* time: 0.09985780715942383
10 2.316355e+03 3.119615e-02
* time: 0.10481095314025879
11 2.316355e+03 6.215513e-03
* time: 0.10923600196838379
12 2.316355e+03 8.313206e-04
* time: 0.11429691314697266FittedPumasModel
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.