I am using the following code from the Pumas tutorials website for fitting a Bayesian model for theophylline

## Code

```
using Pumas, MCMCChains
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
df = CSV.read("/mnt/data/code/tutorial/THEOPP.csv", DataFrame;
missingstring = ["."])
data = read_pumas(df, covariates = [:SEX,:WT])
param = init_params(theopmodel_bayes)
result = fit(theopmodel_bayes, data, param, Pumas.BayesMCMC();
nsamples=2000, nadapts=1000)
```

I have the following questions regarding post-processing and qualification of the model:

- Is there a direct way to get predicted concentrations to qualify the model using individual fits? - this is equivalent to running the
`inspect(fpm)`

in the maximum likelihood estimation (MLE) approach - For covariate inclusion, I would like to test the Bayes factor approach. Is there a way to calculate the said factor in Pumas using the two models (one with the covariate effect and another without the covariate effect)? - this is an equivalent of the log-likelihood ratio in the MLE approach.
- How to generate the mode and 95% credible intervals for a) individual parameters and b) observations?
- What information does the
`result.chains[chain_num][iteration_num]`

object contain? In this case, it seems to be a vector with 47 numbers but I do not know what those numbers are representing. - I keep seeing the warning when fitting -
`Warning: the current proposal will be rejected due to numerical error(s)`

. When I ran a fit with`nsamples=500, nadapts=250`

, I saw this warning for approximately 100 times until the end. - A last general question - are there simple tutorials to understand the NUTS sampling, MCMC (or HMC) stuff for people without a background in Bayesian population PK ?