[R] Fitting compartmental model with nls and lsoda?
Michael A. Miller
mmiller3 at iupui.edu
Thu Jan 22 22:01:35 CET 2004
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
More information about the R-help
mailing list