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.027987957000732422
1 2.343899e+02 1.747348e+03
* time: 0.9532840251922607
2 9.696232e+01 1.198088e+03
* time: 0.9572150707244873
3 -7.818699e+01 5.538151e+02
* time: 0.9602971076965332
4 -1.234803e+02 2.462514e+02
* time: 0.9632959365844727
5 -1.372888e+02 2.067458e+02
* time: 0.9658751487731934
6 -1.410579e+02 1.162950e+02
* time: 0.9683570861816406
7 -1.434754e+02 5.632816e+01
* time: 0.970829963684082
8 -1.453401e+02 7.859270e+01
* time: 0.9737300872802734
9 -1.498185e+02 1.455606e+02
* time: 0.976639986038208
10 -1.534371e+02 1.303682e+02
* time: 0.9794831275939941
11 -1.563557e+02 5.975474e+01
* time: 0.9823429584503174
12 -1.575052e+02 9.308611e+00
* time: 0.985008955001831
13 -1.579357e+02 1.234484e+01
* time: 0.9877109527587891
14 -1.581874e+02 7.478196e+00
* time: 0.9904570579528809
15 -1.582981e+02 2.027162e+00
* time: 0.9931659698486328
16 -1.583375e+02 5.578262e+00
* time: 0.9959421157836914
17 -1.583556e+02 4.727050e+00
* time: 0.9986741542816162
18 -1.583644e+02 2.340173e+00
* time: 1.0014240741729736
19 -1.583680e+02 7.738100e-01
* time: 1.0041749477386475
20 -1.583696e+02 3.300689e-01
* time: 1.0069961547851562
21 -1.583704e+02 3.641985e-01
* time: 1.009817123413086
22 -1.583707e+02 4.365901e-01
* time: 1.012686014175415
23 -1.583709e+02 3.887800e-01
* time: 1.0155270099639893
24 -1.583710e+02 2.766977e-01
* time: 1.0183279514312744
25 -1.583710e+02 1.758029e-01
* time: 1.0208020210266113
26 -1.583710e+02 1.133947e-01
* time: 1.0235049724578857
27 -1.583710e+02 7.922544e-02
* time: 1.025979995727539
28 -1.583710e+02 5.954998e-02
* time: 1.0288059711456299
29 -1.583710e+02 4.157079e-02
* time: 1.0317561626434326
30 -1.583710e+02 4.295447e-02
* time: 1.0346770286560059
31 -1.583710e+02 5.170753e-02
* time: 1.0375621318817139
32 -1.583710e+02 2.644383e-02
* time: 1.0413401126861572
33 -1.583710e+02 4.548993e-03
* time: 1.045241117477417
34 -1.583710e+02 2.501804e-02
* time: 1.0490081310272217
35 -1.583710e+02 3.763440e-02
* time: 1.0519189834594727
36 -1.583710e+02 3.206026e-02
* time: 1.0547480583190918
37 -1.583710e+02 1.003698e-02
* time: 1.057569980621338
38 -1.583710e+02 2.209089e-02
* time: 1.06050705909729
39 -1.583710e+02 4.954172e-03
* time: 1.0633809566497803
40 -1.583710e+02 1.609373e-02
* time: 1.0672039985656738
41 -1.583710e+02 1.579802e-02
* time: 1.0700960159301758
42 -1.583710e+02 1.014113e-03
* time: 1.0729761123657227
43 -1.583710e+02 6.050644e-03
* time: 1.076789140701294
44 -1.583710e+02 1.354412e-02
* time: 1.0796821117401123
45 -1.583710e+02 4.473248e-03
* time: 1.0825550556182861
46 -1.583710e+02 4.644735e-03
* time: 1.0854289531707764
47 -1.583710e+02 9.829910e-03
* time: 1.0883369445800781
48 -1.583710e+02 1.047561e-03
* time: 1.0912880897521973
49 -1.583710e+02 8.366895e-03
* time: 1.0942180156707764
50 -1.583710e+02 7.879055e-04
* time: 1.0971190929412842FittedPumasModel
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: 7.605552673339844e-5
1 -7.022088e+02 1.707063e+02
* time: 0.3500559329986572
2 -7.314067e+02 2.903269e+02
* time: 0.5477230548858643
3 -8.520591e+02 2.285888e+02
* time: 0.8486390113830566
4 -1.120191e+03 3.795410e+02
* time: 1.1841399669647217
5 -1.178784e+03 2.323978e+02
* time: 1.3605549335479736
6 -1.218320e+03 9.699907e+01
* time: 1.535857915878296
7 -1.223641e+03 5.862105e+01
* time: 1.675379991531372
8 -1.227620e+03 1.831403e+01
* time: 1.8478319644927979
9 -1.228381e+03 2.132323e+01
* time: 2.0683600902557373
10 -1.230098e+03 2.921228e+01
* time: 2.2204511165618896
11 -1.230854e+03 2.029662e+01
* time: 2.3679521083831787
12 -1.231116e+03 5.229097e+00
* time: 2.5296781063079834
13 -1.231179e+03 1.689232e+00
* time: 2.7500901222229004
14 -1.231187e+03 1.215379e+00
* time: 2.8995139598846436
15 -1.231188e+03 2.770380e-01
* time: 3.0236079692840576
16 -1.231188e+03 1.636653e-01
* time: 3.144789934158325
17 -1.231188e+03 2.701133e-01
* time: 3.3010470867156982
18 -1.231188e+03 3.163363e-01
* time: 3.406359910964966
19 -1.231188e+03 1.505149e-01
* time: 3.5092921257019043
20 -1.231188e+03 2.484999e-02
* time: 3.6186721324920654
21 -1.231188e+03 8.446863e-04
* time: 3.7400400638580322FittedPumasModel
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.700920104980469e-5
1 -9.197817e+02 9.927951e+02
* time: 0.30173587799072266
2 -1.372640e+03 2.054986e+02
* time: 0.6254029273986816
3 -1.446326e+03 1.543987e+02
* time: 0.9318819046020508
4 -1.545570e+03 1.855028e+02
* time: 1.251338005065918
5 -1.581449e+03 1.713157e+02
* time: 1.7304089069366455
6 -1.639433e+03 1.257382e+02
* time: 2.0164828300476074
7 -1.695964e+03 7.450539e+01
* time: 2.319591999053955
8 -1.722243e+03 5.961044e+01
* time: 2.657282829284668
9 -1.736883e+03 7.320921e+01
* time: 2.933048963546753
10 -1.753547e+03 7.501938e+01
* time: 3.249656915664673
11 -1.764053e+03 6.185661e+01
* time: 3.5680859088897705
12 -1.778991e+03 4.831033e+01
* time: 3.9332008361816406
13 -1.791492e+03 4.943278e+01
* time: 4.31971001625061
14 -1.799847e+03 2.871410e+01
* time: 4.701333999633789
15 -1.805374e+03 7.520791e+01
* time: 5.085456848144531
16 -1.816260e+03 2.990621e+01
* time: 5.483062028884888
17 -1.818252e+03 2.401915e+01
* time: 5.775821924209595
18 -1.822988e+03 2.587225e+01
* time: 6.109920978546143
19 -1.824653e+03 1.550517e+01
* time: 6.40765905380249
20 -1.826074e+03 1.788927e+01
* time: 6.738783836364746
21 -1.826821e+03 1.888389e+01
* time: 7.086375951766968
22 -1.827900e+03 1.432840e+01
* time: 7.381964921951294
23 -1.828511e+03 9.422041e+00
* time: 7.712573051452637
24 -1.828754e+03 5.363442e+00
* time: 8.071256875991821
25 -1.828862e+03 4.916159e+00
* time: 8.347540855407715
26 -1.829007e+03 4.695755e+00
* time: 8.67992091178894
27 -1.829358e+03 1.090249e+01
* time: 8.991798877716064
28 -1.829830e+03 1.451325e+01
* time: 9.341611862182617
29 -1.830201e+03 1.108715e+01
* time: 9.717772960662842
30 -1.830360e+03 2.891223e+00
* time: 10.034783840179443
31 -1.830390e+03 1.695557e+00
* time: 10.410444021224976
32 -1.830404e+03 1.601712e+00
* time: 10.806506872177124
33 -1.830432e+03 2.823385e+00
* time: 11.308536052703857
34 -1.830477e+03 4.060617e+00
* time: 11.72184705734253
35 -1.830528e+03 5.133499e+00
* time: 12.04262399673462
36 -1.830593e+03 2.830970e+00
* time: 12.412673950195312
37 -1.830616e+03 3.342835e+00
* time: 12.738509893417358
38 -1.830622e+03 3.708884e+00
* time: 13.068030834197998
39 -1.830625e+03 2.062934e+00
* time: 13.440907001495361
40 -1.830627e+03 1.278569e+00
* time: 13.71698784828186
41 -1.830628e+03 1.832895e+00
* time: 14.023622035980225
42 -1.830628e+03 3.768840e-01
* time: 14.270959854125977
43 -1.830629e+03 3.152895e-01
* time: 14.533400058746338
44 -1.830630e+03 4.871060e-01
* time: 14.800099849700928
45 -1.830630e+03 3.110627e-01
* time: 15.075872898101807
46 -1.830630e+03 2.687758e-02
* time: 15.327050924301147
47 -1.830630e+03 4.694018e-03
* time: 15.556539058685303
48 -1.830630e+03 8.272969e-03
* time: 15.758961915969849
49 -1.830630e+03 8.249151e-03
* time: 16.019502878189087
50 -1.830630e+03 8.245562e-03
* time: 16.368154048919678
51 -1.830630e+03 8.240030e-03
* time: 16.672675848007202
52 -1.830630e+03 8.240030e-03
* time: 17.047104835510254
53 -1.830630e+03 8.240030e-03
* time: 17.422839879989624FittedPumasModel
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 | 17.423 | 3.74 |
| 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: 7.796287536621094e-5
1 -8.682982e+02 1.000199e+03
* time: 0.36192989349365234
2 -1.381870e+03 5.008081e+02
* time: 0.7124800682067871
3 -1.551053e+03 6.833490e+02
* time: 1.1106960773468018
4 -1.680887e+03 1.834586e+02
* time: 1.4103178977966309
5 -1.726118e+03 8.870274e+01
* time: 1.764415979385376
6 -1.761023e+03 1.162036e+02
* time: 2.0781610012054443
7 -1.786619e+03 1.114552e+02
* time: 2.434468984603882
8 -1.863556e+03 9.914305e+01
* time: 2.77034592628479
9 -1.882942e+03 5.342676e+01
* time: 3.134916067123413
10 -1.888020e+03 2.010181e+01
* time: 3.4656410217285156
11 -1.889832e+03 1.867263e+01
* time: 3.836688995361328
12 -1.891649e+03 1.668512e+01
* time: 4.203073024749756
13 -1.892615e+03 1.820701e+01
* time: 4.510792016983032
14 -1.893453e+03 1.745195e+01
* time: 4.867301940917969
15 -1.894760e+03 1.850174e+01
* time: 5.179548978805542
16 -1.895647e+03 1.773939e+01
* time: 5.524960041046143
17 -1.896597e+03 1.143462e+01
* time: 5.84778904914856
18 -1.897114e+03 9.720097e+00
* time: 6.198713064193726
19 -1.897373e+03 6.054321e+00
* time: 6.517086982727051
20 -1.897498e+03 3.985954e+00
* time: 6.8590168952941895
21 -1.897571e+03 4.262464e+00
* time: 7.215123891830444
22 -1.897633e+03 4.010234e+00
* time: 7.509742021560669
23 -1.897714e+03 4.805375e+00
* time: 7.863183975219727
24 -1.897802e+03 3.508706e+00
* time: 8.164668083190918
25 -1.897865e+03 3.691477e+00
* time: 8.510861873626709
26 -1.897900e+03 2.982720e+00
* time: 8.806504011154175
27 -1.897928e+03 2.563790e+00
* time: 9.13526201248169
28 -1.897968e+03 3.261485e+00
* time: 9.438282012939453
29 -1.898013e+03 3.064690e+00
* time: 9.77065896987915
30 -1.898040e+03 1.636525e+00
* time: 10.076997995376587
31 -1.898051e+03 1.439997e+00
* time: 10.409209966659546
32 -1.898057e+03 1.436504e+00
* time: 10.710906028747559
33 -1.898069e+03 1.881529e+00
* time: 11.043354034423828
34 -1.898095e+03 3.253165e+00
* time: 11.396759033203125
35 -1.898142e+03 4.257942e+00
* time: 11.693532943725586
36 -1.898199e+03 3.685241e+00
* time: 12.053384065628052
37 -1.898245e+03 2.567364e+00
* time: 12.37402892112732
38 -1.898246e+03 2.561591e+00
* time: 12.92623496055603
39 -1.898251e+03 2.530888e+00
* time: 13.408470869064331
40 -1.898298e+03 2.673696e+00
* time: 13.712632894515991
41 -1.898300e+03 2.794639e+00
* time: 14.137242078781128
42 -1.898337e+03 3.751590e+00
* time: 14.577944040298462
43 -1.898421e+03 4.878407e+00
* time: 14.942821025848389
44 -1.898433e+03 4.391719e+00
* time: 15.393600940704346
45 -1.898437e+03 4.216518e+00
* time: 15.871919870376587
46 -1.898442e+03 4.108397e+00
* time: 16.4161958694458
47 -1.898446e+03 3.934902e+00
* time: 16.939286947250366
48 -1.898449e+03 3.769838e+00
* time: 17.46850609779358
49 -1.898450e+03 3.739486e+00
* time: 17.91516089439392
50 -1.898450e+03 3.712049e+00
* time: 18.455398082733154
51 -1.898457e+03 3.623436e+00
* time: 18.93416690826416
52 -1.898471e+03 2.668312e+00
* time: 19.26013994216919
53 -1.898479e+03 2.302438e+00
* time: 19.616926908493042
54 -1.898480e+03 2.386566e-01
* time: 19.953322887420654
55 -1.898480e+03 7.802040e-01
* time: 20.30431294441223
56 -1.898480e+03 7.369786e-01
* time: 20.823550939559937
57 -1.898480e+03 5.113191e-01
* time: 21.213505029678345
58 -1.898480e+03 3.067709e-01
* time: 21.53412699699402
59 -1.898480e+03 3.076791e-01
* time: 21.819719076156616
60 -1.898480e+03 3.102066e-01
* time: 22.13746190071106
61 -1.898480e+03 3.102066e-01
* time: 22.524199962615967
62 -1.898480e+03 3.102069e-01
* time: 22.997545957565308
63 -1.898480e+03 3.102071e-01
* time: 23.736550092697144
64 -1.898480e+03 3.102074e-01
* time: 24.44553303718567
65 -1.898480e+03 3.102076e-01
* time: 25.174213886260986
66 -1.898480e+03 3.102079e-01
* time: 25.88929009437561
67 -1.898480e+03 3.102081e-01
* time: 26.58518695831299
68 -1.898480e+03 3.102081e-01
* time: 27.301661014556885
69 -1.898480e+03 3.102081e-01
* time: 28.0117609500885
70 -1.898480e+03 3.102082e-01
* time: 28.73285698890686
71 -1.898480e+03 3.102082e-01
* time: 29.470452070236206
72 -1.898480e+03 3.102082e-01
* time: 30.202401876449585
73 -1.898480e+03 3.102102e-01
* time: 30.903527975082397
74 -1.898480e+03 3.102102e-01
* time: 31.427478075027466
75 -1.898480e+03 3.102096e-01
* time: 31.978260040283203
76 -1.898480e+03 3.102096e-01
* time: 32.512479066848755
77 -1.898480e+03 3.125688e-01
* time: 33.00318193435669
78 -1.898480e+03 3.125640e-01
* time: 33.45245289802551
79 -1.898480e+03 3.125618e-01
* time: 33.92856287956238
80 -1.898480e+03 3.125615e-01
* time: 34.46659588813782
81 -1.898480e+03 3.125612e-01
* time: 35.013099908828735
82 -1.898480e+03 3.125610e-01
* time: 35.50523591041565
83 -1.898480e+03 3.125609e-01
* time: 36.11022090911865
84 -1.898480e+03 3.125604e-01
* time: 36.63760805130005
85 -1.898480e+03 3.125602e-01
* time: 37.15481090545654
86 -1.898480e+03 3.125602e-01
* time: 37.664392948150635
87 -1.898480e+03 3.125602e-01
* time: 38.403444051742554
88 -1.898480e+03 3.125602e-01
* time: 39.20531988143921
89 -1.898480e+03 3.125602e-01
* time: 39.98598098754883
90 -1.898480e+03 3.125602e-01
* time: 40.75041389465332
91 -1.898480e+03 3.125602e-01
* time: 41.57505202293396
92 -1.898480e+03 3.125602e-01
* time: 42.0564661026001
93 -1.898480e+03 3.125602e-01
* time: 42.61239194869995
94 -1.898480e+03 3.125602e-01
* time: 43.13966488838196
95 -1.898480e+03 1.387453e-01
* time: 43.48903489112854
96 -1.898480e+03 1.387453e-01
* time: 43.91609001159668
97 -1.898480e+03 1.387453e-01
* time: 44.453840017318726
98 -1.898480e+03 1.387453e-01
* time: 45.04751992225647
99 -1.898480e+03 1.387453e-01
* time: 45.414485931396484FittedPumasModel
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.