COVID-19 Model Fitting Alternatives to FME Package

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")))