Проблема с поиском параметров для модели SIR в R

Question

Я пытаюсь смоделировать скорость заражения коронавирусом, используя модель SIR в R. Однако мои значения бета (контролирует переход между Susceptibles и Infected) и gamma (контролирует переход между Infected и Recovered) значения неверны, как вы можете видеть .

      beta     gamma 
 1.0000000 0.8407238

Вот мой код:

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

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

Затем я нахожу 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

Затем я нахожу бета и гамма

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

Что пошло не так с моим кодом и как я могу это исправить? Заранее спасибо!

Answer 1

Оптимизация часто бывает сложной, есть несколько вещей, которые можно улучшить:

• модель очень проста, вы можете вместо нее использовать модель SEIR (будучи все еще слишком простой),
• расширить ограничения блока, это помогло в вашем случае,
• используйте другой оптимизатор, например, из мета-пакета, такого как FME , который также имеет дополнительную диагностику,
• перемасштабируйте переменные до более удобного диапазона, например, разделив совокупность на 1e3,
• ИЛИ, в качестве альтернативы, адаптируйте допуски (atol), соответствующие исходным шкалам.

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

Базовая c модель SEIR и дополнительные ссылки можно найти здесь . Однако обратите внимание, что не только модель слишком проста, но и данные. Количество подтвержденных случаев составляет , а не количество инфицированных, это зависит от того, как они были протестированы, и, вероятно, намного меньше, чем количество зараженных.