I have been using method deSolve::ode45 which has been working until I made a few necessary changes to my equations. Does anyone know why the ODE solver is not working? I have tried running with ode45 as well as the default ode method and neither work. Please let me know if any further explanation would be helpful.
I have checked over the differential equations and I am confident they are correct.
The equations used are as follows:
CCHFModel = function(t,x,params)
{
# get SIR values
SH <- x[1]
EH <- x[2]
IA <- x[3]
IS <- x[4]
RH <- x[5]
ST <- x[6]
IT <- x[7]
SC <- x[9]
IC <- x[10]
RC <- x[11]
# Load values ----
# Beta values
betaHHA = params["betaHHA"]
betaHHS = params["betaHHS"]
betaTH = params["betaTH"]
betaCH = params["betaCH"]
betaTC = params["betaTC"]
betaCT = params["betaCT"]
betaTT = params["betaTT"]
# Gamma value
gamma = params["gamma"]
# death rates
muH = params["muH"]
muT = params["muT"]
muC = params["muC"]
# birth rates
piH = params["piH"]
piT = params["piT"]
piC = params["piC"]
# incubation
deltaHS = params["deltaHS"]
deltaHA = params["deltaHA"]
# recovery rate
alphaA = params["alphaA"]
alphaS = params["alphaS"]
alphaC = params["alphaC"]
# total population
NH = (SH + IA + IS + EH + RH) + (piH * SH) - (muH * SH)
NT = (ST + IT) + (piT * ST) - (muT * ST)
NC = (SC + IC + RC) + (piC * SC) - (muH * SC)
# tick carrying Capacity
# KT = NC * 130 # 130 ticks per carrier max
#computations ----
dSHdt <- (piH * NH) - (betaHHA * IA + betaHHS * IS + betaCH * IC + betaTH * IT)*(SH/NH) - (muH * SH)
dEHdt <- (betaHHA * IA + betaHHS * IS + betaCH * IC + betaTH * IT)*(SH/NH) - ((deltaHA + muH)*EH)
dIAdt <- (deltaHA * EH) - ((alphaA + muH + deltaHS) * IA)
dISdt <- (deltaHS * IA) - ((alphaS + muH + gamma) * IS)
dRHdt <- alphaA * IA + alphaS * IS - muH*RH
dSTdt <- (piT * NT) - (betaTT * IT + betaCT * IC)*(ST/NT) - (muT * ST)
dITdt <- (betaTT * IT + betaCT * IC)*(ST/NT) - (muT * IT)
dSCdt <- (piC * NC) - (betaTC * IT)*(SC/NC) - (muC * SC)
dICdt <- (betaTC * IT)*(SC/NC) - ((alphaC +muC) * IC)
dRCdt <- (alphaC * IC) - (muC * RC)
# return results
list(c(dSHdt, dEHdt, dIAdt, dISdt, dRHdt, dSTdt, dITdt, dSCdt, dICdt, dRCdt))
}
I run the ODE solver using:
defaultParms = c(betaHHA = .0413,
betaHHS = .0413,
betaTH = .2891,
betaCH = .0826,
betaTC = (1/365),
betaCT = 59/365,
betaTT = ((1/(365 * 2)) * .04) * 280,
gamma = 1/10,
muH = (1/(365 * 73)),
muT = (1/(365 * 2)),
muC = (1/(11 * 365)),
piH = 1.25/(73 * 365),
piT = 4.5/730,
piC = 1/(11 * 365),
deltaHS = 1/3,
deltaHA = 1/2,
alphaA = 1/17,
alphaS = 1/17,
alphaC = 1/7)
# time to start solution
t = seq(from = 0, to = 365, by = 0.1)
#initialize initial conditions
initialConditions = c(SH = 10000, EH = 5, IA = 5, IS = 10, RH = 2, ST = 80000, IT = 50, SC = 30000, IC = 5, RC = 1)
dataSet = ode(y = initialConditions, times = t, func = CCHFModel, parms = defaultParms)%>%
as.data.frame()
After running this all the output following the initial conditions is NA.
This is due to a typo - you misnumbered the translation of input values in the first section of your code (i.e., you skipped x[8]. I will go through two (hopefully) useful exercises, first explaining how I debugged this and then showing how to rewrite your function to make it less error-prone ...
t=0, x=<initial conditions>:CCHFModel(0,initialConditions, defaultParms)
## piH betaHHA deltaHA deltaHS alphaA piT
## -15.02882327 12.62349834 0.53902803 0.07805607 0.88227788 385.31052332
## betaTT piC betaTC alphaC
## 0.85526763 NA NA NA
Hmm, we already see we have a problem. Why are the last three elements of the computed gradients NA?
add browser() near the end of the function (before the dsCdt <- ... line) so we can take a closer look. Redefine the function, and try computing the gradient again.
When we get there and print out some of the quantities involved in the computation we see that both NC and RC are NA ... we can also see that an NA value of RC will cause NC to be NA, so let's check the definition of RC ...
aha! RC is defined as x[11], but length(initialConditions) is only 10 ... and a closer look shows that we missed x[8]. Redefining properly gives non-NA values throughout (I don't know if they're correct, but at least they're not NA).
Although using [] or [[]] to extract elements of a vector usually give equivalent answers, you should always use [[]] when you want to extract a single element (scalar) from a vector. Here's why:
initialConditions[11] ## NA
initialConditions[[11]] ## Error in x[[11]] : subscript out of bounds
If you use [], the NA propagates through your code and you have to hunt down the original source. If you use [[]], R fails right away and tells you where the problem is. An additional benefit is that [] propagates the names of the vector elements in a way that doesn't usually make sense (take a look at the names of the output in "debugging/1" above ...)
You can avoid all of the tedious and error-prone unpacking of the parameter and state vectors by replacing the unpacking code (everything before the computation of total populations) with
comb <- c(as.list(x), as.list(params))
attach(comb)
on.exit(detach(comb))
Provided that your parameter and state vectors are properly named (and there are no names that overlap between them), this will create a named list and allow looking up of the elements by name within your function; on.exit(detach(comb)) makes sure that everything gets cleaned up properly at the end. (You will see recommendations to use with() to do this; I prefer the strategy here because it makes debugging within the function [if necessary] easier. But as @tpetzoldt notes in comments, you should always pair attach(...) with on.exit(detach(...)); otherwise things get very confusing and messy ...)
At the end of the function I would use
g <- c(dSHdt, dEHdt, dIAdt, dISdt, dRHdt, dSTdt, dITdt, dSCdt, dICdt, dRCdt)
names(g) <- names(x)
list(g)
to make sure the gradient vector is properly labeled, which makes troubleshooting easier.