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.02411794662475586
1 3.110315e+03 9.706222e+02
* time: 0.952639102935791
2 2.831659e+03 7.817006e+02
* time: 1.0428869724273682
3 2.405281e+03 2.923696e+02
* time: 1.109450101852417
4 2.370406e+03 3.032286e+02
* time: 1.1532480716705322
5 2.313631e+03 3.126188e+02
* time: 1.265902042388916
6 2.263986e+03 2.946697e+02
* time: 1.3093760013580322
7 2.160182e+03 1.917599e+02
* time: 1.3746769428253174
8 2.112467e+03 1.588288e+02
* time: 1.4323899745941162
9 2.090339e+03 1.108334e+02
* time: 1.4750690460205078
10 2.078171e+03 8.108027e+01
* time: 1.5443511009216309
11 2.074517e+03 7.813928e+01
* time: 1.579435110092163
12 2.066270e+03 7.044632e+01
* time: 1.6142029762268066
13 2.049660e+03 1.062584e+02
* time: 1.6493699550628662
14 2.021965e+03 1.130570e+02
* time: 1.6917510032653809
15 1.994936e+03 7.825801e+01
* time: 1.7338879108428955
16 1.979337e+03 5.318263e+01
* time: 1.775886058807373
17 1.972141e+03 6.807046e+01
* time: 1.8587419986724854
18 1.967973e+03 7.896361e+01
* time: 1.8970139026641846
19 1.962237e+03 8.343757e+01
* time: 1.9355950355529785
20 1.952791e+03 5.565304e+01
* time: 1.9759609699249268
21 1.935857e+03 3.923284e+01
* time: 2.0152089595794678
22 1.926254e+03 5.749643e+01
* time: 2.0536739826202393
23 1.922144e+03 4.306225e+01
* time: 2.0903079509735107
24 1.911566e+03 4.810496e+01
* time: 2.1287219524383545
25 1.906464e+03 4.324267e+01
* time: 2.166182041168213
26 1.905339e+03 1.207954e+01
* time: 2.2008299827575684
27 1.905092e+03 1.093479e+01
* time: 2.234447956085205
28 1.904957e+03 1.057034e+01
* time: 2.2679390907287598
29 1.904875e+03 1.060882e+01
* time: 2.299513101577759
30 1.904459e+03 1.031525e+01
* time: 2.36447811126709
31 1.903886e+03 1.041179e+01
* time: 2.398550033569336
32 1.903313e+03 1.135672e+01
* time: 2.432558059692383
33 1.903057e+03 1.075683e+01
* time: 2.465585947036743
34 1.902950e+03 1.091295e+01
* time: 2.4984729290008545
35 1.902887e+03 1.042409e+01
* time: 2.529660940170288
36 1.902640e+03 9.213043e+00
* time: 2.562127113342285
37 1.902364e+03 9.519321e+00
* time: 2.594928026199341
38 1.902156e+03 5.590984e+00
* time: 2.6276400089263916
39 1.902094e+03 1.811898e+00
* time: 2.6595749855041504
40 1.902086e+03 1.644722e+00
* time: 2.6909871101379395
41 1.902084e+03 1.653520e+00
* time: 2.7211129665374756
42 1.902074e+03 1.720184e+00
* time: 2.752492904663086
43 1.902055e+03 2.619061e+00
* time: 2.7838709354400635
44 1.902015e+03 3.519729e+00
* time: 2.841356039047241
45 1.901962e+03 3.403372e+00
* time: 2.8741159439086914
46 1.901924e+03 1.945644e+00
* time: 2.9066500663757324
47 1.901914e+03 6.273342e-01
* time: 2.9380569458007812
48 1.901913e+03 5.374557e-01
* time: 2.9691050052642822
49 1.901913e+03 5.710007e-01
* time: 2.9985239505767822
50 1.901913e+03 6.091390e-01
* time: 3.028083086013794
51 1.901912e+03 6.692417e-01
* time: 3.057373046875
52 1.901909e+03 7.579620e-01
* time: 3.0865209102630615
53 1.901903e+03 8.798211e-01
* time: 3.1161949634552
54 1.901889e+03 1.002981e+00
* time: 3.146406888961792
55 1.901864e+03 1.495512e+00
* time: 3.1769280433654785
56 1.901837e+03 1.664621e+00
* time: 3.207563877105713
57 1.901819e+03 8.601118e-01
* time: 3.2396790981292725
58 1.901815e+03 4.525470e-02
* time: 3.2990291118621826
59 1.901815e+03 1.294280e-02
* time: 3.330566883087158
60 1.901815e+03 2.876568e-03
* time: 3.359036922454834
61 1.901815e+03 2.876568e-03
* time: 3.439333915710449
62 1.901815e+03 2.876568e-03
* time: 3.520004987716675FittedPumasModel
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: 5.5789947509765625e-5
1 2.846050e+03 1.128657e+03
* time: 0.07049679756164551
2 2.472982e+03 7.008264e+02
* time: 0.11585378646850586
3 2.180690e+03 2.312704e+02
* time: 0.16035985946655273
4 2.125801e+03 1.862929e+02
* time: 0.20130681991577148
5 2.103173e+03 1.320946e+02
* time: 0.24557995796203613
6 2.091136e+03 1.103035e+02
* time: 0.28888583183288574
7 2.081443e+03 1.091133e+02
* time: 0.3306598663330078
8 2.071793e+03 7.197675e+01
* time: 0.37203097343444824
9 2.062706e+03 7.623310e+01
* time: 0.41521787643432617
10 2.057515e+03 6.885476e+01
* time: 0.4979557991027832
11 2.051133e+03 6.368504e+01
* time: 0.5344059467315674
12 2.038626e+03 7.730243e+01
* time: 0.5680649280548096
13 2.019352e+03 1.136864e+02
* time: 0.6045629978179932
14 1.997136e+03 1.005748e+02
* time: 0.6465268135070801
15 1.983023e+03 6.831478e+01
* time: 0.6868469715118408
16 1.977700e+03 5.649783e+01
* time: 0.7268428802490234
17 1.974583e+03 6.357642e+01
* time: 0.7651247978210449
18 1.967292e+03 7.658974e+01
* time: 0.8060839176177979
19 1.951045e+03 6.130573e+01
* time: 0.8543078899383545
20 1.935868e+03 4.845839e+01
* time: 0.8974628448486328
21 1.929356e+03 6.325782e+01
* time: 0.944774866104126
22 1.925187e+03 3.142245e+01
* time: 0.9893219470977783
23 1.923733e+03 4.623400e+01
* time: 1.091113805770874
24 1.918498e+03 5.347738e+01
* time: 1.1302859783172607
25 1.912383e+03 5.849125e+01
* time: 1.173386812210083
26 1.905510e+03 3.254038e+01
* time: 1.2157728672027588
27 1.903629e+03 2.905618e+01
* time: 1.2542328834533691
28 1.902833e+03 2.907696e+01
* time: 1.2932329177856445
29 1.902447e+03 2.746037e+01
* time: 1.328449010848999
30 1.899399e+03 1.930949e+01
* time: 1.367475986480713
31 1.898705e+03 1.186361e+01
* time: 1.405216932296753
32 1.898505e+03 1.050402e+01
* time: 1.4435269832611084
33 1.898474e+03 1.042186e+01
* time: 1.4772989749908447
34 1.897992e+03 1.238729e+01
* time: 1.514068841934204
35 1.897498e+03 1.729368e+01
* time: 1.5503950119018555
36 1.896954e+03 1.472554e+01
* time: 1.6143620014190674
37 1.896744e+03 5.852825e+00
* time: 1.6487228870391846
38 1.896705e+03 1.171353e+00
* time: 1.6798489093780518
39 1.896704e+03 1.216117e+00
* time: 1.7116169929504395
40 1.896703e+03 1.230336e+00
* time: 1.7428979873657227
41 1.896698e+03 1.250715e+00
* time: 1.7746479511260986
42 1.896688e+03 1.449552e+00
* time: 1.8065249919891357
43 1.896666e+03 2.533300e+00
* time: 1.839094877243042
44 1.896631e+03 3.075536e+00
* time: 1.8709299564361572
45 1.896599e+03 2.139797e+00
* time: 1.9042038917541504
46 1.896587e+03 6.636031e-01
* time: 1.9358198642730713
47 1.896585e+03 6.303517e-01
* time: 1.9690618515014648
48 1.896585e+03 5.995265e-01
* time: 2.0017249584198
49 1.896584e+03 5.844159e-01
* time: 2.036440849304199
50 1.896583e+03 6.083858e-01
* time: 2.099020004272461
51 1.896579e+03 8.145326e-01
* time: 2.131791830062866
52 1.896570e+03 1.293490e+00
* time: 2.1635940074920654
53 1.896549e+03 1.877705e+00
* time: 2.1954538822174072
54 1.896513e+03 2.217391e+00
* time: 2.227212905883789
55 1.896482e+03 1.658147e+00
* time: 2.2595958709716797
56 1.896466e+03 5.207215e-01
* time: 2.2927098274230957
57 1.896463e+03 1.177468e-01
* time: 2.3250069618225098
58 1.896463e+03 1.678937e-02
* time: 2.3542568683624268
59 1.896463e+03 2.666635e-03
* time: 2.3813729286193848
60 1.896463e+03 2.666635e-03
* time: 2.444844961166382FittedPumasModel
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.602836608886719e-5
1 3.572050e+03 1.302046e+03
* time: 0.09483599662780762
2 3.266947e+03 5.384244e+02
* time: 0.15069818496704102
3 3.150635e+03 1.918079e+02
* time: 0.1995849609375
4 3.108122e+03 1.277799e+02
* time: 0.2444901466369629
5 3.091358e+03 8.838080e+01
* time: 0.37560200691223145
6 3.082997e+03 8.321689e+01
* time: 0.41967201232910156
7 3.076379e+03 8.167702e+01
* time: 0.46193504333496094
8 3.067422e+03 1.177822e+02
* time: 0.5070691108703613
9 3.048580e+03 2.526969e+02
* time: 0.5581271648406982
10 2.993096e+03 6.325396e+02
* time: 0.6385340690612793
11 2.965744e+03 7.634718e+02
* time: 0.8041901588439941
12 2.921212e+03 9.704020e+02
* time: 0.9080960750579834
13 2.553649e+03 6.642510e+02
* time: 0.9877150058746338
14 2.319495e+03 3.204711e+02
* time: 1.1601121425628662
15 2.183040e+03 2.174905e+02
* time: 1.2414591312408447
16 2.157621e+03 3.150983e+02
* time: 1.2905490398406982
17 2.132395e+03 2.847614e+02
* time: 1.336780071258545
18 2.084799e+03 1.563370e+02
* time: 1.3857669830322266
19 2.071497e+03 1.006429e+02
* time: 1.4347870349884033
20 2.064983e+03 9.753313e+01
* time: 1.4840660095214844
21 2.048289e+03 9.230405e+01
* time: 1.5349860191345215
22 2.020938e+03 1.292359e+02
* time: 1.5854620933532715
23 1.983888e+03 2.284990e+02
* time: 1.6389570236206055
24 1.962132e+03 1.220188e+02
* time: 1.7285079956054688
25 1.945947e+03 1.035894e+02
* time: 1.7719080448150635
26 1.917782e+03 8.246698e+01
* time: 1.819336175918579
27 1.905967e+03 5.408054e+01
* time: 1.865328073501587
28 1.898569e+03 2.172222e+01
* time: 1.9125261306762695
29 1.897473e+03 1.689350e+01
* time: 1.956010103225708
30 1.897019e+03 1.666689e+01
* time: 2.001110076904297
31 1.896796e+03 1.699751e+01
* time: 2.048341989517212
32 1.896559e+03 1.645900e+01
* time: 2.09432315826416
33 1.896223e+03 1.415504e+01
* time: 2.140557050704956
34 1.895773e+03 1.630081e+01
* time: 2.186915159225464
35 1.895309e+03 1.723930e+01
* time: 2.235424041748047
36 1.895004e+03 1.229983e+01
* time: 2.2896201610565186
37 1.894871e+03 5.385102e+00
* time: 2.3724701404571533
38 1.894827e+03 3.465463e+00
* time: 2.414201021194458
39 1.894816e+03 3.387474e+00
* time: 2.451479196548462
40 1.894807e+03 3.295388e+00
* time: 2.4886181354522705
41 1.894786e+03 3.089194e+00
* time: 2.5283091068267822
42 1.894737e+03 2.928080e+00
* time: 2.5677080154418945
43 1.894624e+03 3.088723e+00
* time: 2.6093339920043945
44 1.894413e+03 3.493791e+00
* time: 2.653275966644287
45 1.894127e+03 3.142865e+00
* time: 2.695571184158325
46 1.893933e+03 2.145253e+00
* time: 2.7432689666748047
47 1.893888e+03 2.172800e+00
* time: 2.785860061645508
48 1.893880e+03 2.180509e+00
* time: 2.82647705078125
49 1.893876e+03 2.134257e+00
* time: 2.8667690753936768
50 1.893868e+03 2.032354e+00
* time: 2.9428551197052
51 1.893846e+03 1.760874e+00
* time: 2.9810099601745605
52 1.893796e+03 1.779016e+00
* time: 3.021088123321533
53 1.893694e+03 2.018956e+00
* time: 3.0594120025634766
54 1.893559e+03 2.366854e+00
* time: 3.0998220443725586
55 1.893474e+03 3.690142e+00
* time: 3.1418380737304688
56 1.893446e+03 3.675109e+00
* time: 3.1813600063323975
57 1.893439e+03 3.426419e+00
* time: 3.2237601280212402
58 1.893429e+03 3.183164e+00
* time: 3.2659499645233154
59 1.893398e+03 2.695171e+00
* time: 3.308475971221924
60 1.893328e+03 2.753548e+00
* time: 3.3512370586395264
61 1.893169e+03 3.589748e+00
* time: 3.398681163787842
62 1.892920e+03 3.680718e+00
* time: 3.4443411827087402
63 1.892667e+03 2.568107e+00
* time: 3.5280351638793945
64 1.892514e+03 1.087910e+00
* time: 3.5668530464172363
65 1.892493e+03 3.287296e-01
* time: 3.6040070056915283
66 1.892492e+03 2.967465e-01
* time: 3.6408300399780273
67 1.892492e+03 3.020682e-01
* time: 3.6783230304718018
68 1.892491e+03 3.034704e-01
* time: 3.714675188064575
69 1.892491e+03 3.091846e-01
* time: 3.752261161804199
70 1.892491e+03 3.224170e-01
* time: 3.7928550243377686
71 1.892490e+03 6.494197e-01
* time: 3.8346540927886963
72 1.892488e+03 1.115188e+00
* time: 3.875972032546997
73 1.892483e+03 1.838833e+00
* time: 3.920262098312378
74 1.892472e+03 2.765371e+00
* time: 3.9651811122894287
75 1.892452e+03 3.463807e+00
* time: 4.009852170944214
76 1.892431e+03 2.805270e+00
* time: 4.054888963699341
77 1.892411e+03 5.758916e-01
* time: 4.132069110870361
78 1.892410e+03 1.434041e-01
* time: 4.169824123382568
79 1.892409e+03 1.639246e-01
* time: 4.208522081375122
80 1.892409e+03 1.145856e-01
* time: 4.245826005935669
81 1.892409e+03 3.966861e-02
* time: 4.284521102905273
82 1.892409e+03 3.550808e-02
* time: 4.322102069854736
83 1.892409e+03 3.456241e-02
* time: 4.357668161392212
84 1.892409e+03 3.114018e-02
* time: 4.398128032684326
85 1.892409e+03 4.080806e-02
* time: 4.436830997467041
86 1.892409e+03 6.722726e-02
* time: 4.475823163986206
87 1.892409e+03 1.006791e-01
* time: 4.5151190757751465
88 1.892409e+03 1.303988e-01
* time: 4.556627035140991
89 1.892409e+03 1.228919e-01
* time: 4.595091104507446
90 1.892409e+03 6.433813e-02
* time: 4.6344521045684814
91 1.892409e+03 1.314164e-02
* time: 4.708472013473511
92 1.892409e+03 4.929931e-04
* time: 4.742544174194336FittedPumasModel
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.486343383789062e-5
1 3.641387e+03 1.432080e+03
* time: 0.09223389625549316
2 3.290450e+03 5.274782e+02
* time: 0.1540229320526123
3 3.185512e+03 2.173676e+02
* time: 0.21428990364074707
4 3.143009e+03 1.479653e+02
* time: 0.4192969799041748
5 3.128511e+03 8.980031e+01
* time: 0.4674489498138428
6 3.123188e+03 5.033125e+01
* time: 0.5134379863739014
7 3.120794e+03 4.279722e+01
* time: 0.5595099925994873
8 3.118627e+03 3.971051e+01
* time: 0.6061639785766602
9 3.115300e+03 8.456587e+01
* time: 0.6528770923614502
10 3.109353e+03 1.350354e+02
* time: 0.701387882232666
11 3.095894e+03 1.998258e+02
* time: 0.7615790367126465
12 2.988214e+03 4.366433e+02
* time: 0.8531849384307861
13 2.896081e+03 5.505943e+02
* time: 1.1419909000396729
14 2.652467e+03 7.300323e+02
* time: 2.3982529640197754
15 2.560937e+03 6.973661e+02
* time: 2.6945900917053223
16 2.254941e+03 2.740033e+02
* time: 2.76336407661438
17 2.222509e+03 2.034303e+02
* time: 2.8268868923187256
18 2.171255e+03 2.449580e+02
* time: 2.9419350624084473
19 2.024532e+03 1.121511e+02
* time: 2.9967410564422607
20 1.993723e+03 1.042814e+02
* time: 3.0463850498199463
21 1.985113e+03 8.079014e+01
* time: 3.093998908996582
22 1.976757e+03 7.054196e+01
* time: 3.142443895339966
23 1.969970e+03 6.070322e+01
* time: 3.189178943634033
24 1.961095e+03 6.810782e+01
* time: 3.235183000564575
25 1.947983e+03 8.116920e+01
* time: 3.2858269214630127
26 1.930371e+03 8.530051e+01
* time: 3.3434228897094727
27 1.910209e+03 6.993170e+01
* time: 3.402833938598633
28 1.899107e+03 3.362640e+01
* time: 3.462662935256958
29 1.898022e+03 2.642220e+01
* time: 3.5187509059906006
30 1.897055e+03 1.213144e+01
* time: 3.610280990600586
31 1.896596e+03 7.773239e+00
* time: 3.6574530601501465
32 1.896538e+03 7.997039e+00
* time: 3.702409029006958
33 1.896451e+03 8.160909e+00
* time: 3.7460479736328125
34 1.896283e+03 8.237721e+00
* time: 3.790138006210327
35 1.895903e+03 1.520219e+01
* time: 3.835097074508667
36 1.895272e+03 2.358916e+01
* time: 3.881589889526367
37 1.894536e+03 2.461296e+01
* time: 3.928760051727295
38 1.893995e+03 1.546128e+01
* time: 3.97455096244812
39 1.893858e+03 6.976137e+00
* time: 4.0244300365448
40 1.893833e+03 6.019466e+00
* time: 4.073564052581787
41 1.893786e+03 3.827201e+00
* time: 4.124521017074585
42 1.893714e+03 3.323412e+00
* time: 4.173871994018555
43 1.893592e+03 3.215150e+00
* time: 4.262006998062134
44 1.893435e+03 6.534965e+00
* time: 4.308404922485352
45 1.893286e+03 7.424154e+00
* time: 4.353322982788086
46 1.893190e+03 5.552627e+00
* time: 4.397260904312134
47 1.893139e+03 3.222316e+00
* time: 4.440274953842163
48 1.893120e+03 3.015339e+00
* time: 4.483469009399414
49 1.893107e+03 3.244809e+00
* time: 4.525104999542236
50 1.893080e+03 6.163100e+00
* time: 4.567095994949341
51 1.893027e+03 9.824713e+00
* time: 4.611407041549683
52 1.892912e+03 1.390100e+01
* time: 4.658844947814941
53 1.892734e+03 1.510937e+01
* time: 4.706913948059082
54 1.892561e+03 1.008563e+01
* time: 4.754807949066162
55 1.892485e+03 3.730668e+00
* time: 4.805198907852173
56 1.892471e+03 3.380261e+00
* time: 4.882783889770508
57 1.892463e+03 3.167904e+00
* time: 4.925992012023926
58 1.892441e+03 4.152065e+00
* time: 4.968106031417847
59 1.892391e+03 7.355996e+00
* time: 5.010879039764404
60 1.892268e+03 1.195397e+01
* time: 5.054686069488525
61 1.892026e+03 1.640783e+01
* time: 5.099786996841431
62 1.891735e+03 1.593576e+01
* time: 5.143990993499756
63 1.891569e+03 8.316423e+00
* time: 5.186830043792725
64 1.891494e+03 3.948212e+00
* time: 5.22967791557312
65 1.891481e+03 3.911593e+00
* time: 5.275418043136597
66 1.891457e+03 3.875559e+00
* time: 5.322772979736328
67 1.891405e+03 3.811247e+00
* time: 5.3690879344940186
68 1.891262e+03 3.657045e+00
* time: 5.450901031494141
69 1.890930e+03 4.957405e+00
* time: 5.495784044265747
70 1.890317e+03 6.657726e+00
* time: 5.540163040161133
71 1.889660e+03 6.086302e+00
* time: 5.584562063217163
72 1.889303e+03 2.270929e+00
* time: 5.628176927566528
73 1.889253e+03 7.695301e-01
* time: 5.667533874511719
74 1.889252e+03 7.382144e-01
* time: 5.706567049026489
75 1.889251e+03 7.187898e-01
* time: 5.744649887084961
76 1.889251e+03 7.215047e-01
* time: 5.78246808052063
77 1.889250e+03 7.235155e-01
* time: 5.823878049850464
78 1.889249e+03 7.246818e-01
* time: 5.869546890258789
79 1.889244e+03 7.257796e-01
* time: 5.914895057678223
80 1.889233e+03 7.198190e-01
* time: 5.959076881408691
81 1.889204e+03 1.089029e+00
* time: 6.041380882263184
82 1.889142e+03 1.801601e+00
* time: 6.084789037704468
83 1.889043e+03 2.967917e+00
* time: 6.128619909286499
84 1.888889e+03 2.965856e+00
* time: 6.169795989990234
85 1.888705e+03 5.933557e-01
* time: 6.212154865264893
86 1.888655e+03 9.577696e-01
* time: 6.253253936767578
87 1.888582e+03 1.498494e+00
* time: 6.29434609413147
88 1.888533e+03 1.502753e+00
* time: 6.3361029624938965
89 1.888490e+03 1.184664e+00
* time: 6.37846302986145
90 1.888480e+03 6.684517e-01
* time: 6.422600030899048
91 1.888476e+03 3.680034e-01
* time: 6.46841287612915
92 1.888476e+03 4.720040e-01
* time: 6.513005971908569
93 1.888476e+03 4.768646e-01
* time: 6.55768609046936
94 1.888475e+03 4.736673e-01
* time: 6.640733003616333
95 1.888475e+03 4.552764e-01
* time: 6.682302951812744
96 1.888474e+03 5.193733e-01
* time: 6.7237629890441895
97 1.888473e+03 8.850112e-01
* time: 6.763649940490723
98 1.888468e+03 1.461600e+00
* time: 6.802826881408691
99 1.888458e+03 2.209127e+00
* time: 6.841408014297485
100 1.888437e+03 2.961239e+00
* time: 6.881191968917847
101 1.888407e+03 2.978463e+00
* time: 6.922518014907837
102 1.888384e+03 1.707201e+00
* time: 6.963021993637085
103 1.888381e+03 6.199354e-01
* time: 7.0045599937438965
104 1.888380e+03 5.170909e-01
* time: 7.04862904548645
105 1.888378e+03 1.037408e-01
* time: 7.093430042266846
106 1.888378e+03 8.473247e-02
* time: 7.137336015701294
107 1.888378e+03 8.364978e-02
* time: 7.210047960281372
108 1.888378e+03 8.080446e-02
* time: 7.2502760887146
109 1.888378e+03 7.873905e-02
* time: 7.288784980773926
110 1.888378e+03 7.798398e-02
* time: 7.325792074203491
111 1.888378e+03 7.788170e-02
* time: 7.36114501953125
112 1.888378e+03 7.776461e-02
* time: 7.398515939712524
113 1.888378e+03 9.023662e-02
* time: 7.435950040817261
114 1.888378e+03 1.631390e-01
* time: 7.474936008453369
115 1.888378e+03 2.768731e-01
* time: 7.513963937759399
116 1.888377e+03 4.462386e-01
* time: 7.554758071899414
117 1.888377e+03 6.643297e-01
* time: 7.599667072296143
118 1.888375e+03 8.433397e-01
* time: 7.644263982772827
119 1.888374e+03 7.596814e-01
* time: 7.689495086669922
120 1.888373e+03 3.638119e-01
* time: 7.734737873077393
121 1.888372e+03 8.306034e-02
* time: 7.813082933425903
122 1.888372e+03 2.084513e-02
* time: 7.852784872055054
123 1.888372e+03 2.056419e-02
* time: 7.890388011932373
124 1.888372e+03 2.044080e-02
* time: 7.9268410205841064
125 1.888372e+03 2.035196e-02
* time: 7.962079048156738
126 1.888372e+03 2.021264e-02
* time: 7.997268915176392
127 1.888372e+03 1.998163e-02
* time: 8.033145904541016
128 1.888372e+03 3.161638e-02
* time: 8.070631980895996
129 1.888372e+03 5.509218e-02
* time: 8.109816074371338
130 1.888372e+03 9.275848e-02
* time: 8.14736795425415
131 1.888372e+03 1.528749e-01
* time: 8.188378095626831
132 1.888372e+03 2.461766e-01
* time: 8.230078935623169
133 1.888372e+03 3.799362e-01
* time: 8.272496938705444
134 1.888371e+03 5.311665e-01
* time: 8.346170902252197
135 1.888369e+03 6.019039e-01
* time: 8.387784957885742
136 1.888367e+03 4.664896e-01
* time: 8.428416967391968
137 1.888366e+03 1.404934e-01
* time: 8.4686119556427
138 1.888365e+03 8.513331e-02
* time: 8.50783896446228
139 1.888364e+03 1.244077e-01
* time: 8.54677700996399
140 1.888364e+03 1.028085e-01
* time: 8.58587098121643
141 1.888364e+03 5.162231e-02
* time: 8.626079082489014
142 1.888364e+03 5.149630e-02
* time: 8.666070938110352
143 1.888364e+03 3.147284e-02
* time: 8.705656051635742
144 1.888364e+03 2.104766e-02
* time: 8.747345924377441
145 1.888364e+03 6.539151e-03
* time: 8.787769079208374
146 1.888364e+03 2.540196e-03
* time: 8.828561067581177
147 1.888364e+03 4.362661e-03
* time: 8.868603944778442
148 1.888364e+03 3.034416e-03
* time: 8.942184925079346
149 1.888364e+03 5.953892e-04
* time: 8.978913068771362FittedPumasModel
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)17×2 DataFrame
| Row | Metric | Value |
|---|---|---|
| String | Any | |
| 1 | Successful | true |
| 2 | Estimation Time | 3.52 |
| 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.445 |
| 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.743 |
| 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.979 |
| 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.52 | 2.445 | 4.743 | 8.979 |
| 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 : Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}(Optim.NewtonTrustRegion{Float64}(1.0, 100.0, 1.4901161193847656e-8, 0.1, 0.25, 0.75, false), Optim.Options(x_abstol = 0.0, x_reltol = 0.0, f_abstol = 0.0, f_reltol = 0.0, g_abstol = 1.0e-5, g_reltol = 1.0e-8, outer_x_abstol = 0.0, outer_x_reltol = 0.0, outer_f_abstol = 0.0, outer_f_reltol = 0.0, outer_g_abstol = 1.0e-8, outer_g_reltol = 1.0e-8, f_calls_limit = 0, g_calls_limit = 0, h_calls_limit = 0, allow_f_increases = false, allow_outer_f_increases = true, successive_f_tol = 1, iterations = 1000, outer_iterations = 1000, store_trace = false, trace_simplex = false, show_trace = false, extended_trace = false, show_warnings = true, show_every = 1, time_limit = NaN, )
)And now we pass inspect_covariate_model_wt_crcl_sex to the goodness_of_fit plotting function:
goodness_of_fit(inspect_covariate_model_wt_crcl_sex)
The idea is that the population predictions (preds) capture the general tendency of the observations while the individual predictions (ipreds) should coincide much more closely with the observations. That is exactly what we are observing in the top row subplots in our goodness of fit plot.
Regarding the bottom row, on the left we have the weighted population residuals (wres) against time, and on the right we have the weighted individual residuals (iwres) against ipreds. Here we should not see any perceived pattern, which indicates that the error model in the model has a mean 0 and constant variance. Like before, this seems to be what we are observing in our goodness of fit plot.
Hence, our covariate model with allometric scaling and different sex creatinine clearance effectw on clearance is a good candidate for our final model.
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.004964113235473633
2 2.406265e+03 1.884408e+02
* time: 0.009858131408691406
3 2.344175e+03 7.741610e+01
* time: 0.01456904411315918
4 2.323153e+03 2.907642e+01
* time: 0.01942610740661621
5 2.318222e+03 2.273295e+01
* time: 0.02414107322692871
6 2.316833e+03 1.390527e+01
* time: 0.028913021087646484
7 2.316425e+03 4.490883e+00
* time: 0.03368711471557617
8 2.316362e+03 9.374519e-01
* time: 0.038297176361083984
9 2.316356e+03 1.928785e-01
* time: 0.04306817054748535
10 2.316355e+03 3.119615e-02
* time: 0.14829206466674805
11 2.316355e+03 6.215513e-03
* time: 0.15422296524047852
12 2.316355e+03 8.313206e-04
* time: 0.158919095993042FittedPumasModel
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.