[R] Fitting compartmental model with nls and lsoda?

[Peter Dalgaard]

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