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 uncertainty based on the asymptotic variance-covariance formulas (robust or not).
  7. (Optionally) proceeding with more advanced techniques like bootstrapping or SIR for robust uncertainty quantification.

In previous tutorials, we already set up the data and performed a fit. We also obtained some parameter uncertainty estimates. In this tutorial, we will go more into depth with parameter uncertainty calculations using different methods. 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
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, 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 Depots1Central1

    @derived begin
        """
        Concentration (ng/mL)
        """
        conc ~ @. CombinedNormal(cp, σ_add, σ_prop)
    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.

warfarin_data = dataset("pumas/warfarin_pumas")

# Transform the data in a single chain of operations
warfarin_data_scales = @chain warfarin_data begin
    @rtransform begin
        # Scaling factors
        :FSZV = :wtbl / 70            # volume scaling
        :FSZCL = (:wtbl / 70)^0.75     # clearance scaling (allometric)
    end
end
330×12 DataFrame
305 rows omitted
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:

pop_pk = read_pumas(
    warfarin_data_scales;
    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
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.

3.1 Obtaining fit results

Following the examples in previous tutorials, we perform a fit. We need the output of the fit function call to perform inference to obtain parameter uncertainty estimates.

# 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_pk, 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.209064e+04     1.489225e+04
 * time: 0.03820610046386719
     1     2.643772e+03     3.167516e+03
 * time: 2.7492880821228027
     2     1.836601e+03     2.118430e+03
 * time: 2.773488998413086
     3     9.351337e+02     8.722439e+02
 * time: 2.7970900535583496
     4     6.402300e+02     4.199225e+02
 * time: 2.8161370754241943
     5     5.103664e+02     1.642121e+02
 * time: 2.834259033203125
     6     4.760464e+02     5.453749e+01
 * time: 2.852337121963501
     7     4.703757e+02     3.643518e+01
 * time: 2.8720290660858154
     8     4.699019e+02     3.135992e+01
 * time: 2.893193006515503
     9     4.697614e+02     2.953531e+01
 * time: 3.0455751419067383
    10     4.693153e+02     2.463233e+01
 * time: 3.0631890296936035
    11     4.685743e+02     2.580427e+01
 * time: 3.0804741382598877
    12     4.675133e+02     3.864937e+01
 * time: 3.0978469848632812
    13     4.666775e+02     5.495470e+01
 * time: 3.114856004714966
    14     4.661197e+02     5.692101e+01
 * time: 3.1317479610443115
    15     4.656782e+02     4.770992e+01
 * time: 3.148664951324463
    16     4.651802e+02     3.087698e+01
 * time: 3.165487051010132
    17     4.645523e+02     1.184834e+01
 * time: 3.182361125946045
    18     4.641447e+02     1.162249e+01
 * time: 3.1993019580841064
    19     4.639978e+02     1.125144e+01
 * time: 3.2171480655670166
    20     4.639307e+02     1.156463e+01
 * time: 3.2367141246795654
    21     4.638001e+02     1.312870e+01
 * time: 3.258301019668579
    22     4.635282e+02     1.480920e+01
 * time: 3.278942108154297
    23     4.630353e+02     2.169377e+01
 * time: 3.364474058151245
    24     4.623847e+02     4.478029e+01
 * time: 3.38191294670105
    25     4.617426e+02     6.468975e+01
 * time: 3.3993849754333496
    26     4.610293e+02     7.776996e+01
 * time: 3.4165079593658447
    27     4.597628e+02     8.785260e+01
 * time: 3.434009075164795
    28     4.566753e+02     9.769803e+01
 * time: 3.4521079063415527
    29     4.490421e+02     1.008838e+02
 * time: 3.4710381031036377
    30     4.391868e+02     9.978816e+01
 * time: 3.490103006362915
    31     4.130704e+02     5.917685e+01
 * time: 3.5099220275878906
    32     4.055780e+02     3.852824e+01
 * time: 3.528865098953247
    33     4.023118e+02     3.889618e+01
 * time: 3.5485939979553223
    34     4.012516e+02     3.694778e+01
 * time: 3.5887320041656494
    35     4.004391e+02     2.061948e+01
 * time: 3.6309640407562256
    36     3.983040e+02     3.508423e+01
 * time: 3.650810956954956
    37     3.969705e+02     3.841039e+01
 * time: 3.6711740493774414
    38     3.965462e+02     3.738343e+01
 * time: 3.6905710697174072
    39     3.950409e+02     3.064789e+01
 * time: 3.7112059593200684
    40     3.945750e+02     2.876429e+01
 * time: 3.7304770946502686
    41     3.937725e+02     2.571438e+01
 * time: 3.7512009143829346
    42     3.933955e+02     2.436112e+01
 * time: 3.7707459926605225
    43     3.927564e+02     2.051069e+01
 * time: 3.788702964782715
    44     3.916020e+02     1.629035e+01
 * time: 3.807025909423828
    45     3.886991e+02     2.689824e+01
 * time: 3.842358112335205
    46     3.870054e+02     2.298582e+01
 * time: 3.8615570068359375
    47     3.853691e+02     2.614992e+01
 * time: 3.8810341358184814
    48     3.841730e+02     2.207557e+01
 * time: 3.900136947631836
    49     3.825113e+02     2.204399e+01
 * time: 3.918447971343994
    50     3.808880e+02     2.444784e+01
 * time: 3.9377589225769043
    51     3.800407e+02     1.250611e+01
 * time: 3.9571239948272705
    52     3.798092e+02     1.167926e+01
 * time: 3.9763410091400146
    53     3.797789e+02     1.162382e+01
 * time: 4.010581016540527
    54     3.797069e+02     1.152441e+01
 * time: 4.028964042663574
    55     3.794424e+02     1.132717e+01
 * time: 4.048362970352173
    56     3.788131e+02     2.006438e+01
 * time: 4.067728042602539
    57     3.771525e+02     3.584695e+01
 * time: 4.087167978286743
    58     3.731299e+02     5.697249e+01
 * time: 4.1069841384887695
    59     3.658671e+02     6.542042e+01
 * time: 4.127027988433838
    60     3.604194e+02     4.036489e+01
 * time: 4.160676956176758
    61     3.532841e+02     1.574006e+01
 * time: 4.180357933044434
    62     3.520181e+02     1.393300e+01
 * time: 4.199341058731079
    63     3.517984e+02     6.701188e+00
 * time: 4.218188047409058
    64     3.517541e+02     3.503978e+00
 * time: 4.237472057342529
    65     3.516436e+02     8.720957e+00
 * time: 4.2715160846710205
    66     3.511845e+02     1.406200e+01
 * time: 4.29706597328186
    67     3.510647e+02     2.540378e+00
 * time: 4.315419912338257
    68     3.510209e+02     3.157201e+00
 * time: 4.333264112472534
    69     3.509959e+02     3.045642e+00
 * time: 4.350785970687866
    70     3.509765e+02     2.673143e+00
 * time: 4.379575967788696
    71     3.509751e+02     2.603975e+00
 * time: 4.397067070007324
    72     3.509724e+02     2.505719e+00
 * time: 4.414217948913574
    73     3.509666e+02     2.379768e+00
 * time: 4.431554079055786
    74     3.509504e+02     3.572030e+00
 * time: 4.448642015457153
    75     3.509123e+02     6.006350e+00
 * time: 4.476531028747559
    76     3.508288e+02     8.822995e+00
 * time: 4.494671106338501
    77     3.506944e+02     9.708012e+00
 * time: 4.512567043304443
    78     3.505767e+02     6.092631e+00
 * time: 4.530461072921753
    79     3.505358e+02     1.734431e+00
 * time: 4.548381090164185
    80     3.505314e+02     6.749379e-01
 * time: 4.576959133148193
    81     3.505313e+02     6.721982e-01
 * time: 4.5938990116119385
    82     3.505312e+02     6.699487e-01
 * time: 4.610193967819214
    83     3.505307e+02     6.606824e-01
 * time: 4.627099990844727
    84     3.505298e+02     6.413909e-01
 * time: 4.6441969871521
    85     3.505274e+02     9.083363e-01
 * time: 4.672111988067627
    86     3.505222e+02     1.339147e+00
 * time: 4.690190076828003
    87     3.505129e+02     1.608661e+00
 * time: 4.7081239223480225
    88     3.505026e+02     1.293164e+00
 * time: 4.725723028182983
    89     3.504973e+02     5.140504e-01
 * time: 4.743263006210327
    90     3.504963e+02     6.340189e-02
 * time: 4.772479057312012
    91     3.504963e+02     3.137914e-03
 * time: 4.789734125137329
    92     3.504963e+02     5.681551e-04
 * time: 4.805747032165527

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(
    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 also known as the robust variance-covariance 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 calculated using finite difference derivatives of the gradient calculated using automatic differentiation.

An example usage:

inference_results = infer(foce_fit; level = 0.95)
[ 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.116236834944694

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

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

cor_from_cov = cov2cor(vcov(inference_results))
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.116236834945

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:

  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.

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(;
    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/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
┌ Warning: Terminated early due to NaN in gradient.
└ @ Optim ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Optim/mv9zc/src/multivariate/optimize/optimize.jl:136
[ Info: Bootstrap inference finished.

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.018171    [ 0.11966  ; 0.1498  ]
pop_Vc     8.0535     0.26286     [ 7.7221   ; 8.5931  ]
pop_tabs   0.55061    0.18111     [ 0.17059  ; 0.74686 ]
pop_lag    0.87158    0.11801     [ 0.59048  ; 0.99579 ]
pk_Ω₁,₁    0.070642   7.3982e7    [ 0.027809 ; 0.13895 ]
pk_Ω₂,₂    0.018302   0.0052297   [ 0.0093465; 0.025817]
pk_Ω₃,₃    0.91326    812.14      [ 0.23052  ; 3.4542  ]
σ_prop     0.090096   0.029047    [ 0.066047 ; 0.11104 ]
σ_add      0.39115    0.086755    [ 0.19536  ; 0.48012 ]
---------------------------------------------------------
Unique resampled populations: 50 out of 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}}:
   0.000330177   -0.00280934     0.00141662  …  -0.000403363    0.00120112
  -0.00280934     0.0690958     -0.0093858       0.00396455    -0.0108105
   0.00141662    -0.0093858      0.0327998      -0.00170107     0.00711319
  -0.00047939     0.000628082   -0.0106          0.000866542   -0.000830135
  -1.2041e6       1.2716e7      -5.24339e6       1.92546e6     -3.67426e6
   4.89836e-5    -0.000265036    8.73559e-5  …  -7.07644e-5     0.000164492
 -13.2178       139.573        -57.6131         21.1412       -40.3219
  -0.000403363    0.00396455    -0.00170107      0.000843706   -0.000897462
   0.00120112    -0.0108105      0.00711319     -0.000897462    0.00752647

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

DataFrame(bootstrap_results.vcov)
50×9 DataFrame
25 rows omitted
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.131319 7.97105 0.582471 0.871871 0.0654118 0.0154608 1.35437 0.0905466 0.432506
2 0.14124 8.16302 0.481671 0.895266 0.0857694 0.0251178 1.05073 0.0915899 0.436168
3 0.122711 8.22518 0.322153 0.974309 0.028851 0.0191817 1.86095 0.0940816 0.2567
4 0.140951 8.19232 0.683915 0.877617 0.0850955 0.015726 0.611018 0.0813591 0.377295
5 0.118769 7.94039 0.248041 0.88402 0.0330799 0.0126374 0.23487 0.069426 0.18106
6 0.124215 7.96966 0.259253 1.34455 0.0337108 0.019698 0.789186 0.0767908 0.274182
7 0.126638 8.09169 0.624474 0.858893 0.027506 0.017869 0.416183 0.0806154 0.253326
8 0.133532 7.9767 0.260299 0.90792 0.0690409 0.0158791 9.63311e-8 0.0720204 0.298322
9 0.13441 7.96385 0.514761 0.503782 0.105138 0.0258584 1.11628 0.0820345 0.312115
10 0.133249 8.06085 0.5699 0.853994 0.06081 0.016405 0.644981 0.0566908 0.31005
11 0.136863 8.01597 0.62057 0.856669 0.0960049 0.0240731 1.29848 0.103265 0.442277
12 0.136192 8.06985 0.723992 0.899982 0.0719125 0.0138075 1.0677 0.102343 0.471828
13 0.127349 8.12303 0.490345 0.869875 0.0315313 0.0099133 0.940576 0.1069 0.289426
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
39 0.146749 8.35753 0.67614 0.817728 0.0926397 0.0143825 0.769897 0.0895647 0.408862
40 0.137525 8.18238 0.403133 0.844472 0.0825998 0.0251334 0.582921 0.0651941 0.322687
41 0.136465 8.03314 0.284041 0.953259 0.0608299 0.0224049 1.4598 0.106112 0.349143
42 0.145939 8.1475 0.560135 0.815112 0.0925801 0.021901 1.0703 0.0829728 0.417988
43 0.0189348 9.27279 5.43263e-8 1.0 5.23131e8 8.92681e-9 5743.73 0.27279 2.9721e-83
44 0.123052 8.36788 0.455028 0.848899 0.0598354 0.0175145 1.02986 0.0768723 0.362932
45 0.130948 8.17343 0.543095 0.832897 0.0867702 0.0133435 0.547069 0.0875057 0.35223
46 0.129798 7.9 0.560727 0.944459 0.0373703 0.015551 0.762204 0.10551 0.369161
47 0.144096 8.21206 0.675196 0.890738 0.0780203 0.0176547 1.55166 0.111114 0.481656
48 0.12765 7.89345 0.493122 0.751222 0.0605918 0.0169541 0.748374 0.0738874 0.32278
49 0.135769 7.98325 0.537905 0.850689 0.0783345 0.0256738 1.06406 0.0853467 0.379724
50 0.125322 8.05624 0.663164 0.83901 0.0540173 0.0180795 0.841809 0.0898713 0.37739

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.

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

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.0057385   [ 0.12358  ; 0.14764 ]
pop_Vc     8.0535     0.21109     [ 7.6534   ; 8.4406  ]
pop_tabs   0.55061    0.14948     [ 0.26274  ; 0.8637  ]
pop_lag    0.87158    0.03964     [ 0.78326  ; 0.9408  ]
pk_Ω₁,₁    0.070642   0.016758    [ 0.043798 ; 0.11006 ]
pk_Ω₂,₂    0.018302   0.0053721   [ 0.0098975; 0.028853]
pk_Ω₃,₃    0.91326    0.29542     [ 0.47764  ; 1.5919  ]
σ_prop     0.090096   0.0074943   [ 0.076386 ; 0.10385 ]
σ_add      0.39115    0.03792     [ 0.33994  ; 0.48248 ]
---------------------------------------------------------

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.29306e-5    0.000352248  -7.95701e-5   …  -4.4585e-6     4.57126e-5
  0.000352248   0.0445579     0.00241686       0.000120633  -0.00294428
 -7.95701e-5    0.00241686    0.0223433        5.95397e-5    0.000635898
 -3.7593e-6    -0.000388664  -0.00228619       2.50357e-5    6.14858e-5
  5.28056e-5   -0.000202562   8.81238e-5      -3.40801e-5    0.000336331
  1.98446e-6    0.000331548  -0.000107792  …  -5.56025e-6   -2.55938e-6
 -2.23731e-5   -0.00120656   -0.0215181        0.000240516   0.0021967
 -4.4585e-6     0.000120633   5.95397e-5       5.61651e-5   -2.55616e-5
  4.57126e-5   -0.00294428    0.000635898     -2.55616e-5    0.00143791

and

DataFrame(sir_results.vcov)
200×9 DataFrame
175 rows omitted
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.13327 8.23871 0.293479 0.903606 0.0850852 0.0366954 1.32912 0.0862742 0.40325
2 0.136486 8.09315 0.29119 0.876149 0.0957617 0.0154954 1.4708 0.0972086 0.386221
3 0.141388 8.03441 0.597863 0.88173 0.0919856 0.0338386 0.810543 0.0809403 0.432556
4 0.131327 7.81122 0.307676 0.911364 0.0756018 0.024983 1.60084 0.100176 0.386168
5 0.145072 8.35343 0.666395 0.86537 0.0932398 0.0226547 0.512596 0.0855722 0.453321
6 0.131932 7.71444 0.364918 0.878752 0.0826624 0.0228079 1.09324 0.0925351 0.46702
7 0.134163 8.23599 0.557435 0.864556 0.0818525 0.0255342 0.626125 0.093631 0.417867
8 0.131234 7.92992 0.627997 0.861541 0.0632066 0.0220061 0.664695 0.0815402 0.380316
9 0.131944 8.28891 0.539065 0.900762 0.0563799 0.0231639 1.3193 0.0903693 0.38811
10 0.132837 7.75604 0.630664 0.899542 0.067858 0.00960606 0.651845 0.0927538 0.381406
11 0.140125 8.04356 0.580479 0.849856 0.097623 0.0234293 1.3118 0.0785681 0.441195
12 0.134239 8.12596 0.65032 0.848233 0.0735381 0.0257521 0.844522 0.094949 0.441152
13 0.126604 7.77325 0.544872 0.908284 0.0521211 0.0139695 1.44381 0.0936241 0.424924
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
189 0.129437 8.03758 0.622224 0.888079 0.0672195 0.0149515 1.11192 0.0864444 0.44112
190 0.130455 7.6828 0.620931 0.861799 0.0684046 0.0226199 1.01437 0.0843362 0.442657
191 0.141967 8.13256 0.502556 0.87208 0.0803494 0.024851 0.940759 0.087864 0.367553
192 0.13356 8.22674 0.702079 0.829688 0.0438016 0.0217966 0.378632 0.102704 0.299375
193 0.129523 7.71078 0.86102 0.731426 0.0690416 0.0108373 0.477005 0.0921341 0.380976
194 0.132373 8.08544 0.876484 0.785955 0.0939284 0.0144205 0.890306 0.0974734 0.454959
195 0.128896 8.10251 0.596206 0.864424 0.0467019 0.0227244 0.953614 0.0912747 0.384391
196 0.140927 8.1968 0.396393 0.914915 0.0571474 0.026926 1.62209 0.110174 0.395227
197 0.142251 7.9103 0.364357 0.902547 0.0850153 0.0112359 1.08105 0.09387 0.463089
198 0.13961 7.97986 0.66482 0.820203 0.100055 0.0183952 0.402179 0.0867695 0.441889
199 0.136167 8.19218 0.560142 0.842924 0.0732617 0.0169752 0.904549 0.079643 0.392207
200 0.148767 8.44018 0.531465 0.896339 0.0891934 0.0201906 0.708347 0.0964724 0.372812

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.

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.

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.