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.016474008560180664
1 2.343899e+02 1.747348e+03
* time: 0.36257505416870117
2 9.696232e+01 1.198088e+03
* time: 0.36348915100097656
3 -7.818699e+01 5.538151e+02
* time: 0.36414408683776855
4 -1.234803e+02 2.462514e+02
* time: 0.36480116844177246
5 -1.372888e+02 2.067458e+02
* time: 0.3654670715332031
6 -1.410579e+02 1.162950e+02
* time: 0.3663029670715332
7 -1.434754e+02 5.632816e+01
* time: 0.3669471740722656
8 -1.453401e+02 7.859270e+01
* time: 0.36760711669921875
9 -1.498185e+02 1.455606e+02
* time: 0.36853814125061035
10 -1.534371e+02 1.303682e+02
* time: 0.3692209720611572
11 -1.563557e+02 5.975474e+01
* time: 0.3699159622192383
12 -1.575052e+02 9.308611e+00
* time: 0.370743989944458
13 -1.579357e+02 1.234484e+01
* time: 0.3713710308074951
14 -1.581874e+02 7.478196e+00
* time: 0.3719780445098877
15 -1.582981e+02 2.027162e+00
* time: 0.3726069927215576
16 -1.583375e+02 5.578262e+00
* time: 0.37340712547302246
17 -1.583556e+02 4.727050e+00
* time: 0.37400007247924805
18 -1.583644e+02 2.340173e+00
* time: 0.37459397315979004
19 -1.583680e+02 7.738100e-01
* time: 0.3753960132598877
20 -1.583696e+02 3.300689e-01
* time: 0.3759880065917969
21 -1.583704e+02 3.641985e-01
* time: 0.3765881061553955
22 -1.583707e+02 4.365901e-01
* time: 0.37717509269714355
23 -1.583709e+02 3.887800e-01
* time: 0.37798094749450684
24 -1.583710e+02 2.766977e-01
* time: 0.3785829544067383
25 -1.583710e+02 1.758029e-01
* time: 0.37918615341186523
26 -1.583710e+02 1.133947e-01
* time: 0.3800079822540283
27 -1.583710e+02 7.922544e-02
* time: 0.38063716888427734
28 -1.583710e+02 5.954998e-02
* time: 0.3815591335296631
29 -1.583710e+02 4.157079e-02
* time: 0.3822770118713379
30 -1.583710e+02 4.295447e-02
* time: 0.38298606872558594
31 -1.583710e+02 5.170753e-02
* time: 0.3837759494781494
32 -1.583710e+02 2.644383e-02
* time: 0.4224059581756592
33 -1.583710e+02 4.548993e-03
* time: 0.423267126083374
34 -1.583710e+02 2.501804e-02
* time: 0.42408108711242676
35 -1.583710e+02 3.763440e-02
* time: 0.4247009754180908
36 -1.583710e+02 3.206026e-02
* time: 0.425321102142334
37 -1.583710e+02 1.003698e-02
* time: 0.42595696449279785
38 -1.583710e+02 2.209089e-02
* time: 0.42659997940063477
39 -1.583710e+02 4.954172e-03
* time: 0.4272031784057617
40 -1.583710e+02 1.609373e-02
* time: 0.42801594734191895
41 -1.583710e+02 1.579802e-02
* time: 0.4287240505218506
42 -1.583710e+02 1.014113e-03
* time: 0.42972707748413086
43 -1.583710e+02 6.050644e-03
* time: 0.43060803413391113
44 -1.583710e+02 1.354412e-02
* time: 0.4312169551849365
45 -1.583710e+02 4.473248e-03
* time: 0.4318251609802246
46 -1.583710e+02 4.644735e-03
* time: 0.4324300289154053
47 -1.583710e+02 9.829910e-03
* time: 0.43302297592163086
48 -1.583710e+02 1.047561e-03
* time: 0.4336249828338623
49 -1.583710e+02 8.366895e-03
* time: 0.43424010276794434
50 -1.583710e+02 7.879055e-04
* time: 0.43486905097961426FittedPumasModel
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: 6.508827209472656e-5
1 -7.022088e+02 1.707063e+02
* time: 0.13785696029663086
2 -7.314067e+02 2.903269e+02
* time: 0.19669699668884277
3 -8.520591e+02 2.285888e+02
* time: 0.25473594665527344
4 -1.120191e+03 3.795410e+02
* time: 0.38806891441345215
5 -1.178784e+03 2.323978e+02
* time: 0.44876909255981445
6 -1.218320e+03 9.699907e+01
* time: 0.5057520866394043
7 -1.223641e+03 5.862105e+01
* time: 0.5603671073913574
8 -1.227620e+03 1.831403e+01
* time: 0.6150479316711426
9 -1.228381e+03 2.132323e+01
* time: 0.6667420864105225
10 -1.230098e+03 2.921228e+01
* time: 0.7204370498657227
11 -1.230854e+03 2.029662e+01
* time: 0.7731819152832031
12 -1.231116e+03 5.229097e+00
* time: 0.8244979381561279
13 -1.231179e+03 1.689232e+00
* time: 0.8745429515838623
14 -1.231187e+03 1.215379e+00
* time: 0.9231269359588623
15 -1.231188e+03 2.770380e-01
* time: 0.9690670967102051
16 -1.231188e+03 1.636653e-01
* time: 0.996406078338623
17 -1.231188e+03 2.701133e-01
* time: 1.0377490520477295
18 -1.231188e+03 3.163363e-01
* time: 1.0798170566558838
19 -1.231188e+03 1.505149e-01
* time: 1.1220951080322266
20 -1.231188e+03 2.484999e-02
* time: 1.1476020812988281
21 -1.231188e+03 8.446863e-04
* time: 1.1843700408935547FittedPumasModel
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.700920104980469e-5
1 -9.197817e+02 9.927951e+02
* time: 0.11124801635742188
2 -1.372640e+03 2.054986e+02
* time: 0.21738004684448242
3 -1.446326e+03 1.543987e+02
* time: 0.3207859992980957
4 -1.545570e+03 1.855028e+02
* time: 0.4361591339111328
5 -1.581449e+03 1.713157e+02
* time: 0.590252161026001
6 -1.639433e+03 1.257382e+02
* time: 0.686316967010498
7 -1.695964e+03 7.450539e+01
* time: 0.8039040565490723
8 -1.722243e+03 5.961044e+01
* time: 0.921605110168457
9 -1.736883e+03 7.320921e+01
* time: 1.0178210735321045
10 -1.753547e+03 7.501938e+01
* time: 1.1178109645843506
11 -1.764053e+03 6.185661e+01
* time: 1.2193760871887207
12 -1.778991e+03 4.831033e+01
* time: 1.3418869972229004
13 -1.791492e+03 4.943278e+01
* time: 1.455388069152832
14 -1.799847e+03 2.871410e+01
* time: 1.589142084121704
15 -1.805374e+03 7.520789e+01
* time: 1.711332082748413
16 -1.816260e+03 2.990621e+01
* time: 1.8279120922088623
17 -1.818252e+03 2.401915e+01
* time: 1.929955005645752
18 -1.822988e+03 2.587225e+01
* time: 2.0492351055145264
19 -1.824653e+03 1.550517e+01
* time: 2.14900803565979
20 -1.826074e+03 1.788927e+01
* time: 2.248828172683716
21 -1.826821e+03 1.888389e+01
* time: 2.352494955062866
22 -1.827900e+03 1.432840e+01
* time: 2.466905117034912
23 -1.828511e+03 9.422040e+00
* time: 2.5694069862365723
24 -1.828754e+03 5.363444e+00
* time: 2.6756041049957275
25 -1.828862e+03 4.916167e+00
* time: 2.7779250144958496
26 -1.829007e+03 4.695750e+00
* time: 2.931036949157715
27 -1.829358e+03 1.090245e+01
* time: 3.0617079734802246
28 -1.829830e+03 1.451320e+01
* time: 3.1672520637512207
29 -1.830201e+03 1.108694e+01
* time: 3.2763471603393555
30 -1.830360e+03 2.892320e+00
* time: 3.3947880268096924
31 -1.830390e+03 1.699267e+00
* time: 3.495103120803833
32 -1.830404e+03 1.602159e+00
* time: 3.5915441513061523
33 -1.830432e+03 2.822847e+00
* time: 3.6898021697998047
34 -1.830475e+03 4.105639e+00
* time: 3.7928719520568848
35 -1.830527e+03 5.093685e+00
* time: 3.9115381240844727
36 -1.830592e+03 2.697353e+00
* time: 4.018035173416138
37 -1.830615e+03 3.468352e+00
* time: 4.1215500831604
38 -1.830623e+03 2.594581e+00
* time: 4.225497007369995
39 -1.830625e+03 1.770110e+00
* time: 4.337977170944214
40 -1.830627e+03 1.040041e+00
* time: 4.430233001708984
41 -1.830628e+03 1.124112e+00
* time: 4.52751898765564
42 -1.830628e+03 3.547258e-01
* time: 4.619502067565918
43 -1.830629e+03 4.070652e-01
* time: 4.710325002670288
44 -1.830630e+03 7.177755e-01
* time: 4.804986000061035
45 -1.830630e+03 5.121219e-01
* time: 4.911261081695557
46 -1.830630e+03 9.584921e-02
* time: 5.001054048538208
47 -1.830630e+03 1.039362e-02
* time: 5.081762075424194
48 -1.830630e+03 5.788836e-03
* time: 5.159712076187134
49 -1.830630e+03 6.037375e-03
* time: 5.236691951751709
50 -1.830630e+03 7.187933e-03
* time: 5.313377141952515
51 -1.830630e+03 7.187933e-03
* time: 5.411728143692017
52 -1.830630e+03 7.187933e-03
* time: 5.510728120803833
53 -1.830630e+03 7.212102e-03
* time: 5.599266052246094
54 -1.830630e+03 7.214659e-03
* time: 5.7053680419921875
55 -1.830630e+03 7.217015e-03
* time: 5.794765949249268
56 -1.830630e+03 7.219581e-03
* time: 5.883728981018066
57 -1.830630e+03 7.222151e-03
* time: 5.972840070724487
58 -1.830630e+03 7.224724e-03
* time: 6.062087059020996
59 -1.830630e+03 7.227302e-03
* time: 6.153977155685425
60 -1.830630e+03 7.227561e-03
* time: 6.262532949447632
61 -1.830630e+03 7.227819e-03
* time: 6.365787982940674
62 -1.830630e+03 7.228077e-03
* time: 6.464324951171875
63 -1.830630e+03 7.228103e-03
* time: 6.561105966567993
64 -1.830630e+03 7.228129e-03
* time: 6.658561944961548
65 -1.830630e+03 7.228155e-03
* time: 6.7691810131073
66 -1.830630e+03 7.228180e-03
* time: 6.86415696144104
67 -1.830630e+03 7.228206e-03
* time: 6.960355997085571
68 -1.830630e+03 7.228209e-03
* time: 7.075157165527344
69 -1.830630e+03 7.228211e-03
* time: 7.1735711097717285
70 -1.830630e+03 7.228214e-03
* time: 7.271648168563843
71 -1.830630e+03 7.228217e-03
* time: 7.371112108230591
72 -1.830630e+03 7.228217e-03
* time: 7.487725019454956
73 -1.830630e+03 7.228217e-03
* time: 7.589398145675659
74 -1.830630e+03 7.228217e-03
* time: 7.6918439865112305
75 -1.830630e+03 7.228218e-03
* time: 7.808434963226318
76 -1.830630e+03 7.228218e-03
* time: 7.9098920822143555
77 -1.830630e+03 7.228218e-03
* time: 8.03056812286377
78 -1.830630e+03 7.228218e-03
* time: 8.136335134506226
79 -1.830630e+03 7.228218e-03
* time: 8.251532077789307
80 -1.830630e+03 7.214381e-03
* time: 8.335033178329468
81 -1.830630e+03 7.214381e-03
* time: 8.45929503440857
82 -1.830630e+03 7.214381e-03
* time: 8.580814123153687
83 -1.830630e+03 7.214381e-03
* time: 8.715829133987427
84 -1.830630e+03 7.214381e-03
* time: 8.853985071182251FittedPumasModel
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 | 8.854 | 1.184 |
| 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: 7.605552673339844e-5
1 -8.682982e+02 1.000199e+03
* time: 0.13599491119384766
2 -1.381870e+03 5.008081e+02
* time: 0.2697129249572754
3 -1.551053e+03 6.833490e+02
* time: 0.3884439468383789
4 -1.680887e+03 1.834586e+02
* time: 0.5140509605407715
5 -1.726118e+03 8.870274e+01
* time: 0.6200079917907715
6 -1.761023e+03 1.162036e+02
* time: 0.7532620429992676
7 -1.786619e+03 1.114552e+02
* time: 0.8607020378112793
8 -1.863556e+03 9.914305e+01
* time: 0.9983069896697998
9 -1.882942e+03 5.342676e+01
* time: 1.1105499267578125
10 -1.888020e+03 2.010181e+01
* time: 1.2267680168151855
11 -1.889832e+03 1.867263e+01
* time: 1.3595318794250488
12 -1.891649e+03 1.668512e+01
* time: 1.4751739501953125
13 -1.892615e+03 1.820701e+01
* time: 1.6027698516845703
14 -1.893453e+03 1.745195e+01
* time: 1.713608980178833
15 -1.894760e+03 1.850174e+01
* time: 1.8437080383300781
16 -1.895647e+03 1.773939e+01
* time: 1.9528310298919678
17 -1.896597e+03 1.143462e+01
* time: 2.0976200103759766
18 -1.897114e+03 9.720097e+00
* time: 2.206820011138916
19 -1.897373e+03 6.054321e+00
* time: 2.3219170570373535
20 -1.897498e+03 3.985954e+00
* time: 2.447504997253418
21 -1.897571e+03 4.262464e+00
* time: 2.5569310188293457
22 -1.897633e+03 4.010234e+00
* time: 2.683661937713623
23 -1.897714e+03 4.805375e+00
* time: 2.79059100151062
24 -1.897802e+03 3.508706e+00
* time: 2.9232399463653564
25 -1.897865e+03 3.691477e+00
* time: 3.029683828353882
26 -1.897900e+03 2.982720e+00
* time: 3.138169050216675
27 -1.897928e+03 2.563790e+00
* time: 3.2622878551483154
28 -1.897968e+03 3.261485e+00
* time: 3.3667449951171875
29 -1.898013e+03 3.064690e+00
* time: 3.477653980255127
30 -1.898040e+03 1.636525e+00
* time: 3.6012959480285645
31 -1.898051e+03 1.439997e+00
* time: 3.7081968784332275
32 -1.898057e+03 1.436504e+00
* time: 3.8361499309539795
33 -1.898069e+03 1.881529e+00
* time: 3.940809965133667
34 -1.898095e+03 3.253165e+00
* time: 4.0505828857421875
35 -1.898142e+03 4.257942e+00
* time: 4.175104856491089
36 -1.898199e+03 3.685241e+00
* time: 4.281255006790161
37 -1.898245e+03 2.567364e+00
* time: 4.413496017456055
38 -1.898246e+03 2.561588e+00
* time: 4.573415040969849
39 -1.898251e+03 2.530898e+00
* time: 4.728299856185913
40 -1.898298e+03 2.673689e+00
* time: 4.86061692237854
41 -1.898300e+03 2.795097e+00
* time: 4.991982936859131
42 -1.898337e+03 3.699678e+00
* time: 5.149942874908447
43 -1.898435e+03 4.257964e+00
* time: 5.288661003112793
44 -1.898437e+03 4.246435e+00
* time: 5.47176194190979
45 -1.898448e+03 3.821693e+00
* time: 5.618358850479126
46 -1.898477e+03 2.785717e+00
* time: 5.747454881668091
47 -1.898477e+03 2.730327e+00
* time: 5.893364906311035
48 -1.898478e+03 2.516907e+00
* time: 6.04137396812439
49 -1.898479e+03 1.935688e+00
* time: 6.201519012451172
50 -1.898479e+03 1.258213e+00
* time: 6.333572864532471
51 -1.898480e+03 9.887950e-01
* time: 6.478545904159546
52 -1.898480e+03 1.281807e+01
* time: 6.622905969619751
53 -1.898480e+03 2.556060e-01
* time: 6.7776939868927
54 -1.898480e+03 1.118572e-01
* time: 6.920666933059692
55 -1.898480e+03 1.118585e-01
* time: 7.070641040802002
56 -1.898480e+03 1.118595e-01
* time: 7.2318830490112305
57 -1.898480e+03 7.041556e-02
* time: 7.359047889709473
58 -1.898480e+03 6.966559e-02
* time: 7.476991891860962
59 -1.898480e+03 6.385055e-02
* time: 7.573364973068237
60 -1.898480e+03 7.499060e-02
* time: 7.66584587097168
61 -1.898480e+03 7.742415e-02
* time: 7.799831867218018
62 -1.898480e+03 7.742304e-02
* time: 7.920961856842041
63 -1.898480e+03 7.742255e-02
* time: 8.053854942321777
64 -1.898480e+03 7.742192e-02
* time: 8.174669981002808
65 -1.898480e+03 7.742191e-02
* time: 8.312906980514526
66 -1.898480e+03 7.742191e-02
* time: 8.458142042160034
67 -1.898480e+03 7.742191e-02
* time: 8.591567039489746
68 -1.898480e+03 7.979269e-02
* time: 8.728330850601196
69 -1.898480e+03 7.979256e-02
* time: 8.872411966323853
70 -1.898480e+03 7.979253e-02
* time: 9.010485887527466
71 -1.898480e+03 7.979253e-02
* time: 9.213069915771484
72 -1.898480e+03 7.979253e-02
* time: 9.370067834854126
73 -1.898480e+03 7.979253e-02
* time: 9.493891954421997FittedPumasModel
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.