using Pumas
using PumasUtilities
using NCA
using NCAUtilities
A Comprehensive Introduction to Pumas
This tutorial provides a comprehensive introduction to a modeling and simulation workflow in Pumas. The idea is not to get into the details of Pumas specifics, but instead provide a narrative on the lines of a regular workflow in our day-to-day work, with brevity where required to allow a broad overview. Wherever possible, cross-references will be provided to documentation and detailed examples that provide deeper insight into a particular topic.
As part of this workflow, you will be introduced to various aspects such as:
- Data wrangling in Julia
- Exploratory analysis in Julia
- Continuous data non-linear mixed effects modeling in Pumas
- Model comparison routines, post-processing, validation etc.
1 The Study and Design
CTMNopain is a novel anti-inflammatory agent under preliminary investigation. A dose-ranging trial was conducted comparing placebo with 3 doses of CTMNopain (5mg, 20mg and 80 mg QD). The maximum tolerated dose is 160 mg per day. Plasma concentrations (mg/L) of the drug were measured at 0, 0.5, 1, 1.5, 2, 2.5, 3-8 hours.
Pain score (0=no pain, 1=mild, 2=moderate, 3=severe) were obtained at time points when plasma concentration was collected. A pain score of 2 or more is considered as no pain relief.
The subjects can request for remedication if pain relief is not achieved after 2 hours post dose. Some subjects had remedication before 2 hours if they were not able to bear the pain. The time to remedication and the remedication status is available for subjects.
The pharmacokinetic dataset can be accessed using PharmaDatasets.jl.
2 Setup
2.1 Load libraries
These libraries provide the workhorse functionality in the Pumas ecosystem:
In addition, libraries below are good add-on’s that provide ancillary functionality:
using GLM: lm, @formula
using Random
using CSV
using DataFramesMeta
using CairoMakie
using PharmaDatasets2.2 Data Wrangling
We start by reading in the dataset and making some quick summaries.
Tip
If you want to learn more about data wrangling, don’t forget to check our Data Wrangling in Julia tutorials!
pkpain_df = dataset("pk_painrelief")
first(pkpain_df, 5)5×7 DataFrame
| Row | Subject | Time | Conc | PainRelief | PainScore | RemedStatus | Dose |
|---|---|---|---|---|---|---|---|
| Int64 | Float64 | Float64 | Int64 | Int64 | Int64 | String7 | |
| 1 | 1 | 0.0 | 0.0 | 0 | 3 | 1 | 20 mg |
| 2 | 1 | 0.5 | 1.15578 | 1 | 1 | 0 | 20 mg |
| 3 | 1 | 1.0 | 1.37211 | 1 | 0 | 0 | 20 mg |
| 4 | 1 | 1.5 | 1.30058 | 1 | 0 | 0 | 20 mg |
| 5 | 1 | 2.0 | 1.19195 | 1 | 1 | 0 | 20 mg |
Let’s filter out the placebo data as we don’t need that for the PK analysis.
pkpain_noplb_df = @rsubset pkpain_df :Dose != "Placebo";
first(pkpain_noplb_df, 5)5×7 DataFrame
| Row | Subject | Time | Conc | PainRelief | PainScore | RemedStatus | Dose |
|---|---|---|---|---|---|---|---|
| Int64 | Float64 | Float64 | Int64 | Int64 | Int64 | String7 | |
| 1 | 1 | 0.0 | 0.0 | 0 | 3 | 1 | 20 mg |
| 2 | 1 | 0.5 | 1.15578 | 1 | 1 | 0 | 20 mg |
| 3 | 1 | 1.0 | 1.37211 | 1 | 0 | 0 | 20 mg |
| 4 | 1 | 1.5 | 1.30058 | 1 | 0 | 0 | 20 mg |
| 5 | 1 | 2.0 | 1.19195 | 1 | 1 | 0 | 20 mg |
3 Analysis
3.1 Non-compartmental analysis
Let’s begin by performing a quick NCA of the concentration time profiles and view the exposure changes across doses. The input data specification for NCA analysis requires the presence of a :route column and an :amt column that specifies the dose. So, let’s add that in:
@rtransform! pkpain_noplb_df begin
:route = "ev"
:Dose = parse(Int, chop(:Dose; tail = 3))
endWe also need to create an :amt column:
@rtransform! pkpain_noplb_df :amt = :Time == 0 ? :Dose : missingNow, we map the data variables to the read_nca function that prepares the data for NCA analysis.
pkpain_nca = read_nca(
pkpain_noplb_df;
id = :Subject,
time = :Time,
amt = :amt,
observations = :Conc,
group = [:Dose],
route = :route,
)NCAPopulation (120 subjects):
Group: [["Dose" => 5], ["Dose" => 20], ["Dose" => 80]]
Number of missing observations: 0
Number of blq observations: 0Now that we mapped the data in, let’s visualize the concentration vs time plots for a few individuals. When paginate is set to true, a vector of plots are returned and below we display the first element with 9 individuals.
f = observations_vs_time(
pkpain_nca;
paginate = true,
axis = (; xlabel = "Time (hr)", ylabel = "CTMNoPain Concentration (ng/mL)"),
facet = (; combinelabels = true),
)
f[1]
or you can view the summary curves by dose group as passed in to the group argument in read_nca
summary_observations_vs_time(
pkpain_nca,
figure = (; fontsize = 22, resolution = (800, 1000)),
color = "black",
linewidth = 3,
axis = (; xlabel = "Time (hr)", ylabel = "CTMX Concentration (μg/mL)"),
facet = (; combinelabels = true, linkaxes = true),
)
A full NCA Report is now obtained for completeness purposes using the run_nca function, but later we will only extract a couple of key metrics of interest.
pk_nca = run_nca(pkpain_nca; sigdigits = 3)We can look at the NCA fits for some subjects. Here f is a vector or figures. We’ll showcase the first image by indexing f:
f = subject_fits(
pk_nca,
paginate = true,
axis = (; xlabel = "Time (hr)", ylabel = "CTMX Concentration (μg/mL)"),
facet = (; combinelabels = true, linkaxes = true),
)
f[1]
As CTMNopain’s effect maybe mainly related to maximum concentration (cmax) or area under the curve (auc), we present some summary statistics using the summarize function from NCA.
strata = [:Dose]1-element Vector{Symbol}:
:Doseparams = [:cmax, :aucinf_obs]2-element Vector{Symbol}:
:cmax
:aucinf_obsoutput = summarize(pk_nca; stratify_by = strata, parameters = params)6×10 DataFrame
| Row | Dose | parameters | numsamples | minimum | maximum | mean | std | geomean | geostd | geomeanCV |
|---|---|---|---|---|---|---|---|---|---|---|
| Int64 | String | Int64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | 5 | cmax | 40 | 0.19 | 0.539 | 0.356075 | 0.0884129 | 0.345104 | 1.2932 | 26.1425 |
| 2 | 5 | aucinf_obs | 40 | 0.914 | 3.4 | 1.5979 | 0.490197 | 1.53373 | 1.32974 | 29.0868 |
| 3 | 20 | cmax | 40 | 0.933 | 2.7 | 1.4737 | 0.361871 | 1.43408 | 1.2633 | 23.6954 |
| 4 | 20 | aucinf_obs | 40 | 2.77 | 14.1 | 6.377 | 2.22239 | 6.02031 | 1.41363 | 35.6797 |
| 5 | 80 | cmax | 40 | 3.3 | 8.47 | 5.787 | 1.31957 | 5.64164 | 1.25757 | 23.2228 |
| 6 | 80 | aucinf_obs | 40 | 13.7 | 49.1 | 29.5 | 8.68984 | 28.2954 | 1.34152 | 30.0258 |
The statistics printed above are the default, but you can pass in your own statistics using the stats = [] argument to the summarize function.
We can look at a few parameter distribution plots.
parameters_vs_group(
pk_nca,
parameter = :cmax,
axis = (; xlabel = "Dose (mg)", ylabel = "Cₘₐₓ (ng/mL)"),
figure = (; fontsize = 18),
)
Dose normalized PK parameters, cmax and aucinf were essentially dose proportional between for 5 mg, 20 mg and 80 mg doses. You can perform a simple regression to check the impact of dose on cmax:
dp = NCA.DoseLinearityPowerModel(pk_nca, :cmax; level = 0.9)Dose Linearity Power Model
Variable: cmax
Model: log(cmax) ~ log(α) + β × log(dose)
────────────────────────────────────
Estimate low CI 90% high CI 90%
────────────────────────────────────
β 1.00775 0.97571 1.0398
────────────────────────────────────Here’s a visualization for the dose linearity using a power model for cmax:
power_model(dp)
We can also visualize a dose proportionality results with respect to a specific endpoint in a NCA Report; for example cmax and aucinf_obs:
dose_vs_dose_normalized(pk_nca, :cmax)
dose_vs_dose_normalized(pk_nca, :aucinf_obs)
Based on visual inspection of the concentration time profiles as seen earlier, CTMNopain exhibited monophasic decline, and perhaps a one compartment model best fits the PK data.
3.2 Pharmacokinetic modeling
As seen from the plots above, the concentrations decline monoexponentially. We will evaluate both one and two compartment structural models to assess best fit. Further, different residual error models will also be tested.
We will use the results from NCA to provide us good initial estimates.
3.2.1 Data preparation for modeling
PumasNDF requires the presence of :evid and :cmt columns in the dataset.
@rtransform! pkpain_noplb_df begin
:evid = :Time == 0 ? 1 : 0
:cmt = :Time == 0 ? 1 : 2
:cmt2 = 1 # for zero order absorption
endFurther, observations at time of dosing, i.e., when evid = 1 have to be missing
@rtransform! pkpain_noplb_df :Conc = :evid == 1 ? missing : :ConcThe dataframe will now be converted to a Population using read_pumas. Note that both observations and covariates are required to be an array even if it is one element.
pkpain_noplb = read_pumas(
pkpain_noplb_df;
id = :Subject,
time = :Time,
amt = :amt,
observations = [:Conc],
covariates = [:Dose],
evid = :evid,
cmt = :cmt,
)Population
Subjects: 120
Covariates: Dose
Observations: ConcNow that the data is transformed to a Population of subjects, we can explore different models.
3.2.2 One-compartment model
Note
If you are not familiar yet with the @model blocks and syntax, please check our documentation.
pk_1cmp = @model begin
@metadata begin
desc = "One Compartment Model"
timeu = u"hr"
end
@param begin
"""
Clearance (L/hr)
"""
tvcl ∈ RealDomain(; lower = 0, init = 3.2)
"""
Volume (L)
"""
tvv ∈ RealDomain(; lower = 0, init = 16.4)
"""
Absorption rate constant (h-1)
"""
tvka ∈ RealDomain(; lower = 0, init = 3.8)
"""
- ΩCL
- ΩVc
- ΩKa
"""
Ω ∈ PDiagDomain(init = [0.04, 0.04, 0.04])
"""
Proportional RUV
"""
σ_p ∈ RealDomain(; lower = 0.0001, init = 0.2)
end
@random begin
η ~ MvNormal(Ω)
end
@covariates begin
"""
Dose (mg)
"""
Dose
end
@pre begin
CL = tvcl * exp(η[1])
Vc = tvv * exp(η[2])
Ka = tvka * exp(η[3])
end
@dynamics Depots1Central1
@derived begin
cp := @. Central / Vc
"""
CTMx Concentration (ng/mL)
"""
Conc ~ @. Normal(cp, abs(cp) * σ_p)
end
endPumasModel
Parameters: tvcl, tvv, tvka, Ω, σ_p
Random effects: η
Covariates: Dose
Dynamical variables: Depot, Central
Derived: Conc
Observed: Conc
Tip
Note that the local assignment := can be used to define intermediate statements that will not be carried outside of the block. This means that all the resulting data workflows from this model will not contain the intermediate variables defined with :=. We use this when we want to suppress the variable from any further output.
The idea behind := is for performance reasons. If you are not carrying the variable defined with := outside of the block, then it is not necessary to store it in the resulting data structures. Not only will your model run faster, but the resulting data structures will also be smaller.
Before going to fit the model, let’s evaluate some helpful steps via simulation to check appropriateness of data and model
# zero out the random effects
etas = zero_randeffs(pk_1cmp, init_params(pk_1cmp), pkpain_noplb)Above, we are generating a vector of η’s of the same length as the number of subjects to zero out the random effects. We do this as we are evaluating the trajectories of the concentrations at the initial set of parameters at a population level. Other helper functions here are sample_randeffs and init_randeffs. Please refer to the documentation.
simpk_iparams = simobs(pk_1cmp, pkpain_noplb, init_params(pk_1cmp), etas)Simulated population (Vector{<:Subject})
Simulated subjects: 120
Simulated variables: Concsim_plot(
pk_1cmp,
simpk_iparams;
observations = [:Conc],
figure = (; fontsize = 18),
axis = (;
xlabel = "Time (hr)",
ylabel = "Observed/Predicted \n CTMx Concentration (ng/mL)",
),
)
Our NCA based initial guess on the parameters seem to work well.
Lets change the initial estimate of a couple of the parameters to evaluate the sensitivity.
pkparam = (; init_params(pk_1cmp)..., tvka = 2, tvv = 10)(tvcl = 3.2,
tvv = 10,
tvka = 2,
Ω = [0.04 0.0 0.0; 0.0 0.04 0.0; 0.0 0.0 0.04],
σ_p = 0.2,)simpk_changedpars = simobs(pk_1cmp, pkpain_noplb, pkparam, etas)Simulated population (Vector{<:Subject})
Simulated subjects: 120
Simulated variables: Concsim_plot(
pk_1cmp,
simpk_changedpars;
observations = [:Conc],
figure = (; fontsize = 18),
axis = (
xlabel = "Time (hr)",
ylabel = "Observed/Predicted \n CTMx Concentration (ng/mL)",
),
)
Changing the tvka and decreasing the tvv seemed to make an impact and observations go through the simulated lines.
To get a quick ballpark estimate of your PK parameters, we can do a NaivePooled analysis.
3.2.2.1 NaivePooled
pkfit_np = fit(pk_1cmp, pkpain_noplb, init_params(pk_1cmp), NaivePooled(); omegas = (:Ω,))[ 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.744356e+02 3.715711e+03
* time: 0.017796039581298828
1 2.343899e+02 1.747348e+03
* time: 0.3347740173339844
2 9.696232e+01 1.198088e+03
* time: 0.3362419605255127
3 -7.818699e+01 5.538151e+02
* time: 0.33734583854675293
4 -1.234803e+02 2.462514e+02
* time: 0.33852386474609375
5 -1.372888e+02 2.067458e+02
* time: 0.3396739959716797
6 -1.410579e+02 1.162950e+02
* time: 0.3807668685913086
7 -1.434754e+02 5.632816e+01
* time: 0.3816859722137451
8 -1.453401e+02 7.859270e+01
* time: 0.3825080394744873
9 -1.498185e+02 1.455606e+02
* time: 0.38326287269592285
10 -1.534371e+02 1.303682e+02
* time: 0.3840329647064209
11 -1.563557e+02 5.975474e+01
* time: 0.38487982749938965
12 -1.575052e+02 9.308611e+00
* time: 0.3857858180999756
13 -1.579357e+02 1.234484e+01
* time: 0.38680601119995117
14 -1.581874e+02 7.478196e+00
* time: 0.3878359794616699
15 -1.582981e+02 2.027162e+00
* time: 0.3888590335845947
16 -1.583375e+02 5.578262e+00
* time: 0.3899378776550293
17 -1.583556e+02 4.727050e+00
* time: 0.39102888107299805
18 -1.583644e+02 2.340173e+00
* time: 0.3921089172363281
19 -1.583680e+02 7.738100e-01
* time: 0.3932008743286133
20 -1.583696e+02 3.300689e-01
* time: 0.3943350315093994
21 -1.583704e+02 3.641985e-01
* time: 0.3954808712005615
22 -1.583707e+02 4.365901e-01
* time: 0.3966329097747803
23 -1.583709e+02 3.887800e-01
* time: 0.3977530002593994
24 -1.583710e+02 2.766977e-01
* time: 0.3990499973297119
25 -1.583710e+02 1.758029e-01
* time: 0.40010905265808105
26 -1.583710e+02 1.133947e-01
* time: 0.40117788314819336
27 -1.583710e+02 7.922544e-02
* time: 0.4022488594055176
28 -1.583710e+02 5.954998e-02
* time: 0.40346193313598633
29 -1.583710e+02 4.157079e-02
* time: 0.4046509265899658
30 -1.583710e+02 4.295447e-02
* time: 0.40584397315979004
31 -1.583710e+02 5.170754e-02
* time: 0.4069998264312744
32 -1.583710e+02 2.644385e-02
* time: 0.4085419178009033
33 -1.583710e+02 4.548999e-03
* time: 0.41005396842956543
34 -1.583710e+02 2.501802e-02
* time: 0.4115719795227051
35 -1.583710e+02 3.763441e-02
* time: 0.4127190113067627
36 -1.583710e+02 3.206026e-02
* time: 0.4138638973236084
37 -1.583710e+02 1.003695e-02
* time: 0.41501903533935547
38 -1.583710e+02 2.209094e-02
* time: 0.4161999225616455
39 -1.583710e+02 4.954210e-03
* time: 0.4173769950866699
40 -1.583710e+02 1.609377e-02
* time: 0.44985198974609375
41 -1.583710e+02 1.579798e-02
* time: 0.4507119655609131
42 -1.583710e+02 1.014086e-03
* time: 0.45177793502807617
43 -1.583710e+02 6.050530e-03
* time: 0.4528930187225342
44 -1.583710e+02 1.354438e-02
* time: 0.4536750316619873
45 -1.583710e+02 4.473256e-03
* time: 0.4544379711151123
46 -1.583710e+02 4.644149e-03
* time: 0.4552340507507324
47 -1.583710e+02 9.831413e-03
* time: 0.455996036529541
48 -1.583710e+02 1.047835e-03
* time: 0.4567749500274658
49 -1.583710e+02 8.361009e-03
* time: 0.4575538635253906
50 -1.583710e+02 7.909388e-04
* time: 0.4583289623260498FittedPumasModel
Successful minimization: true
Likelihood approximation: NaivePooled
Log-likelihood value: 158.37103
Number of subjects: 120
Number of parameters: Fixed Optimized
1 6
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
------------------
Estimate
------------------
tvcl 3.0054
tvv 14.089
tvka 44.228
Ω₁,₁ 0.0
Ω₂,₂ 0.0
Ω₃,₃ 0.0
σ_p 0.32999
------------------coefficients_table(pkfit_np)7×3 DataFrame
| Row | Parameter | Description | Estimate |
|---|---|---|---|
| String | Abstract… | Float64 | |
| 1 | tvcl | Clearance (L/hr)\n | 3.005 |
| 2 | tvv | Volume (L)\n | 14.089 |
| 3 | tvka | Absorption rate constant (h-1)\n | 44.228 |
| 4 | Ω₁,₁ | ΩCL | 0.0 |
| 5 | Ω₂,₂ | ΩVc | 0.0 |
| 6 | Ω₃,₃ | ΩKa | 0.0 |
| 7 | σ_p | Proportional RUV\n | 0.33 |
The final estimates from the NaivePooled approach seem reasonably close to our initial guess from NCA, except for the tvka parameter. We will stick with our initial guess.
One way to be cautious before going into a complete fitting routine is to evaluate the likelihood of the individual subjects given the initial parameter values and see if any subject(s) pops out as unreasonable. There are a few ways of doing this:
- check the
loglikelihoodsubject wise - check if there any influential subjects
Below, we are basically checking if the initial estimates for any subject are way off that we are unable to compute the initial loglikelihood.
lls = []
for subj in pkpain_noplb
push!(lls, loglikelihood(pk_1cmp, subj, pkparam, FOCE()))
end
# the plot below is using native CairoMakie `hist`
hist(lls; bins = 10, normalization = :none, color = (:black, 0.5), x_gap = 0)
The distribution of the loglikelihood’s suggest no extreme outliers.
A more convenient way is to use the findinfluential function that provides a list of k top influential subjects by showing the normalized (minus) loglikelihood for each subject. As you can see below, the minus loglikelihood in the range of 16 agrees with the histogram plotted above.
influential_subjects = findinfluential(pk_1cmp, pkpain_noplb, pkparam, FOCE())5-element Vector{NamedTuple{(:id, :nll), Tuple{String, Float64}}}:
(id = "148", nll = 16.6596588568447)
(id = "135", nll = 16.648985190076335)
(id = "156", nll = 15.959069556607497)
(id = "159", nll = 15.441218240496482)
(id = "149", nll = 14.715134644119514)3.2.2.2 FOCE
Now that we have a good handle on our data, lets go ahead and fit a population model with FOCE:
pkfit_1cmp = fit(pk_1cmp, pkpain_noplb, pkparam, FOCE(); constantcoef = (; tvka = 2))[ 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 -5.935351e+02 5.597318e+02
* time: 6.794929504394531e-5
1 -7.022088e+02 1.707063e+02
* time: 0.16498708724975586
2 -7.314067e+02 2.903269e+02
* time: 0.22076416015625
3 -8.520591e+02 2.285888e+02
* time: 0.27546000480651855
4 -1.120191e+03 3.795410e+02
* time: 0.39693212509155273
5 -1.178784e+03 2.323978e+02
* time: 0.4663691520690918
6 -1.218320e+03 9.699907e+01
* time: 0.5187010765075684
7 -1.223641e+03 5.862105e+01
* time: 0.5688879489898682
8 -1.227620e+03 1.831402e+01
* time: 0.619326114654541
9 -1.228381e+03 2.132323e+01
* time: 0.6671350002288818
10 -1.230098e+03 2.921228e+01
* time: 0.7165780067443848
11 -1.230854e+03 2.029661e+01
* time: 0.7649531364440918
12 -1.231116e+03 5.229098e+00
* time: 0.8119821548461914
13 -1.231179e+03 1.689231e+00
* time: 0.8577511310577393
14 -1.231187e+03 1.215379e+00
* time: 0.9022819995880127
15 -1.231188e+03 2.770380e-01
* time: 0.94405198097229
16 -1.231188e+03 1.636650e-01
* time: 0.9807579517364502
17 -1.231188e+03 2.701138e-01
* time: 1.005876064300537
18 -1.231188e+03 3.163347e-01
* time: 1.0426220893859863
19 -1.231188e+03 1.505241e-01
* time: 1.0807011127471924
20 -1.231188e+03 2.484090e-02
* time: 1.115065097808838
21 -1.231188e+03 8.344982e-04
* time: 1.1475889682769775FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Log-likelihood value: 1231.188
Number of subjects: 120
Number of parameters: Fixed Optimized
1 6
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
-------------------
Estimate
-------------------
tvcl 3.1642
tvv 13.288
tvka 2.0
Ω₁,₁ 0.08494
Ω₂,₂ 0.048568
Ω₃,₃ 5.5811
σ_p 0.10093
-------------------infer(pkfit_1cmp)[ Info: Calculating: variance-covariance matrix.
[ Info: Done.Asymptotic inference results using sandwich estimator
Successful minimization: true
Likelihood approximation: FOCE
Log-likelihood value: 1231.188
Number of subjects: 120
Number of parameters: Fixed Optimized
1 6
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
-------------------------------------------------------------------
Estimate SE 95.0% C.I.
-------------------------------------------------------------------
tvcl 3.1642 0.086619 [ 2.9944 ; 3.334 ]
tvv 13.288 0.27481 [12.749 ; 13.827 ]
tvka 2.0 NaN [ NaN ; NaN ]
Ω₁,₁ 0.08494 0.011022 [ 0.063338; 0.10654 ]
Ω₂,₂ 0.048568 0.0063501 [ 0.036122; 0.061014]
Ω₃,₃ 5.5811 1.2194 [ 3.1911 ; 7.9711 ]
σ_p 0.10093 0.0057196 [ 0.089718; 0.11214 ]
-------------------------------------------------------------------Notice that tvka is fixed to 2 as we don’t have a lot of information before tmax. From the results above, we see that the parameter precision for this model is reasonable.
3.2.3 Two-compartment model
Just to be sure, let’s fit a 2-compartment model and evaluate:
pk_2cmp = @model begin
@param begin
"""
Clearance (L/hr)
"""
tvcl ∈ RealDomain(; lower = 0, init = 3.2)
"""
Central Volume (L)
"""
tvv ∈ RealDomain(; lower = 0, init = 16.4)
"""
Peripheral Volume (L)
"""
tvvp ∈ RealDomain(; lower = 0, init = 10)
"""
Distributional Clearance (L/hr)
"""
tvq ∈ RealDomain(; lower = 0, init = 2)
"""
Absorption rate constant (h-1)
"""
tvka ∈ RealDomain(; lower = 0, init = 1.3)
"""
- ΩCL
- ΩVc
- ΩKa
- ΩVp
- ΩQ
"""
Ω ∈ PDiagDomain(init = [0.04, 0.04, 0.04, 0.04, 0.04])
"""
Proportional RUV
"""
σ_p ∈ RealDomain(; lower = 0.0001, init = 0.2)
end
@random begin
η ~ MvNormal(Ω)
end
@covariates begin
"""
Dose (mg)
"""
Dose
end
@pre begin
CL = tvcl * exp(η[1])
Vc = tvv * exp(η[2])
Ka = tvka * exp(η[3])
Vp = tvvp * exp(η[4])
Q = tvq * exp(η[5])
end
@dynamics Depots1Central1Periph1
@derived begin
cp := @. Central / Vc
"""
CTMx Concentration (ng/mL)
"""
Conc ~ @. Normal(cp, cp * σ_p)
end
endPumasModel
Parameters: tvcl, tvv, tvvp, tvq, tvka, Ω, σ_p
Random effects: η
Covariates: Dose
Dynamical variables: Depot, Central, Peripheral
Derived: Conc
Observed: Conc3.2.3.1 FOCE
pkfit_2cmp =
fit(pk_2cmp, pkpain_noplb, init_params(pk_2cmp), FOCE(); constantcoef = (; tvka = 2))[ 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 -6.302369e+02 1.021050e+03
* time: 6.103515625e-5
1 -9.197817e+02 9.927951e+02
* time: 0.12894105911254883
2 -1.372640e+03 2.054986e+02
* time: 0.25276613235473633
3 -1.446326e+03 1.543987e+02
* time: 0.3791351318359375
4 -1.545570e+03 1.855028e+02
* time: 0.49924802780151367
5 -1.581449e+03 1.713157e+02
* time: 0.6875019073486328
6 -1.639433e+03 1.257382e+02
* time: 0.8002951145172119
7 -1.695964e+03 7.450539e+01
* time: 0.9222500324249268
8 -1.722243e+03 5.961044e+01
* time: 1.0352609157562256
9 -1.736883e+03 7.320921e+01
* time: 1.1567559242248535
10 -1.753547e+03 7.501938e+01
* time: 1.2730131149291992
11 -1.764053e+03 6.185661e+01
* time: 1.4009671211242676
12 -1.778991e+03 4.831033e+01
* time: 1.5344030857086182
13 -1.791492e+03 4.943278e+01
* time: 1.6653239727020264
14 -1.799847e+03 2.871410e+01
* time: 1.8146750926971436
15 -1.805374e+03 7.520789e+01
* time: 1.9650909900665283
16 -1.816260e+03 2.990621e+01
* time: 2.111330032348633
17 -1.818252e+03 2.401915e+01
* time: 2.238595962524414
18 -1.822988e+03 2.587225e+01
* time: 2.360322952270508
19 -1.824653e+03 1.550517e+01
* time: 2.4876959323883057
20 -1.826074e+03 1.788927e+01
* time: 2.6025099754333496
21 -1.826821e+03 1.888389e+01
* time: 2.7291760444641113
22 -1.827900e+03 1.432840e+01
* time: 2.8555569648742676
23 -1.828511e+03 9.422040e+00
* time: 2.976840019226074
24 -1.828754e+03 5.363445e+00
* time: 3.10815691947937
25 -1.828862e+03 4.916168e+00
* time: 3.224151134490967
26 -1.829007e+03 4.695750e+00
* time: 3.353307008743286
27 -1.829358e+03 1.090244e+01
* time: 3.4827399253845215
28 -1.829830e+03 1.451320e+01
* time: 3.6071200370788574
29 -1.830201e+03 1.108694e+01
* time: 3.7400460243225098
30 -1.830360e+03 2.892316e+00
* time: 3.8706181049346924
31 -1.830390e+03 1.699262e+00
* time: 3.989469051361084
32 -1.830404e+03 1.602221e+00
* time: 4.109914064407349
33 -1.830432e+03 2.823439e+00
* time: 4.225950002670288
34 -1.830475e+03 4.118415e+00
* time: 4.351835012435913
35 -1.830527e+03 5.082915e+00
* time: 4.473376989364624
36 -1.830591e+03 2.670079e+00
* time: 4.608508110046387
37 -1.830615e+03 3.512024e+00
* time: 4.736472129821777
38 -1.830623e+03 2.286718e+00
* time: 4.857501029968262
39 -1.830625e+03 1.670870e+00
* time: 4.981276035308838
40 -1.830627e+03 9.659338e-01
* time: 5.087990045547485
41 -1.830628e+03 9.247684e-01
* time: 5.207504034042358
42 -1.830628e+03 3.479743e-01
* time: 5.312866926193237
43 -1.830629e+03 4.506560e-01
* time: 5.414853096008301
44 -1.830630e+03 6.781412e-01
* time: 5.531583070755005
45 -1.830630e+03 4.430775e-01
* time: 5.637785911560059
46 -1.830630e+03 8.918801e-02
* time: 5.737255096435547
47 -1.830630e+03 2.405868e-03
* time: 5.838191032409668
48 -1.830630e+03 1.870239e-03
* time: 5.923939943313599
49 -1.830630e+03 1.873124e-03
* time: 6.014481067657471
50 -1.830630e+03 1.856009e-03
* time: 6.132394075393677
51 -1.830630e+03 1.856654e-03
* time: 6.251924991607666
52 -1.830630e+03 1.856385e-03
* time: 6.356801986694336
53 -1.830630e+03 1.853241e-03
* time: 6.4754719734191895
54 -1.830630e+03 1.853241e-03
* time: 6.61060094833374
55 -1.830630e+03 1.853241e-03
* time: 6.745827913284302FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Log-likelihood value: 1830.6305
Number of subjects: 120
Number of parameters: Fixed Optimized
1 10
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
-------------------
Estimate
-------------------
tvcl 2.8138
tvv 11.005
tvvp 5.54
tvq 1.5159
tvka 2.0
Ω₁,₁ 0.10267
Ω₂,₂ 0.060776
Ω₃,₃ 1.2012
Ω₄,₄ 0.42349
Ω₅,₅ 0.24473
σ_p 0.048405
-------------------3.3 Comparing One- versus Two-compartment models
The 2-compartment model has a much lower objective function compared to the 1-compartment. Let’s compare the estimates from the 2 models using the compare_estimates function.
compare_estimates(; pkfit_1cmp, pkfit_2cmp)11×3 DataFrame
| Row | parameter | pkfit_1cmp | pkfit_2cmp |
|---|---|---|---|
| String | Float64? | Float64? | |
| 1 | tvcl | 3.1642 | 2.81378 |
| 2 | tvv | 13.288 | 11.0046 |
| 3 | tvka | 2.0 | 2.0 |
| 4 | Ω₁,₁ | 0.0849405 | 0.102669 |
| 5 | Ω₂,₂ | 0.0485682 | 0.0607756 |
| 6 | Ω₃,₃ | 5.58107 | 1.20116 |
| 7 | σ_p | 0.100928 | 0.0484049 |
| 8 | tvvp | missing | 5.53997 |
| 9 | tvq | missing | 1.51591 |
| 10 | Ω₄,₄ | missing | 0.423493 |
| 11 | Ω₅,₅ | missing | 0.244732 |
We perform a likelihood ratio test to compare the two nested models. The test statistic and the \(p\)-value clearly indicate that a 2-compartment model should be preferred.
lrtest(pkfit_1cmp, pkfit_2cmp)Statistic: 1200.0
Degrees of freedom: 4
P-value: 0.0We should also compare the other metrics and statistics, such ηshrinkage, ϵshrinkage, aic, and bic using the metrics_table function.
@chain metrics_table(pkfit_2cmp) begin
leftjoin(metrics_table(pkfit_1cmp); on = :Metric, makeunique = true)
rename!(:Value => :pk2cmp, :Value_1 => :pk1cmp)
end20×3 DataFrame
| Row | Metric | pk2cmp | pk1cmp |
|---|---|---|---|
| String | Any | Any | |
| 1 | Successful | true | true |
| 2 | Estimation Time | 6.746 | 1.148 |
| 3 | Subjects | 120 | 120 |
| 4 | Fixed Parameters | 1 | 1 |
| 5 | Optimized Parameters | 10 | 6 |
| 6 | Conc Active Observations | 1320 | 1320 |
| 7 | Conc Missing Observations | 0 | 0 |
| 8 | Total Active Observations | 1320 | 1320 |
| 9 | Total Missing Observations | 0 | 0 |
| 10 | Likelihood Approximation | Pumas.FOCE | Pumas.FOCE |
| 11 | LogLikelihood (LL) | 1830.63 | 1231.19 |
| 12 | -2LL | -3661.26 | -2462.38 |
| 13 | AIC | -3641.26 | -2450.38 |
| 14 | BIC | -3589.41 | -2419.26 |
| 15 | (η-shrinkage) η₁ | 0.037 | 0.016 |
| 16 | (η-shrinkage) η₂ | 0.047 | 0.04 |
| 17 | (η-shrinkage) η₃ | 0.516 | 0.733 |
| 18 | (ϵ-shrinkage) Conc | 0.185 | 0.105 |
| 19 | (η-shrinkage) η₄ | 0.287 | missing |
| 20 | (η-shrinkage) η₅ | 0.154 | missing |
We next generate some goodness of fit plots to compare which model is performing better. To do this, we first inspect the diagnostics of our model fit.
res_inspect_1cmp = inspect(pkfit_1cmp)[ 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()res_inspect_2cmp = inspect(pkfit_2cmp)[ 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()gof_1cmp = goodness_of_fit(res_inspect_1cmp; figure = (; fontsize = 12))
gof_2cmp = goodness_of_fit(res_inspect_2cmp; figure = (; fontsize = 12))
These plots clearly indicate that the 2-compartment model is a better fit compared to the 1-compartment model.
We can look at selected sample of individual plots.
fig_subject_fits = subject_fits(
res_inspect_2cmp;
separate = true,
paginate = true,
facet = (; combinelabels = true),
figure = (; fontsize = 18),
axis = (; xlabel = "Time (hr)", ylabel = "CTMx Concentration (ng/mL)"),
)
fig_subject_fits[1]
There a lot of important plotting functions you can use for your standard model diagnostics. Please make sure to read the documentation for plotting. Below, we are checking the distribution of the empirical Bayes estimates.
empirical_bayes_dist(res_inspect_2cmp; zeroline_color = :red)
empirical_bayes_vs_covariates(
res_inspect_2cmp;
categorical = [:Dose],
figure = (; resolution = (600, 800)),
)
Clearly, our guess at tvka seems off-target. Let’s try and estimate tvka instead of fixing it to 2:
pkfit_2cmp_unfix_ka = fit(pk_2cmp, pkpain_noplb, init_params(pk_2cmp), 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 -3.200734e+02 1.272671e+03
* time: 7.891654968261719e-5
1 -8.682982e+02 1.000199e+03
* time: 0.19209694862365723
2 -1.381870e+03 5.008081e+02
* time: 0.35257697105407715
3 -1.551053e+03 6.833490e+02
* time: 0.5508878231048584
4 -1.680887e+03 1.834586e+02
* time: 0.6971509456634521
5 -1.726118e+03 8.870274e+01
* time: 0.8735129833221436
6 -1.761023e+03 1.162036e+02
* time: 1.0399467945098877
7 -1.786619e+03 1.114552e+02
* time: 1.2018718719482422
8 -1.863556e+03 9.914305e+01
* time: 1.3825809955596924
9 -1.882942e+03 5.342676e+01
* time: 1.56693696975708
10 -1.888020e+03 2.010181e+01
* time: 1.7427558898925781
11 -1.889832e+03 1.867262e+01
* time: 1.911958932876587
12 -1.891649e+03 1.668510e+01
* time: 2.078805923461914
13 -1.892615e+03 1.820707e+01
* time: 2.256208896636963
14 -1.893453e+03 1.745193e+01
* time: 2.4218368530273438
15 -1.894760e+03 1.850174e+01
* time: 2.589305877685547
16 -1.895647e+03 1.773921e+01
* time: 2.75596284866333
17 -1.896597e+03 1.143421e+01
* time: 2.9222638607025146
18 -1.897114e+03 9.720034e+00
* time: 3.0920848846435547
19 -1.897373e+03 6.054160e+00
* time: 3.257701873779297
20 -1.897498e+03 3.985923e+00
* time: 3.4223248958587646
21 -1.897571e+03 4.262502e+00
* time: 3.5868899822235107
22 -1.897633e+03 4.010316e+00
* time: 3.7477738857269287
23 -1.897714e+03 4.805389e+00
* time: 3.892836809158325
24 -1.897802e+03 3.508614e+00
* time: 4.0554609298706055
25 -1.897865e+03 3.691472e+00
* time: 4.213900804519653
26 -1.897900e+03 2.982676e+00
* time: 4.372879981994629
27 -1.897928e+03 2.563863e+00
* time: 4.534641981124878
28 -1.897968e+03 3.261530e+00
* time: 4.692546844482422
29 -1.898013e+03 3.064695e+00
* time: 4.8509509563446045
30 -1.898040e+03 1.636456e+00
* time: 5.014995813369751
31 -1.898051e+03 1.439998e+00
* time: 5.176599979400635
32 -1.898057e+03 1.436505e+00
* time: 5.320152997970581
33 -1.898069e+03 1.881592e+00
* time: 5.478768825531006
34 -1.898095e+03 3.253228e+00
* time: 5.640621900558472
35 -1.898142e+03 4.257954e+00
* time: 5.814818859100342
36 -1.898199e+03 3.685153e+00
* time: 5.975803852081299
37 -1.898245e+03 2.567367e+00
* time: 6.14474081993103
38 -1.898246e+03 2.561623e+00
* time: 6.386682987213135
39 -1.898251e+03 2.530923e+00
* time: 6.5984838008880615
40 -1.898298e+03 2.674070e+00
* time: 6.765326023101807
41 -1.898300e+03 2.795248e+00
* time: 6.952793836593628
42 -1.898337e+03 3.729579e+00
* time: 7.193566799163818
43 -1.898428e+03 4.552315e+00
* time: 7.360801935195923
44 -1.898441e+03 4.064954e+00
* time: 7.566473007202148
45 -1.898444e+03 3.946564e+00
* time: 7.810468912124634
46 -1.898445e+03 3.887196e+00
* time: 8.041023015975952
47 -1.898447e+03 3.843752e+00
* time: 8.272588014602661
48 -1.898455e+03 1.529300e+02
* time: 8.512478828430176
49 -1.898544e+03 1.505372e+01
* time: 8.678820848464966
50 -1.898861e+03 3.981722e+00
* time: 8.842283010482788
51 -1.898892e+03 4.038144e+00
* time: 9.001397848129272
52 -1.898910e+03 4.249461e+00
* time: 9.144452810287476
53 -1.898976e+03 4.383672e+00
* time: 9.311044931411743
54 -1.899066e+03 6.170489e+00
* time: 9.469858884811401
55 -1.899168e+03 5.312543e+00
* time: 9.630230903625488
56 -1.899211e+03 2.138479e+00
* time: 9.79216194152832
57 -1.899217e+03 3.130196e-01
* time: 9.952345848083496
58 -1.899218e+03 4.118967e-02
* time: 10.102561950683594
59 -1.899218e+03 1.213456e-02
* time: 10.235992908477783
60 -1.899218e+03 3.026041e-03
* time: 10.377192974090576
61 -1.899218e+03 3.725927e-04
* time: 10.51627492904663FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Log-likelihood value: 1899.2177
Number of subjects: 120
Number of parameters: Fixed Optimized
0 11
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
------------------
Estimate
------------------
tvcl 2.6384
tvv 11.36
tvvp 8.1963
tvq 1.3182
tvka 4.8575
Ω₁,₁ 0.12921
Ω₂,₂ 0.06038
Ω₃,₃ 0.40714
Ω₄,₄ 0.14066
Ω₅,₅ 0.25355
σ_p 0.04881
------------------compare_estimates(; pkfit_2cmp, pkfit_2cmp_unfix_ka)11×3 DataFrame
| Row | parameter | pkfit_2cmp | pkfit_2cmp_unfix_ka |
|---|---|---|---|
| String | Float64? | Float64? | |
| 1 | tvcl | 2.81378 | 2.63839 |
| 2 | tvv | 11.0046 | 11.3604 |
| 3 | tvvp | 5.53997 | 8.19634 |
| 4 | tvq | 1.51591 | 1.31818 |
| 5 | tvka | 2.0 | 4.8575 |
| 6 | Ω₁,₁ | 0.102669 | 0.129205 |
| 7 | Ω₂,₂ | 0.0607756 | 0.0603797 |
| 8 | Ω₃,₃ | 1.20116 | 0.407139 |
| 9 | Ω₄,₄ | 0.423493 | 0.140658 |
| 10 | Ω₅,₅ | 0.244732 | 0.253546 |
| 11 | σ_p | 0.0484049 | 0.0488095 |
Let’s revaluate the goodness of fits and η distribution plots.
Not much change in the general gof plots
res_inspect_2cmp_unfix_ka = inspect(pkfit_2cmp_unfix_ka)[ 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()goodness_of_fit(res_inspect_2cmp_unfix_ka; figure = (; fontsize = 12))
But you can see a huge improvement in the ηka, (η₃) distribution which is now centered around zero
empirical_bayes_vs_covariates(
res_inspect_2cmp_unfix_ka;
categorical = [:Dose],
ebes = [:η₃],
figure = (; resolution = (600, 800)),
)
Finally looking at some individual plots for the same subjects as earlier:
fig_subject_fits2 = subject_fits(
res_inspect_2cmp_unfix_ka;
separate = true,
paginate = true,
facet = (; combinelabels = true, linkyaxes = false),
figure = (; fontsize = 18),
axis = (; xlabel = "Time (hr)", ylabel = "CTMx Concentration (ng/mL)"),
)
fig_subject_fits2[6]
The randomly sampled individual fits don’t seem good in some individuals, but we can evaluate this via a vpc to see how to go about.
3.4 Visual Predictive Checks (VPC)
We can now perform a vpc to check. The default plots provide a 80% prediction interval and a 95% simulated CI (shaded area) around each of the quantiles
pk_vpc = vpc(
pkfit_2cmp_unfix_ka,
200;
observations = [:Conc],
stratify_by = [:Dose],
ensemblealg = EnsembleThreads(), # multi-threading
)[ Info: Continuous VPCVisual Predictive Check
Type of VPC: Continuous VPC
Simulated populations: 200
Subjects in data: 40
Stratification variable(s): [:Dose]
Confidence level: 0.95
VPC lines: quantiles ([0.1, 0.5, 0.9])vpc_plot(
pk_2cmp,
pk_vpc;
rows = 1,
columns = 3,
figure = (; resolution = (1400, 1000), fontsize = 22),
axis = (;
xlabel = "Time (hr)",
ylabel = "Observed/Predicted\n CTMx Concentration (ng/mL)",
),
facet = (; combinelabels = true),
)
The visual predictive check suggests that the model captures the data well across all dose levels.
4 Additional Help
If you have questions regarding this tutorial, please post them on our discourse site.