Dear Karline,

Thank you very much for the response and for working through a simplified version of the code. Your trick solved the problem.

Regards,
Emmanuel

From: r-sig-dynamic-models-bounces@r-project.org [mailto:r-sig-dynamic-models-bounces@r-project.org] On Behalf Of Soetaert, Karline
Sent: Tuesday, September 20, 2011 11:51 PM
To: Special Interest Group for Dynamic Simulation Models in R
Subject: Re: [R-sig-dyn-mod] ODE models

Emmanuel,

It looks as if the solver does not work for negative time steps. The easiest solution is to be pragmatic, and rather than having times run from 2009 to 2000 (with a step = -1), you run it from -2009 to -2000 (with positive step = 1).

This solves without problems. Once solved, you can then multiply the first column of the output matrix with -1 to obain positive output times.


Also, your code can be significantly simplified.

Here is my version of your small problem,

library(deSolve)
base_year   <- - 2009
dt.step     <- 1
time.window <- 9
start_year <- base_year
times <- seq(from = start_year, to = (start_year+time.window), by = dt.step)
#State variables
Y <- c(S = 94.85755, I = 2.287803, R = 2.048281)
#Loading parameters
params <- c(theta = 0.01, alpha = 0.02, beta = 0.03,
             mu = 0.001, gamma = 0.8, epsilon = 0.9)
#Model
MyModel <- function(t,Y,params){
  with(as.list(c(Y, params)), {
  dS = theta + epsilon*R + beta*I - (mu + alpha)*S
  dI = alpha*S - (mu + beta + gamma)*I
  dR = gamma*I - (mu + epsilon)*R
  return(list(c(dS = dS, dI = dI, dR = dR)))
  })
}
#Running model
Y.out <- ode(y = Y, times = times, parms = params, func = MyModel)
Y.out[,1] <- -1*Y.out[,1]
print(Y.out)
plot(Y.out)


Hope this helps,

Karline
________________________________
From: r-sig-dynamic-models-bounces@r-project.org on behalf of Drabo, Emmanuel
Sent: Wed 21/09/2011 00:54
To: 'r-sig-dynamic-models@r-project.org'
Subject: [R-sig-dyn-mod] ODE models

Dear all,

I am using the ode package to run a backward integration, but have been facing some inconsistencies in the output, which I believe to be a stability issue. I have pasted the code below with the error, as well as a different testing code for a simple SIR model that I wrote. The testing code does both the forward and backward integrations correctly, yet, with similar approach, the backward integration with the big code is failing. I would appreciate much if anyone has experienced this before and could assist with understanding what is going on and what can be done about it.

Regards,
Emmanuel


########### Warning / error message  ########
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step
     (H = step size). Solver will continue anyway.
      In above,  R1 =  0.2008048127633D+04   R2 = -0.8726792269790D-13
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step
     (H = step size). Solver will continue anyway.
      In above,  R1 =  0.2008048127633D+04   R2 = -0.8726792269790D-13
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step
     (H = step size). Solver will continue anyway.
      In above,  R1 =  0.2008048127633D+04   R2 = -0.8726792269790D-13
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step
     (H = step size). Solver will continue anyway.
      In above,  R1 =  0.2008048127633D+04   R2 = -0.8726792269790D-13
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step
     (H = step size). Solver will continue anyway.
      In above,  R1 =  0.2008048127624D+04   R2 = -0.7576576187189D-13
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step
     (H = step size). Solver will continue anyway.
      In above,  R1 =  0.2008048127624D+04   R2 = -0.7576576187189D-13
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step
     (H = step size). Solver will continue anyway.
      In above,  R1 =  0.2008048127624D+04   R2 = -0.7576576187189D-13
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step
     (H = step size). Solver will continue anyway.
      In above,  R1 =  0.2008048127624D+04   R2 = -0.9125459544063D-13
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step
     (H = step size). Solver will continue anyway.
      In above,  R1 =  0.2008048127624D+04   R2 = -0.9125459544063D-13
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step
     (H = step size). Solver will continue anyway.
      In above,  R1 =  0.2008048127623D+04   R2 = -0.1011825678787D-12
DLSODA-  Above warning has been issued I1 times.
     It will not be issued again for this problem.
      In above message,  I1 =        10
DLSODA-  At current T (=R1), MXSTEP (=I1) steps
      taken on this call before reaching TOUT
      In above message,  I1 =      5000
      In above message,  R1 =  0.2008048127618D+04
Warning in lsoda(y, times, func, parms, ...) :
  an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxsteps
Calls: as.data.frame -> ode -> lsoda -> .Call
Warning in lsoda(y, times, func, parms, ...) :
  Returning early. Results are accurate, as far as they go
