[R] estimating the starting value within a ODE using nls and lsoda

[Conti, David]

All-

I am interested in estimating a parameter that is the starting value for an ODE model. 

That is, in the typical combined fitting procedure using nls and lsoda (alternatively rk4), I first defined the ODE model:

minmod <- function(t, y, parms) {
	G <- y[1]
	X <- y[2]
  with(as.list(parms),{
    I_t <- approx(time, I.input, t)$y
    dG <- -1*(p1 + X)*G +p1*G_b
    dX <- -1*p2*X + p3*(I_t-I_b)
    list(c(dG, dX))
  })
}

Then I estimated the parameters of the model using nls:

fit.rk4 <- nls(noisy ~ rk4(p4, time, minmod, parms=c(p1,p2,p3))

However, my goal is to not only estimate p1, p2 and p3, but also to estimate the starting value p4 from the data.

I am currently using the following procedure using optim to accomplish this:

step 1 (define a DevFun to optimize):
DevFunc <- function(p4) {
	fit.rk4 <- nls(noisy ~ rk4(p4, time, ODE.model, params)
	r <- deviance(fit.rk4)
	r
}

step 2 (optimize with optim):
fit <- optim(inits, DevFunc, method="L-BFGS-B", hessian=T, lower=lower, upper=upper)

I have explored this via simulation and with real data and while the procedure works for most cases, there are a few cases that cause difficulty.  I think it has to do with estimating p1, p2 and p3 conditional on p4 as opposed to estimating all parameters jointly.

Does anyone know of a way in [R] to estimate the starting value jointly with the parameters within a ODE model?  I have included more detailed code below to "spin a perfect cycle" to give a better idea of what I am trying to accomplish.

Thanks in advance,
Dave

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

### ODE model
minmod <- function(t, y, parms) {
	G <- y[1]
	X <- y[2]
  with(as.list(parms),{
    I_t <- approx(time, I.input, t)$y
    dG <- -1*(p1 + X)*G +p1*G_b
    dX <- -1*p2*X + p3*(I_t-I_b)
    list(c(dG, dX))
  })
}

##### simulate data
I.input <- c(5.8,7.8,11.5,10.1,9.9,9.8,12.1,10.7,11.1,11.5,11.2,10.3,171.0,179.0,145.6,147.0,105.0,20.2,15.3,14.7,12.3,10.4,7.7,7.0,7.5,4.9,5.6,5.8)
time <- c(0,1,3,4,6,7,8,10,12,14,16,19,22,23,24,25,27,40,50,60,70,80,90,100,120,140,160,180)
G_b <- 100
I_b <- 5.8
parms <- c(p1=0.012584, p2=0.037708, p3=0.000012524)
xstart <- c(G=273.92,X=0)
sim.data <- as.data.frame(rk4(y=xstart, time=time, func=minmod, parms=parms))  # note: i use rk4 because the standard software for estimating this model uses a Runge-Kutta approach
sim.data$noisy <- sim.data$G + rnorm(nrow(sim.data), mean=0, sd=0.01*sim.data$G)

###### estimation
Weight <- c(0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
xstart.est <- c(G=mean(sim.data$noisy[1:4]), X=0)

# inital fit of nls
fit.rk4 <- nls(noisy ~ rk4(xstart.est, time, minmod, c(p1=p1, p2=p2, p3=p3))[,2],
               data=sim.data,
               start=list(p1=0.025, p2=0.02, p3=0.00004),
               trace=F,
               control=list(tol=1e-4, minFactor=1/1024/1024, warnOnly=T),
               weight=(1/(G*0.015))*Weight
               )

# define deviance function to optimize               
DevFunc <- function(p4) {
	fit.rk4 <- nls(noisy ~ rk4(c(G=p4, X=0), time, minmod, c(p1=p1, p2=p2, p3=p3))[,2],
               data=sim.data,
               start=list(p1=summary(fit.rk4)$coef[1,1], p2=summary(fit.rk4)$coef[2,1], p3=summary(fit.rk4)$coef[3,1]),
               trace=F,
               control=list(tol=1e-4, minFactor=1/1024/1024, warnOnly=T),
               weight=(1/(G*0.015))*Weight
               )     
  		r <- deviance(fit.rk4)
		r
}

##### estimation via optim
inits=c(mean(sim.data$noisy[1:5]))
lower=G_b
upper=500
fit <- optim(inits, DevFunc, method="L-BFGS-B", hessian=T, lower=lower, upper=upper)

#extract p1, p2, and p3 conditional on estimated p4:
fit.rk4 <- nls(noisy ~ rk4(c(G=fit$par, X=0), time, minmod, c(p1=p1, p2=p2, p3=p3))[,2],
               data=sim.data,
               start=list(p1=summary(fit.rk4)$coef[1,1], p2=summary(fit.rk4)$coef[2,1], p3=summary(fit.rk4)$coef[3,1]),
               trace=F,
               control=list(tol=1e-4, minFactor=1/1024/1024, warnOnly=T),
               weight=(1/(G*0.015))*Weight
               )

# get all four parameters:
p <- list(p1=summary(fit.rk4)$coef[1,1], p2=summary(fit.rk4)$coef[2,1], p3=summary(fit.rk4)$coef[3,1], p4=fit$par)
p






David V. Conti, Ph.D.
Associate Professor

University of Southern California
Zilkha Neurogenetics Institute
Keck School of Medicine
Department of Preventive Medicine
Division of Biostatistics

Mailing Address:
University of Southern California
1501 San Pablo Street
ZNI 101, MC2821
Los Angeles, CA 90089-2821

Physical Address:
USC/Harlyne Norris Research Tower
1450 Biggy Street
NRT 2509L
Los Angeles, CA 90033

Tel: 323-442-3140
Fax: 323-442-2448
Email: dconti at usc.edu