Fitting and Inferring a Compartmental PK Model

Author

Vijay Ivaturi

1 Introduction

A typical workflow for fitting a Pumas model and deriving parameter precision typically involves:

  1. Preparing the data and the model.
  2. Checking model-data compatibility.
  3. Obtaining initial parameter estimates.
  4. Fitting the model via a chosen estimation method.
  5. Interpreting the fit results.
  6. Computing parameter precision.
  7. (Optionally) proceeding with more advanced techniques like bootstrapping or SIR for robust uncertainty quantification.

The following sections will walk through these steps using a two-compartment PK model for Warfarin, focusing on the PK aspects only. Exploratory data analysis (EDA), although extremely important, is out of scope here. Readers interested in EDA are encouraged to consult other tutorials.

2 Model and Data

2.1 Model Definition

Below is the PK model, named warfarin_pk_model, defined in Pumas. This model contains:

  • Fixed effects (population parameters): pop_CL, pop_Vc, pop_Q, pop_Vp, pop_tabs, pop_lag
  • Inter-individual variability (IIV) components: pk_Ω, lag_ω
  • Residual error model parameters: σ_prop, σ_add
  • Covariates for scaling: FSZCL and FSZV
  • Differential equations describing the PK behavior in the compartments.
using Pumas
using PharmaDatasets
using DataFramesMeta
using CairoMakie
using AlgebraOfGraphics
using PumasUtilities
using Markdown
warfarin_pk_model = @model begin

    @metadata begin
        desc = "Warfarin 2-compartment PK model (PD removed)"
        timeu = u"hr"
    end

    @param begin
        # PK parameters
        """
        Clearance (L/hr)
        """
        pop_CL  RealDomain(lower = 0.0, init = 0.134)
        """
        Central volume (L)
        """
        pop_Vc  RealDomain(lower = 0.0, init = 8.11)
        """
        Inter-compartmental clearance (L/hr)
        """
        pop_Q  RealDomain(lower = 0.0, init = 0.5)
        """
        Peripheral volume (L)
        """
        pop_Vp  RealDomain(lower = 0.0, init = 20.0)
        """
        Absorption lag time (hr)
        """
        pop_tabs  RealDomain(lower = 0.0, init = 0.523)
        """
        Lag time (hr)
        """
        pop_lag  RealDomain(lower = 0.0, init = 0.1)

        # Inter-individual variability
        """
          - ΩCL: Clearance
          - ΩVc: Central volume
          - Ωtabs: Absorption lag time
        """
        pk_Ω  PDiagDomain([0.01, 0.01, 0.01])
        # Residual variability
        """
        σ_prop: Proportional error
        """
        σ_prop  RealDomain(lower = 0.0, init = 0.00752)
        """
        σ_add: Additive error
        """
        σ_add  RealDomain(lower = 0.0, init = 0.0661)
    end

    @random begin
        pk_η ~ MvNormal(pk_Ω)    # mean = 0, covariance = pk_Ω
    end

    @covariates begin
        """
        FSZCL: Clearance scaling factor
        """
        FSZCL
        """
        FSZV: Volume scaling factor
        """
        FSZV
    end

    @pre begin
        CL = FSZCL * pop_CL * exp(pk_η[1])
        Vc = FSZV * pop_Vc * exp(pk_η[2])
        Q = FSZCL * pop_Q
        Vp = FSZV * pop_Vp

        tabs = pop_tabs * exp(pk_η[3])
        Ka = log(2) / tabs
    end

    @dosecontrol begin
        lags = (Depot = pop_lag,)
    end

    @vars begin
        cp := Central / Vc
    end

    @dynamics begin
        Depot' = -Ka * Depot
        Central' =
            Ka * Depot - (CL / Vc) * Central - (Q / Vc) * Central + (Q / Vp) * Peripheral
        Peripheral' = (Q / Vc) * Central - (Q / Vp) * Peripheral
    end

    @derived begin
        """
        Concentration (ng/mL)
        """
        conc ~ @. Normal(cp, sqrt((σ_prop * cp)^2 + σ_add^2))
    end
end
PumasModel
  Parameters: pop_CL, pop_Vc, pop_Q, pop_Vp, pop_tabs, pop_lag, pk_Ω, σ_prop, σ_add
  Random effects: pk_η
  Covariates: FSZCL, FSZV
  Dynamical system variables: Depot, Central, Peripheral
  Dynamical system type: Matrix exponential
  Derived: conc
  Observed: conc

2.2 Data Preparation

The Warfarin data used in this tutorial is pulled from PharmaDatasets for demonstration purposes. Note how the code reshapes and prepares the data in “wide” format for reading into Pumas. Only the conc column is treated as observations for the PK model.


warfarin_data = dataset("paganz2024/warfarin_long")

# Step 2: Fix Duplicate Time Points
# -------------------------------
# Some subjects have duplicate time points for DVID = 1
# For this dataset, the triple (ID, TIME, DVID) should define
# a row uniquely, but
nrow(warfarin_data)
nrow(unique(warfarin_data, ["ID", "TIME", "DVID"]))

# We can identify the problematic rows by grouping on the index variables
@chain warfarin_data begin
    @groupby :ID :TIME :DVID
    @transform :tmp = length(:ID)
    @rsubset :tmp > 1
end

# It is important to understand the reason for the duplicate values.
# Sometimes the duplication is caused by recording errors, sometimes
# it is a data processing error, e.g. when joining tables, or it can
# be genuine records, e.g. when samples have been analyzed in multiple
# labs. The next step depends on which of the causes are behind the
# duplications.
#
# In this case, we will assume that both values are informative and
# we will therefore just adjust the time stamp a bit for the second
# observation.
warfarin_data_processed = @chain warfarin_data begin
    @groupby :ID :TIME :DVID
    @transform :tmp = 1:length(:ID)
    @rtransform :TIME = :tmp == 1 ? :TIME : :TIME + 1e-6
    @select Not(:tmp)
end

# Transform the data in a single chain of operations
warfarin_data_wide = @chain warfarin_data_processed begin
    @rsubset !contains(:ID, "#")
    @rtransform begin
        # Scaling factors
        :FSZV = :WEIGHT / 70            # volume scaling
        :FSZCL = (:WEIGHT / 70)^0.75     # clearance scaling (allometric)
        # Column name for the DV
        :DVNAME = "DV$(:DVID)"
        # Dosing indicator columns
        :CMT = ismissing(:AMOUNT) ? missing : 1
        :EVID = ismissing(:AMOUNT) ? 0 : 1
    end
    unstack(Not([:DVID, :DVNAME, :DV]), :DVNAME, :DV)
    rename!(:DV1 => :conc, :DV2 => :pca)
