A Comprehensive Introduction to Pumas

Authors

Vijay Ivaturi

Jose Storopoli

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:

  1. Data wrangling in Julia
  2. Exploratory analysis in Julia
  3. Continuous data non-linear mixed effects modeling in Pumas
  4. 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:

using Pumas
using PumasUtilities
using NCA
using NCAUtilities

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.

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))
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.

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: 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.

f = observations_vs_time(
    pkpain_nca;
    paginate = true,
    axis = (; xlabel = "Time (hr)", ylabel = "CTMNoPain Concentration (ng/mL)"),
    facet = (; combinelabels = true),
)
f[1]

An observations versus time profile for all subjects

Observations versus Time

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),
)

An observations versus time profile for all subjects in a summarized manner

Summary Observations versus Time

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]

Trend plot with observations for all individual subjects over time

Subject Fits

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}:
 :Dose
params = [:cmax, :aucinf_obs]
2-element Vector{Symbol}:
 :cmax
 :aucinf_obs
output = 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),
)

A violin plot for the Cmax distribution for each dose group

Cmax for each Dose Group

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)

A dose linearity power model plot for Cmax

Dose Linearity Plot

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)

A dose proportionality plot for Cmax

Dose Proportionality Plot
dose_vs_dose_normalized(pk_nca, :aucinf_obs)

A dose proportionality plot for AUC

Dose Proportionality Plot

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.

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: Conc

Now that the data is transformed to a Population of subjects, we can explore different models.

3.2.2 One-compartment model

Note

If you are not familiar yet with the @model blocks and syntax, please check our documentation.

pk_1cmp = @model begin

    @metadata begin
        desc = "One Compartment Model"
        timeu = u"hr"
    end

    @param begin
        """
        Clearance (L/hr)
        """
        tvcl  RealDomain(; lower = 0, init = 3.2)
        """
        Volume (L)
        """
        tvv  RealDomain(; lower = 0, init = 16.4)
        """
        Absorption rate constant (h-1)
        """
        tvka  RealDomain(; lower = 0, init = 3.8)
        """
          - ΩCL
          - ΩVc
          - ΩKa
        """
        Ω  PDiagDomain(init = [0.04, 0.04, 0.04])
        """
        Proportional RUV
        """
        σ_p  RealDomain(; lower = 0.0001, init = 0.2)
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates begin
        """
        Dose (mg)
        """
        Dose
    end

    @pre begin
        CL = tvcl * exp(η[1])
        Vc = tvv * exp(η[2])
        Ka = tvka * exp(η[3])
    end

    @dynamics Depots1Central1

    @derived begin
        cp := @. Central / Vc
        """
        CTMx Concentration (ng/mL)
        """
        Conc ~ @. Normal(cp, abs(cp) * σ_p)
    end

end
PumasModel
  Parameters: tvcl, tvv, tvka, Ω, σ_p
  Random effects: η
  Covariates: Dose
  Dynamical variables: Depot, Central
  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: Conc
sim_plot(
    pk_1cmp,
    simpk_iparams;
    observations = [:Conc],
    figure = (; fontsize = 18),
    axis = (;
        xlabel = "Time (hr)",
        ylabel = "Observed/Predicted \n CTMx Concentration (ng/mL)",
    ),
)

A simulated observations versus time plot overlaid with the scatter plot of the observed observations

Simulated Observations Plot

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: Conc
sim_plot(
    pk_1cmp,
    simpk_changedpars;
    observations = [:Conc],
    figure = (; fontsize = 18),
    axis = (
        xlabel = "Time (hr)",
        ylabel = "Observed/Predicted \n CTMx Concentration (ng/mL)",
    ),
)

A simulated observations versus time plot overlaid with the scatter plot of the observed observations

