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