Skip to content

Fitting compartmental model with nls and lsoda?

2 messages · Michael A. Miller, Peter Dalgaard

#
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...

Regards, Mike


##====================================================
require(odesolve)

## Simple one compartment blood flow model:
one.compartment.model <- function(t, x, parms) {
  C1 <- x[1] # compartment
  with(as.list(parms),{
    input <- approx(signal$time, signal$input, t)$y
    dC1 <- K1 * input - k2 * C1
    list(c(dC1))
  })
}

## vector of timesteps
time <- seq(0, 100)

## external signal with rectangle impulse
signal <- as.data.frame(list(time=time,
                             input=rep(0,length(time))))
signal$input[signal$time >= 10 & signal$time <=40] <- 0.2

## Parameters for steady state conditions
parms <- c(K1=0.5, k2=0.5)

## Start values for steady state
xstart <- c(C1=0)

## calculate C1 with lsoda:
C1.lsoda <- as.data.frame(lsoda(xstart, time, one.compartment.model, parms))

## Add some noise to the output curve so I can try to fit it...
C1.lsoda$noisy <- C1.lsoda$C1 + rnorm(nrow(C1.lsoda), sd=0.05*C1.lsoda$C1)

## Plot what I've got so far:
plot(input ~ time, data=signal, type='b')
points(C1.lsoda, type='b', pch=16)
points(noisy ~ time, data=C1.lsoda, col='forestgreen')

## See if I can run a fit to find the parameters that I started with...
require(nls)
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
           )
##====================================================
_                
platform i386-pc-linux-gnu
arch     i386             
os       linux-gnu        
system   i386, linux-gnu  
status                    
major    1                
minor    8.1              
year     2003             
month    11               
day      21               
language R   

Package: nls
Version: 1.8.1

Package: odesolve
Version: 0.5-8
#
mmiller3 at iupui.edu (Michael A. Miller) writes:
....
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.