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.025233030319213867
1 3.110315e+03 9.706222e+02
* time: 0.5633931159973145
2 2.831659e+03 7.817006e+02
* time: 0.592919111251831
3 2.405281e+03 2.923696e+02
* time: 0.6443531513214111
4 2.370406e+03 3.032286e+02
* time: 0.6607520580291748
5 2.313631e+03 3.126188e+02
* time: 0.6776890754699707
6 2.263986e+03 2.946697e+02
* time: 0.7137610912322998
7 2.160182e+03 1.917599e+02
* time: 0.7349460124969482
8 2.112467e+03 1.588288e+02
* time: 0.7543799877166748
9 2.090339e+03 1.108334e+02
* time: 0.7697141170501709
10 2.078171e+03 8.108027e+01
* time: 0.7948000431060791
11 2.074517e+03 7.813928e+01
* time: 0.8086459636688232
12 2.066270e+03 7.044632e+01
* time: 0.8227269649505615
13 2.049660e+03 1.062584e+02
* time: 0.836712121963501
14 2.021965e+03 1.130570e+02
* time: 0.8599691390991211
15 1.994936e+03 7.825801e+01
* time: 0.8746011257171631
16 1.979337e+03 5.318263e+01
* time: 0.8891339302062988
17 1.972141e+03 6.807046e+01
* time: 0.9036159515380859
18 1.967973e+03 7.896361e+01
* time: 0.926630973815918
19 1.962237e+03 8.343757e+01
* time: 0.9409420490264893
20 1.952791e+03 5.565304e+01
* time: 0.9556450843811035
21 1.935857e+03 3.923284e+01
* time: 0.970782995223999
22 1.926254e+03 5.749643e+01
* time: 0.9940619468688965
23 1.922144e+03 4.306225e+01
* time: 1.0088369846343994
24 1.911566e+03 4.810496e+01
* time: 1.023353099822998
25 1.906464e+03 4.324267e+01
* time: 1.0381391048431396
26 1.905339e+03 1.207954e+01
* time: 1.0517759323120117
27 1.905092e+03 1.093479e+01
* time: 1.0735750198364258
28 1.904957e+03 1.057034e+01
* time: 1.0867300033569336
29 1.904875e+03 1.060882e+01
* time: 1.099985122680664
30 1.904459e+03 1.031525e+01
* time: 1.1138129234313965
31 1.903886e+03 1.041179e+01
* time: 1.1362731456756592
32 1.903313e+03 1.135672e+01
* time: 1.149960994720459
33 1.903057e+03 1.075683e+01
* time: 1.163370132446289
34 1.902950e+03 1.091295e+01
* time: 1.1768701076507568
35 1.902887e+03 1.042409e+01
* time: 1.1987149715423584
36 1.902640e+03 9.213043e+00
* time: 1.2123329639434814
37 1.902364e+03 9.519321e+00
* time: 1.2256169319152832
38 1.902156e+03 5.590984e+00
* time: 1.2392480373382568
39 1.902094e+03 1.811898e+00
* time: 1.2526030540466309
40 1.902086e+03 1.644722e+00
* time: 1.2745039463043213
41 1.902084e+03 1.653520e+00
* time: 1.2874081134796143
42 1.902074e+03 1.720184e+00
* time: 1.3006670475006104
43 1.902055e+03 2.619061e+00
* time: 1.3138909339904785
44 1.902015e+03 3.519729e+00
* time: 1.3361921310424805
45 1.901962e+03 3.403372e+00
* time: 1.3501780033111572
46 1.901924e+03 1.945644e+00
* time: 1.3634469509124756
47 1.901914e+03 6.273342e-01
* time: 1.3768539428710938
48 1.901913e+03 5.374557e-01
* time: 1.389941930770874
49 1.901913e+03 5.710007e-01
* time: 1.4112370014190674
50 1.901913e+03 6.091390e-01
* time: 1.423902988433838
51 1.901912e+03 6.692417e-01
* time: 1.436927080154419
52 1.901909e+03 7.579620e-01
* time: 1.4499571323394775
53 1.901903e+03 8.798211e-01
* time: 1.4633090496063232
54 1.901889e+03 1.002981e+00
* time: 1.4852571487426758
55 1.901864e+03 1.495512e+00
* time: 1.4990370273590088
56 1.901837e+03 1.664621e+00
* time: 1.5119471549987793
57 1.901819e+03 8.601119e-01
* time: 1.5250351428985596
58 1.901815e+03 4.525470e-02
* time: 1.5471720695495605
59 1.901815e+03 1.294280e-02
* time: 1.5595319271087646
60 1.901815e+03 2.876567e-03
* time: 1.5715320110321045
61 1.901815e+03 2.876567e-03
* time: 1.5984001159667969
62 1.901815e+03 2.876567e-03
* time: 1.6265349388122559FittedPumasModel
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: 8.893013000488281e-5
1 2.846050e+03 1.128657e+03
* time: 0.025747060775756836
2 2.472982e+03 7.008264e+02
* time: 0.044059038162231445
3 2.180690e+03 2.312704e+02
* time: 0.08494710922241211
4 2.125801e+03 1.862929e+02
* time: 0.10081696510314941
5 2.103173e+03 1.320946e+02
* time: 0.11586594581604004
6 2.091136e+03 1.103035e+02
* time: 0.1310410499572754
7 2.081443e+03 1.091133e+02
* time: 0.1656498908996582
8 2.071793e+03 7.197675e+01
* time: 0.18050503730773926
9 2.062706e+03 7.623310e+01
* time: 0.19486188888549805
10 2.057515e+03 6.885476e+01
* time: 0.2091810703277588
11 2.051133e+03 6.368504e+01
* time: 0.24118804931640625
12 2.038626e+03 7.730243e+01
* time: 0.2556910514831543
13 2.019352e+03 1.136864e+02
* time: 0.26987600326538086
14 1.997136e+03 1.005748e+02
* time: 0.28415393829345703
15 1.983023e+03 6.831478e+01
* time: 0.30750298500061035
16 1.977700e+03 5.649783e+01
* time: 0.32238101959228516
17 1.974583e+03 6.357642e+01
* time: 0.336406946182251
18 1.967292e+03 7.658974e+01
* time: 0.3510620594024658
19 1.951045e+03 6.130573e+01
* time: 0.3674759864807129
20 1.935868e+03 4.845839e+01
* time: 0.391002893447876
21 1.929356e+03 6.325782e+01
* time: 0.4063599109649658
22 1.925187e+03 3.142245e+01
* time: 0.4205820560455322
23 1.923733e+03 4.623400e+01
* time: 0.4348750114440918
24 1.918498e+03 5.347738e+01
* time: 0.4583580493927002
25 1.912383e+03 5.849125e+01
* time: 0.47402501106262207
26 1.905510e+03 3.254038e+01
* time: 0.4892909526824951
27 1.903629e+03 2.905618e+01
* time: 0.5034539699554443
28 1.902833e+03 2.907696e+01
* time: 0.5262720584869385
29 1.902447e+03 2.746037e+01
* time: 0.5397360324859619
30 1.899399e+03 1.930949e+01
* time: 0.5538539886474609
31 1.898705e+03 1.186361e+01
* time: 0.567889928817749
32 1.898505e+03 1.050402e+01
* time: 0.5910990238189697
33 1.898474e+03 1.042186e+01
* time: 0.6043200492858887
34 1.897992e+03 1.238729e+01
* time: 0.6176431179046631
35 1.897498e+03 1.729368e+01
* time: 0.6310689449310303
36 1.896954e+03 1.472554e+01
* time: 0.6445930004119873
37 1.896744e+03 5.852825e+00
* time: 0.6665370464324951
38 1.896705e+03 1.171353e+00
* time: 0.6791660785675049
39 1.896704e+03 1.216117e+00
* time: 0.6919090747833252
40 1.896703e+03 1.230336e+00
* time: 0.7047131061553955
41 1.896698e+03 1.250715e+00
* time: 0.7261760234832764
42 1.896688e+03 1.449552e+00
* time: 0.7393569946289062
43 1.896666e+03 2.533300e+00
* time: 0.7522430419921875
44 1.896631e+03 3.075537e+00
* time: 0.7651491165161133
45 1.896599e+03 2.139797e+00
* time: 0.7781789302825928
46 1.896587e+03 6.636030e-01
* time: 0.8000369071960449
47 1.896585e+03 6.303517e-01
* time: 0.813201904296875
48 1.896585e+03 5.995265e-01
* time: 0.8258650302886963
49 1.896584e+03 5.844159e-01
* time: 0.8385980129241943
50 1.896583e+03 6.083858e-01
* time: 0.8606688976287842
51 1.896579e+03 8.145327e-01
* time: 0.874345064163208
52 1.896570e+03 1.293490e+00
* time: 0.8874599933624268
53 1.896549e+03 1.877705e+00
* time: 0.9004569053649902
54 1.896513e+03 2.217392e+00
* time: 0.9135398864746094
55 1.896482e+03 1.658148e+00
* time: 0.9356570243835449
56 1.896466e+03 5.207218e-01
* time: 0.9493510723114014
57 1.896463e+03 1.177468e-01
* time: 0.9623169898986816
58 1.896463e+03 1.678937e-02
* time: 0.9745950698852539
59 1.896463e+03 2.666636e-03
* time: 0.9865350723266602
60 1.896463e+03 2.666636e-03
* time: 1.0213100910186768
61 1.896463e+03 2.666636e-03
* time: 1.0457019805908203FittedPumasModel
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.081031799316406e-5
1 3.572050e+03 1.302046e+03
* time: 0.05837583541870117
2 3.266947e+03 5.384244e+02
* time: 0.07997393608093262
3 3.150635e+03 1.918079e+02
* time: 0.09864187240600586
4 3.108122e+03 1.277799e+02
* time: 0.11684703826904297
5 3.091358e+03 8.838080e+01
* time: 0.1560819149017334
6 3.082997e+03 8.321689e+01
* time: 0.17387604713439941
7 3.076379e+03 8.167702e+01
* time: 0.1948988437652588
8 3.067422e+03 1.177822e+02
* time: 0.2595999240875244
9 3.048580e+03 2.526969e+02
* time: 0.27942800521850586
10 2.993096e+03 6.325396e+02
* time: 0.3072960376739502
11 2.965744e+03 7.634718e+02
* time: 0.37059688568115234
12 2.921212e+03 9.704020e+02
* time: 0.41510581970214844
13 2.553649e+03 6.642510e+02
* time: 0.4413149356842041
14 2.319495e+03 3.204711e+02
* time: 0.49010181427001953
15 2.183040e+03 2.174905e+02
* time: 0.5197699069976807
16 2.157621e+03 3.150983e+02
* time: 0.5373950004577637
17 2.132395e+03 2.847614e+02
* time: 0.5684218406677246
18 2.084799e+03 1.563370e+02
* time: 0.5853660106658936
19 2.071497e+03 1.006429e+02
* time: 0.6019878387451172
20 2.064983e+03 9.753313e+01
* time: 0.6284220218658447
21 2.048289e+03 9.230405e+01
* time: 0.6460068225860596
22 2.020938e+03 1.292359e+02
* time: 0.6630580425262451
23 1.983888e+03 2.284990e+02
* time: 0.6896569728851318
24 1.962132e+03 1.220188e+02
* time: 0.7067289352416992
25 1.945947e+03 1.035894e+02
* time: 0.7228419780731201
26 1.917782e+03 8.246698e+01
* time: 0.7394688129425049
27 1.905967e+03 5.408054e+01
* time: 0.7659568786621094
28 1.898569e+03 2.172222e+01
* time: 0.7830579280853271
29 1.897473e+03 1.689350e+01
* time: 0.7992818355560303
30 1.897019e+03 1.666689e+01
* time: 0.8245558738708496
31 1.896796e+03 1.699751e+01
* time: 0.840648889541626
32 1.896559e+03 1.645900e+01
* time: 0.8564648628234863
33 1.896223e+03 1.415504e+01
* time: 0.8721208572387695
34 1.895773e+03 1.630081e+01
* time: 0.8972899913787842
35 1.895309e+03 1.723930e+01
* time: 0.9136209487915039
36 1.895004e+03 1.229983e+01
* time: 0.9300539493560791
37 1.894871e+03 5.385102e+00
* time: 0.9548499584197998
38 1.894827e+03 3.465463e+00
* time: 0.970836877822876
39 1.894816e+03 3.387474e+00
* time: 0.9860689640045166
40 1.894807e+03 3.295388e+00
* time: 1.0014638900756836
41 1.894786e+03 3.089194e+00
* time: 1.0260870456695557
42 1.894737e+03 2.928080e+00
* time: 1.0418708324432373
43 1.894624e+03 3.088723e+00
* time: 1.0586419105529785
44 1.894413e+03 3.493791e+00
* time: 1.083225965499878
45 1.894127e+03 3.142865e+00
* time: 1.0993478298187256
46 1.893933e+03 2.145253e+00
* time: 1.1153810024261475
47 1.893888e+03 2.172800e+00
* time: 1.131152868270874
48 1.893880e+03 2.180509e+00
* time: 1.1560368537902832
49 1.893876e+03 2.134257e+00
* time: 1.1715550422668457
50 1.893868e+03 2.032354e+00
* time: 1.1869208812713623
51 1.893846e+03 1.760874e+00
* time: 1.2025270462036133
52 1.893796e+03 1.779016e+00
* time: 1.2274470329284668
53 1.893694e+03 2.018956e+00
* time: 1.2430570125579834
54 1.893559e+03 2.366854e+00
* time: 1.2588679790496826
55 1.893474e+03 3.690142e+00
* time: 1.2835140228271484
56 1.893446e+03 3.675109e+00
* time: 1.2992708683013916
57 1.893439e+03 3.426419e+00
* time: 1.3145458698272705
58 1.893429e+03 3.183164e+00
* time: 1.3296568393707275
59 1.893398e+03 2.695171e+00
* time: 1.354038953781128
60 1.893328e+03 2.753548e+00
* time: 1.3702998161315918
61 1.893169e+03 3.589748e+00
* time: 1.386132001876831
62 1.892920e+03 3.680718e+00
* time: 1.4110639095306396
63 1.892667e+03 2.568107e+00
* time: 1.427513837814331
64 1.892514e+03 1.087910e+00
* time: 1.4436569213867188
65 1.892493e+03 3.287296e-01
* time: 1.459183931350708
66 1.892492e+03 2.967465e-01
* time: 1.4836020469665527
67 1.892492e+03 3.020682e-01
* time: 1.498896837234497
68 1.892491e+03 3.034704e-01
* time: 1.5133769512176514
69 1.892491e+03 3.091846e-01
* time: 1.528303861618042
70 1.892491e+03 3.224170e-01
* time: 1.552422046661377
71 1.892490e+03 6.494197e-01
* time: 1.571465015411377
72 1.892488e+03 1.115188e+00
* time: 1.5872108936309814
73 1.892483e+03 1.838833e+00
* time: 1.6115539073944092
74 1.892472e+03 2.765371e+00
* time: 1.6278200149536133
75 1.892452e+03 3.463807e+00
* time: 1.6435439586639404
76 1.892431e+03 2.805270e+00
* time: 1.6590609550476074
77 1.892411e+03 5.758916e-01
* time: 1.6837358474731445
78 1.892410e+03 1.434041e-01
* time: 1.6991708278656006
79 1.892409e+03 1.639246e-01
* time: 1.7142529487609863
80 1.892409e+03 1.145856e-01
* time: 1.7291069030761719
81 1.892409e+03 3.966861e-02
* time: 1.7529828548431396
82 1.892409e+03 3.550808e-02
* time: 1.767482042312622
83 1.892409e+03 3.456241e-02
* time: 1.7817158699035645
84 1.892409e+03 3.114018e-02
* time: 1.804823875427246
85 1.892409e+03 4.080806e-02
* time: 1.8200600147247314
86 1.892409e+03 6.722726e-02
* time: 1.8348100185394287
87 1.892409e+03 1.006791e-01
* time: 1.8495209217071533
88 1.892409e+03 1.303988e-01
* time: 1.8735628128051758
89 1.892409e+03 1.228919e-01
* time: 1.8887019157409668
90 1.892409e+03 6.433813e-02
* time: 1.9034550189971924
91 1.892409e+03 1.314164e-02
* time: 1.9183900356292725
92 1.892409e+03 4.929931e-04
* time: 1.9418559074401855FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -1892.409
Number of subjects: 30
Number of parameters: Fixed Optimized
0 10
Observation records: Active Missing
DV: 540 0
Total: 540 0
-----------------------
Estimate
-----------------------
tvvc 3.9757
tvq 0.042177
tvvp 3.6434
tvcl_hep 0.058572
tvcl_ren 0.1337
Ω₁,₁ 0.18299
Ω₂,₂ 0.081353
σ_add 0.032174
σ_prop 0.06101
dCRCL 1.1331
-----------------------As before, our loglikelihood is higher (implying lower objective function value). Furthermore, our additive and combined error values remain almost the same while our Ω on CL, Ω₁,₁, decreased. This implies that the CRCL covariate with an estimated exponent dCRCL is definitely assisting in a better model fit.
Finally let’s include a separated CL model based on sex as a third covariate model:
covariate_model_wt_crcl_sex = @model begin
@param begin
tvvc ∈ RealDomain(; lower = 0)
tvq ∈ RealDomain(; lower = 0)
tvvp ∈ RealDomain(; lower = 0)
tvcl_hep_M ∈ RealDomain(; lower = 0)
tvcl_hep_F ∈ RealDomain(; lower = 0)
tvcl_ren_M ∈ RealDomain(; lower = 0)
tvcl_ren_F ∈ RealDomain(; lower = 0)
Ω ∈ PDiagDomain(2)
σ_add ∈ RealDomain(; lower = 0)
σ_prop ∈ RealDomain(; lower = 0)
dCRCL_M ∈ RealDomain()
dCRCL_F ∈ RealDomain()
end
@random begin
η ~ MvNormal(Ω)
end
@covariates begin
WT
CRCL
SEX
end
@pre begin
tvcl_hep = ifelse(SEX == "M", tvcl_hep_M, tvcl_hep_F)
tvcl_ren = ifelse(SEX == "M", tvcl_ren_M, tvcl_ren_F)
dCRCL = ifelse(SEX == "M", dCRCL_M, dCRCL_F)
hepCL = tvcl_hep * (WT / 70)^0.75
renCL = tvcl_ren * (CRCL / 100)^dCRCL
CL = (hepCL + renCL) * exp(η[1])
Vc = tvvc * (WT / 70) * exp(η[2])
Q = tvq
Vp = tvvp
end
@dynamics Central1Periph1
@derived begin
cp := @. Central / Vc
DV ~ @. Normal(cp, sqrt(cp^2 * σ_prop^2 + σ_add^2))
end
endPumasModel
Parameters: tvvc, tvq, tvvp, tvcl_hep_M, tvcl_hep_F, tvcl_ren_M, tvcl_ren_F, Ω, σ_add, σ_prop, dCRCL_M, dCRCL_F
Random effects: η
Covariates: WT, CRCL, SEX
Dynamical system variables: Central, Peripheral
Dynamical system type: Closed form
Derived: DV
Observed: DVIn the covariate_model_wt_crcl_sex model we are keeping our allometric scaling on WT, the CRCL effect on renCL, and the breakdown of CL into hepCL and renCL from before. However we are separating them with different values by sex. Hence, we have a new covariate SEX and six new parameters in the @param block by expanding tvcl_hep, tvcl_ren, and dCRCL into male (suffix M) and female (suffix F).
This is a good example on how to add create binary values based on covariate values such as SEX inside the @pre block with the ifelse function. Now, let’s fit this model. Note that we need a new initial parameters values’ list since the previous one we used had a single tvcl_hep, tvcl_ren, and dCRCL:
iparams_covariate_model_wt_crcl_sex = (;
tvvc = 5,
tvcl_hep_M = 0.01,
tvcl_hep_F = 0.01,
tvcl_ren_M = 0.01,
tvcl_ren_F = 0.01,
tvq = 0.01,
tvvp = 10,
Ω = Diagonal([0.01, 0.01]),
σ_add = 0.1,
σ_prop = 0.1,
dCRCL_M = 0.9,
dCRCL_F = 0.9,
)(tvvc = 5,
tvcl_hep_M = 0.01,
tvcl_hep_F = 0.01,
tvcl_ren_M = 0.01,
tvcl_ren_F = 0.01,
tvq = 0.01,
tvvp = 10,
Ω = [0.01 0.0; 0.0 0.01],
σ_add = 0.1,
σ_prop = 0.1,
dCRCL_M = 0.9,
dCRCL_F = 0.9,)fit_covariate_model_wt_crcl_sex =
fit(covariate_model_wt_crcl_sex, pop, iparams_covariate_model_wt_crcl_sex, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 8.554152e+03 5.650279e+03
* time: 7.295608520507812e-5
1 3.641387e+03 1.432080e+03
* time: 0.03085184097290039
2 3.290450e+03 5.274782e+02
* time: 0.07684087753295898
3 3.185512e+03 2.173676e+02
* time: 0.09740591049194336
4 3.143009e+03 1.479653e+02
* time: 0.11724495887756348
5 3.128511e+03 8.980031e+01
* time: 0.1582047939300537
6 3.123188e+03 5.033125e+01
* time: 0.1774768829345703
7 3.120794e+03 4.279722e+01
* time: 0.20621085166931152
8 3.118627e+03 3.971051e+01
* time: 0.22561097145080566
9 3.115300e+03 8.456587e+01
* time: 0.24505400657653809
10 3.109353e+03 1.350354e+02
* time: 0.27469778060913086
11 3.095894e+03 1.998258e+02
* time: 0.29598283767700195
12 2.988214e+03 4.366433e+02
* time: 0.3340299129486084
13 2.896081e+03 5.505943e+02
* time: 0.4085547924041748
14 2.652467e+03 7.300323e+02
* time: 0.8503079414367676
15 2.560937e+03 6.973661e+02
* time: 0.9559969902038574
16 2.254941e+03 2.740033e+02
* time: 0.9795048236846924
17 2.222509e+03 2.034303e+02
* time: 1.0100939273834229
18 2.171255e+03 2.449580e+02
* time: 1.0313007831573486
19 2.024532e+03 1.121511e+02
* time: 1.0613298416137695
20 1.993723e+03 1.042814e+02
* time: 1.0815958976745605
21 1.985113e+03 8.079014e+01
* time: 1.1004078388214111
22 1.976757e+03 7.054196e+01
* time: 1.1284708976745605
23 1.969970e+03 6.070322e+01
* time: 1.1475379467010498
24 1.961095e+03 6.810782e+01
* time: 1.16597580909729
25 1.947983e+03 8.116920e+01
* time: 1.1945569515228271
26 1.930371e+03 8.530051e+01
* time: 1.213709831237793
27 1.910209e+03 6.993170e+01
* time: 1.2425739765167236
28 1.899107e+03 3.362640e+01
* time: 1.262436866760254
29 1.898022e+03 2.642220e+01
* time: 1.280707836151123
30 1.897055e+03 1.213144e+01
* time: 1.3081879615783691
31 1.896596e+03 7.773239e+00
* time: 1.3272387981414795
32 1.896538e+03 7.997039e+00
* time: 1.345228910446167
33 1.896451e+03 8.160909e+00
* time: 1.3728687763214111
34 1.896283e+03 8.237721e+00
* time: 1.391319990158081
35 1.895903e+03 1.520219e+01
* time: 1.4093759059906006
36 1.895272e+03 2.358916e+01
* time: 1.437647819519043
37 1.894536e+03 2.461296e+01
* time: 1.4570119380950928
38 1.893995e+03 1.546128e+01
* time: 1.4855828285217285
39 1.893858e+03 6.976137e+00
* time: 1.5043158531188965
40 1.893833e+03 6.019466e+00
* time: 1.5223608016967773
41 1.893786e+03 3.827201e+00
* time: 1.5500338077545166
42 1.893714e+03 3.323412e+00
* time: 1.5691449642181396
43 1.893592e+03 3.215150e+00
* time: 1.5886638164520264
44 1.893435e+03 6.534965e+00
* time: 1.6169419288635254
45 1.893286e+03 7.424154e+00
* time: 1.6363658905029297
46 1.893190e+03 5.552627e+00
* time: 1.6549489498138428
47 1.893139e+03 3.222316e+00
* time: 1.683459997177124
48 1.893120e+03 3.015339e+00
* time: 1.7016489505767822
49 1.893107e+03 3.244809e+00
* time: 1.7194859981536865
50 1.893080e+03 6.163100e+00
* time: 1.7465438842773438
51 1.893027e+03 9.824713e+00
* time: 1.7645368576049805
52 1.892912e+03 1.390100e+01
* time: 1.7921888828277588
53 1.892734e+03 1.510937e+01
* time: 1.8107588291168213
54 1.892561e+03 1.008563e+01
* time: 1.8285629749298096
55 1.892485e+03 3.730668e+00
* time: 1.8560519218444824
56 1.892471e+03 3.380261e+00
* time: 1.8743538856506348
57 1.892463e+03 3.167904e+00
* time: 1.8919718265533447
58 1.892441e+03 4.152065e+00
* time: 1.919480800628662
59 1.892391e+03 7.355996e+00
* time: 1.9374818801879883
60 1.892268e+03 1.195397e+01
* time: 1.955068826675415
61 1.892026e+03 1.640783e+01
* time: 1.9830129146575928
62 1.891735e+03 1.593576e+01
* time: 2.00213885307312
63 1.891569e+03 8.316423e+00
* time: 2.0302469730377197
64 1.891494e+03 3.948212e+00
* time: 2.0486087799072266
65 1.891481e+03 3.911593e+00
* time: 2.0666329860687256
66 1.891457e+03 3.875559e+00
* time: 2.094381809234619
67 1.891405e+03 3.811247e+00
* time: 2.112484931945801
68 1.891262e+03 3.657045e+00
* time: 2.129974842071533
69 1.890930e+03 4.957405e+00
* time: 2.157970905303955
70 1.890317e+03 6.657726e+00
* time: 2.177428960800171
71 1.889660e+03 6.086302e+00
* time: 2.196016788482666
72 1.889303e+03 2.270929e+00
* time: 2.2243077754974365
73 1.889253e+03 7.695301e-01
* time: 2.243053913116455
74 1.889252e+03 7.382144e-01
* time: 2.261218786239624
75 1.889251e+03 7.187898e-01
* time: 2.289094924926758
76 1.889251e+03 7.215047e-01
* time: 2.307013988494873
77 1.889250e+03 7.235155e-01
* time: 2.3255558013916016
78 1.889249e+03 7.246818e-01
* time: 2.3531899452209473
79 1.889244e+03 7.257796e-01
* time: 2.3713319301605225
80 1.889233e+03 7.198190e-01
* time: 2.3993239402770996
81 1.889204e+03 1.089029e+00
* time: 2.4178009033203125
82 1.889142e+03 1.801601e+00
* time: 2.4356319904327393
83 1.889043e+03 2.967917e+00
* time: 2.462601900100708
84 1.888889e+03 2.965856e+00
* time: 2.481362819671631
85 1.888705e+03 5.933554e-01
* time: 2.4997448921203613
86 1.888655e+03 9.577699e-01
* time: 2.527549982070923
87 1.888582e+03 1.498494e+00
* time: 2.5458779335021973
88 1.888533e+03 1.502750e+00
* time: 2.564161777496338
89 1.888490e+03 1.184664e+00
* time: 2.59330677986145
90 1.888480e+03 6.684513e-01
* time: 2.611515998840332
91 1.888476e+03 3.680030e-01
* time: 2.629363775253296
92 1.888476e+03 4.720039e-01
* time: 2.6566929817199707
93 1.888476e+03 4.768646e-01
* time: 2.6746628284454346
94 1.888475e+03 4.736674e-01
* time: 2.7018067836761475
95 1.888475e+03 4.552766e-01
* time: 2.7193689346313477
96 1.888474e+03 5.193719e-01
* time: 2.7368509769439697
97 1.888473e+03 8.850088e-01
* time: 2.763478994369507
98 1.888468e+03 1.461597e+00
* time: 2.781486988067627
99 1.888458e+03 2.209123e+00
* time: 2.7990989685058594
100 1.888437e+03 2.961234e+00
* time: 2.826258897781372
101 1.888407e+03 2.978462e+00
* time: 2.8448219299316406
102 1.888384e+03 1.707197e+00
* time: 2.8627548217773438
103 1.888381e+03 6.198730e-01
* time: 2.8900928497314453
104 1.888380e+03 5.171201e-01
* time: 2.9080889225006104
105 1.888378e+03 1.037261e-01
* time: 2.9253768920898438
106 1.888378e+03 8.473257e-02
* time: 2.9518918991088867
107 1.888378e+03 8.364956e-02
* time: 2.9692459106445312
108 1.888378e+03 8.080438e-02
* time: 2.9861629009246826
109 1.888378e+03 7.873896e-02
* time: 3.0139708518981934
110 1.888378e+03 7.798398e-02
* time: 3.031312942504883
111 1.888378e+03 7.788171e-02
* time: 3.047755002975464
112 1.888378e+03 7.776461e-02
* time: 3.0744469165802
113 1.888378e+03 9.023533e-02
* time: 3.092615842819214
114 1.888378e+03 1.631356e-01
* time: 3.10988187789917
115 1.888378e+03 2.768664e-01
* time: 3.1370038986206055
116 1.888377e+03 4.462262e-01
* time: 3.154770851135254
117 1.888377e+03 6.643078e-01
* time: 3.173435926437378
118 1.888375e+03 8.433023e-01
* time: 3.200929880142212
119 1.888374e+03 7.596239e-01
* time: 3.218766927719116
120 1.888373e+03 3.637667e-01
* time: 3.2460968494415283
121 1.888372e+03 8.304667e-02
* time: 3.2636868953704834
122 1.888372e+03 2.084518e-02
* time: 3.281003952026367
123 1.888372e+03 2.056414e-02
* time: 3.3073649406433105
124 1.888372e+03 2.044078e-02
* time: 3.324521780014038
125 1.888372e+03 2.035197e-02
* time: 3.340725898742676
126 1.888372e+03 2.021268e-02
* time: 3.367154836654663
127 1.888372e+03 1.998172e-02
* time: 3.384230852127075
128 1.888372e+03 3.162406e-02
* time: 3.400624990463257
129 1.888372e+03 5.510549e-02
* time: 3.426881790161133
130 1.888372e+03 9.278088e-02
* time: 3.443943977355957
131 1.888372e+03 1.529116e-01
* time: 3.460805892944336
132 1.888372e+03 2.462349e-01
* time: 3.4872758388519287
133 1.888372e+03 3.800236e-01
* time: 3.504758834838867
134 1.888371e+03 5.312831e-01
* time: 3.5222268104553223
135 1.888369e+03 6.020265e-01
* time: 3.549909830093384
136 1.888367e+03 4.665657e-01
* time: 3.5678648948669434
137 1.888366e+03 1.404905e-01
* time: 3.585958957672119
138 1.888365e+03 8.513244e-02
* time: 3.6129918098449707
139 1.888364e+03 1.244427e-01
* time: 3.630300998687744
140 1.888364e+03 1.028331e-01
* time: 3.6475398540496826
141 1.888364e+03 5.164076e-02
* time: 3.6747729778289795
142 1.888364e+03 5.147918e-02
* time: 3.691927909851074
143 1.888364e+03 3.147222e-02
* time: 3.708711862564087
144 1.888364e+03 2.104481e-02
* time: 3.735074996948242
145 1.888364e+03 6.543267e-03
* time: 3.751828908920288
146 1.888364e+03 2.537332e-03
* time: 3.7684438228607178
147 1.888364e+03 4.361311e-03
* time: 3.794477939605713
148 1.888364e+03 3.035139e-03
* time: 3.8108878135681152
149 1.888364e+03 5.966636e-04
* time: 3.8278708457946777FittedPumasModel
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.627 |
| 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.046 |
| 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 | 1.942 |
| 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 | 3.828 |
| 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.627 | 1.046 | 1.942 | 3.828 |
| 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: 6.008148193359375e-5
1 2.994747e+03 1.083462e+03
* time: 0.0060040950775146484
2 2.406265e+03 1.884408e+02
* time: 0.07446503639221191
3 2.344175e+03 7.741610e+01
* time: 0.07858705520629883
4 2.323153e+03 2.907642e+01
* time: 0.08254194259643555
5 2.318222e+03 2.273295e+01
* time: 0.08660387992858887
6 2.316833e+03 1.390527e+01
* time: 0.09092402458190918
7 2.316425e+03 4.490883e+00
* time: 0.09606099128723145
8 2.316362e+03 9.374519e-01
* time: 0.10081911087036133
9 2.316356e+03 1.928785e-01
* time: 0.10485696792602539
10 2.316355e+03 3.119615e-02
* time: 0.10899901390075684
11 2.316355e+03 6.215513e-03
* time: 0.11344504356384277
12 2.316355e+03 8.313206e-04
* time: 0.11805391311645508FittedPumasModel
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.