stiff ODEs and population PK models

[Linas Mockus]

Hi,

Are there any improvements in stiff ODE solver in 2.12 as compared to 2.11?

Does anybody have experience with population PK models (essentially hierarchical models) that are based on system of ODEs? My specific question is whether to assume that PK parameters for each individual are distributed normally around corresponding population parameter - individual_parameter~N(population parameter,sd) - or it is better to assume log normal - log(individual_parameter)~N(log(population parameter),sd) - or something else? In addition I would love to hear comments on PBPK (physiologically based PK) population models but any comments on traditional compartmental PK population models would be a good starting point.

Thanks,
Linas

[Bob Carpenter]

Not sure what you're looking for. We fit these all the time,
including lots of hierarchical models. The hierarchical model
can usually be done two ways---directly on whatever constrained
scale you have with priors that match the constraints or on
an unconstrained scale with transforms to the constrained scale.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[lmockus]

Bob,

For some reason I thought that posterior estimates are different. Maybe
I was wrong.

As I understand NONMEM uses log transform and the log of individual
volume of distribution (Vd) is distributed normally around the
population Vd (as an example). However I guess I saw examples when Vd of
individual~N(Vd of population). I wander what is best practice and what
is most efficient in stan.

Linas

>You received this message because you are subscribed to a topic in the

>Google Groups "Stan users mailing list" group.

>To unsubscribe from this topic, visit
>https://groups.google.com/d/topic/stan-users/ghKrcnKd6Hg/unsubscribe.
>To unsubscribe from this group and all its topics, send an email to

[Bob Carpenter]

You want to look at the noise in the model itself.

For instance, in the soil carbon models, respirated carbon,
which is positive constrained, tends to have lognormal measurement
noise (noise is proportional to measurement rather than being
a constant scale additive term).

There's not a huge difference in efficiency between declaring
a positive value and giving it a lognormal distribution and
declaring an unconstrained value, giving it a normal distribution,
and then exp() transforming it to positive.

- Bob

[Linas Mockus]

Bob,

Thanks, I think that PK models also have lognormal measurement noise.

Another question I originally had is are there any improvements with respect to solving ODEs in stan cmd 2.12 vs 2.11. Currently I am generating jacobian using dmod and code libraries but I am not sure if the process of linking jacobians to executable is the same.

Linas

[Daniel Lee]

Whatever you're doing for CmdStan 2.11 should work for 2.12 and now 2.13. The CmdStan code base hasn't changed much.

Daniel

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Devin Pastoor]

Linus,

Not to be nitpicky but be careful distinguishing 'measurement' noise - which you might attribute to assay error, which will be normally distributed (your $SIGMA) in nonmem. The lognormal distribution of PK parameters is driven mostly by physiology. Given the etas are normally distributed, you generate the lognormal parameter distribution via

TVP*EXP(nP) 

And FYI, for bayesian models, nonmem will use a normal distribution for variance around thetas and an inverse wichtart distribution for omegas, if you want to match that with STAN.

[Bill Gillespie]

Hi Linas,

You may want to look at the work we're doing to develop a suite of Stan functions for pharmacometric applications, particularly PKPD models. Here's a couple relevant links:

• https://github.com/charlesm93/example-models/blob/feature/issue-70-PKPDexamples-torsten/PKPD/torsten/torstenManual.pdf
  
• https://github.com/charlesm93/example-models/blob/feature/issue-70-PKPDexamples-torsten/PKPD/torsten/README.md
  

The user manual provides some popPK and popPKPD examples. The code and data for the examples is also available in the repo.

We usually describe inter-individual variation in the model parameters as lognormal, though usually coded in terms of a multivariate normal distribution of log transformed parameters. Of course this is because the parameters typically make sense only if they are positive -- but there are sometimes exceptions.

Cheers,

Bill

[Bob Carpenter]

We have our own recommendations for covariance. Many
systems were stuck with inverse Wisharts for covariance
(or Wisharts for precision matrices) because of the
conjugacy. We prefer scaled LKJ priors as they make more
sense from a modeling perspective in most cases.

- Bob

[Devin Pastoor]

Hi Bob,

Do you have any recommendations reading-wise, or at least pointers in the direction, to back up the scaled LKJ priors suggestion? I have a feeling that one of the first questions that will be posed is why not wishart, and I have a feeling "because people smarter than me recommended it" won't cut it :-)

