Calculating Parameter Uncertainty

Author

Patrick Kofod Mogensen

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 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 and FSZV
  • Differential equations describing the PK behavior in the compartments.
using Pumas
using PharmaDatasets
using DataFramesMeta
using PumasUtilities
using StatsBase # for cov2cor
warfarin_pk_model = @model begin
    @metadata begin
        desc = "Warfarin 1-compartment PK model (PD removed)"
        timeu = u"hr"
    end

    @param begin
        # PK parameters
        """
        Clearance (L/hr)
        """
        pop_CL  RealDomain(lower = 0.0, upper = 13.4, init = 0.134)
        """
        Central volume (L)
        """
        pop_Vc  RealDomain(lower = 0.0, upper = 81.1, init = 8.11)
        """
        Absorption lag time (hr)
        """
        pop_tabs  RealDomain(lower = 0.0, upper = 5.23, init = 0.523)
        """
        Lag time (hr)
        """
        pop_lag  RealDomain(lower = 0.0, upper = 5.0, init = 0.1)

        # Inter-individual variability
        """
          - ΩCL: Clearance
        """
        pk_Ω  PDiagDomain([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
        Vc = FSZV * pop_Vc * exp(pk_η[1])

        tabs = pop_tabs
        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 ~ @. Normal(cp, sqrt((σ_prop * cp)^2 + σ_add^2))
    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: 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.

Also note, that parameter inference can be expensive and for that reason we simplified the model for this tutorial to decrease overall runtime.

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

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.12,
    pop_Vc = 7.3,
    pop_tabs = 0.523,
    pop_lag = 0.5,
    pk_Ω = Diagonal([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     2.754619e+04     3.345035e+04
 * time: 0.024158000946044922
     1     5.618146e+03     1.386157e+04
 * time: 0.9305980205535889
     2     2.806727e+03     6.164132e+03
 * time: 0.9508500099182129
     3     1.814316e+03     3.717848e+03
 * time: 0.9698288440704346
     4     9.458168e+02     1.261428e+03
 * time: 0.9876370429992676
     5     6.506311e+02     6.846058e+02
 * time: 1.0048530101776123
     6     5.067237e+02     3.407504e+02
 * time: 1.0217719078063965
     7     4.543150e+02     1.615765e+02
 * time: 1.0379738807678223
     8     4.390128e+02     1.001202e+02
 * time: 1.0538949966430664
     9     4.359539e+02     1.004875e+02
 * time: 1.0689358711242676
    10     4.350899e+02     9.448070e+01
 * time: 1.083570957183838
    11     4.340047e+02     8.066223e+01
 * time: 1.098559856414795
    12     4.322893e+02     5.164762e+01
 * time: 1.1137919425964355
    13     4.307161e+02     2.264796e+01
 * time: 1.1290059089660645
    14     4.302010e+02     1.080854e+01
 * time: 1.143723964691162
    15     4.301247e+02     1.195746e+01
 * time: 1.1579129695892334
    16     4.301086e+02     9.445190e+00
 * time: 1.1723220348358154
    17     4.300988e+02     7.960007e+00
 * time: 1.2312710285186768
    18     4.300587e+02     7.716936e+00
 * time: 1.2456250190734863
    19     4.300051e+02     7.117247e+00
 * time: 1.2601749897003174
    20     4.299453e+02     7.208415e+00
 * time: 1.2747459411621094
    21     4.299169e+02     7.127430e+00
 * time: 1.2890698909759521
    22     4.299077e+02     7.022170e+00
 * time: 1.3030400276184082
    23     4.298994e+02     7.386201e+00
 * time: 1.3175418376922607
    24     4.298766e+02     8.313301e+00
 * time: 1.3312029838562012
    25     4.298227e+02     9.481826e+00
 * time: 1.3450238704681396
    26     4.296924e+02     1.085355e+01
 * time: 1.3595008850097656
    27     4.294405e+02     1.706433e+01
 * time: 1.37431001663208
    28     4.291427e+02     1.848276e+01
 * time: 1.3894219398498535
    29     4.289303e+02     1.058023e+01
 * time: 1.4050629138946533
    30     4.288371e+02     9.815220e+00
 * time: 1.420529842376709
    31     4.288319e+02     9.597527e+00
 * time: 1.4349720478057861
    32     4.288091e+02     8.868287e+00
 * time: 1.4499528408050537
    33     4.287649e+02     7.800905e+00
 * time: 1.4650509357452393
    34     4.286492e+02     8.240185e+00
 * time: 1.4800300598144531
    35     4.284193e+02     1.118489e+01
 * time: 1.4956419467926025
    36     4.280494e+02     1.142225e+01
 * time: 1.5110909938812256
    37     4.277503e+02     6.478284e+00
 * time: 1.5263378620147705
    38     4.276843e+02     3.593549e+00
 * time: 1.5418128967285156
    39     4.276750e+02     2.901808e+00
 * time: 1.5566258430480957
    40     4.276737e+02     2.925309e+00
 * time: 1.570786952972412
    41     4.276724e+02     2.942991e+00
 * time: 1.583894968032837
    42     4.276678e+02     2.975008e+00
 * time: 1.597606897354126
    43     4.276570e+02     3.010425e+00
 * time: 1.6108229160308838
    44     4.276278e+02     3.489417e+00
 * time: 1.6250059604644775
    45     4.275555e+02     4.423545e+00
 * time: 1.6657419204711914
    46     4.273850e+02     5.582714e+00
 * time: 1.6805288791656494
    47     4.270496e+02     6.557308e+00
 * time: 1.6958160400390625
    48     4.266090e+02     6.537839e+00
 * time: 1.7161118984222412
    49     4.263506e+02     3.908843e+00
 * time: 1.7330639362335205
    50     4.262968e+02     2.885389e+00
 * time: 1.7507109642028809
    51     4.262907e+02     2.842439e+00
 * time: 1.7680978775024414
    52     4.262900e+02     2.821199e+00
 * time: 1.784026861190796
    53     4.262894e+02     2.809832e+00
 * time: 1.799576997756958
    54     4.262879e+02     2.794691e+00
 * time: 1.8161790370941162
    55     4.262842e+02     2.781348e+00
 * time: 1.8322770595550537
    56     4.262743e+02     3.009687e+00
 * time: 1.848479986190796
    57     4.262477e+02     5.430536e+00
 * time: 1.8648090362548828
    58     4.261705e+02     9.974222e+00
 * time: 1.8814430236816406
    59     4.258983e+02     2.055347e+01
 * time: 1.8988220691680908
    60     4.256807e+02     2.645976e+01
 * time: 1.9267208576202393
    61     4.254650e+02     3.195354e+01
 * time: 1.9541189670562744
    62     4.250777e+02     4.156365e+01
 * time: 2.037363052368164
    63     4.250303e+02     5.260486e+01
 * time: 2.0557448863983154
    64     4.240447e+02     4.793859e+01
 * time: 2.073246955871582
    65     4.222991e+02     3.988756e+01
 * time: 2.1197099685668945
    66     4.175763e+02     1.012857e+01
 * time: 2.138514995574951
    67     4.156021e+02     1.547767e+01
 * time: 2.1644489765167236
    68     4.143567e+02     2.010443e+01
 * time: 2.1912620067596436
    69     4.130607e+02     1.768446e+01
 * time: 2.208889961242676
    70     4.116886e+02     8.431093e+00
 * time: 2.227548837661743
    71     4.112484e+02     7.131364e+00
 * time: 2.246048927307129
    72     4.109556e+02     4.884649e+00
 * time: 2.2643978595733643
    73     4.106945e+02     2.263126e+00
 * time: 2.2823939323425293
    74     4.105934e+02     3.304597e+00
 * time: 2.2998549938201904
    75     4.105660e+02     9.482860e-01
 * time: 2.317352056503296
    76     4.105529e+02     3.362095e-01
 * time: 2.334207057952881
    77     4.105426e+02     3.015047e-01
 * time: 2.3513100147247314
    78     4.105395e+02     1.804386e-01
 * time: 2.367568016052246
    79     4.105392e+02     3.818348e-02
 * time: 2.383018970489502
    80     4.105391e+02     4.591111e-02
 * time: 2.398068904876709
    81     4.105391e+02     8.641460e-03
 * time: 2.412148952484131
    82     4.105391e+02     2.258521e-03
 * time: 2.4255380630493164
    83     4.105391e+02     2.645649e-04
 * time: 2.437397003173828
FittedPumasModel

Successful minimization:                      true

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

Log-likelihood value:                   -410.53908
Number of subjects:                             31
Number of parameters:         Fixed      Optimized
                                  0              7
Observation records:         Active        Missing
    conc:                       239             47
    Total:                      239             47

-----------------------
             Estimate
-----------------------
pop_CL        0.12758
pop_Vc        8.3876
pop_tabs      0.48983
pop_lag       0.72171
pk_Ω₁,₁       0.010524
σ_prop        0.19957
σ_add         0.55153
-----------------------

2.5 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; 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:                   -410.53908
Number of subjects:                             31
Number of parameters:         Fixed      Optimized
                                  0              7
Observation records:         Active        Missing
    conc:                       239             47
    Total:                      239             47

---------------------------------------------------------------------
            Estimate           SE                      95.0% C.I.
---------------------------------------------------------------------
pop_CL       0.12758         0.0046779        [ 0.11842  ; 0.13675 ]
pop_Vc       8.3876          0.22464          [ 7.9474   ; 8.8279  ]
pop_tabs     0.48983         0.17922          [ 0.13856  ; 0.8411  ]
pop_lag      0.72171         0.16464          [ 0.39902  ; 1.0444  ]
pk_Ω₁,₁      0.010524        0.0071909        [-0.0035698; 0.024618]
σ_prop       0.19957         0.044554         [ 0.11225  ; 0.28689 ]
σ_add        0.55153         0.2842           [-0.0054881; 1.1085  ]
---------------------------------------------------------------------

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)
7×7 Symmetric{Float64, Matrix{Float64}}:
  2.1883e-5     0.000332493   5.45952e-6   …   3.30925e-5    0.000379372
  0.000332493   0.050461      0.00233924       0.00117512   -0.00391389
  5.45952e-6    0.00233924    0.032121        -0.000413843   0.00462368
  0.00019032   -0.0071644    -0.0160153       -0.00259753    0.0226855
 -1.83205e-6   -9.94307e-5   -0.000353466     -0.000234415   0.00119
  3.30925e-5    0.00117512   -0.000413843  …   0.00198502   -0.0102392
  0.000379372  -0.00391389    0.00462368      -0.0102392     0.0807681

and to get the condition number of the correlation matrix implied by vcov use

cond(inference_results)
44.47139507390104

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

It is also possible to calculate the correlation matrix from the vcov output using the cov2cor function from StatsBase

cor_from_cov = cov2cor(Matrix(vcov(inference_results)))
7×7 Matrix{Float64}:
  1.0          0.31641     0.00651188  …  -0.0544625   0.158779    0.285359
  0.31641      1.0         0.0581035      -0.0615543   0.117415   -0.0613072
  0.00651188   0.0581035   1.0            -0.274263   -0.0518272   0.0907765
  0.247112    -0.193717   -0.542756        0.550551   -0.354113    0.484835
 -0.0544625   -0.0615543  -0.274263        1.0        -0.731675    0.582294
  0.158779     0.117415   -0.0518272   …  -0.731675    1.0        -0.808653
  0.285359    -0.0613072   0.0907765       0.582294   -0.808653    1.0

And we see that the cond call above matches the condition number of the correlation matrix

cond(cor_from_cov)
44.47139507390104

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

  1. There are parameters very close to a bound (often 0)
  2. The parameter vector does not represent a local minimum (optimization failed)
  3. 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.

2.5.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(;
    rng = Random.default_rng,
    samples = 200,
    stratify_by = nothing,
    ensemblealg = EnsembleThreads(),
)

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.

bootstrap_results = infer(foce_fit, Bootstrap(samples = 50); level = 0.95)
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/ZhuZN/src/multivariate/optimize/optimize.jl:98
┌ Info: Bootstrap inference finished.
│   Total resampled fits = 50
│   Success rate = 1.0
└   Unique resampled populations = 50
Bootstrap inference results

Successful minimization:                      true

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

Log-likelihood value:                   -410.53908
Number of subjects:                             31
Number of parameters:         Fixed      Optimized
                                  0              7
Observation records:         Active        Missing
    conc:                       239             47
    Total:                      239             47

---------------------------------------------------------------------
            Estimate           SE                      95.0% C.I.
---------------------------------------------------------------------
pop_CL       0.12758         0.0091084       [0.11329    ;  0.13746]
pop_Vc       8.3876          1.9246          [3.9022     ; 10.166  ]
pop_tabs     0.48983         0.33162         [1.4535e-6  ;  0.95859]
pop_lag      0.72171         0.21772         [0.48501    ;  1.3232 ]
pk_Ω₁,₁      0.010524    11845.0             [4.8222e-37 ; 14.277  ]
σ_prop       0.19957         0.2432          [0.088629   ;  0.57126]
σ_add        0.55153         0.33856         [2.8645e-110;  1.0101 ]
---------------------------------------------------------------------
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)
7×7 Matrix{Float64}:
   8.29633e-5      -0.00168494    -9.86676e-5  …   4.96516e-5     0.000103152
  -0.00168494       3.70415        0.0272816       0.0523728      0.0200379
  -9.86676e-5       0.0272816      0.109974       -0.0237614      0.0578509
   0.000121945     -0.09738       -0.0388545       0.00924125    -0.0106758
 -10.4906       16547.1         -329.951          64.4109      -846.206
   4.96516e-5       0.0523728     -0.0237614   …   0.059145      -0.033538
   0.000103152      0.0200379      0.0578509      -0.033538       0.11462

and we can even get a DataFrame that includes all the estimated parameters from the sampled population fits

DataFrame(bootstrap_results.vcov)
50×7 DataFrame
25 rows omitted
Row pop_CL pop_Vc pop_tabs pop_lag pk_Ω₁,₁ σ_prop σ_add
Float64 Float64 Float64 Float64 Float64 Float64 Float64
1 0.134016 8.65287 0.484542 0.891717 0.0200174 0.162405 0.633772
2 0.130861 8.20114 0.254805 0.939278 0.00622446 0.208666 0.56751
3 0.125724 8.30701 0.393502 0.825934 0.00903701 0.178083 0.628499
4 0.129568 8.79904 0.389889 0.715993 0.0320396 0.112779 0.793014
5 0.129804 8.65905 0.385626 0.628387 0.0113616 0.197784 0.535303
6 0.119671 8.07233 0.415196 0.898305 0.0209074 0.167609 0.610805
7 0.132512 1.97477 2.4284e-5 0.999993 15.2185 0.276848 2.00509e-84
8 0.125363 8.51168 0.812658 0.48172 1.94724e-7 0.189389 0.628456
9 0.122847 8.35063 8.58455e-7 0.999999 2.14318e-36 0.235924 2.55822e-82
10 0.121942 3.24819 0.031055 0.98655 0.768065 0.219011 9.30696e-87
11 0.12435 7.66617 0.292211 0.826527 0.0174121 0.106904 0.686635
12 0.126239 8.10416 0.909928 0.810124 3.97662e-11 0.123419 1.08996
13 0.119051 8.28455 0.609158 0.505811 0.0488787 0.0998134 0.80045
39 0.120257 18.0464 0.225274 0.500656 83756.9 0.252149 2.2317e-91
40 0.127118 8.73291 0.57322 0.590384 0.00776682 0.187224 0.499899
41 0.123256 9.05071 0.656273 0.496347 0.0239417 0.164002 0.5845
42 0.137551 8.43927 0.511881 0.790102 0.0146792 0.129186 1.03698
43 0.134054 8.28148 0.332856 0.803674 0.0186361 0.189244 0.448263
44 0.130228 8.4935 0.679309 0.624027 0.013717 0.161494 0.757152
45 0.125817 6.45824 0.00206356 0.995432 0.000203833 0.209588 1.2731e-109
46 0.132371 8.37512 0.533469 0.74652 0.0119179 0.211287 0.559197
47 0.126551 8.14899 0.533174 0.821946 0.00317238 0.18244 0.469497
48 0.119491 8.3544 0.624988 0.623705 0.0201748 0.100483 0.792862
49 0.175267 8.88947 2.73607e-7 0.99936 8.74275e-5 0.65674 2.25566e-97
50 0.124357 8.67008 3.50293e-6 0.999998 4.3041e-248 0.214766 5.05989e-85

This is very useful for histogram plotting of parameter distributions.

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

sir_results = infer(foce_fit, SIR(samples = 1000, resamples = 200); level = 0.95)
[ Info: Calculating: variance-covariance matrix.
[ Info: Done.
[ Info: Running SIR.
[ Info: Resampling.
Simulated inference results

Successful minimization:                      true

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

Log-likelihood value:                   -410.53908
Number of subjects:                             31
Number of parameters:         Fixed      Optimized
                                  0              7
Observation records:         Active        Missing
    conc:                       239             47
    Total:                      239             47

---------------------------------------------------------------------
            Estimate           SE                      95.0% C.I.
---------------------------------------------------------------------
pop_CL       0.12758         0.0027897         [0.12254  ; 0.13311 ]
pop_Vc       8.3876          0.23385           [7.9652   ; 8.8762  ]
pop_tabs     0.48983         0.1317            [0.2235   ; 0.7372  ]
pop_lag      0.72171         0.099031          [0.50517  ; 0.85837 ]
pk_Ω₁,₁      0.010524        0.0047239         [0.0033439; 0.020362]
σ_prop       0.19957         0.028544          [0.14064  ; 0.2582  ]
σ_add        0.55153         0.21143           [0.060896 ; 0.89448 ]
---------------------------------------------------------------------

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)
7×7 Matrix{Float64}:
  7.78242e-6    0.000131615   4.20902e-5  …   1.88317e-6    0.000107411
  0.000131615   0.0546839     0.00167402      0.000498246  -0.00423702
  4.20902e-5    0.00167402    0.0173439      -0.00122985    0.0104688
  1.2793e-5    -0.00388864   -0.00856774     -0.000243919   0.00150448
 -6.24208e-7    4.99346e-5    5.84688e-5     -6.15249e-5    0.000348246
  1.88317e-6    0.000498246  -0.00122985  …   0.000814775  -0.00514636
  0.000107411  -0.00423702    0.0104688      -0.00514636    0.0447009

and

DataFrame(sir_results.vcov)
200×7 DataFrame
175 rows omitted
Row pop_CL pop_Vc pop_tabs pop_lag pk_Ω₁,₁ σ_prop σ_add
Float64 Float64 Float64 Float64 Float64 Float64 Float64
1 0.118713 8.36464 0.465706 0.585337 0.00454954 0.270517 0.00116172
2 0.126531 8.92586 0.462745 0.615701 0.00736712 0.234534 0.12411
3 0.128235 8.60987 0.193538 0.812456 0.0188321 0.243148 0.0272635
4 0.127477 8.39428 0.223876 0.856864 0.00334827 0.231972 0.135339
5 0.127826 8.31589 0.412528 0.625932 0.0178256 0.186454 0.484356
6 0.131338 9.12558 0.562954 0.736443 0.0124856 0.22665 0.545595
7 0.131355 8.27218 0.680454 0.616491 0.0106911 0.195922 0.655791
8 0.130895 8.51689 0.574083 0.458172 0.0163588 0.241413 0.435062
9 0.12912 8.3895 0.54489 0.679992 0.0155591 0.204949 0.500991
10 0.128168 8.57043 0.548686 0.622847 0.0152219 0.20032 0.666749
11 0.123633 8.87755 0.203274 0.806741 0.00421935 0.219183 0.227295
12 0.125555 8.34988 0.208634 0.692567 0.0117479 0.262446 0.0154857
13 0.126851 8.06447 0.455413 0.514951 0.00374792 0.260789 0.0308789
189 0.125441 8.55567 0.286595 0.822721 0.00979064 0.190246 0.428465
190 0.127533 8.28288 0.4076 0.810833 0.0123393 0.219448 0.338617
191 0.129956 8.64986 0.424785 0.758608 0.0114686 0.186786 0.688531
192 0.131063 8.22629 0.582104 0.705057 0.00971109 0.17166 0.777585
193 0.128001 8.22615 0.705439 0.606682 0.0062057 0.187526 0.68666
194 0.126105 8.1764 0.275277 0.847257 0.0131304 0.234261 0.394563
195 0.1274 8.64415 0.537182 0.693211 0.0038667 0.216172 0.552178
196 0.12978 8.55739 0.590816 0.672415 0.00474161 0.234714 0.340018
197 0.126596 8.24084 0.429096 0.776939 0.0195613 0.174611 0.544439
198 0.131805 8.2133 0.383215 0.816528 0.0116719 0.190063 0.663302
199 0.133768 8.65198 0.416143 0.759388 0.0173328 0.215892 0.636105
200 0.128892 8.12137 0.573466 0.69629 0.00300946 0.17674 0.624796

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

inference_results = infer(foce_fit, MarginalMCMC(); level = 0.95)

As sampling based inference can be computationally intensive we exclude the actual invocation of this method from this tutorial.

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