Simulated Observations Plot

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.019459962844848633
     1     2.343899e+02     1.747348e+03
 * time: 0.41583800315856934
     2     9.696232e+01     1.198088e+03
 * time: 0.4171018600463867
     3    -7.818699e+01     5.538151e+02
 * time: 0.41799187660217285
     4    -1.234803e+02     2.462514e+02
 * time: 0.41909003257751465
     5    -1.372888e+02     2.067458e+02
 * time: 0.42035794258117676
     6    -1.410579e+02     1.162950e+02
 * time: 0.42162203788757324
     7    -1.434754e+02     5.632816e+01
 * time: 0.42294788360595703
     8    -1.453401e+02     7.859270e+01
 * time: 0.4241809844970703
     9    -1.498185e+02     1.455606e+02
 * time: 0.425400972366333
    10    -1.534371e+02     1.303682e+02
 * time: 0.42666101455688477
    11    -1.563557e+02     5.975474e+01
 * time: 0.4278700351715088
    12    -1.575052e+02     9.308611e+00
 * time: 0.4291059970855713
    13    -1.579357e+02     1.234484e+01
 * time: 0.4303319454193115
    14    -1.581874e+02     7.478196e+00
 * time: 0.43155694007873535
    15    -1.582981e+02     2.027162e+00
 * time: 0.43279194831848145
    16    -1.583375e+02     5.578262e+00
 * time: 0.43403196334838867
    17    -1.583556e+02     4.727050e+00
 * time: 0.4352738857269287
    18    -1.583644e+02     2.340173e+00
 * time: 0.4365348815917969
    19    -1.583680e+02     7.738100e-01
 * time: 0.43784594535827637
    20    -1.583696e+02     3.300689e-01
 * time: 0.43922996520996094
    21    -1.583704e+02     3.641985e-01
 * time: 0.4405529499053955
    22    -1.583707e+02     4.365901e-01
 * time: 0.4418959617614746
    23    -1.583709e+02     3.887800e-01
 * time: 0.4432210922241211
    24    -1.583710e+02     2.766977e-01
 * time: 0.44457101821899414
    25    -1.583710e+02     1.758029e-01
 * time: 0.4459190368652344
    26    -1.583710e+02     1.133947e-01
 * time: 0.4824390411376953
    27    -1.583710e+02     7.922544e-02
 * time: 0.4833869934082031
    28    -1.583710e+02     5.954998e-02
 * time: 0.484234094619751
    29    -1.583710e+02     4.157079e-02
 * time: 0.48508405685424805
    30    -1.583710e+02     4.295447e-02
 * time: 0.4859309196472168
    31    -1.583710e+02     5.170753e-02
 * time: 0.48670101165771484
    32    -1.583710e+02     2.644383e-02
 * time: 0.4879450798034668
    33    -1.583710e+02     4.548993e-03
 * time: 0.48917198181152344
    34    -1.583710e+02     2.501804e-02
 * time: 0.490354061126709
    35    -1.583710e+02     3.763440e-02
 * time: 0.49118804931640625
    36    -1.583710e+02     3.206026e-02
 * time: 0.49208807945251465
    37    -1.583710e+02     1.003698e-02
 * time: 0.49289393424987793
    38    -1.583710e+02     2.209089e-02
 * time: 0.4936549663543701
    39    -1.583710e+02     4.954172e-03
 * time: 0.4944279193878174
    40    -1.583710e+02     1.609373e-02
 * time: 0.4954850673675537
    41    -1.583710e+02     1.579802e-02
 * time: 0.4962799549102783
    42    -1.583710e+02     1.014113e-03
 * time: 0.49739694595336914
    43    -1.583710e+02     6.050644e-03
 * time: 0.4984629154205322
    44    -1.583710e+02     1.354412e-02
 * time: 0.4992549419403076
    45    -1.583710e+02     4.473248e-03
 * time: 0.5000360012054443
    46    -1.583710e+02     4.644735e-03
 * time: 0.5007989406585693
    47    -1.583710e+02     9.829910e-03
 * time: 0.501579999923706
    48    -1.583710e+02     1.047561e-03
 * time: 0.5023620128631592
    49    -1.583710e+02     8.366895e-03
 * time: 0.5032360553741455
    50    -1.583710e+02     7.879055e-04
 * time: 0.5042920112609863
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:              NaivePooled
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 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), x_gap = 0)

