using Pumas
using PumasUtilities
using NCA
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 a column that specifies the dose. The dose is already there as :Dose, we only have to turn it into a number, so let’s add the route and do that:
@rtransform! pkpain_noplb_df begin
:route = "ev"
:Dose = parse(Int, chop(:Dose; tail = 3))
endNow, we map the data variables to the read_nca function that prepares the data for NCA analysis. Every row of a profile carries its dose here, which is the layout that read_nca expects; data where the dose sits on a single row per profile and is missing elsewhere, such as a Pumas simulation table, is read with read_nca_pumas instead.
pkpain_nca = read_nca(
pkpain_noplb_df;
id = :Subject,
time = :Time,
dose = :Dose,
concentration = :Conc,
group = [:Dose],
route = :route,
)1440×12 DataFrame
1415 rows omitted
| Row | Dose | id | profile | time | nominal_time | dose_time | nominal_dose_time | concentration | dose | route | exclusion | exclusion_mask |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Int64 | Int64 | Int64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | String | Missing | Bool | |
| 1 | 20 | 1 | 1 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 20.0 | ev | missing | false |
| 2 | 20 | 1 | 1 | 0.5 | 0.5 | 0.0 | 0.0 | 1.15578 | 20.0 | ev | missing | false |
| 3 | 20 | 1 | 1 | 1.0 | 1.0 | 0.0 | 0.0 | 1.37211 | 20.0 | ev | missing | false |
| 4 | 20 | 1 | 1 | 1.5 | 1.5 | 0.0 | 0.0 | 1.30058 | 20.0 | ev | missing | false |
| 5 | 20 | 1 | 1 | 2.0 | 2.0 | 0.0 | 0.0 | 1.19195 | 20.0 | ev | missing | false |
| 6 | 20 | 1 | 1 | 2.5 | 2.5 | 0.0 | 0.0 | 1.13602 | 20.0 | ev | missing | false |
| 7 | 20 | 1 | 1 | 3.0 | 3.0 | 0.0 | 0.0 | 0.873224 | 20.0 | ev | missing | false |
| 8 | 20 | 1 | 1 | 4.0 | 4.0 | 0.0 | 0.0 | 0.739963 | 20.0 | ev | missing | false |
| 9 | 20 | 1 | 1 | 5.0 | 5.0 | 0.0 | 0.0 | 0.600143 | 20.0 | ev | missing | false |
| 10 | 20 | 1 | 1 | 6.0 | 6.0 | 0.0 | 0.0 | 0.425624 | 20.0 | ev | missing | false |
| 11 | 20 | 1 | 1 | 7.0 | 7.0 | 0.0 | 0.0 | 0.363418 | 20.0 | ev | missing | false |
| 12 | 20 | 1 | 1 | 8.0 | 8.0 | 0.0 | 0.0 | 0.304177 | 20.0 | ev | missing | false |
| 13 | 80 | 2 | 2 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 80.0 | ev | missing | false |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| 1429 | 80 | 159 | 120 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 80.0 | ev | missing | false |
| 1430 | 80 | 159 | 120 | 0.5 | 0.5 | 0.0 | 0.0 | 5.52603 | 80.0 | ev | missing | false |
| 1431 | 80 | 159 | 120 | 1.0 | 1.0 | 0.0 | 0.0 | 7.12166 | 80.0 | ev | missing | false |
| 1432 | 80 | 159 | 120 | 1.5 | 1.5 | 0.0 | 0.0 | 6.89178 | 80.0 | ev | missing | false |
| 1433 | 80 | 159 | 120 | 2.0 | 2.0 | 0.0 | 0.0 | 6.42015 | 80.0 | ev | missing | false |
| 1434 | 80 | 159 | 120 | 2.5 | 2.5 | 0.0 | 0.0 | 5.70904 | 80.0 | ev | missing | false |
| 1435 | 80 | 159 | 120 | 3.0 | 3.0 | 0.0 | 0.0 | 5.21105 | 80.0 | ev | missing | false |
| 1436 | 80 | 159 | 120 | 4.0 | 4.0 | 0.0 | 0.0 | 4.27633 | 80.0 | ev | missing | false |
| 1437 | 80 | 159 | 120 | 5.0 | 5.0 | 0.0 | 0.0 | 3.64571 | 80.0 | ev | missing | false |
| 1438 | 80 | 159 | 120 | 6.0 | 6.0 | 0.0 | 0.0 | 3.1888 | 80.0 | ev | missing | false |
| 1439 | 80 | 159 | 120 | 7.0 | 7.0 | 0.0 | 0.0 | 2.86869 | 80.0 | ev | missing | false |
| 1440 | 80 | 159 | 120 | 8.0 | 8.0 | 0.0 | 0.0 | 2.33494 | 80.0 | ev | missing | false |
Now 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)"),
)
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, size = (800, 1000)),
color = "black",
linewidth = 3,
axis = (; xlabel = "Time (hr)", ylabel = "CTMX Concentration (μg/mL)"),
)
All NCA endpoints that apply to this data are now computed with the run_nca function, but later we will only extract a couple of key metrics of interest.
pk_nca = run_nca(pkpain_nca)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(
pkpain_nca,
paginate = true,
axis = (; xlabel = "Time (hr)", ylabel = "CTMX Concentration (μg/mL)"),
# Legend options
legend = (; position = :bottom),
)
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×11 DataFrame
| Row | Dose | parameters | numsamples | nonmissing | minimum | maximum | mean | std | geomean | geostd | geomeanCV |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Int64 | String | Int64 | Int64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | 20 | cmax | 40 | 40 | 0.933113 | 2.70061 | 1.47354 | 0.361599 | 1.43399 | 1.26304 | 23.674 |
| 2 | 20 | aucinf_obs | 40 | 40 | 2.77495 | 14.1139 | 6.3764 | 2.22274 | 6.01975 | 1.41356 | 35.6743 |
| 3 | 80 | cmax | 40 | 40 | 3.29685 | 8.47195 | 5.78671 | 1.31987 | 5.64125 | 1.25771 | 23.2336 |
| 4 | 80 | aucinf_obs | 40 | 40 | 13.7102 | 49.122 | 29.5023 | 8.69558 | 28.2973 | 1.34149 | 30.0237 |
| 5 | 5 | cmax | 40 | 40 | 0.190482 | 0.538566 | 0.356087 | 0.0883522 | 0.345132 | 1.29294 | 26.1215 |
| 6 | 5 | aucinf_obs | 40 | 40 | 0.9137 | 3.39947 | 1.59819 | 0.490386 | 1.53397 | 1.32989 | 29.0992 |
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.0077 0.97567 1.03973
────────────────────────────────────Here’s a visualization for the dose linearity using a power model for cmax:
power_model(dp; legend = (; position = :bottom))
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, :cmt and :amt columns in the dataset, where :amt holds the dose on the dosing rows and is missing elsewhere.
@rtransform! pkpain_noplb_df begin
:evid = :Time == 0 ? 1 : 0
:cmt = :Time == 0 ? 1 : 2
:cmt2 = 1 # for zero order absorption
:amt = :Time == 0 ? :Dose : missing
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 = "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
@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 ~ @. ProportionalNormal(cp, σ_p)
end
endPumasModel
Parameters: tvcl, tvv, tvka, Ω, σ_p
Random effects: η
Covariates:
Dynamical system variables: Depot, Central
Dynamical system type: Closed form
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, pkpain_noplb, init_params(pk_1cmp))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
pkparam_np = (; init_params(pk_1cmp)..., Ω = Diagonal(zeros(3)))
pkfit_np = fit(pk_1cmp, pkpain_noplb, pkparam_np, NaivePooled(); constantcoef = (:Ω,))[ 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.0511012077331543 1 7.651412e+02 1.939231e+03 * time: 3.4250881671905518 2 7.079893e+02 1.675715e+03 * time: 3.4265999794006348 3 2.791962e+02 1.782931e+03 * time: 3.4275450706481934 4 4.841092e+01 1.134072e+03 * time: 3.4288010597229004 5 -4.109061e+01 1.279327e+03 * time: 3.429651975631714 6 -1.062281e+02 3.565024e+02 * time: 3.4305059909820557 7 -1.298325e+02 2.245444e+02 * time: 3.431339979171753 8 -1.432572e+02 1.031979e+02 * time: 3.43216609954834 9 -1.449990e+02 5.226590e+01 * time: 3.4330291748046875 10 -1.462807e+02 4.410109e+01 * time: 3.4338581562042236 11 -1.488553e+02 1.089712e+02 * time: 3.4347779750823975 12 -1.528708e+02 1.191595e+02 * time: 3.4356470108032227 13 -1.565497e+02 4.858221e+01 * time: 3.436462163925171 14 -1.576249e+02 1.255252e+01 * time: 3.437309980392456 15 -1.578681e+02 2.886131e+01 * time: 3.438129186630249 16 -1.580766e+02 3.050114e+01 * time: 3.4389421939849854 17 -1.582541e+02 1.790554e+01 * time: 3.4398090839385986 18 -1.583290e+02 3.831173e+00 * time: 3.440667152404785 19 -1.583498e+02 2.133704e+00 * time: 3.441495180130005 20 -1.583602e+02 2.667446e+00 * time: 3.4423232078552246 21 -1.583666e+02 1.429682e+00 * time: 3.4431450366973877 22 -1.583691e+02 6.661105e-01 * time: 3.443970203399658 23 -1.583700e+02 1.121904e+00 * time: 3.444784164428711 24 -1.583705e+02 9.437259e-01 * time: 3.4456610679626465 25 -1.583708e+02 5.280225e-01 * time: 3.44649600982666 26 -1.583709e+02 1.853276e-01 * time: 3.447331190109253 27 -1.583710e+02 1.341436e-01 * time: 3.4481630325317383 28 -1.583710e+02 8.843462e-02 * time: 3.4489831924438477 29 -1.583710e+02 1.332545e-01 * time: 3.449794054031372 30 -1.583710e+02 5.096336e-02 * time: 3.45062518119812 31 -1.583710e+02 1.096148e-01 * time: 3.4518091678619385 32 -1.583710e+02 1.369709e-01 * time: 3.4526891708374023 33 -1.583710e+02 5.987000e-02 * time: 3.4535770416259766 34 -1.583710e+02 2.466081e-02 * time: 3.454425096511841 35 -1.583710e+02 2.778113e-02 * time: 3.4552412033081055 36 -1.583710e+02 1.821778e-02 * time: 3.456057071685791 37 -1.583710e+02 3.099380e-03 * time: 3.457205057144165 38 -1.583710e+02 1.157716e-02 * time: 3.4583911895751953 39 -1.583710e+02 5.585459e-02 * time: 3.459207057952881 40 -1.583710e+02 1.894124e-02 * time: 3.460024118423462 41 -1.583710e+02 1.328404e-02 * time: 3.4608261585235596 42 -1.583710e+02 2.380617e-02 * time: 3.4617021083831787 43 -1.583710e+02 9.069363e-03 * time: 3.462522029876709 44 -1.583710e+02 2.126942e-02 * time: 3.46333909034729 45 -1.583710e+02 1.156032e-03 * time: 3.4641740322113037 46 -1.583710e+02 7.539635e-03 * time: 3.4649910926818848 47 -1.583710e+02 1.669066e-02 * time: 3.4658100605010986 48 -1.583710e+02 4.231777e-03 * time: 3.466658115386963 49 -1.583710e+02 4.547109e-03 * time: 3.4675910472869873 50 -1.583710e+02 1.437776e-02 * time: 3.4684700965881348 51 -1.583710e+02 2.130662e-04 * time: 3.469325065612793
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 120
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
Number of parameters: Constant Optimized
3 4
Likelihood approximation: NaivePooled
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: 158.37103
------------------
Estimate
------------------
tvcl 3.0057
tvv 14.09
tvka 43.37
† Ω₁,₁ 0.0
† Ω₂,₂ 0.0
† Ω₃,₃ 0.0
σ_p 0.32999
------------------
† indicates constant parameterscoefficients_table(pkfit_np)7×4 DataFrame
| Row | Parameter | Description | Constant | Estimate |
|---|---|---|---|---|
| String | String | Bool | Float64 | |
| 1 | tvcl | Clearance (L/hr) | false | 3.006 |
| 2 | tvv | Volume (L) | false | 14.09 |
| 3 | tvka | Absorption rate constant (h-1) | false | 43.37 |
| 4 | Ω₁,₁ | ΩCL | true | 0.0 |
| 5 | Ω₂,₂ | ΩVc | true | 0.0 |
| 6 | Ω₃,₃ | ΩKa | true | 0.0 |
| 7 | σ_p | Proportional RUV | false | 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 = [loglikelihood(pk_1cmp, subj, pkparam, FOCE()) for subj in pkpain_noplb]
# the plot below is using native CairoMakie `hist`
hist(lls; bins = 10, normalization = :none, color = (:black, 0.5))
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())120-element Vector{@NamedTuple{id::String, nll::Float64}}:
(id = "148", nll = 16.659658856844796)
(id = "135", nll = 16.648985190076335)
(id = "156", nll = 15.959069556607494)
(id = "159", nll = 15.441218240496484)
(id = "149", nll = 14.715134644119509)
(id = "88", nll = 13.09709837464614)
(id = "16", nll = 12.982280521931443)
(id = "61", nll = 12.652182902303682)
(id = "71", nll = 12.500330088085484)
(id = "59", nll = 12.24151025480522)
⋮
(id = "57", nll = -22.797674232534312)
(id = "93", nll = -22.836900711478208)
(id = "12", nll = -23.007742339519247)
(id = "123", nll = -23.292751843079234)
(id = "41", nll = -23.425412534960515)
(id = "99", nll = -23.535214841901052)
(id = "29", nll = -24.025959868383097)
(id = "52", nll = -24.164757842493685)
(id = "24", nll = -25.57209232565844)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,))[ 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: 3.2901763916015625e-5 1 -7.019247e+02 1.748155e+02 * time: 6.636348009109497 2 -7.308960e+02 2.936116e+02 * time: 6.702507019042969 3 -9.540611e+02 1.972709e+02 * time: 7.076305866241455 4 -1.110598e+03 2.654265e+02 * time: 7.233352899551392 5 -1.200293e+03 9.910898e+01 * time: 7.284558057785034 6 -1.216709e+03 1.484081e+02 * time: 7.3378660678863525 7 -1.225139e+03 2.202360e+01 * time: 7.39548397064209 8 -1.226781e+03 1.766571e+01 * time: 7.453456878662109 9 -1.228206e+03 2.073312e+01 * time: 7.610852956771851 10 -1.229694e+03 2.425736e+01 * time: 7.659614086151123 11 -1.230765e+03 1.519853e+01 * time: 7.708623886108398 12 -1.231081e+03 4.521041e+00 * time: 7.7539238929748535 13 -1.231178e+03 1.723965e+00 * time: 7.800030946731567 14 -1.231188e+03 9.412196e-01 * time: 7.841705083847046 15 -1.231188e+03 2.205522e-01 * time: 7.883446931838989 16 -1.231188e+03 5.487003e-02 * time: 7.962429046630859 17 -1.231188e+03 4.802975e-02 * time: 7.992333889007568 18 -1.231188e+03 9.339034e-02 * time: 8.024332046508789 19 -1.231188e+03 5.586507e-02 * time: 8.051910877227783 20 -1.231188e+03 1.076107e-02 * time: 8.08040189743042 21 -1.231188e+03 7.633227e-04 * time: 8.105845928192139
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 120
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
Number of parameters: Constant Optimized
1 6
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: 1231.188
-------------------
Estimate
-------------------
tvcl 3.1642
tvv 13.288
† tvka 2.0
Ω₁,₁ 0.08494
Ω₂,₂ 0.048569
Ω₃,₃ 5.581
σ_p 0.10093
-------------------
† indicates constant parametersinfer(pkfit_1cmp)[ Info: Calculating: variance-covariance matrix. [ Info: Done.
Asymptotic inference results using sandwich estimator
Dynamical system type: Closed form
Number of subjects: 120
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
Number of parameters: Constant Optimized
1 6
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: 1231.188
---------------------------------------------------------
Estimate SE 95.0% C.I.
---------------------------------------------------------
tvcl 3.1642 0.08662 [ 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.048569 0.0063501 [ 0.036123; 0.061015]
Ω₃,₃ 5.581 1.2188 [ 3.1922 ; 7.9699 ]
σ_p 0.10093 0.0057196 [ 0.089718; 0.11214 ]
---------------------------------------------------------
† indicates constant parametersNotice 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
@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 ~ @. ProportionalNormal(cp, σ_p)
end
endPumasModel
Parameters: tvcl, tvv, tvvp, tvq, tvka, Ω, σ_p
Random effects: η
Covariates:
Dynamical system variables: Depot, Central, Peripheral
Dynamical system type: Closed form
Derived: Conc
Observed: Conc3.2.3.1 FOCE
pkfit_2cmp = fit(
pk_2cmp,
pkpain_noplb,
(; init_params(pk_2cmp)..., tvka = 2),
FOCE();
constantcoef = (:tvka,),
)[ 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: 4.00543212890625e-5 1 -9.197817e+02 9.927951e+02 * time: 11.594305992126465 2 -1.372640e+03 2.054986e+02 * time: 11.693562984466553 3 -1.446326e+03 1.543987e+02 * time: 11.794328927993774 4 -1.545570e+03 1.855028e+02 * time: 11.890448093414307 5 -1.581346e+03 1.625556e+02 * time: 12.165079116821289 6 -1.639868e+03 1.274103e+02 * time: 12.256431102752686 7 -1.695755e+03 9.009168e+01 * time: 12.337693929672241 8 -1.722327e+03 5.933190e+01 * time: 12.421711921691895 9 -1.736765e+03 7.225508e+01 * time: 12.547975063323975 10 -1.753684e+03 7.288038e+01 * time: 12.630341053009033 11 -1.764130e+03 6.011627e+01 * time: 12.717697143554688 12 -1.779290e+03 4.959305e+01 * time: 12.82308292388916 13 -1.791334e+03 4.832595e+01 * time: 12.96094298362732 14 -1.800062e+03 2.768105e+01 * time: 13.07219409942627 15 -1.810796e+03 2.472791e+01 * time: 13.20137095451355 16 -1.812985e+03 2.690187e+01 * time: 13.292884111404419 17 -1.817454e+03 5.103483e+01 * time: 13.41114592552185 18 -1.820664e+03 3.560949e+01 * time: 13.502815961837769 19 -1.823411e+03 1.709308e+01 * time: 13.619204998016357 20 -1.824134e+03 1.659633e+01 * time: 13.702304124832153 21 -1.824937e+03 1.543537e+01 * time: 13.811867952346802 22 -1.825374e+03 1.842823e+01 * time: 13.90004301071167 23 -1.825634e+03 5.318790e+00 * time: 13.989026069641113 24 -1.825811e+03 6.286861e+00 * time: 14.113507986068726 25 -1.825982e+03 6.154970e+00 * time: 14.199043035507202 26 -1.826330e+03 9.237152e+00 * time: 14.313848972320557 27 -1.826668e+03 8.860375e+00 * time: 14.402260065078735 28 -1.826861e+03 5.853416e+00 * time: 14.5143301486969 29 -1.827033e+03 4.447191e+00 * time: 14.601380109786987 30 -1.827093e+03 4.010670e+00 * time: 14.707279920578003 31 -1.827256e+03 5.055045e+00 * time: 14.797250986099243 32 -1.827365e+03 4.933784e+00 * time: 14.908875942230225 33 -1.827475e+03 3.602187e+00 * time: 14.993454933166504 34 -1.827518e+03 1.391397e+00 * time: 15.074494123458862 35 -1.827529e+03 8.330800e-01 * time: 15.180057048797607 36 -1.827533e+03 2.767235e-01 * time: 15.258908987045288 37 -1.827533e+03 2.613588e-01 * time: 15.354691982269287 38 -1.827533e+03 2.710603e-01 * time: 15.426382064819336 39 -1.827535e+03 4.294332e-01 * time: 15.535665035247803 40 -1.827535e+03 3.836718e-01 * time: 15.62889313697815 41 -1.827536e+03 1.534356e-01 * time: 15.7155179977417 42 -1.827536e+03 2.708580e-02 * time: 15.819788932800293 43 -1.827536e+03 1.834493e-03 * time: 15.93691611289978 44 -1.827536e+03 6.451297e-04 * time: 16.050787925720215
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 120
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
Number of parameters: Constant Optimized
1 10
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: 1827.5359
-------------------
Estimate
-------------------
tvcl 2.8451
tvv 10.988
tvvp 5.0546
tvq 1.5282
† tvka 2.0
Ω₁,₁ 0.10042
Ω₂,₂ 0.061795
Ω₃,₃ 1.1334
Ω₄,₄ 0.46411
Ω₅,₅ 0.24145
σ_p 0.048443
-------------------
† indicates constant parameters3.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.84515 |
| 2 | tvv | 13.288 | 10.9884 |
| 3 | tvka | 2.0 | 2.0 |
| 4 | Ω₁,₁ | 0.0849404 | 0.100422 |
| 5 | Ω₂,₂ | 0.0485685 | 0.0617954 |
| 6 | Ω₃,₃ | 5.58105 | 1.13336 |
| 7 | σ_p | 0.100928 | 0.0484426 |
| 8 | tvvp | missing | 5.05458 |
| 9 | tvq | missing | 1.52817 |
| 10 | Ω₄,₄ | missing | 0.464107 |
| 11 | Ω₅,₅ | missing | 0.241455 |
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: 1190.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)
end19×3 DataFrame
| Row | Metric | pk2cmp | pk1cmp |
|---|---|---|---|
| String | Any | Any | |
| 1 | Estimation Time | 16.051 | 8.106 |
| 2 | Subjects | 120 | 120 |
| 3 | Fixed Parameters | 1 | 1 |
| 4 | Optimized Parameters | 10 | 6 |
| 5 | Conc Active Observations | 1320 | 1320 |
| 6 | Conc Missing Observations | 0 | 0 |
| 7 | Total Active Observations | 1320 | 1320 |
| 8 | Total Missing Observations | 0 | 0 |
| 9 | Likelihood Approximation | FOCE | FOCE |
| 10 | LogLikelihood (LL) | 1827.54 | 1231.19 |
| 11 | -2LL | -3655.07 | -2462.38 |
| 12 | AIC | -3635.07 | -2450.38 |
| 13 | BIC | -3583.22 | -2419.26 |
| 14 | (η-shrinkage) η₁ | 0.038 | 0.016 |
| 15 | (η-shrinkage) η₂ | 0.047 | 0.04 |
| 16 | (η-shrinkage) η₃ | 0.511 | 0.733 |
| 17 | (ϵ-shrinkage) Conc | 0.185 | 0.105 |
| 18 | (η-shrinkage) η₄ | 0.266 | missing |
| 19 | (η-shrinkage) η₅ | 0.19 | 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 dose control parameters. [ Info: Evaluating individual parameters. [ Info: Done.
FittedPumasModelInspection
Likelihood approximation used for weighted residuals: FOCEres_inspect_2cmp = inspect(pkfit_2cmp)[ 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: FOCEgof_1cmp = goodness_of_fit(
res_inspect_1cmp;
figure = (; fontsize = 12),
legend = (; position = :bottom),
)
gof_2cmp = goodness_of_fit(
res_inspect_2cmp;
figure = (; fontsize = 12),
legend = (; position = :bottom),
)
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,
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 = (; size = (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: 1.9073486328125e-5 1 -8.682982e+02 1.000199e+03 * time: 1.5660550594329834 2 -1.381870e+03 5.008081e+02 * time: 1.6824331283569336 3 -1.551053e+03 6.833490e+02 * time: 1.825402021408081 4 -1.680887e+03 1.834586e+02 * time: 4.601530075073242 5 -1.726118e+03 8.870274e+01 * time: 4.754997968673706 6 -1.761023e+03 1.162036e+02 * time: 4.855561971664429 7 -1.786619e+03 1.114552e+02 * time: 4.956387996673584 8 -1.863556e+03 9.914305e+01 * time: 5.0818421840667725 9 -1.882942e+03 5.342676e+01 * time: 5.203715085983276 10 -1.888020e+03 2.010181e+01 * time: 5.330049991607666 11 -1.889832e+03 1.867263e+01 * time: 5.5873260498046875 12 -1.891649e+03 1.668512e+01 * time: 5.690200090408325 13 -1.892615e+03 1.820701e+01 * time: 5.791347026824951 14 -1.893453e+03 1.745195e+01 * time: 5.890608072280884 15 -1.894760e+03 1.850174e+01 * time: 5.998337984085083 16 -1.895647e+03 1.773939e+01 * time: 6.220883131027222 17 -1.896597e+03 1.143462e+01 * time: 6.328296184539795 18 -1.897114e+03 9.720097e+00 * time: 6.442536115646362 19 -1.897373e+03 6.054321e+00 * time: 6.611508131027222 20 -1.897498e+03 3.985954e+00 * time: 6.714517116546631 21 -1.897571e+03 4.262464e+00 * time: 6.811750173568726 22 -1.897633e+03 4.010234e+00 * time: 6.921353101730347 23 -1.897714e+03 4.805376e+00 * time: 7.062586069107056 24 -1.897802e+03 3.508705e+00 * time: 7.1640331745147705 25 -1.897865e+03 3.691477e+00 * time: 7.30308198928833 26 -1.897900e+03 2.982719e+00 * time: 7.399280071258545 27 -1.897928e+03 2.563789e+00 * time: 7.526124954223633 28 -1.897968e+03 3.261484e+00 * time: 7.622564077377319 29 -1.898013e+03 3.064689e+00 * time: 7.7281270027160645 30 -1.898040e+03 1.636525e+00 * time: 7.885051965713501 31 -1.898051e+03 1.439997e+00 * time: 7.982959985733032 32 -1.898057e+03 1.436504e+00 * time: 8.11266303062439 33 -1.898069e+03 1.881528e+00 * time: 8.209388017654419 34 -1.898095e+03 3.253165e+00 * time: 8.340306043624878 35 -1.898142e+03 4.257941e+00 * time: 8.439387083053589 36 -1.898199e+03 3.685241e+00 * time: 8.572865962982178 37 -1.898245e+03 2.567364e+00 * time: 8.676421165466309 38 -1.898246e+03 2.525090e+00 * time: 9.069570064544678 39 -1.898277e+03 2.322951e+00 * time: 9.211320161819458 40 -1.898288e+03 2.250254e+00 * time: 9.391138076782227 41 -1.898389e+03 4.741678e+00 * time: 9.555300951004028 42 -1.898399e+03 4.764969e+00 * time: 9.83967113494873 43 -1.898458e+03 3.847089e+00 * time: 9.938238143920898 44 -1.898479e+03 2.349070e+00 * time: 10.056939125061035 45 -1.898479e+03 1.423477e+00 * time: 10.303730964660645 46 -1.898478e+03 5.093603e+00 * time: 14.122579097747803
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 120
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
Number of parameters: Constant Optimized
0 11
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: NoXChange
Log-likelihood value: 1898.4776
-----------------
Estimate
-----------------
tvcl 2.6193
tvv 11.379
tvvp 8.4617
tvq 1.3162
tvka 4.8941
Ω₁,₁ 0.13246
Ω₂,₂ 0.059667
Ω₃,₃ 0.41569
Ω₄,₄ 0.080891
Ω₅,₅ 0.24988
σ_p 0.049098
-----------------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.84515 | 2.61934 |
| 2 | tvv | 10.9884 | 11.3791 |
| 3 | tvvp | 5.05458 | 8.46174 |
| 4 | tvq | 1.52817 | 1.31618 |
| 5 | tvka | 2.0 | 4.89405 |
| 6 | Ω₁,₁ | 0.100422 | 0.132462 |
| 7 | Ω₂,₂ | 0.0617954 | 0.0596667 |
| 8 | Ω₃,₃ | 1.13336 | 0.415689 |
| 9 | Ω₄,₄ | 0.464107 | 0.0808911 |
| 10 | Ω₅,₅ | 0.241455 | 0.249884 |
| 11 | σ_p | 0.0484426 | 0.0490982 |
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 dose control parameters. [ Info: Evaluating individual parameters. [ Info: Done.
FittedPumasModelInspection
Likelihood approximation used for weighted residuals: FOCEgoodness_of_fit(
res_inspect_2cmp_unfix_ka;
figure = (; fontsize = 12),
legend = (; position = :bottom),
)
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 = (; size = (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 = (; 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])[ Info: Continuous VPC
Visual Predictive Check
Type of VPC: Continuous VPC
Simulated populations: 200
Subjects in data: 120
Stratification variable(s): [:Dose]
Confidence level: 0.95
VPC lines: quantiles ([0.1, 0.5, 0.9])vpc_plot(
pk_2cmp,
pk_vpc;
figure = (; size = (1400, 1000), fontsize = 22),
axis = (;
xlabel = "Time (hr)",
ylabel = "Observed/Predicted\n CTMx Concentration (ng/mL)",
),
)
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.