[Bob Carpenter]

Check out this to develop some intuitions:

http://www.stat.columbia.edu/~gelman/research/unpublished/Visualization.pdf

What do you tell people now when they ask you to justify an
(inverse) Wishart?

What the inverse Wishart has going for it is that it's the conjugate
prior for covariance matrices in a multivariate normal model. As
such, it can be represented as prior multivariate normal data.

The LKJ prior is like a beta prior in some ways. The
bigger the parameter, the more concentrated the prior is around the
unit correlation matrix. So it lets you directly parameterize
correlation.

Then you can control the scales any way you want to allow for both long
and short-tailed behavior or zero-avoiding behavior (e.g., Cauchy, normal,
lognormal).

In the 2D case this is like putting a beta prior stretched from (0,1)
to (-1, 1) on the correlation rho, and putting an independent
prior on the scales sigma[1] and sigma[2]:

Sigma = sigma[1] * sigma[1] rho * sigma[1] * sigma[2]
rho * sigma[1] * sigma[2] sigma[2] * sigma[2]

You see people do this all the time in the 2D case and the LKJ
prior is just a natural generalization to higher dimensions.

To gain some degree of separation between the scale priors and
the correlation priors, Gelman and Hill hacked something similar
in their regression book by scaling a covariance matrix with an
inverse Wishart prior. It's much easier to get the effect they were
after with a scaled LKJ prior.

- Bob

[Devin Pastoor]

awesome. Thanks for taking the time to respond. 

To answer your question - people don't ask at this point. Mainly because of the contraints of what NONMEM gives you (basically only inverse wishart unless you want to jump through a lot of hoops).

Also, other than the very small (passionate) minority, the additional hurdles bayesian analysis poses to the status quo of model development makes it so very few people are developing models in that fashion. It's really only under situations where either

a) you a priori know you want to make quantification and handling of uncertainity 'first class'

b) all other options have failed

are people turning to bayesian techniques.

For example, if you look at the ddmore model repository of ~200 models http://repository.ddmore.eu/ as well as my personal interactions with the model bases at Merck, the FDA, and our own group pushing that number to probably 300+ models, I think less than 5% of the models used full bayes, especially for PK models. You see it more when doing like model based meta analysis or some systems biology/PBPK modeling though, and there only metrum and some select individuals seem to make an consistent noise in the space.

Anyway - I don't want to turn this into an exposition, however if you are interested in the why of lack of full bayes (much of which I bet you can guess at) or what stan can (continue) doing to make a dent in that, I'm happy to expand on the state of things.

Thanks again for the pointers!

Devin

[Bob Carpenter]

> On Dec 19, 2016, at 5:30 PM, Devin Pastoor <devin....@gmail.com> wrote:
>
> awesome. Thanks for taking the time to respond.
>
> To answer your question - people don't ask at this point. Mainly because of the contraints of what NONMEM gives you (basically only inverse wishart unless you want to jump through a lot of hoops).

I figured that'd be the answer.

> Also, other than the very small (passionate) minority, the additional hurdles bayesian analysis poses to the status quo of model development

You mean in publishing? Because for the most part, Bayesian
modeling is easier because posterior means are so much better
behaved than posterior modes (e.g., we don't have to limit models
like NONMEM [or lme4 does] in order to fit them with max marginal
likelihood).

> makes it so very few people are developing models in that fashion.

We must know them all :-)

> It's really only under situations where either
>
> a) you a priori know you want to make quantification and handling of uncertainity 'first class'
>
> b) all other options have failed
>
> are people turning to bayesian techniques.

I think as people get more sophisticated in thinking about uncertainty,
(a) will be come a priority, and as they want to build more
sophisticated models, (b) is critical. We've come into projects largely
because of (b) --- people had models to fit that they couldn't fit
with NONMEM or other off-the-shelf tools.

> For example, if you look at the ddmore model repository of ~200 models http://repository.ddmore.eu/ as well as my personal interactions with the model bases at Merck, the FDA, and our own group pushing that number to probably 300+ models, I think less than 5% of the models used full bayes, especially for PK models. You see it more when doing like model based meta analysis or some systems biology/PBPK modeling though, and there only metrum and some select individuals seem to make an consistent noise in the space.

