Introduction to Bayesian Models in Pumas

Authors

Jose Storopoli

Mohamed Tarek 

using Pumas

1 Introduction to Bayesian Models in Pumas

This notebook will provide the basics on how to specify and fit Bayesian models in Pumas.

1.1 Bayesian Statistics

Bayesian statistics is a data analysis approach based on Bayes’ theorem where available knowledge about the parameters of a statistical model is updated with the information of observed data (Gelman et al., 2013; McElreath, 2020).

Previous knowledge is expressed as a prior distribution and combined with the observed data in the form of a likelihood function to generate a posterior distribution. The posterior can also be used to make predictions about future events.

1.2 Advantages of Bayesian analysis

The Bayesian workflow allows analysts to:

  1. Incorporate domain knowledge and insights from previous studies using prior distributions.
  2. Quantify the epistemic uncertainty in the model parameters’ values. Parameter values can be uncertain in reality due to model non-identifiability or practical non-identifiability issues because of lack of data.

This epistemic uncertainty can be further propagated forward via simulation to get a distribution of predictions which when used to make decisions, e.g. dosing decisions, can make those decisions more robust to changes in the parameter values. The above are advantages of Bayesian analysis which the traditional frequentist workflow typically doesn’t have a satisfactory answer for. Using bootstrapping or asymptotic estimates of standard errors can be considered somewhat ad-hoc and assumptions methods to quantify uncertainty in the parameter estimates when Bayesian inference uses the established theory of probability to more rigorously quantify said uncertainty with fewer assumptions about the model.

It is important to note that one can still reap the second benefit of Bayesian analysis due to its flexibility even when little to no domain knowledge is imposed on the model. So the Bayesian workflow doesn’t force you to incorporate your domain knowledge but it empowers you to if you want.

Let’s formalize some of the advantages.

First, Bayesian Statistics uses probabilistic statements in terms of:

  • one or more parameters \(\theta\)
  • unobserved data \(\tilde{y}\)

These statements are conditioned on the observed values of \(y\):

  • \(P(\theta \mid y)\)
  • \(P(\tilde{y} \mid y)\)

We also, implicitly, condition on the observed data from any covariate \(x\)

Generally, we are interested in:

  • expected response of a new subject to a drug, e.g. \(\operatorname{E}[\hat{y} \mid y]\)
  • the probability of drug effect is higher than zero, e.g. \(P(\theta > 0 \mid y) \geq 0.95\)

1.3 How to specify Bayesian models in Pumas?

Let’s revisit the 1-compartment model we specified in the Introduction to Pumas tutorial:

@model begin
    @param begin
        tvcl ∈ RealDomain(lower = 0, init = 3.2)
        tvv ∈ RealDomain(lower = 0, init = 16.4)
        tvka ∈ RealDomain(lower = 0, init = 3.8)
        Ω ∈ PDiagDomain(init = [0.04, 0.04, 0.04])
        σ_p ∈ RealDomain(lower = 0.0001, init = 0.2)
    end
    @random begin
        η ~ MvNormal(Ω)
    end
    @covariates begin
        Dose
    end
    @pre begin
        CL = tvcl * exp(η[1])
        Vc = tvv * exp(η[2])
        Ka = tvka * exp(η[3])
    end
    @dynamics Depots1Central1
    @derived begin
        cp := @. Central / Vc
        Conc ~ @. Normal(cp, abs(cp) * σ_p)
    end
end

Bayesian models in Pumas are quite easy: you just need to specify priors for all parameters in the @param block.

Instead of having:

@param begin
    tvcl ∈ RealDomain(lower = 0, init = 3.2)
    ...
end

We would say that tvcl has a prior:

@param begin
    tvcl ~ LogNormal(log(1.5), 1)
    ...
end

where tvcl (typical value for clearance) is the parameter that we want the model to estimate and LogNormal(log(1.5), 1) is the prior we are giving to tvcl.

Note

