using Pumas
using PharmaDatasets
using DataFramesMeta
using CairoMakie
using AlgebraOfGraphics
using PumasUtilities
using Markdown
Interpreting Results from Pumas Fits
1 Introduction
After following a typical workflow for fitting a Pumas model it is important to assess and interpret the results. This includes steps such as
- Looking at log-likelihood related quantities such as
loglikelihood,aicandbicto compare the fit to other existing fits. - Assessing convergence of the maximum likelihood optimization.
- Calculating the condition number of the correlation matrix implied by the estimated variance-covariance matrix.
- Calculating shrinkage for random effects and the dependent variable.
The following sections will walk through these steps using a one-compartment PK model for Warfarin, focusing on the PK aspects only.
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_tabs, pop_lag - Inter-individual variability (IIV) components:
pk_Ω - Residual error model parameters:
σ_prop, σ_add - Covariates for scaling:
FSZCLandFSZV - Differential equations describing the PK behavior in the compartments.
warfarin_pk_model = @model begin
@metadata begin
desc = "Warfarin 1-compartment PK model (PD removed)"
timeu = "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)
"""
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])
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
end
@derived begin
"""
Concentration (ng/mL)
"""
conc ~ @. CombinedNormal(cp, σ_add, σ_prop)
end
endPumasModel
Parameters: pop_CL, pop_Vc, pop_tabs, pop_lag, pk_Ω, σ_prop, σ_add
Random effects: pk_η
Covariates: FSZCL, FSZV
Dynamical system variables: Depot, Central
Dynamical system type: Matrix exponential
Derived: conc
Observed: conc2.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("pumas/warfarin_nonmem")
# 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
# Transform the data in a single chain of operations
warfarin_data_wide = @chain warfarin_data begin
@rtransform begin
# Scaling factors
:FSZV = :wtbl / 70 # volume scaling
:FSZCL = (:wtbl / 70)^0.75 # clearance scaling (allometric)
# Column name for the dv
:dvname = "dv$(:dvid)"
# Dosing indicator columns
:cmt = ismissing(:amt) ? missing : 1
:evid = iszero(:amt) ? 0 : 1
end
unstack(Not([:dvid, :dvname, :dv]), :dvname, :dv)
rename!(:dv1 => :conc, :dv2 => :pca)
select!(Not([:dv0]))
# Concentrations of exactly zero are below the limit of quantification
@rtransform @passmissing :conc = iszero(:conc) ? missing : :conc
end
first(warfarin_data_wide, 10)10×12 DataFrame
| Row | id | time | evid | amt | wtbl | age | sex | FSZV | FSZCL | cmt | conc | pca |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Int64 | Float64 | Int64 | Float64 | Float64 | Int64 | String1 | Float64 | Float64 | Int64 | Float64? | Float64? | |
| 1 | 1 | 0.0 | 1 | 100.0 | 66.7 | 50 | M | 0.952857 | 0.96443 | 1 | missing | missing |
| 2 | 1 | 0.5 | 0 | 0.0 | 66.7 | 50 | M | 0.952857 | 0.96443 | 1 | missing | missing |
| 3 | 1 | 1.0 | 0 | 0.0 | 66.7 | 50 | M | 0.952857 | 0.96443 | 1 | 1.9 | missing |
| 4 | 1 | 2.0 | 0 | 0.0 | 66.7 | 50 | M | 0.952857 | 0.96443 | 1 | 3.3 | missing |
| 5 | 1 | 3.0 | 0 | 0.0 | 66.7 | 50 | M | 0.952857 | 0.96443 | 1 | 6.6 | missing |
| 6 | 1 | 6.0 | 0 | 0.0 | 66.7 | 50 | M | 0.952857 | 0.96443 | 1 | 9.1 | missing |
| 7 | 1 | 9.0 | 0 | 0.0 | 66.7 | 50 | M | 0.952857 | 0.96443 | 1 | 10.8 | missing |
| 8 | 1 | 12.0 | 0 | 0.0 | 66.7 | 50 | M | 0.952857 | 0.96443 | 1 | 8.6 | missing |
| 9 | 1 | 24.0 | 0 | 0.0 | 66.7 | 50 | M | 0.952857 | 0.96443 | 1 | 5.6 | 44.0 |
| 10 | 1 | 36.0 | 0 | 0.0 | 66.7 | 50 | M | 0.952857 | 0.96443 | 1 | 4.0 | 27.0 |
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 = :amt,
cmt = :cmt,
evid = :evid,
covariates = [:sex, :wtbl, :FSZV, :FSZCL],
observations = [:conc],
)Population
Subjects: 32
Covariates: sex, wtbl, FSZV, FSZCL
Observations: conc# A named tuple of parameter values
param_vals = (
pop_CL = 0.134,
pop_Vc = 8.11,
pop_tabs = 0.523,
pop_lag = 0.1,
pk_Ω = Diagonal([0.01, 0.01, 0.01]),
σ_prop = 0.00752,
σ_add = 0.0661,
)
foce_fit = fit(warfarin_pk_model, pop, param_vals, 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 1.030400e+04 1.181292e+04 * time: 0.07335996627807617 1 2.262448e+03 2.457021e+03 * time: 9.045576810836792 2 1.606175e+03 1.689206e+03 * time: 9.07464599609375 3 8.249479e+02 6.944903e+02 * time: 9.102895975112915 4 5.830633e+02 3.309739e+02 * time: 9.130812883377075 5 4.835343e+02 1.258919e+02 * time: 9.159070014953613 6 4.597452e+02 3.446181e+01 * time: 9.186965942382812 7 4.564004e+02 2.157414e+01 * time: 9.214821815490723 8 4.561766e+02 1.841633e+01 * time: 9.241413831710815 9 4.560933e+02 1.718302e+01 * time: 9.27746295928955 10 4.557723e+02 1.676719e+01 * time: 9.322022914886475 11 4.552137e+02 2.081976e+01 * time: 9.368091821670532 12 4.542913e+02 2.597069e+01 * time: 9.400276899337769 13 4.534692e+02 3.527716e+01 * time: 9.426272869110107 14 4.529383e+02 3.402880e+01 * time: 9.453419923782349 15 4.526116e+02 2.490635e+01 * time: 9.479892015457153 16 4.523614e+02 1.340594e+01 * time: 9.505985021591187 17 4.520146e+02 1.295774e+01 * time: 9.532482862472534 18 4.517070e+02 1.823845e+01 * time: 9.560969829559326 19 4.514823e+02 2.416090e+01 * time: 9.589457988739014 20 4.512324e+02 2.583177e+01 * time: 9.617363929748535 21 4.506901e+02 2.632776e+01 * time: 9.645208835601807 22 4.493304e+02 3.098682e+01 * time: 9.673577785491943 23 4.457731e+02 4.060488e+01 * time: 9.702276945114136 24 4.374540e+02 5.129320e+01 * time: 9.732040882110596 25 4.278696e+02 5.289924e+01 * time: 9.763497829437256 26 4.221945e+02 3.807050e+01 * time: 9.794572830200195 27 4.126997e+02 2.293280e+01 * time: 9.829111814498901 28 4.085329e+02 2.226988e+01 * time: 9.86014699935913 29 4.073530e+02 1.516843e+01 * time: 10.080879926681519 30 4.066443e+02 1.009682e+01 * time: 10.109987020492554 31 4.064325e+02 7.665676e+00 * time: 10.137465000152588 32 4.063169e+02 7.945666e+00 * time: 10.165079832077026 33 4.062598e+02 8.057281e+00 * time: 10.191873788833618 34 4.061917e+02 8.096741e+00 * time: 10.218864917755127 35 4.060744e+02 1.001791e+01 * time: 10.245957851409912 36 4.058146e+02 1.313770e+01 * time: 10.273317813873291 37 4.052323e+02 1.489612e+01 * time: 10.301275968551636 38 4.040793e+02 2.371901e+01 * time: 10.329726934432983 39 4.022553e+02 3.961615e+01 * time: 10.35873794555664 40 4.000227e+02 4.701290e+01 * time: 10.389318943023682 41 3.977591e+02 4.815312e+01 * time: 10.41731595993042 42 3.965739e+02 4.653027e+01 * time: 10.44510293006897 43 3.916661e+02 4.293389e+01 * time: 10.473394870758057 44 3.895667e+02 4.702429e+01 * time: 10.500940799713135 45 3.889001e+02 5.508621e+01 * time: 10.530346870422363 46 3.868670e+02 5.359373e+01 * time: 10.559207916259766 47 3.857187e+02 5.483763e+01 * time: 10.58902883529663 48 3.801079e+02 6.169025e+01 * time: 10.618764877319336 49 3.758914e+02 5.242790e+01 * time: 10.646840810775757 50 3.704875e+02 3.681725e+01 * time: 10.677296876907349 51 3.647525e+02 6.405499e+01 * time: 10.707692861557007 52 3.593661e+02 5.044040e+01 * time: 10.739197969436646 53 3.525545e+02 2.168782e+01 * time: 10.771612882614136 54 3.519477e+02 9.794489e+00 * time: 10.802923917770386 55 3.518024e+02 1.071209e+01 * time: 10.894738912582397 56 3.517873e+02 1.052820e+01 * time: 10.922974824905396 57 3.516873e+02 8.722340e+00 * time: 10.95199990272522 58 3.515999e+02 7.655533e+00 * time: 10.981554985046387 59 3.513182e+02 3.886802e+00 * time: 11.008091926574707 60 3.511058e+02 4.956859e+00 * time: 11.033994913101196 61 3.510219e+02 3.013130e+00 * time: 11.060414791107178 62 3.510122e+02 2.597003e+00 * time: 11.0867338180542 63 3.510101e+02 2.262248e+00 * time: 11.113844871520996 64 3.510069e+02 2.230919e+00 * time: 11.14003586769104 65 3.509975e+02 2.157109e+00 * time: 11.166553974151611 66 3.509747e+02 2.135766e+00 * time: 11.192803859710693 67 3.509195e+02 3.319996e+00 * time: 11.219040870666504 68 3.508069e+02 4.888964e+00 * time: 11.246124982833862 69 3.506399e+02 5.090855e+00 * time: 11.273002862930298 70 3.505135e+02 2.778974e+00 * time: 11.300212860107422 71 3.504733e+02 9.941507e-01 * time: 11.32697582244873 72 3.504659e+02 2.186575e-01 * time: 11.353760957717896 73 3.504655e+02 1.956737e-01 * time: 11.383034944534302 74 3.504655e+02 1.910197e-01 * time: 11.426628828048706 75 3.504654e+02 1.883550e-01 * time: 11.469931840896606 76 3.504654e+02 1.890016e-01 * time: 11.510186910629272 77 3.504653e+02 1.881134e-01 * time: 11.53519582748413 78 3.504651e+02 1.815387e-01 * time: 11.582655906677246 79 3.504646e+02 1.749731e-01 * time: 11.608130931854248 80 3.504636e+02 2.187541e-01 * time: 11.635742902755737 81 3.504626e+02 1.812098e-01 * time: 11.662891864776611 82 3.504620e+02 7.622078e-02 * time: 11.688638925552368 83 3.504619e+02 1.141440e-02 * time: 11.714066982269287 84 3.504619e+02 9.572856e-04 * time: 11.739869832992554
FittedPumasModel
Dynamical system type: Matrix exponential
Number of subjects: 32
Observation records: Active Missing
conc: 247 51
Total: 247 51
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.46191
--------------------
Estimate
--------------------
pop_CL 0.13469
pop_Vc 8.048
pop_tabs 0.5449
pop_lag 0.87702
pk_Ω₁,₁ 0.070708
pk_Ω₂,₂ 0.018349
pk_Ω₃,₃ 0.92441
σ_prop 0.087933
σ_add 0.41375
--------------------2.4 Evaluating the termination of the fit
Fitting non-linear mixed effects model is not always a simple task. The model has to be appropriate to the data, data must be accurately sampled and recorded, no information should be truncated, and so on for the fits to run reliably. Even in the most perfect version of data collection and modeling there can still be numerical challenges in evaluating various expressions, calculating empirical bayes estimates in an inner loop, handling vanishing variances, local optima and more.
In Pumas we use Optim.jl (Mogensen and Riseth 2018) to perform the minimization of any FO, FOCE, or LaplaceI fits with out without MAP for maximum aposterori fitting. The returned FittedPumasModel contains the optimization results structure from Optim.jl in the optim field. This can be useful when assessing convergence.
Convergence in non-linear optimization is not always a simply task to assess:
- We can look at the gradient element sizes and we often look at the infinity norm of the gradient to this end. The infinity norm of a vector picks out the largest absolute value of the elements. By default we require this to be less than
1e-3to declare convergence, but some problems may be more or less numerically stable so that number may be large or small depending on the circumstances.
NoteGradient Elements and Infinity Norm of the Gradient
- Gradient Elements
- The gradient is a vector that tells us which direction to move to increase a function’s value the fastest
- Each element in this vector corresponds to how much the function changes when we adjust one parameter
- “Gradient element sizes” refers to the magnitude (absolute value) of each of these directional changes
To illustrate with a simple example, consider a model with 3 parameters (like CL, V, and Ka), a gradient might look like:
gradient = [0.0001, -0.002, 0.0005]3-element Vector{Float64}:
0.0001
-0.002
0.0005Here:
- 0.0001 tells us how much the function would change if we adjust the first parameter (
CL) - -0.002 for the second parameter (
V) - 0.0005 for the third parameter (
Ka)
- Infinity Norm of the Gradient
- The infinity norm is simply the largest absolute value among all gradient elements
- In our example above:
|0.0001| = 0.0001
|-0.002| = 0.002 # <- This is the largest!
|0.0005| = 0.0005Therefore, infinity norm = 0.002
Why This Matters for Convergence:
- When optimizing a model (finding the best parameters), we want all gradient elements to be close to zero
- A large gradient element means we can still improve our solution by moving in that direction
- A small gradient (all elements close to zero) suggests we’re at or near the optimal solution
- By convention in Pumas, if the infinity norm is less than 0.001 (1e-3), we consider the optimization “converged”
Think of it like hiking to find the highest point:
- The gradient tells you which direction to walk
- The size of each gradient element tells you how steep it is in each direction
- When you’re at the top, it should be flat in all directions (all gradient elements close to zero)
- The infinity norm tells you about the steepest direction you could still go
- If the steepest direction is very gentle (< 0.001), you’re probably at the peak!
- We do not stop according to small changes in the parameters or objective function value by default. This is because it is typically not a direct sign of convergence but rather that progress is stalled. That said, many researches set a relative change tolerance in the parameters in the same way that is done in NONMEM because often you will see many final iterations with next to no improvement. This can be done by setting
optim_options = (x_reltol = 1e-5)infitto match NONMEM for example.
Let us consider a bunch of different cases and how you can find and interpret the information in the optim field. This looks something like the following.
julia> foce_fit.optim
L1 * Status: success
L2
L3 * Candidate solution
L4 Final objective value: 3.199344e+02
L5
L6 * Found with
L7 Algorithm: BFGS
L8
L9 * Convergence measures
L10 |x - x'| = 4.36e-04 ≰ 0.0e+00
L11 |x - x'|/|x'| = 1.10e-04 ≰ 0.0e+00
L12 |f(x) - f(x')| = 3.57e-06 ≰ 0.0e+00
L13 |f(x) - f(x')|/|f(x')| = 1.12e-08 ≰ 0.0e+00
L14 |g(x)| = 5.62e-04 ≤ 1.0e-03
L15
L16 * Work counters
L17 Seconds run: 1 (vs limit Inf)
L18 Iterations: 95
L19 f(x) calls: 101
L20 ∇f(x) calls: 97We will use the line numbers to the left below to refer to the relevant part of the output. While interpreting the output it may be worth remembering that the optimization is cast as the problem of minimizing the negative log of the marginal likelihood.
2.4.1 Case 1: gradient convergence
Since a (local) minimum of a smooth optimization function satisfies that the partial derivatives are zero we will often look at the norm of the gradient to declare convergence. The default value for the tolerance set as the g_tol option in optim_options in Pumas is 1e-3. The final state of the measure is found in L14 in the output above where it is seen that the tolerance was not met. This tolerance is set to match what is often possible to achieve in well-specified models, but smaller is always better. To derive a theoretically well-reasoned value we would need to know the Hessian at the optimum as this is not possible to know before running a fit and very expensive to calculate exactly after a fit. Several factors affect the level that can be reached in practice including the accuracy of the solution of the embedded dynamical system and the curvature of the objective function surface. If the gradient norm is lower than 1e-3 we will in general accept the fit as having converged, but “near convergence” is often accepted as well if there are not other indications of issues with the model.
2.4.2 Case 2: no change in parameters
If x_abstol or x_reltol was not manually set by the user we will only terminate the optimization procedure if no change in parameters could be made at all from one iteration to the next. The final state of these termination tests are found in L10 and L11 respectively. This can indicate a few things but essentially indicates that the procedure that scales the search direction to find the next value of the parameters failed. If this occurs with a relatively small gradient norm we may accept the solution though the termination reason cannot be used to conclude convergence in isolation. The reasons why no step led to an improvement of the objective function value can be many. Some reasons include:
- The objective function is not smooth due to numerical noise from calculations in the objective function including model solutions.
- The inverse Hessian approximation in the BFGS algorithm is far from the true inverse Hessian of the objective function. Restarting the fit may help.
As mentioned above, NONMEM sets x_reltol=1e-5 by default but it is 0 in Pumas.
2.4.3 Case 3: small change in parameters
As mentioned above, Pumas defaults to only terminate due to the size of the change in parameters from one iteration to the next is exactly zero. If the relative change tolerance x_reltol is set to a number such as 1e-5 to match other software such as NONMEM then the user may well accept the solution since they deliberately used this criterion. In general, we advise using this setting with caution as it can terminate with quite large gradient norm.
Again, this information is seen in L10 and L11 but the limit after the inequality will change to the number provided.
2.4.4 Case 4: no change in objective function
This case is related to Case 2. We do not allow for termination due to the absence of improvement in the objective function, so in general this should be seen as a failure. The rest of the comments about reasons for stopping, acceptance based on gradient norms etc applies here as well.
The information about these termination tests can be found in L12 and L13.
2.4.5 Case 5: small change in objective function
Just as Case 4 is related to Case 2, this case is related to Case 3. A very small improvement is often not a valid reason to accept the returned solution as the maximum likelihood estimates. If the user has good reason to actively set one of the options f_abstol or f_reltol, the result may be accepted.
As in Case 4, the information can be found in L12 and L13, but the right-hand side of the inequalities are updated as in Case 3.
2.4.6 Case 6: “maximum number of iterations reached”
The pre-specified limit on the number of iterations was reached. This can be increased in optim_options by setting the iterations keyword. This defaults to 1000 which should generally be plenty so termination according to this criterion is not a sign of convergence. The output may be accepted if the solution seems to cycle between parameter values with a gradient norm close to the prescribed tolerance.
The number iterations can be found in L18.
2.4.7 Case 7: line search failed
This case is related to Case 2, 3, 4, and 5. Given the search direction it was not possible to find a step that sufficiently decreases the objective function value. The advice found in Case 2 applies here as well.
In the output this will be indicates in L1.
2.4.8 Case 8: “user provided callback termination”
This is an advanced use-case where a user-provided callback function returned the value true to indicate that the optimization should terminate. Can only happen if the user supplied a custom callback function as part of optim_options.
In the output this will be indicated in L1.
2.4.9 Case 9: “time limit reached”
The optimization routine allows for a soft time limit that is checked after every iteration. The default allows for infinite run-time, but if the time_limit option in optim_options is set this termination reason may occur.
In the output this will be indicated in L17.
2.4.10 Case 10: “objective evaluation limit reached
The optimization routine allows for a limit to be placed on the total number of objective function evaluations through the f_calls_limit option in optim_options. It is not set by default in Pumas. The interpretation is similar to Case 6.
In the output this will be indicated in L19.
2.4.11 Case 11: “gradient evaluation limit reached
The optimization routine allows for a limit to be placed on the total number of gradient function evaluations through the g_calls_limit option in optim_options. It is not set by default in Pumas. The interpretation is similar to Case 6.
In the output this will be indicated in L20.
2.4.12 Case 12: “hessian evaluation limit reached
The optimization routine allows for a limit to be placed on the total number of hessian function evaluations through the h_calls_limit option in optim_options. It is not set by default in Pumas. The interpretation is similar to Case 6.
In the output this would be placed after L20 if it was present, but it is not because the BFGS algorithm only maintains an approximation to the inverse Hessian it doesn’t actually compute and utilize the real Hessian.
2.4.13 Case 13: “objective function increased”
This can happen if the allow_f_increases option is set to false in optim_options. The default is true. The causes are the same as in Case 2.
This would be indicated in L1 if it failed.
2.5 Optimization didn’t succeed - then what?
If the optimization did not success according to the description above it is always important to think about what to do to improve the situation if possible. If the model is or initial parameters are grossly wrong there may not be much to do about the failed optimization. It may always be a good idea to try to reduce the complexity of the model by removing random effects that cannot be supported by the data or fixing parameters that may be know with high certainty from the relevant literature.
Before changing the model too much it may be worth trying a few different things.
2.5.1 Try different starting parameters
Given the complexity and high number of parameters in many NLME models it is always important to initialize the fit from a set of starting parameters instead of only one parameter NamedTuple. It is quite likely that the optimization can end up in a local minimum that is not the global minimum, so we need to check the robustness of our maximum likelihood estimate candidates. Performing this robustness step is explained more in detail in a different tutorial, but just to show the effect we can try with a different initial fixed effect guess for volume of distribution.
foce_fit2 = fit(
warfarin_pk_model,
pop,
merge(param_vals, (; pop_Vc = 12.0)),
FOCE();
optim_options = (show_trace = false,),
)[ Info: Checking the initial parameter values. [ Info: The initial negative log-likelihood and its gradient are finite. Check passed.
FittedPumasModel
Dynamical system type: Matrix exponential
Number of subjects: 32
Observation records: Active Missing
conc: 247 51
Total: 247 51
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.46191
--------------------
Estimate
--------------------
pop_CL 0.13469
pop_Vc 8.048
pop_tabs 0.5449
pop_lag 0.87702
pk_Ω₁,₁ 0.070709
pk_Ω₂,₂ 0.018349
pk_Ω₃,₃ 0.92439
σ_prop 0.087933
σ_add 0.41375
--------------------2.5.2 Resetting the optimization procedure
Whenever a model fit concludes we will have ended up with a final set of fixed effect estimates (population level parameters) and well as a set of empirical bayes estimates. We also end up with a final inverse Hessian approximation. There are ways to restart the optimization at the final fixed effects parameters while resetting the random effects to the central value of the prior and to reset the optimization procedure.
To restart the optimization procedure and use the final fixed effects and empirical bayes estimates use the following.
foce_fit3 = fit(foce_fit2; reset_optim = true)[ Info: Checking the initial parameter values. [ Info: The initial negative log-likelihood and its gradient are finite. Check passed.
FittedPumasModel
Dynamical system type: Matrix exponential
Number of subjects: 32
Observation records: Active Missing
conc: 247 51
Total: 247 51
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.46191
--------------------
Estimate
--------------------
pop_CL 0.13469
pop_Vc 8.048
pop_tabs 0.5449
pop_lag 0.87702
pk_Ω₁,₁ 0.070709
pk_Ω₂,₂ 0.018349
pk_Ω₃,₃ 0.92439
σ_prop 0.087933
σ_add 0.41375
--------------------2.5.3 Resetting the optimization procedure and the empirical Bayes estimates
To restart both, we can simply invoke fit as normal with coef(foce_fit3) as the initial fixed effects.
foce_fit4 = fit(
warfarin_pk_model,
pop,
coef(foce_fit3),
FOCE();
optim_options = (show_trace = false,),
)[ Info: Checking the initial parameter values. [ Info: The initial negative log-likelihood and its gradient are finite. Check passed.
FittedPumasModel
Dynamical system type: Matrix exponential
Number of subjects: 32
Observation records: Active Missing
conc: 247 51
Total: 247 51
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.46191
--------------------
Estimate
--------------------
pop_CL 0.13469
pop_Vc 8.048
pop_tabs 0.5449
pop_lag 0.87702
pk_Ω₁,₁ 0.070709
pk_Ω₂,₂ 0.018349
pk_Ω₃,₃ 0.92439
σ_prop 0.087933
σ_add 0.41375
--------------------Restarting the optimization is often a good way to avoid problems with gradient convergence.
2.5.4 A computationally expensive alternative
If gradient convergence below a certain threshold is very important it can often be achieved with GradientDescent. Gradient Descent has very good asymptotic (for iterations going to infinity) properties and avoids potential numerical issues with the inverse Hessian approximation in situations with numerical noise in the objective that could for example come from numerical integration of the dynamical system. However, since it is not using (an approximation to the) second order information, it may also take quite large steps that lead to extreme values of the parameters causing warnings to appear on screen from the ODE model solvers and it may well use a lot of iterations to reach the prescribed tolerance. Consider the following.
foce_fit5 = fit(
warfarin_pk_model,
pop,
coef(foce_fit3),
FOCE();
optim_alg = Pumas.Optim.GradientDescent(),
optim_options = (g_tol = 1e-5, show_trace = false),
)
foce_fit5.optim[ Info: Checking the initial parameter values. [ Info: The initial negative log-likelihood and its gradient are finite. Check passed.
* Status: success
* Candidate solution
Final objective value: 3.504619e+02
* Found with
Algorithm: Gradient Descent
* Convergence measures
|x - x'| = 1.56e-08 ≰ 0.0e+00
|x - x'|/|x'| = 3.91e-09 ≰ 0.0e+00
|f(x) - f(x')| = 2.81e-06 ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = 8.01e-09 ≰ 0.0e+00
|g(x)| = 8.91e-06 ≤ 1.0e-05
* Work counters
Seconds run: 9 (vs limit Inf)
Iterations: 298
f(x) calls: 340
∇f(x) calls: 299
∇f(x)ᵀv calls: 339We were indeed able to get a very low gradient tolerance, but notice the many iterations used! This may simply not be feasible for very large data sets or computationally very expensive models.
3 Calculating Some Key Diagnostics
3.1 Model Fit and Comparison
When we specify FOCE(), LaplaceI() or FO() we are doing Maximum Likelihood Estimation. This also means that the most natural number we may have to gauge and compare model fit is the log-likelihood. If we want to extract the log-likelihood at the maximum likelihood estimates we can simply use the loglikelihood function that has a method for just a FittedPumasModel.
loglikelihood(foce_fit5)-350.4619068208We see that this is the same number as in the printed output although it was obviously rounded there. We can then compare this number to the log-likelihood of another fitted model, and if they are nested we can even do a log-likelihood ratio test with lrtest. Please consult the docstring for how to use this function ?lrtest.
Other more general approaches of model comparison involves information criteria such as aic, aicc, bic, and more. These all have methods for FittedPumasModel so to calculate them simply use them as follows
aic_fit, bic_fit = aic(foce_fit5), bic(foce_fit5)(718.9238136416, 750.5083086712517)3.2 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, correlation matrices, and condition numbers.
An example usage:
inference_results = infer(foce_fit4; level = 0.95)[ Info: Calculating: variance-covariance matrix. [ Info: Done.
Asymptotic inference results using sandwich estimator
Dynamical system type: Matrix exponential
Number of subjects: 32
Observation records: Active Missing
conc: 247 51
Total: 247 51
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.46191
---------------------------------------------------------
Estimate SE 95.0% C.I.
---------------------------------------------------------
pop_CL 0.13469 0.0066655 [ 0.12163 ; 0.14776 ]
pop_Vc 8.048 0.2209 [ 7.6151 ; 8.481 ]
pop_tabs 0.5449 0.18347 [ 0.1853 ; 0.9045 ]
pop_lag 0.87702 0.052956 [ 0.77323 ; 0.98082 ]
pk_Ω₁,₁ 0.070709 0.024568 [ 0.022555 ; 0.11886 ]
pk_Ω₂,₂ 0.018349 0.0051547 [ 0.0082457; 0.028452]
pk_Ω₃,₃ 0.92439 0.40817 [ 0.12439 ; 1.7244 ]
σ_prop 0.087933 0.014414 [ 0.059682 ; 0.11619 ]
σ_add 0.41375 0.068273 [ 0.27994 ; 0.54756 ]
---------------------------------------------------------This augments the printed table from above with the uncertainty information. A useful function to know is the coefficients_table from PumasUtilities. It needs the model and inference objects as inputs and returns a DataFrame with model parameter names, descriptions from the model annotations, estimates, standard errors, relative standard errors, and confidence interval limits.
coefficients_table(foce_fit5, inference_results)9×8 DataFrame
| Row | Parameter | Description | Constant | Estimate | Standard Error | Relative_Se | CI Lower | CI Upper |
|---|---|---|---|---|---|---|---|---|
| String | String | Bool | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | pop_CL | Clearance (L/hr) | false | 0.135 | 0.007 | 0.049 | 0.122 | 0.148 |
| 2 | pop_Vc | Central volume (L) | false | 8.048 | 0.221 | 0.027 | 7.615 | 8.481 |
| 3 | pop_tabs | Absorption lag time (hr) | false | 0.545 | 0.183 | 0.337 | 0.185 | 0.904 |
| 4 | pop_lag | Lag time (hr) | false | 0.877 | 0.053 | 0.06 | 0.773 | 0.981 |
| 5 | pk_Ω₁,₁ | ΩCL: Clearance | false | 0.071 | 0.025 | 0.347 | 0.023 | 0.119 |
| 6 | pk_Ω₂,₂ | ΩVc: Central volume | false | 0.018 | 0.005 | 0.281 | 0.008 | 0.028 |
| 7 | pk_Ω₃,₃ | Ωtabs: Absorption lag time | false | 0.924 | 0.408 | 0.442 | 0.124 | 1.724 |
| 8 | σ_prop | σ_prop: Proportional error | false | 0.088 | 0.014 | 0.164 | 0.06 | 0.116 |
| 9 | σ_add | σ_add: Additive error | false | 0.414 | 0.068 | 0.165 | 0.28 | 0.548 |
In another tutorial we will look closer at different ways of calculating uncertainty for the maximum likelihood estimates.
3.3 Calculate the Condition Number
Some researchers like to include the condition number of the correlation matrix implied by the variance-covarinace matrix that was estimated by infer. This is simply calculated by using the cond function on the output of infer above.
cond(inference_results)47.21414791748396The researcher may set a limit as to how large this number can be. Models with large condition numbers are sometimes labeled “unstable” and models with small condition numbers may be labeled “stable”. This categorization is contested and discussed among practitioners and the condition number should be reported and interpreted with care outside of the linear regression context where the interpretation was originally proposed.
3.4 Shrinkage
A Non-linear Mixed Effects Model includes random effects to represent individual specific variability and well as sampling noise. In cases where the model is over-parameterized compared to the true model or to the available data we may observe what has been called shrinkage in the literature. The thing that is shrinking is the distribution of empirical bayes estimates (η-shrinkage) or in the individual weighted residuals (ϵ-shrinkage). If there is insufficient data or the model is over-parameterized it will be possible to find exact or almost exact fits for each individual. There is not enough data to move the parameters away from the prior on the random effects and the individual distributions collapse into the population distribution.
3.4.1 Random Effect (eta) shrinkage
η-shrinkage is defined as \[ ηsh = 1-\frac{sd(\eta_i)}{\omega_\eta} \] Such that if there is agreement between the prior standard deviation (\(\omega_\eta\)) and the standard deviation of \(\eta_i\) we get no shrinkage or a shrinkage of 0. As \(sd(\eta_i)\) shrinks to 0 we get shrinkages closer and closer to 1 indicating an over-parameterization or lack individual data.
In Pumas, this number is simply calculated by invoking ηshrinkage on the fitted model.
ηshrinkage(foce_fit)(pk_η = [0.008033743165338225, 0.13591927715921814, 0.4109126857830827],)This returns a NamedTuple with one element per ~ in the @random-block.
3.4.2 Error Model (sigma) shrinkage
Sometimes a number called ϵ-shrinkage is reported as well. The definition is analogous to the above. It divides a standard deviation by another standard deviation and subtracts it from zero to gauge how well they match. The problem for the derived variables is that we often have heteroscedasticity built into our models. This means that the standard deviation is a function of the model prediction such as in the proportional or combined gaussian error models. Then, there is not a single standard deviation to compare to. For this reason, we compute the individual weighted residuals (IWRES) and calculate the standard deviation of that. Since IWRES is already weighted we simply subtract the standard deviation from 1. \[
ϵsh = 1-sd(IWRES_i)
\] The interpretation is as above. We can calculate this using the ϵshrinkage in Pumas.
ϵshrinkage(foce_fit)(conc = 0.14069266539378955,)Will return a NamedTuple with one element per ~ that defines a Normal distribution in the @derived-block.
3.5 Reporting Metrics
PumasUtilities contains a convenience function called metrics_table that can be used to report a number of metrics in a compact format.
metrics_table(foce_fit)17×2 DataFrame
| Row | Metric | Value |
|---|---|---|
| String | Any | |
| 1 | Estimation Time | 11.74 |
| 2 | Subjects | 32 |
| 3 | Fixed Parameters | 0 |
| 4 | Optimized Parameters | 9 |
| 5 | conc Active Observations | 247 |
| 6 | conc Missing Observations | 51 |
| 7 | Total Active Observations | 247 |
| 8 | Total Missing Observations | 51 |
| 9 | Likelihood Approximation | FOCE |
| 10 | LogLikelihood (LL) | -350.462 |
| 11 | -2LL | 700.924 |
| 12 | AIC | 718.924 |
| 13 | BIC | 750.508 |
| 14 | (η-shrinkage) pk_η₁ | 0.008 |
| 15 | (η-shrinkage) pk_η₂ | 0.136 |
| 16 | (η-shrinkage) pk_η₃ | 0.411 |
| 17 | (ϵ-shrinkage) conc | 0.141 |
This will return a DataFrame with the metrics.
4 Concluding Remarks
This tutorial showcased a subset of the steps users want to perform after a model fit:
- Calculate uncertainty of the estimated parameters
- Calculate log-likelihood and information criteria for fit assessment and model comparison purposes
- Calculate shrinkage numbers to assess over-parametrization or excessive sparseness of data
Other tutorials will cover many more diagnostics steps to be done to asses the quality of the fitted model.
References
Mogensen, P, and A Riseth. 2018. “Optim: A Mathematical Optimization Package for Julia.” Journal of Open Source Software 3 (24).