[R] question about model formula

[Peter Dalgaard BSA]

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