We have an ONR grant to make PK/PD easier with Metrum. They're
working on NONMEM-style interfaces for "events" (dosing and measurement)
and on better solvers (just added a matrix exponential).

> Anyway - I don't want to turn this into an exposition, however if you are interested in the why of lack of full bayes (much of which I bet you can guess at) or what stan can (continue) doing to make a dent in that, I'm happy to expand on the state of things.

I don't care so much about why people don't do full Bayes, but
I am curious about what you think Stan can do.

- Bob

[Devin Pastoor]

> You mean in publishing?  Because for the most part, Bayesian 

modeling is easier because posterior means are so much better 
behaved than posterior modes (e.g., we don't have to limit models 
like NONMEM [or lme4 does] in order to fit them with max marginal 
likelihood). 

Hurdles namely in terms of implementation details (eg. I know i want a 2 compartment model with linear absorption with covariates of weight and gender, but how do I code that and then look at model diagnostics), and from an acceptance perspective - 'matching' up to historical comparisons. Eg, its very easy at this point to say, model development was performed using step-wise covariate modeling using a forward inclusion criteria of p < 0.05 and backwards step of p < 0.001. Here are x, y, z diagnostic plots. You can point to 20 other publications with similar diagnostics and interpretations. For example, the 2 primary plots shown to demonstrate model adequacy is conditional weighted residuals  and observed vs individual predicted concentrations, both shown here: http://devinpastoor.com/simplest_bayes/foce-analysis.html#visual-summaries 

These frequentist diagnostics by and large go out the window for bayesian analysis. I don't say this as a good or bad thing, just that the current status quo heavily favors coloring by those numbers to provide historical backing to the modeling strategy used.

> I don't care so much about why people don't do full Bayes, but I am curious about what you think Stan can do. 

I believe they are reasonably tightly coupled in a lot of situations, somewhat touched on above.

First, for the remaining discussion we can set aside models that only stan/winbugs/etc can handle, and instead focus on what can be done such that people who would normally have turned to nonmem would instead grab Stan. 

1) Push out some more (simple, but holistic) examples. 

This is going to be the largest hurdle once the ODE implementation settles. For better or worse, the majority of people will want to 'copy-paste' their way through scenarios they are familiar with to wrap their head around. Many companies have model libraries of 100's of models with 'tricks' on how to implement things. 

For example, I'm picking away at an extremely simple example: http://devinpastoor.com/simplest_bayes/foce-analysis.html#visual-summaries and the associated code for the project: https://github.com/dpastoor/simplest_bayes

Such case studies would cover:

- what the data looks like

- model setup with some information

- diagnostic outputs with some minimal interpretation

I know documentation is so low on priorities for the core team, but without some more examples, people such as myself are stuck being unable to directly contribute because it can be difficult to even get off the ground. 2-3 different case-study-based examples would give enough overlapping surface area to provide a feedback loop of being able to generate more community-driven examples off those templates. Similarly, as I mentioned at the beginning, there needs to be at least some guidance in terms of interpretation. This might be comparing a false model to the true model and noting some of the things, or even just a note here or there about what one might look for. 

2) Provide some 'expert opinion' or commentary from Andrew or others that can allay some question to approaches to model selection. Namely companies often waste time debating over covariate selection. Giving a nod to reasonable approach(es) would be a huge boon. This is made more tricky by the current gold standard of stepwise covariate modeling via change in objective function values does not have a clear a bayesian equivalent?.  Metrum, and myself personally, are using a full model approach to essentially select for covariates that can be deemed clinically relevant based on their posterior distribution http://metrumrg.com/assets/pubs/GastonguayPAGE2011.pdf  

2.1) If accomplishing this requires pre/post-processing, providing some tooling around accomplishing this goal. The #1 way via MLE right now is via step-wise modeling, and much of that is due to it being 'good enough' statistically, and dead simple due to it being implemented as little more than a config file + base model structure http://psn.sourceforge.net/pdfdocs/scm_userguide.pdf I think the above full model approach diagnostics are almost completely implemented from a visually in jonah's bayesplot, but could use some additional helpers for continuous covariates (as shown on slide 27 above).

