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.0407719612121582 1 3.110315e+03 9.706222e+02 * time: 1.6153340339660645 2 2.831659e+03 7.817006e+02 * time: 1.651737928390503 3 2.405281e+03 2.923696e+02 * time: 1.8357789516448975 4 2.370406e+03 3.032286e+02 * time: 1.858201026916504 5 2.313631e+03 3.126188e+02 * time: 1.8816471099853516 6 2.263986e+03 2.946697e+02 * time: 1.9025020599365234 7 2.160182e+03 1.917599e+02 * time: 1.9308528900146484 8 2.112467e+03 1.588288e+02 * time: 1.9569180011749268 9 2.090339e+03 1.108334e+02 * time: 1.9762170314788818 10 2.078171e+03 8.108027e+01 * time: 1.9957959651947021 11 2.074517e+03 7.813928e+01 * time: 2.0137319564819336 12 2.066270e+03 7.044632e+01 * time: 2.0318539142608643 13 2.049660e+03 1.062584e+02 * time: 2.112204074859619 14 2.021965e+03 1.130570e+02 * time: 2.1312460899353027 15 1.994936e+03 7.825801e+01 * time: 2.1499290466308594 16 1.979337e+03 5.318263e+01 * time: 2.168313980102539 17 1.972141e+03 6.807046e+01 * time: 2.1864969730377197 18 1.967973e+03 7.896361e+01 * time: 2.204461097717285 19 1.962237e+03 8.343757e+01 * time: 2.2227120399475098 20 1.952791e+03 5.565304e+01 * time: 2.2415170669555664 21 1.935857e+03 3.923284e+01 * time: 2.27559494972229 22 1.926254e+03 5.749643e+01 * time: 2.2947399616241455 23 1.922144e+03 4.306225e+01 * time: 2.3143510818481445 24 1.911566e+03 4.810496e+01 * time: 2.333578109741211 25 1.906464e+03 4.324267e+01 * time: 2.3524200916290283 26 1.905339e+03 1.207954e+01 * time: 2.3695240020751953 27 1.905092e+03 1.093479e+01 * time: 2.3863229751586914 28 1.904957e+03 1.057034e+01 * time: 2.403101921081543 29 1.904875e+03 1.060882e+01 * time: 2.43349289894104 30 1.904459e+03 1.031525e+01 * time: 2.451401948928833 31 1.903886e+03 1.041179e+01 * time: 2.469212055206299 32 1.903313e+03 1.135672e+01 * time: 2.487020969390869 33 1.903057e+03 1.075683e+01 * time: 2.5039889812469482 34 1.902950e+03 1.091295e+01 * time: 2.5209779739379883 35 1.902887e+03 1.042409e+01 * time: 2.5503830909729004 36 1.902640e+03 9.213043e+00 * time: 2.5679678916931152 37 1.902364e+03 9.519321e+00 * time: 2.5853381156921387 38 1.902156e+03 5.590984e+00 * time: 2.6025450229644775 39 1.902094e+03 1.811898e+00 * time: 2.6192378997802734 40 1.902086e+03 1.644722e+00 * time: 2.6470460891723633 41 1.902084e+03 1.653520e+00 * time: 2.6636810302734375 42 1.902074e+03 1.720184e+00 * time: 2.680561065673828 43 1.902055e+03 2.619061e+00 * time: 2.69730806350708 44 1.902015e+03 3.519729e+00 * time: 2.7246201038360596 45 1.901962e+03 3.403372e+00 * time: 2.7416889667510986 46 1.901924e+03 1.945644e+00 * time: 2.758613109588623 47 1.901914e+03 6.273342e-01 * time: 2.7754220962524414 48 1.901913e+03 5.374557e-01 * time: 2.802412986755371 49 1.901913e+03 5.710007e-01 * time: 2.8202269077301025 50 1.901913e+03 6.091390e-01 * time: 2.8367180824279785 51 1.901912e+03 6.692417e-01 * time: 2.8533570766448975 52 1.901909e+03 7.579620e-01 * time: 2.8803489208221436 53 1.901903e+03 8.798211e-01 * time: 2.8974030017852783 54 1.901889e+03 1.002981e+00 * time: 2.9143130779266357 55 1.901864e+03 1.495512e+00 * time: 2.931148052215576 56 1.901837e+03 1.664621e+00 * time: 2.9576690196990967 57 1.901819e+03 8.601119e-01 * time: 2.974764108657837 58 1.901815e+03 4.525470e-02 * time: 2.9916510581970215 59 1.901815e+03 1.294280e-02 * time: 3.007736921310425 60 1.901815e+03 2.876567e-03 * time: 3.033768892288208 61 1.901815e+03 2.876567e-03 * time: 3.069301128387451 62 1.901815e+03 2.876567e-03 * time: 3.1114120483398438
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 30
Observation records: Active Missing
DV: 540 0
Total: 540 0
Number of parameters: Constant Optimized
0 8
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: NoObjectiveChange
Log-likelihood value: -1901.815
------------------
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: 1.71661376953125e-5 1 2.846050e+03 1.128657e+03 * time: 0.44837498664855957 2 2.472982e+03 7.008264e+02 * time: 0.4706740379333496 3 2.180690e+03 2.312704e+02 * time: 0.5692121982574463 4 2.125801e+03 1.862929e+02 * time: 0.5887610912322998 5 2.103173e+03 1.320946e+02 * time: 0.608051061630249 6 2.091136e+03 1.103035e+02 * time: 0.6265091896057129 7 2.081443e+03 1.091133e+02 * time: 0.6438171863555908 8 2.071793e+03 7.197675e+01 * time: 0.6618061065673828 9 2.062706e+03 7.623310e+01 * time: 0.6801130771636963 10 2.057515e+03 6.885476e+01 * time: 0.731827974319458 11 2.051133e+03 6.368504e+01 * time: 0.7495851516723633 12 2.038626e+03 7.730243e+01 * time: 0.7672750949859619 13 2.019352e+03 1.136864e+02 * time: 0.7850940227508545 14 1.997136e+03 1.005748e+02 * time: 0.8027920722961426 15 1.983023e+03 6.831478e+01 * time: 0.820976972579956 16 1.977700e+03 5.649783e+01 * time: 0.8505160808563232 17 1.974583e+03 6.357642e+01 * time: 0.8686120510101318 18 1.967292e+03 7.658974e+01 * time: 0.8874011039733887 19 1.951045e+03 6.130573e+01 * time: 0.9097490310668945 20 1.935868e+03 4.845839e+01 * time: 0.9287080764770508 21 1.929356e+03 6.325782e+01 * time: 0.9585230350494385 22 1.925187e+03 3.142245e+01 * time: 0.9771041870117188 23 1.923733e+03 4.623400e+01 * time: 0.9949309825897217 24 1.918498e+03 5.347738e+01 * time: 1.0229079723358154 25 1.912383e+03 5.849125e+01 * time: 1.0429420471191406 26 1.905510e+03 3.254038e+01 * time: 1.0627000331878662 27 1.903629e+03 2.905618e+01 * time: 1.0805912017822266 28 1.902833e+03 2.907696e+01 * time: 1.1081430912017822 29 1.902447e+03 2.746037e+01 * time: 1.1250250339508057 30 1.899399e+03 1.930949e+01 * time: 1.1431150436401367 31 1.898705e+03 1.186361e+01 * time: 1.160823106765747 32 1.898505e+03 1.050402e+01 * time: 1.1875929832458496 33 1.898474e+03 1.042186e+01 * time: 1.2037789821624756 34 1.897992e+03 1.238729e+01 * time: 1.2202579975128174 35 1.897498e+03 1.729368e+01 * time: 1.2373640537261963 36 1.896954e+03 1.472554e+01 * time: 1.26387619972229 37 1.896744e+03 5.852824e+00 * time: 1.2808270454406738 38 1.896705e+03 1.171353e+00 * time: 1.2964940071105957 39 1.896704e+03 1.216117e+00 * time: 1.312488079071045 40 1.896703e+03 1.230336e+00 * time: 1.3376569747924805 41 1.896698e+03 1.250715e+00 * time: 1.3541309833526611 42 1.896688e+03 1.449552e+00 * time: 1.3704571723937988 43 1.896666e+03 2.533300e+00 * time: 1.3868000507354736 44 1.896631e+03 3.075536e+00 * time: 1.4141981601715088 45 1.896599e+03 2.139797e+00 * time: 1.4308011531829834 46 1.896587e+03 6.636030e-01 * time: 1.447242021560669 47 1.896585e+03 6.303517e-01 * time: 1.4636600017547607 48 1.896585e+03 5.995265e-01 * time: 1.4888989925384521 49 1.896584e+03 5.844159e-01 * time: 1.5049970149993896 50 1.896583e+03 6.083858e-01 * time: 1.5211551189422607 51 1.896579e+03 8.145327e-01 * time: 1.5375070571899414 52 1.896570e+03 1.293490e+00 * time: 1.5633721351623535 53 1.896549e+03 1.877705e+00 * time: 1.579888105392456 54 1.896513e+03 2.217392e+00 * time: 1.596285104751587 55 1.896482e+03 1.658148e+00 * time: 1.6124210357666016 56 1.896466e+03 5.207217e-01 * time: 1.6380131244659424 57 1.896463e+03 1.177468e-01 * time: 1.6547160148620605 58 1.896463e+03 1.678937e-02 * time: 1.6702990531921387 59 1.896463e+03 2.666635e-03 * time: 1.6850600242614746 60 1.896463e+03 2.666635e-03 * time: 1.7247881889343262 61 1.896463e+03 2.666635e-03 * time: 1.748979091644287
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 30
Observation records: Active Missing
DV: 540 0
Total: 540 0
Number of parameters: Constant Optimized
0 8
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: NoXChange
Log-likelihood value: -1896.4632
------------------
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: 1.1920928955078125e-5 1 3.572050e+03 1.302046e+03 * time: 0.5227758884429932 2 3.266947e+03 5.384244e+02 * time: 0.6842260360717773 3 3.150635e+03 1.918079e+02 * time: 0.7098879814147949 4 3.108122e+03 1.277799e+02 * time: 0.7342989444732666 5 3.091358e+03 8.838080e+01 * time: 0.7580199241638184 6 3.082997e+03 8.321689e+01 * time: 0.78171706199646 7 3.076379e+03 8.167702e+01 * time: 0.8054249286651611 8 3.067422e+03 1.177822e+02 * time: 0.8293530941009521 9 3.048580e+03 2.526969e+02 * time: 0.9102180004119873 10 2.993096e+03 6.325396e+02 * time: 0.9483928680419922 11 2.965744e+03 7.634718e+02 * time: 1.018122911453247 12 2.921212e+03 9.704020e+02 * time: 1.060434103012085 13 2.553649e+03 6.642510e+02 * time: 1.112104892730713 14 2.319495e+03 3.204711e+02 * time: 1.1616430282592773 15 2.183040e+03 2.174905e+02 * time: 1.2002639770507812 16 2.157621e+03 3.150983e+02 * time: 1.2239129543304443 17 2.132395e+03 2.847614e+02 * time: 1.2604029178619385 18 2.084799e+03 1.563370e+02 * time: 1.2838809490203857 19 2.071497e+03 1.006429e+02 * time: 1.3066260814666748 20 2.064983e+03 9.753313e+01 * time: 1.3294000625610352 21 2.048289e+03 9.230405e+01 * time: 1.3661549091339111 22 2.020938e+03 1.292359e+02 * time: 1.3896009922027588 23 1.983888e+03 2.284990e+02 * time: 1.413412094116211 24 1.962132e+03 1.220188e+02 * time: 1.4365239143371582 25 1.945947e+03 1.035894e+02 * time: 1.4712440967559814 26 1.917782e+03 8.246698e+01 * time: 1.4968490600585938 27 1.905967e+03 5.408054e+01 * time: 1.5209400653839111 28 1.898569e+03 2.172222e+01 * time: 1.555495023727417 29 1.897473e+03 1.689350e+01 * time: 1.5768640041351318 30 1.897019e+03 1.666689e+01 * time: 1.5980639457702637 31 1.896796e+03 1.699751e+01 * time: 1.6286659240722656 32 1.896559e+03 1.645900e+01 * time: 1.6498360633850098 33 1.896223e+03 1.415504e+01 * time: 1.6710150241851807 34 1.895773e+03 1.630081e+01 * time: 1.702873945236206 35 1.895309e+03 1.723930e+01 * time: 1.7250139713287354 36 1.895004e+03 1.229983e+01 * time: 1.7471859455108643 37 1.894871e+03 5.385102e+00 * time: 1.7692198753356934 38 1.894827e+03 3.465463e+00 * time: 1.8017210960388184 39 1.894816e+03 3.387474e+00 * time: 1.8227190971374512 40 1.894807e+03 3.295388e+00 * time: 1.8429710865020752 41 1.894786e+03 3.089194e+00 * time: 1.874772071838379 42 1.894737e+03 2.928080e+00 * time: 1.8959078788757324 43 1.894624e+03 3.088723e+00 * time: 1.9167308807373047 44 1.894413e+03 3.493791e+00 * time: 1.9480400085449219 45 1.894127e+03 3.142865e+00 * time: 1.9696869850158691 46 1.893933e+03 2.145253e+00 * time: 1.9927690029144287 47 1.893888e+03 2.172800e+00 * time: 2.0143768787384033 48 1.893880e+03 2.180509e+00 * time: 2.0466489791870117 49 1.893876e+03 2.134257e+00 * time: 2.0672359466552734 50 1.893868e+03 2.032354e+00 * time: 2.088555097579956 51 1.893846e+03 1.760874e+00 * time: 2.1195690631866455 52 1.893796e+03 1.779016e+00 * time: 2.1403579711914062 53 1.893694e+03 2.018956e+00 * time: 2.161504030227661 54 1.893559e+03 2.366854e+00 * time: 2.1926639080047607 55 1.893474e+03 3.690142e+00 * time: 2.2139930725097656 56 1.893446e+03 3.675109e+00 * time: 2.234433889389038 57 1.893439e+03 3.426419e+00 * time: 2.2547478675842285 58 1.893429e+03 3.183164e+00 * time: 2.2856509685516357 59 1.893398e+03 2.695171e+00 * time: 2.3070390224456787 60 1.893328e+03 2.753548e+00 * time: 2.3288168907165527 61 1.893169e+03 3.589748e+00 * time: 2.3617630004882812 62 1.892920e+03 3.680718e+00 * time: 2.383953094482422 63 1.892667e+03 2.568107e+00 * time: 2.4067020416259766 64 1.892514e+03 1.087910e+00 * time: 2.4404640197753906 65 1.892493e+03 3.287296e-01 * time: 2.4631099700927734 66 1.892492e+03 2.967465e-01 * time: 2.486008882522583 67 1.892492e+03 3.020682e-01 * time: 2.5074350833892822 68 1.892491e+03 3.034704e-01 * time: 2.539201021194458 69 1.892491e+03 3.091846e-01 * time: 2.5601160526275635 70 1.892491e+03 3.224170e-01 * time: 2.580622911453247 71 1.892490e+03 6.494197e-01 * time: 2.6130878925323486 72 1.892488e+03 1.115188e+00 * time: 2.6345410346984863 73 1.892483e+03 1.838833e+00 * time: 2.6560049057006836 74 1.892472e+03 2.765371e+00 * time: 2.6888110637664795 75 1.892452e+03 3.463807e+00 * time: 2.7110350131988525 76 1.892431e+03 2.805270e+00 * time: 2.733109951019287 77 1.892411e+03 5.758916e-01 * time: 2.7551889419555664 78 1.892410e+03 1.434041e-01 * time: 2.7888309955596924 79 1.892409e+03 1.639246e-01 * time: 2.810991048812866 80 1.892409e+03 1.145856e-01 * time: 2.8322560787200928 81 1.892409e+03 3.966861e-02 * time: 2.864665985107422 82 1.892409e+03 3.550808e-02 * time: 2.885540008544922 83 1.892409e+03 3.456241e-02 * time: 2.905478000640869 84 1.892409e+03 3.114018e-02 * time: 2.9252779483795166 85 1.892409e+03 4.080806e-02 * time: 2.956547975540161 86 1.892409e+03 6.722726e-02 * time: 2.9769420623779297 87 1.892409e+03 1.006791e-01 * time: 2.9985098838806152 88 1.892409e+03 1.303988e-01 * time: 3.029465913772583 89 1.892409e+03 1.228919e-01 * time: 3.0491480827331543 90 1.892409e+03 6.433813e-02 * time: 3.0688600540161133 91 1.892409e+03 1.314164e-02 * time: 3.0889980792999268 92 1.892409e+03 4.929931e-04 * time: 3.1202969551086426
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 30
Observation records: Active Missing
DV: 540 0
Total: 540 0
Number of parameters: Constant Optimized
0 10
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -1892.409
--------------------
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: 1.5974044799804688e-5 1 3.641387e+03 1.432080e+03 * time: 0.5076448917388916 2 3.290450e+03 5.274782e+02 * time: 0.5392189025878906 3 3.185512e+03 2.173676e+02 * time: 0.5688400268554688 4 3.143009e+03 1.479653e+02 * time: 0.5977778434753418 5 3.128511e+03 8.980031e+01 * time: 0.749561071395874 6 3.123188e+03 5.033125e+01 * time: 0.7765090465545654 7 3.120794e+03 4.279722e+01 * time: 0.8030259609222412 8 3.118627e+03 3.971051e+01 * time: 0.8293209075927734 9 3.115300e+03 8.456587e+01 * time: 0.8555068969726562 10 3.109353e+03 1.350354e+02 * time: 0.8827428817749023 11 3.095894e+03 1.998258e+02 * time: 1.2391459941864014 12 2.988214e+03 4.366433e+02 * time: 1.2809689044952393 13 2.896081e+03 5.505943e+02 * time: 1.3751299381256104 14 2.652467e+03 7.300323e+02 * time: 3.2114789485931396 15 2.560937e+03 6.973661e+02 * time: 3.3451929092407227 16 2.254941e+03 2.740033e+02 * time: 3.3839869499206543 17 2.222509e+03 2.034303e+02 * time: 3.419456958770752 18 2.171255e+03 2.449580e+02 * time: 3.496535062789917 19 2.024532e+03 1.121511e+02 * time: 3.5247888565063477 20 1.993723e+03 1.042814e+02 * time: 3.550813913345337 21 1.985113e+03 8.079014e+01 * time: 3.5763280391693115 22 1.976757e+03 7.054196e+01 * time: 3.6015419960021973 23 1.969970e+03 6.070322e+01 * time: 3.628135919570923 24 1.961095e+03 6.810782e+01 * time: 3.6541318893432617 25 1.947983e+03 8.116920e+01 * time: 3.680279016494751 26 1.930371e+03 8.530051e+01 * time: 3.705137014389038 27 1.910209e+03 6.993170e+01 * time: 3.7351338863372803 28 1.899107e+03 3.362640e+01 * time: 3.768244981765747 29 1.898022e+03 2.642220e+01 * time: 3.8156919479370117 30 1.897055e+03 1.213144e+01 * time: 3.841033935546875 31 1.896596e+03 7.773239e+00 * time: 3.8664188385009766 32 1.896538e+03 7.997039e+00 * time: 3.890878915786743 33 1.896451e+03 8.160909e+00 * time: 3.914691925048828 34 1.896283e+03 8.237721e+00 * time: 3.938375949859619 35 1.895903e+03 1.520219e+01 * time: 3.962601900100708 36 1.895272e+03 2.358916e+01 * time: 3.987252950668335 37 1.894536e+03 2.461296e+01 * time: 4.012096881866455 38 1.893995e+03 1.546128e+01 * time: 4.054764032363892 39 1.893858e+03 6.976137e+00 * time: 4.079833030700684 40 1.893833e+03 6.019466e+00 * time: 4.1041810512542725 41 1.893786e+03 3.827201e+00 * time: 4.129745960235596 42 1.893714e+03 3.323412e+00 * time: 4.15383505821228 43 1.893592e+03 3.215150e+00 * time: 4.17781400680542 44 1.893435e+03 6.534965e+00 * time: 4.2021989822387695 45 1.893286e+03 7.424154e+00 * time: 4.242321968078613 46 1.893190e+03 5.552627e+00 * time: 4.266967058181763 47 1.893139e+03 3.222316e+00 * time: 4.291357040405273 48 1.893120e+03 3.015339e+00 * time: 4.315032958984375 49 1.893107e+03 3.244809e+00 * time: 4.338006019592285 50 1.893080e+03 6.163100e+00 * time: 4.360895872116089 51 1.893027e+03 9.824713e+00 * time: 4.398411989212036 52 1.892912e+03 1.390100e+01 * time: 4.422744035720825 53 1.892734e+03 1.510937e+01 * time: 4.4470930099487305 54 1.892561e+03 1.008563e+01 * time: 4.471181869506836 55 1.892485e+03 3.730668e+00 * time: 4.494951963424683 56 1.892471e+03 3.380261e+00 * time: 4.53122091293335 57 1.892463e+03 3.167904e+00 * time: 4.555018901824951 58 1.892441e+03 4.152065e+00 * time: 4.578194856643677 59 1.892391e+03 7.355996e+00 * time: 4.6013689041137695 60 1.892268e+03 1.195397e+01 * time: 4.637695074081421 61 1.892026e+03 1.640783e+01 * time: 4.662276983261108 62 1.891735e+03 1.593576e+01 * time: 4.686686992645264 63 1.891569e+03 8.316423e+00 * time: 4.721685886383057 64 1.891494e+03 3.948212e+00 * time: 4.745671987533569 65 1.891481e+03 3.911593e+00 * time: 4.7691810131073 66 1.891457e+03 3.875559e+00 * time: 4.802548885345459 67 1.891405e+03 3.811247e+00 * time: 4.826221942901611 68 1.891262e+03 3.657045e+00 * time: 4.8593809604644775 69 1.890930e+03 4.957405e+00 * time: 4.883656024932861 70 1.890317e+03 6.657726e+00 * time: 4.907949924468994 71 1.889660e+03 6.086302e+00 * time: 4.942239999771118 72 1.889303e+03 2.270929e+00 * time: 4.96658992767334 73 1.889253e+03 7.695301e-01 * time: 5.000185012817383 74 1.889252e+03 7.382144e-01 * time: 5.023742914199829 75 1.889251e+03 7.187898e-01 * time: 5.047044992446899 76 1.889251e+03 7.215047e-01 * time: 5.08021092414856 77 1.889250e+03 7.235155e-01 * time: 5.103350877761841 78 1.889249e+03 7.246818e-01 * time: 5.127575874328613 79 1.889244e+03 7.257796e-01 * time: 5.161609888076782 80 1.889233e+03 7.198190e-01 * time: 5.18543004989624 81 1.889204e+03 1.089029e+00 * time: 5.219203948974609 82 1.889142e+03 1.801601e+00 * time: 5.24354887008667 83 1.889043e+03 2.967917e+00 * time: 5.267734050750732 84 1.888889e+03 2.965856e+00 * time: 5.301977872848511 85 1.888705e+03 5.933555e-01 * time: 5.3263280391693115 86 1.888655e+03 9.577698e-01 * time: 5.350190877914429 87 1.888582e+03 1.498494e+00 * time: 5.3846330642700195 88 1.888533e+03 1.502751e+00 * time: 5.409010887145996 89 1.888490e+03 1.184664e+00 * time: 5.4431469440460205 90 1.888480e+03 6.684515e-01 * time: 5.467205047607422 91 1.888476e+03 3.680032e-01 * time: 5.491019010543823 92 1.888476e+03 4.720040e-01 * time: 5.52454686164856 93 1.888476e+03 4.768646e-01 * time: 5.548112869262695 94 1.888475e+03 4.736674e-01 * time: 5.581001043319702 95 1.888475e+03 4.552765e-01 * time: 5.604321002960205 96 1.888474e+03 5.193725e-01 * time: 5.628641843795776 97 1.888473e+03 8.850098e-01 * time: 5.662456035614014 98 1.888468e+03 1.461598e+00 * time: 5.686336040496826 99 1.888458e+03 2.209124e+00 * time: 5.709831953048706 100 1.888437e+03 2.961236e+00 * time: 5.743630886077881 101 1.888407e+03 2.978462e+00 * time: 5.767673969268799 102 1.888384e+03 1.707199e+00 * time: 5.801056861877441 103 1.888381e+03 6.198984e-01 * time: 5.824913024902344 104 1.888380e+03 5.171082e-01 * time: 5.84883189201355 105 1.888378e+03 1.037321e-01 * time: 5.882760047912598 106 1.888378e+03 8.473253e-02 * time: 5.905645847320557 107 1.888378e+03 8.364965e-02 * time: 5.928377866744995 108 1.888378e+03 8.080442e-02 * time: 5.961735963821411 109 1.888378e+03 7.873900e-02 * time: 5.984438896179199 110 1.888378e+03 7.798398e-02 * time: 6.018709897994995 111 1.888378e+03 7.788170e-02 * time: 6.040844917297363 112 1.888378e+03 7.776461e-02 * time: 6.063107013702393 113 1.888378e+03 9.023586e-02 * time: 6.095097064971924 114 1.888378e+03 1.631370e-01 * time: 6.1188600063323975 115 1.888378e+03 2.768691e-01 * time: 6.141672849655151 116 1.888377e+03 4.462313e-01 * time: 6.1743528842926025 117 1.888377e+03 6.643167e-01 * time: 6.197537899017334 118 1.888375e+03 8.433175e-01 * time: 6.230587005615234 119 1.888374e+03 7.596473e-01 * time: 6.253756999969482 120 1.888373e+03 3.637851e-01 * time: 6.277091026306152 121 1.888372e+03 8.305224e-02 * time: 6.309720039367676 122 1.888372e+03 2.084516e-02 * time: 6.332524061203003 123 1.888372e+03 2.056416e-02 * time: 6.355186939239502 124 1.888372e+03 2.044079e-02 * time: 6.386999845504761 125 1.888372e+03 2.035197e-02 * time: 6.408705949783325 126 1.888372e+03 2.021266e-02 * time: 6.44010591506958 127 1.888372e+03 1.998168e-02 * time: 6.46229100227356 128 1.888372e+03 3.162094e-02 * time: 6.484605073928833 129 1.888372e+03 5.510007e-02 * time: 6.516520023345947 130 1.888372e+03 9.277177e-02 * time: 6.538651943206787 131 1.888372e+03 1.528967e-01 * time: 6.560559034347534 132 1.888372e+03 2.462112e-01 * time: 6.5921900272369385 133 1.888372e+03 3.799880e-01 * time: 6.615468978881836 134 1.888371e+03 5.312357e-01 * time: 6.637573957443237 135 1.888369e+03 6.019766e-01 * time: 6.669717073440552 136 1.888367e+03 4.665348e-01 * time: 6.692209959030151 137 1.888366e+03 1.404917e-01 * time: 6.724050998687744 138 1.888365e+03 8.513280e-02 * time: 6.747042894363403 139 1.888364e+03 1.244285e-01 * time: 6.769166946411133 140 1.888364e+03 1.028231e-01 * time: 6.801823854446411 141 1.888364e+03 5.163326e-02 * time: 6.824615955352783 142 1.888364e+03 5.148616e-02 * time: 6.846876859664917 143 1.888364e+03 3.147247e-02 * time: 6.879116058349609 144 1.888364e+03 2.104597e-02 * time: 6.902534008026123 145 1.888364e+03 6.541596e-03 * time: 6.934670925140381 146 1.888364e+03 2.538502e-03 * time: 6.95622992515564 147 1.888364e+03 4.361857e-03 * time: 6.977724075317383 148 1.888364e+03 3.034845e-03 * time: 7.008401870727539 149 1.888364e+03 5.961460e-04 * time: 7.029558897018433
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 30
Observation records: Active Missing
DV: 540 0
Total: 540 0
Number of parameters: Constant Optimized
0 13
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -1888.3638
-----------------------
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.
5 Step 4 - Model Comparison
Now that we’ve fitted all of our models we need to compare them and choose one for our final model.
We begin by analyzing the model metrics. This can be done with the metrics_table function:
metrics_table(fit_base_model)WARNING: using deprecated binding Distributions.MatrixReshaped in Pumas.
, use Distributions.ReshapedDistribution{2, S, D} where D<:Distributions.Distribution{Distributions.ArrayLikeVariate{1}, S} where S<:Distributions.ValueSupport instead.17×2 DataFrame
| Row | Metric | Value |
|---|---|---|
| Any | Any | |
| 1 | Successful | true |
| 2 | Estimation Time | 3.112 |
| 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 | RawTypst("LogLikelihood (LL)") | -1901.82 |
| 12 | RawTypst("-2LL") | 3803.63 |
| 13 | AIC | 3819.63 |
| 14 | BIC | 3853.96 |
| 15 | RawTypst("(η-shrinkage) η₁") | -0.015 |
| 16 | RawTypst("(η-shrinkage) η₂") | -0.013 |
| 17 | RawTypst("(ϵ-shrinkage) DV") | 0.056 |
metrics_table(fit_covariate_model_wt)17×2 DataFrame
| Row | Metric | Value |
|---|---|---|
| Any | Any | |
| 1 | Successful | true |
| 2 | Estimation Time | 1.749 |
| 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 | RawTypst("LogLikelihood (LL)") | -1896.46 |
| 12 | RawTypst("-2LL") | 3792.93 |
| 13 | AIC | 3808.93 |
| 14 | BIC | 3843.26 |
| 15 | RawTypst("(η-shrinkage) η₁") | -0.014 |
| 16 | RawTypst("(η-shrinkage) η₂") | -0.012 |
| 17 | RawTypst("(ϵ-shrinkage) DV") | 0.056 |
metrics_table(fit_covariate_model_wt_crcl)17×2 DataFrame
| Row | Metric | Value |
|---|---|---|
| Any | Any | |
| 1 | Successful | true |
| 2 | Estimation Time | 3.12 |
| 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 | RawTypst("LogLikelihood (LL)") | -1892.41 |
| 12 | RawTypst("-2LL") | 3784.82 |
| 13 | AIC | 3804.82 |
| 14 | BIC | 3847.73 |
| 15 | RawTypst("(η-shrinkage) η₁") | -0.014 |
| 16 | RawTypst("(η-shrinkage) η₂") | -0.012 |
| 17 | RawTypst("(ϵ-shrinkage) DV") | 0.056 |
metrics_table(fit_covariate_model_wt_crcl_sex)17×2 DataFrame
| Row | Metric | Value |
|---|---|---|
| Any | Any | |
| 1 | Successful | true |
| 2 | Estimation Time | 7.03 |
| 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 | RawTypst("LogLikelihood (LL)") | -1888.36 |
| 12 | RawTypst("-2LL") | 3776.73 |
| 13 | AIC | 3802.73 |
| 14 | BIC | 3858.52 |
| 15 | RawTypst("(η-shrinkage) η₁") | -0.013 |
| 16 | RawTypst("(η-shrinkage) η₂") | -0.012 |
| 17 | RawTypst("(ϵ-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 |
|---|---|---|---|---|---|
| Any | Any | Any | Any | Any | |
| 1 | Successful | true | true | true | true |
| 2 | Estimation Time | 3.112 | 1.749 | 3.12 | 7.03 |
| 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 | RawTypst("LogLikelihood (LL)") | -1901.82 | -1896.46 | -1892.41 | -1888.36 |
| 12 | RawTypst("-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 | RawTypst("(η-shrinkage) η₁") | -0.015 | -0.014 | -0.014 | -0.013 |
| 16 | RawTypst("(η-shrinkage) η₂") | -0.013 | -0.012 | -0.012 | -0.012 |
| 17 | RawTypst("(ϵ-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.6299849528186aic(fit_covariate_model_wt)3808.9264607805967aic(fit_covariate_model_wt_crcl)3804.8179473717055aic(fit_covariate_model_wt_crcl_sex)3802.7275243740673We 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 dose control parameters. [ Info: Evaluating individual parameters. [ Info: Done.
FittedPumasModelInspection
Likelihood approximation used for weighted residuals: FOCEAnd now we pass inspect_covariate_model_wt_crcl_sex to the goodness_of_fit plotting function:
goodness_of_fit(inspect_covariate_model_wt_crcl_sex)
The idea is that the population predictions (preds) capture the general tendency of the observations while the individual predictions (ipreds) should coincide much more closely with the observations. That is exactly what we are observing in the top row subplots in our goodness of fit plot.
Regarding the bottom row, on the left we have the weighted population residuals (wres) against time, and on the right we have the weighted individual residuals (iwres) against ipreds. Here we should not see any perceived pattern, which indicates that the error model in the model has a mean 0 and constant variance. Like before, this seems to be what we are observing in our goodness of fit plot.
Hence, our covariate model with allometric scaling and different sex creatinine clearance effectw on clearance is a good candidate for our final model.
6 Time-Varying Covariates
Pumas can handle time-varying covariates. This happens automatically if, when parsing a dataset, read_pumas detects that covariate values change over time.
6.1 painord Dataset
Here’s an example with an ordinal regression using the painord dataset from PharmaDatasets.jl. :painord is our observations measuring the perceived pain in a scale from 0 to 3, which we need to have the range shifted by 1 (1 to 4). Additionally, we’ll use the concentration in plasma, :conc as a covariate. Of course, :conc varies with time, thus, it is a time-varying covariate:
painord = dataset("pumas/pain_remed")
first(painord, 5)5×8 DataFrame
| Row | id | arm | dose | time | conc | painord | dv | remed |
|---|---|---|---|---|---|---|---|---|
| Int64 | Int64 | Int64 | Float64 | Float64 | Int64 | Int64 | Int64 | |
| 1 | 1 | 2 | 20 | 0.0 | 0.0 | 3 | 0 | 0 |
| 2 | 1 | 2 | 20 | 0.5 | 1.15578 | 1 | 1 | 0 |
| 3 | 1 | 2 | 20 | 1.0 | 1.37211 | 0 | 1 | 0 |
| 4 | 1 | 2 | 20 | 1.5 | 1.30058 | 0 | 1 | 0 |
| 5 | 1 | 2 | 20 | 2.0 | 1.19195 | 1 | 1 | 0 |
@rtransform! painord :painord = :painord + 1;
describe(painord, :mean, :std, :first, :last, :eltype)8×6 DataFrame
| Row | variable | mean | std | first | last | eltype |
|---|---|---|---|---|---|---|
| Symbol | Float64 | Float64 | Real | Real | DataType | |
| 1 | id | 80.5 | 46.1992 | 1 | 160 | Int64 |
| 2 | arm | 1.5 | 1.11833 | 2 | 0 | Int64 |
| 3 | dose | 26.25 | 31.9017 | 20 | 0 | Int64 |
| 4 | time | 3.375 | 2.5183 | 0.0 | 8.0 | Float64 |
| 5 | conc | 0.93018 | 1.49902 | 0.0 | 0.0 | Float64 |
| 6 | painord | 2.50208 | 0.863839 | 4 | 4 | Int64 |
| 7 | dv | 0.508333 | 0.500061 | 0 | 0 | Int64 |
| 8 | remed | 0.059375 | 0.236387 | 0 | 0 | Int64 |
unique(painord.dose)4-element Vector{Int64}:
20
80
0
5As we can see we have 160 subjects were given either 0, 5, 20, or 80 units of a certain painkiller drug.
:conc is the drug concentration in plasma and :painord is the perceived pain in a scale from 1 to 4.
First, we’ll parse the painord dataset into a Population. Note that we’ll be using event_data=false since we do not have any dosing rows:
pop_ord =
read_pumas(painord; observations = [:painord], covariates = [:conc], event_data = false)
Note
We won’t be going into the details of the ordinal regression model in this tutorial. We highly encourage you to take a look at the Ordinal Regression Pumas Tutorial for an in-depth explanation.
We’ll build an ordinal regression model declaring :conc as a covariate. In the @derived block we’ll state the the likelihood of :painord follows a Categorical distribution:
ordinal_model = @model begin
@param begin
b₁ ∈ RealDomain(; init = 0)
b₂ ∈ RealDomain(; init = 1)
b₃ ∈ RealDomain(; init = 1)
slope ∈ RealDomain(; init = 0)
end
@covariates conc # time-varying
@pre begin
effect = slope * conc
# Logit of cumulative probabilities
lge₁ = b₁ + effect
lge₂ = lge₁ - b₂
lge₃ = lge₂ - b₃
# Probabilities of >=1 and >=2 and >=3
pge₁ = logistic(lge₁)
pge₂ = logistic(lge₂)
pge₃ = logistic(lge₃)
# Probabilities of Y=1,2,3,4
p₁ = 1.0 - pge₁
p₂ = pge₁ - pge₂
p₃ = pge₂ - pge₃
p₄ = pge₃
end
@derived begin
painord ~ @. Categorical(p₁, p₂, p₃, p₄)
end
endPumasModel
Parameters: b₁, b₂, b₃, slope
Random effects:
Covariates: conc
Dynamical system variables:
Dynamical system type: No dynamical model
Derived: painord
Observed: painordFinally we’ll fit our model using NaivePooled estimation method since it does not have any random-effects, i.e. no @random block:
ordinal_fit = fit(ordinal_model, pop_ord, init_params(ordinal_model), NaivePooled())[ Info: Checking the initial parameter values. [ Info: The initial negative log likelihood and its gradient are finite. Check passed. Iter Function value Gradient norm 0 3.103008e+03 7.031210e+02 * time: 3.1948089599609375e-5 1 2.994747e+03 1.083462e+03 * time: 0.4444308280944824 2 2.406265e+03 1.884408e+02 * time: 0.4495198726654053 3 2.344175e+03 7.741610e+01 * time: 0.454571008682251 4 2.323153e+03 2.907642e+01 * time: 0.45897984504699707 5 2.318222e+03 2.273295e+01 * time: 0.4631359577178955 6 2.316833e+03 1.390527e+01 * time: 0.46710681915283203 7 2.316425e+03 4.490883e+00 * time: 0.47116684913635254 8 2.316362e+03 9.374519e-01 * time: 0.4751608371734619 9 2.316356e+03 1.928785e-01 * time: 0.4793968200683594 10 2.316355e+03 3.119615e-02 * time: 0.4836399555206299 11 2.316355e+03 6.215513e-03 * time: 0.4879460334777832 12 2.316355e+03 8.313206e-04 * time: 0.4928100109100342
FittedPumasModel
Dynamical system type: No dynamical model
Number of subjects: 160
Observation records: Active Missing
painord: 1920 0
Total: 1920 0
Number of parameters: Constant Optimized
0 4
Likelihood approximation: NaivePooled
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -2316.3554
-----------------
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.