Bayesian Estimation with Pumas

PumasAI
October, 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

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 be estimated. For this tutorial we will retrict the focus to a simple pharmacokinetics (PK) model. The Bayesian estimation procedures also works for more complicated models but such models might take a long time to estimate with an MCMC based procedure.

The prior for 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 normal but constrained to avoid negative and too extreme draws. The variance of the random effects vector η is modeled as an InverseWishart which is parameterized with degrees of freedom parmeter ν and scale matrix Ψ. It is sometimes more intuitive use the mode to specify priors and which is our choice. Since the mode of the InverseWishart is Ψ/(ν + p + 1), the model of our prior is diagm(fill(0.2, 3)) which corresponds to a between 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 so the posterior wont be forced away from zero. The scale parameter is then also equal to the mean which makes the interpretation easy.

Fitting models

To fit the model, we use the fit function. It requires a model, a population, a named tuple of parameters and a likelihood approximation method, since we want to use Bayesian estimation with Markov Chain Monte Carlo we pass BayesMCMC as the likelihood approximation argument.

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 estimate the model 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.8803    0.5259     0.0118    0.0127   1067.8560    0.9995  
          θ₂    0.0415    0.0076     0.0002    0.0003    699.0790    0.9998  
          θ₃    0.4712    0.0590     0.0013    0.0016    867.9547    1.0005  
          θ₄    1.8017    0.4868     0.0109    0.0093   1016.2725    1.0000  
        Ω₁,₁    0.5802    0.2587     0.0058    0.0060   1192.9478    1.0004  
        Ω₂,₁   -0.0494    0.1154     0.0026    0.0039    996.3132    0.9995  
        Ω₃,₁    0.0045    0.0888     0.0020    0.0017   1065.6074    1.0014  
        Ω₂,₂    0.2569    0.1141     0.0026    0.0032   1049.3889    0.9998  
        Ω₃,₂    0.0282    0.0625     0.0014    0.0018   1142.6894    0.9995  
        Ω₃,₃    0.1761    0.0815     0.0018    0.0026    977.5103    0.9996  
           σ    0.7054    0.0654     0.0015    0.0021    777.7266    1.0020  

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%  
      Symbol   Float64   Float64   Float64   Float64   Float64  
                                                                
          θ₁    1.0235    1.5123    1.8104    2.1939    3.0681  
          θ₂    0.0310    0.0369    0.0408    0.0452    0.0547  
          θ₃    0.3725    0.4319    0.4672    0.5061    0.5976  
          θ₄    1.0351    1.4508    1.7494    2.0726    2.8994  
        Ω₁,₁    0.2572    0.3957    0.5207    0.7015    1.2060  
        Ω₂,₁   -0.2971   -0.1081   -0.0427    0.0174    0.1787  
        Ω₃,₁   -0.1723   -0.0467    0.0039    0.0552    0.1887  
        Ω₂,₂    0.1179    0.1781    0.2323    0.3058    0.5564  
        Ω₃,₂   -0.0888   -0.0072    0.0245    0.0595    0.1646  
        Ω₃,₃    0.0857    0.1245    0.1574    0.2019    0.3633  
           σ    0.6091    0.6650    0.7011    0.7377    0.8197

The show method for the fit will print various summary statistics useful for evaluating the fit. The summary output is based on the the postprocessing by the MCMCChains.jl package. Hence, 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])