using Pumas
Posterior Postprocessing - Plots and Queries
After you fit your Bayesian Pumas model, there are a number of functions and plots you can call on the output of the fit
or Pumas.truncate
function.
In this tutorial we’ll showcase some of them. As always, please refer to the Bayesian Workflow Pumas documentation for more details.
Let’s recall the model from the Introduction to Bayesian Models in Pumas tutorial:
= @model begin
pk_1cmp @param begin
~ LogNormal(log(3.2), 1)
tvcl ~ LogNormal(log(16.4), 1)
tvv ~ LogNormal(log(3.8), 1)
tvka ~ LogNormal(log(0.04), 0.25)
ω²cl ~ LogNormal(log(0.04), 0.25)
ω²v ~ LogNormal(log(0.04), 0.25)
ω²ka ∈ LogNormal(log(0.2), 0.25)
σ_p end
@random begin
~ Normal(0, sqrt(ω²cl))
ηcl ~ Normal(0, sqrt(ω²v))
ηv ~ Normal(0, sqrt(ω²ka))
ηka end
@covariates begin
Doseend
@pre begin
= tvcl * exp(ηcl)
CL = tvv * exp(ηv)
Vc = tvka * exp(ηka)
Ka end
@dynamics Depots1Central1
@derived begin
:= @. Central / Vc
cp ~ @. Normal(cp, abs(cp) * σ_p)
Conc 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
using PharmaDatasets
= dataset("pk_painrelief");
pkpain_df using DataFramesMeta
@chain pkpain_df begin
@rsubset! :Dose != "Placebo"
@rtransform! begin
:amt = :Time == 0 ? parse(Int, chop(:Dose; tail = 3)) : missing
:evid = :Time == 0 ? 1 : 0
:cmt = :Time == 0 ? 1 : 2
end
@rtransform! :Conc = :evid == 1 ? missing : :Conc
end;
first(pkpain_df, 5)
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 |
=
pop = read_pumas(
pkpain_noplb
pkpain_df,= :Subject,
id = :Time,
time = :amt,
amt = [:Conc],
observations = [:Dose],
covariates = :evid,
evid = :cmt,
cmt )
Population
Subjects: 120
Covariates: Dose
Observations: Conc
= fit(
pk_1cmp_fit
pk_1cmp,
pop,init_params(pk_1cmp),
BayesMCMC(; nsamples = 2_000, nadapts = 1_000, constantcoef = (; tvka = 2)),
);= Pumas.truncate(pk_1cmp_fit; burnin = 1_000) pk_1cmp_tfit
[ 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 = 1381.4 seconds
Compute duration = 1380.11 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.1896 0.0818 0.0068 147.5797 265.4966 1.0338 ⋯
tvv 13.2402 0.2683 0.0115 542.7381 995.9879 1.0038 ⋯
ω²cl 0.0738 0.0086 0.0004 575.2162 1117.9990 1.0013 ⋯
ω²v 0.0460 0.0054 0.0002 1116.2422 2080.7402 1.0018 ⋯
ω²ka 1.0977 0.1436 0.0022 4227.6316 3132.0972 1.0002 ⋯
σ_p 0.1040 0.0024 0.0000 5223.4984 2882.0876 1.0001 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
tvcl 3.0343 3.1338 3.1882 3.2430 3.3575
tvv 12.7207 13.0620 13.2322 13.4108 13.7975
ω²cl 0.0588 0.0680 0.0731 0.0787 0.0926
ω²v 0.0366 0.0423 0.0457 0.0493 0.0579
ω²ka 0.8339 0.9957 1.0915 1.1887 1.4040
σ_p 0.0994 0.1023 0.1040 0.1056 0.1088
1 Predictive Checks
In the Bayesian workflow, we generally perform predictive checks both in our prior selection and also in our posterior estimation (Gelman et al., 2013; McElreath, 2020).
1.1 Prior Predictive Checks
In a very simple way, prior predictive checks consists in simulate parameter values based on prior distribution without conditioning on any data or employing any likelihood function.
Independent of the level of information specified in the priors, it is always important to perform a prior sensitivity analysis in order to have a deep understanding of the prior influence onto the posterior.
To perform a prior predictive check in Pumas, we need to use the simobs
function with:
positional arguments:
- model
Population
orSubject
keyword arguments:
samples
: an integer; how many samples to drawsimulate_error
a boolean; whether to sample from the error model in the@derived block
(true
) or to simply return the expected value of the error distribution (false
)
In the case of a prior predictive check, we’ll use the following options in simobs
:
= simobs(pk_1cmp, pop; samples = 200, simulate_error = true) prior_sims
We get back a Vector
of SimulatedObservations
:
typeof(prior_sims)
Vector{SimulatedObservations{PumasModel{(tvcl = 1, tvv = 1, tvka = 1, ω²cl = 1, ω²v = 1, ω²ka = 1, σ_p = 1), 3, (:Depot, :Central), (:Ka, :CL, :Vc), ParamSet{NamedTuple{(:tvcl, :tvv, :tvka, :ω²cl, :ω²v, :ω²ka, :σ_p), NTuple{7, LogNormal{Float64}}}}, var"#1#7", TimeDispatcher{var"#2#8", var"#3#9"}, Nothing, var"#4#10", Depots1Central1, var"#5#11", var"#6#12", Nothing, PumasModelOptions}, Subject{NamedTuple{(:Conc,), Tuple{Vector{Union{Missing, Float64}}}}, ConstantCovar{NamedTuple{(:Dose,), Tuple{String7}}}, Vector{Event{Float64, Float64, Float64, Float64, Float64, Float64, Symbol}}, Vector{Float64}}, Vector{Float64}, NamedTuple{(:tvcl, :tvv, :tvka, :ω²cl, :ω²v, :ω²ka, :σ_p), NTuple{7, Float64}}, NamedTuple{(:ηcl, :ηv, :ηka), Tuple{Float64, Float64, Float64}}, NamedTuple{(:Dose,), Tuple{Vector{String7}}}, NamedTuple{(:callback, :continuity, :saveat, :alg), Tuple{Nothing, Symbol, Vector{Float64}, CompositeAlgorithm{Tuple{Tsit5{typeof(trivial_limiter!), typeof(trivial_limiter!), False}, Rosenbrock23{1, true, GenericLUFactorization{RowMaximum}, typeof(DEFAULT_PRECS), Val{:forward}, true, nothing}}, AutoSwitch{Tsit5{typeof(trivial_limiter!), typeof(trivial_limiter!), False}, Rosenbrock23{1, true, GenericLUFactorization{RowMaximum}, typeof(DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}}, PKPDAnalyticalSolution{SVector{2, Float64}, 2, Vector{SVector{2, Float64}}, Vector{Float64}, Vector{SVector{2, Float64}}, Vector{SVector{2, Float64}}, Returns{NamedTuple{(:Ka, :CL, :Vc), Tuple{Float64, Float64, Float64}}}, AnalyticalPKPDProblem{SVector{2, Float64}, Float64, false, Depots1Central1, Vector{Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}, Returns{NamedTuple{(:Ka, :CL, :Vc), Tuple{Float64, Float64, Float64}}}}}, NamedTuple{(), Tuple{}}, NamedTuple{(:CL, :Vc, :Ka), Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}}, NamedTuple{(:Conc,), Tuple{Vector{Float64}}}}} (alias for Array{Pumas.SimulatedObservations{PumasModel{(tvcl = 1, tvv = 1, tvka = 1, ω²cl = 1, ω²v = 1, ω²ka = 1, σ_p = 1), 3, (:Depot, :Central), (:Ka, :CL, :Vc), ParamSet{NamedTuple{(:tvcl, :tvv, :tvka, :ω²cl, :ω²v, :ω²ka, :σ_p), NTuple{7, LogNormal{Float64}}}}, var"#1#7", Pumas.TimeDispatcher{var"#2#8", var"#3#9"}, Nothing, var"#4#10", Depots1Central1, var"#5#11", var"#6#12", Nothing, Pumas.PumasModelOptions}, Subject{NamedTuple{(:Conc,), Tuple{Array{Union{Missing, Float64}, 1}}}, Pumas.ConstantCovar{NamedTuple{(:Dose,), Tuple{InlineStrings.String7}}}, Array{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Symbol}, 1}, Array{Float64, 1}}, Array{Float64, 1}, NamedTuple{(:tvcl, :tvv, :tvka, :ω²cl, :ω²v, :ω²ka, :σ_p), NTuple{7, Float64}}, NamedTuple{(:ηcl, :ηv, :ηka), Tuple{Float64, Float64, Float64}}, NamedTuple{(:Dose,), Tuple{Array{InlineStrings.String7, 1}}}, NamedTuple{(:callback, :continuity, :saveat, :alg), Tuple{Nothing, Symbol, Array{Float64, 1}, CompositeAlgorithm{Tuple{OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{1, true, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitch{OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{1, true, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}}, Pumas.PKPDAnalyticalSolution{StaticArraysCore.SArray{Tuple{2}, Float64, 1, 2}, 2, Array{StaticArraysCore.SArray{Tuple{2}, Float64, 1, 2}, 1}, Array{Float64, 1}, Array{StaticArraysCore.SArray{Tuple{2}, Float64, 1, 2}, 1}, Array{StaticArraysCore.SArray{Tuple{2}, Float64, 1, 2}, 1}, Pumas.Returns{NamedTuple{(:Ka, :CL, :Vc), Tuple{Float64, Float64, Float64}}}, Pumas.AnalyticalPKPDProblem{StaticArraysCore.SArray{Tuple{2}, Float64, 1, 2}, Float64, false, Depots1Central1, Array{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}, 1}, Array{Float64, 1}, Pumas.Returns{NamedTuple{(:Ka, :CL, :Vc), Tuple{Float64, Float64, Float64}}}}}, NamedTuple{(), Tuple{}}, NamedTuple{(:CL, :Vc, :Ka), Tuple{Array{Float64, 1}, Array{Float64, 1}, Array{Float64, 1}}}, NamedTuple{(:Conc,), Tuple{Array{Float64, 1}}}}, 1})
which we can then perform a visual predictive check (VPC):
To use the VPC plots, you have to load the PumasUtilities
package.
= vpc(prior_sims; stratify_by = [:Dose]); prior_vpc
using PumasUtilities
Most of Pumas’ plots can have the axes linked by using the facet
keyword argument:
= (; linkaxes = true) facet
You can also control which axis is linked by using:
linkxaxes
: only the X-axislinkyaxes
: only the Y-axislinkaxes
: both X- and Y-axis
vpc_plot(prior_vpc; facet = (; linkaxes = true,))
1.2 Posterior Predictive Checks
We need to make sure that the MCMC-approximated posterior distribution of our dependent variable (DV) can capture all the nuances of the real distribution density/mass of DV.
This procedure is called posterior predictive check, and it is generally carried on by a visual inspection of the real density/mass of our DV against generated samples of DV by the Bayesian model.
We also perform other inspections, see the Model Comparison with Crossvalidation tutorial.
To perform a posterior predictive check in Pumas, we need to use the simobs
function with:
positional arguments:
- result from a
fit
/Pumas.truncate
using a Bayesian Pumas model
- result from a
keyword arguments:
samples
: an integer; how many samples to drawsimulate_error
a boolean; whether to sample from the error model in the@derived block
(true
) or to simply return the expected value of the error distribution (false
).
In the case of a posterior predictive check, we’ll use the following options in simobs
:
= simobs(pk_1cmp_tfit; samples = 200, simulate_error = true) posterior_sims
= vpc(posterior_sims; stratify_by = [:Dose]); posterior_vpc
vpc_plot(posterior_vpc; facet = (; linkaxes = true,))
2 Posterior Queries
Often you want to execute posterior queries. A common posterior query is, for example, whether a certain parameter \(\theta\) is higher than \(0\):
\[\operatorname{E}[\theta > 0 \mid \text{data}]\]
The way we can do this is using the mean
function with a convenient do
operator:
mean(pk_1cmp_tfit) do p
>= 3
p.tvcl end
0.996
It also works for between-subject variability parameters:
mean(pk_1cmp_tfit) do p
>= 0.27^2
p.ω²cl end
0.5115
Or random effects in a specific subject:
mean(pk_1cmp_tfit; subject = 1) do p
<= 0
p.ηcl end
1.0
You can use any summarizing function in a posterior query. This means any function that takes a vector of numbers and outputs a single number; even user-defined functions.
3 Visualizations
There are a number of plots which one can use to gain insights into how the sampler and model are performing.
To use the following plots, you have to load the PumasUtilities
package.
The Bayesian visualization plots can have the axes linked by using the following keyword arguments:
linkxaxes
: only the X-axislinkyaxes
: only the Y-axis
For all of those you can use the following options:
:all
: link all axes:minimal
: link only within columns or rows:none
: unlink all axes.
The trace plot of a parameter shows the value of the parameter in each iteration of the MCMC algorithm. A good trace plot is one that:
- Is noisy, not an increasing or decreasing line for example.
- Has a fixed mean.
- Has a fixed variance.
- Shows all chains overlapping with each other, i.e. chain mixing.
You can plot a trace plot for the population parameters :tvcl
and :tvv
using:
trace_plot(pk_1cmp_tfit; parameters = [:tvcl, :tvv], linkyaxes = :none)
You can also plot the trace plots of the subject-specific parameters for a selection of subjects using:
trace_plot(pk_1cmp_tfit; subjects = [1, 2], linkyaxes = :none)
The density plot of a parameter shows a smoothed version of the histogram of the parameter values, giving an approximate probability density function for the marginal posterior of the parameter considered. This helps us visualize the shape of the marginal posterior of each parameter.
You can plot a density plot for the population parameter :tvcl
using:
density_plot(pk_1cmp_tfit; parameters = [:tvcl])
You can also plot the density plots of the subject-specific parameters for a selection of subjects using:
density_plot(pk_1cmp_tfit; subjects = [1, 2], linkxaxes = :none, linkyaxes = :none)
Another plot that can be used for parameter estimates is the ridgeline plot, which outputs a single density, even if you are using multiple chains, along with relevant statistical information about your parameter. The information that it outputs is mean, median, 80% quantiles (Q0.1
and Q0.9
), along with 95% and 80% highest posterior density interval(HPDI). HPDIs are the narrowest interval that capture certain percentage of the posterior density. You can plot ridge plots with the function ridgeline_plot
, which has a similar syntax as density_plot
.
Here’s an example with the population parameters :tvcl
:
ridgeline_plot(pk_1cmp_tfit; parameters = [:tvcl])
MCMC chains are prone to auto-correlation between the samples because each sample in the chain is a function of the previous sample.
You can plot an auto-correlation plot for the population parameter :tvcl
using:
autocor_plot(pk_1cmp_tfit; parameters = [:tvcl])
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.