Calls: as.data.frame -> ode -> lsoda -> .Call
#################################################



################# Original Code #################
#Loading parameters
#Loading libraries
require(Hmisc)
library(tgp)
library(MASS)
library(deSolve)
library(matlab)

#Switches
diff.method <-"backward"
extra.output <-TRUE

#Setting integration variables
if(toupper(diff.method)=="BACKWARD"){
  base_year <- 2009
  dt.step <- -1
  time.window <- -9
 }else{
  base_year <- 2000
  dt.step <- 1
  time.window <- 9
}
start_year <- base_year
times <- round(seq(from=start_year,to=(start_year+time.window),by=dt.step),2)

#Loading initial values
if (toupper(diff.method)=="FORWARD") {
Y <- c(1.457760e+05, 3.091247e+02, 1.063831e+04, 9.361904e+03, 1.416520e+03, 2.124780e+03, 3.953649e+03, 1.598546e+04, 9.889438e+00, 3.403382e+02, 2.995036e+02, 4.531695e+01, 6.797542e+01, 1.264841e+02, 5.114029e+02, 0, 0, 0, 0, 0, 0)
}else{
Y <- c(161466.8, 119.1497,  2709.332, 3122.580, 2642.434, 10380.388,  667.9746,  9552.27, 0.8302728,  18.46581,  39.71779,  84.39928, 1217.89197,  58.69279, 1606.9190, 15173.821, 68.36607, 19455.452, 23677.615, 32736.754, 28105.102)
}

names(Y) <- c("S","Ps","Is","Js","Es","Ts","As","TAs","Pr","Ir","Jr","Er","Tr","Ar","TAr","TotInfs","TotInfr","AIDS","Tested","Treatment","Deaths")

if (extra.output == FALSE){
Y <- Y[1:15]
}else{
Y <- Y
}

#loading parameters
params <- c(1.029718828,0.200396625,16.10975417,0.319044033,0.087139753,0.169026706,0.37111486,0.312412041,0.683488044,0.062658294,0.121984533,1.80173656,1.831660659,0.12382535,0.066790796,0.000307645,3.429975016,5.795929898,0.163094284,0.493413633,0.014542249,0.021519464,0.108808315,5.916216994,0.273186851,0.071536041,0.022121953,0.022088591,0.000628439,0.101163446,0.014902887)

extraparams <-params[25:length(params)]*params[21]
params <- c(params,extraparams)

names(params) <-c("pi","ageing","rho","omega","omega.A","chi","nu","gamma.Es","gamma.Er","gamma.Ts","gamma.Tr","sigma","sigma.A","g","g.A","mu","gamma.As","gamma.Ar","gamma.TAs","gamma.TAr","r","q","h.resistant","C.mix","beta.P","beta.I","beta.J","beta.E","beta.T","beta.A","beta.TA","beta.Pr","beta.Ir","beta.Jr","beta.Er","beta.Tr","beta.Ar","beta.TAr")

