# Individual PK Parameters following Bayesian Estimation

**URL:** <https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665>\
**Category:** bayesian\
**Created:** [April 7, 2022, 5:44pm UTC](https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665 "2022-04-07T17:44:42Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![ahmed.salem](https://avatars.discourse-cdn.com/v4/letter/a/858c86/32.png) [@ahmed.salem](https://discourse.pumas.ai/u/ahmed.salem)\
**Post date:** [April 7, 2022, 5:44pm UTC](https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665/1 "2022-04-07T17:44:42Z")

</div>

Hello,  
I am trying to run inspect function following Bayesian estimation to get the individual random effects but I am getting the following error

```auto
ERROR: type BayesMCMCResults has no field approx
Stacktrace:
 [1] getproperty(x::Pumas.BayesMCMCResults{Pumas.BayesLogDensity{PumasModel{ParamSet{NamedTuple{(:θ, :Ω, :σ), Tuple{Constrained{DiagNormal, VectorDomain{Vector{Float64}, Vector{Float64}, Vector{Float64}}}, InverseWishart{Float64, PDMats.PDMat{Float64, Matrix{Float64}}}, Gamma{Float64}}}}, var"#5#17", var"#6#18", var"#8#20", var"#10#22", Central1, var"#11#23", var"#14#26", Nothing}, Vector{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, T2, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}} where T2}, Vector{Float64}, ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(Pumas.logdensity), Float64}, Float64, 8, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(Pumas.logdensity), Float64}, Float64, 8}}}, DiffResults.MutableDiffResult{1, Float64, Tuple{Vector{Float64}}}, NamedTuple{(), Tuple{}}}, Vector{Vector{Vector{Float64}}}, Vector{Vector{NamedTuple}}}, f::Symbol)
   @ Base ./Base.jl:33
 [2] inspect(fpm::Pumas.BayesMCMCResults{Pumas.BayesLogDensity{PumasModel{ParamSet{NamedTuple{(:θ, :Ω, :σ), Tuple{Constrained{DiagNormal, VectorDomain{Vector{Float64}, Vector{Float64}, Vector{Float64}}}, InverseWishart{Float64, PDMats.PDMat{Float64, Matrix{Float64}}}, Gamma{Float64}}}}, var"#5#17", var"#6#18", var"#8#20", var"#10#22", Central1, var"#11#23", var"#14#26", Nothing}, Vector{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, T2, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}} where T2}, Vector{Float64}, ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(Pumas.logdensity), Float64}, Float64, 8, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(Pumas.logdensity), Float64}, Float64, 8}}}, DiffResults.MutableDiffResult{1, Float64, Tuple{Vector{Float64}}}, NamedTuple{(), Tuple{}}}, Vector{Vector{Vector{Float64}}}, Vector{Vector{NamedTuple}}})
   @ Pumas /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/HDuXQ/src/estimation/diagnostics.jl:1096
 [3] top-level scope
   @ REPL[13]:1

```

Is there a different way to run inspect function specific to Bayesian estimation ?

---

<div class="post-metadata">

**Author:** ![storopoli](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/storopoli/32/190_2.png) [@storopoli](https://discourse.pumas.ai/u/storopoli)\
**Post date:** [April 9, 2022, 11:40am UTC](https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665/2 "2022-04-09T11:40:36Z")

</div>

Could you please post your model and some simulated data for us to debug the error?

---

<div class="post-metadata">

**Author:** ![ahmed.salem](https://avatars.discourse-cdn.com/v4/letter/a/858c86/32.png) [@ahmed.salem](https://discourse.pumas.ai/u/ahmed.salem)\
**Post date:** [April 10, 2022, 7:42pm UTC](https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665/3 "2022-04-10T19:42:17Z")

</div>

Hi Dr. Storopoli,  
Thanks for the reply, Please find attached the code

```auto
using Pumas, MCMCChains

model_bayes = @model begin
  @param begin
    θ ~ Constrained(
      MvNormal(
        [3.0, 5.0],
        Diagonal(ones(2))
      ),
      lower = zeros(2),
      upper = fill(50.0, 2),
      init = [3.0,5.0])
    Ω ~ InverseWishart(6, diagm(fill(0.2, 2)) .* (6 + 3 + 1))
    σ ~ Gamma(1.0, 0.5)
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @pre begin
    CL = θ[1]*(WT/50)^0.75*exp(η[1])
    Vc = θ[2]*(WT/50)*exp(η[2])
  end

  @covariates WT

  @dynamics Central1

  @derived begin
    μ := @. Central / Vc
    CONC ~ @. Normal(μ, σ)
  end
end

df = CSV.read("./df_PUMAS.csv",DataFrame)

estimation = read_pumas(df,
    observations=[:CONC], 
    id=:ID,
    time=:TIME,
    evid=:EVID,
    mdv = :MDV,
    amt=:AMT,
    rate = :RATE,
    cmt = :CMT,
    covariates = [:WT]
    )

param = (θ = [3.0, 5.0], Ω = [1.0 0.0; 0.0 1.0], σ = 0.5)

result = fit(model_bayes, estimation, param, Pumas.BayesMCMC();
  nsamples=2000, nadapts=1000)

```

---

<div class="post-metadata">

**Author:** ![ahmed.salem](https://avatars.discourse-cdn.com/v4/letter/a/858c86/32.png) [@ahmed.salem](https://discourse.pumas.ai/u/ahmed.salem)\
**Post date:** [April 10, 2022, 7:46pm UTC](https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665/4 "2022-04-10T19:46:21Z")

</div>

and the dataset  
[dataset](https://docs.google.com/spreadsheets/d/1Dz9RGnogIYkzTBhN2VuAqfzKtLV83j8H7EG17_8EixY/edit?usp=sharing)

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/mohamed82008/32/127_2.png) [@mohamed82008](https://discourse.pumas.ai/u/mohamed82008)\
**Post date:** [April 11, 2022, 7:54pm UTC](https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665/5 "2022-04-11T19:54:08Z")

</div>

@ahmed.salem it’s not too easy to get the individual random effects now. I would wait until the next release where we have many more Bayesian features coming out that let you access samples from the posterior or posterior predictive distribution of any quantity in the model.

---

<div class="post-metadata">

**Author:** ![ahmed.salem](https://avatars.discourse-cdn.com/v4/letter/a/858c86/32.png) [@ahmed.salem](https://discourse.pumas.ai/u/ahmed.salem)\
**Post date:** [April 23, 2022, 11:38pm UTC](https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665/6 "2022-04-23T23:38:07Z")

</div>

Hi @mohamed82008 , Thanks for the reply. I have a followup question. In the fit output the following information is reported

```auto
Iterations = 1:2000
Thinning interval = 1
Chains = 1
Samples per chain = 2000

```

My question is about how to specify how many chains I would like to use (say for example I want to run 4 chains). How can I specify this? Is there an argument inside fit function to define the chain or I need to do certain maneuver? . Same also for thinning.

I am also a little bit unfamiliar with nadapts argument… is this similar to the burn in phase ? if no, how can I specify the number of samples to be in the burn in phase

Finally, I am getting warning messages like these during the fitting process

```auto
┌ Warning: The current proposal will be rejected due to numerical error(s).
│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /builds/PumasAI/PumasSystemImages-jl/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:47
┌ Warning: The current proposal will be rejected due to numerical error(s).
│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /builds/PumasAI/PumasSystemImages-jl/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:47
┌ Warning: The current proposal will be rejected due to numerical error(s).
│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /builds/PumasAI/PumasSystemImages-jl/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:47
┌ Warning: timestamp of type Missing unknown
└ @ MCMCChains /builds/PumasAI/PumasSystemImages-jl/.julia/packages/MCMCChains/VwH27/src/chains.jl:364
┌ Warning: timestamp of type Missing unknown
└ @ MCMCChains /builds/PumasAI/PumasSystemImages-jl/.julia/packages/MCMCChains/VwH27/src/chains.jl:364
Chains MCMC chain (2000×6×1 Array{Float64, 3}):

```

Do these warnings indicate a problem with the model specification ?

Thanks!

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/mohamed82008/32/127_2.png) [@mohamed82008](https://discourse.pumas.ai/u/mohamed82008)\
**Post date:** [April 24, 2022, 2:37am UTC](https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665/7 "2022-04-24T02:37:35Z")

</div>

> how to specify how many chains I would like to use (say for example I want to run 4 chains).

There is an `nchains` keyword argument you can specify in the `fit` function. In the next release, it will be 4 by default.

> I am also a little bit unfamiliar with nadapts argument… is this similar to the burn in phase ?

No, `nadapts` is different from burn-in. In Pumas, we use the No-U-Turn Sampler (NUTS) algorithm which is a variant of the Hamiltonian Monte Carlo sampling algorithm that’s more robust and was popularised by Stan. Our NUTS implementation and the Stan implementation are pretty much identical. The way the NUTS algorithm works involves 2 phases: a sampling + adaptation phase and a sampling with no adaptation phase. In the first phase, the algorithm tries to fine tune its parameters, e.g. the so-called step size and mass matrix, to achieve a certain target acceptance ratio. Once it’s done fine tuning the parameters, adaptation stops. `nadapts` describes how many of the `nsamples` iterations are going to be spent doing adaptation next to sampling. You can choose to get rid of all the adaptation steps as “burn-in” but that might be too much to throw away.

> if no, how can I specify the number of samples to be in the burn in phase

In the next release, there is a function that lets you do burn-in and thinning as a post-processing step.

> Finally, I am getting warning messages like these during the fitting process  
> Do these warnings indicate a problem with the model specification ?

The timestamp warning is fixed in the next release. The numerical error warning is not a big deal if it only shows up a few times. Basically, during sampling sometimes the sampler can wander off to parameter values that either cause the differential equation solver to not converge, or cause the likelihood to be 0 (log likelihood = -Inf). In these cases, you get numerical errors and the NUTS sampler will “reject” this problematic parameter value. These numerical errors can more commonly happen at the first part of the sampling rather than the later parts of the sampling for well-defined models. Sometimes good initialisation can also help avoid these numerical errors. If you are getting too many such warnings, then this is a sign your model “blows up” too often so you can try restricting the parameter space by using stronger priors centred around the parameter values that are known to “work”, as in they have non-zero likelihoods and the ODE solver converges. Alternatively, start with a simpler model and stat building the complexity up step by step until sampling breaks and then try to understand why it breaks.

---

<div class="post-metadata">

**Author:** ![ahmed.salem](https://avatars.discourse-cdn.com/v4/letter/a/858c86/32.png) [@ahmed.salem](https://discourse.pumas.ai/u/ahmed.salem)\
**Post date:** [April 24, 2022, 2:51am UTC](https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665/8 "2022-04-24T02:51:09Z")

</div>

Great thanks for the detailed reply! The warning is happening just two to three times at the beginning of the sampling process and the parameter estimates makes so I do not think it is a problem.

---

<div class="post-metadata">

**Author:** ![ahmed.salem](https://avatars.discourse-cdn.com/v4/letter/a/858c86/32.png) [@ahmed.salem](https://discourse.pumas.ai/u/ahmed.salem)\
**Post date:** [May 18, 2022, 10:10am UTC](https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665/9 "2022-05-18T10:10:03Z")

</div>

Hi @mohamed82008 ,

I currently have access to Pumas v2.1. Could you please let me what functions to use to get the Bayesian individual posteriors?

Thanks

---

<div class="post-metadata">

**Author:** ![storopoli](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/storopoli/32/190_2.png) [@storopoli](https://discourse.pumas.ai/u/storopoli)\
**Post date:** [May 18, 2022, 11:59am UTC](https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665/10 "2022-05-18T11:59:48Z")

</div>

@ahmed.salem, I think you can find more at the docs, [Estimating Parameters of Pumas Models · Pumas](https://docs.pumas.ai/dev/basics/estimation/#BayesMCMCResults).

You would probably need to convert the output to a `MCMCChains.jl` output, e.g.

```julia
my_chain = Chains(fit_result.chain)

```

---

<div class="post-metadata">

**Author:** ![ahmed.salem](https://avatars.discourse-cdn.com/v4/letter/a/858c86/32.png) [@ahmed.salem](https://discourse.pumas.ai/u/ahmed.salem)\
**Post date:** [May 18, 2022, 12:46pm UTC](https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665/11 "2022-05-18T12:46:11Z")

</div>

> [@storopoli](#):
>
> `my_chain = Chains(fit_result.chain)`

Hi @storopoli  
I tried the function you provided but it gives the following error

```auto
julia> my_chain = Chains(result_1.chain)
ERROR: type BayesMCMCResults has no field chain
Stacktrace:
 [1] getproperty(x::Pumas.BayesMCMCResults{Vector{Pumas.ThreadedBayesLogDensity{PumasModel{(θ = 2, Ω = 3, σ = 1), 2, ParamSet{NamedTuple{(:θ, :Ω, :σ), Tuple{Constrained{FullNormal, VectorDomain{Vector{Float64}, Vector{Int64}, Vector{Float64}}}, InverseWishart{Float64, PDMats.PDMat{Float64, Matrix{Float64}}}, Gamma{Float64}}}}, var"#5#17", var"#6#18", var"#8#20", var"#10#22", Central1, var"#11#23", var"#14#26", Nothing}, NamedTuple{(:θ, :Ω, :σ), Tuple{Vector{Float64}, Matrix{Float64}, Float64}}, Vector{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, T2, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}} where T2}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Vector{ForwardDiff.GradientConfig{ForwardDiff.Tag{Pumas.Tag1, Float64}, Float64, 8, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Pumas.Tag1, Float64}, Float64, 8}}}}, Vector{DiffResults.MutableDiffResult{1, Float64, Tuple{Vector{Float64}}}}, NamedTuple{(), Tuple{}}}}, Vector{Vector{Vector{Float64}}}, Vector{AbstractMCMC.SamplingStats}}, f::Symbol)
   @ Base ./Base.jl:42
 [2] top-level scope
   @ ~/data/code/BAYESIAN_PUMAS/Code.jl:60

```

Also, I think the document has not been updated yet to include information about getting individual parameters or individual posteriors

---

<div class="post-metadata">

**Author:** ![ahmed.salem](https://avatars.discourse-cdn.com/v4/letter/a/858c86/32.png) [@ahmed.salem](https://discourse.pumas.ai/u/ahmed.salem)\
**Post date:** [May 23, 2022, 11:34am UTC](https://discourse.pumas.ai/t/individual-pk-parameters-following-bayesian-estimation/665/12 "2022-05-23T11:34:53Z")

</div>

Hi @mohamed82008, I am following up on my last response since I have access to Pumas 2.1 now. Could you please let me know the new feature that allow accessing the individual posteriors of PKPD parameters of each subject included in the analysis?