# Bayesian pharmacometrics modeling with Stan and Torsten preprint

**URL:** <https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537>\
**Category:** Events\
**Created:** [September 22, 2021, 7:22pm UTC](https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537 "2021-09-22T19:22:43Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [September 22, 2021, 7:22pm UTC](https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537/1 "2021-09-22T19:22:43Z")

</div>

Dear Stan community,  
@yizhang , @billg, and I released a preprint this week demonstrating the use of Stan and Metrum’s homebrewed extension, Torsten, for pharmacometrics modeling.

> **[Flexible and efficient Bayesian pharmacometrics modeling using Stan and...](https://arxiv.org/abs/2109.10184)**
>
> Stan is an open-source probabilistic programing language, primarily designed
> to do Bayesian data analysis. Its main inference algorithm is an adaptive
> Hamiltonian Monte Carlo sampler, supported by state of the art gradient
> computation. Stan's...

This is the first part of a two-parts tutorial, targeted at pharmacometricians. It includes detailed examples and code. We wrote this paper as Torsten and Stan developers, and also instructors who have taught workshops on Bayesian pharmacometrics modeling: we devoted a lot of effort to clarifying various Bayesian concepts and emphasize some of the more subtle points which sometimes get lost. Part II of the tutorial is in progress and will demonstrate more advanced features of Stan and Torsten.

I’m always excited to write about Torsten because this is the project that introduced me to Bayesian statistics and ODE-based models!

---

<div class="post-metadata">

**Author:** ![Funko\_Unko](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/funko_unko/32/10786_2.png) [@Funko\_Unko](https://discourse.mc-stan.org/u/Funko_Unko)\
**Post date:** [September 23, 2021, 3:13pm UTC](https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537/2 "2021-09-23T15:13:32Z")

</div>

So how would eg for the last nonlinear example compare NONMEM, Torstan and Stan without Torsten? Let’s say in terms of usability, accuracy of the uncertainty quantification and efficiency?

---

<div class="post-metadata">

**Author:** ![yizhang](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/yizhang/32/15673_2.png) [@yizhang](https://discourse.mc-stan.org/u/yizhang)\
**Post date:** [September 23, 2021, 5:51pm UTC](https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537/3 "2021-09-23T17:51:47Z")

</div>

Good question.

- The stan/torsten vs nonmem benchmark is in progress. On the hand, Stan’s rich diagnostics is definitely missing in nonmem’s BAYES.
- The last example is actually easier to code using just Stan, because it involves only a single patient. It’s there as a footstep to the part 2 of the tutorial, where population models with various likelihood structures will be explored. Generally when it comes to population model, coding various dosing events becomes less pleasant in Stan, and that’s one of the major motivations of torsten.
- We choose to follow nonmem events schedule nomenclature and this has both
  - pros
    - it’s friendly to PMX practitioners due to nonmem prevalence.
    - it’s easier to be integrated into existing workflows

  - and cons:
    - ODE definition must comply with old Stan style.
    - function signature could be complicated (but we are working on simplified versions)

---

<div class="post-metadata">

**Author:** ![yizhang](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/yizhang/32/15673_2.png) [@yizhang](https://discourse.mc-stan.org/u/yizhang)\
**Post date:** [September 23, 2021, 6:44pm UTC](https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537/4 "2021-09-23T18:44:36Z")

</div>

Another issue of using nonmem-style events is that when some of the events are Stan `parameters` we have to define the entire array as parameter, e.g on p.26 of the paper

```stan
real<lower = 0> parms[9];
parms = {CL , Q, V1, V2, ka, mtt , circ0 , gamma , alpha };

```

becomes inefficient when, say only clearance `CL` is a parameter. This can be resolved when tuple becomes available.

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [September 23, 2021, 9:56pm UTC](https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537/5 "2021-09-23T21:56:50Z")

</div>

> [@yizhang](#):
>
> The last example is actually easier to code using just Stan, because it involves only a single patient.

Even though there’s only a single patient, you still need to handle multiple doses. What’s more, Torsten supports a mixed solver, meaning you can solve the PK analytically, while numerically solving the remaining PD equations. Coding this directly in Stan is possible, but long and error prone. Not only do you need to write the analytical solution to the PK system, it also takes a bit of book-keeping to then write the ODE which gets solved numerically.

---

<div class="post-metadata">

**Author:** ![Funko\_Unko](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/funko_unko/32/10786_2.png) [@Funko\_Unko](https://discourse.mc-stan.org/u/Funko_Unko)\
**Post date:** [September 24, 2021, 8:24am UTC](https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537/6 "2021-09-24T08:24:59Z")

</div>

Thank you both @charlesm93 and @yizhang for responding. I can’t wait for part 2!

---

<div class="post-metadata">

**Author:** ![Funko\_Unko](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/funko_unko/32/10786_2.png) [@Funko\_Unko](https://discourse.mc-stan.org/u/Funko_Unko)\
**Post date:** [October 7, 2021, 6:27am UTC](https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537/7 "2021-10-07T06:27:03Z")

</div>

So inspired by [this thread on the adjoint ODE solver](https://discourse.mc-stan.org/t/adjoint-ode-behaviour/24703), what are your experiences with the adjoint solver? Or with the BDF solver? BTW, out of curiosity I’ve reimplemented the nonlinear model from section 4 and it looks like this implementation using Stan’s regular rk45 is quite a bit (~30%) faster than the specialized `pmx_solve_twocpt_rk45`, is this to be expected? Also, don’t you guys have some kind of adaptive (stiff-\>non-stiff) and parallel (combining draws across chains for cov-estimation) warm-up. How much faster does this make the problem?

---

<div class="post-metadata">

**Author:** ![yizhang](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/yizhang/32/15673_2.png) [@yizhang](https://discourse.mc-stan.org/u/yizhang)\
**Post date:** [October 7, 2021, 3:41pm UTC](https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537/8 "2021-10-07T15:41:59Z")

</div>

I haven’t used the adj solver, BDF on the other hand on a regular basis. It’d be great to see your implementation as I do notice that sometimes torsten’s version of rk45 performs differently for better or worse. We’ve discussed stiffness/nonstiffness detection solvers but in general the technique is both involved and not very robust. Torsten does have cross-chain implemented and we plan to discuss that in the coming part 2. Recently I did a benchmark for the ACoP12 conference and can post it here next month.

---

<div class="post-metadata">

**Author:** ![Funko\_Unko](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/funko_unko/32/10786_2.png) [@Funko\_Unko](https://discourse.mc-stan.org/u/Funko_Unko)\
**Post date:** [October 7, 2021, 7:14pm UTC](https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537/9 "2021-10-07T19:14:58Z")

</div>

I would actually have thought that the BDF solver could work better for the linked nonlinear problem, as the stepsize of rk45 is bounded by the first compartment due to stability concerns and for the other compartments not much is really going on. But I tried to drop it in, and it just did nothing… 🤷

---

<div class="post-metadata">

**Author:** ![arya](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/arya/32/1626_2.png) [@arya](https://discourse.mc-stan.org/u/arya)\
**Post date:** [October 22, 2021, 12:34am UTC](https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537/10 "2021-10-22T00:34:20Z")

</div>

Hey guys cool paper! Looking forward to part two. I’ve been using Torsten every day as of late. I had a couple of questions.

1. Does Torsten currently support the adjoint solver in Stan?

2. Are the group integrator functions wrappers for parallelization? Is parallelization only over a cluster or is it possible to parallelize across cores in a multicore machine? I have access to a 48-core machine and my naive plan was to just wrapping the `pmx_solve_bdf()` function so that I can call it using Stan’s `reduce_sum` functions. Would that be possible in Torsten?

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [October 1, 2022, 6:12pm UTC](https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537/11 "2022-10-01T18:12:12Z")

</div>

Hey @arya,  
Sorry it took me so long to get back to you on this, this slipped through the cracks.

> [@arya](#):
>
> - Does Torsten currently support the adjoint solver in Stan?

No (unless this was recently added…). As far as I can tell, the size of typical PKPD ODEs doesn’t warrant using the adjoint method but maybe you have some counterexamples in mind. If so, we can open an issue on GitHub. The algebraic solvers, on the other hand, all use adjoint methods, and that’s something you get in Torsten.

> [@arya](#):
>
> Are the group integrator functions wrappers for parallelization? Is parallelization only over a cluster or is it possible to parallelize across cores in a multicore machine?

You can certainly use `reduce_sum`, but Torsten supports a population ODE solver, which handles the parallelization more seamlessly (see [General ODE-based Population Model Function - Torsten](https://metrumresearchgroup.github.io/Torsten/function/ode-pop/)). So my recommendation is to use the `pmx_solver_group` function.

The upcoming part II of our Torsten Tutorial should cover these topics ([Part I](https://ascpt.onlinelibrary.wiley.com/doi/10.1002/psp4.12812) just appeared in CPT: pharmacometrics)

@yizhang can probably provide additional insights.

---

<div class="post-metadata">

**Author:** ![arya](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/arya/32/1626_2.png) [@arya](https://discourse.mc-stan.org/u/arya)\
**Post date:** [October 1, 2022, 11:32pm UTC](https://discourse.mc-stan.org/t/bayesian-pharmacometrics-modeling-with-stan-and-torsten-preprint/24537/12 "2022-10-01T23:32:42Z")

</div>

Thanks @charlesm93! Yeah I’ve since played around with both `pmx_solver_group` and MPI and also just wrapping `pmx_solve_twocpt` in a `reduce_sum`. You’re right that the former technique is seamless in terms of cleaner Stan code, but the latter can be useful at times. For example, I ran into a bug that only shows up in the MPI version, but I made an issue on Github. Also the MPI installation can be a little daunting, although I do have a Bash script I can share that I use to install it now on fresh Metworx instances.