using Pumas
using PharmaDatasets
using DataFramesMeta
using PumasUtilities
Calculating Parameter Uncertainty
1 Introduction
A typical workflow for fitting a Pumas model and deriving parameter precision typically involves:
- Preparing the data and the model.
- Checking model-data compatibility.
- Obtaining initial parameter estimates.
- Fitting the model via a chosen estimation method.
- Interpreting the fit results.
- Computing parameter precision.
- (Optionally) proceeding with more advanced techniques like bootstrapping or SIR for robust uncertainty quantification.
The following sections will walk through these steps using a one-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_tabs, pop_lag
- Inter-individual variability (IIV) components:
pk_Ω
- Residual error model parameters:
σ_prop
,σ_add
- Covariates for scaling:
FSZCL
andFSZV
- Differential equations describing the PK behavior in the compartments.
= @model begin
warfarin_pk_model @metadata begin
= "Warfarin 1-compartment PK model (PD removed)"
desc = u"hr"
timeu end
@param begin
# PK parameters
"""
Clearance (L/hr)
"""
∈ RealDomain(lower = 0.0, init = 0.134)
pop_CL """
Central volume (L)
"""
∈ RealDomain(lower = 0.0, init = 8.11)
pop_Vc """
Absorption lag time (hr)
"""
∈ RealDomain(lower = 0.0, init = 0.523)
pop_tabs """
Lag time (hr)
"""
∈ RealDomain(lower = 0.0, init = 0.1)
pop_lag # Inter-individual variability
"""
- ΩCL: Clearance
- ΩVc: Central volume
- Ωtabs: Absorption lag time
"""
∈ PDiagDomain([0.01, 0.01, 0.01])
pk_Ω # Residual variability
"""
σ_prop: Proportional error
"""
∈ RealDomain(lower = 0.0, init = 0.00752)
σ_prop """
σ_add: Additive error
"""
∈ RealDomain(lower = 0.0, init = 0.0661)
σ_add end
@random begin
~ MvNormal(pk_Ω) # mean = 0, covariance = pk_Ω
pk_η end
@covariates begin
"""
FSZCL: Clearance scaling factor
"""
FSZCL"""
FSZV: Volume scaling factor
"""
FSZVend
@pre begin
= FSZCL * pop_CL * exp(pk_η[1])
CL = FSZV * pop_Vc * exp(pk_η[2])
Vc = pop_tabs * exp(pk_η[3])
tabs = log(2) / tabs
Ka end
@dosecontrol begin
= (Depot = pop_lag,)
lags end
@vars begin
:= Central / Vc
cp end
@dynamics Depots1Central1
@derived begin
"""
Concentration (ng/mL)
"""
~ @. Normal(cp, sqrt((σ_prop * cp)^2 + σ_add^2))
conc end
end
PumasModel
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: Closed form
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.
= dataset("pumas/warfarin_pumas")
warfarin_data
# Transform the data in a single chain of operations
= @chain warfarin_data begin
warfarin_data_wide @rtransform begin
# Scaling factors
:FSZV = :wtbl / 70 # volume scaling
:FSZCL = (:wtbl / 70)^0.75 # clearance scaling (allometric)
end
end
Row | id | time | evid | amt | cmt | conc | pca | wtbl | age | sex | FSZV | FSZCL |
---|---|---|---|---|---|---|---|---|---|---|---|---|
Int64 | Float64 | Int64 | Float64? | Int64? | Float64? | Float64? | Float64 | Int64 | String1 | Float64 | Float64 | |
1 | 1 | 0.0 | 1 | 100.0 | 1 | missing | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
2 | 1 | 0.5 | 0 | missing | missing | 0.0 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
3 | 1 | 1.0 | 0 | missing | missing | 1.9 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
4 | 1 | 2.0 | 0 | missing | missing | 3.3 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
5 | 1 | 3.0 | 0 | missing | missing | 6.6 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
6 | 1 | 6.0 | 0 | missing | missing | 9.1 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
7 | 1 | 9.0 | 0 | missing | missing | 10.8 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
8 | 1 | 12.0 | 0 | missing | missing | 8.6 | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
9 | 1 | 24.0 | 0 | missing | missing | 5.6 | 44.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
10 | 1 | 36.0 | 0 | missing | missing | 4.0 | 27.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
11 | 1 | 48.0 | 0 | missing | missing | 2.7 | 28.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
12 | 1 | 72.0 | 0 | missing | missing | 0.8 | 31.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
13 | 1 | 96.0 | 0 | missing | missing | missing | 60.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
319 | 32 | 48.0 | 0 | missing | missing | 6.9 | 24.0 | 62.0 | 21 | M | 0.885714 | 0.912999 |
320 | 32 | 72.0 | 0 | missing | missing | 4.4 | 23.0 | 62.0 | 21 | M | 0.885714 | 0.912999 |
321 | 32 | 96.0 | 0 | missing | missing | 3.5 | 20.0 | 62.0 | 21 | M | 0.885714 | 0.912999 |
322 | 32 | 120.0 | 0 | missing | missing | 2.5 | 22.0 | 62.0 | 21 | M | 0.885714 | 0.912999 |
323 | 33 | 0.0 | 1 | 100.0 | 1 | missing | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
324 | 33 | 0.0 | 0 | missing | missing | missing | 100.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
325 | 33 | 24.0 | 0 | missing | missing | 9.2 | 49.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
326 | 33 | 36.0 | 0 | missing | missing | 8.5 | 32.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
327 | 33 | 48.0 | 0 | missing | missing | 6.4 | 26.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
328 | 33 | 72.0 | 0 | missing | missing | 4.8 | 22.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
329 | 33 | 96.0 | 0 | missing | missing | 3.1 | 28.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
330 | 33 | 120.0 | 0 | missing | missing | 2.5 | 33.0 | 66.7 | 50 | M | 0.952857 | 0.96443 |
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:
= read_pumas(
pop_pk
warfarin_data_wide;= :id,
id = :time,
time = :amt,
amt = :cmt,
cmt = :evid,
evid = [:sex, :wtbl, :FSZV, :FSZCL],
covariates = [:conc],
observations )
Population
Subjects: 32
Covariates: sex, wtbl, FSZV, FSZCL
Observations: conc
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.
Also note, that parameter inference can be expensive and for that reason we simplified the model for this tutorial to decrease overall runtime.
3.1 Checking Model-Data Compatibility
Before performing any fit, it is recommended to verify whether the defined model can handle the provided dataset. Pumas offers functions such as loglikelihood
and findinfluential
for these checks.
3.1.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 includeFOCE()
,FO()
,LaplaceI()
, etc.
For example, one might write:
# A named tuple of parameter values
= (
param_vals = 0.134,
pop_CL = 8.11,
pop_Vc = 0.523,
pop_tabs = 0.1,
pop_lag = Diagonal([0.01, 0.01, 0.01]),
pk_Ω = 0.00752,
σ_prop = 0.0661,
σ_add
)= fit(warfarin_pk_model, pop_pk, param_vals, FOCE();) foce_fit
[ 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.209064e+04 1.489225e+04 * time: 0.03675198554992676 1 2.643772e+03 3.167516e+03 * time: 1.4747381210327148 2 1.836601e+03 2.118430e+03 * time: 1.7526400089263916 3 9.351337e+02 8.722439e+02 * time: 1.7702059745788574 4 6.402300e+02 4.199225e+02 * time: 1.7873950004577637 5 5.103664e+02 1.642121e+02 * time: 1.8039441108703613 6 4.760464e+02 5.453749e+01 * time: 1.820162057876587 7 4.703757e+02 3.643518e+01 * time: 1.8363800048828125 8 4.699019e+02 3.135992e+01 * time: 1.8520629405975342 9 4.697614e+02 2.953531e+01 * time: 1.8669240474700928 10 4.693153e+02 2.463233e+01 * time: 1.8822100162506104 11 4.685743e+02 2.580427e+01 * time: 1.8974249362945557 12 4.675133e+02 3.864937e+01 * time: 1.9128310680389404 13 4.666775e+02 5.495470e+01 * time: 1.9279429912567139 14 4.661197e+02 5.692101e+01 * time: 1.943000078201294 15 4.656782e+02 4.770992e+01 * time: 1.9578680992126465 16 4.651802e+02 3.087698e+01 * time: 1.9735970497131348 17 4.645523e+02 1.184834e+01 * time: 1.991873025894165 18 4.641447e+02 1.162249e+01 * time: 2.0895190238952637 19 4.639978e+02 1.125144e+01 * time: 2.1046030521392822 20 4.639307e+02 1.156463e+01 * time: 2.120753049850464 21 4.638001e+02 1.312870e+01 * time: 2.1357409954071045 22 4.635282e+02 1.480920e+01 * time: 2.1507019996643066 23 4.630353e+02 2.169377e+01 * time: 2.1660830974578857 24 4.623847e+02 4.478029e+01 * time: 2.1816020011901855 25 4.617426e+02 6.468975e+01 * time: 2.1972689628601074 26 4.610293e+02 7.776996e+01 * time: 2.2128751277923584 27 4.597628e+02 8.785260e+01 * time: 2.2286999225616455 28 4.566753e+02 9.769803e+01 * time: 2.245326042175293 29 4.490421e+02 1.008838e+02 * time: 2.2656331062316895 30 4.391868e+02 9.978816e+01 * time: 2.286550998687744 31 4.130704e+02 5.917685e+01 * time: 2.3081319332122803 32 4.055780e+02 3.852824e+01 * time: 2.3519699573516846 33 4.023118e+02 3.889618e+01 * time: 2.369731903076172 34 4.012516e+02 3.694778e+01 * time: 2.3940329551696777 35 4.004391e+02 2.061948e+01 * time: 2.4114699363708496 36 3.983040e+02 3.508423e+01 * time: 2.4288880825042725 37 3.969705e+02 3.841039e+01 * time: 2.4470369815826416 38 3.965462e+02 3.738343e+01 * time: 2.464261054992676 39 3.950409e+02 3.064789e+01 * time: 2.4817941188812256 40 3.945750e+02 2.876429e+01 * time: 2.4992449283599854 41 3.937725e+02 2.571438e+01 * time: 2.5170140266418457 42 3.933955e+02 2.436112e+01 * time: 2.5357470512390137 43 3.927564e+02 2.051069e+01 * time: 2.5753109455108643 44 3.916020e+02 1.629035e+01 * time: 2.5926859378814697 45 3.886991e+02 2.689824e+01 * time: 2.610496997833252 46 3.870054e+02 2.298582e+01 * time: 2.62937593460083 47 3.853691e+02 2.614992e+01 * time: 2.6474111080169678 48 3.841730e+02 2.207557e+01 * time: 2.665019989013672 49 3.825113e+02 2.204399e+01 * time: 2.683074951171875 50 3.808880e+02 2.444784e+01 * time: 2.7004339694976807 51 3.800407e+02 1.250611e+01 * time: 2.717205047607422 52 3.798092e+02 1.167926e+01 * time: 2.7336010932922363 53 3.797789e+02 1.162382e+01 * time: 2.7675321102142334 54 3.797069e+02 1.152441e+01 * time: 2.7833831310272217 55 3.794424e+02 1.132717e+01 * time: 2.799989938735962 56 3.788131e+02 2.006438e+01 * time: 2.8164401054382324 57 3.771525e+02 3.584695e+01 * time: 2.833604097366333 58 3.731299e+02 5.697249e+01 * time: 2.8507890701293945 59 3.658671e+02 6.542042e+01 * time: 2.8686940670013428 60 3.604194e+02 4.036489e+01 * time: 2.8864059448242188 61 3.532841e+02 1.574006e+01 * time: 2.918376922607422 62 3.520181e+02 1.393300e+01 * time: 2.935529947280884 63 3.517984e+02 6.701188e+00 * time: 2.9518630504608154 64 3.517541e+02 3.503978e+00 * time: 2.9681029319763184 65 3.516436e+02 8.720957e+00 * time: 2.9857680797576904 66 3.511845e+02 1.406200e+01 * time: 3.0082690715789795 67 3.510647e+02 2.540378e+00 * time: 3.038135051727295 68 3.510209e+02 3.157201e+00 * time: 3.054983139038086 69 3.509959e+02 3.045642e+00 * time: 3.0705180168151855 70 3.509765e+02 2.673143e+00 * time: 3.0861690044403076 71 3.509751e+02 2.603975e+00 * time: 3.1010830402374268 72 3.509724e+02 2.505719e+00 * time: 3.117516040802002 73 3.509666e+02 2.379768e+00 * time: 3.144541025161743 74 3.509504e+02 3.572030e+00 * time: 3.160053014755249 75 3.509123e+02 6.006350e+00 * time: 3.175511121749878 76 3.508288e+02 8.822995e+00 * time: 3.1909849643707275 77 3.506944e+02 9.708012e+00 * time: 3.2066640853881836 78 3.505767e+02 6.092631e+00 * time: 3.222661018371582 79 3.505358e+02 1.734431e+00 * time: 3.2509260177612305 80 3.505314e+02 6.749379e-01 * time: 3.2661449909210205 81 3.505313e+02 6.721982e-01 * time: 3.2807891368865967 82 3.505312e+02 6.699487e-01 * time: 3.295156955718994 83 3.505307e+02 6.606824e-01 * time: 3.3096511363983154 84 3.505298e+02 6.413909e-01 * time: 3.3243799209594727 85 3.505274e+02 9.083363e-01 * time: 3.351072072982788 86 3.505222e+02 1.339147e+00 * time: 3.366420030593872 87 3.505129e+02 1.608661e+00 * time: 3.3817520141601562 88 3.505026e+02 1.293164e+00 * time: 3.396914005279541 89 3.504973e+02 5.140504e-01 * time: 3.4122350215911865 90 3.504963e+02 6.340189e-02 * time: 3.4274189472198486 91 3.504963e+02 3.137914e-03 * time: 3.454767942428589 92 3.504963e+02 5.681551e-04 * time: 3.4686319828033447
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 32
Observation records: Active Missing
conc: 251 47
Total: 251 47
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.49625
--------------------
Estimate
--------------------
pop_CL 0.13465
pop_Vc 8.0535
pop_tabs 0.55061
pop_lag 0.87158
pk_Ω₁,₁ 0.070642
pk_Ω₂,₂ 0.018302
pk_Ω₃,₃ 0.91326
σ_prop 0.090096
σ_add 0.39115
--------------------
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, and correlation matrices.
The signature for infer
often looks like:
infer(
::FittedPumasModel;
fpm= 0.95,
level ::Bool = false,
rethrow_error::Bool = true,
sandwich_estimator )
where:
fpm::FittedPumasModel
: The result offit
(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 parametersrethrow_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 totrue
(the default value), the sandwich estimator will be used. If set tofalse
, the standard error will be calculated using the inverse of the Hessian matrix.
An example usage:
= infer(foce_fit; level = 0.95) inference_results
[ Info: Calculating: variance-covariance matrix. [ Info: Done.
Asymptotic inference results using sandwich estimator
Dynamical system type: Closed form
Number of subjects: 32
Observation records: Active Missing
conc: 251 47
Total: 251 47
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.49625
---------------------------------------------------------
Estimate SE 95.0% C.I.
---------------------------------------------------------
pop_CL 0.13465 0.0066546 [ 0.12161 ; 0.1477 ]
pop_Vc 8.0535 0.22108 [ 7.6201 ; 8.4868 ]
pop_tabs 0.55061 0.18702 [ 0.18406 ; 0.91717 ]
pop_lag 0.87158 0.056687 [ 0.76048 ; 0.98269 ]
pk_Ω₁,₁ 0.070642 0.024577 [ 0.022472 ; 0.11881 ]
pk_Ω₂,₂ 0.018302 0.0051549 [ 0.0081988; 0.028406]
pk_Ω₃,₃ 0.91326 0.40637 [ 0.11678 ; 1.7097 ]
σ_prop 0.090096 0.014521 [ 0.061636 ; 0.11856 ]
σ_add 0.39115 0.065398 [ 0.26297 ; 0.51932 ]
---------------------------------------------------------
This is the usual asymptotic variance-covariance estimator and we already saw this previous tutorials.
To get a matrix representation of this, use vcov()
vcov(inference_results)
9×9 Symmetric{Float64, Matrix{Float64}}:
4.42841e-5 0.000217445 0.000302094 … 3.99916e-6 0.00019736
0.000217445 0.0488775 0.00571323 -0.000846166 -0.0056657
0.000302094 0.00571323 0.0349767 0.000227818 0.00412692
-7.40855e-5 -0.00207014 -0.00450616 0.000458813 0.000494683
0.000120614 5.09406e-5 0.00164596 -9.1424e-5 0.000734901
2.90008e-7 0.000292148 -0.000131446 … 3.99746e-6 1.80866e-5
-0.000263152 -0.023877 -0.0275659 0.00328879 0.0126135
3.99916e-6 -0.000846166 0.000227818 0.000210856 0.000518153
0.00019736 -0.0056657 0.00412692 0.000518153 0.00427687
and to get the condition number of the correlation matrix implied by vcov
use
cond(inference_results)
50.11623683487897
Some users request the condition number of the covariance matrix itself and even if the use is often misguided it can be calculated as well.
cond(inference_results; correlation = false)
13082.75373321667
It is also possible to calculate the correlation matrix from the vcov
output using the cov2cor
function
= cov2cor(vcov(inference_results)) cor_from_cov
9×9 Symmetric{Float64, Matrix{Float64}}:
1.0 0.147799 0.242733 … -0.0973098 0.0413859 0.453494
0.147799 1.0 0.138178 -0.265766 -0.263578 -0.391865
0.242733 0.138178 1.0 -0.362707 0.083889 0.337422
-0.196394 -0.165183 -0.425047 0.555027 0.557394 0.133439
0.737483 0.00937536 0.358102 -0.28125 -0.25618 0.45724
0.00845409 0.256348 -0.136345 … 0.315212 0.0534038 0.0536508
-0.0973098 -0.265766 -0.362707 1.0 0.557335 0.47462
0.0413859 -0.263578 0.083889 0.557335 1.0 0.545635
0.453494 -0.391865 0.337422 0.47462 0.545635 1.0
And we see that the cond
call above matches the condition number of the correlation matrix
cond(cor_from_cov)
50.1162368348788
3.2.1 Failure of the asymptotic variance-covariance matrix
It is well-known that the asymptotic variance-covariance matrix can sometimes fail to compute. This can happen for a variety of reasons including:
- There are parameters very close to a bound (often 0)
- The parameter vector does not represent a local minimum (optimization failed)
- The parameter vector does represent a local minimum but it’s not the global solution
The first one is often easy to check. The problematic parameters can be ones than have lower or upper bounds set. Often this will be a variance of standard deviation that has moved very close to the lower boundary. Consider removing the associated random effect if the problematic parameter is a variance in its specification or error model component if a combined additive and proportional error model is used and a standard deviation is close to zero.
It is also possible that the parameters do not represent a local minimum. In other words, they come from a failed fit. In this case, it can often be hard to perform the associated calculations in a stable way, but most importantly the results would not be interpretable even if they can be calculated in this case. The formulas are only valid for parameters that represent the actual (maximum likelihood) estimates. Please try to restart the optimization at different starting points in this case.
If you have reasons to believe that the model should indeed be a combined error model or if the random effect should be present it is also possible that the model converged to a local minimum that is not the global minimum. If the optimization happened to move to a region of the parameter state space that is hard to get out of you will often have to restart the fit at different parameter values. It is not possible to verify if the minimum is global in general, but it is always advised to try out more than one set of initial parameters when fitting models.
3.2.2 Bootstrap
Sometimes it is appropriate to use a different method to calculate estimate the uncertainty of the estimated parameters. Bootstrapping is a very popular approach that is simply but can often come a quite significant computational cost. Researchers often perform a bootstrapping step if their computational budget allows it or if the asymptotic variance-covariance estimator fails. Bootstrapping is advantageous because it does not rely on any invertability of matrices and it cannot produce negative variance confidence intervals because the resampled estimator respects the bounds of the model.
The signature for bootstrapping in infer
looks as follows.
infer(fpm::FittedPumasModel, bts::Bootstrap; level = 0.95)
This does not help much before also looking at the interface for Bootstrap
itself.
Bootstrap(;
= Random.default_rng,
rng = 200,
samples = nothing,
stratify_by = EnsembleThreads(),
ensemblealg )
Bootstrap
accepts a random number generator rng
, the number of resampled datasets to produce samples
, if sampling should be stratified according to the covariates in stratify_by
, and finally the ensemble algorithm to control parallelization across fits. On the JuliaHub platform this can be used together with distributed computing to perform many resampled estimations in a short time.
= infer(foce_fit, Bootstrap(samples = 50); level = 0.95) bootstrap_results
┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Warning: Terminated early due to NaN in gradient. └ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/7krni/src/multivariate/optimize/optimize.jl:117 ┌ Info: Bootstrap inference finished. │ Total resampled fits = 50 │ Success rate = 1.0 └ Unique resampled populations = 50
Bootstrap inference results
Dynamical system type: Closed form
Number of subjects: 32
Observation records: Active Missing
conc: 251 47
Total: 251 47
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.49625
---------------------------------------------------------
Estimate SE 95.0% C.I.
---------------------------------------------------------
pop_CL 0.13465 0.0090551 [ 0.12428 ; 0.14912 ]
pop_Vc 8.0535 0.26504 [ 7.6274 ; 8.5159 ]
pop_tabs 0.55061 0.22486 [ 0.21093 ; 1.085 ]
pop_lag 0.87158 0.14698 [ 0.58479 ; 1.2821 ]
pk_Ω₁,₁ 0.070642 561220.0 [ 0.028144 ; 0.11522 ]
pk_Ω₂,₂ 0.018302 0.004948 [ 0.0095219; 0.027023]
pk_Ω₃,₃ 0.91326 0.79202 [ 0.083747 ; 3.0749 ]
σ_prop 0.090096 0.014425 [ 0.070673 ; 0.11906 ]
σ_add 0.39115 0.10914 [ 0.036122 ; 0.47577 ]
---------------------------------------------------------
Successful fits: 50 out of 50
Unique resampled populations: 50
No stratification.
Again, we can calculate a covariance matrix based on the samples with vcov
vcov(bootstrap_results)
9×9 Symmetric{Float64, Matrix{Float64}}:
8.19945e-5 -0.000336711 … -3.50055e-5 0.000565455
-0.000336711 0.0702444 -0.000894007 -0.0176072
0.000558086 0.0178959 -0.000197232 0.00387331
-0.000111521 -0.00466242 0.000512402 0.000357513
-2951.78 56078.7 2709.21 -27779.4
-8.40975e-6 0.000347322 … 1.4446e-5 -7.80706e-5
-0.000171646 -0.0827463 0.00520728 0.0267877
-3.50055e-5 -0.000894007 0.000208079 2.8332e-5
0.000565455 -0.0176072 2.8332e-5 0.0119106
and we can even get a DataFrame
that includes all the estimated parameters from the sampled population fits
DataFrame(bootstrap_results.vcov)
Row | pop_CL | pop_Vc | pop_tabs | pop_lag | pk_Ω₁,₁ | pk_Ω₂,₂ | pk_Ω₃,₃ | σ_prop | σ_add |
---|---|---|---|---|---|---|---|---|---|
Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
1 | 0.13323 | 8.20728 | 0.791375 | 0.561339 | 0.0823593 | 0.0136195 | 0.318462 | 0.100694 | 0.232389 |
2 | 0.135569 | 8.44134 | 0.92093 | 0.84089 | 0.0878481 | 0.0223436 | 0.434904 | 0.0793278 | 0.379577 |
3 | 0.133214 | 8.15451 | 0.209688 | 0.946414 | 0.0548924 | 0.0191834 | 1.43885 | 0.0821152 | 0.303382 |
4 | 0.146243 | 8.09525 | 0.680989 | 0.848553 | 0.0948266 | 0.0112764 | 0.664869 | 0.0992746 | 0.379791 |
5 | 0.124129 | 7.73684 | 0.441653 | 0.842237 | 0.0444963 | 0.0117525 | 1.15752 | 0.0921815 | 0.397297 |
6 | 0.12592 | 8.25691 | 0.635714 | 0.836178 | 0.0550814 | 0.014361 | 0.458466 | 0.0753128 | 0.298527 |
7 | 0.126899 | 7.78716 | 0.439159 | 0.972769 | 0.0385816 | 0.0224259 | 1.9789 | 0.114202 | 0.420057 |
8 | 0.130988 | 8.07576 | 0.311589 | 0.98284 | 0.0282472 | 0.0253044 | 3.4725 | 0.112748 | 0.294359 |
9 | 0.136263 | 8.37749 | 0.454351 | 0.870117 | 0.0382458 | 0.0208256 | 0.727583 | 0.095338 | 0.240551 |
10 | 0.134155 | 7.54576 | 0.279517 | 0.973944 | 0.0564424 | 0.0197273 | 3.07544 | 0.108975 | 0.469758 |
11 | 0.130906 | 8.2373 | 0.564182 | 0.862791 | 0.0450941 | 0.0180138 | 0.364845 | 0.0973565 | 0.213472 |
12 | 0.143213 | 8.47163 | 0.552519 | 0.889574 | 0.0876238 | 0.0155491 | 0.39629 | 0.061615 | 0.338408 |
13 | 0.0971723 | 8.77843 | 0.675515 | 0.710536 | 3.96844e6 | 0.025277 | 0.0204085 | 0.126042 | 1.10943e-82 |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
39 | 0.140619 | 8.4742 | 1.09293 | 0.787401 | 0.0878795 | 0.0235782 | 0.393446 | 0.0779579 | 0.383531 |
40 | 0.128871 | 8.43283 | 0.21521 | 0.91707 | 0.0429032 | 0.0176422 | 0.668267 | 0.0780908 | 0.186376 |
41 | 0.125697 | 7.91828 | 0.480461 | 0.919592 | 0.0336388 | 0.0149118 | 0.664353 | 0.0994008 | 0.316135 |
42 | 0.140551 | 7.78494 | 0.549585 | 0.928 | 0.0883875 | 0.0147415 | 0.898622 | 0.10708 | 0.424709 |
43 | 0.144365 | 8.08573 | 0.5153 | 0.87089 | 0.0971503 | 0.014082 | 0.488065 | 0.0844205 | 0.398055 |
44 | 0.133604 | 8.38586 | 0.260778 | 0.972484 | 0.050185 | 0.0255931 | 2.89587 | 0.101091 | 0.290101 |
45 | 0.127517 | 7.91961 | 0.428351 | 0.976971 | 0.0371348 | 0.0272688 | 3.07325 | 0.119089 | 0.441089 |
46 | 0.124923 | 8.51307 | 0.640704 | 0.855135 | 0.0286727 | 0.0176239 | 0.342263 | 0.0818859 | 0.160004 |
47 | 0.135456 | 8.18322 | 0.742447 | 0.904681 | 0.0676346 | 0.0180688 | 0.496355 | 0.0966218 | 0.440981 |
48 | 0.135845 | 7.69767 | 0.536743 | 0.915432 | 0.0840448 | 0.0181668 | 0.948525 | 0.0969929 | 0.424724 |
49 | 0.126358 | 8.19263 | 0.682 | 0.841587 | 0.0470215 | 0.0145047 | 1.02964 | 0.0787489 | 0.358142 |
50 | 0.128161 | 7.76568 | 0.436558 | 0.546353 | 0.0614043 | 0.0164034 | 3.73374e-7 | 0.0736364 | 0.307656 |
This is very useful for histogram plotting of parameter distributions.
3.2.3 Sampling Importance Re-sampling
Pumas has support for inference through Sampling Importance Re-sampling through the SIR()
input to infer
. The signature for SIR in infer
looks as follows.
infer(fpm::FittedPumasModel, sir::SIR; level = 0.95, ensemblealg = EnsembleThreads())
This performs sampling importance re-sampling for the population in fpm
. The confidence intervals are calculated as the (1-level)/2
and (1+level)/2
quantiles of the sampled parameters. ensemblealg
can be EnsembleThreads() (the default value) to use multi-threading or EnsembleSerial()
to use a single thread.
The signature for the SIR
specification is
SIR(; rng, samples, resamples)
SIR
accepts a random number generator rng
, the number of samples from the proposal, samples
, can be set and to complete the specification the resample
has to be set. It is suggested that samples
is at least 5 times larger than resamples
in practice to have sufficient samples to resample from.
= infer(foce_fit, SIR(samples = 1000, resamples = 200); level = 0.95) sir_results
[ Info: Calculating: variance-covariance matrix. [ Info: Done. [ Info: Running SIR. [ Info: Resampling.
Simulated inference results
Dynamical system type: Closed form
Number of subjects: 32
Observation records: Active Missing
conc: 251 47
Total: 251 47
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.49625
--------------------------------------------------------
Estimate SE 95.0% C.I.
--------------------------------------------------------
pop_CL 0.13465 0.0055232 [ 0.12578 ; 0.14599 ]
pop_Vc 8.0535 0.21377 [ 7.6726 ; 8.4502 ]
pop_tabs 0.55061 0.15057 [ 0.29822 ; 0.89115 ]
pop_lag 0.87158 0.042053 [ 0.78352 ; 0.94544 ]
pk_Ω₁,₁ 0.070642 0.016553 [ 0.045459; 0.10976 ]
pk_Ω₂,₂ 0.018302 0.0048786 [ 0.010063; 0.028944]
pk_Ω₃,₃ 0.91326 0.29955 [ 0.50762 ; 1.6117 ]
σ_prop 0.090096 0.0085777 [ 0.074035; 0.10566 ]
σ_add 0.39115 0.036412 [ 0.33811 ; 0.47912 ]
--------------------------------------------------------
Notice, that SIR
bases its first samples
number of samples from a truncated multivariate normal distribution with mean of the maximum likelihood population level parameters and covariance matrix that is the asymptotic matrix calculated by infer(fpm)
. This means that to use SIR
the matrix is question has to be successfully calculated by infer(fpm)
under the hood.
The methods for vcov
and DataFrame(sir_results.vcov)
that we saw for Bootstrap
also applies here
vcov(sir_results)
9×9 Symmetric{Float64, Matrix{Float64}}:
3.05057e-5 0.00017376 3.16741e-5 … -1.08916e-6 3.89227e-5
0.00017376 0.0456974 0.00910612 3.82711e-7 -0.002038
3.16741e-5 0.00910612 0.0226718 9.46417e-5 0.000396189
-1.51387e-5 -0.000610178 -0.00206591 2.29116e-5 0.000309257
4.8686e-5 0.000280504 0.000245383 -5.19439e-5 0.000251573
-2.90886e-6 0.000189743 4.79939e-5 … 7.45685e-7 8.72613e-6
-0.000366786 -0.00720379 -0.0174358 0.000261838 0.00316718
-1.08916e-6 3.82711e-7 9.46417e-5 7.35768e-5 2.44909e-5
3.89227e-5 -0.002038 0.000396189 2.44909e-5 0.00132585
and
DataFrame(sir_results.vcov)
Row | pop_CL | pop_Vc | pop_tabs | pop_lag | pk_Ω₁,₁ | pk_Ω₂,₂ | pk_Ω₃,₃ | σ_prop | σ_add |
---|---|---|---|---|---|---|---|---|---|
Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
1 | 0.135756 | 8.18513 | 0.457613 | 0.862083 | 0.0676902 | 0.0162232 | 0.904281 | 0.0978225 | 0.369858 |
2 | 0.127934 | 8.42382 | 0.952258 | 0.774151 | 0.0795748 | 0.0221921 | 0.513181 | 0.0886247 | 0.346316 |
3 | 0.127604 | 8.08186 | 0.641844 | 0.93448 | 0.083258 | 0.0276672 | 1.52666 | 0.0910196 | 0.490284 |
4 | 0.145313 | 8.2139 | 0.566593 | 0.85655 | 0.0744819 | 0.0291473 | 0.935252 | 0.0888907 | 0.351748 |
5 | 0.126793 | 7.97982 | 0.419822 | 0.743315 | 0.0623115 | 0.018158 | 1.12384 | 0.0913996 | 0.389919 |
6 | 0.135247 | 7.95475 | 0.609286 | 0.896347 | 0.0525052 | 0.0225094 | 1.39235 | 0.0825916 | 0.360687 |
7 | 0.126216 | 7.73592 | 0.530339 | 0.89951 | 0.0732214 | 0.0207412 | 1.3877 | 0.0737389 | 0.462568 |
8 | 0.139774 | 8.23812 | 0.530381 | 0.910542 | 0.090941 | 0.017896 | 1.45225 | 0.0895194 | 0.435131 |
9 | 0.13722 | 8.16592 | 0.621373 | 0.853793 | 0.0759836 | 0.0242841 | 0.618417 | 0.0849662 | 0.338282 |
10 | 0.141234 | 7.66084 | 0.596256 | 0.88745 | 0.103689 | 0.01412 | 0.637426 | 0.0820216 | 0.381504 |
11 | 0.144751 | 7.87291 | 0.497025 | 0.87971 | 0.0737385 | 0.018409 | 1.12892 | 0.0864012 | 0.415998 |
12 | 0.129264 | 7.97477 | 0.522049 | 0.869903 | 0.0670709 | 0.00955507 | 0.571405 | 0.0896627 | 0.406787 |
13 | 0.138135 | 7.82616 | 0.527258 | 0.867346 | 0.0738438 | 0.0171755 | 0.750477 | 0.102415 | 0.36373 |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
189 | 0.137869 | 8.19329 | 0.299244 | 0.945433 | 0.0801471 | 0.0172385 | 0.799525 | 0.0948349 | 0.382412 |
190 | 0.131889 | 7.76434 | 0.233769 | 0.86404 | 0.0485432 | 0.0195741 | 1.79979 | 0.101398 | 0.415065 |
191 | 0.139349 | 8.23077 | 0.338273 | 0.91902 | 0.1022 | 0.0267065 | 1.07595 | 0.0723795 | 0.39755 |
192 | 0.128974 | 8.43291 | 0.699697 | 0.899561 | 0.0809745 | 0.0251594 | 0.843238 | 0.0754693 | 0.370244 |
193 | 0.143735 | 7.95856 | 0.624922 | 0.896066 | 0.0918681 | 0.0124163 | 0.626936 | 0.0887354 | 0.454116 |
194 | 0.132335 | 7.9153 | 0.455948 | 0.911961 | 0.0513684 | 0.0209757 | 1.27935 | 0.0910339 | 0.410635 |
195 | 0.141326 | 8.16733 | 0.656509 | 0.878783 | 0.0713747 | 0.0233366 | 1.2275 | 0.100199 | 0.421034 |
196 | 0.121969 | 8.01575 | 0.546848 | 0.918101 | 0.0591082 | 0.0100144 | 0.668178 | 0.0836779 | 0.373118 |
197 | 0.135153 | 7.91359 | 0.718934 | 0.806777 | 0.0720071 | 0.00968338 | 0.986095 | 0.105033 | 0.450284 |
198 | 0.136877 | 7.79438 | 0.522326 | 0.865159 | 0.0891616 | 0.0128566 | 0.901763 | 0.0854831 | 0.440135 |
199 | 0.143798 | 8.01544 | 0.548749 | 0.899228 | 0.084835 | 0.0167247 | 0.608901 | 0.090777 | 0.373111 |
200 | 0.139225 | 7.96604 | 0.35164 | 0.870921 | 0.0849718 | 0.0160526 | 0.834434 | 0.0928043 | 0.392189 |
3.2.4 Marginal MCMC
An alternative to Bootstrap
and SIR
is to simply use the MarginalMCMC
sampler which is a Hamiltonian Monte Carlo (HMC) No-U-Turn Sampler (NUTS) that will sample from the marginal loglikelihood. This means that individual effects are marginalized out and then we sample the population level parameters. This does not resample populations like Bootstrap
so inference may be more stable if many resampled populations lead to extreme estimates and it differs from SIR
in that it does not need the asymptotic covariance matrix to be calculated and sampled from.
This method requires slightly more understanding from the user when setting the options that can be found through the docstring of MarginalMCMC
. Some knowledge of Bayesian inference is advised.
= infer(foce_fit, MarginalMCMC(); level = 0.95) inference_results
As sampling based inference can be computationally intensive we exclude the actual invocation of this method from this tutorial.
4 Concluding Remarks
This tutorial showcased a typical Pumas workflow for parameter inference in Pumas models. We showed the different methods supported by Pumas for calculating parameter uncertainty.