Bill Gillespie, of Metrum, is giving a tutorial next week at ACoP:

- Getting Started with Bayesian PK/PD Modeling Using Stan: Practical use of Stan and R for PK/PD applications

Thursday 8 October 2015, 8 AM — 5 PM, Crystal City, VA

This is super cool for us, because Bill’s not one of our core developers and has created this tutorial without the core development team’s help. Having said that, we’ve learned a lot from Bill and colleagues on our mailing lists as we were designing ODE solvers for Stan (an ongoing issue—see below for future plans).

Bill’s tutorial is up against a 2-day Monolix tutorial and a 2-day tutorial on R by Devin Pastoor, who’s also been active on our mailing lists recently.

**Why Stan for PK/PD?**

In case you’re wondering why people would use Stan for this instead of something more specialized like Monolix or NONMEM, it’s because of the modeling flexiblity provided by the Stan language and the effectiveness of NUTS for MCMC. So far, though, we’re in the hole in not having a stiff ODE solver in place. Or a good NONMEM-like event data language on top.

Maybe Bill will jump in with some other motivations.

**What’s in Store for Stan’s ODE Solvers?**

There’s been lots of behind-the-scenes activity on our ODE solvers—we’re really just getting ~~burned in~~ warmed up.

The next minor release of Stan (2.9) should stop the freezing issue when parameters wander into regions of parameter space that lead to stiff ODEs. And we’ve really sped up the Jacobian calculations when Michael Betancourt realized we were doing a lot of redundant calculation and he and I put a patch in to fix it. We should also allow user-defined control of absolute and relative tolerances.

Next, hopefully by Stan 2.10, we’ll have a stiff solver and maybe a way for users to supply analytic coupled-system gradients and Jacobians. Stay tuned. These new designs are largely being guided by Sebastian Weber and Wenping Wang at Novartis. And of course, by Michael Betancourt working out all the math and Daniel, Michael, and I working out the code with Sebastian’s and Wenping’s input.

We also need to evaluate how well variational inference works for ODE problems. Our early trials are very promising. Then we could replace the max marginal likelihood approach of NONMEM with a very speedy variational inference mechanism allowing much more general models.

There’s more in the works, but the above are the top of our to-do list.

Nice post, Bob, and an excellent illustration of why we prefer “warm up” to “burn in”!

P.S. And we have plans to implement maximum marginal likelihood in Stan too.

I can imagine variational inference working well for many ODEs. Often the probabilistic error model in an ODE system is just “plus wiggle room because we know it’s not going to go through all the data points” and if you can find a reasonable and fast to calculate approximation to the already approximate “plus wiggle room” you are doing well.

Hi Bob,

Thanks for the plug.

As you say the primary motivations for using Stan for pharmacometric applications fall under the headings of flexibility and the effectiveness/efficiency of the sampler. In my workshop I have a slide headed “Why Stan?” with the following response:

Flexibility

* Flexible w.r.t. stochastic structure

** Any number of levels variability

** Large selection of built-in probability distributions

** Permits sub-models with very different stochastic hierarchies

* Flexible w.r.t. deterministic structure

** Control structures: if-then-else, for and while loops

** Large collection of built-in functions

** Operators and functions for vector and matrix calculations

* Computational efficiency

** Often faster than Gibbs or Metropolis-Hastings

*** Measured in terms of time / (effective sample size)

* Also includes optimization and variational inference methods for rapid approximate Bayesian analysis

For example when comparing Stan to BUGS variants on one hand and more specialized pharmacometric modeling tools like NONMEM and Monolix on the other, I ask the question “What can you do with Stan that you can’t do with those tools. My answer is models containing submodels with different stochastic structures and complex deterministic structures. BUGS can do the first but not the second and vice versa for NONMEM and Monolix.

Stan also distinguishes itself as freely available, open source software and by the very active ongoing development program that means there are only better things to come.