end
317×13 DataFrame
292 rows omitted
Row ID TIME WEIGHT AGE SEX AMOUNT FSZV FSZCL CMT EVID DV0 pca conc
String3 Float64 Float64 Int64 Int64 Float64? Float64 Float64 Int64? Int64 Float64? Float64? Float64?
1 1 0.0 66.7 50 1 100.0 0.952857 0.96443 1 1 missing missing missing
2 1 0.0 66.7 50 1 missing 0.952857 0.96443 missing 0 missing 100.0 missing
3 1 24.0 66.7 50 1 missing 0.952857 0.96443 missing 0 missing 49.0 9.2
4 1 36.0 66.7 50 1 missing 0.952857 0.96443 missing 0 missing 32.0 8.5
5 1 48.0 66.7 50 1 missing 0.952857 0.96443 missing 0 missing 26.0 6.4
6 1 72.0 66.7 50 1 missing 0.952857 0.96443 missing 0 missing 22.0 4.8
7 1 96.0 66.7 50 1 missing 0.952857 0.96443 missing 0 missing 28.0 3.1
8 1 120.0 66.7 50 1 missing 0.952857 0.96443 missing 0 missing 33.0 2.5
9 2 0.0 66.7 31 1 100.0 0.952857 0.96443 1 1 missing missing missing
10 2 0.0 66.7 31 1 missing 0.952857 0.96443 missing 0 missing 100.0 missing
11 2 0.5 66.7 31 1 missing 0.952857 0.96443 missing 0 missing missing 0.0
12 2 2.0 66.7 31 1 missing 0.952857 0.96443 missing 0 missing missing 8.4
13 2 3.0 66.7 31 1 missing 0.952857 0.96443 missing 0 missing missing 9.7
306 31 48.0 83.3 24 1 missing 1.19 1.13936 missing 0 missing 24.0 6.4
307 31 72.0 83.3 24 1 missing 1.19 1.13936 missing 0 missing 22.0 4.5
308 31 96.0 83.3 24 1 missing 1.19 1.13936 missing 0 missing 28.0 3.4
309 31 120.0 83.3 24 1 missing 1.19 1.13936 missing 0 missing 42.0 2.5
310 32 0.0 62.0 21 1 93.0 0.885714 0.912999 1 1 missing missing missing
311 32 0.0 62.0 21 1 missing 0.885714 0.912999 missing 0 missing 100.0 missing
312 32 24.0 62.0 21 1 missing 0.885714 0.912999 missing 0 missing 36.0 8.9
313 32 36.0 62.0 21 1 missing 0.885714 0.912999 missing 0 missing 27.0 7.7
314 32 48.0 62.0 21 1 missing 0.885714 0.912999 missing 0 missing 24.0 6.9
315 32 72.0 62.0 21 1 missing 0.885714 0.912999 missing 0 missing 23.0 4.4
316 32 96.0 62.0 21 1 missing 0.885714 0.912999 missing 0 missing 20.0 3.5
317 32 120.0 62.0 21 1 missing 0.885714 0.912999 missing 0 missing 22.0 2.5

2.3 Creating a Pumas Population

Below is the creation of a population object in Pumas using read_pumas. Only the conc data are treated as the observation variable:

pop = read_pumas(
    warfarin_data_wide;
    id = :ID,
    time = :TIME,
    amt = :AMOUNT,
    cmt = :CMT,
    evid = :EVID,
    covariates = [:SEX, :WEIGHT, :FSZV, :FSZCL],
    observations = [:conc],
)
Population
  Subjects: 31
  Covariates: SEX, WEIGHT, FSZV, FSZCL
  Observations: conc
Note

The same data can contain multiple endpoints or PD observations. In this tutorial, the focus is solely on PK fitting. PKPD modeling on this warfarin dataset will be introduced later.

2.4 Checking Model-Data Compatibility

Before performing any fit, it is recommended to verify whether the defined model is consistent with the provided dataset. Pumas offers functions such as loglikelihood and findinfluential for these checks.

2.4.1 The loglikelihood Function

The loglikelihood function computes the log-likelihood of the model given data and parameters. In Pumas, the signature typically looks like:

loglikelihood(model, population, params, approx)

where:

  • model: The Pumas model definition (e.g., warfarin_pk_model).
  • population: A Pumas population object (e.g., pop).
  • params: A named tuple or dictionary containing parameter values.
  • approx: The approximation method to use. Common options include FOCE(), FO(), LaplaceI(), etc.

For example, one might write:

# A named tuple of parameter values
param_vals = (
    pop_CL = 0.134,
    pop_Vc = 8.11,
    pop_Q = 0.5,
    pop_Vp = 20.0,
    pop_tabs = 0.523,
    pop_lag = 0.1,
    pk_Ω = Diagonal([0.01, 0.01, 0.01]),
    σ_prop = 0.00752,
    σ_add = 0.0661,
)

ll_value = loglikelihood(warfarin_pk_model, pop, param_vals, FOCE())
-91123.78983612436

The initial loglikelihood of the warfarin PK model given the data and parameter values is `LL = -91123.78983612436 `.

If the model and data are incompatible (e.g., missing doses for compartments, or out-of-range parameter values), loglikelihood might return an error or produce a warning. A successful computation is a good sign that the model can handle the data.

2.4.2 The findinfluential Function

The findinfluential function helps identify observations that disproportionately influence the fit. It can be used before or after fitting, but even with initial guesses, it can highlight potentially problematic data points. The most notable case is when the loglikelihood cannot be evaluated in which case the data returned for the problematic subject does not return a finite value.

Typical usage is:

findinfluential(warfarin_pk_model, pop, params; approx)

This will list the individual loglikelihoods for each subject. Readers can then inspect these observations further, potentially re-checking the data quality or adjusting the model assumptions if needed for each subject.

fi = findinfluential(warfarin_pk_model, pop, param_vals, FOCE())
31-element Vector{@NamedTuple{id::String, nll::Float64}}:
 (id = "12", nll = 7035.6308301243125)
 (id = "8", nll = 6811.854099061583)
 (id = "25", nll = 5140.036812994075)
 (id = "27", nll = 4720.678341527095)
 (id = "11", nll = 3856.5709209433408)
 (id = "16", nll = 3853.3494408795573)
 (id = "22", nll = 3616.0674596218623)
 (id = "7", nll = 3592.014215601264)
 (id = "3", nll = 3434.8872372376954)
 (id = "26", nll = 3358.499378697223)
 ⋮
 (id = "18", nll = 2084.1283016469442)
 (id = "9", nll = 1950.7850340954703)
 (id = "24", nll = 1717.087929029517)
 (id = "29", nll = 1381.3498671515454)
 (id = "17", nll = 1350.9070098668137)
 (id = "6", nll = 1319.2234098688702)
 (id = "21", nll = 1072.85052503673)
 (id = "19", nll = 989.3840625094336)
 (id = "5", nll = 984.368482951681)