3) Stay the course with Metrum implementing an events system and general ODE implementation details. I imagine there will also need to be some wrappers/templating around other constructs that might not fall under that initial objective (lag time, bioavailability fractions, transit compartments)?

4) longer term (that I hesitate to even pose at this point due to how many more important fundamental constructs need to be implemented first) an additional abstraction layer to give more nonmem style control-file that removes some constructs like explicitly coding to loop through individuals. For example, you can see here how 1 cmt model with analytical solution is defined for a population. https://github.com/dpastoor/simplest_bayes/blob/master/models/one_cmt_iv.cpp vs the current stan implementation https://github.com/charlesm93/example-models/blob/feature/issue-70-PKPDexamples-torsten/PKPD/torsten/TwoCptModelPopulation/TwoCptModelPopulation.stan

By and large, it seems people running models in pharma seem to fall under opposite ends of the spectrum: 1) your formally trained statisticians, and 2) (myself included) have background in biological sciences, some math and have been picking up stats on the way, but virtually no formal training, including 'simple' things like digesting proofs. For those of us in camp #2, there is this huge gap between the baby step of here is bayes rule and why its awesome, to 'introductory' texts where the expectation is still at minimum undergrad courses in inference, probability, and general comfort digesting proofs and formula.   An intermediate ground of applied examples could allow people to more directly connect with the domain while also providing pointers to the theory such that they can digest as they become more comfortable. 

Hopefully not too long winded. 

And as a bonus, here, for example, might be an excellent case study:

Which dose should we choose to maximize the benefit/risk profile in a target population, given rich data from a small subset of individuals studied in a previous trial? Let's say we've done a trial and have data for ~100 individuals at 4 escalating dose levels, each with 4-6 samples. To answer this we will need to do the following 

1) Build and evaluate structural and residual error model(s) - is it one cmt/two cmt, and is the residual error additive, exponential, some combination

    - what are the diagnostic plots/calculated values used to compare such models

2) Determine the effect of covariates on model parameters to determine which to retain, as well as to get a feel for which will need to be most closely monitored for potential dose adjustment in later phases (or to kill the drug all together). 

    - how do you address covariate inclusion/exclusion and overall building

3) Simulate from the final model and compare against some criteria (such as % of individuals in a given exposure window). 

    - how to simulate under different designs

3.5) perform sensitivity analysis to determine the impact of being 'wrong' on a given parameter - say the population variability on clearance - what if it was 50% greater, lower, etc. This can be relevant when moving to different subpopulation, eg studied in healthy now going to sick patients with additional comorbidities that may effect drug clearance.

4) bonus: use the above model as a prior to fit some new data.

Hope this helps.

Devin

[Bob Carpenter]

Thanks. That's very much in line with what others have
told us. We are teaching such classes with such case
studies---Bill Gillespie of Metrum's taught a bunch and
we just taught a class in Paris this Fall and will probably
be doing more. We just need to make some of the case studies
public.

To think that this is how some pharmacologists choose covariates
is frightening. We tend to include all the covariates we have
and let the priors and hierarchical modeling sort them out.
Or we'll do that, then select out the ones that aren't doing
predictively after the model is fit (Aki Vehtari has some papers
on this approach).

Those 100 model PK/PD libraries can (and are) being built on
top of Stan. Our target audience is really the people writing those
libraries rather than the users of those libraries. That is,
it's another level on top of Stan, like the RStanArm package.
Hopefully this level of user knows more stats. What we'd like to
do is bring more people up to this level. That's why doc's so
important. The whole point of working directly in Stan (or BUGS
or JAGS or NONMEM) is that you can make more reasonable data
and modeling choices for your scientific process and measurement
process.

Just to clear up one misconception: doc is *not* a low
priority for us! It is, unfortunately, a very low priority
for most programs in most funding agencies.

- Bob

[Devin Pastoor]

Good points - and yeah there is a lot of baggage on covariate model building. There is an intense fear of 

multiple comparisons so techniques like scm + some reduction are in place. The other, historically, difficult part 

about letting things sort, is convergence issues with more complex models, especially to get standard errors. 

Thankfully, this is extremely natural in the bayesian context. 

I've pinged Charles and by proxy Bill to help out with the case study generation, so hopefully I can alleviate some

of the time for them on that end and let them focus on implementations.

Thanks again for the pointers. 

[lmockus]