You can use any prior you want. Here we are giving tvcl a LogNormal prior since it has a positive-only domain and allows us to give more weight to certain values while also making it possible for higher values due to the wide-tail nature.

There is no one prior that fits all cases and generally it might be a good practice to follow a previous similar study’s priors where a good reference can be found.

However, if you have the task of choosing good prior distributions for an all-new model, it will generally be a multi-step process consisting of:

  1. Deciding the support of the prior. The support of the prior distribution must match the the domain of the parameter. For example, different priors can be used for positive parameters than those for parameters between 0 and 1.
  2. Deciding the center of the prior, e.g. mean, median or mode.
  3. Deciding the strength (aka informativeness) of the prior. This is often controlled by a standard deviation or scale parameter in the distribution constructor. A small standard deviation or scale parameter implies low uncertainty in the parameter value which leads to a stronger (aka more informative) prior. A large standard deviation or scale parameter implies high uncertainty in the parameter value which leads to a weak (aka less informative) prior. You should study each prior distribution you are considering before using it to make sure the strength of the prior reflects your confidence level in the parameter values.
  4. Deciding the shape of the probability density function (PDF) of the prior. Some distributions are left skewed, others are right skewed and some are symmetric. Some have heavier tails than others. You should make sure the shape of the PDF reflects your knowledge about the parameter value prior to observing the data.

1.3.0.1 Data preparation for modeling

We’ll use the PharmaDatasets.jl package to load the data:

using PharmaDatasets
pkpain_df = dataset("pk_painrelief")
first(pkpain_df, 5)
5×7 DataFrame
Row Subject Time Conc PainRelief PainScore RemedStatus Dose
Int64 Float64 Float64 Int64 Int64 Int64 String7
1 1 0.0 0.0 0 3 1 20 mg
2 1 0.5 1.15578 1 1 0 20 mg
3 1 1.0 1.37211 1 0 0 20 mg
4 1 1.5 1.30058 1 0 0 20 mg
5 1 2.0 1.19195 1 1 0 20 mg

Let’s filter out the placebo data as we don’t need that for the PK analysis:

using DataFramesMeta
pkpain_noplb_df = @rsubset pkpain_df :Dose != "Placebo";
Tip

If you want to learn more about data wrangling, don’t forget to check our Data Wrangling in Julia tutorials!

Also we need to add the :amt column:

@rtransform! pkpain_noplb_df :amt = :Time == 0 ? parse(Int, chop(:Dose; tail = 3)) : missing;

PumasNDF requires the presence of :evid and :cmt columns in the dataset:

@rtransform! pkpain_noplb_df begin
    :evid = :Time == 0 ? 1 : 0
    :cmt = :Time == 0 ? 1 : 2
end;

Further, observations at time of dosing, i.e., when evid = 1 have to be missing:

@rtransform! pkpain_noplb_df :Conc = :evid == 1 ? missing : :Conc;

Here is the final result:

first(pkpain_noplb_df, 10)
10×10 DataFrame
Row Subject Time Conc PainRelief PainScore RemedStatus Dose amt evid cmt
Int64 Float64 Float64? Int64 Int64 Int64 String7 Int64? Int64 Int64
1 1 0.0 missing 0 3 1 20 mg 20 1 1
2 1 0.5 1.15578 1 1 0 20 mg missing 0 2
3 1 1.0 1.37211 1 0 0 20 mg missing 0 2
4 1 1.5 1.30058 1 0 0 20 mg missing 0 2
5 1 2.0 1.19195 1 1 0 20 mg missing 0 2
6 1 2.5 1.13602 1 1 0 20 mg missing 0 2
7 1 3.0 0.873224 1 0 0 20 mg missing 0 2
8 1 4.0 0.739963 1 1 0 20 mg missing 0 2
9 1 5.0 0.600143 0 2 0 20 mg missing 0 2
10 1 6.0 0.425624 1 1 0 20 mg missing 0 2