A histogram of the individual loglikelihoods

Histogram of Loglikelihoods

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())
5-element Vector{NamedTuple{(:id, :nll), Tuple{String, Float64}}}:
 (id = "148", nll = 16.65965885684475)
 (id = "135", nll = 16.64898519007633)
 (id = "156", nll = 15.9590695566075)
 (id = "159", nll = 15.441218240496482)
 (id = "149", nll = 14.715134644119512)
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.009506225585938e-5
     1    -7.022088e+02     1.707063e+02
 * time: 0.1868741512298584
     2    -7.314067e+02     2.903269e+02
 * time: 0.2677011489868164
     3    -8.520591e+02     2.285888e+02
 * time: 0.3528010845184326
     4    -1.120191e+03     3.795410e+02
 * time: 0.504932165145874
     5    -1.178784e+03     2.323978e+02
 * time: 0.5758771896362305
     6    -1.218320e+03     9.699907e+01
 * time: 0.6426181793212891
     7    -1.223641e+03     5.862105e+01
 * time: 0.7062301635742188
     8    -1.227620e+03     1.831402e+01
 * time: 0.7696301937103271
     9    -1.228381e+03     2.132323e+01
 * time: 0.8294341564178467
    10    -1.230098e+03     2.921228e+01
 * time: 0.8918039798736572
    11    -1.230854e+03     2.029662e+01
 * time: 0.9690201282501221
    12    -1.231116e+03     5.229099e+00
 * time: 1.0285770893096924
    13    -1.231179e+03     1.689231e+00
 * time: 1.0871729850769043
    14    -1.231187e+03     1.215379e+00
 * time: 1.1435470581054688
    15    -1.231188e+03     2.770381e-01
 * time: 1.1796331405639648
    16    -1.231188e+03     1.636651e-01
 * time: 1.2269890308380127
    17    -1.231188e+03     2.701087e-01
 * time: 1.2760560512542725
    18    -1.231188e+03     3.163373e-01
 * time: 1.3244211673736572
    19    -1.231188e+03     1.505097e-01
 * time: 1.372650146484375
    20    -1.231188e+03     2.485456e-02
 * time: 1.4019861221313477
    21    -1.231188e+03     8.381370e-04
 * time: 1.443122148513794
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
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
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
end
PumasModel
  Parameters: tvcl, tvv, tvvp, tvq, tvka, Ω, σ_p
  Random effects: η
  Covariates: Dose
  Dynamical variables: Depot, Central, Peripheral
  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: 7.700920104980469e-5
     1    -9.197817e+02     9.927951e+02
 * time: 0.12532997131347656
     2    -1.372640e+03     2.054986e+02
 * time: 0.2588388919830322
     3    -1.446326e+03     1.543987e+02
 * time: 0.3678269386291504
     4    -1.545570e+03     1.855028e+02
 * time: 0.4759979248046875
     5    -1.581449e+03     1.713157e+02
 * time: 0.6520099639892578
     6    -1.639433e+03     1.257382e+02
 * time: 0.7558748722076416
     7    -1.695964e+03     7.450539e+01
 * time: 0.8623418807983398
     8    -1.722243e+03     5.961044e+01
 * time: 0.9806869029998779
     9    -1.736883e+03     7.320921e+01
 * time: 1.0837979316711426
    10    -1.753547e+03     7.501938e+01
 * time: 1.1941108703613281
    11    -1.764053e+03     6.185661e+01
 * time: 1.3044540882110596
    12    -1.778991e+03     4.831033e+01
 * time: 1.4509179592132568
    13    -1.791492e+03     4.943278e+01
 * time: 1.591062068939209
    14    -1.799847e+03     2.871410e+01
 * time: 1.753920078277588
    15    -1.805374e+03     7.520792e+01
 * time: 1.8918890953063965
    16    -1.816260e+03     2.990621e+01
 * time: 2.0262789726257324
    17    -1.818252e+03     2.401915e+01
 * time: 2.134427070617676
    18    -1.822988e+03     2.587225e+01
 * time: 2.247175931930542
    19    -1.824653e+03     1.550517e+01
 * time: 2.3696448802948
    20    -1.826074e+03     1.788927e+01
 * time: 2.501286029815674
    21    -1.826821e+03     1.888389e+01
 * time: 2.606563091278076
    22    -1.827900e+03     1.432840e+01
 * time: 2.716836929321289
    23    -1.828511e+03     9.422042e+00
 * time: 2.8403279781341553
    24    -1.828754e+03     5.363448e+00
 * time: 2.9545300006866455
    25    -1.828862e+03     4.916153e+00
 * time: 3.063361883163452
    26    -1.829007e+03     4.695757e+00
 * time: 3.18571400642395
    27    -1.829358e+03     1.090248e+01
 * time: 3.29803204536438
    28    -1.829830e+03     1.451324e+01
 * time: 3.411407947540283
    29    -1.830201e+03     1.108690e+01
 * time: 3.5418200492858887
    30    -1.830360e+03     2.892290e+00
 * time: 3.6551480293273926
    31    -1.830390e+03     1.698884e+00
 * time: 3.7673840522766113
    32    -1.830404e+03     1.602129e+00
 * time: 3.9007480144500732
    33    -1.830432e+03     2.823129e+00
 * time: 4.019996881484985
    34    -1.830475e+03     4.105407e+00
 * time: 4.129487991333008
    35    -1.830527e+03     5.093470e+00
 * time: 4.243844032287598
    36    -1.830592e+03     2.695875e+00
 * time: 4.371697902679443
    37    -1.830615e+03     3.470118e+00
 * time: 4.482125997543335
    38    -1.830623e+03     2.577088e+00
 * time: 4.594481945037842
    39    -1.830625e+03     1.765655e+00
 * time: 4.701636075973511
    40    -1.830627e+03     1.036472e+00
 * time: 4.816457033157349
    41    -1.830628e+03     1.124524e+00
 * time: 4.92054009437561
    42    -1.830628e+03     3.545350e-01
 * time: 5.018733024597168
    43    -1.830629e+03     4.125688e-01
 * time: 5.118269920349121
    44    -1.830630e+03     7.227016e-01
 * time: 5.231782913208008
    45    -1.830630e+03     5.071624e-01
 * time: 5.345696926116943
    46    -1.830630e+03     1.064027e-01
 * time: 5.442035913467407
    47    -1.830630e+03     5.753961e-03
 * time: 5.528441905975342
    48    -1.830630e+03     3.093135e-03
 * time: 5.610526084899902
    49    -1.830630e+03     2.462622e-03
 * time: 5.692231893539429
    50    -1.830630e+03     1.882491e-03
 * time: 5.772938966751099
    51    -1.830630e+03     1.882491e-03
 * time: 5.888859033584595
    52    -1.830630e+03     1.882491e-03
 * time: 6.013195991516113
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
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.423493
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.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
20×3 DataFrame
Row Metric pk2cmp pk1cmp
String Any Any
1 Successful true true
2 Estimation Time 6.013 1.443
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 Pumas.FOCE
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()
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()
gof_1cmp = goodness_of_fit(res_inspect_1cmp; figure = (; fontsize = 12))

