Derivative w.r.t lag times

[charlesm93]

Summary:

A bolus dosing event causes a discontinuity in amounts w.r.t time. This leads to potential issues with NUTS because the gradient of the log posterior cannot be properly computed.

Description:

A bolus dosing event causes a discontinuity in X (the returned matrix that contains drug amounts) w.r.t time and time-dependent parameters. This is problematic when we try to estimate lag times. Consider a bolus dosing at time t = 0 with amt = 1000 and cmt = 1. The derivatives of X w.r.t t or tlag 1 at t = 0 is ill-defined because the function is discontinuous.

This issue propagates to all events and all compartments (I'm not sure why the latter happens). I have not yet tested derivatives with respect to time.

Reproducible Steps:

One can observe this issue with the unit test prim/PKModelOneCpt_test.cpp by setting discontinuous = false in util_torsten.hpp before the if (discontinuous == false) statement, line 131.

Current Output:

The unit tests do not check derivatives at discontinuous points. They do not compute derivatives w.r.t lag times in dosing compartments. We need to figure out how this impacts NUTS. We may need to revise our algorithm.

Current Version:

Stan v2.13.1
Torsten v0.82 - beta

[charlesm93]

This issue propagates to all events and all compartments (I'm not sure why the latter happens). I have not yet tested derivatives with respect to time.

What happens in other compartments is a different issue. It is specific to the case where tlag = 0 and has to do with the if (tlag != 0) conditional statement in Event.hpp. Auto-diff does not handle the change in the function (when tlag changes from 0) that occurs because the conditional statement becomes verified.

I'm discussing the problem on discourse: http://discourse.mc-stan.org/t/autodiff-w-r-t-time-lags/191/4

[charlesm93]

Regarding unit tests for torsten:

• yes, they do not test for gradients w.r.t tlags at a dosing event (evid = 1 or 4).
• yes there are issues when tlag = 0, but when this is the case, it is unlikely tlag will be a parameter.
• furthermore, the unit tests verify the gradients at all other events, and these pass. This means that even though the gradient fails to compute properly at the dosing event, autodiff still successfully computes it at all subsequent events (where there is no discontinuity in X w.r.t time).
  Edit: the unit tests skip entries when evaluating gradients w.r.t to lag times in dosing compartment. See two comments down for more details.

An implication is that if the likelihood does NOT depend X at dosing events, then the gradient of the log posterior should be properly computed. We haven't explored yet what happens when the likelihood does depend on X at a dosing event.

[charlesm93]

In any case, this matter needs to be thoroughly discussed in the user manual.

[charlesm93]

I was confronted with lag times. Playing with a one compartment model (see unit tests and util_* files), I found:

• when tlag = 0, the jacobian with respect to lag times goes to 0. This is simply because of the if (lagtime != 0) statement. Autodiff disregards lagtimes, because of this conditional.
• when tlag is nonzero, we observe disagreements between finite and auto-diff, at events which occur at the same time as the dosing, after accounting for the lag time.
• However, if the dosing does not coincide with any event (i.e time + lagtime does not coincide with the time of another event), the jacobian evaluated by autodiff agrees with that evaluated by finite diff. This is rather good news, since this last case seems the most plausible, when tlag is a parameter we are trying to evaluate.