Finally, we’ll import pkpain_noplb_df DataFrame into a Population with read_pumas:

pkpain_noplb = read_pumas(
    pkpain_noplb_df,
    id = :Subject,
    time = :Time,
    amt = :amt,
    observations = [:Conc],
    covariates = [:Dose],
    evid = :evid,
    cmt = :cmt,
)
Population
  Subjects: 120
  Covariates: Dose
  Observations: Conc

1.3.1 Pumas Bayesian 1-compartment Model

Let’s then add the priors in the @param block.

Note

Since the original model in the Introduction to Pumas tutorial has Ω as a PDiagDomain, we’ll be using scalars and not a matrix for our between-subject variability parameters. Hence, we’ll define ω² for every diagonal element of the original Ω.

We’ll cover covariance matrices in the Random Effects in Bayesian Models tutorial.

Now we fit our model by using the fit function. Instead of using estimation method as FOCE or NaivePooled, we will be using BayesMCMC which uses No-U-Turn Sampler (NUTS) – a very fast and state-of-the-art Hamiltonian Monte Carlo (HMC) based sampler:

Note

If you want to learn more about Markov chain Monte Carlo (MCMC) methods, Hamiltonian Monte Carlo (HMC), and No-U-Turn Sampler (NUTS); please check the Pumas Bayesian Workflow Documentation at docs.pumas.ai.

pk_1cmp = @model begin
    @param begin
        tvcl ~ LogNormal(log(3.2), 1)
        tvv ~ LogNormal(log(16.4), 1)
        tvka ~ LogNormal(log(3.8), 1)
        ω²cl ~ LogNormal(log(0.04), 0.25)
        ω²v ~ LogNormal(log(0.04), 0.25)
        ω²ka ~ LogNormal(log(0.04), 0.25)
        σ_p ∈ LogNormal(log(0.2), 0.25)
    end
    @random begin
        ηcl ~ Normal(0, sqrt(ω²cl))
        ηv ~ Normal(0, sqrt(ω²v))
        ηka ~ Normal(0, sqrt(ω²ka))
    end
    @covariates begin
        Dose
    end
    @pre begin
        CL = tvcl * exp(ηcl)
        Vc = tvv * exp(ηv)
        Ka = tvka * exp(ηka)
    end
    @dynamics Depots1Central1
    @derived begin
        cp := @. Central / Vc
        Conc ~ @. Normal(cp, abs(cp) * σ_p)
    end
end
┌ Warning: Covariate Dose is not used in the model.
└ @ Pumas ~/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/aZRyj/src/dsl/model_macro.jl:2856
PumasModel
  Parameters: tvcl, tvv, tvka, ω²cl, ω²v, ω²ka, σ_p
  Random effects: ηcl, ηv, ηka
  Covariates: Dose
  Dynamical system variables: Depot, Central
  Dynamical system type: Closed form
  Derived: Conc
  Observed: Conc
Tip

By default BayesMCMC uses nsamples = 2_000 as the number of MCMC samples to generate (including the adaptation samples); and nadapts = 1_000 as the number of adaptation steps in the NUTS algorithm which must be less than nsamples. Additionally, BayesMCMC by default uses 4 chains with nchains = 4. By default, multithreading is enabled (ensemblealg = EnsembleThreads): If Julia runs multithreaded, sampling of the chains is performed in parallel (parallel_chains = true), and if the number of threads is at least twice the number of chains, then parallelization is extended to subjects (parallel_subjects = true).

Here we’ll use nsamples = 1_000 and nadapts = 500.

Caution

Similar to the original model fit in the Introduction to Pumas tutorial, we’ll fix tvka to 2 as we don’t have a lot of information before tmax.

