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.023007869720458984
1 2.343899e+02 1.747348e+03
* time: 0.4593689441680908
2 9.696232e+01 1.198088e+03
* time: 0.4613628387451172
3 -7.818699e+01 5.538151e+02
* time: 0.4623379707336426
4 -1.234803e+02 2.462514e+02
* time: 0.4633939266204834
5 -1.372888e+02 2.067458e+02
* time: 0.4644429683685303
6 -1.410579e+02 1.162950e+02
* time: 0.4659128189086914
7 -1.434754e+02 5.632816e+01
* time: 0.5086348056793213
8 -1.453401e+02 7.859270e+01
* time: 0.5098237991333008
9 -1.498185e+02 1.455606e+02
* time: 0.5108017921447754
10 -1.534371e+02 1.303682e+02
* time: 0.5117688179016113
11 -1.563557e+02 5.975474e+01
* time: 0.5127229690551758
12 -1.575052e+02 9.308611e+00
* time: 0.5136549472808838
13 -1.579357e+02 1.234484e+01
* time: 0.5145819187164307
14 -1.581874e+02 7.478196e+00
* time: 0.5155220031738281
15 -1.582981e+02 2.027162e+00
* time: 0.5165379047393799
16 -1.583375e+02 5.578262e+00
* time: 0.5175778865814209
17 -1.583556e+02 4.727050e+00
* time: 0.5190989971160889
18 -1.583644e+02 2.340173e+00
* time: 0.5202310085296631
19 -1.583680e+02 7.738100e-01
* time: 0.5213217735290527
20 -1.583696e+02 3.300689e-01
* time: 0.5228948593139648
21 -1.583704e+02 3.641985e-01
* time: 0.5239639282226562
22 -1.583707e+02 4.365901e-01
* time: 0.5249907970428467
23 -1.583709e+02 3.887800e-01
* time: 0.5261509418487549
24 -1.583710e+02 2.766977e-01
* time: 0.5276229381561279
25 -1.583710e+02 1.758029e-01
* time: 0.528580904006958
26 -1.583710e+02 1.133947e-01
* time: 0.5295398235321045
27 -1.583710e+02 7.922544e-02
* time: 0.5309319496154785
28 -1.583710e+02 5.954998e-02
* time: 0.5318548679351807
29 -1.583710e+02 4.157079e-02
* time: 0.5328149795532227
30 -1.583710e+02 4.295447e-02
* time: 0.5342378616333008
31 -1.583710e+02 5.170753e-02
* time: 0.5352089405059814
32 -1.583710e+02 2.644383e-02
* time: 0.5365688800811768
33 -1.583710e+02 4.548993e-03
* time: 0.5383257865905762
34 -1.583710e+02 2.501804e-02
* time: 0.539543867111206
35 -1.583710e+02 3.763440e-02
* time: 0.5405049324035645
36 -1.583710e+02 3.206026e-02
* time: 0.5419049263000488
37 -1.583710e+02 1.003698e-02
* time: 0.5428338050842285
38 -1.583710e+02 2.209089e-02
* time: 0.5439188480377197
39 -1.583710e+02 4.954172e-03
* time: 0.5454268455505371
40 -1.583710e+02 1.609373e-02
* time: 0.5473079681396484
41 -1.583710e+02 1.579802e-02
* time: 0.5487899780273438
42 -1.583710e+02 1.014113e-03
* time: 0.5503959655761719
43 -1.583710e+02 6.050644e-03
* time: 0.5522170066833496
44 -1.583710e+02 1.354412e-02
* time: 0.5539968013763428
45 -1.583710e+02 4.473248e-03
* time: 0.5548789501190186
46 -1.583710e+02 4.644735e-03
* time: 0.588284969329834
47 -1.583710e+02 9.829910e-03
* time: 0.5896649360656738
48 -1.583710e+02 1.047561e-03
* time: 0.5906839370727539
49 -1.583710e+02 8.366895e-03
* time: 0.5916528701782227
50 -1.583710e+02 7.879055e-04
* time: 0.5926249027252197FittedPumasModel
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: 7.104873657226562e-5
1 -7.022088e+02 1.707063e+02
* time: 0.20298099517822266
2 -7.314067e+02 2.903269e+02
* time: 0.3285560607910156
3 -8.520591e+02 2.285888e+02
* time: 0.4285571575164795
4 -1.120191e+03 3.795410e+02
* time: 0.660750150680542
5 -1.178784e+03 2.323978e+02
* time: 0.7673070430755615
6 -1.218320e+03 9.699907e+01
* time: 0.8654439449310303
7 -1.223641e+03 5.862105e+01
* time: 0.9589660167694092
8 -1.227620e+03 1.831403e+01
* time: 1.0542960166931152
9 -1.228381e+03 2.132323e+01
* time: 1.1419241428375244
10 -1.230098e+03 2.921228e+01
* time: 1.2335381507873535
11 -1.230854e+03 2.029662e+01
* time: 1.3255980014801025
12 -1.231116e+03 5.229097e+00
* time: 1.4139759540557861
13 -1.231179e+03 1.689232e+00
* time: 1.485464096069336
14 -1.231187e+03 1.215379e+00
* time: 1.5669829845428467
15 -1.231188e+03 2.770380e-01
* time: 1.6419999599456787
16 -1.231188e+03 1.636653e-01
* time: 1.7093250751495361
17 -1.231188e+03 2.701133e-01
* time: 1.7767159938812256
18 -1.231188e+03 3.163363e-01
* time: 1.8306660652160645
19 -1.231188e+03 1.505149e-01
* time: 1.8950231075286865
20 -1.231188e+03 2.484999e-02
* time: 1.9575550556182861
21 -1.231188e+03 8.446863e-04
* time: 2.016414165496826FittedPumasModel
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: 8.106231689453125e-5
1 -9.197817e+02 9.927951e+02
* time: 0.24565410614013672
2 -1.372640e+03 2.054986e+02
* time: 0.4350099563598633
3 -1.446326e+03 1.543987e+02
* time: 0.6477639675140381
4 -1.545570e+03 1.855028e+02
* time: 0.8503780364990234
5 -1.581449e+03 1.713157e+02
* time: 1.1719300746917725
6 -1.639433e+03 1.257382e+02
* time: 1.3612439632415771
7 -1.695964e+03 7.450539e+01
* time: 1.556236982345581
8 -1.722243e+03 5.961044e+01
* time: 1.7617769241333008
9 -1.736883e+03 7.320921e+01
* time: 1.9566221237182617
10 -1.753547e+03 7.501938e+01
* time: 2.1680891513824463
11 -1.764053e+03 6.185661e+01
* time: 2.3746581077575684
12 -1.778991e+03 4.831033e+01
* time: 2.5966410636901855
13 -1.791492e+03 4.943278e+01
* time: 2.8218400478363037
14 -1.799847e+03 2.871410e+01
* time: 3.073413133621216
15 -1.805374e+03 7.520789e+01
* time: 3.3226089477539062
16 -1.816260e+03 2.990621e+01
* time: 3.5513241291046143
17 -1.818252e+03 2.401915e+01
* time: 3.769308090209961
18 -1.822988e+03 2.587225e+01
* time: 3.9771780967712402
19 -1.824653e+03 1.550517e+01
* time: 4.188854932785034
20 -1.826074e+03 1.788927e+01
* time: 4.410035133361816
21 -1.826821e+03 1.888389e+01
* time: 4.616788148880005
22 -1.827900e+03 1.432840e+01
* time: 4.814429998397827
23 -1.828511e+03 9.422040e+00
* time: 5.027518033981323
24 -1.828754e+03 5.363444e+00
* time: 5.238345146179199
25 -1.828862e+03 4.916167e+00
* time: 5.444849967956543
26 -1.829007e+03 4.695750e+00
* time: 5.644551038742065
27 -1.829358e+03 1.090245e+01
* time: 5.863084077835083
28 -1.829830e+03 1.451320e+01
* time: 6.075309991836548
29 -1.830201e+03 1.108694e+01
* time: 6.297608137130737
30 -1.830360e+03 2.892320e+00
* time: 6.509427070617676
31 -1.830390e+03 1.699267e+00
* time: 6.7196900844573975
32 -1.830404e+03 1.602159e+00
* time: 6.908579111099243
33 -1.830432e+03 2.822847e+00
* time: 7.098942995071411
34 -1.830475e+03 4.105639e+00
* time: 7.31061315536499
35 -1.830527e+03 5.093685e+00
* time: 7.513188123703003
36 -1.830592e+03 2.697353e+00
* time: 7.736372947692871
37 -1.830615e+03 3.468352e+00
* time: 7.939083099365234
38 -1.830623e+03 2.594581e+00
* time: 8.155632019042969
39 -1.830625e+03 1.770110e+00
* time: 8.34770393371582
40 -1.830627e+03 1.040041e+00
* time: 8.538861989974976
41 -1.830628e+03 1.124112e+00
* time: 8.726274967193604
42 -1.830628e+03 3.547258e-01
* time: 8.90208911895752
43 -1.830629e+03 4.070652e-01
* time: 9.073261022567749
44 -1.830630e+03 7.177755e-01
* time: 9.265002965927124
45 -1.830630e+03 5.121219e-01
* time: 9.44118595123291
46 -1.830630e+03 9.584921e-02
* time: 9.6127450466156
47 -1.830630e+03 1.039362e-02
* time: 9.769767999649048
48 -1.830630e+03 5.788836e-03
* time: 9.924812078475952
49 -1.830630e+03 6.037375e-03
* time: 10.064406156539917
50 -1.830630e+03 7.187933e-03
* time: 10.204432010650635
51 -1.830630e+03 7.187933e-03
* time: 10.401621103286743
52 -1.830630e+03 7.187933e-03
* time: 10.585646152496338
53 -1.830630e+03 7.212102e-03
* time: 10.762042045593262
54 -1.830630e+03 7.214659e-03
* time: 10.93213701248169
55 -1.830630e+03 7.217015e-03
* time: 11.101531982421875
56 -1.830630e+03 7.219581e-03
* time: 11.282581090927124
57 -1.830630e+03 7.222151e-03
* time: 11.453551054000854
58 -1.830630e+03 7.224724e-03
* time: 11.631006002426147
59 -1.830630e+03 7.227302e-03
* time: 11.818910121917725
60 -1.830630e+03 7.227561e-03
* time: 12.003399133682251
61 -1.830630e+03 7.227819e-03
* time: 12.20097303390503
62 -1.830630e+03 7.228077e-03
* time: 12.382674932479858
63 -1.830630e+03 7.228103e-03
* time: 12.560586929321289
64 -1.830630e+03 7.228129e-03
* time: 12.746459007263184
65 -1.830630e+03 7.228155e-03
* time: 12.93652892112732
66 -1.830630e+03 7.228180e-03
* time: 13.117639064788818
67 -1.830630e+03 7.228206e-03
* time: 13.301226139068604
68 -1.830630e+03 7.228209e-03
* time: 13.482129096984863
69 -1.830630e+03 7.228211e-03
* time: 13.662544012069702
70 -1.830630e+03 7.228214e-03
* time: 13.837778091430664
71 -1.830630e+03 7.228217e-03
* time: 14.031696081161499
72 -1.830630e+03 7.228217e-03
* time: 14.214507102966309
73 -1.830630e+03 7.228217e-03
* time: 14.391113996505737
74 -1.830630e+03 7.228217e-03
* time: 14.577691078186035
75 -1.830630e+03 7.228218e-03
* time: 14.761384010314941
76 -1.830630e+03 7.228218e-03
* time: 14.93655014038086
77 -1.830630e+03 7.228218e-03
* time: 15.130031108856201
78 -1.830630e+03 7.228218e-03
* time: 15.34204912185669
79 -1.830630e+03 7.228218e-03
* time: 15.566895008087158
80 -1.830630e+03 7.214381e-03
* time: 15.713245153427124
81 -1.830630e+03 7.214381e-03
* time: 15.927448987960815
82 -1.830630e+03 7.214381e-03
* time: 16.15268301963806
83 -1.830630e+03 7.214381e-03
* time: 16.375805139541626
84 -1.830630e+03 7.214381e-03
* time: 16.597398042678833FittedPumasModel
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 | 16.598 | 2.017 |
| 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.918212890625e-5
1 -8.682982e+02 1.000199e+03
* time: 0.31469011306762695
2 -1.381870e+03 5.008081e+02
* time: 0.5582900047302246
3 -1.551053e+03 6.833490e+02
* time: 0.8199670314788818
4 -1.680887e+03 1.834586e+02
* time: 1.071230173110962
5 -1.726118e+03 8.870274e+01
* time: 1.3435111045837402
6 -1.761023e+03 1.162036e+02
* time: 1.6267061233520508
7 -1.786619e+03 1.114552e+02
* time: 1.9000999927520752
8 -1.863556e+03 9.914305e+01
* time: 2.1820931434631348
9 -1.882942e+03 5.342676e+01
* time: 2.4329450130462646
10 -1.888020e+03 2.010181e+01
* time: 2.696493148803711
11 -1.889832e+03 1.867263e+01
* time: 2.9631760120391846
12 -1.891649e+03 1.668512e+01
* time: 3.2173070907592773
13 -1.892615e+03 1.820701e+01
* time: 3.474558115005493
14 -1.893453e+03 1.745195e+01
* time: 3.7264840602874756
15 -1.894760e+03 1.850174e+01
* time: 3.994093179702759
16 -1.895647e+03 1.773939e+01
* time: 4.230128049850464
17 -1.896597e+03 1.143462e+01
* time: 4.476845026016235
18 -1.897114e+03 9.720097e+00
* time: 4.728129148483276
19 -1.897373e+03 6.054321e+00
* time: 4.9801270961761475
20 -1.897498e+03 3.985954e+00
* time: 5.240940093994141
21 -1.897571e+03 4.262464e+00
* time: 5.503745079040527
22 -1.897633e+03 4.010234e+00
* time: 5.746827125549316
23 -1.897714e+03 4.805375e+00
* time: 5.984598159790039
24 -1.897802e+03 3.508706e+00
* time: 6.2309770584106445
25 -1.897865e+03 3.691477e+00
* time: 6.476570129394531
26 -1.897900e+03 2.982720e+00
* time: 6.715179204940796
27 -1.897928e+03 2.563790e+00
* time: 6.952192068099976
28 -1.897968e+03 3.261485e+00
* time: 7.184490203857422
29 -1.898013e+03 3.064690e+00
* time: 7.420509099960327
30 -1.898040e+03 1.636525e+00
* time: 7.674176216125488
31 -1.898051e+03 1.439997e+00
* time: 7.9276649951934814
32 -1.898057e+03 1.436504e+00
* time: 8.159252166748047
33 -1.898069e+03 1.881529e+00
* time: 8.410726070404053
34 -1.898095e+03 3.253165e+00
* time: 8.63944411277771
35 -1.898142e+03 4.257942e+00
* time: 8.883410215377808
36 -1.898199e+03 3.685241e+00
* time: 9.119894027709961
37 -1.898245e+03 2.567364e+00
* time: 9.38542103767395
38 -1.898246e+03 2.561588e+00
* time: 9.736809015274048
39 -1.898251e+03 2.530898e+00
* time: 10.035003185272217
40 -1.898298e+03 2.673689e+00
* time: 10.262229204177856
41 -1.898300e+03 2.795097e+00
* time: 10.562350034713745
42 -1.898337e+03 3.699678e+00
* time: 10.891691207885742
43 -1.898435e+03 4.257964e+00
* time: 11.136443138122559
44 -1.898437e+03 4.246435e+00
* time: 11.469350099563599
45 -1.898448e+03 3.821693e+00
* time: 11.776835203170776
46 -1.898477e+03 2.785717e+00
* time: 11.998000144958496
47 -1.898477e+03 2.730327e+00
* time: 12.299758195877075
48 -1.898478e+03 2.516907e+00
* time: 12.568533182144165
49 -1.898479e+03 1.935688e+00
* time: 12.861083030700684
50 -1.898479e+03 1.258213e+00
* time: 13.129094123840332
51 -1.898480e+03 9.887950e-01
* time: 13.39158010482788
52 -1.898480e+03 1.281807e+01
* time: 13.641710042953491
53 -1.898480e+03 2.556060e-01
* time: 13.932107210159302
54 -1.898480e+03 1.118572e-01
* time: 14.21106219291687
55 -1.898480e+03 1.118585e-01
* time: 14.549192190170288
56 -1.898480e+03 1.118595e-01
* time: 14.893805027008057
57 -1.898480e+03 7.041556e-02
* time: 15.146278142929077
58 -1.898480e+03 6.966559e-02
* time: 15.3181471824646
59 -1.898480e+03 6.385055e-02
* time: 15.507527112960815
60 -1.898480e+03 7.499060e-02
* time: 15.663819074630737
61 -1.898480e+03 7.742415e-02
* time: 15.876546144485474
62 -1.898480e+03 7.742304e-02
* time: 16.13555121421814
63 -1.898480e+03 7.742255e-02
* time: 16.382019996643066
64 -1.898480e+03 7.742192e-02
* time: 16.631980180740356
65 -1.898480e+03 7.742191e-02
* time: 16.89803719520569
66 -1.898480e+03 7.742191e-02
* time: 17.16580605506897
67 -1.898480e+03 7.742191e-02
* time: 17.40164804458618
68 -1.898480e+03 7.979269e-02
* time: 17.650196075439453
69 -1.898480e+03 7.979256e-02
* time: 17.929628133773804
70 -1.898480e+03 7.979253e-02
* time: 18.222684144973755
71 -1.898480e+03 7.979253e-02
* time: 18.56340718269348
72 -1.898480e+03 7.979253e-02
* time: 18.897646188735962
73 -1.898480e+03 7.979253e-02
* time: 19.111366033554077FittedPumasModel
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.