It’s great to hear that ODE solvers in Stan suitable for pharmacometric applications are near at hand. At Metrum we are developing Stan functionality for dealing with event schedules, e.g., dosing, that may be coupled with such solvers. That combination will provide functionality similar to the PREDPP component of NONMEM. Between us or others like Sebastian and Wenping working along similar lines, the pharmacometric user community will likely see Stan become a logical alternative to NONMEM in a reasonably short time.

Thanks,

Bill

I don’t work in Pharmacokinetics but Chemical Kinetics & almost every problem I work on has stiff ODEs. Simply because rate constants span such a large range of magnitudes in any system of pathways.

Is the situation similar in PK / PD? Or is there work in PK/PD that doesn’t encounter stiffness?

Out of curiosity, how well do things like regularized collocation by interpolator functions work? Example;

you have a spring/mass system with high spring stiffness of one spring and low stiffness of another so that you have two widely separated timescales. You’re interested in the long time scale behavior but the fast oscillations of the short timescale behavior make long timesteps impossible. Instead, you approximate the displacement behavior by two radial basis functions with certain constraints on the centers and smoothness parameters. You then turn the differential equation into an algebraic equation requiring the approximate solution (which is now an analytic RBF) to satisfy the differential equation at a finite set of points within the time domain of interest. The regularization of the RBF prevents it from being able to “wiggle” in a way that fits the fast timescale behavior.

I could imagine this kind of technique could be appropriately applied to such systems, though I’ve more heard it applied in the spatial domain of PDEs with the time-stepping being more of a standard time-stepping algorithm.

Question also for Stan guys: does Stan have a root finding algorithm that could be used to solve nonlinear algebraic equations? How about a linear solver that can be used to efficiently solve linear systems without calculating the inverse? That would be quite useful for utilizing algebraic model structure in defining the posterior. trivial nonlinear example:

a,b are parameters, x,y are data, k a local variable, I am making this up out of whole cloth for illustration purposes only:

a ~ normal(0,1);

b ~ normal(0,1);

solve(x-(2*k*a^2-b^2/k) = 0,k);

y ~ normal(k,1);

No, we don’t have any tools for efficient root finding yet. But it’s on our longer-term to-do list.

Is it necessary to be able to differentiate this kind of thing wrt. each parameter for Stan to work? I’d imagine that would make root-finding algorithms quite challenging.

@Daniel

I don’t know about the implementation details much. I use the Stiff ODE Solvers in Matlab or the Sundials Toolbox and those seem to work fine for me.

But try a non-stiff solver on the same problem and it usually fails miserably.

There’s Stochastic Solvers like the Gillispie based ones that can also be used but even they have some serious problems when they hit a stiff system.

Yes, many of them are one- or two-compartment models with very limited stiffness. The problem we’ve been having is that the systems aren’t stiff in the neighborhood of high posterior mass, but can be stiff out in the tail.

How much do informative priors help in your example models? Are the models expressed in dimensionless form? The choice of scaling can help a lot in my experience.

Rahul, have you seen GNU MCSim (https://www.gnu.org/software/mcsim/)? That’s a much earlier MCMC system developed by Frederic Bois and Don Maszle under Andrew’s guidance, as I understand it. Frederic uses and supports it today. MCSim is based on MH, not HMC, but it can handle stiff ODEs and events. Frederic has consulted on the ODE solver approach in Stan, too.

Bill:

Yes, much of MCSim’s capabilities now exist in Stan (and of course Stan can do lots and lots of things that MCSim can’t). We’re working with Frederic and others to get stiff solvers working in Stan too.

I gave up after making Stan freeze with some pretty trivial systems, glad to hear there’s some development happening.

Maybe if you posted your cases / code it’d help the developers improve on the bugs?

+1 on that. We love getting these reports on the Stan users list. It’s not always bugs in Stan or in user code, it can also be models that are not buggy but can be written more efficiently, or models that Stan an’t currently fit efficiently, in which case we’d like to know. Of course if Stan is actually freezing the computer, then, yes, that’s a bug!

This is a great post. The link to ACoP tutorial is broken though. I am looking for example Stan code for the PK/PD models. I can’t seem to find a good example with accompanying code anywhere.