R help - Jun 2007 - Using odesolve to produce non-negative solutions

Hello,

I am using odesolve to simulate a group of people moving through time and
transmitting infections to one another.

In Matlab, there is a NonNegative option which tells the Matlab solver to keep
the vector elements of the ODE solution non-negative at all times. What is the
right way to do this in R?

Thanks,
Jeremy

P.S., Below is a simplified version of the code I use to try to do this, but I
am not sure that it is theoretically right

dynmodel <- function(t,y,p) 
{ 
## Initialize parameter values

	birth <- p$mybirth(t)
	death <- p$mydeath(t)
	recover <- p$myrecover
	beta <- p$mybeta
	vaxeff <- p$myvaxeff
	vaccinated <- p$myvax(t)

	vax <- vaxeff*vaccinated/100

## If the state currently has negative quantities (shouldn't have), then
reset to reasonable values for computing meaningful derivatives

	for (i in 1:length(y)) {
		if (y[i]<0) {
			y[i] <- 0
		}
	}

	S <- y[1]
	I <- y[2]
	R <- y[3]
	N <- y[4]

	shat <- (birth*(1-vax)) - (death*S) - (beta*S*I/N)
	ihat <- (beta*S*I/N) - (death*I) - (recover*I)
	rhat <- (birth*(vax)) + (recover*I) - (death*R)

## Do we overshoot into negative space, if so shrink derivative to bring state
to 0
## then rescale the components that take the derivative negative

	if (shat+S<0) {
		shat_old <- shat
		shat <- -1*S
		scaled_transmission <- (shat/shat_old)*(beta*S*I/N)
		ihat <- scaled_transmission - (death*I) - (recover*I)
		
	}	
	if (ihat+I<0) {
		ihat_old <- ihat
		ihat <- -1*I
		scaled_recovery <- (ihat/ihat_old)*(recover*I)
		rhat <- scaled_recovery +(birth*(vax)) - (death*R)
	
	}	
	if (rhat+R<0) {
		rhat <- -1*R
	}	

	nhat <- shat + ihat + rhat

	if (nhat+N<0) {
		nhat <- -1*N	
	}	

## return derivatives

	list(c(shat,ihat,rhat,nhat),c(0))

}

Hello,

I am using odesolve to simulate a group of people moving through time and
transmitting infections to one another.

In Matlab, there is a NonNegative option which tells the Matlab solver to keep
the vector elements of the ODE solution non-negative at all times. What is the
right way to do this in R?

Thanks,
Jeremy

P.S., Below is a simplified version of the code I use to try to do this, but I
am not sure that it is theoretically right

dynmodel <- function(t,y,p) 
{ 
## Initialize parameter values

	birth <- p$mybirth(t)
	death <- p$mydeath(t)
	recover <- p$myrecover
	beta <- p$mybeta
	vaxeff <- p$myvaxeff
	vaccinated <- p$myvax(t)

	vax <- vaxeff*vaccinated/100

## If the state currently has negative quantities (shouldn't have), then
reset to reasonable values for computing meaningful derivatives

	for (i in 1:length(y)) {
		if (y[i]<0) {
			y[i] <- 0
		}
	}

	S <- y[1]
	I <- y[2]
	R <- y[3]
	N <- y[4]

	shat <- (birth*(1-vax)) - (death*S) - (beta*S*I/N)
	ihat <- (beta*S*I/N) - (death*I) - (recover*I)
	rhat <- (birth*(vax)) + (recover*I) - (death*R)

## Do we overshoot into negative space, if so shrink derivative to bring state
to 0
## then rescale the components that take the derivative negative

	if (shat+S<0) {
		shat_old <- shat
		shat <- -1*S
		scaled_transmission <- (shat/shat_old)*(beta*S*I/N)
		ihat <- scaled_transmission - (death*I) - (recover*I)
		
	}	
	if (ihat+I<0) {
		ihat_old <- ihat
		ihat <- -1*I
		scaled_recovery <- (ihat/ihat_old)*(recover*I)
		rhat <- scaled_recovery +(birth*(vax)) - (death*R)
	
	}	
	if (rhat+R<0) {
		rhat <- -1*R
	}	

	nhat <- shat + ihat + rhat

	if (nhat+N<0) {
		nhat <- -1*N	
	}	

## return derivatives

	list(c(shat,ihat,rhat,nhat),c(0))

}

Dear Jeremy,

a few notes about your model: The equations of your derivatives are 
designed in a way that can lead to negative state variables with certain 
parameter combinations. In order to avoid this, you are using "if 
constructions" which are intended to correct this. This method is 
however (as far as I have learned from theory and own experience ;-) a 
bad idea and should be strictly avoided.

This trick may violate assumptions of the ODE solvers and in many cases 
also mass-balances (or similar).

Instead of this, processes should be modeled in a way that avoids 
"crossing zero", e.g. exponential decay (x, k > 0):

dx = - k * x  (1)

and not a linear decay like

dx = -k       (2)

which by its nature can lead to negative state values.

Case (1) can be managed almost perfectly by lsoda with his automatic 
internal time step algorithm. hmax is intended for non-autonomous models 
to ensure that external signals are not skipped by the automatism, which 
may be appropriate in your case because p seems to contain time 
dependent functions.

As far as I can see, your equations belong to type (1) and should not 
need any of the if and for constructs, as long as your parameters (which 
are not given in your post) do have the correct sign and range (for 
example, vax <= 1, death >= 0  etc.).

If you perform optimization, you must ensure that parameters stay in the 
plausible range is met using transformations of the parameters or 
constraints in the optimization procedure.

Thomas

PS: another question, what is the purpose of the state variable N?
I guess it can be derived from the other states.


Jeremy Goldhaber-Fiebert wrote:
> Hello,
> 
> I am using odesolve to simulate a group of people moving through time
> and transmitting infections to one another.
> 
> In Matlab, there is a NonNegative option which tells the Matlab
> solver to keep the vector elements of the ODE solution non-negative
> at all times. What is the right way to do this in R?
> 
> Thanks, Jeremy
[... code deleted, see original post ...]

If you didn't get this solved.

I have done parameter estimation with models
defined by ODE's where negative solutions are a problem
and one can only avoid them with great difficulty if the
standard explicit methods for solving the ODE are used. I found that
using implicit methods could be a great help.

For example in the exponential case

     dN/dt = -k*N

the simple finite difference approximation is

     N_{t+1}-N+t
     --------------  = -k*N_t ,  k>=0
        h

   or
     N_{t+1} = N_t -k*h*N_t

   and if  k*h gets too large N_{t+1} goes negative and you are in trouble.

   Consider instead the implicit formulation where the second
   N_t on the RHS is replaced by N_{t+1}  and one gets

     N_{t+1} = N_t/(1+k*h)

    which is correct  for k*h=0 and as k*h--> infinity

   For a more complicated example see

    http://otter-rsch.com/admodel/cc4.html

   for something I called "semi-implicit".
   I hope these ideas will be useful for your problem.


         Cheers,

          Dave





-- 
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com