[R] Fitting compartmental model with nls and lsoda?

[Michael A. Miller]

>>>>> "DivineSAAM" == DivineSAAM  <DivineSAAM at aol.com> writes:

    > Dear Colleagues,
    > Our group is also working on implementing the use of R for
    > pharmacokinetic compartmental analysis. Perhaps I have
    > missed something, but 

    > > 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
    > +            )
    > Error in eval(as.name(varName), data) : Object "C1.lsoda" not found

    > What part of the e-mail did I miss? I would like to get
    > this problem up an running. 

Oscar,

There were several problems with the code in my previous posting.
I'll append an example that does work.  While this example often
works, there are cases when nls fails by reducing the step factor
below minFactor.  It is even less stable for fitting less trivial
examples with real data sets.  I'll be very interested to keep in
touch as we make progress on this problem.

Regards, Mike


Here's my (more) correct example:

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, 1)

## 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
C1.lsoda$noisy <- C1.lsoda$C1 + rnorm(nrow(C1.lsoda), sd=0.15*C1.lsoda$C1)

## 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=K1, k2=k2))[,2],
           data=C1.lsoda,
           start=list(K1=0.3, k2=0.7),
           trace=T,
           control=list(tol=1e-2,
             minFactor=1/1024/1024)
           )

fit.rk4 <- nls(noisy ~ rk4(xstart, time, one.compartment.model, c(K1=K1, k2=k2))[,2],
               data=C1.lsoda,
               start=list(K1=0.3, k2=0.7),
               trace=T,
               control=list(tol=1e-2,
                 minFactor=1/1024/1024)
               )

## Plot what I've got so far:
par(mfrow=c(2,2))
plot(noisy ~ time, data=C1.lsoda, main='Input, C1, C1+noise', col='forestgreen')
points(input ~ time, data=signal, type='b')
points(C1 ~ time, data=C1.lsoda, type='b', pch=16)

t <- seq(0,100,0.1)
plot(noisy ~ time, data=C1.lsoda, main='Input, C1+noise, lsoda',
     col='forestgreen')
lines(t, predict(fit, list(time=t)), col='red',type='l')

plot(noisy ~ time, data=C1.lsoda, main='Input, C1+noise, rk4',
     col='forestgreen')
lines(t, predict(fit.rk4, list(time=t)), col='blue', type='l')

plot(t, predict(fit.rk4, list(time=t))-predict(fit, list(time=t)),
     type='l')
abline(h=0)


print('Coefficients for lsoda solution:')
print(coef(fit))
print(vcov(fit))

print('Coefficients for rk4 solution:')
print(coef(fit.rk4))
print(vcov(fit.rk4))

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