From: Ravi Varadhan <rvaradhan_at_jhmi.edu>

Date: Mon, 11 Jun 2007 15:55:38 -0400

p.s. Your example seems not to be self contained. If I could have easily copied it from your email and run it myself, I might have been able to offer more useful suggestions.

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.

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. Received on Mon 11 Jun 2007 - 20:38:18 GMT

Date: Mon, 11 Jun 2007 15:55:38 -0400

Jeremy,

You should examine the steady-state solution to your system of equations, by setting the time-derivatives to zero and then solving/analyzing the resulting algebraic equations. This should give you some insights.

Let us say you have 3 groups, A,B, and C, with initial conditions: N_A(t=0) = N_{A0}, N_B(t=0) = N_{B0}, and N_C(t=0) = N_{C0}, and that people transition in and out of these 3 states (one of the states could even be absorbing, e.g. death), but it is true for any time t that N_A(t) + N_B(t) + N_C(t) = N_{A0} + N_{B0} + N_{C0}. Furthermore, you have 3 diff-equations that describe the rate of change of N_A, N_B, and N_C, for t > 0. If it happens that one of the N's, say N_A, becomes negative, you could set it equal to zero. But then you have to figure out how to re-adjust N_B and N_C so that they add up to the initial total count. After re-adjustment, you also have to think about whether the resulting system of equations are valid, when there are no A people.

Ravi.

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan_at_jhmi.edu

Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

>From your response, I am not sure if I asked my question clearly.

I am modeling a group of people (in a variety of health states) moving through time (and getting infected with an infectious disease). This means that the count of the number of people in each state should be positive at all times.

What appears to happen is that lsoda asks for a derivative at a given point in time t and then adjusts the state of the population. However, perhaps due to numerical instability, it occasionally lower the population count below 0 for one of the health states (perhaps because it's step size is too big or something).

I have tried both the logarithm trick and also changing the relative and absolute tolerance inputs but I still get the problem for certain combinations of parameters and initial conditions.

It occurs both under MS Windows XP Service Pack 2 and on a Linux cluster so I am pretty sure it is not platform specific.

My real question to the group is if there is not a work around in lsoda are there other ode solvers in R that will allow the constraint of solutions to the ODEs remain non-negative?

Best regards,

Jeremy

>>> Spencer Graves <spencer.graves@pdf.com> 6/8/2007 9:51 AM >>>

On the 'lsoda' help page, I did not see any option to force some
or all parameters to be nonnegative.

Have you considered replacing the parameters that must be nonnegative with their logarithms? This effective moves the 0 lower limit to (-Inf) and seems to have worked well for me in the past. Often, it can even make the log likelihood or sum of squares surface more elliptical, which means that the standard normal approximation for the sampling distribution of parameter estimates will likely be more accurate.

Hope this helps. Spencer Graves

p.s. Your example seems not to be self contained. If I could have easily copied it from your email and run it myself, I might have been able to offer more useful suggestions.

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
**>
**> 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))
**>
**> }
**>
**> ______________________________________________
**> 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.

*>
*

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.

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. Received on Mon 11 Jun 2007 - 20:38:18 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Mon 11 Jun 2007 - 21:32:20 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*