CODEWITHSUNDEEP
r
moistbananacake
moistbananacake

Reputation: 25

Problem in finding parameters for SIR model in R

I am trying to model the coronavirus infectious rate using SIR model in R. However, my beta (controls the transition between Susceptibles and Infected) and gamma (controls the transition between Infected and Recovered) values are not correct as you can see. 

      beta     gamma 
 1.0000000 0.8407238

Here is my code:

library(deSolve)


Infected <- c(994, 307, 329, 553, 587, 843, 983, 1750, 2950, 4569, 5632, 4848, 9400, 10311, 11166, 13451, 17388, 18743, 19452, 20065, 20732,24914)
Day <- 1:(length(Infected))
N <- 331002651 # population of the us

This is the SIR function

SIR <- function(time, state, parameters) {
  par <- as.list(c(state, parameters))
  with(par, {
    dS <- -beta/N * I * S
    dI <- beta/N * I * S - gamma * I
    dR <- gamma * I
    list(c(dS, dI, dR))
  })
}

Then I find the RSS

init <- c(S = N-Infected[1], I = Infected[1], R = 0)
RSS <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- ode(y = init, times = Day, func = SIR, parms = parameters)
  fit <- out[ , 3]
  sum((Infected - fit)^2)
}

Opt <- optim(c(0.5, 0.5), RSS, method = "L-BFGS-B", lower = c(0, 0), upper = c(1, 1)) # optimize with some sensible conditions
Opt$message

Then I find beta and gamma

Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par
##      beta     gamma 
## 1.0000000 0.8407238

What went wrong with my code and how can I fix it? Thank you in advance!

Upvotes: 2

Views: 422

Answers (1)

tpetzoldt
tpetzoldt

Reputation: 5813

Optimization is often tricky, there are several things that may be improved:

  • the model is very simple, you may consinder to use a SEIR model instead (being still to simple),
  • extend the box constraints, this helped in your case,
  • use another optimizer, e.g. from a meta-package like FME, that has also additional diagnostics,
  • re-scale variables to a more convenient range, e.g. by dividing the population by 1e3,
  • OR, as an alternative, adapt tolerances (atol) appropriate to the original scales.
library(deSolve)

Infected <- c(994, 307, 329, 553, 587, 843, 983, 1750, 2950, 4569, 5632, 4848, 
              9400, 10311, 11166, 13451, 17388, 18743, 19452, 20065, 20732,24914)
Day <- 1:(length(Infected))
N <- 331002651 # population of the us

SIR <- function(time, state, parameters) {
  par <- as.list(c(state, parameters))
  with(par, {
    dS <- -beta/N * I * S
    dI <- beta/N * I * S - gamma * I
    dR <- gamma * I
    list(c(dS, dI, dR))
  })
}

init <- c(S = N-Infected[1], I = Infected[1], R = 0)
RSS <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- ode(y = init, times = Day, func = SIR, parms = parameters)
  fit <- out[ , 3]
  sum((Infected - fit)^2)
}

# optimize with some sensible conditions
Opt <- optim(c(0.5, 0.5), RSS, method = "L-BFGS-B", 
             lower = c(0, 0), upper = c(10, 10))
Opt$message

Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par
# beta     gamma 
# 0.6512503 0.4920399 

out <- ode(y = init, times = Day, func = SIR, parms = Opt_par)

plot(out)
plot(out, obs=data.frame(time=Day, I=Infected))

A basic SEIR model and additional links can be found here. Note however that not only the model is too simple, but also the data. The number of confirmed cases is not the number of infected, it depends on how it was tested and is probably much less than the number of infected.

Upvotes: 1

Related Questions