# Model function
MyModel <- function(t,Y,params,extra.output) {
  with(as.list(c(Y,params)), {

  Xs  <- Y[2:8]
  Xr <- Y[9:15]

  # Assigning state variables
  S  <- Y["S"]
  Ps  <- Y["Ps"]
  Is  <- Y["Is"]
  Js  <- Y["Js"]
  Es  <- Y["Es"]
  Ts  <- Y["Ts"]
  As  <- Y["As"]
  TAs <- Y["TAs"]
  Pr <- Y["Pr"]
  Ir <- Y["Ir"]
  Jr <- Y["Jr"]
  Er <- Y["Er"]
  Tr <- Y["Tr"]
  Ar <- Y["Ar"]
  TAr <- Y["TAr"]

  beta.Xs <- c(beta.P,beta.I,beta.J,beta.E,beta.T,beta.A,beta.TA)
  names(beta.Xs) <- X.states

  beta.Xr <- c(beta.Pr,beta.Ir,beta.Jr,beta.Er,beta.Tr,beta.Ar,beta.TAr)
  names(beta.Xr) <- Xr.states

  # Creating force of infection for each strata
  N <- S+Ps+Is+Js+Es+Ts+As+TAs+Pr+Ir+Jr+Er+Tr+Ar+TAr
  lambda <- as.numeric(sum(C.mix*beta.Xs*Xs)/N)
  lambda.r <- as.numeric(sum(C.mix*beta.Xr*Xr)/N)


  dS = pi - (mu + lambda + lambda.r)*S

  dPs = lambda*S - (mu + rho)*Ps
  dIs = rho*Ps - (mu + chi + omega + omega.A)*Is
  dJs = omega*Is - (mu + nu)*Js
  dEs = chi*Is + nu*Js + g*Ts + q*Er - (mu + gamma.Es + sigma)*Es
  dTs = sigma*Es - (mu + gamma.Ts + g + r)*Ts
  dAs = omega.A*Is + gamma.Es*Es + g.A*TAs +q*Ar -(mu + gamma.As + sigma.A)*As
  dTAs = gamma.Ts*Ts + sigma.A*As - (mu + gamma.TAs + g.A + r)*TAs

  dPr = lambda.r*S - (mu + rho)*Pr
  dIr = rho*Pr - (mu + chi + omega + omega.A)*Ir
  dJr = omega*Ir - (mu + nu)*Jr
  dEr = chi*Ir + nu*Jr + g*Tr - (mu + gamma.Er + sigma + q)*Er
  dTr = r*Ts + sigma*Er - (mu + gamma.Tr + g)*Tr
  dAr = omega.A*Ir + gamma.Er*Er + g*TAr  - (mu + gamma.Ar + sigma.A + q)*Ar
  dTAr = r*TAs + gamma.Tr*Tr + sigma.A*Ar  - (mu + gamma.TAr + g)*TAr

  RR <- as.numeric(c(dS,dPs,dIs,dJs,dEs,dTs,dAs,dTAs,dPr,dIr,dJr,dEr,dTr,dAr,dTAr))

  if(extra.output==TRUE){
    dTotInfs = lambda*S
    dTotInfr = lambda.r*S
    dAIDS = gamma.Es*Es+gamma.Er*Er+gamma.Ts*Ts+gamma.Tr*Tr+omega.A*(Is+Ir)
    dTested = (omega+omega.A+chi)*(Is+Ir)
    dTreatment = sigma*(Es+Er)+sigma.A*(As+Ar)-g*(Ts+Tr)+g.A*(TAs+TAr)
    dDeaths = gamma.As*As + gamma.TAs*TAs + gamma.Ar*Ar + gamma.TAr*TAr

    RR <- c(RR,as.numeric(c(dTotInfs,dTotInfr,dAIDS,dTested,dTreatment,dDeaths)))
  }

  names(RR) <- paste("d",names(Y),sep="")
  return(list(RR))
  })
}

#Running model
Y.out <- as.data.frame(ode(y=Y, times = times, parms = params, func = MyModel, hmax=0.025, extra.output = extra.output))
names(Y.out) <- c("times",names(Y))
print(Y.out)

#################################################



################ Test model ################
#Loading parameters
#Loading libraries
require(Hmisc)
library(tgp)
library(MASS)
library(deSolve)
library(matlab)


#Switches
diff.method <-"backward"

#Setting integration variables
if(toupper(diff.method)=="BACKWARD"){
  base_year <- 2009
  dt.step <- -1
  time.window <- -9
 }else{
  base_year <- 2000
  dt.step <- 1
  time.window <- 9
}
start_year <- base_year
times <- round(seq(from=start_year,to=(start_year+time.window),by=dt.step),2)

#State variables
Y.states <- c("S","I","R")

#Loading initial values
if (toupper(diff.method)=="FORWARD") {
Y <- c(90,10,0)
}else{
Y <- c(94.85755,2.287803,2.048281)
}

names(Y) <- c("S","I","R")
#rownames(Y) <- c(times)


#Loading parameters
data <- NULL
params <- c(0.01,0.02,0.03,0.001,0.8,0.9)
params <- cbind(data,params)
names(params) <- c("theta","alpha","beta","mu","gamma","epsilon")
params <- c(params)

#Model
MyModel <- function(t,Y,params){
with(as.list(c(Y,params)), {
# Assigning state variables
S <- Y["S"]
I <- Y["I"]
R <- Y["R"]

dS = theta + epsilon*R + beta*I - (mu + alpha)*S
dI = alpha*S - (mu + beta + gamma)*I
dR = gamma*I - (mu + epsilon)*R

RR <- as.numeric(c(dS,dI,dR))
names(RR) <- paste("d",colnames(Y),sep="")
return(list(RR))
})
}

#Running model
Y.out <- as.data.frame(ode(y=Y, times = times, parms = params, func = MyModel, hmax=0.025))
names(Y.out) <- c("times",Y.states)
print(Y.out)
################################################




Emmanuel F. Drabo
Doctoral Fellow
Pardee RAND Graduate School
RAND Corporation
1776 Main ST MSC M1N
SANTA MONICA, CA 90401-3208
Phone: 310.393.0411 x7231
Mobile: 718.419.4508
Fax: 310.260.8142



__________________________________________________________________________

This email message is for the sole use of the intended r...{{dropped:18}}

