Posterior Postprocessing - Plots and Queries

Authors

Jose Storopoli

Mohamed Tarek 

using Pumas

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:

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
using PharmaDatasets
pkpain_df = dataset("pk_painrelief");
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)
5×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 
pop =
    pkpain_noplb = read_pumas(
        pkpain_df,
        id = :Subject,
        time = :Time,
        amt = :amt,
        observations = [:Conc],
        covariates = [:Dose],
        evid = :evid,
        cmt = :cmt,
    )
Population
  Subjects: 120
  Covariates: Dose
  Observations: Conc
pk_1cmp_fit = fit(
    pk_1cmp,
    pop,
    init_params(pk_1cmp),
    BayesMCMC(; nsamples = 2_000, nadapts = 1_000, constantcoef = (; tvka = 2)),
);
pk_1cmp_tfit = Pumas.truncate(pk_1cmp_fit; burnin = 1_000)
[ 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     = 1084.52 seconds
Compute duration  = 1083.62 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.1870    0.0790    0.0052    229.0134    452.5206    1.0147   ⋯
         tvv   13.2594    0.2660    0.0109    598.0363   1091.9517    1.0063   ⋯
        ω²cl    0.0736    0.0079    0.0003    663.2039   1202.4435    1.0059   ⋯
         ω²v    0.0459    0.0055    0.0002   1159.2336   1832.3081    1.0017   ⋯
        ω²ka    1.0986    0.1426    0.0025   3293.5732   3578.0865    1.0003   ⋯
         σ_p    0.1040    0.0023    0.0000   4861.6934   3240.1004    1.0013   ⋯
                                                                1 column omitted

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

        tvcl    3.0316    3.1330    3.1880    3.2414    3.3407
         tvv   12.7290   13.0871   13.2629   13.4350   13.7741
        ω²cl    0.0595    0.0680    0.0731    0.0785    0.0906
         ω²v    0.0359    0.0421    0.0457    0.0493    0.0580
        ω²ka    0.8503    1.0001    1.0866    1.1824    1.4091
         σ_p    0.0995    0.1024    0.1039    0.1056    0.1087

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:

    1. model
    2. Population or Subject

  • keyword arguments:

    • samples: an integer; how many samples to draw
    • simulate_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:

prior_sims = simobs(pk_1cmp, pop; samples = 200, simulate_error = true)

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):

Note

To use the VPC plots, you have to load the PumasUtilities package.

prior_vpc = vpc(prior_sims; stratify_by = [:Dose]);
using PumasUtilities
Tip

Most of Pumas’ plots can have the axes linked by using the facet keyword argument:

facet = (; linkaxes = true)

You can also control which axis is linked by using:

  • linkxaxes: only the X-axis
  • linkyaxes: only the Y-axis
  • linkaxes: 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.

Note

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:

    1. result from a fit/Pumas.truncate using a Bayesian Pumas model

  • keyword arguments:

    • samples: an integer; how many samples to draw
    • simulate_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:

posterior_sims = simobs(pk_1cmp_tfit; samples = 200, simulate_error = true)
posterior_vpc = vpc(posterior_sims; stratify_by = [:Dose]);
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
    p.tvcl >= 3
end
0.98875

It also works for between-subject variability parameters:

mean(pk_1cmp_tfit) do p
    p.ω²cl >= 0.27^2
end
0.51

Or random effects in a specific subject:

mean(pk_1cmp_tfit; subject = 1) do p
    p.ηcl <= 0
end
1.0
Tip

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.

Note

To use the following plots, you have to load the PumasUtilities package.

Tip

The Bayesian visualization plots can have the axes linked by using the following keyword arguments:

  • linkxaxes: only the X-axis
  • linkyaxes: 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.