The default output is a vector of NamedTuples, which can easily be converted into a DataFrame for further analysis. One can then plot the distribution of the individual loglikelihoods to identify any potential problematic subjects.

fidf = DataFrame(fi)
hist(
    fidf.nll,
    axis = (xlabel = "Log-likelihood", ylabel = "Frequency"),
    color = (:blue, 0.5),
)

The plot above shows that there are some subjects with very high initial log-likelihoods, which might be worth investigating. Notice, that a high initial log-likelihood does not necessarily imply a model or data problem. It can simply be that for the chosen initial population parameters the model does not fit the observations well for the subject in question.

2.5 Getting Good Initial Estimates Using Naive Pooled

A reliable starting point for nonlinear mixed-effects modeling is to get reasonable initial parameter estimates. Pumas provides approaches such as:

  • Naive Pooled: Treats all data as if there were no inter-individual variability.
  • NCA (Non-compartmental Analysis): Uses PK metrics (CL, V, etc.) from a non-compartmental approach to seed the parameter values (not showcased in this tutorial).

Here, the Naive Pooled approach is demonstrated. This method estimates population parameters by ignoring inter-individual differences, effectively pooling data as if from a single “super-subject.”

Below is the signature and documentation for the Pumas fit function. For more details, see the Fitting in Pumas Documentation.

fit(
  model::PumasModel,
  population::Population,
  param::NamedTuple,
  approx::Union{LikelihoodApproximation, MAP};
  optim_alg = Optim.BFGS(linesearch = Optim.LineSearches.BackTracking(), initial_stepnorm = 1.0),
  optim_options::NamedTuple = NamedTuple(),
  optimize_fn = nothing,
  constantcoef::Tuple = (),
  ensemblealg::SciMLBase.EnsembleAlgorithm = EnsembleThreads(),
  checkidentification = true,
  diffeq_options = NamedTuple(),
  verbose = true,
  ignore_numerical_error = true,
)

Fit the Pumas model model to the dataset population with starting values param using the estimation method approx. Currently supported values for the `approx` argument are `FO`, `FOCE`,
`LaplaceI`, `NaivePooled`, and `BayesMCMC`. See the online documentation for more details about the different methods.

