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.
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 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.
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.
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.02519702911376953
1 3.110315e+03 9.706222e+02
* time: 1.0291941165924072
2 2.831659e+03 7.817006e+02
* time: 1.2014200687408447
3 2.405281e+03 2.923696e+02
* time: 1.2706952095031738
4 2.370406e+03 3.032286e+02
* time: 1.3183560371398926
5 2.313631e+03 3.126188e+02
* time: 1.3666820526123047
6 2.263986e+03 2.946697e+02
* time: 1.4109220504760742
7 2.160182e+03 1.917599e+02
* time: 1.4767122268676758
8 2.112467e+03 1.588288e+02
* time: 1.5379841327667236
9 2.090339e+03 1.108334e+02
* time: 1.5829980373382568
10 2.078171e+03 8.108027e+01
* time: 1.630803108215332
11 2.074517e+03 7.813928e+01
* time: 1.6732990741729736
12 2.066270e+03 7.044632e+01
* time: 1.716184139251709
13 2.049660e+03 1.062584e+02
* time: 1.8025610446929932
14 2.021965e+03 1.130570e+02
* time: 1.8425121307373047
15 1.994936e+03 7.825801e+01
* time: 1.8816630840301514
16 1.979337e+03 5.318263e+01
* time: 1.9202711582183838
17 1.972141e+03 6.807046e+01
* time: 1.9589650630950928
18 1.967973e+03 7.896361e+01
* time: 1.9974100589752197
19 1.962237e+03 8.343757e+01
* time: 2.0364010334014893
20 1.952791e+03 5.565304e+01
* time: 2.076489210128784
21 1.935857e+03 3.923284e+01
* time: 2.1177151203155518
22 1.926254e+03 5.749643e+01
* time: 2.158841133117676
23 1.922144e+03 4.306225e+01
* time: 2.196493148803711
24 1.911566e+03 4.810496e+01
* time: 2.2374460697174072
25 1.906464e+03 4.324267e+01
* time: 2.282677173614502
26 1.905339e+03 1.207954e+01
* time: 2.345985174179077
27 1.905092e+03 1.093479e+01
* time: 2.381989002227783
28 1.904957e+03 1.057034e+01
* time: 2.417689085006714
29 1.904875e+03 1.060882e+01
* time: 2.451439142227173
30 1.904459e+03 1.031525e+01
* time: 2.4874050617218018
31 1.903886e+03 1.041179e+01
* time: 2.5234901905059814
32 1.903313e+03 1.135672e+01
* time: 2.5597522258758545
33 1.903057e+03 1.075683e+01
* time: 2.5943992137908936
34 1.902950e+03 1.091295e+01
* time: 2.6303930282592773
35 1.902887e+03 1.042409e+01
* time: 2.6642260551452637
36 1.902640e+03 9.213043e+00
* time: 2.699651002883911
37 1.902364e+03 9.519321e+00
* time: 2.735473155975342
38 1.902156e+03 5.590984e+00
* time: 2.771583080291748
39 1.902094e+03 1.811898e+00
* time: 2.811385154724121
40 1.902086e+03 1.644722e+00
* time: 2.8692212104797363
41 1.902084e+03 1.653520e+00
* time: 2.902998208999634
42 1.902074e+03 1.720184e+00
* time: 2.9375991821289062
43 1.902055e+03 2.619061e+00
* time: 2.971679210662842
44 1.902015e+03 3.519729e+00
* time: 3.0056681632995605
45 1.901962e+03 3.403372e+00
* time: 3.0402510166168213
46 1.901924e+03 1.945644e+00
* time: 3.0752391815185547
47 1.901914e+03 6.273342e-01
* time: 3.1091771125793457
48 1.901913e+03 5.374557e-01
* time: 3.143541097640991
49 1.901913e+03 5.710007e-01
* time: 3.1751811504364014
50 1.901913e+03 6.091390e-01
* time: 3.2068140506744385
51 1.901912e+03 6.692417e-01
* time: 3.2397491931915283
52 1.901909e+03 7.579620e-01
* time: 3.2724180221557617
53 1.901903e+03 8.798211e-01
* time: 3.3060731887817383
54 1.901889e+03 1.002981e+00
* time: 3.3650760650634766
55 1.901864e+03 1.495512e+00
* time: 3.399322032928467
56 1.901837e+03 1.664621e+00
* time: 3.4314310550689697
57 1.901819e+03 8.601118e-01
* time: 3.4640581607818604
58 1.901815e+03 4.525470e-02
* time: 3.4967610836029053
59 1.901815e+03 1.294280e-02
* time: 3.5273921489715576
60 1.901815e+03 2.876568e-03
* time: 3.5562450885772705
61 1.901815e+03 2.876568e-03
* time: 3.6422221660614014
62 1.901815e+03 2.876568e-03
* time: 3.7290191650390625FittedPumasModel
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
---------------------4 Step 3 - Covariate Model
The third step of our covariate model building workflow is to actually develop one or more covariate models. Let’s develop three covariate models:
- allometric scaling based on weight
- clearance effect based on creatinine clearance
- separated error model based on sex
To include covariates in a Pumas model we need to first include them in the @covariates block. Then, we are free to use them inside the @pre block
Here’s our covariate model with allometric scaling based on weight:
Tip
When building covariate models, unlike in NONMEM, it is highly recommended to derive covariates or perform any required transformations or scaling as a data wrangling step and pass the derived/transformed into read_pumas as a covariate. In this particular case, it is easy to create two columns in the original data as:
@rtransform! pkdata :ASWT_CL = (:WT / 70)^0.75
@rtransform! pkdata :ASWT_Vc = (:WT / 70)Then, you bring these covariates into read_pumas and use them directly on CL and Vc. This will avoid computations of your transformations during the model fit.
covariate_model_wt = @model begin
@param begin
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvq ∈ RealDomain(; lower = 0)
tvvp ∈ RealDomain(; lower = 0)
Ω ∈ PDiagDomain(2)
σ_add ∈ RealDomain(; lower = 0)
σ_prop ∈ RealDomain(; lower = 0)
end
@random begin
η ~ MvNormal(Ω)
end
@covariates begin
WT
end
@pre begin
CL = tvcl * (WT / 70)^0.75 * exp(η[1])
Vc = tvvc * (WT / 70) * exp(η[2])
Q = tvq
Vp = tvvp
end
@dynamics Central1Periph1
@derived begin
cp := @. Central / Vc
DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
end
endPumasModel
Parameters: tvcl, tvvc, tvq, tvvp, Ω, σ_add, σ_prop
Random effects: η
Covariates: WT
Dynamical system variables: Central, Peripheral
Dynamical system type: Closed form
Derived: DV
Observed: DVOnce we included the WT covariate in the @covariates block we can use the WT values inside the @pre block. For both clearance (CL) and volume of the central compartment (Vc), we are allometric scaling by the WT value by the mean weight 70 and, in the case of CL using an allometric exponent with value 0.75.
Let’s fit our covariate_model_wt. Notice that we have not added any new parameters inside the @param block, thus, we can use the same iparams_base_model initial parameters values’ list:
fit_covariate_model_wt = fit(covariate_model_wt, pop, iparams_base_model, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 7.695401e+03 4.898919e+03
* time: 6.389617919921875e-5
1 2.846050e+03 1.128657e+03
* time: 0.07277798652648926
2 2.472982e+03 7.008264e+02
* time: 0.1244349479675293
3 2.180690e+03 2.312704e+02
* time: 0.1765739917755127
4 2.125801e+03 1.862929e+02
* time: 0.22098588943481445
5 2.103173e+03 1.320946e+02
* time: 0.32357287406921387
6 2.091136e+03 1.103035e+02
* time: 0.3651449680328369
7 2.081443e+03 1.091133e+02
* time: 0.4029698371887207
8 2.071793e+03 7.197675e+01
* time: 0.44199490547180176
9 2.062706e+03 7.623310e+01
* time: 0.48282384872436523
10 2.057515e+03 6.885476e+01
* time: 0.5209488868713379
11 2.051133e+03 6.368504e+01
* time: 0.5583789348602295
12 2.038626e+03 7.730243e+01
* time: 0.5958988666534424
13 2.019352e+03 1.136864e+02
* time: 0.6345798969268799
14 1.997136e+03 1.005748e+02
* time: 0.6733617782592773
15 1.983023e+03 6.831478e+01
* time: 0.7144818305969238
16 1.977700e+03 5.649783e+01
* time: 0.7557239532470703
17 1.974583e+03 6.357642e+01
* time: 0.7970879077911377
18 1.967292e+03 7.658974e+01
* time: 0.8828198909759521
19 1.951045e+03 6.130573e+01
* time: 0.9298779964447021
20 1.935868e+03 4.845839e+01
* time: 0.9701778888702393
21 1.929356e+03 6.325782e+01
* time: 1.0121707916259766
22 1.925187e+03 3.142245e+01
* time: 1.049914836883545
23 1.923733e+03 4.623400e+01
* time: 1.0875787734985352
24 1.918498e+03 5.347738e+01
* time: 1.1261248588562012
25 1.912383e+03 5.849125e+01
* time: 1.1693439483642578
26 1.905510e+03 3.254038e+01
* time: 1.2107958793640137
27 1.903629e+03 2.905618e+01
* time: 1.2480947971343994
28 1.902833e+03 2.907696e+01
* time: 1.2871358394622803
29 1.902447e+03 2.746037e+01
* time: 1.3243048191070557
30 1.899399e+03 1.930949e+01
* time: 1.3665308952331543
31 1.898705e+03 1.186361e+01
* time: 1.4306278228759766
32 1.898505e+03 1.050402e+01
* time: 1.4674649238586426
33 1.898474e+03 1.042186e+01
* time: 1.5010387897491455
34 1.897992e+03 1.238729e+01
* time: 1.5354578495025635
35 1.897498e+03 1.729368e+01
* time: 1.570779800415039
36 1.896954e+03 1.472554e+01
* time: 1.6064488887786865
37 1.896744e+03 5.852825e+00
* time: 1.6416618824005127
38 1.896705e+03 1.171353e+00
* time: 1.6737008094787598
39 1.896704e+03 1.216117e+00
* time: 1.7071528434753418
40 1.896703e+03 1.230336e+00
* time: 1.7390708923339844
41 1.896698e+03 1.250715e+00
* time: 1.7724409103393555
42 1.896688e+03 1.449552e+00
* time: 1.8059988021850586
43 1.896666e+03 2.533300e+00
* time: 1.8424530029296875
44 1.896631e+03 3.075536e+00
* time: 1.8815388679504395
45 1.896599e+03 2.139797e+00
* time: 1.9444868564605713
46 1.896587e+03 6.636031e-01
* time: 1.9803948402404785
47 1.896585e+03 6.303517e-01
* time: 2.013828992843628
48 1.896585e+03 5.995265e-01
* time: 2.0461058616638184
49 1.896584e+03 5.844159e-01
* time: 2.078763008117676
50 1.896583e+03 6.083858e-01
* time: 2.1110479831695557
51 1.896579e+03 8.145326e-01
* time: 2.144243001937866
52 1.896570e+03 1.293490e+00
* time: 2.1772937774658203
53 1.896549e+03 1.877705e+00
* time: 2.210641860961914
54 1.896513e+03 2.217391e+00
* time: 2.2436258792877197
55 1.896482e+03 1.658147e+00
* time: 2.2770118713378906
56 1.896466e+03 5.207215e-01
* time: 2.3108839988708496
57 1.896463e+03 1.177468e-01
* time: 2.345505952835083
58 1.896463e+03 1.678937e-02
* time: 2.3790199756622314
59 1.896463e+03 2.666635e-03
* time: 2.43227481842041
60 1.896463e+03 2.666635e-03
* time: 2.4980459213256836FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -1896.4632
Number of subjects: 30
Number of parameters: Fixed Optimized
0 8
Observation records: Active Missing
DV: 540 0
Total: 540 0
---------------------
Estimate
---------------------
tvcl 0.13915
tvvc 3.9754
tvq 0.041988
tvvp 3.5722
Ω₁,₁ 0.23874
Ω₂,₂ 0.081371
σ_add 0.032174
σ_prop 0.061012
---------------------We can definitely see that, despite not increasing the parameters of the model, our loglikelihood is higher (implying lower objective function value). Furthermore, our additive and combined error values remain almost the same while our Ωs decreased for CL and Vc. This implies that the WT covariate is definitely assisting in a better model fit by capturing the variability in CL and Vc. We’ll compare models in the next step.
Let’s now try to incorporate into this model creatinine clearance (CRCL) effect on clearance (CL):
Tip
Like the tip above, you can derive covariates or perform any required transformations or scaling as a data wrangling step and pass the derived/transformed into read_pumas as a covariate. In this particular case, it is easy to create three more columns in the original data as:
@rtransform! pkdata :CRCL_CL = :CRCL / 100
@rtransform! pkdata :ASWT_CL = (:WT / 70)^0.75
@rtransform! pkdata :ASWT_Vc = (:WT / 70)Then, you bring these covariates into read_pumas and use them directly on renCL, CL and Vc. This will avoid computations of your transformations during the model fit.
covariate_model_wt_crcl = @model begin
@param begin
tvvc ∈ RealDomain(; lower = 0)
tvq ∈ RealDomain(; lower = 0)
tvvp ∈ RealDomain(; lower = 0)
tvcl_hep ∈ RealDomain(; lower = 0)
tvcl_ren ∈ RealDomain(; lower = 0)
Ω ∈ PDiagDomain(2)
σ_add ∈ RealDomain(; lower = 0)
σ_prop ∈ RealDomain(; lower = 0)
dCRCL ∈ RealDomain()
end
@random begin
η ~ MvNormal(Ω)
end
@covariates begin
WT
CRCL
end
@pre begin
hepCL = tvcl_hep * (WT / 70)^0.75
renCL = tvcl_ren * (CRCL / 100)^dCRCL
CL = (hepCL + renCL) * exp(η[1])
Vc = tvvc * (WT / 70) * exp(η[2])
Q = tvq
Vp = tvvp
end
@dynamics Central1Periph1
@derived begin
cp := @. Central / Vc
DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
end
endPumasModel
Parameters: tvvc, tvq, tvvp, tvcl_hep, tvcl_ren, Ω, σ_add, σ_prop, dCRCL
Random effects: η
Covariates: WT, CRCL
Dynamical system variables: Central, Peripheral
Dynamical system type: Closed form
Derived: DV
Observed: DVIn the covariate_model_wt_crcl model we are keeping our allometric scaling on WT from before. But we are also adding a new covariate creatinine clearance (CRCL), dividing clearance (CL) into hepatic (hepCL) and renal clearance (renCL), along with a new parameter dCRCL.
dCRCL is the exponent of the power function for the effect of creatinine clearance on renal clearance. In some models this parameter is fixed, however we’ll allow the model to estimate it.
This is a good example on how to add covariate coefficients such as dCRCL in any Pumas covariate model. Now, let’s fit this model. Note that we need a new initial parameters values’ list since the previous one we used doesn’t include dCRCL, tvcl_hep or tvcl_ren:
iparams_covariate_model_wt_crcl = (;
tvvc = 5,
tvcl_hep = 0.01,
tvcl_ren = 0.01,
tvq = 0.01,
tvvp = 10,
Ω = Diagonal([0.01, 0.01]),
σ_add = 0.1,
σ_prop = 0.1,
dCRCL = 0.9,
)(tvvc = 5,
tvcl_hep = 0.01,
tvcl_ren = 0.01,
tvq = 0.01,
tvvp = 10,
Ω = [0.01 0.0; 0.0 0.01],
σ_add = 0.1,
σ_prop = 0.1,
dCRCL = 0.9,)fit_covariate_model_wt_crcl =
fit(covariate_model_wt_crcl, pop, iparams_covariate_model_wt_crcl, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 8.554152e+03 5.650279e+03
* time: 5.793571472167969e-5
1 3.572050e+03 1.302046e+03
* time: 0.21841692924499512
2 3.266947e+03 5.384244e+02
* time: 0.27515196800231934
3 3.150635e+03 1.918079e+02
* time: 0.324815034866333
4 3.108122e+03 1.277799e+02
* time: 0.37119293212890625
5 3.091358e+03 8.838080e+01
* time: 0.4176170825958252
6 3.082997e+03 8.321689e+01
* time: 0.46323490142822266
7 3.076379e+03 8.167702e+01
* time: 0.5127871036529541
8 3.067422e+03 1.177822e+02
* time: 0.5668690204620361
9 3.048580e+03 2.526969e+02
* time: 0.6254889965057373
10 2.993096e+03 6.325396e+02
* time: 0.7074980735778809
11 2.965744e+03 7.634718e+02
* time: 0.935020923614502
12 2.921212e+03 9.704020e+02
* time: 1.0268890857696533
13 2.553649e+03 6.642510e+02
* time: 1.0974819660186768
14 2.319495e+03 3.204711e+02
* time: 1.210387945175171
15 2.183040e+03 2.174905e+02
* time: 1.296699047088623
16 2.157621e+03 3.150983e+02
* time: 1.3495030403137207
17 2.132395e+03 2.847614e+02
* time: 1.40122389793396
18 2.084799e+03 1.563370e+02
* time: 1.4501760005950928
19 2.071497e+03 1.006429e+02
* time: 1.531275987625122
20 2.064983e+03 9.753313e+01
* time: 1.576456069946289
21 2.048289e+03 9.230405e+01
* time: 1.6211769580841064
22 2.020938e+03 1.292359e+02
* time: 1.6662020683288574
23 1.983888e+03 2.284990e+02
* time: 1.7130451202392578
24 1.962132e+03 1.220188e+02
* time: 1.7582299709320068
25 1.945947e+03 1.035894e+02
* time: 1.8013570308685303
26 1.917782e+03 8.246698e+01
* time: 1.851668119430542
27 1.905967e+03 5.408054e+01
* time: 1.9043910503387451
28 1.898569e+03 2.172222e+01
* time: 1.9569690227508545
29 1.897473e+03 1.689350e+01
* time: 2.0059900283813477
30 1.897019e+03 1.666689e+01
* time: 2.0531671047210693
31 1.896796e+03 1.699751e+01
* time: 2.103271007537842
32 1.896559e+03 1.645900e+01
* time: 2.1837799549102783
33 1.896223e+03 1.415504e+01
* time: 2.2232420444488525
34 1.895773e+03 1.630081e+01
* time: 2.263493061065674
35 1.895309e+03 1.723930e+01
* time: 2.305633068084717
36 1.895004e+03 1.229983e+01
* time: 2.344667911529541
37 1.894871e+03 5.385102e+00
* time: 2.3823230266571045
38 1.894827e+03 3.465463e+00
* time: 2.4238240718841553
39 1.894816e+03 3.387474e+00
* time: 2.4637880325317383
40 1.894807e+03 3.295388e+00
* time: 2.5031309127807617
41 1.894786e+03 3.089194e+00
* time: 2.5441110134124756
42 1.894737e+03 2.928080e+00
* time: 2.5860989093780518
43 1.894624e+03 3.088723e+00
* time: 2.6284871101379395
44 1.894413e+03 3.493791e+00
* time: 2.671765089035034
45 1.894127e+03 3.142865e+00
* time: 2.7440340518951416
46 1.893933e+03 2.145253e+00
* time: 2.782705068588257
47 1.893888e+03 2.172800e+00
* time: 2.819520950317383
48 1.893880e+03 2.180509e+00
* time: 2.8553431034088135
49 1.893876e+03 2.134257e+00
* time: 2.8918190002441406
50 1.893868e+03 2.032354e+00
* time: 2.9276700019836426
51 1.893846e+03 1.760874e+00
* time: 2.9648890495300293
52 1.893796e+03 1.779016e+00
* time: 3.005415916442871
53 1.893694e+03 2.018956e+00
* time: 3.0469419956207275
54 1.893559e+03 2.366854e+00
* time: 3.0884580612182617
55 1.893474e+03 3.690142e+00
* time: 3.131103038787842
56 1.893446e+03 3.675109e+00
* time: 3.172192096710205
57 1.893439e+03 3.426419e+00
* time: 3.2164580821990967
58 1.893429e+03 3.183164e+00
* time: 3.292294979095459
59 1.893398e+03 2.695171e+00
* time: 3.3285279273986816
60 1.893328e+03 2.753548e+00
* time: 3.3655591011047363
61 1.893169e+03 3.589748e+00
* time: 3.403580904006958
62 1.892920e+03 3.680718e+00
* time: 3.4412479400634766
63 1.892667e+03 2.568107e+00
* time: 3.478679895401001
64 1.892514e+03 1.087910e+00
* time: 3.5157461166381836
65 1.892493e+03 3.287296e-01
* time: 3.557209014892578
66 1.892492e+03 2.967465e-01
* time: 3.596877098083496
67 1.892492e+03 3.020682e-01
* time: 3.6355700492858887
68 1.892491e+03 3.034704e-01
* time: 3.672713041305542
69 1.892491e+03 3.091846e-01
* time: 3.712996006011963
70 1.892491e+03 3.224170e-01
* time: 3.753592014312744
71 1.892490e+03 6.494197e-01
* time: 3.829068899154663
72 1.892488e+03 1.115188e+00
* time: 3.865631103515625
73 1.892483e+03 1.838833e+00
* time: 3.9024550914764404
74 1.892472e+03 2.765371e+00
* time: 3.9382588863372803
75 1.892452e+03 3.463807e+00
* time: 3.974879026412964
76 1.892431e+03 2.805270e+00
* time: 4.011466026306152
77 1.892411e+03 5.758916e-01
* time: 4.048983097076416
78 1.892410e+03 1.434041e-01
* time: 4.08971905708313
79 1.892409e+03 1.639246e-01
* time: 4.132880926132202
80 1.892409e+03 1.145856e-01
* time: 4.175347089767456
81 1.892409e+03 3.966861e-02
* time: 4.2162721157073975
82 1.892409e+03 3.550808e-02
* time: 4.255378007888794
83 1.892409e+03 3.456241e-02
* time: 4.292864084243774
84 1.892409e+03 3.114018e-02
* time: 4.331734895706177
85 1.892409e+03 4.080806e-02
* time: 4.405987977981567
86 1.892409e+03 6.722726e-02
* time: 4.44135308265686
87 1.892409e+03 1.006791e-01
* time: 4.475115060806274
88 1.892409e+03 1.303988e-01
* time: 4.51045298576355
89 1.892409e+03 1.228919e-01
* time: 4.545026063919067
90 1.892409e+03 6.433813e-02
* time: 4.581047058105469
91 1.892409e+03 1.314164e-02
* time: 4.619441986083984
92 1.892409e+03 4.929931e-04
* time: 4.658466100692749FittedPumasModel
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.198883056640625e-5
1 3.641387e+03 1.432080e+03
* time: 0.09166121482849121
2 3.290450e+03 5.274782e+02
* time: 0.15293121337890625
3 3.185512e+03 2.173676e+02
* time: 0.20676517486572266
4 3.143009e+03 1.479653e+02
* time: 0.259350061416626
5 3.128511e+03 8.980031e+01
* time: 0.3171060085296631
6 3.123188e+03 5.033125e+01
* time: 0.3750782012939453
7 3.120794e+03 4.279722e+01
* time: 0.43281006813049316
8 3.118627e+03 3.971051e+01
* time: 0.49019598960876465
9 3.115300e+03 8.456587e+01
* time: 0.5485281944274902
10 3.109353e+03 1.350354e+02
* time: 0.6797389984130859
11 3.095894e+03 1.998258e+02
* time: 0.7335090637207031
12 2.988214e+03 4.366433e+02
* time: 0.8147280216217041
13 2.896081e+03 5.505943e+02
* time: 1.0103600025177002
14 2.652467e+03 7.300323e+02
* time: 2.310288190841675
15 2.560937e+03 6.973661e+02
* time: 2.634766101837158
16 2.254941e+03 2.740033e+02
* time: 2.6938750743865967
17 2.222509e+03 2.034303e+02
* time: 2.7452781200408936
18 2.171255e+03 2.449580e+02
* time: 2.7963030338287354
19 2.024532e+03 1.121511e+02
* time: 2.846710205078125
20 1.993723e+03 1.042814e+02
* time: 2.896397113800049
21 1.985113e+03 8.079014e+01
* time: 2.948674201965332
22 1.976757e+03 7.054196e+01
* time: 3.001394033432007
23 1.969970e+03 6.070322e+01
* time: 3.0525200366973877
24 1.961095e+03 6.810782e+01
* time: 3.132786989212036
25 1.947983e+03 8.116920e+01
* time: 3.179230213165283
26 1.930371e+03 8.530051e+01
* time: 3.2245709896087646
27 1.910209e+03 6.993170e+01
* time: 3.2706332206726074
28 1.899107e+03 3.362640e+01
* time: 3.31832218170166
29 1.898022e+03 2.642220e+01
* time: 3.3621280193328857
30 1.897055e+03 1.213144e+01
* time: 3.4061031341552734
31 1.896596e+03 7.773239e+00
* time: 3.4501500129699707
32 1.896538e+03 7.997039e+00
* time: 3.492859125137329
33 1.896451e+03 8.160909e+00
* time: 3.536785125732422
34 1.896283e+03 8.237721e+00
* time: 3.5858981609344482
35 1.895903e+03 1.520219e+01
* time: 3.6400492191314697
36 1.895272e+03 2.358916e+01
* time: 3.7324600219726562
37 1.894536e+03 2.461296e+01
* time: 3.7794620990753174
38 1.893995e+03 1.546128e+01
* time: 3.823970079421997
39 1.893858e+03 6.976137e+00
* time: 3.8665771484375
40 1.893833e+03 6.019466e+00
* time: 3.9074130058288574
41 1.893786e+03 3.827201e+00
* time: 3.9493181705474854
42 1.893714e+03 3.323412e+00
* time: 3.993798017501831
43 1.893592e+03 3.215150e+00
* time: 4.0382421016693115
44 1.893435e+03 6.534965e+00
* time: 4.081248044967651
45 1.893286e+03 7.424154e+00
* time: 4.1263182163238525
46 1.893190e+03 5.552627e+00
* time: 4.171922206878662
47 1.893139e+03 3.222316e+00
* time: 4.217095136642456
48 1.893120e+03 3.015339e+00
* time: 4.261669158935547
49 1.893107e+03 3.244809e+00
* time: 4.336937189102173
50 1.893080e+03 6.163100e+00
* time: 4.378269195556641
51 1.893027e+03 9.824713e+00
* time: 4.419520139694214
52 1.892912e+03 1.390100e+01
* time: 4.460760116577148
53 1.892734e+03 1.510937e+01
* time: 4.5030601024627686
54 1.892561e+03 1.008563e+01
* time: 4.545773029327393
55 1.892485e+03 3.730668e+00
* time: 4.588163137435913
56 1.892471e+03 3.380261e+00
* time: 4.628920078277588
57 1.892463e+03 3.167904e+00
* time: 4.66928505897522
58 1.892441e+03 4.152065e+00
* time: 4.717846155166626
59 1.892391e+03 7.355996e+00
* time: 4.76901912689209
60 1.892268e+03 1.195397e+01
* time: 4.821491003036499
61 1.892026e+03 1.640783e+01
* time: 4.874075174331665
62 1.891735e+03 1.593576e+01
* time: 4.958032131195068
63 1.891569e+03 8.316423e+00
* time: 5.001077175140381
64 1.891494e+03 3.948212e+00
* time: 5.041632175445557
65 1.891481e+03 3.911593e+00
* time: 5.080954074859619
66 1.891457e+03 3.875559e+00
* time: 5.12056303024292
67 1.891405e+03 3.811247e+00
* time: 5.160194158554077
68 1.891262e+03 3.657045e+00
* time: 5.200317144393921
69 1.890930e+03 4.957405e+00
* time: 5.242178201675415
70 1.890317e+03 6.657726e+00
* time: 5.285262107849121
71 1.889660e+03 6.086302e+00
* time: 5.334060192108154
72 1.889303e+03 2.270929e+00
* time: 5.38079309463501
73 1.889253e+03 7.695301e-01
* time: 5.4264280796051025
74 1.889252e+03 7.382144e-01
* time: 5.498983144760132
75 1.889251e+03 7.187898e-01
* time: 5.541349172592163
76 1.889251e+03 7.215047e-01
* time: 5.581266164779663
77 1.889250e+03 7.235155e-01
* time: 5.621249198913574
78 1.889249e+03 7.246818e-01
* time: 5.6608030796051025
79 1.889244e+03 7.257796e-01
* time: 5.701889991760254
80 1.889233e+03 7.198190e-01
* time: 5.7435760498046875
81 1.889204e+03 1.089029e+00
* time: 5.785771131515503
82 1.889142e+03 1.801601e+00
* time: 5.829092025756836
83 1.889043e+03 2.967917e+00
* time: 5.871860027313232
84 1.888889e+03 2.965856e+00
* time: 5.9192280769348145
85 1.888705e+03 5.933557e-01
* time: 5.966817140579224
86 1.888655e+03 9.577696e-01
* time: 6.012159109115601
87 1.888582e+03 1.498494e+00
* time: 6.093837022781372
88 1.888533e+03 1.502753e+00
* time: 6.136113166809082
89 1.888490e+03 1.184664e+00
* time: 6.177707195281982
90 1.888480e+03 6.684517e-01
* time: 6.217783212661743
91 1.888476e+03 3.680034e-01
* time: 6.261268138885498
92 1.888476e+03 4.720040e-01
* time: 6.300936222076416
93 1.888476e+03 4.768646e-01
* time: 6.339308023452759
94 1.888475e+03 4.736673e-01
* time: 6.377235174179077
95 1.888475e+03 4.552764e-01
* time: 6.41481614112854
96 1.888474e+03 5.193733e-01
* time: 6.453266143798828
97 1.888473e+03 8.850112e-01
* time: 6.497563123703003
98 1.888468e+03 1.461600e+00
* time: 6.541508197784424
99 1.888458e+03 2.209127e+00
* time: 6.585536003112793
100 1.888437e+03 2.961239e+00
* time: 6.668373107910156
101 1.888407e+03 2.978463e+00
* time: 6.7123730182647705
102 1.888384e+03 1.707201e+00
* time: 6.754968166351318
103 1.888381e+03 6.199354e-01
* time: 6.7993152141571045
104 1.888380e+03 5.170909e-01
* time: 6.841914176940918
105 1.888378e+03 1.037408e-01
* time: 6.884415149688721
106 1.888378e+03 8.473247e-02
* time: 6.923729181289673
107 1.888378e+03 8.364978e-02
* time: 6.963803052902222
108 1.888378e+03 8.080446e-02
* time: 7.003874063491821
109 1.888378e+03 7.873905e-02
* time: 7.04256010055542
110 1.888378e+03 7.798398e-02
* time: 7.085312128067017
111 1.888378e+03 7.788170e-02
* time: 7.12665319442749
112 1.888378e+03 7.776461e-02
* time: 7.168719053268433
113 1.888378e+03 9.023662e-02
* time: 7.248385190963745
114 1.888378e+03 1.631390e-01
* time: 7.289265155792236
115 1.888378e+03 2.768731e-01
* time: 7.33092999458313
116 1.888377e+03 4.462386e-01
* time: 7.3718321323394775
117 1.888377e+03 6.643297e-01
* time: 7.41212010383606
118 1.888375e+03 8.433397e-01
* time: 7.451935052871704
119 1.888374e+03 7.596814e-01
* time: 7.4936230182647705
120 1.888373e+03 3.638119e-01
* time: 7.534501075744629
121 1.888372e+03 8.306034e-02
* time: 7.5742270946502686
122 1.888372e+03 2.084513e-02
* time: 7.6138200759887695
123 1.888372e+03 2.056419e-02
* time: 7.655054092407227
124 1.888372e+03 2.044080e-02
* time: 7.695799112319946
125 1.888372e+03 2.035196e-02
* time: 7.7358880043029785
126 1.888372e+03 2.021264e-02
* time: 7.776039123535156
127 1.888372e+03 1.998163e-02
* time: 7.845039129257202
128 1.888372e+03 3.161638e-02
* time: 7.884169101715088
129 1.888372e+03 5.509218e-02
* time: 7.922773122787476
130 1.888372e+03 9.275848e-02
* time: 7.9626100063323975
131 1.888372e+03 1.528749e-01
* time: 8.001838207244873
132 1.888372e+03 2.461766e-01
* time: 8.041743040084839
133 1.888372e+03 3.799362e-01
* time: 8.082329034805298
134 1.888371e+03 5.311665e-01
* time: 8.122904062271118
135 1.888369e+03 6.019039e-01
* time: 8.163978099822998
136 1.888367e+03 4.664896e-01
* time: 8.205352067947388
137 1.888366e+03 1.404934e-01
* time: 8.249914169311523
138 1.888365e+03 8.513331e-02
* time: 8.29580807685852
139 1.888364e+03 1.244077e-01
* time: 8.341844081878662
140 1.888364e+03 1.028085e-01
* time: 8.422417163848877
141 1.888364e+03 5.162231e-02
* time: 8.464368104934692
142 1.888364e+03 5.149630e-02
* time: 8.504476070404053
143 1.888364e+03 3.147284e-02
* time: 8.544013023376465
144 1.888364e+03 2.104766e-02
* time: 8.58292818069458
145 1.888364e+03 6.539151e-03
* time: 8.620529174804688
146 1.888364e+03 2.540196e-03
* time: 8.657837152481079
147 1.888364e+03 4.362661e-03
* time: 8.695493221282959
148 1.888364e+03 3.034416e-03
* time: 8.733221054077148
149 1.888364e+03 5.953892e-04
* time: 8.76960301399231FittedPumasModel
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.
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)WARNING: using deprecated binding Distributions.MatrixReshaped in Pumas.
, use Distributions.ReshapedDistribution{2, S, D} where D<:Distributions.Distribution{Distributions.ArrayLikeVariate{1}, S} where S<:Distributions.ValueSupport instead.17×2 DataFrame
| Row | Metric | Value |
|---|---|---|
| String | Any | |
| 1 | Successful | true |
| 2 | Estimation Time | 3.729 |
| 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.498 |
| 3 | Subjects | 30 |
| 4 | Fixed Parameters | 0 |
| 5 | Optimized Parameters | 8 |
| 6 | DV Active Observations | 540 |
| 7 | DV Missing Observations | 0 |
| 8 | Total Active Observations | 540 |
| 9 | Total Missing Observations | 0 |
| 10 | Likelihood Approximation | Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}} |
| 11 | LogLikelihood (LL) | -1896.46 |
| 12 | -2LL | 3792.93 |
| 13 | AIC | 3808.93 |
| 14 | BIC | 3843.26 |
| 15 | (η-shrinkage) η₁ | -0.014 |
| 16 | (η-shrinkage) η₂ | -0.012 |
| 17 | (ϵ-shrinkage) DV | 0.056 |
metrics_table(fit_covariate_model_wt_crcl)17×2 DataFrame
| Row | Metric | Value |
|---|---|---|
| String | Any | |
| 1 | Successful | true |
| 2 | Estimation Time | 4.659 |
| 3 | Subjects | 30 |
| 4 | Fixed Parameters | 0 |
| 5 | Optimized Parameters | 10 |
| 6 | DV Active Observations | 540 |
| 7 | DV Missing Observations | 0 |
| 8 | Total Active Observations | 540 |
| 9 | Total Missing Observations | 0 |
| 10 | Likelihood Approximation | Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}} |
| 11 | LogLikelihood (LL) | -1892.41 |
| 12 | -2LL | 3784.82 |
| 13 | AIC | 3804.82 |
| 14 | BIC | 3847.73 |
| 15 | (η-shrinkage) η₁ | -0.014 |
| 16 | (η-shrinkage) η₂ | -0.012 |
| 17 | (ϵ-shrinkage) DV | 0.056 |
metrics_table(fit_covariate_model_wt_crcl_sex)17×2 DataFrame
| Row | Metric | Value |
|---|---|---|
| String | Any | |
| 1 | Successful | true |
| 2 | Estimation Time | 8.77 |
| 3 | Subjects | 30 |
| 4 | Fixed Parameters | 0 |
| 5 | Optimized Parameters | 13 |
| 6 | DV Active Observations | 540 |
| 7 | DV Missing Observations | 0 |
| 8 | Total Active Observations | 540 |
| 9 | Total Missing Observations | 0 |
| 10 | Likelihood Approximation | Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}} |
| 11 | LogLikelihood (LL) | -1888.36 |
| 12 | -2LL | 3776.73 |
| 13 | AIC | 3802.73 |
| 14 | BIC | 3858.52 |
| 15 | (η-shrinkage) η₁ | -0.013 |
| 16 | (η-shrinkage) η₂ | -0.012 |
| 17 | (ϵ-shrinkage) DV | 0.056 |
metrics_table outputs all of the model metrics we might be interested with respect to a certain model. That includes metadata such as estimation time, number of subjects, how many parameters were optimized and fixed, and number of observations. It also includes common model metrics like AIC, BIC, objective function value with constant (-2 loglikelihood), and so on.
We can also do an innerjoin (check our Data Wrangling Tutorials) to get all metrics into a single DataFrame:
all_metrics = innerjoin(
metrics_table(fit_base_model),
metrics_table(fit_covariate_model_wt),
metrics_table(fit_covariate_model_wt_crcl),
metrics_table(fit_covariate_model_wt_crcl_sex);
on = :Metric,
makeunique = true,
);
rename!(
all_metrics,
:Value => :Base_Model,
:Value_1 => :Covariate_Model_WT,
:Value_2 => :Covariate_Model_WT_CRCL,
:Value_3 => :Covariate_Model_WT_CRCL_SEX,
)17×5 DataFrame
| Row | Metric | Base_Model | Covariate_Model_WT | Covariate_Model_WT_CRCL | Covariate_Model_WT_CRCL_SEX |
|---|---|---|---|---|---|
| String | Any | Any | Any | Any | |
| 1 | Successful | true | true | true | true |
| 2 | Estimation Time | 3.729 | 2.498 | 4.659 | 8.77 |
| 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.
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: FOCEAnd 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.
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.
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: 6.794929504394531e-5
1 2.994747e+03 1.083462e+03
* time: 0.004851818084716797
2 2.406265e+03 1.884408e+02
* time: 0.009682893753051758
3 2.344175e+03 7.741610e+01
* time: 0.014475822448730469
4 2.323153e+03 2.907642e+01
* time: 0.019198894500732422
5 2.318222e+03 2.273295e+01
* time: 0.023793935775756836
6 2.316833e+03 1.390527e+01
* time: 0.028393030166625977
7 2.316425e+03 4.490883e+00
* time: 0.03313803672790527
8 2.316362e+03 9.374519e-01
* time: 0.03775501251220703
9 2.316356e+03 1.928785e-01
* time: 0.04253196716308594
10 2.316355e+03 3.119615e-02
* time: 0.16691899299621582
11 2.316355e+03 6.215513e-03
* time: 0.17255091667175293
12 2.316355e+03 8.313206e-04
* time: 0.17697501182556152FittedPumasModel
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.
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.
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.
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.