Hello! I'm working on COVID-19 outbreak in Colorado and am currently using the FME package to determine the best-fitting model to the data. However, FME is very slow, often taking hours to run enough iterations using the pseudo-random method to come up with good parameter estimates. If I instead use Marquardt, it takes only a minute or two, but comes up with horrible parameter estimates because it only runs a few iterations.

All the research I've done on infectious disease modeling supports use of the FME package, but I'd like to know if there are alternative packages for modeling ODE's. Does anyone have any ideas on how I can accomplish the same task using a more efficient R package? Code pasted below.

Thank you in advance!

```
## Load deSolve package
library(deSolve)
library(rootSolve)
library(coda)
library(FME)
Cp <- 5840795 ## state population
##cases <- 29 ## county cases
gam <- 1/8 ## recovery rate (1/average length of infection)
alpha <- 5.1 ## incubation period
hosp <- 0.07 ## percent of symptomatic cases hospitalized
cc <- 0.009166 ## percent of cases requiring critical care
## Define parameters that will change
parameters <- list( ## mu = 0, ## pop import/export rate - not neccessary if we assume same import/export risk
beta = 0.3696, ## transmission rate
gam = gam, ## recovery rate (1/average length of infection)
alpha = alpha, ## incubation period
ef = 0, ## effectiveness of social distancing (vary from 0.5 - 1)
pS = 0.5, ## proportion of infectious individuals symptomatic
pID = 0.3703, ## proportion of symptomatic individuals identified
siI = .399,## proportion of symptomatic individuals self isolate
lambda = 1.575, ## difference in infectiousness symptomatic/asymptomatic
hosp = hosp, cc = cc, ## percent of symptomatic cases hospitalized, ## percent of cases requiring critical care
dc = 0.5, ## proportion of ICU patients that die
mag = 0) ## percent magnitude of effectiveness of social distancing interventions
N <- Cp
trial <- function (parameters, times = seq(0, 500, 1)) {
seir1 <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
ef <- ifelse (t < 64, 0, mag) ## vary the magnitude of social distancing, from 0% before time 64 to
## defined percentages (35%, 65%, etc.) at time 64 and after
siI <- ifelse (t < 42, 0, siI) ## vary the proportion of symptomatic individuals self-isolating,
## from 0% before time 42 to defined percentages at time 42 and after
## generate differential equations for the model
dS <- - (I)*(beta*lambda*S*(1-siI)*(1-ef))/N - (beta*S*A*(1-ef))/N ## susceptible
dE <- - E/alpha + (I)*(beta*lambda*S*(1-siI)*(1-ef))/N + (beta*S*A*(1-ef))/N ## exposed (not infectious)
dI <- (E*pS)/alpha - I*(gam) ## infectious (symptomatic)
dIh <- I*hosp*gam - Ih*1/8 ## hospitalized (non-critical care)
dIc <- I*cc*gam - dc*Ic*(1/10) ## hospitalized (critical care)
dA <- (E*(1-pS))/alpha - A*gam ## infectious (asymptomatic)
dR <- I*(gam*(1-hosp+cc)) + A*gam ## recovered (implied immune)
dRh <- Ih*1/8 ## recovered from hospital
dRc <- Ic*1/10 ## recovered from critical care
dD <- dc*Ic*(1/10) ## died from COVID-19
return(list(c(dS, dE, dI, dIh, dIc, dA, dR, dRh, dRc, dD),
Id = (I+Ih+Ic)*pID,
Iht = Ih+Ic))
})
}
state <- c(S = Cp-(1), E = 0, I = 1, Ih = 0, Ic = 0, A = 0, R = 0, Rh = 0, Rc = 0, D = 0)
return(ode(y=state, times=times, func=seir1, parms=parameters))
}
out <- trial(parameters)
## change to data frame
out <- as.data.frame(out)
## read in and subset COVID-19 raw data
epi2 <- read.csv("COcases 5_04.csv", header = TRUE)
fit <- subset(epi2, select=c(time, Iht))
## determine the deviance of the model from the data
cost <- function(P){
parameters["beta"]<-P[1]
parameters["lambda"]<-P[2]
parameters["mag"]<-P[3]
parameters["siI"]<-P[4]
out <- trial(parameters)
return(modCost(out,fit))
}
## use the functions build into the R package FME
## modFit will run multiple iterations with the given lower and upper limits for these parameters
## pseudo random method will run over 1,000 iterations but can take hours
## marq method is faster, running a few iterations but generating poor parameter estimates
(Fit<-modFit(p = c(beta = 0.5, lambda = 1.6, mag = 0.8, siI = 0.4),
f = cost
,lower=c( 0.2, 1.05, 0.4, 0.1),
upper =c( 0.7, 3.00, 0.99, 0.7),
method = ("Marq")))
```