charlesm93 / example-models

Example models for Stan
http://mc-stan.org/
5 stars 1 forks source link

Feature Request: Initial Conditions for ODE Solvers #8

Open billdenney opened 7 years ago

billdenney commented 7 years ago

For generalOdeModel, PKModelCpt, and linOdeModel, it would be nice to be able to set the initial conditions. That will enable simpler handling of endogenous compounds when at time 0, the amount is nonzero.

As a thought, it could be developed within the current function definitions so that if nCmt is an integer it is replicated with zeros; if it is a real vector, then the initial conditions are nCmt and the number of compartments is the size of nCmt.

charlesm93 commented 7 years ago

Great idea. The solution you propose works well for generalOdeModel and mixOdeModel but the other functions will need an extra argument. The signature could look like:

PKModelOneCpt(time, amt, rate, ii, evid, cmt, addl, ss, theta, F, tlag, init)

with init at the end. In future releases, the last arguments (F, tlag, and init) will be optional.

Implementing initial conditions should be straightforward. I think you can currently hack it by using multiple dosing events at time t = 0, which is of course not ideal. Maybe @billgillespie knows other ways to implement initial conditions.

billdenney commented 7 years ago

Great!

The work-around of a dose at time 0 generally works, but if the first measurement time starts before 0, that has to be detected, set, and then the system can start. (It's error prone because time = 0 may not be the first measurement.)

As an interface thought: It would be nice if the calls had identical signatures or as near identical as possible. That way, changes between methods would be simpler. You've had that until now, and I don't want this to introduce an incompatibility or a signature mismatch (mainly with argument order changes). I'd like to be sensitive between adding new features quickly and long term usability (what works in version 0.83 will hopefully work in version 1.0).

My concern with the interface is that as you add many optional arguments keeping track of which optional arguments are in effect can become hard for the user (aka me) and I may want to have different arguments in effect at different times such as F and init but not tlag or init and tlag but not F.

billgillespie commented 7 years ago

Agree that I would like the signatures as similar as possible across the Torsten functions. If we introduce initial values as an argument I am inclined to make it a required argument. In fact I am leaning toward keeping F and tlag as required arguments too.

I am leaning toward adding initial values as a separate argument for the following reasons:

The only negative that immediately comes to mind is increasing th number of arguments...

billdenney commented 7 years ago

With @billgillespie's comments, I'll weigh in that I think making all arguments mandatory is probably the best path. Having a mix of mandatory and optional arguments will make life more complex (and to harp on it, error prone).

Perhaps the most useful for users (aka my way of working) would be to have two signatures: standard minimal case (similar to what is in the system now) and full complexity case (initial values, tlag, F, ...).

While thinking of this interface, I can imagine a few different signatures within each of those:

Scalar vs. Vector Scalar values get applied across all measurements. Vector values get applied specifically to the row where it is observed. The use case here is modeling multiple drugs simultaneously. When modeling more than one drug, I will want to have flexibility to set parameters (e.g. tlag, F) separately for each drug.

Standard Values Most of the time, I will not be fitting some of the optional values in the full complexity case. The function signatures should be forgiving when I'm not using the full complexity. As an example, F is a real (not integer) value. But, if I accidentally give F=1 as an integer (the default), the code should not complain. Values other than the default (e.g. F=0 as an integer), should require the correct type.

Consistency Checks (I'm a Torsten novice at this point, so this may already be in the code.) When two values are inconsistent, as informative as possible an error should be given. As an example, if the differential equations are for a 3 compartment model and the initial conditions are for a 2 compartment model, please tell me specifically that rather than something like "Differential equation result size incorrect." (Again, I've not checked-- you may already be good there.)

charlesm93 commented 7 years ago

Perhaps the most useful for users (aka my way of working) would be to have two signatures: standard minimal case (similar to what is in the system now) and full complexity case (initial values, tlag, F, ...).

Yes, this is exactly what I had in mind: a minimal and a full complexity case. Just realized I didn't make this clear at all.

Scalars vs Vectors [...]

We currently allow users to pass theta (the ODE parameters), F, and tlag as either arrays or 2D arrays. In the first case, there is one element per compartment, which applies to all events. In the second case, each array element applies to one event.

I'm not sure how easy this makes it to test multiple drugs. You can have some events correspond to one drug. That said, I'm guessing different patients receive different drugs? If that's true, you need to pass the relevant theta, F, and tlag for each patient. Recall users call a torsten function for each patient, see section 3.1. of the manual.

Standard Values [...]

This would require modifying Stan's parser. The cost/benefit ratio isn't that good. If the user specified F as an integer instead of a real, the code will (unnecessarily) complain but it's a straightforward fix.

Consistency Check [...]

Yes, we indeed have some consistency checks. For instance, time and amt are expected to be arrays of the same length. The example you give is implemented (through Stan's built-in ODE integrator). I'll make a list of the consistency checks, to make sure I did not miss important cases. In addition, I need to write unit tests for the error messages the code produces.