[R] Fitting compartmental model with nls and lsoda?

[Michael A. Miller]

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
           )
##====================================================


> R.version
         _                
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



-- 

Michael A. Miller                               mmiller3 at iupui.edu
  Imaging Sciences, Department of Radiology, IU School of Medicine