R help - Mar 2003 - question about model formula

Dear R Gang,

I'm interested in using R and the nls package for fitting kinetic
models.  I'm having some difficulty getting a model specified for
nls though.  The math for the model that I want to fit is

   dg(t)/dt = K1 f(t) - k2 g(t)

where g(t) and f(t) are measured data at a sequence of times t.
K1 and k2 are the parameters of the model.  If I solve this, the
solution is

  g(t) = K1 \int_0^t f(t') \exp(-k2(t-t')) dt'

and I'm not sure how to write a formula in R that I can pass to
nls that will handle both the implicit loop over the data values
of t and the interpolation and numeric integration over t'.

Can anyone help me to get this properly coded for R?

Mike

-- 
Michael A. Miller                               mmiller3 at iupui.edu
  Imaging Sciences, Department of Radiology, IU School of Medicine

mmiller3 at iupui.edu (Michael A. Miller) writes:
> Dear R Gang,
> 
> I'm interested in using R and the nls package for fitting kinetic
> models.  I'm having some difficulty getting a model specified for
> nls though.  The math for the model that I want to fit is
> 
>    dg(t)/dt = K1 f(t) - k2 g(t)
> 
> where g(t) and f(t) are measured data at a sequence of times t.
> K1 and k2 are the parameters of the model.  If I solve this, the
> solution is
> 
>   g(t) = K1 \int_0^t f(t') \exp(-k2(t-t')) dt'
> 
> and I'm not sure how to write a formula in R that I can pass to
> nls that will handle both the implicit loop over the data values
> of t and the interpolation and numeric integration over t'.
> 
> Can anyone help me to get this properly coded for R?
This stuff gets tricky. The way that you *want* to solve it is to feed
the differential equation into an automatic solver like lsoda or rk4
from the odesolve package, and then the result into nls. However, even
in dead simple cases, the adaptive stepsize of the integrators makes
the likelihood discontinuous enough to confuse the optimizer in nls.
The Nelder-Mead algorithm in optim() did actually converge on our
examples, but you'll need to get that wired into nls. I have a PhD
student planning to dig into just that shortly.

One trick that is known to work is to fix the integration scheme, this
may cost you in precision, but stability is much better because you use
the same time steps and the same interpolations throughout the fit.

However, you'll probably have to decide whether you can assume f to be
known, or whether the true f(t) should be considered a parameter of a
model with a bivariate response. Then it gets *really* tricky.

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907

I did not know about stepsize being an issue.  I had thought that
problems with convergence in this case were due to bad approximations of
the finite difference gradient.  I guessed that around the optimum,
numerical errors would come to dominate the gradient calculations,
causing convergence to fail.

I've found that the issue about confusing optimizers that use gradients
can sometimes be fixed by augmenting the original system of odes with
what I believe engineers call the sensitivity equations.  If your
original equation is
dg/dt = f(t, b), where b is a parameter to be estimated,
then include

d^2g/dtdb (be sure to remember the chain rule when doing this!)

with the original equation.  With some regularity assumptions, this
integrates to dg/db, which can be used to give nls a gradient to work
with.

R. Woodrow Setzer, Jr.                        Phone: (919) 541-0128
Experimental Toxicology Division             Fax:  (919) 541-4284
Pharmacokinetics Branch
NHEERL B143-05; US EPA; RTP, NC 27711