Fitting compartmental model with nls and lsoda?

mmiller3 at iupui.edu (Michael A. Miller) writes:

A while back I started looking into using R to fit kinetic models
and I'm finally getting back to it.  It was suggested that I use
lsoda (thanks Peter Dalgaard).  So I've tried that, even though
Peter warned me that it'd be tricky.  I've put the following
example together from the lsoda help pages, but I'm not sure that
this is the correct/best way to fit a model like this.  Or maybe
this is just what Peter warned me about?  When I run the
following code, I get this error message:

  Error in qr.qty(QR, resid) : qr and y must have the same number of rows

I'm not sure where to go from there...

....

fit <- nls(noisy ~ lsoda(xstart, time, one.compartment.model, c(K1=0.5, k2=0.5)),
           data=C1.lsoda,
           start=list(K1=0.3, k2=0.7),
           trace=T
           )

Well, for a start, 

nls(noisy ~ lsoda(xstart, time, one.compartment.model, c(K1=K1, k2=k2))[,2],
etc.

works better. Notice that lsoda returns a matrix, not just the vector
of means. Also, you were supplying parameters both held constant at
0.5. If you do that, the thing at least moves:

0.3746248 :  0.3 0.7
0.007367827 :  0.3530283 0.3598175
0.007210442 :  0.3556961 0.3510272
0.002582476 :  0.4579374 0.4585350
0.002476932 :  0.4774264 0.4775360
0.002476923 :  0.4773889 0.4775243
0.002476915 :  0.4773416 0.4775221
0.002476910 :  0.4773417 0.4775218
0.00247685 :  0.4773748 0.4775564
0.002476846 :  0.4773745 0.4775566
0.002476845 :  0.4773744 0.4775566

Error in nls(noisy ~ lsoda(xstart, time, one.compartment.model, c(K1 = K1,  :
        step factor 0.000488281 reduced below `minFactor' of 0.000976562

...but now you run into a problem which is probably of the kind I
warned you against before. Be careful about those
singularities in the input function! 

Things that I know of which might help:

A. Split the region into pieces where the input is continuous to at
least 2nd order (meaning that I don't know if that is enough in
general, but first-order continuity is clearly not enough). You have 3
piecewise constant regions, so that should be easy.

B. Use a less intelligent integrator, e.g. rk4 (on a finer grid!)

C. Modify the RHS to supply a gradient with respect to the parameters
(easy, by implicit differentiation).

D. Adjust tolerances in both lsoda and nls.

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