pk_1cmp_fit = fit(
    pk_1cmp,
    pkpain_noplb,
    init_params(pk_1cmp),
    BayesMCMC(; nsamples = 1_000, nadapts = 500, constantcoef = (; tvka = 2)),
)
[ Info: Checking the initial parameter values.
[ Info: The initial log probability and its gradient are finite. Check passed.
[ Info: Checking the initial parameter values.
[ Info: The initial log probability and its gradient are finite. Check passed.
[ Info: Checking the initial parameter values.
[ Info: The initial log probability and its gradient are finite. Check passed.
[ Info: Checking the initial parameter values.
[ Info: The initial log probability and its gradient are finite. Check passed.
Chains MCMC chain (1000×6×4 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 4
Samples per chain = 1000
Wall duration     = 484.37 seconds
Compute duration  = 1816.82 seconds
parameters        = tvcl, tvv, ω²cl, ω²v, ω²ka, σ_p

Summary Statistics
  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ⋯
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64   ⋯

        tvcl    3.1965    0.0820    0.0061    171.1794    228.9697    1.0130   ⋯
         tvv   13.2594    0.2949    0.0163    294.3668    784.5622    1.0185   ⋯
        ω²cl    0.0726    0.0090    0.0005    386.8321    560.5880    1.0047   ⋯
         ω²v    0.0458    0.0056    0.0002    944.7915   1594.0172    1.0001   ⋯
        ω²ka    1.0875    0.1748    0.0058   2087.6109   1141.6907    1.0012   ⋯
         σ_p    0.1054    0.0215    0.0007   1629.6766    862.5197    1.0011   ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

        tvcl    3.0534    3.1441    3.1914    3.2462    3.3626
         tvv   12.7116   13.0851   13.2619   13.4419   13.7885
        ω²cl    0.0564    0.0671    0.0724    0.0783    0.0907
         ω²v    0.0358    0.0420    0.0455    0.0493    0.0575
        ω²ka    0.8205    0.9914    1.0859    1.1842    1.4075
         σ_p    0.0997    0.1024    0.1040    0.1057    0.1094
Note

Pumas Bayesian models returns a different object that of other Pumas models which can be:

  • exported to a DataFrame containing all of the MCMC samples with DataFrame(Chains(fit_result))
  • summarized to a DataFrame containing summary statistics of the MCMC samples with DataFrame(summarystats(fit_result))
  • summarized to a DataFrame containing quantiles of the MCMC samples with DataFrame(quantile(fit_result))

Here are our estimates both as summary statistics with summarystats and quantiles with quantile:

DataFrame(summarystats(pk_1cmp_fit))
6×8 DataFrame
Row parameters mean std mcse ess_bulk ess_tail rhat ess_per_sec
Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64
1 tvcl 3.19648 0.0819599 0.00609143 171.179 228.97 1.01302 0.0942191
2 tvv 13.2594 0.294855 0.0163149 294.367 784.562 1.0185 0.162023
3 ω²cl 0.0726201 0.0089747 0.000464749 386.832 560.588 1.00468 0.212917
4 ω²v 0.0458332 0.00559946 0.000182312 944.792 1594.02 1.00013 0.520024
5 ω²ka 1.08754 0.17484 0.00575527 2087.61 1141.69 1.00115 1.14904
6 σ_p 0.10538 0.0215159 0.000720165 1629.68 862.52 1.00112 0.896992 
DataFrame(quantile(pk_1cmp_fit))
6×6 DataFrame
Row parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
1 tvcl 3.05343 3.14405 3.19144 3.24625 3.36262
2 tvv 12.7116 13.0851 13.2619 13.4419 13.7885
3 ω²cl 0.056393 0.067066 0.0723587 0.0782869 0.0906654
4 ω²v 0.0357811 0.0420332 0.0455254 0.049348 0.0575088
5 ω²ka 0.820523 0.991375 1.08591 1.18424 1.40755
6 σ_p 0.0996584 0.10245 0.104041 0.105706 0.109388

1.4 References

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis. Chapman and Hall/CRC.

McElreath, R. (2020). Statistical rethinking: A Bayesian course with examples in R and Stan. CRC press.