R help - Aug 2006 - Fitting data with optim or nls--different time scales

Hi,

I have a system of ODE's I can solve with lsoda.

Model=function(t,x,parms) 
{
    #parameter definitions
    lambda=parms[1]; beta=parms[2]; 
    d = parms[3]; delta = parms[4]; 
     p=parms[5];    c=parms[6]
     
      xdot[1] = lambda - (d*x[1])- (beta*x[3]*x[1])
      xdot[2] = (beta*x[3]*x[1]) - (delta*x[2])
      xdot[3] = (p*x[2]) - (c*x[3])
     
    return(list(xdot))
}

I want to fit the output out[,4] to experimental data that is only 
available on days 0, 7, 12, 14, 17, and 20. I don't know how to set up 
optim or nls so that it takes out[,4] on the appropriate day, but still 
runs lsoda on a time scale of 0.01 day.

Below is the function I've been using to run 'optim', at the 
course-grained time scale:

Modelfit=function(s) {
	parms[1:4]=s[1:4]; times=c(0,7,12,14,17,20,25)
	out=lsoda(x0,times,Model,parms)
	mse=mean((log10(out[,4])-log10(i(times)))^2)
#	cat(times)
	return(mse)
}
#x0=c(T0,I0,V0)
x0=c(2249,0,1)
#parms(lambda, beta, d, delta, p, c)
parms[5:6]=c(1.0,23)

s0=c(49994,8456,6.16E-8,0.012) #initial values

fit=optim(s0,Modelfit)

Right now, lsoda is being run on too course-grained a time scale in the 
function Modelfit. Most examples of optim and nls I have found compare 
two data sets at the same times, and run lsoda on the time scale the 
data is available at, but I would like to run lsoda at a finer scale, and 
only compare the appropriate time points with the experiment.  I have also 
tried using nls, but I have the same problem. Does anyone have 
suggestions?

Thank you very much,

Leslie

<see in line>

Leslie Chavez wrote:
> Hi,
> 
> I have a system of ODE's I can solve with lsoda.
> 
> Model=function(t,x,parms) 
> {
>     #parameter definitions
>     lambda=parms[1]; beta=parms[2]; 
>     d = parms[3]; delta = parms[4]; 
>      p=parms[5];    c=parms[6]
>      
>       xdot[1] = lambda - (d*x[1])- (beta*x[3]*x[1])
>       xdot[2] = (beta*x[3]*x[1]) - (delta*x[2])
>       xdot[3] = (p*x[2]) - (c*x[3])
>      
>     return(list(xdot))
> }
> 
> I want to fit the output out[,4] to experimental data that is only 
> available on days 0, 7, 12, 14, 17, and 20. I don't know how to set up 
> optim or nls so that it takes out[,4] on the appropriate day, but still 
> runs lsoda on a time scale of 0.01 day.
> 
> Below is the function I've been using to run 'optim', at the 
> course-grained time scale:
> 
SG:  What about the following:

  Modelfit=function(s) {
  	parms[1:4]=s[1:4];
	times=c(0,7,12,14,17,20,25)
	lsodaTimes <- seq(min(times),max(times), by=0.01)
  	out=lsoda(x0,lsodaTimes,Model,parms)
	obsTimes <- (100*times-1)
  	mse=mean((log10(out[obsTimes,4])-log10(i(times)))^2)
  #	cat(times)
  	return(mse)
  }

	  Your example is not self contained, so obviously I haven't tried this 
with it.  However, something of this nature should work fine, I believe. 
  Something similar but different should also work, I believe, with 
'nls';  this would give you access to many helper functions (see 
"methods(class='nls')").  If 'nls' bombed on me,
I'd then try 'optim' as
it is less brittle.  Then I might use the output of 'optim' as initial 
values for 'nls' to get confidence intervals etc.

	  hope this helps.
	  Spencer Graves
> #x0=c(T0,I0,V0)
> x0=c(2249,0,1)
> #parms(lambda, beta, d, delta, p, c)
> parms[5:6]=c(1.0,23)
> 
> s0=c(49994,8456,6.16E-8,0.012) #initial values
> 
> fit=optim(s0,Modelfit)
> 
> Right now, lsoda is being run on too course-grained a time scale in the 
> function Modelfit. Most examples of optim and nls I have found compare 
> two data sets at the same times, and run lsoda on the time scale the 
> data is available at, but I would like to run lsoda at a finer scale, and 
> only compare the appropriate time points with the experiment.  I have also 
> tried using nls, but I have the same problem. Does anyone have 
> suggestions?
> 
> Thank you very much,
> 
> Leslie
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.