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, size = (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
end┌ Warning: Covariate Dose is not used in the model.
└ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/aZRyj/src/dsl/model_macro.jl:2856PumasModel
Parameters: tvcl, tvv, tvka, Ω, σ_p
Random effects: η
Covariates: Dose
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, 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.024767160415649414
1 2.343899e+02 1.747348e+03
* time: 0.8609800338745117
2 9.696232e+01 1.198088e+03
* time: 0.8638792037963867
3 -7.818699e+01 5.538151e+02
* time: 0.8661110401153564
4 -1.234803e+02 2.462514e+02
* time: 0.8687350749969482
5 -1.372888e+02 2.067458e+02
* time: 0.8711390495300293
6 -1.410579e+02 1.162950e+02
* time: 0.8735032081604004
7 -1.434754e+02 5.632816e+01
* time: 0.8761460781097412
8 -1.453401e+02 7.859270e+01
* time: 0.8787841796875
9 -1.498185e+02 1.455606e+02
* time: 0.8812761306762695
10 -1.534371e+02 1.303682e+02
* time: 0.883613109588623
11 -1.563557e+02 5.975474e+01
* time: 0.8858611583709717
12 -1.575052e+02 9.308611e+00
* time: 0.8879201412200928
13 -1.579357e+02 1.234484e+01
* time: 0.8903341293334961
14 -1.581874e+02 7.478196e+00
* time: 0.8921821117401123
15 -1.582981e+02 2.027162e+00
* time: 0.8944389820098877
16 -1.583375e+02 5.578262e+00
* time: 0.8962001800537109
17 -1.583556e+02 4.727050e+00
* time: 0.8984131813049316
18 -1.583644e+02 2.340173e+00
* time: 0.9006831645965576
19 -1.583680e+02 7.738100e-01
* time: 0.9024100303649902
20 -1.583696e+02 3.300689e-01
* time: 0.9047331809997559
21 -1.583704e+02 3.641985e-01
* time: 0.906508207321167
22 -1.583707e+02 4.365901e-01
* time: 0.908811092376709
23 -1.583709e+02 3.887800e-01
* time: 0.9106271266937256
24 -1.583710e+02 2.766977e-01
* time: 0.9129171371459961
25 -1.583710e+02 1.758029e-01
* time: 0.9152061939239502
26 -1.583710e+02 1.133947e-01
* time: 0.9169390201568604
27 -1.583710e+02 7.922544e-02
* time: 0.9192841053009033
28 -1.583710e+02 5.954998e-02
* time: 0.9210011959075928
29 -1.583710e+02 4.157079e-02
* time: 0.9232480525970459
30 -1.583710e+02 4.295447e-02
* time: 0.9254801273345947
31 -1.583710e+02 5.170753e-02
* time: 0.9272100925445557
32 -1.583710e+02 2.644383e-02
* time: 0.9301230907440186
33 -1.583710e+02 4.548993e-03
* time: 0.932481050491333
34 -1.583710e+02 2.501804e-02
* time: 0.9355301856994629
35 -1.583710e+02 3.763440e-02
* time: 0.9377779960632324
36 -1.583710e+02 3.206026e-02
* time: 0.9394772052764893
37 -1.583710e+02 1.003698e-02
* time: 0.9418141841888428
38 -1.583710e+02 2.209089e-02
* time: 0.9442269802093506
39 -1.583710e+02 4.954172e-03
* time: 0.9461190700531006
40 -1.583710e+02 1.609373e-02
* time: 0.9491031169891357
41 -1.583710e+02 1.579802e-02
* time: 0.9508721828460693
42 -1.583710e+02 1.014113e-03
* time: 0.9531991481781006
43 -1.583710e+02 6.050644e-03
* time: 0.9561731815338135
44 -1.583710e+02 1.354412e-02
* time: 0.9579191207885742
45 -1.583710e+02 4.473248e-03
* time: 0.960237979888916
46 -1.583710e+02 4.644735e-03
* time: 0.9619970321655273
47 -1.583710e+02 9.829910e-03
* time: 0.9643549919128418
48 -1.583710e+02 1.047561e-03
* time: 0.966576099395752
49 -1.583710e+02 8.366895e-03
* time: 0.9683160781860352
50 -1.583710e+02 7.879055e-04
* time: 0.9706981182098389FittedPumasModel
Successful minimization: true
Likelihood approximation: NaivePooled
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
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))
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.65965885684477)
(id = "135", nll = 16.648985190076335)
(id = "156", nll = 15.959069556607496)
(id = "159", nll = 15.441218240496484)
(id = "149", nll = 14.71513464411951)
(id = "88", nll = 13.09709837464614)
(id = "16", nll = 12.98228052193144)
(id = "61", nll = 12.652182902303679)
(id = "71", nll = 12.500330088085505)
(id = "59", nll = 12.241510254805235)
⋮
(id = "57", nll = -22.79767423253431)
(id = "93", nll = -22.836900711478208)
(id = "12", nll = -23.007742339519247)
(id = "123", nll = -23.292751843079234)
(id = "41", nll = -23.425412534960515)
(id = "99", nll = -23.535214841901112)
(id = "29", nll = -24.025959868383083)
(id = "52", nll = -24.164757842493685)
(id = "24", nll = -25.57209232565845)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.318092346191406e-5
1 -7.022088e+02 1.707063e+02
* time: 0.2911660671234131
2 -7.314067e+02 2.903269e+02
* time: 0.4536252021789551
3 -8.520591e+02 2.285888e+02
* time: 0.609976053237915
4 -1.120191e+03 3.795410e+02
* time: 0.9874510765075684
5 -1.178784e+03 2.323978e+02
* time: 1.1393101215362549
6 -1.218320e+03 9.699907e+01
* time: 1.3297491073608398
7 -1.223641e+03 5.862105e+01
* time: 1.452510118484497
8 -1.227620e+03 1.831403e+01
* time: 1.5801920890808105
9 -1.228381e+03 2.132323e+01
* time: 1.708693027496338
10 -1.230098e+03 2.921228e+01
* time: 1.8575990200042725
11 -1.230854e+03 2.029662e+01
* time: 1.9839379787445068
12 -1.231116e+03 5.229097e+00
* time: 2.1377670764923096
13 -1.231179e+03 1.689232e+00
* time: 2.3081750869750977
14 -1.231187e+03 1.215379e+00
* time: 2.432142972946167
15 -1.231188e+03 2.770380e-01
* time: 2.542647123336792
16 -1.231188e+03 1.636653e-01
* time: 2.6455509662628174
17 -1.231188e+03 2.701133e-01
* time: 2.7702670097351074
18 -1.231188e+03 3.163363e-01
* time: 2.862377166748047
19 -1.231188e+03 1.505149e-01
* time: 2.9520931243896484
20 -1.231188e+03 2.484999e-02
* time: 3.0393831729888916
21 -1.231188e+03 8.446863e-04
* time: 3.1454319953918457FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
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
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
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.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.048568 0.0063502 [ 0.036122; 0.061014]
Ω₃,₃ 5.5811 1.2188 [ 3.1922 ; 7.97 ]
σ_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
end┌ Warning: Covariate Dose is not used in the model.
└ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/aZRyj/src/dsl/model_macro.jl:2856PumasModel
Parameters: tvcl, tvv, tvvp, tvq, tvka, Ω, σ_p
Random effects: η
Covariates: Dose
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), 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: 7.390975952148438e-5
1 -9.197817e+02 9.927951e+02
* time: 0.28101396560668945
2 -1.372640e+03 2.054986e+02
* time: 0.5808839797973633
3 -1.446326e+03 1.543987e+02
* time: 0.8994128704071045
4 -1.545570e+03 1.855028e+02
* time: 1.171355962753296
5 -1.581449e+03 1.713157e+02
* time: 1.6162738800048828
6 -1.639433e+03 1.257382e+02
* time: 1.9116718769073486
7 -1.695964e+03 7.450539e+01
* time: 2.1721668243408203
8 -1.722243e+03 5.961044e+01
* time: 2.4609179496765137
9 -1.736883e+03 7.320921e+01
* time: 2.762345790863037
10 -1.753547e+03 7.501938e+01
* time: 3.022596836090088
11 -1.764053e+03 6.185661e+01
* time: 3.3250157833099365
12 -1.778991e+03 4.831033e+01
* time: 3.662087917327881
13 -1.791492e+03 4.943278e+01
* time: 4.002189874649048
14 -1.799847e+03 2.871410e+01
* time: 4.3572258949279785
15 -1.805374e+03 7.520791e+01
* time: 4.762135982513428
16 -1.816260e+03 2.990621e+01
* time: 5.0929179191589355
17 -1.818252e+03 2.401915e+01
* time: 5.4055140018463135
18 -1.822988e+03 2.587225e+01
* time: 5.717067003250122
19 -1.824653e+03 1.550517e+01
* time: 5.969350814819336
20 -1.826074e+03 1.788927e+01
* time: 6.251224994659424
21 -1.826821e+03 1.888389e+01
* time: 6.516048908233643
22 -1.827900e+03 1.432840e+01
* time: 6.797647953033447
23 -1.828511e+03 9.422041e+00
* time: 7.093999862670898
24 -1.828754e+03 5.363442e+00
* time: 7.37075400352478
25 -1.828862e+03 4.916159e+00
* time: 7.6524498462677
26 -1.829007e+03 4.695755e+00
* time: 7.949594974517822
27 -1.829358e+03 1.090249e+01
* time: 8.220670938491821
28 -1.829830e+03 1.451325e+01
* time: 8.522029876708984
29 -1.830201e+03 1.108715e+01
* time: 8.828874826431274
30 -1.830360e+03 2.891223e+00
* time: 9.097190856933594
31 -1.830390e+03 1.695557e+00
* time: 9.576762914657593
32 -1.830404e+03 1.601712e+00
* time: 10.05839991569519
33 -1.830432e+03 2.823385e+00
* time: 10.373932838439941
34 -1.830477e+03 4.060617e+00
* time: 10.721182823181152
35 -1.830528e+03 5.133499e+00
* time: 11.046433925628662
36 -1.830593e+03 2.830970e+00
* time: 11.394709825515747
37 -1.830616e+03 3.342835e+00
* time: 11.726136922836304
38 -1.830622e+03 3.708884e+00
* time: 12.01149296760559
39 -1.830625e+03 2.062934e+00
* time: 12.29072380065918
40 -1.830627e+03 1.278569e+00
* time: 12.56161093711853
41 -1.830628e+03 1.832895e+00
* time: 12.800863981246948
42 -1.830628e+03 3.768840e-01
* time: 13.065598011016846
43 -1.830629e+03 3.152895e-01
* time: 13.292333841323853
44 -1.830630e+03 4.871060e-01
* time: 13.564238786697388
45 -1.830630e+03 3.110627e-01
* time: 13.843886852264404
46 -1.830630e+03 2.687758e-02
* time: 14.067235946655273
47 -1.830630e+03 4.694018e-03
* time: 14.304260969161987
48 -1.830630e+03 8.272969e-03
* time: 14.48662281036377
49 -1.830630e+03 8.249151e-03
* time: 14.748056888580322
50 -1.830630e+03 8.245562e-03
* time: 15.0417799949646
51 -1.830630e+03 8.240030e-03
* time: 15.371457815170288
52 -1.830630e+03 8.240030e-03
* time: 15.729442834854126
53 -1.830630e+03 8.240030e-03
* time: 16.083330869674683FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: 1830.6304
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.4235
Ω₅,₅ 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.0607757 |
| 6 | Ω₃,₃ | 5.58107 | 1.20115 |
| 7 | σ_p | 0.100928 | 0.0484049 |
| 8 | tvvp | missing | 5.53998 |
| 9 | tvq | missing | 1.51591 |
| 10 | Ω₄,₄ | missing | 0.423495 |
| 11 | Ω₅,₅ | missing | 0.244731 |
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 | 16.084 | 3.146 |
| 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{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}} | Pumas.FOCE{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}} |
| 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{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}(Optim.NewtonTrustRegion{Float64}(1.0, 100.0, 1.4901161193847656e-8, 0.1, 0.25, 0.75, false), Optim.Options(x_abstol = 0.0, x_reltol = 0.0, f_abstol = 0.0, f_reltol = 0.0, g_abstol = 1.0e-5, g_reltol = 1.0e-8, outer_x_abstol = 0.0, outer_x_reltol = 0.0, outer_f_abstol = 0.0, outer_f_reltol = 0.0, outer_g_abstol = 1.0e-8, outer_g_reltol = 1.0e-8, f_calls_limit = 0, g_calls_limit = 0, h_calls_limit = 0, allow_f_increases = false, allow_outer_f_increases = true, successive_f_tol = 1, iterations = 1000, outer_iterations = 1000, store_trace = false, trace_simplex = false, show_trace = false, extended_trace = false, show_warnings = true, show_every = 1, time_limit = NaN, )
)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{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}(Optim.NewtonTrustRegion{Float64}(1.0, 100.0, 1.4901161193847656e-8, 0.1, 0.25, 0.75, false), Optim.Options(x_abstol = 0.0, x_reltol = 0.0, f_abstol = 0.0, f_reltol = 0.0, g_abstol = 1.0e-5, g_reltol = 1.0e-8, outer_x_abstol = 0.0, outer_x_reltol = 0.0, outer_f_abstol = 0.0, outer_f_reltol = 0.0, outer_g_abstol = 1.0e-8, outer_g_reltol = 1.0e-8, f_calls_limit = 0, g_calls_limit = 0, h_calls_limit = 0, allow_f_increases = false, allow_outer_f_increases = true, successive_f_tol = 1, iterations = 1000, outer_iterations = 1000, store_trace = false, trace_simplex = false, show_trace = false, extended_trace = false, show_warnings = true, show_every = 1, time_limit = NaN, )
)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 = (; 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: 8.296966552734375e-5
1 -8.682982e+02 1.000199e+03
* time: 0.3516368865966797
2 -1.381870e+03 5.008081e+02
* time: 0.6976919174194336
3 -1.551053e+03 6.833490e+02
* time: 1.0851240158081055
4 -1.680887e+03 1.834586e+02
* time: 1.3713510036468506
5 -1.726118e+03 8.870274e+01
* time: 1.7111878395080566
6 -1.761023e+03 1.162036e+02
* time: 2.0172178745269775
7 -1.786619e+03 1.114552e+02
* time: 2.3584039211273193
8 -1.863556e+03 9.914305e+01
* time: 2.688314914703369
9 -1.882942e+03 5.342676e+01
* time: 3.038245916366577
10 -1.888020e+03 2.010181e+01
* time: 3.403244972229004
11 -1.889832e+03 1.867263e+01
* time: 3.7195088863372803
12 -1.891649e+03 1.668512e+01
* time: 4.078242063522339
13 -1.892615e+03 1.820701e+01
* time: 4.398313999176025
14 -1.893453e+03 1.745195e+01
* time: 4.741631031036377
15 -1.894760e+03 1.850174e+01
* time: 5.061720848083496
16 -1.895647e+03 1.773939e+01
* time: 5.399042844772339
17 -1.896597e+03 1.143462e+01
* time: 5.762431859970093
18 -1.897114e+03 9.720097e+00
* time: 6.065305948257446
19 -1.897373e+03 6.054321e+00
* time: 6.41254186630249
20 -1.897498e+03 3.985954e+00
* time: 6.727407932281494
21 -1.897571e+03 4.262464e+00
* time: 7.0630738735198975
22 -1.897633e+03 4.010234e+00
* time: 7.372957944869995
23 -1.897714e+03 4.805375e+00
* time: 7.708884000778198
24 -1.897802e+03 3.508706e+00
* time: 8.086956024169922
25 -1.897865e+03 3.691477e+00
* time: 8.379258871078491
26 -1.897900e+03 2.982720e+00
* time: 8.713791847229004
27 -1.897928e+03 2.563790e+00
* time: 9.012503862380981
28 -1.897968e+03 3.261485e+00
* time: 9.343833923339844
29 -1.898013e+03 3.064690e+00
* time: 9.640743017196655
30 -1.898040e+03 1.636525e+00
* time: 9.96767783164978
31 -1.898051e+03 1.439997e+00
* time: 10.278594017028809
32 -1.898057e+03 1.436504e+00
* time: 10.59801197052002
33 -1.898069e+03 1.881529e+00
* time: 10.936486959457397
34 -1.898095e+03 3.253165e+00
* time: 11.224890947341919
35 -1.898142e+03 4.257942e+00
* time: 11.564955949783325
36 -1.898199e+03 3.685241e+00
* time: 11.869658946990967
37 -1.898245e+03 2.567364e+00
* time: 12.213368892669678
38 -1.898246e+03 2.561591e+00
* time: 12.715200901031494
39 -1.898251e+03 2.530888e+00
* time: 13.108697891235352
40 -1.898298e+03 2.673696e+00
* time: 13.453145027160645
41 -1.898300e+03 2.794639e+00
* time: 13.888885974884033
42 -1.898337e+03 3.751590e+00
* time: 14.30534291267395
43 -1.898421e+03 4.878407e+00
* time: 14.6543869972229
44 -1.898433e+03 4.391719e+00
* time: 15.042754888534546
45 -1.898437e+03 4.216518e+00
* time: 15.531946897506714
46 -1.898442e+03 4.108397e+00
* time: 16.035634994506836
47 -1.898446e+03 3.934902e+00
* time: 16.539177894592285
48 -1.898449e+03 3.769838e+00
* time: 17.043819904327393
49 -1.898450e+03 3.739486e+00
* time: 17.480304956436157
50 -1.898450e+03 3.712049e+00
* time: 17.99446988105774
51 -1.898457e+03 3.623436e+00
* time: 18.428550004959106
52 -1.898471e+03 2.668312e+00
* time: 18.744637966156006
53 -1.898479e+03 2.302438e+00
* time: 19.086650848388672
54 -1.898480e+03 2.386566e-01
* time: 19.396994829177856
55 -1.898480e+03 7.802040e-01
* time: 19.720026969909668
56 -1.898480e+03 7.369786e-01
* time: 20.181758880615234
57 -1.898480e+03 5.113191e-01
* time: 20.55393695831299
58 -1.898480e+03 3.067709e-01
* time: 20.868165969848633
59 -1.898480e+03 3.076791e-01
* time: 21.207254886627197
60 -1.898480e+03 3.102066e-01
* time: 21.47196102142334
61 -1.898480e+03 3.102066e-01
* time: 21.881314039230347
62 -1.898480e+03 3.102069e-01
* time: 22.356987953186035
63 -1.898480e+03 3.102071e-01
* time: 23.031436920166016
64 -1.898480e+03 3.102074e-01
* time: 23.705273866653442
65 -1.898480e+03 3.102076e-01
* time: 24.39507794380188
66 -1.898480e+03 3.102079e-01
* time: 25.07362985610962
67 -1.898480e+03 3.102081e-01
* time: 25.761672019958496
68 -1.898480e+03 3.102081e-01
* time: 26.46256685256958
69 -1.898480e+03 3.102081e-01
* time: 27.166467905044556
70 -1.898480e+03 3.102082e-01
* time: 27.866047859191895
71 -1.898480e+03 3.102082e-01
* time: 28.56844997406006
72 -1.898480e+03 3.102082e-01
* time: 29.275399923324585
73 -1.898480e+03 3.102102e-01
* time: 29.94611096382141
74 -1.898480e+03 3.102102e-01
* time: 30.50294303894043
75 -1.898480e+03 3.102096e-01
* time: 31.010547876358032
76 -1.898480e+03 3.102096e-01
* time: 31.535022974014282
77 -1.898480e+03 3.125688e-01
* time: 32.054117918014526
78 -1.898480e+03 3.125640e-01
* time: 32.485291957855225
79 -1.898480e+03 3.125618e-01
* time: 32.952473878860474
80 -1.898480e+03 3.125615e-01
* time: 33.46808195114136
81 -1.898480e+03 3.125612e-01
* time: 33.99279189109802
82 -1.898480e+03 3.125610e-01
* time: 34.47139883041382
83 -1.898480e+03 3.125609e-01
* time: 35.0406060218811
84 -1.898480e+03 3.125604e-01
* time: 35.54052495956421
85 -1.898480e+03 3.125602e-01
* time: 36.038084983825684
86 -1.898480e+03 3.125602e-01
* time: 36.58007884025574
87 -1.898480e+03 3.125602e-01
* time: 37.2889039516449
88 -1.898480e+03 3.125602e-01
* time: 38.0173180103302
89 -1.898480e+03 3.125602e-01
* time: 38.768173933029175
90 -1.898480e+03 3.125602e-01
* time: 39.49185800552368
91 -1.898480e+03 3.125602e-01
* time: 40.23722696304321
92 -1.898480e+03 3.125602e-01
* time: 40.75749206542969
93 -1.898480e+03 3.125602e-01
* time: 41.25179100036621
94 -1.898480e+03 3.125602e-01
* time: 41.751808881759644
95 -1.898480e+03 1.387453e-01
* time: 42.123191833496094
96 -1.898480e+03 1.387453e-01
* time: 42.564483880996704
97 -1.898480e+03 1.387453e-01
* time: 43.05810189247131
98 -1.898480e+03 1.387453e-01
* time: 43.64898705482483
99 -1.898480e+03 1.387453e-01
* time: 44.03232002258301FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: 1898.4797
Number of subjects: 120
Number of parameters: Fixed Optimized
0 11
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
-------------------
Estimate
-------------------
tvcl 2.6191
tvv 11.378
tvvp 8.453
tvq 1.3164
tvka 4.8926
Ω₁,₁ 0.13243
Ω₂,₂ 0.059669
Ω₃,₃ 0.41581
Ω₄,₄ 0.080679
Ω₅,₅ 0.24996
σ_p 0.049097
-------------------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.6191 |
| 2 | tvv | 11.0046 | 11.3784 |
| 3 | tvvp | 5.53998 | 8.45297 |
| 4 | tvq | 1.51591 | 1.31637 |
| 5 | tvka | 2.0 | 4.89257 |
| 6 | Ω₁,₁ | 0.102669 | 0.132432 |
| 7 | Ω₂,₂ | 0.0607757 | 0.0596693 |
| 8 | Ω₃,₃ | 1.20115 | 0.415811 |
| 9 | Ω₄,₄ | 0.423495 | 0.0806789 |
| 10 | Ω₅,₅ | 0.244731 | 0.249956 |
| 11 | σ_p | 0.0484049 | 0.0490975 |
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{Optim.NewtonTrustRegion{Float64}, Optim.Options{Float64, Nothing}}(Optim.NewtonTrustRegion{Float64}(1.0, 100.0, 1.4901161193847656e-8, 0.1, 0.25, 0.75, false), Optim.Options(x_abstol = 0.0, x_reltol = 0.0, f_abstol = 0.0, f_reltol = 0.0, g_abstol = 1.0e-5, g_reltol = 1.0e-8, outer_x_abstol = 0.0, outer_x_reltol = 0.0, outer_f_abstol = 0.0, outer_f_reltol = 0.0, outer_g_abstol = 1.0e-8, outer_g_reltol = 1.0e-8, f_calls_limit = 0, g_calls_limit = 0, h_calls_limit = 0, allow_f_increases = false, allow_outer_f_increases = true, successive_f_tol = 1, iterations = 1000, outer_iterations = 1000, store_trace = false, trace_simplex = false, show_trace = false, extended_trace = false, show_warnings = true, show_every = 1, time_limit = NaN, )
)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 = (; 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 = (; 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: 120
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 = (; size = (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.