[R-sig-dyn-mod] ODE models
Drabo, Emmanuel
edrabo at prgs.edu
Wed Sep 21 00:54:09 CEST 2011
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:6}}
More information about the R-sig-dynamic-models
mailing list