Bayesian Estimation with Pumas
PumasAI
November, 2020
Fitting a PK model
In this tutorial we will go through the steps for bayesian parameter estimation of the Theophylline model in Pumas.jl. During the tutorial we use the following packages.
using Pumas, MCMCChains, StatsPlots
Important Note: Pumas.jl is a proprietary package. It is free to use for non-commercial academic teaching and research purposes. For commercial users, license fees apply. Further, Pumas may not be used under any license agreement for patient management services in a hospital or clinical settings. Please refer to End User License Agreement (https://juliacomputing.com/eula) for details. Please contact sales@juliacomputing.com for purchase.
The core modeling and estimation functionality is provided by the Pumas package. Postprocessing of the posterior data is handled by the MCMCChains package while the plots are made with the Plots package.
The PK model of drug concentration and elimination
First, we have to set up the pharmacometric model to estimate its parameters. For this tutorial, we will retrict the focus to a simple pharmacokinetics (PK) model. The Bayesian estimation procedure also works for more complicated models but such models might take longer to estimate with an MCMC-based procedure.
The prior distribution of a parameter of the model can be any distribution from the Distributions.jl package and is specified with the tilde (~) symbol.
theopmodel_bayes = @model begin @param begin # Mode at [2.0, 0.2, 0.8, 2.0] θ ~ Constrained( MvNormal( [2.0, 0.2, 0.8, 2.0], Diagonal(ones(4)) ), lower = zeros(4), upper = fill(10.0, 4), init = [2.0, 0.2, 0.8, 2.0]) # Mode at diagm(fill(0.2, 3)) Ω ~ InverseWishart(6, diagm(fill(0.2, 3)) .* (6 + 3 + 1)) # Mean at 0.5 and positive density at 0.0 σ ~ Gamma(1.0, 0.5) end @random begin η ~ MvNormal(Ω) end @pre begin Ka = (SEX == 1 ? θ[1] : θ[4])*exp(η[1]) CL = θ[2]*(WT/70) *exp(η[2]) Vc = θ[3] *exp(η[3]) end @covariates SEX WT @dynamics Depots1Central1 @derived begin # The conditional mean μ := @. Central / Vc # Additive error model dv ~ @. Normal(μ, σ) end end
PumasModel Parameters: θ, Ω, σ Random effects: η Covariates: SEX, WT Dynamical variables: Depot, Central Derived: dv Observed: dv
The joint prior distribution of θ is specified as the normal distribution but constrained to avoid negative and extreme draws. In the Bayesian framework, the prior distribution represents our initial belief about the value of the parameter prior to observing any data. The covariance matrix of the random effects vector η is given a prior distribution of InverseWishart which is parameterized with degrees of freedom parameter ν and scale matrix Ψ. It is sometimes more intuitive to use the mode to specify priors which is our choice. Since the mode of InverseWishart is Ψ/(ν + p + 1), the mode of our prior is diagm(fill(0.2, 3)) which corresponds to an inter-subject variability peaking at 20%. Finally, the prior additive error component σ is modeled as a Gamma(1.0, 0.5). Notice that Gamma is parameterized with a shape and a scale parameter and not a rate. Setting the shape parameter equal to one has the advantage that the density has support at zero, meaning the value zero has a non-zero probablility, so the posterior distribution won't be forced away from zero. The scale parameter is then also equal to the mean which makes the interpretation easy. In the Bayesian framework, the posterior distribution represents our belief about the value of the parameter after observing some data.
Fitting models
To fit the model, we use the fit function. It requires a model, a population, a named tuple of parameters and an estimation method. Since we want to use Bayesian estimation with Markov Chain Monte Carlo (MCMC), we pass BayesMCMC as the estimation method argument. This will return a sample from the joint posterior distribution of the model's parameters. A sample in the MCMC context is also known as a chain.
First, we load the theophylline dataset and specify SEX and WT as covariates.
data = read_pumas(example_data("event_data/THEOPP"), covariates = [:SEX,:WT])
Population Subjects: 12 Covariates: SEX, WT Observables: dv
Next, we extract the initial parameters from theopmodel_bayes. Alternatively, the initial values could be specified directly as a NamedTuple with keys matching the parameter names used in theopmodel_bayes.
param = init_param(theopmodel_bayes)
(θ = [2.0, 0.2, 0.8, 2.0], Ω = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], σ = 0.5)
We can now infer the model's parameters given the data using the No U-Turn Sampler (NUTS). The number of samples and the number of steps used for adaptation are set with the nsamples and nadapts arguments respectively. The adaptation uses the prodecure introduced by Stan and it is currently not possible to adjust the adaption parameters.
result = fit(theopmodel_bayes, data, param, Pumas.BayesMCMC(); nsamples=2000, nadapts=1000)
Chains MCMC chain (2000×11×1 Array{Float64,3}):
Iterations = 1:2000
Thinning interval = 1
Chains = 1
Samples per chain = 2000
parameters = θ₁, θ₂, θ₃, θ₄, Ω₁,₁, Ω₂,₁, Ω₃,₁, Ω₂,₂, Ω₃,₂, Ω₃,₃, σ
Summary Statistics
parameters mean std naive_se mcse ess rhat
Symbol Float64 Float64 Float64 Float64 Float64 Float64
θ₁ 1.8721 0.4879 0.0109 0.0140 1058.7378 0.9995
θ₂ 0.0417 0.0080 0.0002 0.0003 750.2298 1.0013
θ₃ 0.4752 0.0613 0.0014 0.0018 909.7713 1.0000
θ₄ 1.8088 0.4659 0.0104 0.0143 850.1282 0.9995
Ω₁,₁ 0.5749 0.2584 0.0058 0.0074 1219.7135 0.9995
Ω₂,₁ -0.0546 0.1098 0.0025 0.0045 965.5898 1.0008
Ω₃,₁ -0.0078 0.0941 0.0021 0.0022 1030.1632 0.9997
Ω₂,₂ 0.2573 0.1093 0.0024 0.0034 1222.8134 1.0009
Ω₃,₂ 0.0296 0.0651 0.0015 0.0024 1074.5796 1.0000
Ω₃,₃ 0.1790 0.0824 0.0018 0.0031 825.4327 1.0028
σ 0.7050 0.0627 0.0014 0.0020 756.6869 1.0005
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
θ₁ 1.0609 1.5257 1.8051 2.1496 2.9694
θ₂ 0.0301 0.0372 0.0411 0.0452 0.0557
θ₃ 0.3668 0.4334 0.4709 0.5102 0.6072
θ₄ 1.0529 1.4840 1.7697 2.0877 2.8534
Ω₁,₁ 0.2496 0.3956 0.5226 0.6889 1.2473
Ω₂,₁ -0.2847 -0.1079 -0.0510 0.0072 0.1632
Ω₃,₁ -0.2185 -0.0592 -0.0062 0.0474 0.1793
Ω₂,₂ 0.1238 0.1830 0.2345 0.3045 0.5528
Ω₃,₂ -0.0849 -0.0078 0.0242 0.0597 0.1827
Ω₃,₃ 0.0852 0.1269 0.1607 0.2083 0.3883
σ 0.6128 0.6666 0.7009 0.7361 0.8172
The show method for the result will print various summary statistics useful for evaluating the sample/chain. The summary output is based on some postprocessing by the MCMCChains.jl package. Hence, a similar output will be presented if the fitted model is converted to a Chains object from MCMCChains.
chains = Pumas.Chains(result)
However, the MCMCChains package provides many other out of the box diagnostics and plotting functionality for MCMC chains. E.g. a default plotting method to that generates time series and kernel density plots for each of the parmeters.
plot(chains)
The plots are based on the complete chain, including the very unstable initial phase of the sampling. This can distort both the scale of the of the times series plot and the shape of the kernel density. Hence, it is often useful to exlcude the initial burn-in phase. This is easily done simply by slicing the Chains structure, i.e.
plot(chains[200:end])