[R] Fitting compartmental model with nls and lsoda?
Peter Dalgaard
p.dalgaard at biostat.ku.dk
Thu Jan 22 22:59:28 CET 2004
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
More information about the R-help
mailing list