To control the optimization procedure for finding population parameters (fixed effects), use the `optim_alg` and `optim_options` keyword arguments. In previous versions of Pumas the argument
`optimize_fn` was used, but is now discouraged and will be removed in a later version of Pumas. These options control the optimization of all `approx` methods except BayesMCMC. The default
optimization function uses the quasi-Newton routine BFGS method from the `Optim` package. It can be changed by setting the `optim_alg` to an algorithm implemented in Optim.jl as long as it
does not use second order derivatives. Optimization specific options can be passed in using the `optim_options` keyword and has to be a `NamedTuple` with names and values that match the
`Optim.Options` type. For example, the optimization trace can be disabled and the algorithm can be changed to L-BFGS by passing `optim_alg=Optim.LBFGS(), optim_options = ;(show_trace=false)`
to fit. See [Optim](https://docs.pumas.ai/stable/basics/estimation/#Optimization-option) for more defails.

It is possible to fix one or more parameters of the fit by passing a Tuple of Symbols as the `constantcoef` argument with elements corresponding to the names of the fixed parameters, e.g.
`constantcoef=(:σ,)`.

When models include an `@random` block and fitting with NaivePooled is requested, it is required that the user sets all variability parameters to zero with `constantcoef` such that these can
be ignored in the optimization, e.g. `constantcoef = (:Ω,)` while overwriting the corresponding values in `param` with `(; init_params..., Ω = zeros(3, 3))`.


Parallelization of the optimization is supported for most estimation methods via the ensemble interface of DifferentialEquations.jl. Currently, the only supported options are:

    •  `EnsembleThreads()`: the default. Accelerate fits by using multiple threads.

    •  `EnsembleSerial()`: fit using a single thread.

    •  `EnsembleDistributed()`: fit by using multiple worker processes.

The `fit` function will check if any gradients and throw an exception if any of the elements are exactly zero unless `checkidentification` is set to false.

Further keyword arguments can be passed via the `diffeq_options` argument. This allows for passing arguments to the differential equations solver such as `alg`, `abstol`, and `reltol`. The
default values for these are `AutoVern7(Rodas5P(autodiff=true)), 1e-12, and 1e-8` respectively. See the [documentation](https://docs.pumas.ai/stable/basics/simulation/#Options-and-settings-for-simobs) for more details.

The keyword `verbose` controls if info statements about initial evaluations of the loglikelihood function and gradient should be printed or not. Defaults to true.

Since the numerical optimization routine can sometimes visit extreme regions of the parameter space during it's exploration we automatically handle the situation where the input
parameters result in errors when solving the dynamical system. This allows the algorithm to recover and continue in many situations that would have otherwise stalled early. Sometimes, it
is useful to turn this error handling off when debugging a model fit, and this can be done by setting `ignore_numerical_error = false`.

Below is how to run a Naive Pooled analysis:

naive_fit = fit(warfarin_pk_model, pop, param_vals, NaivePooled(); constantcoef = (:pk_Ω,))
[ 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.188875e+05     5.217691e+05
 * time: 0.026215076446533203
     1     3.248274e+04     5.854052e+04
 * time: 1.0637309551239014
     2     2.414204e+04     4.403968e+04
 * time: 1.0672879219055176
     3     9.800450e+03     1.834949e+04
 * time: 1.0707130432128906
     4     5.218319e+03     9.776442e+03
 * time: 1.0741140842437744
     5     2.657496e+03     4.805323e+03
 * time: 1.0775010585784912
     6     1.508416e+03     2.469328e+03
 * time: 1.0810329914093018
     7     9.524014e+02     1.261463e+03
 * time: 1.0844430923461914
     8     6.999017e+02     6.470385e+02
 * time: 1.0878310203552246
     9     5.922734e+02     3.262900e+02
 * time: 1.091162919998169
    10     5.535691e+02     2.252527e+02
 * time: 1.094351053237915
    11     5.434957e+02     1.979047e+02
 * time: 1.0975310802459717
    12     5.419714e+02     1.852856e+02
 * time: 1.1006379127502441
    13     5.417601e+02     1.809508e+02
 * time: 1.103679895401001
    14     5.415227e+02     1.766418e+02
 * time: 1.1067509651184082
    15     5.408349e+02     1.663266e+02
 * time: 1.1099491119384766
    16     5.392823e+02     1.469356e+02
 * time: 1.1131150722503662
    17     5.357939e+02     1.101891e+02
 * time: 1.1163370609283447
    18     5.295274e+02     6.101717e+01
 * time: 1.119739055633545
    19     5.227089e+02     4.376424e+01
 * time: 1.1231749057769775
    20     5.205584e+02     6.762423e+01
 * time: 1.126629114151001
    21     5.200667e+02     6.460758e+01
 * time: 1.1302299499511719
    22     5.198970e+02     5.159635e+01
 * time: 1.1338250637054443
    23     5.198793e+02     5.168473e+01
 * time: 1.1372630596160889
    24     5.182416e+02     5.493697e+01
 * time: 1.1407029628753662
    25     5.154632e+02     5.617036e+01
 * time: 1.144129991531372
    26     5.085379e+02     6.965501e+01
 * time: 1.1477079391479492
    27     4.992897e+02     8.613012e+01
 * time: 1.1510200500488281
    28     4.952284e+02     6.666806e+01
 * time: 1.154336929321289
    29     4.937816e+02     4.845720e+01
 * time: 1.1576409339904785
    30     4.933235e+02     4.622865e+01
 * time: 1.1609349250793457
    31     4.932968e+02     4.584617e+01
 * time: 1.1642329692840576
    32     4.932703e+02     4.565222e+01
 * time: 1.167543888092041
    33     4.931155e+02     4.490937e+01
 * time: 1.1709809303283691
    34     4.927970e+02     4.552021e+01
 * time: 1.1743929386138916
    35     4.919145e+02     5.125228e+01
 * time: 1.1777949333190918
    36     4.898542e+02     5.869445e+01
 * time: 1.1814560890197754
    37     4.854151e+02     6.533002e+01
 * time: 1.1850640773773193
    38     4.786263e+02     5.894401e+01
 * time: 1.1885490417480469
    39     4.723401e+02     3.679497e+01
 * time: 1.1920900344848633
    40     4.676138e+02     2.454358e+01
 * time: 1.1956191062927246
    41     4.667006e+02     1.298513e+01
 * time: 1.1990959644317627
    42     4.666338e+02     1.353668e+01
 * time: 1.2026050090789795
    43     4.666323e+02     1.384287e+01
 * time: 1.206115961074829
    44     4.666314e+02     1.374470e+01
 * time: 1.2096099853515625
    45     4.666256e+02     1.344668e+01
 * time: 1.213109016418457
    46     4.666098e+02     1.292659e+01
 * time: 1.2166039943695068
    47     4.665692e+02     1.205145e+01
 * time: 1.220120906829834
    48     4.664639e+02     1.210731e+01
 * time: 1.223634958267212
    49     4.662037e+02     1.193827e+01
 * time: 1.2271909713745117
    50     4.656294e+02     1.246811e+01
 * time: 1.2307629585266113
    51     4.646738e+02     1.521581e+01
 * time: 1.234287977218628
    52     4.635500e+02     1.353950e+01
 * time: 1.237696886062622
    53     4.625753e+02     6.648229e+00
 * time: 1.2410330772399902
    54     4.622478e+02     3.557188e+00
 * time: 1.2443010807037354
    55     4.622081e+02     3.487277e+00
 * time: 1.247499942779541
    56     4.622033e+02     3.503679e+00
 * time: 1.250627040863037
    57     4.622028e+02     3.522469e+00
 * time: 1.2539091110229492
    58     4.622027e+02     3.528199e+00
 * time: 1.257232904434204
    59     4.622027e+02     3.533686e+00
 * time: 1.2606050968170166
    60     4.622025e+02     3.543319e+00
 * time: 1.264112949371338
    61     4.622021e+02     3.558678e+00
 * time: 1.2678909301757812
    62     4.622011e+02     3.584789e+00
 * time: 1.2715730667114258
    63     4.621983e+02     3.629313e+00
 * time: 1.2753760814666748
    64     4.621910e+02     3.708490e+00
 * time: 1.2795109748840332
    65     4.621716e+02     3.857090e+00
 * time: 1.283452033996582
    66     4.621190e+02     5.146880e+00
 * time: 1.28751802444458
    67     4.619639e+02     8.785594e+00
 * time: 1.2915899753570557
    68     4.611900e+02     1.480755e+01
 * time: 1.2956769466400146
    69     4.590238e+02     2.147025e+01
 * time: 1.2997329235076904
    70     4.586338e+02     1.072817e+02
 * time: 1.3039319515228271
    71     4.574219e+02     3.842059e+01
 * time: 1.3080439567565918
    72     4.558833e+02     2.153200e+01
 * time: 1.3120388984680176
    73     4.541875e+02     3.129058e+01
 * time: 1.3160150051116943
    74     4.519886e+02     2.114260e+01
 * time: 1.3205509185791016
    75     4.502880e+02     9.897421e+00
 * time: 1.3248710632324219
    76     4.501239e+02     5.986907e+00
 * time: 1.3291599750518799
    77     4.499380e+02     3.686877e+00
 * time: 1.3330209255218506
    78     4.498291e+02     1.392883e+01
 * time: 1.3370709419250488
    79     4.496998e+02     1.236634e+00
 * time: 1.341120958328247
    80     4.496737e+02     5.760379e-01
 * time: 1.3451499938964844
    81     4.496709e+02     5.258588e-01
 * time: 1.3491981029510498
    82     4.496708e+02     5.255864e-01
 * time: 1.3531899452209473
    83     4.496708e+02     5.255561e-01
 * time: 1.357234001159668
    84     4.496708e+02     5.255238e-01
 * time: 1.361232042312622
    85     4.496708e+02     5.254038e-01
 * time: 1.3652620315551758
    86     4.496708e+02     5.252870e-01
 * time: 1.3692851066589355
    87     4.496707e+02     5.252319e-01
 * time: 1.3733150959014893
    88     4.496707e+02     5.252208e-01
 * time: 1.3772668838500977
    89     4.496707e+02     5.252152e-01
 * time: 1.38108491897583
    90     4.496707e+02     5.252180e-01
 * time: 1.3848450183868408
    91     4.496707e+02     5.252473e-01
 * time: 1.3885478973388672
    92     4.496707e+02     5.253567e-01
 * time: 1.3921840190887451
    93     4.496706e+02     5.256970e-01
 * time: 1.3958139419555664
    94     4.496703e+02     5.266758e-01
 * time: 1.4582209587097168
    95     4.496695e+02     7.151059e-01
 * time: 1.461850881576538
    96     4.496675e+02     1.184913e+00
 * time: 1.4653129577636719
    97     4.496620e+02     2.005293e+00
 * time: 1.4687919616699219
    98     4.496459e+02     3.596953e+00
 * time: 1.4722979068756104
    99     4.495867e+02     7.748827e+00
 * time: 1.4758100509643555
   100     4.494486e+02     1.448289e+01
 * time: 1.4798851013183594
   101     4.492033e+02     2.150839e+01
 * time: 1.483901023864746
   102     4.481338e+02     3.873096e+01
 * time: 1.4878599643707275
   103     4.478699e+02     4.381183e+01
 * time: 1.492448091506958
   104     4.471856e+02     6.224688e+01
 * time: 1.496467113494873
   105     4.450564e+02     4.951027e+01
 * time: 1.4998950958251953
   106     4.343470e+02     4.636061e+01
 * time: 1.5033080577850342
   107     4.267316e+02     3.426986e+01
 * time: 1.5072839260101318
   108     4.175694e+02     9.150175e+01
 * time: 1.5107579231262207
   109     4.157174e+02     3.480801e+01
 * time: 1.5147700309753418
   110     4.155775e+02     2.472778e+01
 * time: 1.5182960033416748
   111     4.143712e+02     1.336454e+01
 * time: 1.5217509269714355
   112     4.139564e+02     1.078090e+01
 * time: 1.52512788772583
   113     4.136943e+02     1.025862e+01
 * time: 1.528378963470459
   114     4.133535e+02     1.035907e+01
 * time: 1.5315001010894775
   115     4.129611e+02     1.043919e+01
 * time: 1.5346169471740723
   116     4.125629e+02     6.594579e+00
 * time: 1.5377230644226074
   117     4.122546e+02     2.487983e+00
 * time: 1.5407118797302246
   118     4.121191e+02     3.188012e+00
 * time: 1.5437569618225098
   119     4.120857e+02     3.683605e+00
 * time: 1.5468430519104004
   120     4.119667e+02     3.824309e+00
 * time: 1.5499060153961182
   121     4.117773e+02     6.314328e+00
 * time: 1.5529799461364746
   122     4.115062e+02     5.379015e+00
 * time: 1.556152105331421
   123     4.103579e+02     7.425434e+00
 * time: 1.560072898864746
   124     4.087289e+02     1.410270e+01
 * time: 1.5640521049499512
   125     4.073973e+02     7.419002e+01
 * time: 1.5681769847869873
   126     4.066611e+02     8.314123e+01
 * time: 1.5723309516906738
   127     4.039550e+02     4.224650e+01
 * time: 1.576509952545166
   128     4.009380e+02     6.289836e+02
 * time: 1.5807769298553467
   129     3.873725e+02     1.366039e+02
 * time: 1.5846309661865234
   130     3.850714e+02     2.493050e+03
 * time: 1.5891010761260986
   131     3.782113e+02     5.514808e+02
 * time: 1.592900037765503
   132     3.663733e+02     2.708336e+03
 * time: 1.596893072128296
   133     3.650225e+02     1.533270e+03
 * time: 1.6009070873260498
   134     3.577179e+02     1.856345e+03
 * time: 1.6056311130523682
   135     3.495763e+02     3.463897e+03
 * time: 1.610403060913086
   136     3.419370e+02     7.057854e+03
 * time: 1.6191949844360352
   137     3.407800e+02     6.820066e+03
 * time: 1.8277599811553955
   138     3.387080e+02     1.197776e+04
 * time: 1.8320820331573486
   139     3.320277e+02     2.319497e+05
 * time: 1.8366239070892334
   140     3.059164e+02     5.391851e+04
 * time: 1.8406999111175537
   141     2.991093e+02     6.878031e+04
 * time: 1.8455729484558105
   142     2.953936e+02     9.185108e+04
 * time: 1.8511769771575928
   143     2.932414e+02     2.622618e+05
 * time: 1.8564109802246094
   144     2.888083e+02     3.572167e+05
 * time: 1.861569881439209
   145     2.836077e+02     3.685959e+05
 * time: 1.8668060302734375
   146     2.780634e+02     5.513396e+05
 * time: 1.8723459243774414
   147     2.673610e+02     1.211363e+05
 * time: 1.8776888847351074
   148     2.619248e+02     3.471854e+06
 * time: 1.8838419914245605
   149     2.578706e+02     4.094938e+06
 * time: 1.889328956604004
   150     2.491806e+02     1.613891e+06
 * time: 1.8948719501495361
   151     2.420469e+02     2.814961e+06
 * time: 1.9004158973693848
   152     2.330126e+02     1.276101e+07
 * time: 1.9059638977050781
   153     2.266521e+02     1.767899e+07
 * time: 1.911513090133667
   154     2.159157e+02     1.136918e+07
 * time: 1.9172120094299316
   155     2.077832e+02     2.529660e+07
 * time: 1.922914981842041
   156     1.975289e+02     4.174341e+07
 * time: 1.928692102432251
   157     1.860481e+02     3.810927e+07
 * time: 1.9344770908355713
   158     1.779841e+02     4.879259e+07
 * time: 1.9402339458465576
   159     1.628585e+02     1.895301e+06
 * time: 1.9460809230804443
   160     1.533996e+02     8.832463e+07
 * time: 1.9525830745697021
   161     1.457890e+02     3.666339e+08
 * time: 1.9583048820495605
   162     1.378113e+02     4.345112e+08
 * time: 1.9638359546661377
   163     1.293752e+02     2.686581e+08
 * time: 1.9693210124969482
   164     1.185882e+02     7.011174e+08
 * time: 1.9747869968414307
   165     1.092327e+02     1.053691e+09
 * time: 1.9806289672851562
   166     9.628190e+01     8.925871e+08
 * time: 1.9865880012512207
   167     8.650073e+01     1.463066e+09
 * time: 1.9928040504455566
   168     6.461456e+01     1.542218e+09
 * time: 1.9991240501403809
   169     5.110825e+01     2.702010e+09
 * time: 2.005484104156494
   170     3.190847e+00     6.463555e+09
 * time: 2.012165069580078
   171    -2.866391e+01     9.596926e+09
 * time: 2.0189530849456787
   172    -5.508309e+01     3.429359e+10
 * time: 2.0272860527038574
   173    -5.791486e+01     1.341726e+11
 * time: 2.0361011028289795
   174    -7.457934e+01     3.870780e+11
 * time: 2.0429489612579346
   175    -8.802561e+01     4.146441e+11
 * time: 2.0497610569000244
   176    -9.882043e+01     2.380592e+11
 * time: 2.056593894958496
   177    -1.153795e+02     1.734800e+12
 * time: 2.063647985458374
   178    -1.238138e+02     1.841969e+12
 * time: 2.070713996887207
   179    -1.632583e+02     3.160994e+11
 * time: 2.0769948959350586
   180    -1.706494e+02     3.160994e+11
 * time: 2.1000730991363525
   181    -2.167586e+02     3.160994e+11
 * time: 2.115726947784424
   182    -2.167586e+02     3.160994e+11
 * time: 2.1375889778137207
   183    -2.167586e+02     3.160994e+11
 * time: 2.160944938659668
   184    -2.167586e+02     3.160994e+11
 * time: 2.1853229999542236
   185    -2.167586e+02     3.160994e+11
 * time: 2.2105350494384766
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:              NaivePooled
Likelihood Optimizer:                         BFGS
Dynamical system type:          Matrix exponential

Log-likelihood value:                    216.75859
Number of subjects:                             31
Number of parameters:         Fixed      Optimized
                                  1             10
Observation records:         Active        Missing
    conc:                       239             47
    Total:                      239             47

-------------------------
             Estimate
-------------------------
pop_CL        0.12098
pop_Vc        8.2019
pop_Q         0.013114
pop_Vp        4.6242e12
pop_tabs      2.3621e-12
pop_lag       1.0
pk_Ω₁,₁       0.01
pk_Ω₂,₂       0.01
pk_Ω₃,₃       0.01
σ_prop        0.24211
σ_add         3.0263e-96
-------------------------

As you can see, the constantcoef argument is used to fix the inter-individual variability parameters to zero. The parameter estimates from the naive pooled fit seem out of place for some parameters such as pop_Vp and pop_tabs. Perhaps the initial parameter values are not close to the true values and some adjustments are needed. In Pumas, you can rerun the fit with new initial parameter values very easily. Suppose, we want to run the naive pooled fit again with 4 new initial parameter values that are within 10% of the original values.

We can use the power of the Julia language to generate 4 new parameter values that are within 10% of the original values.

# Generate parameter sets by scaling original values by ±10%
param_sets = [
    NamedTuple(k => v * scale for (k, v) in pairs(param_vals)) for
    scale in [0.9, 0.95, 1.05, 1.1]
]
4-element Vector{@NamedTuple{pop_CL::Float64, pop_Vc::Float64, pop_Q::Float64, pop_Vp::Float64, pop_tabs::Float64, pop_lag::Float64, pk_Ω::Diagonal{Float64, Vector{Float64}}, σ_prop::Float64, σ_add::Float64}}:
 (pop_CL = 0.12060000000000001, pop_Vc = 7.2989999999999995, pop_Q = 0.45, pop_Vp = 18.0, pop_tabs = 0.4707, pop_lag = 0.09000000000000001, pk_Ω = [0.009000000000000001 0.0 0.0; 0.0 0.009000000000000001 0.0; 0.0 0.0 0.009000000000000001], σ_prop = 0.006768, σ_add = 0.05949000000000001)
 (pop_CL = 0.1273, pop_Vc = 7.7044999999999995, pop_Q = 0.475, pop_Vp = 19.0, pop_tabs = 0.49685, pop_lag = 0.095, pk_Ω = [0.0095 0.0 0.0; 0.0 0.0095 0.0; 0.0 0.0 0.0095], σ_prop = 0.007143999999999999, σ_add = 0.062795)
 (pop_CL = 0.14070000000000002, pop_Vc = 8.5155, pop_Q = 0.525, pop_Vp = 21.0, pop_tabs = 0.54915, pop_lag = 0.10500000000000001, pk_Ω = [0.0105 0.0 0.0; 0.0 0.0105 0.0; 0.0 0.0 0.0105], σ_prop = 0.007896, σ_add = 0.06940500000000001)
 (pop_CL = 0.14740000000000003, pop_Vc = 8.921, pop_Q = 0.55, pop_Vp = 22.0, pop_tabs = 0.5753, pop_lag = 0.11000000000000001, pk_Ω = [0.011000000000000001 0.0 0.0; 0.0 0.011000000000000001 0.0; 0.0 0.0 0.011000000000000001], σ_prop = 0.008272, σ_add = 0.07271000000000001)

Now, we can run the fit again with the new parameter values.

new_fits = [
    fit(warfarin_pk_model, pop, param_set, NaivePooled(), constantcoef = (:pk_Ω,)) for
    param_set in param_sets
];

You can compare the results of the new fits with perturbed initial parameter values.

compare_estimates(;
    p1 = new_fits[1],
    p2 = new_fits[2],
    p3 = new_fits[3],
    p4 = new_fits[4],
    original = naive_fit,
)
11×6 DataFrame
Row parameter p1 p2 p3 p4 original
String Float64? Float64? Float64? Float64? Float64?
1 pop_CL 0.0436272 0.000370378 0.0470208 0.127227 0.12098
2 pop_Vc 7.95126 14.3362 9.34242 8.34565 8.20193
3 pop_Q 0.0812342 0.179464 0.0867748 4.31655e-5 0.0131141
4 pop_Vp 181.671 1.92817e16 148.995 1.75913 4.62419e12
5 pop_tabs 1.09629 2.35217e-12 0.0381164 0.37191 2.36208e-12
6 pop_lag 0.5 0.5 0.5 0.764776 1.0
7 pk_Ω₁,₁ 0.009 0.0095 0.0105 0.011 0.01
8 pk_Ω₂,₂ 0.009 0.0095 0.0105 0.011 0.01
9 pk_Ω₃,₃ 0.009 0.0095 0.0105 0.011 0.01
10 σ_prop 0.310588 0.463829 0.256087 0.228482 0.242111
11 σ_add 1.48972e-17 1.39471e-161 6.52968e-11 0.380485 3.0263e-96

The compare_estimates function is a convenience function that is part of the PumasUtilities package. It is used to compare the estimates of the new fits with the original fit. The p4 results seem to be reasonable amongst all the fits and these results will be taken forward for the FOCE() fit.

2.6 Estimating Parameters via Maximum Likelihood (FOCE)

After obtaining initial estimates, the next step is fitting the model to data using a maximum likelihood approach. One of the most commonly used methods in Pumas is the First-Order Conditional Estimation (FOCE) method. This is invoked by specifying FOCE() as the estimation algorithm:

foce_fit = fit(warfarin_pk_model, pop, coef(new_fits[4]), FOCE())

The coef(new_fits[4]) function extracts the parameter estimates from the new_fits[4] fit and passes them to the fit function.

Note

By default, the fit function will use all the cores of your machine to parallelize the optimization. You can control this using the ensemblealg argument. For more details, see the fit options section in the Pumas Documentation.

2.7 Fixing Parameters with constantcoef

Sometimes, it is desirable to hold certain parameters constant during estimation. For instance, suppose one wants to fix pop_Q at a particular value:

foce_fit_fixedQ = fit(
    warfarin_pk_model,
    pop,
    (; param_vals..., pop_Q = 0.5),
    FOCE();
    constantcoef = (:pop_Q,),
);   # This fixes the pop_Q parameter
foce_fit_fixedQ
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
Likelihood Optimizer:                         BFGS
Dynamical system type:          Matrix exponential

Log-likelihood value:                    -321.9637
Number of subjects:                             31
Number of parameters:         Fixed      Optimized
                                  1             10
Observation records:         Active        Missing
    conc:                       239             47
    Total:                      239             47

------------------------
             Estimate
------------------------
pop_CL        0.13128
pop_Vc        8.0921
pop_Q         0.5
pop_Vp        0.0040524
pop_tabs      0.44966
pop_lag       0.92947
pk_Ω₁,₁       0.049875
pk_Ω₂,₂       0.018981
pk_Ω₃,₃       1.7058
σ_prop        0.090433
σ_add         0.33064
------------------------

This approach is useful for sensitivity analysis, exploring different values, or simplifying the model.

2.8 Changing Tolerance During Fitting

The fit function accepts additional keyword arguments to tweak the fitting process, such as tolerances for convergence:

foce_fit_fixedQ = fit(
    warfarin_pk_model,
    pop,
    (; param_vals..., pop_Q = 0.5),
    FOCE();
    constantcoef = (pop_Q = 0.5,),
    optim_options = (; x_reltol = 1e-6, x_abstol = 1e-8, iterations = 200),
);

These arguments (x_reltol, x_abstol, iterations) help fine-tune the solver’s stopping criteria. For more details, see the Optimization Options during Fitting.

2.9 Interpreting the Printed Results

Once the model fit is complete, Pumas will print a summary of the fit to the console. This typically includes:

  1. Objective Function Value (OFV): The minimized negative log-likelihood (or similar objective). Lower is better.
  2. Convergence Status: An indication of whether the solver found a local optimum or if it ran into issues (e.g., maximum iterations reached).
  3. Parameter Estimates: A table listing final estimates of fixed effects, random effects (variances or standard deviations), and residual variability components.
  4. Standard Errors (if computed): By default, these might not be computed until using infer (discussed in the next section).

The printed output might look like:

julia> foce_fit_fixedQ
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
Likelihood Optimizer:                         BFGS
Dynamical system type:          Matrix exponential

Log-likelihood value:                    -321.9637
Number of subjects:                             31
Number of parameters:         Fixed      Optimized
                                  1             10
Observation records:         Active        Missing
    conc:                       239             47
    Total:                      239             47

------------------------
             Estimate
------------------------
pop_CL        0.13128
pop_Vc        8.0919
pop_Q         0.5
pop_Vp        0.0042035
pop_tabs      0.44966
pop_lag       0.92947
pk_Ω₁,₁       0.049875
pk_Ω₂,₂       0.018981
pk_Ω₃,₃       1.7058
σ_prop        0.090433
σ_add         0.33064
------------------------

We can break down the output of the fit function into the following components:

2.9.1 FittedPumasModel

FittedPumasModel

The FittedPumasModel object contains the following fields:

  • model: The original model definition
foce_fit_fixedQ.model
PumasModel
  Parameters: pop_CL, pop_Vc, pop_Q, pop_Vp, pop_tabs, pop_lag, pk_Ω, σ_prop, σ_add
  Random effects: pk_η
  Covariates: FSZCL, FSZV
  Dynamical system variables: Depot, Central, Peripheral
  Dynamical system type: Matrix exponential
  Derived: conc
  Observed: conc
  • data: The population data
foce_fit_fixedQ.data
Population
  Subjects: 31
  Covariates: SEX, WEIGHT, FSZV, FSZCL
  Observations: conc
  • optim: The optimization results
foce_fit_fixedQ.optim
 * Status: success

 * Candidate solution
    Final objective value:     3.219637e+02

 * Found with
    Algorithm:     BFGS

 * Convergence measures
    |x - x'|               = 0.00e+00 ≤ 1.0e-08
    |x - x'|/|x'|          = 0.00e+00 ≤ 1.0e-06
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 1.49e-03 ≰ 1.0e-03

 * Work counters
    Seconds run:   28  (vs limit Inf)
    Iterations:    172
    f(x) calls:    206
    ∇f(x) calls:   173
  • approx: The likelihood approximation method
foce_fit_fixedQ.approx
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_warnings = true, show_every = 1, time_limit = NaN, )
)
  • kwargs: Additional keyword arguments
foce_fit_fixedQ.kwargs
  • fixedparamset: The fixed parameter set
foce_fit_fixedQ.fixedparamset
  • optim_state: The optimization state
foce_fit_fixedQ.optim_state

2.9.2 Minimization Status

Successful minimization:                      true

The Successful minimization field indicates whether the minimization was successful. It depends on the convergence of the optimizer and rules set in the optim_options argument, for example x_reltol and x_abstol.

2.9.3 Likelihood Approximation

  Likelihood approximation:                     FOCE

The Likelihood approximation field indicates the likelihood approximation method used.

2.9.4 Likelihood Optimizer

  Likelihood Optimizer:                         BFGS

The Likelihood Optimizer field indicates the optimizer used which by default is BFGS from the Optim package. A user can change the optimizer by setting the optim_alg argument in the fit function.

2.9.5 Dynamical system type

  Dynamical system type:          Matrix exponential

The Dynamical system type field indicates the type of dynamical system used. In the example here, even though the model is a differential equation model, the Matrix exponential indicates that the model is solved using a matrix exponential solver as the system has been deemed linear. The user can check the linear nature of the system as follows: This example model is: {julia} foce_fit.model.prob. If the user does not want the automatic check for linearity, they can turn it off in the @options block of the model definition.

@model begin
    @param begin
        ...
    end
    @random begin
        ...
    end
    @dynamics begin
        ...
    end
    @pre begin
        ...
    end
    @options begin
        checklinear = false
    end
end

The checklinear option is a boolean that determines whether the solver should check if the system defined in the @dynamics block is linear. If the system is linear, setting this option to true (default) enables calculating the solution through matrix exponentials. If it is set to false or time (t) appears in the @pre block, this optimization is disabled. This option can be useful when the matrix exponential solver is not superior to general numerical integrators or for debugging purposes.

2.9.6 Log-likelihood value

  Log-likelihood value:                    -321.9637

The Log-likelihood value field indicates the value of the log-likelihood function at the maximum likelihood estimate.

2.9.7 Number of subjects

  Number of subjects:                             31

The Number of subjects field indicates the number of subjects from the population data used in the fit.

2.9.8 Number of parameters

  Number of parameters:         Fixed      Optimized
                                  1             10

The Number of parameters field indicates the number of fixed and optimized parameters. The Fixed column indicates the number of parameters that were fixed using the constantcoef argument during the fit and the Optimized column indicates the number of parameters that were optimized.

2.9.9 Observation records

  Observation records:         Active        Missing
    conc:                       239             47
    Total:                      239             47

The Observation records field indicates the number of active and missing observations in the population data. The Active column indicates the number of observations that were used in the fit and the Missing column indicates the number of observations that were missing and not used in the fit. missing in this case is due to the missingness of the conc observations, which is the collection of all the conc observations from all the subjects.

2.9.10 Parameter estimates

------------------------
             Estimate
------------------------
pop_CL        0.13128
pop_Vc        8.0919
pop_Q         0.5
pop_Vp        0.0042035
pop_tabs      0.44966
pop_lag       0.92947
pk_Ω₁,₁       0.049875
pk_Ω₂,₂       0.018981
pk_Ω₃,₃       1.7058
σ_prop        0.090433
σ_add         0.33064
------------------------

The Parameter estimates field indicates the final parameter estimates. The Estimate column indicates the estimated value of the parameter.

2.10 Computing Parameter Precision with infer

The infer function in Pumas estimates the uncertainty (precision) of parameter estimates. Depending on the chosen method, infer can provide standard errors, confidence intervals, and correlation matrices.

The signature for infer often looks like:

infer(
    fpm::FittedPumasModel;
    level = 0.95,
    rethrow_error::Bool = false,
    sandwich_estimator::Bool = true,
)

where:

  • fpm::FittedPumasModel: The result of fit (e.g., foce_fit).
  • level: The confidence interval level (e.g., 0.95). The confidence intervals are calculated as the (1-level)/2 and (1+level)/2 quantiles of the estimated parameters
  • rethrow_error: If rethrow_error is false (the default value), no error will be thrown if the variance-covariance matrix estimator fails. If it is true, an error will be thrown if the estimator fails.
  • sandwich_estimator: Whether to use the sandwich estimator. If set to true (the default value), the sandwich estimator will be used. If set to false, the standard error will be calculated using the inverse of the Hessian matrix.

An example usage:

inference_results = infer(foce_fit_fixedQ; level = 0.95)
[ 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:          Matrix exponential

Log-likelihood value:                    -321.9637
Number of subjects:                             31
Number of parameters:         Fixed      Optimized
                                  1             10
Observation records:         Active        Missing
    conc:                       239             47
    Total:                      239             47

-------------------------------------------------------------------------
             Estimate            SE                      95.0% C.I.
-------------------------------------------------------------------------
pop_CL        0.13128          0.005417         [ 0.12066  ; 0.14189 ]
pop_Vc        8.0921           0.23906          [ 7.6235   ; 8.5606  ]
pop_Q         0.5              NaN              [  NaN     ;  NaN      ]
pop_Vp        0.0040524        0.0098746        [-0.015301 ; 0.023406]
pop_tabs      0.44966          0.22631          [ 0.0060942; 0.89322 ]
pop_lag       0.92947          0.052185         [ 0.82719  ; 1.0318  ]
pk_Ω₁,₁       0.049875         0.01679          [ 0.016967 ; 0.082782]
pk_Ω₂,₂       0.018981         0.0057796        [ 0.0076528; 0.030308]
pk_Ω₃,₃       1.7058           1.1673           [-0.58197  ; 3.9936  ]
σ_prop        0.090433         0.014936         [ 0.061159 ; 0.11971 ]
σ_add         0.33064          0.086237         [ 0.16162  ; 0.49967 ]
-------------------------------------------------------------------------

We can use the sandwich_estimator argument to get a more robust estimate of the standard errors.

inference_results = infer(foce_fit_fixedQ; level = 0.95, sandwich_estimator = false)
[ Info: Calculating: variance-covariance matrix.
[ Info: Done.
Asymptotic inference results using negative inverse Hessian

Successful minimization:                      true

Likelihood approximation:                     FOCE
Likelihood Optimizer:                         BFGS
Dynamical system type:          Matrix exponential

Log-likelihood value:                    -321.9637
Number of subjects:                             31
Number of parameters:         Fixed      Optimized
                                  1             10
Observation records:         Active        Missing
    conc:                       239             47
    Total:                      239             47

------------------------------------------------------------------------
             Estimate            SE                     95.0% C.I.
------------------------------------------------------------------------
pop_CL        0.13128          0.0054515       [ 0.12059  ; 0.14196 ]
pop_Vc        8.0921           0.47585         [ 7.1594   ; 9.0247  ]
pop_Q         0.5              NaN             [  NaN     ;  NaN      ]
pop_Vp        0.0040524        0.41066         [-0.80083  ; 0.80894 ]
pop_tabs      0.44966          0.21485         [ 0.028563 ; 0.87075 ]
pop_lag       0.92947          0.032557        [ 0.86566  ; 0.99328 ]
pk_Ω₁,₁       0.049875         0.013572        [ 0.023273 ; 0.076476]
pk_Ω₂,₂       0.018981         0.0067243       [ 0.0058011; 0.03216 ]
pk_Ω₃,₃       1.7058           0.97858         [-0.21215  ; 3.6238  ]
σ_prop        0.090433         0.0086792       [ 0.073422 ; 0.10744 ]
σ_add         0.33064          0.047455        [ 0.23763  ; 0.42366 ]
------------------------------------------------------------------------

This result above indicates that both with the sandwich estimator and the inverse of the Hessian matrix on the fixed parameter, the inference worked. We will pursue other methods to obtain parameter precision in later tutorials, such as bootstrap and SIR.

3 Concluding Remarks

This tutorial showcased a typical Pumas workflow for PK model fitting using a Warfarin dataset:

  1. Model Definition and Data Preparation.
  2. Checking Compatibility via loglikelihood and findinfluential.
  3. Initial Parameter Estimation using Naive Pooled.
  4. Fitting with FOCE, demonstrating how to fix parameters or change solver tolerances.
  5. Interpreting Fit Results, exploring the components of a fittedpumasmodel.
  6. Computing Precision of estimates with infer using different methods.

Readers are encouraged to refine their understanding by:

  • Performing thorough Exploratory Data Analysis prior to modeling.
  • Exploring Bootstrap or SIR methods for deeper uncertainty quantification (covered in later tutorials).
  • Validating the final model through visual diagnostics (e.g., VPCs, GOF plots).

More advanced topics can be found throughout the Pumas Documentation and Pumas Tutorials.

Tip

Real-world workflows are iterative. Analysts often revisit model assumptions, re-check data, and explore alternative parameterizations before finalizing any one model as “best.”