A typical workflow for fitting a Pumas model and deriving parameter precision typically involves:
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.
Below is the PK model, named warfarin_pk_model, defined in Pumas. This model contains:
pop_CL, pop_Vc, pop_tabs, pop_lagpk_Ωσ_prop,σ_addFSZCL and FSZVusing Pumas
using PharmaDatasets
using DataFramesMeta
using PumasUtilities
using DimensionalData: DimArray, X, Y
using Random
# Bootstrap and SIR draw random numbers, so we seed the random number generator to make the
# results of this tutorial reproducible.
Random.seed!(545)TaskLocalRNG()warfarin_pk_model = @model begin
@metadata begin
desc = "Warfarin 1-compartment PK model (PD removed)"
timeu = "hr"
end
@param begin
# PK parameters
"""
Clearance (L/hr)
"""
pop_CL ∈ RealDomain(lower = 0.0, init = 0.134)
"""
Central volume (L)
"""
pop_Vc ∈ RealDomain(lower = 0.0, init = 8.11)
"""
Absorption lag time (hr)
"""
pop_tabs ∈ RealDomain(lower = 0.0, init = 0.523)
"""
Lag time (hr)
"""
pop_lag ∈ RealDomain(lower = 0.0, init = 0.1)
# Inter-individual variability
"""
- ΩCL: Clearance
- ΩVc: Central volume
- Ωtabs: Absorption lag time
"""
pk_Ω ∈ PDiagDomain([0.01, 0.01, 0.01])
# Residual variability
"""
σ_prop: Proportional error
"""
σ_prop ∈ RealDomain(lower = 0.0, init = 0.00752)
"""
σ_add: Additive error
"""
σ_add ∈ RealDomain(lower = 0.0, init = 0.0661)
end
@random begin
pk_η ~ MvNormal(pk_Ω) # mean = 0, covariance = pk_Ω
end
@covariates begin
"""
FSZCL: Clearance scaling factor
"""
FSZCL
"""
FSZV: Volume scaling factor
"""
FSZV
end
@pre begin
CL = FSZCL * pop_CL * exp(pk_η[1])
Vc = FSZV * pop_Vc * exp(pk_η[2])
tabs = pop_tabs * exp(pk_η[3])
Ka = log(2) / tabs
end
@dosecontrol begin
lags = (Depot = pop_lag,)
end
@vars begin
cp := Central / Vc
end
@dynamics Depots1Central1
@derived begin
"""
Concentration (ng/mL)
"""
conc ~ @. CombinedNormal(cp, σ_add, σ_prop)
end
endPumasModel
Parameters: pop_CL, pop_Vc, pop_tabs, pop_lag, pk_Ω, σ_prop, σ_add
Random effects: pk_η
Covariates: FSZCL, FSZV
Dynamical system variables: Depot, Central
Dynamical system type: Closed form
Derived: conc
Observed: concThe 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.
Four concentrations in this dataset are recorded as exactly zero, 0.5 hours after the dose and hence before the absorption lag, where the model predicts a concentration of exactly zero as well. The likelihood of such an observation is unbounded as the additive standard deviation σ_add approaches zero, so we treat these values as missing: a concentration of exactly zero means that the drug was not quantifiable, not that it was absent. See the BQL case study for a more detailed discussion.
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
# Concentrations of exactly zero are below the limit of quantification
@rtransform @passmissing :conc = iszero(:conc) ? missing : :conc
end| Row | id | time | evid | amt | cmt | conc | pca | wtbl | age | sex | FSZV | FSZCL |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Int64 | Float64 | Int64 | Float64? | Int64? | Float64? | Float64? | Float64 | Int64 | String1 | Float64 | Float64 | |
| 1 | 1 | 0.0 | 1 | 100.0 | 1 | missing | missing | 66.7 | 50 | M | 0.952857 | 0.96443 |
| 2 | 1 | 0.5 | 0 | missing | missing | missing | 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 |
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: concThe 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.
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.030400e+04 1.181292e+04 * time: 0.041338205337524414 1 2.262448e+03 2.457021e+03 * time: 8.838338136672974 2 1.606175e+03 1.689206e+03 * time: 8.85036015510559 3 8.249479e+02 6.944903e+02 * time: 8.862088203430176 4 5.830633e+02 3.309739e+02 * time: 9.14536714553833 5 4.835343e+02 1.258919e+02 * time: 9.154090166091919 6 4.597452e+02 3.446181e+01 * time: 9.162438154220581 7 4.564004e+02 2.157414e+01 * time: 9.170516014099121 8 4.561766e+02 1.841633e+01 * time: 9.178284168243408 9 4.560933e+02 1.718302e+01 * time: 9.185716152191162 10 4.557723e+02 1.676719e+01 * time: 9.193552017211914 11 4.552137e+02 2.081976e+01 * time: 9.201372146606445 12 4.542913e+02 2.597069e+01 * time: 9.209287166595459 13 4.534692e+02 3.527716e+01 * time: 9.217392206192017 14 4.529383e+02 3.402880e+01 * time: 9.225499153137207 15 4.526116e+02 2.490636e+01 * time: 9.233367204666138 16 4.523614e+02 1.340594e+01 * time: 9.241147994995117 17 4.520146e+02 1.295774e+01 * time: 9.248842000961304 18 4.517070e+02 1.823845e+01 * time: 9.25629711151123 19 4.514823e+02 2.416090e+01 * time: 9.263714075088501 20 4.512324e+02 2.583177e+01 * time: 9.271071195602417 21 4.506901e+02 2.632776e+01 * time: 9.278501033782959 22 4.493304e+02 3.098682e+01 * time: 9.286169052124023 23 4.457731e+02 4.060488e+01 * time: 9.294275045394897 24 4.374540e+02 5.129320e+01 * time: 9.302879095077515 25 4.278696e+02 5.289924e+01 * time: 9.312424182891846 26 4.221945e+02 3.807050e+01 * time: 9.321943998336792 27 4.126997e+02 2.293280e+01 * time: 9.332703113555908 28 4.085329e+02 2.226988e+01 * time: 9.34449815750122 29 4.073530e+02 1.516843e+01 * time: 9.355534076690674 30 4.066443e+02 1.009682e+01 * time: 9.366791009902954 31 4.064325e+02 7.665676e+00 * time: 9.377676010131836 32 4.063169e+02 7.945666e+00 * time: 9.388749122619629 33 4.062598e+02 8.057281e+00 * time: 9.39911699295044 34 4.061917e+02 8.096740e+00 * time: 9.409893989562988 35 4.060744e+02 1.001791e+01 * time: 9.420825004577637 36 4.058146e+02 1.313770e+01 * time: 9.431816101074219 37 4.052323e+02 1.489612e+01 * time: 9.529898166656494 38 4.040793e+02 2.371902e+01 * time: 9.539597034454346 39 4.022553e+02 3.961614e+01 * time: 9.549413204193115 40 4.000227e+02 4.701290e+01 * time: 9.559117078781128 41 3.977591e+02 4.815312e+01 * time: 9.568570137023926 42 3.965739e+02 4.653027e+01 * time: 9.57756519317627 43 3.916661e+02 4.293389e+01 * time: 9.586708068847656 44 3.895667e+02 4.702428e+01 * time: 9.595414161682129 45 3.889000e+02 5.508619e+01 * time: 9.604875087738037 46 3.868670e+02 5.359374e+01 * time: 9.614046096801758 47 3.857187e+02 5.483765e+01 * time: 9.623219013214111 48 3.801077e+02 6.169015e+01 * time: 9.633714199066162 49 3.758913e+02 5.242755e+01 * time: 9.643076181411743 50 3.704872e+02 3.681565e+01 * time: 9.652875185012817 51 3.647521e+02 6.405310e+01 * time: 9.662696123123169 52 3.593658e+02 5.044209e+01 * time: 9.672236204147339 53 3.525544e+02 2.170001e+01 * time: 9.68156909942627 54 3.519477e+02 9.793970e+00 * time: 9.690324068069458 55 3.518024e+02 1.071185e+01 * time: 9.698784112930298 56 3.517872e+02 1.052811e+01 * time: 9.706431150436401 57 3.516872e+02 8.722351e+00 * time: 9.714680194854736 58 3.515998e+02 7.654274e+00 * time: 9.722624063491821 59 3.513181e+02 3.885265e+00 * time: 9.732073068618774 60 3.511057e+02 4.956890e+00 * time: 9.7415771484375 61 3.510219e+02 3.011426e+00 * time: 9.75108003616333 62 3.510122e+02 2.597058e+00 * time: 9.760331153869629 63 3.510101e+02 2.262243e+00 * time: 9.769317150115967 64 3.510069e+02 2.230923e+00 * time: 9.778213024139404 65 3.509975e+02 2.157133e+00 * time: 9.814602136611938 66 3.509746e+02 2.136095e+00 * time: 9.822688102722168 67 3.509195e+02 3.322024e+00 * time: 9.83090615272522 68 3.508068e+02 4.891317e+00 * time: 9.839011192321777 69 3.506398e+02 5.092565e+00 * time: 9.847247123718262 70 3.505135e+02 2.778868e+00 * time: 9.855555057525635 71 3.504732e+02 9.934168e-01 * time: 9.863656044006348 72 3.504659e+02 2.186379e-01 * time: 9.871395111083984 73 3.504655e+02 1.956674e-01 * time: 9.878870010375977 74 3.504655e+02 1.910188e-01 * time: 9.885939121246338 75 3.504654e+02 1.883529e-01 * time: 9.892736196517944 76 3.504654e+02 1.889961e-01 * time: 9.899737119674683 77 3.504653e+02 1.880938e-01 * time: 9.906795978546143 78 3.504651e+02 1.815117e-01 * time: 9.914002180099487 79 3.504646e+02 1.748458e-01 * time: 9.921251058578491 80 3.504636e+02 2.185751e-01 * time: 9.928873062133789 81 3.504626e+02 1.810451e-01 * time: 9.936617136001587 82 3.504620e+02 7.615336e-02 * time: 9.943977117538452 83 3.504619e+02 1.139943e-02 * time: 9.95112919807434 84 3.504619e+02 9.569155e-04 * time: 9.958011150360107
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 32
Observation records: Active Missing
conc: 247 51
Total: 247 51
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.46191
--------------------
Estimate
--------------------
pop_CL 0.13469
pop_Vc 8.048
pop_tabs 0.5449
pop_lag 0.87702
pk_Ω₁,₁ 0.070708
pk_Ω₂,₂ 0.018349
pk_Ω₃,₃ 0.92441
σ_prop 0.087933
σ_add 0.41375
--------------------inferThe 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.
In its simplest form infer takes the output of fit and returns the asymptotic variance-covariance based inference results. The keyword arguments that this tutorial relies on are:
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.The complete interface is documented in the Inference section of the Pumas documentation.
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: 247 51
Total: 247 51
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.46191
---------------------------------------------------------
Estimate SE 95.0% C.I.
---------------------------------------------------------
pop_CL 0.13469 0.0066655 [ 0.12163 ; 0.14776 ]
pop_Vc 8.048 0.2209 [ 7.6151 ; 8.481 ]
pop_tabs 0.5449 0.18347 [ 0.1853 ; 0.9045 ]
pop_lag 0.87702 0.052956 [ 0.77323 ; 0.98081 ]
pk_Ω₁,₁ 0.070708 0.024568 [ 0.022555 ; 0.11886 ]
pk_Ω₂,₂ 0.018349 0.0051547 [ 0.0082457; 0.028452]
pk_Ω₃,₃ 0.92441 0.40816 [ 0.12443 ; 1.7244 ]
σ_prop 0.087933 0.014414 [ 0.059682 ; 0.11618 ]
σ_add 0.41375 0.068273 [ 0.27994 ; 0.54756 ]
---------------------------------------------------------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.44286e-5 0.000209426 0.000291729 … 1.71964e-6 0.000214248
0.000209426 0.0487961 0.0052402 -0.000877066 -0.00559521
0.000291729 0.0052402 0.0336614 9.51391e-5 0.00507281
-6.46913e-5 -0.00187707 -0.00376579 0.000470405 0.000401524
0.000120852 4.07894e-5 0.00160924 -9.8414e-5 0.000775274
3.18546e-7 0.000284206 -0.000124859 … 5.47807e-6 8.6199e-6
-0.000257165 -0.0233434 -0.0263984 0.00340652 0.012124
1.71964e-6 -0.000877066 9.51391e-5 0.000207774 0.000517359
0.000214248 -0.00559521 0.00507281 0.000517359 0.00466121and to get the condition number of the correlation matrix implied by vcov use
cond(inference_results)47.213447836476696Some 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)12933.533857380817It 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.142235 0.238552 … -0.0945255 0.0178983 0.4708
0.142235 1.0 0.129297 -0.258905 -0.275451 -0.371001
0.238552 0.129297 1.0 -0.352516 0.0359747 0.404979
-0.183275 -0.160463 -0.387594 0.544998 0.616259 0.111058
0.737986 0.00751592 0.35701 -0.280957 -0.2779 0.462203
0.00927129 0.249597 -0.132024 … 0.30886 0.0737276 0.0244935
-0.0945255 -0.258905 -0.352516 1.0 0.579006 0.435075
0.0178983 -0.275451 0.0359747 0.579006 1.0 0.525711
0.4708 -0.371001 0.404979 0.435075 0.525711 1.0And we see that the cond call above matches the condition number of the correlation matrix
cond(cor_from_cov)47.213447836476526It is well-known that the asymptotic variance-covariance matrix can sometimes fail to compute. This can happen for a variety of reasons including:
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.
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.
To bootstrap, pass a Bootstrap specification as the second argument to infer. Bootstrap accepts the random number generator rng, the number samples of resampled datasets to produce, the covariates stratify_by to stratify the sampling by, and finally the ensemble algorithm ensemblealg to control parallelization across fits. On the JuliaHub platform this can be used together with distributed computing to perform many resampled estimations in parallel. See Bootstrap-based inference for the full interface.
bootstrap_results = infer(foce_fit, Bootstrap(samples = 50); level = 0.95)[ Info: Bootstrap inference finished.
Bootstrap inference results
Dynamical system type: Closed form
Number of subjects: 32
Observation records: Active Missing
conc: 247 51
Total: 247 51
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.46191
---------------------------------------------------------
Estimate SE 95.0% C.I.
---------------------------------------------------------
pop_CL 0.13469 0.00705 [ 0.1258 ; 0.15008 ]
pop_Vc 8.048 0.20482 [ 7.7298 ; 8.5011 ]
pop_tabs 0.5449 0.22249 [ 0.20383 ; 0.8891 ]
pop_lag 0.87702 0.11873 [ 0.71447 ; 1.2905 ]
pk_Ω₁,₁ 0.070708 0.026616 [ 0.027214 ; 0.13497 ]
pk_Ω₂,₂ 0.018349 0.0047874 [ 0.0091463; 0.025145]
pk_Ω₃,₃ 0.92441 0.57267 [ 0.34131 ; 2.2985 ]
σ_prop 0.087933 0.018767 [ 0.049607 ; 0.12258 ]
σ_add 0.41375 0.067692 [ 0.26804 ; 0.5059 ]
---------------------------------------------------------
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}}:
4.97025e-5 0.000196475 0.000229811 … 1.24279e-5 0.000206047
0.000196475 0.041951 0.0184019 -0.00135934 -0.00347977
0.000229811 0.0184019 0.0495035 -0.000464115 0.000868209
-0.000103609 -0.00830307 -0.00788933 0.000699514 0.00266545
0.00012788 0.00109018 0.00165868 -0.000213898 0.000736655
4.54419e-6 -0.000118695 -0.000252352 … -1.06194e-5 9.85166e-5
-0.0003823 -0.0385999 -0.0636734 0.00299463 0.0150342
1.24279e-5 -0.00135934 -0.000464115 0.000352201 0.000231977
0.000206047 -0.00347977 0.000868209 0.000231977 0.00458217and we can even get a DataFrame that includes all the estimated parameters from the sampled population fits
samples_table(bootstrap_results)| Row | pop_CL | pop_Vc | pop_tabs | pop_lag | pk_Ω₁,₁ | pk_Ω₂,₂ | pk_Ω₃,₃ | σ_prop | σ_add | loglikelihood |
|---|---|---|---|---|---|---|---|---|---|---|
| Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | 0.138199 | 7.88269 | 0.261311 | 0.853437 | 0.0652341 | 0.019051 | 0.839352 | 0.07037 | 0.388899 | -359.84 |
| 2 | 0.147409 | 8.06405 | 0.542726 | 0.898654 | 0.0766633 | 0.00943141 | 0.446809 | 0.0865969 | 0.40379 | -356.595 |
| 3 | 0.130045 | 8.32838 | 0.187929 | 0.937123 | 0.0327084 | 0.00915136 | 1.85535 | 0.0917191 | 0.303618 | -369.617 |
| 4 | 0.142971 | 7.99311 | 0.621799 | 0.87429 | 0.0821683 | 0.0207287 | 0.712888 | 0.089067 | 0.430751 | -351.717 |
| 5 | 0.136384 | 8.04222 | 0.659984 | 1.38333 | 0.0756629 | 0.0219668 | 1.18834 | 0.0893547 | 0.505606 | -444.244 |
| 6 | 0.135858 | 7.89788 | 0.208673 | 0.943898 | 0.0652191 | 0.0196336 | 1.92101 | 0.0932864 | 0.378394 | -354.771 |
| 7 | 0.131619 | 8.19803 | 0.375198 | 0.893224 | 0.0458887 | 0.0134883 | 1.09276 | 0.0810651 | 0.355713 | -356.315 |
| 8 | 0.145999 | 8.15935 | 0.663138 | 0.795883 | 0.10125 | 0.0164858 | 0.484316 | 0.0760988 | 0.387434 | -358.184 |
| 9 | 0.148771 | 8.29753 | 0.611214 | 0.860173 | 0.0609043 | 0.0179712 | 1.36196 | 0.0953851 | 0.422294 | -354.566 |
| 10 | 0.146446 | 8.13234 | 0.450113 | 0.83895 | 0.13755 | 0.0248324 | 0.5953 | 0.0733565 | 0.418513 | -358.32 |
| 11 | 0.133372 | 7.75314 | 0.693663 | 0.898928 | 0.0523328 | 0.0169064 | 0.918215 | 0.0938977 | 0.457898 | -353.912 |
| 12 | 0.143182 | 8.06781 | 0.837878 | 0.712793 | 0.0684965 | 0.0101765 | 0.389407 | 0.10978 | 0.296499 | -360.457 |
| 13 | 0.141068 | 8.27269 | 0.773788 | 0.851858 | 0.0878996 | 0.0135328 | 0.353973 | 0.0778773 | 0.384881 | -358.158 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| 39 | 0.145104 | 7.72301 | 0.490232 | 0.96761 | 0.0826272 | 0.0153741 | 2.31155 | 0.1321 | 0.549725 | -374.62 |
| 40 | 0.137929 | 8.50375 | 0.541095 | 0.849905 | 0.0738461 | 0.0175471 | 0.398782 | 0.0694119 | 0.314915 | -369.007 |
| 41 | 0.138132 | 8.04736 | 0.536048 | 0.820793 | 0.105426 | 0.0145463 | 1.05958 | 0.080938 | 0.451293 | -353.848 |
| 42 | 0.137203 | 8.08756 | 0.437777 | 0.873898 | 0.0694537 | 0.0216998 | 0.960834 | 0.0777146 | 0.41107 | -352.07 |
| 43 | 0.151397 | 8.01287 | 0.245036 | 0.911006 | 0.107094 | 0.0252359 | 1.61701 | 0.0878252 | 0.424717 | -356.191 |
| 44 | 0.141688 | 8.00417 | 0.202426 | 0.970769 | 0.0580152 | 0.0209218 | 2.25374 | 0.100388 | 0.419911 | -375.654 |
| 45 | 0.125658 | 7.89791 | 0.500011 | 0.947889 | 0.0273326 | 0.0214628 | 1.48732 | 0.106666 | 0.362837 | -367.292 |
| 46 | 0.139285 | 7.88893 | 0.379051 | 0.944136 | 0.0695366 | 0.0246278 | 1.1083 | 0.0927457 | 0.378043 | -354.556 |
| 47 | 0.131241 | 8.10441 | 0.482215 | 0.862209 | 0.0606528 | 0.0163825 | 0.458794 | 0.0798909 | 0.334757 | -357.276 |
| 48 | 0.135471 | 8.1099 | 0.485273 | 0.885359 | 0.0680653 | 0.0129145 | 0.915354 | 0.0493309 | 0.334917 | -390.426 |
| 49 | 0.12667 | 7.85383 | 0.339741 | 0.927379 | 0.0271793 | 0.0228594 | 1.07212 | 0.104092 | 0.331625 | -366.854 |
| 50 | 0.129295 | 8.10596 | 0.436212 | 0.872782 | 0.0501195 | 0.0127704 | 0.503576 | 0.0773306 | 0.221905 | -382.848 |
This is very useful for histogram plotting of parameter distributions.
The bootstrap estimates are meant to be distributed like the maximum likelihood estimator, and Dosne et al. (Dosne, Niebecker, et al. 2016) proposed a diagnostic for whether they are. The likelihood-ratio statistic of a bootstrap estimate, the difference in objective function value (dOFV) of the pharmacometrics literature, is twice the difference between the log-likelihood of the parameter estimates and the log-likelihood of that bootstrap estimate, both evaluated on the original dataset. It therefore measures how much worse the bootstrap estimate explains those data, and is not the objective function value that its own fit reached. Were the estimates distributed like the estimator, their statistics would asymptotically follow a χ² distribution with one degree of freedom per estimated parameter for an unconstrained fixed effects model. For a nonlinear mixed effects model the degrees of freedom can only be lower. The likelihood approximation of the estimation method and any misspecification of the model influence the statistics as well (Dosne, Niebecker, et al. 2016). likelihood_ratio_ppplot plots the empirical cumulative distribution of the statistics of the bootstrap estimates against the cumulative distribution function of the χ² reference.
likelihood_ratio_ppplot(bootstrap_results)
A curve above the identity line has smaller statistics than the reference and is spread more narrowly than it, and below the line it is spread more widely. The annotated degrees of freedom, obtained by fitting a χ² distribution to the statistics, express that comparison as a single number. Since the degrees of freedom cannot exceed the number of estimated parameters, a curve on or above the identity line is acceptable, whereas one below it describes the estimator badly, the estimates then including parameter vectors that explain the data worse than the reference allows for. Since the reference assumes that the fit converged to the maximum likelihood estimates, the diagnostic is only meaningful for a converged fit.
Pumas supports inference through SIR with the SIR specification. SIR accepts the random number generator rng, the number samples of parameter samples drawn from the proposal distribution, the number resamples of parameter samples that are then drawn from those by weighted sampling without replacement, and finally the scale proposal_scale and the shape proposal_shape of the space that a proposal distribution is fitted on. Neither samples nor resamples has a default value, so both have to be set. Dosne et al. (Dosne, Bergstrand, et al. 2016), who introduced SIR for nonlinear mixed effects models, recommend a samples/resamples ratio of 5 as a starting point so that enough samples are available to resample from; the ratio that turned out to be sufficient varied between 2 and 6 across their examples. Like the bootstrap confidence intervals, and unlike the asymptotic ones, the confidence intervals based on SIR respect the bounds of the model, so for instance the interval of a strictly positive parameter cannot contain negative values.
In each iteration, parameter samples are drawn from a multivariate normal proposal distribution whose mean is the parameter estimates of the fitted model. From the second iteration onwards, the proposal distribution is fitted to the parameter samples resampled in the previous iteration, with its covariance matrix set to their empirical covariance matrix. The proposal distribution of the initial iteration, on the other hand, is built from a covariance matrix that has to be specified. Such a matrix is expressed on the natural scale of the parameters, where a sample can violate the domain of a parameter and is then rejected, so that proposal distribution is a truncated multivariate normal. One choice is to set it to an estimate of the variance-covariance matrix of the parameters, which can be taken from existing inference results:
sir_results = infer(inference_results, SIR(samples = 1000, resamples = 200); level = 0.95)[ Info: Running SIR.
Simulated inference results
Dynamical system type: Closed form
Number of subjects: 32
Observation records: Active Missing
conc: 247 51
Total: 247 51
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.46191
--------------------------------------------------------
Estimate SE 95.0% C.I.
--------------------------------------------------------
pop_CL 0.13469 0.0060247 [ 0.12473 ; 0.14795 ]
pop_Vc 8.048 0.20123 [ 7.651 ; 8.3759 ]
pop_tabs 0.5449 0.13797 [ 0.27665 ; 0.795 ]
pop_lag 0.87702 0.032241 [ 0.81433 ; 0.93974 ]
pk_Ω₁,₁ 0.070708 0.016364 [ 0.044078; 0.10684 ]
pk_Ω₂,₂ 0.018349 0.0046353 [ 0.011472; 0.028963]
pk_Ω₃,₃ 0.92441 0.27962 [ 0.52279 ; 1.5776 ]
σ_prop 0.087933 0.0075669 [ 0.076266; 0.10429 ]
σ_add 0.41375 0.038512 [ 0.36037 ; 0.49908 ]
--------------------------------------------------------Since inference_results was computed with infer(foce_fit; level = 0.95) above, the initial proposal distribution here is based on the robust variance-covariance matrix, which is reused rather than recomputed.
The ensemblealg keyword argument controls the parallelization of the SIR procedure itself: it accepts EnsembleThreads() (the default value) to use multi-threading or EnsembleSerial() to use a single thread. Because the inference results are computed in a separate step, that step can use an ensemble algorithm that the SIR procedure does not support, such as EnsembleDistributed(). See Sampling importance resampling (SIR) for the full interface.
Because the initial proposal distribution above is based on the robust variance-covariance matrix, that SIR procedure fails whenever the matrix fails to be computed or is not positive definite. This is not inherent to SIR: any of the other inference results in this tutorial can serve as its basis instead. Passing the results of an infer call with sandwich_estimator = false bases it on the inverse of the Hessian matrix,
sir_from_hessian =
infer(infer(foce_fit; sandwich_estimator = false), SIR(samples = 1000, resamples = 200))[ 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: 247 51
Total: 247 51
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.46191
---------------------------------------------------------
Estimate SE 95.0% C.I.
---------------------------------------------------------
pop_CL 0.13469 0.0068339 [ 0.12127 ; 0.14765 ]
pop_Vc 8.048 0.24899 [ 7.633 ; 8.5803 ]
pop_tabs 0.5449 0.15659 [ 0.26393 ; 0.91255 ]
pop_lag 0.87702 0.038342 [ 0.79218 ; 0.93918 ]
pk_Ω₁,₁ 0.070708 0.018764 [ 0.043353 ; 0.11691 ]
pk_Ω₂,₂ 0.018349 0.0062882 [ 0.0099289; 0.033066]
pk_Ω₃,₃ 0.92441 0.36868 [ 0.47945 ; 1.8382 ]
σ_prop 0.087933 0.0098462 [ 0.070299 ; 0.11006 ]
σ_add 0.41375 0.047888 [ 0.33664 ; 0.52017 ]
---------------------------------------------------------and passing the results of a Bootstrap, SIR, or MarginalMCMC inference bases it on the covariance matrix of their parameter samples, which requires no matrix to be inverted at all:
sir_from_bootstrap = infer(bootstrap_results, SIR(samples = 1000, resamples = 200))[ Info: Running SIR.
Simulated inference results
Dynamical system type: Closed form
Number of subjects: 32
Observation records: Active Missing
conc: 247 51
Total: 247 51
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.46191
--------------------------------------------------------
Estimate SE 95.0% C.I.
--------------------------------------------------------
pop_CL 0.13469 0.0058976 [ 0.12603 ; 0.14767 ]
pop_Vc 8.048 0.18456 [ 7.6978 ; 8.3988 ]
pop_tabs 0.5449 0.18299 [ 0.23223 ; 0.9267 ]
pop_lag 0.87702 0.054062 [ 0.74943 ; 0.95286 ]
pk_Ω₁,₁ 0.070708 0.022496 [ 0.046921; 0.13495 ]
pk_Ω₂,₂ 0.018349 0.0054718 [ 0.011848; 0.034778]
pk_Ω₃,₃ 0.92441 0.5839 [ 0.48379 ; 2.8955 ]
σ_prop 0.087933 0.010739 [ 0.069692; 0.10951 ]
σ_add 0.41375 0.048408 [ 0.3414 ; 0.53533 ]
--------------------------------------------------------The covariance matrix can also be supplied directly with the proposal_covariance keyword argument. It has to be a DimArray whose rows and columns are labelled with the names of the estimated parameters, so that the matrix is unambiguous regardless of how Pumas vectorizes the parameters internally. These names are the parameter entries of the rows of coeftable for which constant is false, since constant parameters are not estimated:
ct = coeftable(foce_fit)
labels = ct.parameter[.!ct.constant]
sir_from_matrix = infer(
foce_fit,
SIR(samples = 1000, resamples = 200);
proposal_covariance = DimArray(vcov(bootstrap_results), (X(labels), Y(labels))),
)[ Info: Running SIR.
Simulated inference results
Dynamical system type: Closed form
Number of subjects: 32
Observation records: Active Missing
conc: 247 51
Total: 247 51
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.46191
---------------------------------------------------------
Estimate SE 95.0% C.I.
---------------------------------------------------------
pop_CL 0.13469 0.0061292 [ 0.1264 ; 0.15063 ]
pop_Vc 8.048 0.19751 [ 7.6649 ; 8.3943 ]
pop_tabs 0.5449 0.16568 [ 0.2332 ; 0.85807 ]
pop_lag 0.87702 0.052954 [ 0.76716 ; 0.95258 ]
pk_Ω₁,₁ 0.070708 0.020905 [ 0.042124 ; 0.116 ]
pk_Ω₂,₂ 0.018349 0.005012 [ 0.0099553; 0.028281]
pk_Ω₃,₃ 0.92441 0.41815 [ 0.40821 ; 1.9073 ]
σ_prop 0.087933 0.011591 [ 0.069468 ; 0.11815 ]
σ_add 0.41375 0.047599 [ 0.33622 ; 0.51876 ]
---------------------------------------------------------samples and resamples also accept vectors, in which case the SIR procedure is repeated for as many iterations as there are entries. Since every iteration but the first builds its proposal distribution from the parameter samples resampled in the previous one, iterating lets the procedure move away from an initial proposal distribution that approximated the parameter uncertainty poorly. The schedule below is the one Dosne et al. (Dosne et al. 2017) use for the first three iterations of their automated SIR procedure: a large samples/resamples ratio improves the proposal distribution quickly while the number of likelihood evaluations stays small.
sir_iterated =
infer(inference_results, SIR(samples = [1000, 1000, 1000], resamples = [200, 400, 500]))[ Info: Running SIR.
Simulated inference results
Dynamical system type: Closed form
Number of subjects: 32
Observation records: Active Missing
conc: 247 51
Total: 247 51
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -350.46191
--------------------------------------------------------
Estimate SE 95.0% C.I.
--------------------------------------------------------
pop_CL 0.13469 0.0054486 [ 0.12571 ; 0.1472 ]
pop_Vc 8.048 0.20905 [ 7.6207 ; 8.4479 ]
pop_tabs 0.5449 0.15053 [ 0.32467 ; 0.93324 ]
pop_lag 0.87702 0.03415 [ 0.81052 ; 0.94343 ]
pk_Ω₁,₁ 0.070708 0.01799 [ 0.04753 ; 0.11676 ]
pk_Ω₂,₂ 0.018349 0.0062482 [ 0.010781; 0.035804]
pk_Ω₃,₃ 0.92441 0.39666 [ 0.52473 ; 2.0959 ]
σ_prop 0.087933 0.0080058 [ 0.072476; 0.10375 ]
σ_add 0.41375 0.036754 [ 0.35459 ; 0.49593 ]
--------------------------------------------------------In every iteration but the last, resamples has to be greater than the number of estimated parameters, since otherwise the covariance matrix of the proposal distribution of the next iteration is singular.
Whenever a proposal distribution is fitted to parameter samples rather than built from a covariance matrix, the proposal_scale and proposal_shape keyword arguments select the space it is fitted on.
proposal_scale selects the scale of the parameters. Its default value UnconstrainedScale() is the unconstrained parameterization that Pumas uses internally, on which every sample drawn from the proposal distribution corresponds to a valid parameter value, so that no sample is ever rejected. The alternative NaturalScale() is the original domain of the parameters, which is also the scale that covariance matrices are expressed on, and there a sample can violate the domain of a parameter and is then rejected.
proposal_shape selects whether the marginals of the parameter samples are reshaped before the proposal distribution is fitted to them. Its default value nothing leaves them as they are, whereas YeoJohnson() applies a Yeo-Johnson transformation with one parameter per marginal that is estimated from the samples themselves, and rejects the samples that fall outside the image of the transformation. The proposal distribution is multivariate normal on the space it is fitted on, and the unconstrained scale can already reduce the skewness that the domain of a parameter induces, for instance by taking the logarithm of a strictly positive parameter. Reshaping the marginals addresses skewness that remains on the selected scale: SIR(samples = [1000, 1000, 1000], resamples = [200, 400, 500], proposal_shape = YeoJohnson()) runs the schedule above with the proposal distributions of its second and third iteration fitted in this way.
Both keyword arguments apply to the proposal distributions of the second and all subsequent iterations, and to the one of the initial iteration if the SIR procedure is based on the parameter samples of a Bootstrap, SIR, or MarginalMCMC inference, as sir_from_bootstrap is. A proposal distribution that is built from a covariance matrix is used on the natural scale of that matrix, and neither keyword argument affects it, so neither of them has any effect on sir_results, sir_from_hessian, or sir_from_matrix.
How many iterations are needed cannot be known in advance. Dosne et al. (Dosne, Bergstrand, et al. 2016; Dosne et al. 2017) diagnose this with the same distribution of likelihood-ratio statistics as for the bootstrap above, here computed for the parameter samples. likelihood_ratio_ppplot shows it separately for the proposal samples of each iteration and for the samples resampled from them, which the plot labels Posterior.
likelihood_ratio_ppplot(sir_iterated)
As discussed above for the bootstrap, curves below the identity line contain parameter vectors that are less likely under the model and the data than the reference allows for. However, the initial proposal distribution should actually lie at or below the identity line, i.e., it should be at least as wide as the reference. Since the resampled parameters are a subset of the sampled ones, it is much harder for the SIR procedure to widen the spread of the samples than to narrow it, so a proposal distribution above the line is too narrow to start from and further iterations are not a reliable remedy. Dosne et al. (Dosne et al. 2017) inflate the proposal distribution in that case, which here means passing a scaled covariance matrix as proposal_covariance.
The SIR procedure is considered to have converged once the posterior curves of two consecutive iterations agree up to sampling noise (Dosne et al. 2017).
The methods for vcov and samples_table that we saw for Bootstrap also apply here
vcov(sir_results)9×9 Symmetric{Float64, Matrix{Float64}}:
3.62968e-5 0.000348386 2.40861e-5 … 1.7757e-6 5.02633e-5
0.000348386 0.0404937 0.0102667 0.00012874 -0.00115141
2.40861e-5 0.0102667 0.0190346 -3.35534e-5 0.00117162
-1.11595e-5 -0.000128123 -0.00116382 3.22481e-5 -2.8092e-5
5.70055e-5 0.000264664 0.000158533 -2.14363e-5 0.000235953
-3.59692e-6 0.000109225 6.55346e-5 … 4.76438e-7 -6.32331e-6
-0.000146238 -0.011341 -0.0142066 0.000241036 0.00284776
1.7757e-6 0.00012874 -3.35534e-5 5.72577e-5 2.50515e-5
5.02633e-5 -0.00115141 0.00117162 2.50515e-5 0.00148314and
samples_table(sir_results)| Row | iteration | resampled | pop_CL | pop_Vc | pop_tabs | pop_lag | pk_Ω₁,₁ | pk_Ω₂,₂ | pk_Ω₃,₃ | σ_prop | σ_add | loglikelihood |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Int64 | Bool | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | 1 | false | 0.139913 | 8.09898 | 0.270909 | 0.944072 | 0.0602871 | 0.0153502 | 1.71146 | 0.104718 | 0.43459 | -356.354 |
| 2 | 1 | false | 0.126388 | 7.43719 | 0.546844 | 0.913006 | 0.0552564 | 0.0158506 | 1.16432 | 0.103922 | 0.414951 | -359.261 |
| 3 | 1 | false | 0.144884 | 8.11891 | 0.532135 | 0.813961 | 0.110881 | 0.021031 | 0.843298 | 0.0723295 | 0.431184 | -356.796 |
| 4 | 1 | true | 0.134259 | 8.12845 | 0.663904 | 0.873656 | 0.0421814 | 0.0138431 | 1.47412 | 0.104963 | 0.402271 | -355.184 |
| 5 | 1 | false | 0.132217 | 8.20752 | 0.609049 | 0.858437 | 0.0501672 | 0.0121309 | 0.722428 | 0.0847403 | 0.387989 | -353.575 |
| 6 | 1 | false | 0.145849 | 8.0596 | 0.515502 | 0.877354 | 0.0879174 | 0.0186515 | 1.32292 | 0.0999536 | 0.455451 | -354.294 |
| 7 | 1 | true | 0.136892 | 8.09476 | 0.659616 | 0.836869 | 0.0874305 | 0.0280986 | 0.770005 | 0.0844024 | 0.427413 | -352.351 |
| 8 | 1 | false | 0.132131 | 8.02162 | 0.527528 | 0.951634 | 0.0605511 | 0.0244894 | 1.40016 | 0.10706 | 0.499345 | -359.483 |
| 9 | 1 | false | 0.147171 | 8.22759 | 0.628914 | 0.912468 | 0.0890172 | 0.0207625 | 1.48357 | 0.107325 | 0.523869 | -360.517 |
| 10 | 1 | false | 0.128671 | 8.11524 | 0.31241 | 0.961572 | 0.0465465 | 0.0242904 | 1.33969 | 0.100604 | 0.392312 | -375.195 |
| 11 | 1 | false | 0.138694 | 8.19763 | 0.821608 | 0.784441 | 0.11134 | 0.0108094 | 0.178263 | 0.0706693 | 0.408259 | -370.052 |
| 12 | 1 | false | 0.13054 | 8.26618 | 0.48309 | 0.825316 | 0.0715759 | 0.0112498 | 0.37904 | 0.0786218 | 0.319611 | -363.423 |
| 13 | 1 | false | 0.136568 | 8.03981 | 0.713437 | 0.829953 | 0.0765048 | 0.0188543 | 1.16659 | 0.0839245 | 0.465457 | -353.174 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| 989 | 1 | false | 0.12702 | 7.85043 | 0.504996 | 0.916498 | 0.0348412 | 0.0112806 | 0.710457 | 0.106329 | 0.387983 | -360.861 |
| 990 | 1 | false | 0.120026 | 7.73825 | 0.401402 | 0.920782 | 0.0358037 | 0.0331952 | 1.55216 | 0.0880829 | 0.357694 | -364.288 |
| 991 | 1 | false | 0.13984 | 7.58684 | 0.236183 | 0.935877 | 0.101818 | 0.0213419 | 1.30243 | 0.0995688 | 0.495436 | -359.812 |
| 992 | 1 | false | 0.136618 | 8.11007 | 0.652327 | 0.839113 | 0.102084 | 0.0213826 | 0.803312 | 0.0842942 | 0.45614 | -352.538 |
| 993 | 1 | true | 0.132678 | 8.26397 | 0.794173 | 0.896231 | 0.0752803 | 0.0193108 | 1.06417 | 0.0846101 | 0.44992 | -352.808 |
| 994 | 1 | false | 0.137411 | 8.05602 | 0.580874 | 0.858191 | 0.0703926 | 0.0197046 | 0.345727 | 0.0708231 | 0.345134 | -361.836 |
| 995 | 1 | false | 0.126043 | 7.71693 | 0.201182 | 0.942442 | 0.0410586 | 0.0231867 | 0.954445 | 0.0866231 | 0.325292 | -365.09 |
| 996 | 1 | false | 0.135913 | 7.8272 | 0.6272 | 0.904967 | 0.0717663 | 0.0199426 | 1.42365 | 0.104868 | 0.548703 | -359.714 |
| 997 | 1 | false | 0.14539 | 7.98263 | 0.862303 | 0.856531 | 0.0723463 | 0.0169899 | 1.22037 | 0.110906 | 0.518864 | -361.296 |
| 998 | 1 | false | 0.14473 | 8.34529 | 0.62758 | 0.847014 | 0.0893616 | 0.019536 | 0.444219 | 0.0754726 | 0.361905 | -358.389 |
| 999 | 1 | false | 0.133425 | 8.04954 | 0.336651 | 1.00642 | 0.0432638 | 0.0196364 | 1.51013 | 0.110695 | 0.443763 | -485.18 |
| 1000 | 1 | true | 0.13575 | 7.83892 | 0.550872 | 0.782704 | 0.110335 | 0.0196834 | 0.819436 | 0.0890779 | 0.446447 | -354.983 |
An alternative to Bootstrap and SIR is 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 a proposal distribution to be specified.
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.
mcmc_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.
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.