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

@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
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 = 10_000 as the number of MCMC samples to generate (including the adaptation samples); and nadapts = 2_000 as the number of adaptation steps in the NUTS algorithm which must be less than nsamples. Additionally, BayesMCMC by default uses 4 parallels chains with nchains = 4 and parallel_chains = true. Parallelism is also extended to subjects with parallel_subjects = true by default.

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

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 = 2_000, nadapts = 1_000, 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 (2000×6×4 Array{Float64, 3}):

Iterations        = 1:1:2000
Number of chains  = 4
Samples per chain = 2000
Wall duration     = 786.35 seconds
Compute duration  = 785.71 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.1979    0.0827    0.0044    349.7693    716.9232    1.0065   ⋯
tvv   13.2427    0.2779    0.0089    932.2831   1574.6139    1.0065   ⋯
ω²cl    0.0729    0.0087    0.0003    964.1235   1578.2516    1.0042   ⋯
ω²v    0.0460    0.0058    0.0001   1786.4488   3114.4032    1.0011   ⋯
ω²ka    1.0974    0.1586    0.0027   5326.4694   3775.9048    1.0000   ⋯
σ_p    0.1046    0.0095    0.0003   4295.5948   2915.6221    1.0002   ⋯
1 column omitted

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

tvcl    3.0445    3.1424    3.1948    3.2497    3.3736
tvv   12.6994   13.0674   13.2475   13.4226   13.7719
ω²cl    0.0580    0.0672    0.0725    0.0782    0.0914
ω²v    0.0360    0.0420    0.0456    0.0496    0.0582
ω²ka    0.8327    0.9983    1.0902    1.1926    1.4119
σ_p    0.0996    0.1024    0.1040    0.1057    0.1091

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.1979 0.0827039 0.00442662 349.769 716.923 1.00651 0.445163
2 tvv 13.2427 0.277872 0.00894386 932.283 1574.61 1.00648 1.18655
3 ω²cl 0.0729004 0.00870667 0.000282789 964.124 1578.25 1.00423 1.22707
4 ω²v 0.0459729 0.00578891 0.000137998 1786.45 3114.4 1.00108 2.27367
5 ω²ka 1.09741 0.158608 0.00269952 5326.47 3775.9 1.00004 6.77918
6 σ_p 0.104554 0.00953128 0.000274593 4295.59 2915.62 1.00018 5.46715
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.04449 3.14245 3.19476 3.24968 3.37361
2 tvv 12.6994 13.0674 13.2475 13.4226 13.7719
3 ω²cl 0.0579619 0.0671621 0.0724724 0.0782095 0.0914481
4 ω²v 0.0360289 0.041986 0.045597 0.0495636 0.0582459
5 ω²ka 0.832674 0.998338 1.09023 1.19263 1.41191
6 σ_p 0.0996258 0.102373 0.104004 0.105699 0.109082

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