It is a great discussion but I just wanted to chime in with my humble needs...

 

I did comparison with stan cmd 2.13.1 between the performance when jacobian of ODE is not supplied and when jacobian is supplied. Sebastian taught me the trick by using dMod and cOde. My ODE is 34 state and quite stiff. I have lot of parameters but just estimate few of them.

 

1. When jacobian is supplied diagnose gives this result. I am curious how to interpret it but it seems there is big difference between finite differences and model. Can somebody explain why?

TEST GRADIENT MODE

 

 Log probability=-17.7607

 

 param idx           value           model     finite diff           error
         0               0        -40.7104        -40.7104   -1.58568e-008
         1          1.1474        -2.88137       -0.280441        -2.60093
         2         5.39363        -6.77934        -8.89517         2.11583
         3         2.07944       -0.829541         1.06564        -1.89518

2. Variational inference breaks after few iterations because cvode error -4. It might be explained that var inference is not designed to handle such errors.

 

3. Although gradient evaluation for the case when jacobian is supplied is much less the total simulation time is pretty much the same. I am still running second simulations but previous results showed that total time is less when jacobian is not supplied. Strange? Can it be attributed to stiffness?

 

4. Even worse is with convergence. While both converge to the similar posterior the case when jacobian is not supplied shows excellent mixing (visual) while when jacobian is supplied it shows autocorrelation. I still didn't run several disperse chains but can it be explained somehow? Maybe stiffness.

 

Data is simulated so my expectation that supplying jacobian makes posterior estimates easy and fast. All files are attached to make example reproducible. Please note that when generating jacobian I just copy ....stan to ...jac.stan to be consistent with namespace.

 

Thanks for any insights,

Linas

[Sebastian Weber]

It's quite amazing what you did there... would it be possible to first simplify things to the minimum? Let's start with few a ODEs (a single one only if possible) and then build up the complexity gradually. Doing 34 ODEs straight is a tough call.

I hope I find some time to look into this again.

[lmockus]

Bob/Aki,

Could you please give references to papers mentioned in the second
paragraph below?

Thank you,

Linas

------ Original Message ------
From: "Bob Carpenter" <ca...@alias-i.com>
To: stan-...@googlegroups.com
Sent: 12/19/2016 11:04:53 PM
Subject: Re: [stan-users] stiff ODEs and population PK models

[Jacob Leander]

Really interesting Linas. I am really curious about this. Could you provide us with a simpler example?

[lmockus]

Thanks for the interest.

 

I don't know how to set up a simpler example to demonstrate the behavior but I run the model with different priors (less and more tight) and stan code with jacobian still didn't finish yet while the one without jacobian is done for both priors. Will post the results shortly. Below is output for relaxed priors and stan model is attached:

-mixing for 4 chains is good

-pairs shows very good correlation for kid_vmax and kid_km in log scale

-this correlation looks bit like banana for tight priors which probably means that relaxed priors are more convincing. Relaxed priors for theta~lognorn(5,1), tight priors for theta~lognorm(0,1)

-I am quite curious why diagnose for model with jacobian gives such large errors. Any ideas?

-prediction is not perfect but makes sense (blue lines 95% cred int, red - median, black - data for 3 patients pooled together). The next step is to build population model. It is much butter for tight priors - I understand that priors still overwhelm the data. It would be good to have some scaling.

#  Log probability=-379.133
#
#  param idx           value           model     finite diff           error
#          0        0.144682         426.009         426.009   -6.39926e-008
#          1          1.2404        -88.3948        -55.9049          -32.49
#          2       -0.804449         100.322         101.914        -1.59196
#          3        0.684437        -3.35495        -3.50431        0.149352
#          4        -1.30776         215.469         215.469    1.91892e-008

 

   

For tight prior:

   

 

------ Original Message ------

From: "Jacob Leander" <jacob.e...@gmail.com>

To: "Stan users mailing list" <stan-...@googlegroups.com>

Cc: linasm...@gmail.com

Sent: 1/3/2017 9:16:24 AM

Subject: Re: Re[2]: [stan-users] stiff ODEs and population PK models

 

Lookings really interesting. I am really curious about this. Could we set up a simpler example?

[Bob Carpenter]

http://arxiv.org/abs/1510.04813

http://arxiv.org/abs/1503.08650