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 PharmaDatasets
2.2 Data Wrangling
We start by reading in the dataset and making some quick summaries.
If you want to learn more about data wrangling, don’t forget to check our Data Wrangling in Julia tutorials!
= dataset("pk_painrelief")
pkpain_df first(pkpain_df, 5)
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.
= @rsubset pkpain_df :Dose != "Placebo";
pkpain_noplb_df first(pkpain_noplb_df, 5)
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))
end
We also need to create an :amt
column:
@rtransform! pkpain_noplb_df :amt = :Time == 0 ? :Dose : missing
Now, we map the data variables to the read_nca
function that prepares the data for NCA analysis.
= read_nca(
pkpain_nca
pkpain_noplb_df;= :Subject,
id = :Time,
time = :amt,
amt = :Conc,
observations = [:Dose],
group = :route,
route )
NCAPopulation (120 subjects):
Group: [["Dose" => 5], ["Dose" => 20], ["Dose" => 80]]
Number of missing observations: 0
Number of blq observations: 0
Now 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.
= observations_vs_time(
f
pkpain_nca;= true,
paginate = (; xlabel = "Time (hr)", ylabel = "CTMNoPain Concentration (ng/mL)"),
axis = (; combinelabels = true),
facet
)1] f[
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,= (; fontsize = 22, size = (800, 1000)),
figure = "black",
color = 3,
linewidth = (; xlabel = "Time (hr)", ylabel = "CTMX Concentration (μg/mL)"),
axis = (; combinelabels = true, linkaxes = true),
facet )
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.
= run_nca(pkpain_nca; sigdigits = 3) pk_nca
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
:
= subject_fits(
f
pk_nca,= true,
paginate = (; xlabel = "Time (hr)", ylabel = "CTMX Concentration (μg/mL)"),
axis = (; combinelabels = true, linkaxes = true),
facet
)1] f[
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
.
= [:Dose] strata
1-element Vector{Symbol}:
:Dose
= [:cmax, :aucinf_obs] params
2-element Vector{Symbol}:
:cmax
:aucinf_obs
= summarize(pk_nca; stratify_by = strata, parameters = params) output
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,= :cmax,
parameter = (; xlabel = "Dose (mg)", ylabel = "Cₘₐₓ (ng/mL)"),
axis = (; fontsize = 18),
figure )
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
:
= NCA.DoseLinearityPowerModel(pk_nca, :cmax; level = 0.9) dp
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
end
Further, observations at time of dosing, i.e., when evid = 1
have to be missing
@rtransform! pkpain_noplb_df :Conc = :evid == 1 ? missing : :Conc
The 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.
= read_pumas(
pkpain_noplb
pkpain_noplb_df;= :Subject,
id = :Time,
time = :amt,
amt = [:Conc],
observations = [:Dose],
covariates = :evid,
evid = :cmt,
cmt )
Population
Subjects: 120
Covariates: Dose
Observations: Conc
Now that the data is transformed to a Population
of subjects, we can explore different models.
3.2.2 One-compartment model
If you are not familiar yet with the @model
blocks and syntax, please check our documentation.
= @model begin
pk_1cmp
@metadata begin
= "One Compartment Model"
desc = u"hr"
timeu end
@param begin
"""
Clearance (L/hr)
"""
∈ RealDomain(; lower = 0, init = 3.2)
tvcl """
Volume (L)
"""
∈ RealDomain(; lower = 0, init = 16.4)
tvv """
Absorption rate constant (h-1)
"""
∈ RealDomain(; lower = 0, init = 3.8)
tvka """
- ΩCL
- ΩVc
- ΩKa
"""
∈ PDiagDomain(init = [0.04, 0.04, 0.04])
Ω """
Proportional RUV
"""
∈ RealDomain(; lower = 0.0001, init = 0.2)
σ_p end
@random begin
~ MvNormal(Ω)
η end
@covariates begin
"""
Dose (mg)
"""
Doseend
@pre begin
= tvcl * exp(η[1])
CL = tvv * exp(η[2])
Vc = tvka * exp(η[3])
Ka end
@dynamics Depots1Central1
@derived begin
:= @. Central / Vc
cp """
CTMx Concentration (ng/mL)
"""
~ @. Normal(cp, abs(cp) * σ_p)
Conc end
end
┌ Warning: Covariate Dose is not used in the model.
└ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GJfPM/src/dsl/model_macro.jl:2958
PumasModel
Parameters: tvcl, tvv, tvka, Ω, σ_p
Random effects: η
Covariates: Dose
Dynamical system variables: Depot, Central
Dynamical system type: Closed form
Derived: Conc
Observed: Conc
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
= zero_randeffs(pk_1cmp, init_params(pk_1cmp), pkpain_noplb) etas
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.
= simobs(pk_1cmp, pkpain_noplb, init_params(pk_1cmp), etas) simpk_iparams
Simulated population (Vector{<:Subject})
Simulated subjects: 120
Simulated variables: Conc
sim_plot(
pk_1cmp,
simpk_iparams;= [:Conc],
observations = (; fontsize = 18),
figure = (;
axis = "Time (hr)",
xlabel = "Observed/Predicted \n CTMx Concentration (ng/mL)",
ylabel
), )
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.
= (; init_params(pk_1cmp)..., tvka = 2, tvv = 10) pkparam
(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,)
= simobs(pk_1cmp, pkpain_noplb, pkparam, etas) simpk_changedpars
Simulated population (Vector{<:Subject})
Simulated subjects: 120
Simulated variables: Conc
sim_plot(
pk_1cmp,
simpk_changedpars;= [:Conc],
observations = (; fontsize = 18),
figure = (
axis = "Time (hr)",
xlabel = "Observed/Predicted \n CTMx Concentration (ng/mL)",
ylabel
), )
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
= fit(pk_1cmp, pkpain_noplb, init_params(pk_1cmp), NaivePooled(); omegas = (:Ω,)) pkfit_np
[ 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.024713993072509766
1 2.343899e+02 1.747348e+03
* time: 0.8817849159240723
2 9.696232e+01 1.198088e+03
* time: 0.885037899017334
3 -7.818699e+01 5.538151e+02
* time: 0.8872499465942383
4 -1.234803e+02 2.462514e+02
* time: 0.8897709846496582
5 -1.372888e+02 2.067458e+02
* time: 0.8925299644470215
6 -1.410579e+02 1.162950e+02
* time: 0.8953390121459961
7 -1.434754e+02 5.632816e+01
* time: 0.8982620239257812
8 -1.453401e+02 7.859270e+01
* time: 0.900968074798584
9 -1.498185e+02 1.455606e+02
* time: 0.9031789302825928
10 -1.534371e+02 1.303682e+02
* time: 0.9052610397338867
11 -1.563557e+02 5.975474e+01
* time: 0.9073779582977295
12 -1.575052e+02 9.308611e+00
* time: 0.9095370769500732
13 -1.579357e+02 1.234484e+01
* time: 0.911628007888794
14 -1.581874e+02 7.478196e+00
* time: 0.9136559963226318
15 -1.582981e+02 2.027162e+00
* time: 0.9156849384307861
16 -1.583375e+02 5.578262e+00
* time: 0.9176769256591797
17 -1.583556e+02 4.727050e+00
* time: 0.9196810722351074
18 -1.583644e+02 2.340173e+00
* time: 0.9216759204864502
19 -1.583680e+02 7.738100e-01
* time: 0.9236700534820557
20 -1.583696e+02 3.300689e-01
* time: 0.9257779121398926
21 -1.583704e+02 3.641985e-01
* time: 0.9277820587158203
22 -1.583707e+02 4.365901e-01
* time: 0.929771900177002
23 -1.583709e+02 3.887800e-01
* time: 1.1037778854370117
24 -1.583710e+02 2.766977e-01
* time: 1.10697603225708
25 -1.583710e+02 1.758029e-01
* time: 1.1112360954284668
26 -1.583710e+02 1.133947e-01
* time: 1.1137669086456299
27 -1.583710e+02 7.922544e-02
* time: 1.1162428855895996
28 -1.583710e+02 5.954998e-02
* time: 1.1186621189117432
29 -1.583710e+02 4.157079e-02
* time: 1.121211051940918
30 -1.583710e+02 4.295447e-02
* time: 1.1237890720367432
31 -1.583710e+02 5.170753e-02
* time: 1.1265020370483398
32 -1.583710e+02 2.644383e-02
* time: 1.129951000213623
33 -1.583710e+02 4.548993e-03
* time: 1.13372802734375
34 -1.583710e+02 2.501804e-02
* time: 1.1370019912719727
35 -1.583710e+02 3.763440e-02
* time: 1.1394550800323486
36 -1.583710e+02 3.206026e-02
* time: 1.1418681144714355
37 -1.583710e+02 1.003698e-02
* time: 1.1442670822143555
38 -1.583710e+02 2.209089e-02
* time: 1.146756887435913
39 -1.583710e+02 4.954172e-03
* time: 1.149190902709961
40 -1.583710e+02 1.609373e-02
* time: 1.1524770259857178
41 -1.583710e+02 1.579802e-02
* time: 1.154918909072876
42 -1.583710e+02 1.014113e-03
* time: 1.1572670936584473
43 -1.583710e+02 6.050644e-03
* time: 1.1604349613189697
44 -1.583710e+02 1.354412e-02
* time: 1.162842035293579
45 -1.583710e+02 4.473248e-03
* time: 1.165208101272583
46 -1.583710e+02 4.644735e-03
* time: 1.16758394241333
47 -1.583710e+02 9.829910e-03
* time: 1.1700119972229004
48 -1.583710e+02 1.047561e-03
* time: 1.1724090576171875
49 -1.583710e+02 8.366895e-03
* time: 1.1747419834136963
50 -1.583710e+02 7.879055e-04
* time: 1.1771268844604492
FittedPumasModel
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)
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 fit
ting 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
loglikelihood
subject wise - check if there any influential subjects
Below, we are basically checking if the initial estimates for any subject are way off that we are unable to compute the initial loglikelihood
.
= []
lls for subj in pkpain_noplb
push!(lls, loglikelihood(pk_1cmp, subj, pkparam, FOCE()))
end
# the plot below is using native CairoMakie `hist`
hist(lls; bins = 10, normalization = :none, color = (:black, 0.5))
The distribution of the loglikelihood’s suggest no extreme outliers.
A more convenient way is to use the findinfluential
function that provides a list of k
top influential subjects by showing the normalized (minus) loglikelihood for each subject. As you can see below, the minus loglikelihood in the range of 16 agrees with the histogram plotted above.
= findinfluential(pk_1cmp, pkpain_noplb, pkparam, FOCE()) influential_subjects
120-element Vector{@NamedTuple{id::String, nll::Float64}}:
(id = "148", nll = 16.65965885684477)
(id = "135", nll = 16.648985190076335)
(id = "156", nll = 15.959069556607496)
(id = "159", nll = 15.441218240496484)
(id = "149", nll = 14.71513464411951)
(id = "88", nll = 13.09709837464614)
(id = "16", nll = 12.98228052193144)
(id = "61", nll = 12.652182902303679)
(id = "71", nll = 12.500330088085505)
(id = "59", nll = 12.241510254805235)
⋮
(id = "57", nll = -22.79767423253431)
(id = "93", nll = -22.836900711478208)
(id = "12", nll = -23.007742339519247)
(id = "123", nll = -23.292751843079234)
(id = "41", nll = -23.425412534960515)
(id = "99", nll = -23.535214841901112)
(id = "29", nll = -24.025959868383083)
(id = "52", nll = -24.164757842493685)
(id = "24", nll = -25.57209232565845)
3.2.2.2 FOCE
Now that we have a good handle on our data, lets go ahead and fit
a population model with FOCE
:
= fit(pk_1cmp, pkpain_noplb, pkparam, FOCE(); constantcoef = (; tvka = 2)) pkfit_1cmp
[ 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.604194641113281e-5
1 -7.022088e+02 1.707063e+02
* time: 0.3335421085357666
2 -7.314067e+02 2.903269e+02
* time: 0.9243202209472656
3 -8.520591e+02 2.285888e+02
* time: 1.0918550491333008
4 -1.120191e+03 3.795410e+02
* time: 1.4610600471496582
5 -1.178784e+03 2.323978e+02
* time: 1.6336250305175781
6 -1.218320e+03 9.699907e+01
* time: 3.4604082107543945
7 -1.223641e+03 5.862105e+01
* time: 3.5862600803375244
8 -1.227620e+03 1.831403e+01
* time: 3.7190370559692383
9 -1.228381e+03 2.132323e+01
* time: 3.8521080017089844
10 -1.230098e+03 2.921228e+01
* time: 3.9922091960906982
11 -1.230854e+03 2.029662e+01
* time: 4.190099000930786
12 -1.231116e+03 5.229097e+00
* time: 4.306979179382324
13 -1.231179e+03 1.689232e+00
* time: 4.418698072433472
14 -1.231187e+03 1.215379e+00
* time: 4.539494037628174
15 -1.231188e+03 2.770380e-01
* time: 4.691087007522583
16 -1.231188e+03 1.636653e-01
* time: 4.7735700607299805
17 -1.231188e+03 2.701133e-01
* time: 4.857354164123535
18 -1.231188e+03 3.163363e-01
* time: 4.950387001037598
19 -1.231188e+03 1.505149e-01
* time: 5.0461461544036865
20 -1.231188e+03 2.484999e-02
* time: 5.1495420932769775
21 -1.231188e+03 8.446863e-04
* time: 5.21833610534668
FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: 1231.188
Number of subjects: 120
Number of parameters: Fixed Optimized
1 6
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
-------------------
Estimate
-------------------
tvcl 3.1642
tvv 13.288
tvka 2.0
Ω₁,₁ 0.08494
Ω₂,₂ 0.048568
Ω₃,₃ 5.5811
σ_p 0.10093
-------------------
infer(pkfit_1cmp)
[ Info: Calculating: variance-covariance matrix.
[ Info: Done.
Asymptotic inference results using sandwich estimator
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: 1231.188
Number of subjects: 120
Number of parameters: Fixed Optimized
1 6
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
-------------------------------------------------------------------
Estimate SE 95.0% C.I.
-------------------------------------------------------------------
tvcl 3.1642 0.08662 [ 2.9944 ; 3.334 ]
tvv 13.288 0.27481 [12.749 ; 13.827 ]
tvka 2.0 NaN [ NaN ; NaN ]
Ω₁,₁ 0.08494 0.011022 [ 0.063338; 0.10654 ]
Ω₂,₂ 0.048568 0.0063502 [ 0.036122; 0.061014]
Ω₃,₃ 5.5811 1.2188 [ 3.1922 ; 7.97 ]
σ_p 0.10093 0.0057196 [ 0.089718; 0.11214 ]
-------------------------------------------------------------------
Notice that tvka
is fixed to 2 as we don’t have a lot of information before tmax
. From the results above, we see that the parameter precision for this model is reasonable.
3.2.3 Two-compartment model
Just to be sure, let’s fit a 2-compartment model and evaluate:
= @model begin
pk_2cmp
@param begin
"""
Clearance (L/hr)
"""
∈ RealDomain(; lower = 0, init = 3.2)
tvcl """
Central Volume (L)
"""
∈ RealDomain(; lower = 0, init = 16.4)
tvv """
Peripheral Volume (L)
"""
∈ RealDomain(; lower = 0, init = 10)
tvvp """
Distributional Clearance (L/hr)
"""
∈ RealDomain(; lower = 0, init = 2)
tvq """
Absorption rate constant (h-1)
"""
∈ RealDomain(; lower = 0, init = 1.3)
tvka """
- ΩCL
- ΩVc
- ΩKa
- ΩVp
- ΩQ
"""
∈ PDiagDomain(init = [0.04, 0.04, 0.04, 0.04, 0.04])
Ω """
Proportional RUV
"""
∈ RealDomain(; lower = 0.0001, init = 0.2)
σ_p end
@random begin
~ MvNormal(Ω)
η end
@covariates begin
"""
Dose (mg)
"""
Doseend
@pre begin
= tvcl * exp(η[1])
CL = tvv * exp(η[2])
Vc = tvka * exp(η[3])
Ka = tvvp * exp(η[4])
Vp = tvq * exp(η[5])
Q end
@dynamics Depots1Central1Periph1
@derived begin
:= @. Central / Vc
cp """
CTMx Concentration (ng/mL)
"""
~ @. Normal(cp, cp * σ_p)
Conc end
end
┌ Warning: Covariate Dose is not used in the model.
└ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GJfPM/src/dsl/model_macro.jl:2958
PumasModel
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: Conc
3.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: 6.985664367675781e-5
1 -9.197817e+02 9.927951e+02
* time: 0.32958984375
2 -1.372640e+03 2.054986e+02
* time: 0.6624188423156738
3 -1.446326e+03 1.543987e+02
* time: 0.9322929382324219
4 -1.545570e+03 1.855028e+02
* time: 1.2357759475708008
5 -1.581449e+03 1.713157e+02
* time: 1.7014610767364502
6 -1.639433e+03 1.257382e+02
* time: 1.9680500030517578
7 -1.695964e+03 7.450539e+01
* time: 2.270176887512207
8 -1.722243e+03 5.961044e+01
* time: 2.5457990169525146
9 -1.736883e+03 7.320921e+01
* time: 2.8763740062713623
10 -1.753547e+03 7.501938e+01
* time: 3.231657028198242
11 -1.764053e+03 6.185661e+01
* time: 3.526944875717163
12 -1.778991e+03 4.831033e+01
* time: 3.8545429706573486
13 -1.791492e+03 4.943278e+01
* time: 4.201708078384399
14 -1.799847e+03 2.871410e+01
* time: 4.522182941436768
15 -1.805374e+03 7.520791e+01
* time: 4.882382869720459
16 -1.816260e+03 2.990621e+01
* time: 5.233353853225708
17 -1.818252e+03 2.401915e+01
* time: 5.4972710609436035
18 -1.822988e+03 2.587225e+01
* time: 5.817925930023193
19 -1.824653e+03 1.550517e+01
* time: 6.0952370166778564
20 -1.826074e+03 1.788927e+01
* time: 6.386954069137573
21 -1.826821e+03 1.888389e+01
* time: 6.692451000213623
22 -1.827900e+03 1.432840e+01
* time: 6.958230972290039
23 -1.828511e+03 9.422041e+00
* time: 7.2681920528411865
24 -1.828754e+03 5.363442e+00
* time: 7.560417890548706
25 -1.828862e+03 4.916159e+00
* time: 7.850364923477173
26 -1.829007e+03 4.695755e+00
* time: 8.153811931610107
27 -1.829358e+03 1.090249e+01
* time: 8.431594848632812
28 -1.829830e+03 1.451325e+01
* time: 8.74743103981018
29 -1.830201e+03 1.108715e+01
* time: 9.106160879135132
30 -1.830360e+03 2.891223e+00
* time: 9.381738901138306
31 -1.830390e+03 1.695557e+00
* time: 9.683528900146484
32 -1.830404e+03 1.601712e+00
* time: 9.94815707206726
33 -1.830432e+03 2.823385e+00
* time: 10.237279891967773
34 -1.830477e+03 4.060617e+00
* time: 10.55091905593872
35 -1.830528e+03 5.133499e+00
* time: 10.823570966720581
36 -1.830593e+03 2.830970e+00
* time: 11.1370370388031
37 -1.830616e+03 3.342835e+00
* time: 11.420769929885864
38 -1.830622e+03 3.708884e+00
* time: 11.726306915283203
39 -1.830625e+03 2.062934e+00
* time: 12.048168897628784
40 -1.830627e+03 1.278569e+00
* time: 12.330615043640137
41 -1.830628e+03 1.832895e+00
* time: 12.66934084892273
42 -1.830628e+03 3.768840e-01
* time: 12.940047979354858
43 -1.830629e+03 3.152895e-01
* time: 13.20452094078064
44 -1.830630e+03 4.871060e-01
* time: 13.506587982177734
45 -1.830630e+03 3.110627e-01
* time: 13.74018907546997
46 -1.830630e+03 2.687758e-02
* time: 14.016929864883423
47 -1.830630e+03 4.694018e-03
* time: 14.211304903030396
48 -1.830630e+03 8.272969e-03
* time: 14.437211036682129
49 -1.830630e+03 8.249151e-03
* time: 14.661539077758789
50 -1.830630e+03 8.245562e-03
* time: 15.009922981262207
51 -1.830630e+03 8.240030e-03
* time: 15.364125967025757
52 -1.830630e+03 8.240030e-03
* time: 15.681519031524658
53 -1.830630e+03 8.240030e-03
* time: 16.040081024169922
FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: 1830.6304
Number of subjects: 120
Number of parameters: Fixed Optimized
1 10
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
-------------------
Estimate
-------------------
tvcl 2.8138
tvv 11.005
tvvp 5.54
tvq 1.5159
tvka 2.0
Ω₁,₁ 0.10267
Ω₂,₂ 0.060776
Ω₃,₃ 1.2012
Ω₄,₄ 0.4235
Ω₅,₅ 0.24473
σ_p 0.048405
-------------------
3.3 Comparing One- versus Two-compartment models
The 2-compartment model has a much lower objective function compared to the 1-compartment. Let’s compare the estimates from the 2 models using the compare_estimates
function.
compare_estimates(; pkfit_1cmp, pkfit_2cmp)
Row | parameter | pkfit_1cmp | pkfit_2cmp |
---|---|---|---|
String | Float64? | Float64? | |
1 | tvcl | 3.1642 | 2.81378 |
2 | tvv | 13.288 | 11.0046 |
3 | tvka | 2.0 | 2.0 |
4 | Ω₁,₁ | 0.0849405 | 0.102669 |
5 | Ω₂,₂ | 0.0485682 | 0.0607757 |
6 | Ω₃,₃ | 5.58107 | 1.20115 |
7 | σ_p | 0.100928 | 0.0484049 |
8 | tvvp | missing | 5.53998 |
9 | tvq | missing | 1.51591 |
10 | Ω₄,₄ | missing | 0.423495 |
11 | Ω₅,₅ | missing | 0.244731 |
We perform a likelihood ratio test to compare the two nested models. The test statistic and the \(p\)-value clearly indicate that a 2-compartment model should be preferred.
lrtest(pkfit_1cmp, pkfit_2cmp)
Statistic: 1200.0
Degrees of freedom: 4
P-value: 0.0
We 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)
end
Row | Metric | pk2cmp | pk1cmp |
---|---|---|---|
String | Any | Any | |
1 | Successful | true | true |
2 | Estimation Time | 16.04 | 5.218 |
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.
= inspect(pkfit_1cmp) res_inspect_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: FOCE
= inspect(pkfit_2cmp) res_inspect_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: FOCE
= goodness_of_fit(res_inspect_1cmp; figure = (; fontsize = 12)) gof_1cmp
= goodness_of_fit(res_inspect_2cmp; figure = (; fontsize = 12)) gof_2cmp
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.
= subject_fits(
fig_subject_fits
res_inspect_2cmp;= true,
separate = true,
paginate = (; combinelabels = true),
facet = (; fontsize = 18),
figure = (; xlabel = "Time (hr)", ylabel = "CTMx Concentration (ng/mL)"),
axis
)1] fig_subject_fits[
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;= [:Dose],
categorical = (; size = (600, 800)),
figure )
Clearly, our guess at tvka
seems off-target. Let’s try and estimate tvka
instead of fixing it to 2
:
= fit(pk_2cmp, pkpain_noplb, init_params(pk_2cmp), FOCE()) pkfit_2cmp_unfix_ka
[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 -3.200734e+02 1.272671e+03
* time: 7.796287536621094e-5
1 -8.682982e+02 1.000199e+03
* time: 0.4357719421386719
2 -1.381870e+03 5.008081e+02
* time: 0.7661659717559814
3 -1.551053e+03 6.833490e+02
* time: 1.1560449600219727
4 -1.680887e+03 1.834586e+02
* time: 1.5097808837890625
5 -1.726118e+03 8.870274e+01
* time: 1.807318925857544
6 -1.761023e+03 1.162036e+02
* time: 2.161576986312866
7 -1.786619e+03 1.114552e+02
* time: 2.4786858558654785
8 -1.863556e+03 9.914305e+01
* time: 2.8541688919067383
9 -1.882942e+03 5.342676e+01
* time: 3.193096876144409
10 -1.888020e+03 2.010181e+01
* time: 3.555824041366577
11 -1.889832e+03 1.867263e+01
* time: 3.9458370208740234
12 -1.891649e+03 1.668512e+01
* time: 4.260592937469482
13 -1.892615e+03 1.820701e+01
* time: 4.624940872192383
14 -1.893453e+03 1.745195e+01
* time: 4.947556972503662
15 -1.894760e+03 1.850174e+01
* time: 5.318979024887085
16 -1.895647e+03 1.773939e+01
* time: 5.642477035522461
17 -1.896597e+03 1.143462e+01
* time: 6.001314878463745
18 -1.897114e+03 9.720097e+00
* time: 6.329721927642822
19 -1.897373e+03 6.054321e+00
* time: 6.684862852096558
20 -1.897498e+03 3.985954e+00
* time: 7.054429054260254
21 -1.897571e+03 4.262464e+00
* time: 7.360589981079102
22 -1.897633e+03 4.010234e+00
* time: 7.714879989624023
23 -1.897714e+03 4.805375e+00
* time: 8.030707836151123
24 -1.897802e+03 3.508706e+00
* time: 8.386404037475586
25 -1.897865e+03 3.691477e+00
* time: 8.703810930252075
26 -1.897900e+03 2.982720e+00
* time: 9.046013832092285
27 -1.897928e+03 2.563790e+00
* time: 9.365079879760742
28 -1.897968e+03 3.261485e+00
* time: 9.716320991516113
29 -1.898013e+03 3.064690e+00
* time: 10.034571886062622
30 -1.898040e+03 1.636525e+00
* time: 10.372076034545898
31 -1.898051e+03 1.439997e+00
* time: 10.733641862869263
32 -1.898057e+03 1.436504e+00
* time: 11.015001058578491
33 -1.898069e+03 1.881529e+00
* time: 11.351459980010986
34 -1.898095e+03 3.253165e+00
* time: 11.641527891159058
35 -1.898142e+03 4.257942e+00
* time: 11.981464862823486
36 -1.898199e+03 3.685241e+00
* time: 12.290660858154297
37 -1.898245e+03 2.567364e+00
* time: 12.635010957717896
38 -1.898246e+03 2.561591e+00
* time: 13.150048971176147
39 -1.898251e+03 2.530888e+00
* time: 13.545409917831421
40 -1.898298e+03 2.673696e+00
* time: 13.921424865722656
41 -1.898300e+03 2.794639e+00
* time: 14.316861867904663
42 -1.898337e+03 3.751590e+00
* time: 14.80947995185852
43 -1.898421e+03 4.878407e+00
* time: 15.175278902053833
44 -1.898433e+03 4.391719e+00
* time: 15.568634986877441
45 -1.898437e+03 4.216518e+00
* time: 16.10071587562561
46 -1.898442e+03 4.108397e+00
* time: 16.63192105293274
47 -1.898446e+03 3.934902e+00
* time: 17.17805004119873
48 -1.898449e+03 3.769838e+00
* time: 17.62833595275879
49 -1.898450e+03 3.739486e+00
* time: 18.140362977981567
50 -1.898450e+03 3.712049e+00
* time: 18.703497886657715
51 -1.898457e+03 3.623436e+00
* time: 19.13950490951538
52 -1.898471e+03 2.668312e+00
* time: 19.50264883041382
53 -1.898479e+03 2.302438e+00
* time: 19.871415853500366
54 -1.898480e+03 2.386566e-01
* time: 20.18097186088562
55 -1.898480e+03 7.802040e-01
* time: 20.52798104286194
56 -1.898480e+03 7.369786e-01
* time: 20.980095863342285
57 -1.898480e+03 5.113191e-01
* time: 21.410794019699097
58 -1.898480e+03 3.067709e-01
* time: 21.6919949054718
59 -1.898480e+03 3.076791e-01
* time: 22.006651878356934
60 -1.898480e+03 3.102066e-01
* time: 22.34007501602173
61 -1.898480e+03 3.102066e-01
* time: 22.71593189239502
62 -1.898480e+03 3.102069e-01
* time: 23.22360396385193
63 -1.898480e+03 3.102071e-01
* time: 23.95078182220459
64 -1.898480e+03 3.102074e-01
* time: 24.6696879863739
65 -1.898480e+03 3.102076e-01
* time: 25.389315843582153
66 -1.898480e+03 3.102079e-01
* time: 26.11537504196167
67 -1.898480e+03 3.102081e-01
* time: 26.847025871276855
68 -1.898480e+03 3.102081e-01
* time: 27.580785036087036
69 -1.898480e+03 3.102081e-01
* time: 28.316622018814087
70 -1.898480e+03 3.102082e-01
* time: 29.049312829971313
71 -1.898480e+03 3.102082e-01
* time: 29.790189027786255
72 -1.898480e+03 3.102082e-01
* time: 30.540308952331543
73 -1.898480e+03 3.102102e-01
* time: 31.33962392807007
74 -1.898480e+03 3.102102e-01
* time: 31.920801877975464
75 -1.898480e+03 3.102096e-01
* time: 32.4870879650116
76 -1.898480e+03 3.102096e-01
* time: 33.047523021698
77 -1.898480e+03 3.125688e-01
* time: 33.54310584068298
78 -1.898480e+03 3.125640e-01
* time: 34.008383989334106
79 -1.898480e+03 3.125618e-01
* time: 34.50533699989319
80 -1.898480e+03 3.125615e-01
* time: 35.046995878219604
81 -1.898480e+03 3.125612e-01
* time: 35.54880905151367
82 -1.898480e+03 3.125610e-01
* time: 36.08809494972229
83 -1.898480e+03 3.125609e-01
* time: 36.69255185127258
84 -1.898480e+03 3.125604e-01
* time: 37.21844983100891
85 -1.898480e+03 3.125602e-01
* time: 37.7409508228302
86 -1.898480e+03 3.125602e-01
* time: 38.24379301071167
87 -1.898480e+03 3.125602e-01
* time: 38.98719382286072
88 -1.898480e+03 3.125602e-01
* time: 39.76678991317749
89 -1.898480e+03 3.125602e-01
* time: 40.54706692695618
90 -1.898480e+03 3.125602e-01
* time: 41.35654401779175
91 -1.898480e+03 3.125602e-01
* time: 42.15867900848389
92 -1.898480e+03 3.125602e-01
* time: 42.7078058719635
93 -1.898480e+03 3.125602e-01
* time: 43.281699895858765
94 -1.898480e+03 3.125602e-01
* time: 43.82082390785217
95 -1.898480e+03 1.387453e-01
* time: 44.18882083892822
96 -1.898480e+03 1.387453e-01
* time: 44.646087884902954
97 -1.898480e+03 1.387453e-01
* time: 45.26629090309143
98 -1.898480e+03 1.387453e-01
* time: 45.89601397514343
99 -1.898480e+03 1.387453e-01
* time: 46.26556086540222
FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: 1898.4797
Number of subjects: 120
Number of parameters: Fixed Optimized
0 11
Observation records: Active Missing
Conc: 1320 0
Total: 1320 0
-------------------
Estimate
-------------------
tvcl 2.6191
tvv 11.378
tvvp 8.453
tvq 1.3164
tvka 4.8926
Ω₁,₁ 0.13243
Ω₂,₂ 0.059669
Ω₃,₃ 0.41581
Ω₄,₄ 0.080679
Ω₅,₅ 0.24996
σ_p 0.049097
-------------------
compare_estimates(; pkfit_2cmp, pkfit_2cmp_unfix_ka)
Row | parameter | pkfit_2cmp | pkfit_2cmp_unfix_ka |
---|---|---|---|
String | Float64? | Float64? | |
1 | tvcl | 2.81378 | 2.6191 |
2 | tvv | 11.0046 | 11.3784 |
3 | tvvp | 5.53998 | 8.45297 |
4 | tvq | 1.51591 | 1.31637 |
5 | tvka | 2.0 | 4.89257 |
6 | Ω₁,₁ | 0.102669 | 0.132432 |
7 | Ω₂,₂ | 0.0607757 | 0.0596693 |
8 | Ω₃,₃ | 1.20115 | 0.415811 |
9 | Ω₄,₄ | 0.423495 | 0.0806789 |
10 | Ω₅,₅ | 0.244731 | 0.249956 |
11 | σ_p | 0.0484049 | 0.0490975 |
Let’s revaluate the goodness of fits and η distribution plots.
Not much change in the general gof
plots
= inspect(pkfit_2cmp_unfix_ka) res_inspect_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: FOCE
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;= [:Dose],
categorical = [:η₃],
ebes = (; size = (600, 800)),
figure )
Finally looking at some individual plots for the same subjects as earlier:
= subject_fits(
fig_subject_fits2
res_inspect_2cmp_unfix_ka;= true,
separate = true,
paginate = (; combinelabels = true, linkyaxes = false),
facet = (; fontsize = 18),
figure = (; xlabel = "Time (hr)", ylabel = "CTMx Concentration (ng/mL)"),
axis
)6] fig_subject_fits2[
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
= vpc(
pk_vpc
pkfit_2cmp_unfix_ka,200;
= [:Conc],
observations = [:Dose],
stratify_by = EnsembleThreads(), # multi-threading
ensemblealg )
[ Info: Continuous VPC
Visual Predictive Check
Type of VPC: Continuous VPC
Simulated populations: 200
Subjects in data: 120
Stratification variable(s): [:Dose]
Confidence level: 0.95
VPC lines: quantiles ([0.1, 0.5, 0.9])
vpc_plot(
pk_2cmp,
pk_vpc;= 1,
rows = 3,
columns = (; size = (1400, 1000), fontsize = 22),
figure = (;
axis = "Time (hr)",
xlabel = "Observed/Predicted\n CTMx Concentration (ng/mL)",
ylabel
),= (; combinelabels = true),
facet )
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.