A 4-mosaic goodness of fit plot showing the 1-compartment model

Goodness of Fit Plots
gof_2cmp = goodness_of_fit(res_inspect_2cmp; figure = (; fontsize = 12))

Trend plot with observations for all individual subjects over time

Subject Fits

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]

Trend plot with observations for 4 individual subjects over time

Subject Fits for 4 Individuals

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)

A histogram for the empirical Bayes distribution of all subject-specific parameters

Empirical Bayes Distribution
empirical_bayes_vs_covariates(
    res_inspect_2cmp;
    categorical = [:Dose],
    figure = (; resolution = (600, 800)),
)

A histogram for the empirical Bayes distribution of all subject-specific parameters stratified by categorical covariates

Empirical Bayes Distribution Stratified by Covariates

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: 0.0001049041748046875
     1    -8.682982e+02     1.000199e+03
 * time: 0.23105287551879883
     2    -1.381870e+03     5.008081e+02
 * time: 0.4362177848815918
     3    -1.551053e+03     6.833490e+02
 * time: 0.6362729072570801
     4    -1.680887e+03     1.834586e+02
 * time: 0.8405759334564209
     5    -1.726118e+03     8.870274e+01
 * time: 1.0202717781066895
     6    -1.761023e+03     1.162036e+02
 * time: 1.2339417934417725
     7    -1.786619e+03     1.114552e+02
 * time: 1.419010877609253
     8    -1.863556e+03     9.914305e+01
 * time: 1.6381158828735352
     9    -1.882942e+03     5.342676e+01
 * time: 1.8343398571014404
    10    -1.888020e+03     2.010181e+01
 * time: 2.0425679683685303
    11    -1.889832e+03     1.867263e+01
 * time: 2.255631923675537
    12    -1.891649e+03     1.668512e+01
 * time: 2.4641478061676025
    13    -1.892615e+03     1.820701e+01
 * time: 2.6557397842407227
    14    -1.893453e+03     1.745195e+01
 * time: 2.8603739738464355
    15    -1.894760e+03     1.850174e+01
 * time: 3.068542003631592
    16    -1.895647e+03     1.773939e+01
 * time: 3.2737669944763184
    17    -1.896597e+03     1.143462e+01
 * time: 3.463918924331665
    18    -1.897114e+03     9.720097e+00
 * time: 3.666883945465088
    19    -1.897373e+03     6.054321e+00
 * time: 3.8588268756866455
    20    -1.897498e+03     3.985954e+00
 * time: 4.065310955047607
    21    -1.897571e+03     4.262464e+00
 * time: 4.249583959579468
    22    -1.897633e+03     4.010234e+00
 * time: 4.450654983520508
    23    -1.897714e+03     4.805375e+00
 * time: 4.6342408657073975
    24    -1.897802e+03     3.508706e+00
 * time: 4.836848974227905
    25    -1.897865e+03     3.691477e+00
 * time: 5.018253803253174
    26    -1.897900e+03     2.982720e+00
 * time: 5.223111867904663
    27    -1.897928e+03     2.563790e+00
 * time: 5.39683198928833
    28    -1.897968e+03     3.261484e+00
 * time: 5.603308916091919
    29    -1.898013e+03     3.064690e+00
 * time: 5.776801824569702
    30    -1.898040e+03     1.636525e+00
 * time: 5.955432891845703
    31    -1.898051e+03     1.439997e+00
 * time: 6.148314952850342
    32    -1.898057e+03     1.436504e+00
 * time: 6.31996488571167
    33    -1.898069e+03     1.881528e+00
 * time: 6.515344858169556
    34    -1.898095e+03     3.253165e+00
 * time: 6.684129953384399
    35    -1.898142e+03     4.257942e+00
 * time: 6.889653921127319
    36    -1.898199e+03     3.685241e+00
 * time: 7.061599969863892
    37    -1.898245e+03     2.567364e+00
 * time: 7.2769269943237305
    38    -1.898246e+03     2.561590e+00
 * time: 7.585847854614258
    39    -1.898251e+03     2.530885e+00
 * time: 7.869495868682861
    40    -1.898298e+03     2.673694e+00
 * time: 8.071818828582764
    41    -1.898300e+03     2.794535e+00
 * time: 8.324241876602173
    42    -1.898337e+03     3.768309e+00
 * time: 8.616428852081299
    43    -1.898415e+03     5.140597e+00
 * time: 8.798900842666626
    44    -1.898427e+03     4.646687e+00
 * time: 9.039944887161255
    45    -1.898432e+03     4.441103e+00
 * time: 9.340512990951538
    46    -1.898437e+03     4.297436e+00
 * time: 9.624772787094116
    47    -1.898443e+03     4.170006e+00
 * time: 9.924351930618286
    48    -1.898448e+03     3.949707e+00
 * time: 10.216111898422241
    49    -1.898451e+03     3.720254e+00
 * time: 10.501424789428711
    50    -1.898451e+03     3.700859e+00
 * time: 10.786243915557861
    51    -1.898456e+03     3.967298e+00
 * time: 11.049617767333984
    52    -1.898479e+03     2.134121e+00
 * time: 11.246042966842651
    53    -1.898480e+03     7.166192e-01
 * time: 11.459459781646729
    54    -1.898480e+03     7.319824e-01
 * time: 11.64109492301941
    55    -1.898480e+03     1.358764e-01
 * time: 11.829197883605957
    56    -1.898480e+03     1.239822e-01
 * time: 11.99403190612793
    57    -1.898480e+03     1.239752e-01
 * time: 12.19646692276001
    58    -1.898480e+03     1.239711e-01
 * time: 12.484867811203003
    59    -1.898480e+03     1.239711e-01
 * time: 12.765661001205444
    60    -1.898480e+03     1.239711e-01
 * time: 13.101184844970703
    61    -1.898480e+03     1.248189e-01
 * time: 13.438716888427734
    62    -1.898480e+03     1.248189e-01
 * time: 13.767398834228516
    63    -1.898480e+03     1.248189e-01
 * time: 14.108152866363525
    64    -1.898480e+03     1.248189e-01
 * time: 14.480135917663574
    65    -1.898480e+03     1.248189e-01
 * time: 14.818104982376099
    66    -1.898480e+03     1.248189e-01
 * time: 15.198214769363403
    67    -1.898480e+03     1.248189e-01
 * time: 15.556413888931274
    68    -1.898480e+03     1.248189e-01
 * time: 15.901340961456299
    69    -1.898480e+03     1.248189e-01
 * time: 16.276745796203613
    70    -1.898480e+03     1.248189e-01
 * time: 16.521309852600098
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
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.4529
tvq       1.3164
tvka      4.8925
Ω₁,₁      0.13243
Ω₂,₂      0.05967
Ω₃,₃      0.41581
Ω₄,₄      0.080679
Ω₅,₅      0.24996
σ_p       0.049097
-------------------
compare_estimates(; pkfit_2cmp, pkfit_2cmp_unfix_ka)
11×3 DataFrame
Row parameter pkfit_2cmp pkfit_2cmp_unfix_ka
String Float64? Float64?
1 tvcl 2.81378 2.61912
2 tvv 11.0046 11.3783
3 tvvp 5.53998 8.45294
4 tvq 1.51591 1.31637
5 tvka 2.0 4.89251
6 Ω₁,₁ 0.102669 0.132432
7 Ω₂,₂ 0.0607755 0.0596696
8 Ω₃,₃ 1.20116 0.415809
9 Ω₄,₄ 0.423493 0.0806791
10 Ω₅,₅ 0.244732 0.249956
11 σ_p 0.0484049 0.0490975

Let’s revaluate the goodness of fits and η distribution plots.

Not much change in the general gof plots

res_inspect_2cmp_unfix_ka = inspect(pkfit_2cmp_unfix_ka)
[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating individual parameters.
[ Info: Evaluating dose control parameters.
[ Info: Done.
FittedPumasModelInspection

Fitting was successful: true
Likehood approximation used for weighted residuals : Pumas.FOCE()
goodness_of_fit(res_inspect_2cmp_unfix_ka; figure = (; fontsize = 12))

A 4-mosaic goodness of fit plot showing the 2-compartment model

Goodness of Fit Plots

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)),
)

A histogram for the empirical Bayes distribution of all subject-specific parameters stratified by categorical covariates

Empirical Bayes Distribution Stratified by Covariates

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]

Trend plot with observations for 4 individual subjects over time

Subject Fits for 4 Individuals

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 VPC
Visual 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),
)

A visual predictive plot stratified by dose group

Visual Predictive Plots

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.