using Pumas
using PumasUtilities
using NCA
using NCAUtilities
A Comprehensive Introduction to Pumas
This tutorial provides a comprehensive introduction to a modeling and simulation workflow in Pumas. The idea is not to get into the details of Pumas specifics, but instead provide a narrative on the lines of a regular workflow in our day-to-day work, with brevity where required to allow a broad overview. Wherever possible, cross-references will be provided to documentation and detailed examples that provide deeper insight into a particular topic.
As part of this workflow, you will be introduced to various aspects such as:
- Data wrangling in Julia
- Exploratory analysis in Julia
- Continuous data non-linear mixed effects modeling in Pumas
- Model comparison routines, post-processing, validation etc.
1 The Study and Design
CTMNopain is a novel anti-inflammatory agent under preliminary investigation. A dose-ranging trial was conducted comparing placebo with 3 doses of CTMNopain (5mg, 20mg and 80 mg QD). The maximum tolerated dose is 160 mg per day. Plasma concentrations (mg/L) of the drug were measured at 0, 0.5, 1, 1.5, 2, 2.5, 3-8 hours.
Pain score (0=no pain, 1=mild, 2=moderate, 3=severe) were obtained at time points when plasma concentration was collected. A pain score of 2 or more is considered as no pain relief.
The subjects can request for remedication if pain relief is not achieved after 2 hours post dose. Some subjects had remedication before 2 hours if they were not able to bear the pain. The time to remedication and the remedication status is available for subjects.
The pharmacokinetic dataset can be accessed using PharmaDatasets.jl.
2 Setup
2.1 Load libraries
These libraries provide the workhorse functionality in the Pumas ecosystem:
In addition, libraries below are good add-on’s that provide ancillary functionality:
using GLM: lm, @formula
using Random
using CSV
using DataFramesMeta
using CairoMakie
using PharmaDatasets2.2 Data Wrangling
We start by reading in the dataset and making some quick summaries.
Tip
If you want to learn more about data wrangling, don’t forget to check our Data Wrangling in Julia tutorials!
pkpain_df = dataset("pk_painrelief")
first(pkpain_df, 5)5×7 DataFrame
| Row | Subject | Time | Conc | PainRelief | PainScore | RemedStatus | Dose |
|---|---|---|---|---|---|---|---|
| Int64 | Float64 | Float64 | Int64 | Int64 | Int64 | String7 | |
| 1 | 1 | 0.0 | 0.0 | 0 | 3 | 1 | 20 mg |
| 2 | 1 | 0.5 | 1.15578 | 1 | 1 | 0 | 20 mg |
| 3 | 1 | 1.0 | 1.37211 | 1 | 0 | 0 | 20 mg |
| 4 | 1 | 1.5 | 1.30058 | 1 | 0 | 0 | 20 mg |
| 5 | 1 | 2.0 | 1.19195 | 1 | 1 | 0 | 20 mg |
Let’s filter out the placebo data as we don’t need that for the PK analysis.
pkpain_noplb_df = @rsubset pkpain_df :Dose != "Placebo";
first(pkpain_noplb_df, 5)5×7 DataFrame
| Row | Subject | Time | Conc | PainRelief | PainScore | RemedStatus | Dose |
|---|---|---|---|---|---|---|---|
| Int64 | Float64 | Float64 | Int64 | Int64 | Int64 | String7 | |
| 1 | 1 | 0.0 | 0.0 | 0 | 3 | 1 | 20 mg |
| 2 | 1 | 0.5 | 1.15578 | 1 | 1 | 0 | 20 mg |
| 3 | 1 | 1.0 | 1.37211 | 1 | 0 | 0 | 20 mg |
| 4 | 1 | 1.5 | 1.30058 | 1 | 0 | 0 | 20 mg |
| 5 | 1 | 2.0 | 1.19195 | 1 | 1 | 0 | 20 mg |
3 Analysis
3.1 Non-compartmental analysis
Let’s begin by performing a quick NCA of the concentration time profiles and view the exposure changes across doses. The input data specification for NCA analysis requires the presence of a :route column and an :amt column that specifies the dose. So, let’s add that in:
@rtransform! pkpain_noplb_df begin
:route = "ev"
:Dose = parse(Int, chop(:Dose; tail = 3))
endWe also need to create an :amt column:
@rtransform! pkpain_noplb_df :amt = :Time == 0 ? :Dose : missingNow, we map the data variables to the read_nca function that prepares the data for NCA analysis.
pkpain_nca = read_nca(
pkpain_noplb_df;
id = :Subject,
time = :Time,
amt = :amt,
observations = :Conc,
group = [:Dose],
route = :route,
)NCAPopulation (120 subjects):
Group: [["Dose" => 5], ["Dose" => 20], ["Dose" => 80]]
Number of missing observations: 0
Number of blq observations: 0Now that we mapped the data in, let’s visualize the concentration vs time plots for a few individuals. When paginate is set to true, a vector of plots are returned and below we display the first element with 9 individuals.
f = observations_vs_time(
pkpain_nca;
paginate = true,
axis = (; xlabel = "Time (hr)", ylabel = "CTMNoPain Concentration (ng/mL)"),
facet = (; combinelabels = true),
)
f[1]
or you can view the summary curves by dose group as passed in to the group argument in read_nca
summary_observations_vs_time(
pkpain_nca,
figure = (; fontsize = 22, resolution = (800, 1000)),
color = "black",
linewidth = 3,
axis = (; xlabel = "Time (hr)", ylabel = "CTMX Concentration (μg/mL)"),
facet = (; combinelabels = true, linkaxes = true),
)
A full NCA Report is now obtained for completeness purposes using the run_nca function, but later we will only extract a couple of key metrics of interest.
pk_nca = run_nca(pkpain_nca; sigdigits = 3)We can look at the NCA fits for some subjects. Here f is a vector or figures. We’ll showcase the first image by indexing f:
f = subject_fits(
pk_nca,
paginate = true,
axis = (; xlabel = "Time (hr)", ylabel = "CTMX Concentration (μg/mL)"),
facet = (; combinelabels = true, linkaxes = true),
)
f[1]
As CTMNopain’s effect maybe mainly related to maximum concentration (cmax) or area under the curve (auc), we present some summary statistics using the summarize function from NCA.
strata = [:Dose]1-element Vector{Symbol}:
:Doseparams = [:cmax, :aucinf_obs]2-element Vector{Symbol}:
:cmax
:aucinf_obsoutput = summarize(pk_nca; stratify_by = strata, parameters = params)6×10 DataFrame
| Row | Dose | parameters | numsamples | minimum | maximum | mean | std | geomean | geostd | geomeanCV |
|---|---|---|---|---|---|---|---|---|---|---|
| Int64 | String | Int64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | 5 | cmax | 40 | 0.19 | 0.539 | 0.356075 | 0.0884129 | 0.345104 | 1.2932 | 26.1425 |
| 2 | 5 | aucinf_obs | 40 | 0.914 | 3.4 | 1.5979 | 0.490197 | 1.53373 | 1.32974 | 29.0868 |
| 3 | 20 | cmax | 40 | 0.933 | 2.7 | 1.4737 | 0.361871 | 1.43408 | 1.2633 | 23.6954 |
| 4 | 20 | aucinf_obs | 40 | 2.77 | 14.1 | 6.377 | 2.22239 | 6.02031 | 1.41363 | 35.6797 |
| 5 | 80 | cmax | 40 | 3.3 | 8.47 | 5.787 | 1.31957 | 5.64164 | 1.25757 | 23.2228 |
| 6 | 80 | aucinf_obs | 40 | 13.7 | 49.1 | 29.5 | 8.68984 | 28.2954 | 1.34152 | 30.0258 |
The statistics printed above are the default, but you can pass in your own statistics using the stats = [] argument to the summarize function.
We can look at a few parameter distribution plots.
parameters_vs_group(
pk_nca,
parameter = :cmax,
axis = (; xlabel = "Dose (mg)", ylabel = "Cₘₐₓ (ng/mL)"),
figure = (; fontsize = 18),
)
Dose normalized PK parameters, cmax and aucinf were essentially dose proportional between for 5 mg, 20 mg and 80 mg doses. You can perform a simple regression to check the impact of dose on cmax:
dp = NCA.DoseLinearityPowerModel(pk_nca, :cmax; level = 0.9)Dose Linearity Power Model
Variable: cmax
Model: log(cmax) ~ log(α) + β × log(dose)
────────────────────────────────────
Estimate low CI 90% high CI 90%
────────────────────────────────────
β 1.00775 0.97571 1.0398
────────────────────────────────────Here’s a visualization for the dose linearity using a power model for cmax:
power_model(dp)
We can also visualize a dose proportionality results with respect to a specific endpoint in a NCA Report; for example cmax and aucinf_obs:
dose_vs_dose_normalized(pk_nca, :cmax)
dose_vs_dose_normalized(pk_nca, :aucinf_obs)
Based on visual inspection of the concentration time profiles as seen earlier, CTMNopain exhibited monophasic decline, and perhaps a one compartment model best fits the PK data.
3.2 Pharmacokinetic modeling
As seen from the plots above, the concentrations decline monoexponentially. We will evaluate both one and two compartment structural models to assess best fit. Further, different residual error models will also be tested.
We will use the results from NCA to provide us good initial estimates.
3.2.1 Data preparation for modeling
PumasNDF requires the presence of :evid and :cmt columns in the dataset.
@rtransform! pkpain_noplb_df begin
:evid = :Time == 0 ? 1 : 0
:cmt = :Time == 0 ? 1 : 2
:cmt2 = 1 # for zero order absorption
endFurther, observations at time of dosing, i.e., when evid = 1 have to be missing
@rtransform! pkpain_noplb_df :Conc = :evid == 1 ? missing : :ConcThe dataframe will now be converted to a Population using read_pumas. Note that both observations and covariates are required to be an array even if it is one element.
pkpain_noplb = read_pumas(
pkpain_noplb_df;
id = :Subject,
time = :Time,
amt = :amt,
observations = [:Conc],
covariates = [:Dose],
evid = :evid,
cmt = :cmt,
)Population
Subjects: 120
Covariates: Dose
Observations: ConcNow that the data is transformed to a Population of subjects, we can explore different models.
3.2.2 One-compartment model
Note
If you are not familiar yet with the @model blocks and syntax, please check our documentation.
pk_1cmp = @model begin
@metadata begin
desc = "One Compartment Model"
timeu = u"hr"
end
@param begin
"""
Clearance (L/hr)
"""
tvcl ∈ RealDomain(; lower = 0, init = 3.2)
"""
Volume (L)
"""
tvv ∈ RealDomain(; lower = 0, init = 16.4)
"""
Absorption rate constant (h-1)
"""
tvka ∈ RealDomain(; lower = 0, init = 3.8)
"""
- ΩCL
- ΩVc
- ΩKa
"""
Ω ∈ PDiagDomain(init = [0.04, 0.04, 0.04])
"""
Proportional RUV
"""
σ_p ∈ RealDomain(; lower = 0.0001, init = 0.2)
end
@random begin
η ~ MvNormal(Ω)
end
@covariates begin
"""
Dose (mg)
"""
Dose
end
@pre begin
CL = tvcl * exp(η[1])
Vc = tvv * exp(η[2])
Ka = tvka * exp(η[3])
end
@dynamics Depots1Central1
@derived begin
cp := @. Central / Vc
"""
CTMx Concentration (ng/mL)
"""
Conc ~ @. Normal(cp, abs(cp) * σ_p)
end
endPumasModel
Parameters: tvcl, tvv, tvka, Ω, σ_p
Random effects: η
Covariates: Dose
Dynamical 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.026766061782836914
1 2.343899e+02 1.747348e+03
* time: 0.4976530075073242
2 9.696232e+01 1.198088e+03
* time: 0.4996790885925293
3 -7.818699e+01 5.538151e+02
* time: 0.5012938976287842
4 -1.234803e+02 2.462514e+02
* time: 0.5029759407043457
5 -1.372888e+02 2.067458e+02
* time: 0.5046780109405518
6 -1.410579e+02 1.162950e+02
* time: 0.5551629066467285
7 -1.434754e+02 5.632816e+01
* time: 0.5564119815826416
8 -1.453401e+02 7.859270e+01
* time: 0.5575258731842041
9 -1.498185e+02 1.455606e+02
* time: 0.5586040019989014
10 -1.534371e+02 1.303682e+02
* time: 0.5597009658813477
11 -1.563557e+02 5.975474e+01
* time: 0.5607888698577881
12 -1.575052e+02 9.308611e+00
* time: 0.5620739459991455
13 -1.579357e+02 1.234484e+01
* time: 0.5636448860168457
14 -1.581874e+02 7.478196e+00
* time: 0.565234899520874
15 -1.582981e+02 2.027162e+00
* time: 0.5668668746948242
16 -1.583375e+02 5.578262e+00
* time: 0.5684390068054199
17 -1.583556e+02 4.727050e+00
* time: 0.5700418949127197
18 -1.583644e+02 2.340173e+00
* time: 0.5716519355773926
19 -1.583680e+02 7.738100e-01
* time: 0.5735750198364258
20 -1.583696e+02 3.300689e-01
* time: 0.5752630233764648
21 -1.583704e+02 3.641985e-01
* time: 0.5770928859710693
22 -1.583707e+02 4.365901e-01
* time: 0.5788829326629639
23 -1.583709e+02 3.887800e-01
* time: 0.580543041229248
24 -1.583710e+02 2.766977e-01
* time: 0.5821409225463867
25 -1.583710e+02 1.758029e-01
* time: 0.5838029384613037
26 -1.583710e+02 1.133947e-01
* time: 0.5854370594024658
27 -1.583710e+02 7.922544e-02
* time: 0.5869929790496826
28 -1.583710e+02 5.954998e-02
* time: 0.5885839462280273
29 -1.583710e+02 4.157079e-02
* time: 0.5901889801025391
30 -1.583710e+02 4.295447e-02
* time: 0.5918159484863281
31 -1.583710e+02 5.170753e-02
* time: 0.5934419631958008
32 -1.583710e+02 2.644383e-02
* time: 0.5954558849334717
33 -1.583710e+02 4.548993e-03
* time: 0.5973970890045166
34 -1.583710e+02 2.501804e-02
* time: 0.5995140075683594
35 -1.583710e+02 3.763440e-02
* time: 0.601172924041748
36 -1.583710e+02 3.206026e-02
* time: 0.6028480529785156
37 -1.583710e+02 1.003698e-02
* time: 0.604496955871582
38 -1.583710e+02 2.209089e-02
* time: 0.6061689853668213
39 -1.583710e+02 4.954172e-03
* time: 0.6078689098358154
40 -1.583710e+02 1.609373e-02
* time: 0.6099510192871094
41 -1.583710e+02 1.579802e-02
* time: 0.6115579605102539
42 -1.583710e+02 1.014113e-03
* time: 0.6131930351257324
43 -1.583710e+02 6.050644e-03
* time: 0.6152610778808594
44 -1.583710e+02 1.354412e-02
* time: 0.616919994354248
45 -1.583710e+02 4.473248e-03
* time: 0.6558220386505127
46 -1.583710e+02 4.644735e-03
* time: 0.657249927520752
47 -1.583710e+02 9.829910e-03
* time: 0.6583340167999268
48 -1.583710e+02 1.047561e-03
* time: 0.6594228744506836
49 -1.583710e+02 8.366895e-03
* time: 0.6604878902435303
50 -1.583710e+02 7.879055e-04
* time: 0.6616060733795166FittedPumasModel
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), x_gap = 0)
The distribution of the loglikelihood’s suggest no extreme outliers.
A more convenient way is to use the findinfluential function that provides a list of k top influential subjects by showing the normalized (minus) loglikelihood for each subject. As you can see below, the minus loglikelihood in the range of 16 agrees with the histogram plotted above.
influential_subjects = findinfluential(pk_1cmp, pkpain_noplb, pkparam, FOCE())120-element Vector{NamedTuple{(:id, :nll), Tuple{String, Float64}}}:
(id = "148", nll = 16.65965885684477)
(id = "135", nll = 16.64898519007633)
(id = "156", nll = 15.959069556607496)
(id = "159", nll = 15.441218240496486)
(id = "149", nll = 14.715134644119509)
(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: 8.296966552734375e-5
1 -7.022088e+02 1.707063e+02
* time: 0.20783185958862305
2 -7.314067e+02 2.903269e+02
* time: 0.30812597274780273
3 -8.520591e+02 2.285888e+02
* time: 0.4109029769897461
4 -1.120191e+03 3.795410e+02
* time: 0.6144449710845947
5 -1.178784e+03 2.323978e+02
* time: 0.708690881729126
6 -1.218320e+03 9.699907e+01
* time: 0.7945859432220459
7 -1.223641e+03 5.862105e+01
* time: 0.8768119812011719
8 -1.227620e+03 1.831403e+01
* time: 0.9596998691558838
9 -1.228381e+03 2.132323e+01
* time: 1.0385079383850098
10 -1.230098e+03 2.921228e+01
* time: 1.1187078952789307
11 -1.230854e+03 2.029662e+01
* time: 1.2044639587402344
12 -1.231116e+03 5.229097e+00
* time: 1.2802090644836426
13 -1.231179e+03 1.689232e+00
* time: 1.3406670093536377
14 -1.231187e+03 1.215379e+00
* time: 1.4132559299468994
15 -1.231188e+03 2.770380e-01
* time: 1.4812140464782715
16 -1.231188e+03 1.636653e-01
* time: 1.5407640933990479
17 -1.231188e+03 2.701133e-01
* time: 1.6017448902130127
18 -1.231188e+03 3.163363e-01
* time: 1.649095058441162
19 -1.231188e+03 1.505149e-01
* time: 1.7093169689178467
20 -1.231188e+03 2.484999e-02
* time: 1.765408992767334
21 -1.231188e+03 8.446863e-04
* time: 1.8042089939117432FittedPumasModel
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.086619 [ 2.9944 ; 3.334 ]
tvv 13.288 0.27481 [12.749 ; 13.827 ]
tvka 2.0 NaN [ NaN ; NaN ]
Ω₁,₁ 0.08494 0.011022 [ 0.063338; 0.10654 ]
Ω₂,₂ 0.048568 0.0063501 [ 0.036122; 0.061014]
Ω₃,₃ 5.5811 1.2194 [ 3.1911 ; 7.9711 ]
σ_p 0.10093 0.0057196 [ 0.089718; 0.11214 ]
-------------------------------------------------------------------Notice that tvka is fixed to 2 as we don’t have a lot of information before tmax. From the results above, we see that the parameter precision for this model is reasonable.
3.2.3 Two-compartment model
Just to be sure, let’s fit a 2-compartment model and evaluate:
pk_2cmp = @model begin
@param begin
"""
Clearance (L/hr)
"""
tvcl ∈ RealDomain(; lower = 0, init = 3.2)
"""
Central Volume (L)
"""
tvv ∈ RealDomain(; lower = 0, init = 16.4)
"""
Peripheral Volume (L)
"""
tvvp ∈ RealDomain(; lower = 0, init = 10)
"""
Distributional Clearance (L/hr)
"""
tvq ∈ RealDomain(; lower = 0, init = 2)
"""
Absorption rate constant (h-1)
"""
tvka ∈ RealDomain(; lower = 0, init = 1.3)
"""
- ΩCL
- ΩVc
- ΩKa
- ΩVp
- ΩQ
"""
Ω ∈ PDiagDomain(init = [0.04, 0.04, 0.04, 0.04, 0.04])
"""
Proportional RUV
"""
σ_p ∈ RealDomain(; lower = 0.0001, init = 0.2)
end
@random begin
η ~ MvNormal(Ω)
end
@covariates begin
"""
Dose (mg)
"""
Dose
end
@pre begin
CL = tvcl * exp(η[1])
Vc = tvv * exp(η[2])
Ka = tvka * exp(η[3])
Vp = tvvp * exp(η[4])
Q = tvq * exp(η[5])
end
@dynamics Depots1Central1Periph1
@derived begin
cp := @. Central / Vc
"""
CTMx Concentration (ng/mL)
"""
Conc ~ @. Normal(cp, cp * σ_p)
end
endPumasModel
Parameters: tvcl, tvv, tvvp, tvq, tvka, Ω, σ_p
Random effects: η
Covariates: Dose
Dynamical 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.009506225585938e-5
1 -9.197817e+02 9.927951e+02
* time: 0.2316291332244873
2 -1.372640e+03 2.054986e+02
* time: 0.4178810119628906
3 -1.446326e+03 1.543987e+02
* time: 0.6274740695953369
4 -1.545570e+03 1.855028e+02
* time: 0.8226001262664795
5 -1.581449e+03 1.713157e+02
* time: 1.1153671741485596
6 -1.639433e+03 1.257382e+02
* time: 1.314897060394287
7 -1.695964e+03 7.450539e+01
* time: 1.5016250610351562
8 -1.722243e+03 5.961044e+01
* time: 1.7025911808013916
9 -1.736883e+03 7.320921e+01
* time: 1.8853960037231445
10 -1.753547e+03 7.501938e+01
* time: 2.0762460231781006
11 -1.764053e+03 6.185661e+01
* time: 2.2891762256622314
12 -1.778991e+03 4.831033e+01
* time: 2.4988701343536377
13 -1.791492e+03 4.943278e+01
* time: 2.7384450435638428
14 -1.799847e+03 2.871410e+01
* time: 2.9735631942749023
15 -1.805374e+03 7.520789e+01
* time: 3.207710027694702
16 -1.816260e+03 2.990621e+01
* time: 3.419020175933838
17 -1.818252e+03 2.401915e+01
* time: 3.6309170722961426
18 -1.822988e+03 2.587225e+01
* time: 3.8329060077667236
19 -1.824653e+03 1.550517e+01
* time: 4.039568185806274
20 -1.826074e+03 1.788927e+01
* time: 4.233752012252808
21 -1.826821e+03 1.888389e+01
* time: 4.445890188217163
22 -1.827900e+03 1.432840e+01
* time: 4.645512104034424
23 -1.828511e+03 9.422040e+00
* time: 4.860310077667236
24 -1.828754e+03 5.363444e+00
* time: 5.071684122085571
25 -1.828862e+03 4.916167e+00
* time: 5.286909103393555
26 -1.829007e+03 4.695750e+00
* time: 5.4859092235565186
27 -1.829358e+03 1.090245e+01
* time: 5.692968130111694
28 -1.829830e+03 1.451320e+01
* time: 5.9144251346588135
29 -1.830201e+03 1.108694e+01
* time: 6.140575170516968
30 -1.830360e+03 2.892320e+00
* time: 6.346736192703247
31 -1.830390e+03 1.699267e+00
* time: 6.559553146362305
32 -1.830404e+03 1.602159e+00
* time: 6.749340057373047
33 -1.830432e+03 2.822847e+00
* time: 6.945397138595581
34 -1.830475e+03 4.105639e+00
* time: 7.164409160614014
35 -1.830527e+03 5.093685e+00
* time: 7.366429090499878
36 -1.830592e+03 2.697353e+00
* time: 7.584075212478638
37 -1.830615e+03 3.468352e+00
* time: 7.785872220993042
38 -1.830623e+03 2.594581e+00
* time: 7.998565196990967
39 -1.830625e+03 1.770110e+00
* time: 8.206013202667236
40 -1.830627e+03 1.040041e+00
* time: 8.385918140411377
41 -1.830628e+03 1.124112e+00
* time: 8.588779211044312
42 -1.830628e+03 3.547258e-01
* time: 8.770697116851807
43 -1.830629e+03 4.070652e-01
* time: 8.946278095245361
44 -1.830630e+03 7.177755e-01
* time: 9.1458420753479
45 -1.830630e+03 5.121219e-01
* time: 9.324467182159424
46 -1.830630e+03 9.584921e-02
* time: 9.500436067581177
47 -1.830630e+03 1.039362e-02
* time: 9.659080028533936
48 -1.830630e+03 5.788836e-03
* time: 9.809008121490479
49 -1.830630e+03 6.037375e-03
* time: 9.957494020462036
50 -1.830630e+03 7.187933e-03
* time: 10.107605218887329
51 -1.830630e+03 7.187933e-03
* time: 10.320370197296143
52 -1.830630e+03 7.187933e-03
* time: 10.516804218292236
53 -1.830630e+03 7.212102e-03
* time: 10.709197044372559
54 -1.830630e+03 7.214659e-03
* time: 10.886619091033936
55 -1.830630e+03 7.217015e-03
* time: 11.067561149597168
56 -1.830630e+03 7.219581e-03
* time: 11.264463186264038
57 -1.830630e+03 7.222151e-03
* time: 11.442285060882568
58 -1.830630e+03 7.224724e-03
* time: 11.623488187789917
59 -1.830630e+03 7.227302e-03
* time: 11.821527004241943
60 -1.830630e+03 7.227561e-03
* time: 12.005751132965088
61 -1.830630e+03 7.227819e-03
* time: 12.193777084350586
62 -1.830630e+03 7.228077e-03
* time: 12.39099907875061
63 -1.830630e+03 7.228103e-03
* time: 12.579838037490845
64 -1.830630e+03 7.228129e-03
* time: 12.783632040023804
65 -1.830630e+03 7.228155e-03
* time: 12.97727108001709
66 -1.830630e+03 7.228180e-03
* time: 13.183101177215576
67 -1.830630e+03 7.228206e-03
* time: 13.372436046600342
68 -1.830630e+03 7.228209e-03
* time: 13.588906049728394
69 -1.830630e+03 7.228211e-03
* time: 13.79707407951355
70 -1.830630e+03 7.228214e-03
* time: 13.98587417602539
71 -1.830630e+03 7.228217e-03
* time: 14.19002914428711
72 -1.830630e+03 7.228217e-03
* time: 14.385045051574707
73 -1.830630e+03 7.228217e-03
* time: 14.675770998001099
74 -1.830630e+03 7.228217e-03
* time: 15.021639108657837
75 -1.830630e+03 7.228218e-03
* time: 15.29368805885315
76 -1.830630e+03 7.228218e-03
* time: 15.553704023361206
77 -1.830630e+03 7.228218e-03
* time: 15.809790134429932
78 -1.830630e+03 7.228218e-03
* time: 16.060631036758423
79 -1.830630e+03 7.228218e-03
* time: 16.324639081954956
80 -1.830630e+03 7.214381e-03
* time: 16.512296199798584
81 -1.830630e+03 7.214381e-03
* time: 16.77218508720398
82 -1.830630e+03 7.214381e-03
* time: 17.05467200279236
83 -1.830630e+03 7.214381e-03
* time: 17.321096181869507
84 -1.830630e+03 7.214381e-03
* time: 17.585496187210083FittedPumasModel
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.42349
Ω₅,₅ 0.24473
σ_p 0.048405
-------------------3.3 Comparing One- versus Two-compartment models
The 2-compartment model has a much lower objective function compared to the 1-compartment. Let’s compare the estimates from the 2 models using the compare_estimates function.
compare_estimates(; pkfit_1cmp, pkfit_2cmp)11×3 DataFrame
| Row | parameter | pkfit_1cmp | pkfit_2cmp |
|---|---|---|---|
| String | Float64? | Float64? | |
| 1 | tvcl | 3.1642 | 2.81378 |
| 2 | tvv | 13.288 | 11.0046 |
| 3 | tvka | 2.0 | 2.0 |
| 4 | Ω₁,₁ | 0.0849405 | 0.102669 |
| 5 | Ω₂,₂ | 0.0485682 | 0.0607755 |
| 6 | Ω₃,₃ | 5.58107 | 1.20116 |
| 7 | σ_p | 0.100928 | 0.0484049 |
| 8 | tvvp | missing | 5.53998 |
| 9 | tvq | missing | 1.51591 |
| 10 | Ω₄,₄ | missing | 0.423494 |
| 11 | Ω₅,₅ | missing | 0.244732 |
We perform a likelihood ratio test to compare the two nested models. The test statistic and the \(p\)-value clearly indicate that a 2-compartment model should be preferred.
lrtest(pkfit_1cmp, pkfit_2cmp)Statistic: 1200.0
Degrees of freedom: 4
P-value: 0.0We should also compare the other metrics and statistics, such ηshrinkage, ϵshrinkage, aic, and bic using the metrics_table function.
@chain metrics_table(pkfit_2cmp) begin
leftjoin(metrics_table(pkfit_1cmp); on = :Metric, makeunique = true)
rename!(:Value => :pk2cmp, :Value_1 => :pk1cmp)
end20×3 DataFrame
| Row | Metric | pk2cmp | pk1cmp |
|---|---|---|---|
| String | Any | Any | |
| 1 | Successful | true | true |
| 2 | Estimation Time | 17.586 | 1.804 |
| 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_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_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 = (; resolution = (600, 800)),
)
Clearly, our guess at tvka seems off-target. Let’s try and estimate tvka instead of fixing it to 2:
pkfit_2cmp_unfix_ka = fit(pk_2cmp, pkpain_noplb, init_params(pk_2cmp), FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 -3.200734e+02 1.272671e+03
* time: 9.179115295410156e-5
1 -8.682982e+02 1.000199e+03
* time: 0.2515130043029785
2 -1.381870e+03 5.008081e+02
* time: 0.43483400344848633
3 -1.551053e+03 6.833490e+02
* time: 0.6496539115905762
4 -1.680887e+03 1.834586e+02
* time: 0.8517158031463623
5 -1.726118e+03 8.870274e+01
* time: 1.0524239540100098
6 -1.761023e+03 1.162036e+02
* time: 1.256809949874878
7 -1.786619e+03 1.114552e+02
* time: 1.4591748714447021
8 -1.863556e+03 9.914305e+01
* time: 1.6703059673309326
9 -1.882942e+03 5.342676e+01
* time: 1.8825349807739258
10 -1.888020e+03 2.010181e+01
* time: 2.0884039402008057
11 -1.889832e+03 1.867263e+01
* time: 2.3101279735565186
12 -1.891649e+03 1.668512e+01
* time: 2.505246877670288
13 -1.892615e+03 1.820701e+01
* time: 2.710768938064575
14 -1.893453e+03 1.745195e+01
* time: 2.916289806365967
15 -1.894760e+03 1.850174e+01
* time: 3.1224098205566406
16 -1.895647e+03 1.773939e+01
* time: 3.326152801513672
17 -1.896597e+03 1.143462e+01
* time: 3.5344769954681396
18 -1.897114e+03 9.720097e+00
* time: 3.7398529052734375
19 -1.897373e+03 6.054321e+00
* time: 3.9454619884490967
20 -1.897498e+03 3.985954e+00
* time: 4.134361982345581
21 -1.897571e+03 4.262464e+00
* time: 4.333568811416626
22 -1.897633e+03 4.010234e+00
* time: 4.535064935684204
23 -1.897714e+03 4.805375e+00
* time: 4.736138820648193
24 -1.897802e+03 3.508706e+00
* time: 4.937376976013184
25 -1.897865e+03 3.691477e+00
* time: 5.13475489616394
26 -1.897900e+03 2.982720e+00
* time: 5.313155889511108
27 -1.897928e+03 2.563790e+00
* time: 5.510150909423828
28 -1.897968e+03 3.261485e+00
* time: 5.70566987991333
29 -1.898013e+03 3.064690e+00
* time: 5.900296926498413
30 -1.898040e+03 1.636525e+00
* time: 6.081982851028442
31 -1.898051e+03 1.439997e+00
* time: 6.280400991439819
32 -1.898057e+03 1.436504e+00
* time: 6.474470853805542
33 -1.898069e+03 1.881529e+00
* time: 6.667884826660156
34 -1.898095e+03 3.253165e+00
* time: 6.850181818008423
35 -1.898142e+03 4.257942e+00
* time: 7.074469804763794
36 -1.898199e+03 3.685241e+00
* time: 7.273766994476318
37 -1.898245e+03 2.567364e+00
* time: 7.475966930389404
38 -1.898246e+03 2.561588e+00
* time: 7.763505935668945
39 -1.898251e+03 2.530898e+00
* time: 8.0169038772583
40 -1.898298e+03 2.673689e+00
* time: 8.225608825683594
41 -1.898300e+03 2.795097e+00
* time: 8.465049028396606
42 -1.898337e+03 3.699678e+00
* time: 8.750589847564697
43 -1.898435e+03 4.257964e+00
* time: 8.967885971069336
44 -1.898437e+03 4.246435e+00
* time: 9.266204833984375
45 -1.898448e+03 3.821693e+00
* time: 9.526943922042847
46 -1.898477e+03 2.785717e+00
* time: 9.73681092262268
47 -1.898477e+03 2.730327e+00
* time: 10.012123823165894
48 -1.898478e+03 2.516907e+00
* time: 10.26265001296997
49 -1.898479e+03 1.935688e+00
* time: 10.504621982574463
50 -1.898479e+03 1.258213e+00
* time: 10.755369901657104
51 -1.898480e+03 9.887950e-01
* time: 10.998169898986816
52 -1.898480e+03 1.281807e+01
* time: 11.227849006652832
53 -1.898480e+03 2.556060e-01
* time: 11.51050090789795
54 -1.898480e+03 1.118572e-01
* time: 11.765657901763916
55 -1.898480e+03 1.118585e-01
* time: 12.020663976669312
56 -1.898480e+03 1.118595e-01
* time: 12.291797876358032
57 -1.898480e+03 7.041556e-02
* time: 12.520057916641235
58 -1.898480e+03 6.966559e-02
* time: 12.69608187675476
59 -1.898480e+03 6.385055e-02
* time: 12.86320686340332
60 -1.898480e+03 7.499060e-02
* time: 13.031009912490845
61 -1.898480e+03 7.742415e-02
* time: 13.239248991012573
62 -1.898480e+03 7.742304e-02
* time: 13.466570854187012
63 -1.898480e+03 7.742255e-02
* time: 13.686180830001831
64 -1.898480e+03 7.742192e-02
* time: 13.909379005432129
65 -1.898480e+03 7.742191e-02
* time: 14.154387950897217
66 -1.898480e+03 7.742191e-02
* time: 14.389948844909668
67 -1.898480e+03 7.742191e-02
* time: 14.643712997436523
68 -1.898480e+03 7.979269e-02
* time: 14.875096797943115
69 -1.898480e+03 7.979256e-02
* time: 15.108335971832275
70 -1.898480e+03 7.979253e-02
* time: 15.371480941772461
71 -1.898480e+03 7.979253e-02
* time: 15.675747871398926
72 -1.898480e+03 7.979253e-02
* time: 15.964759826660156
73 -1.898480e+03 7.979253e-02
* time: 16.146328926086426FittedPumasModel
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.6192
tvv 11.378
tvvp 8.4522
tvq 1.3163
tvka 4.8925
Ω₁,₁ 0.13243
Ω₂,₂ 0.059668
Ω₃,₃ 0.41582
Ω₄,₄ 0.080661
Ω₅,₅ 0.25001
σ_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.81378 | 2.61918 |
| 2 | tvv | 11.0046 | 11.3783 |
| 3 | tvvp | 5.53998 | 8.45221 |
| 4 | tvq | 1.51591 | 1.31633 |
| 5 | tvka | 2.0 | 4.89247 |
| 6 | Ω₁,₁ | 0.102669 | 0.13243 |
| 7 | Ω₂,₂ | 0.0607755 | 0.0596677 |
| 8 | Ω₃,₃ | 1.20116 | 0.415822 |
| 9 | Ω₄,₄ | 0.423494 | 0.0806615 |
| 10 | Ω₅,₅ | 0.244732 | 0.250008 |
| 11 | σ_p | 0.0484049 | 0.0490976 |
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_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 = (; resolution = (600, 800)),
)
Finally looking at some individual plots for the same subjects as earlier:
fig_subject_fits2 = subject_fits(
res_inspect_2cmp_unfix_ka;
separate = true,
paginate = true,
facet = (; combinelabels = true, linkyaxes = false),
figure = (; fontsize = 18),
axis = (; xlabel = "Time (hr)", ylabel = "CTMx Concentration (ng/mL)"),
)
fig_subject_fits2[6]
The randomly sampled individual fits don’t seem good in some individuals, but we can evaluate this via a vpc to see how to go about.
3.4 Visual Predictive Checks (VPC)
We can now perform a vpc to check. The default plots provide a 80% prediction interval and a 95% simulated CI (shaded area) around each of the quantiles
pk_vpc = vpc(
pkfit_2cmp_unfix_ka,
200;
observations = [:Conc],
stratify_by = [:Dose],
ensemblealg = EnsembleThreads(), # multi-threading
)[ Info: Continuous VPCVisual Predictive Check
Type of VPC: Continuous VPC
Simulated populations: 200
Subjects in data: 40
Stratification variable(s): [:Dose]
Confidence level: 0.95
VPC lines: quantiles ([0.1, 0.5, 0.9])vpc_plot(
pk_2cmp,
pk_vpc;
rows = 1,
columns = 3,
figure = (; resolution = (1400, 1000), fontsize = 22),
axis = (;
xlabel = "Time (hr)",
ylabel = "Observed/Predicted\n CTMx Concentration (ng/mL)",
),
facet = (; combinelabels = true),
)
The visual predictive check suggests that the model captures the data well across all dose levels.
4 Additional Help
If you have questions regarding this tutorial, please post them on our discourse site.