Modelo SIR e Previsão do Caso Brasileiro

Introdução

O avanço da Pandemia causada pelo vírus COVID-19 trouxe o interesse de pesquisadores de distintas áreas em mensurar o avanço e impactos da doença. Nosso interesse neste momento é conhecer um modelo epidemiológico simples que nos permita inferir a trajetória de crescimento da epidemia, bem como a necessidade de leitos que essa venha a apresentar para o país.

Isso posto, vamos utilizar o pacote, recém lançado no R, chamado nCov2019 para a obtenção da série histórica. Além desse, utilizaremos os pacotes dplyr, deSolve, xts e highcharter. O primeiro será utilizado para organizar os dados, o segundo para resolver equações diferenciais e o terceiro para criar séries temporais descontínuas e o último para plotar os gráficos.

library(nCov2019)
library(dplyr)
library(deSolve)
library(xts)
library(highcharter)
library(widgetframe)

### Obtenção dos Dados

Dados         <- load_nCov2019()
Dados_Brasil  <- Dados$global %>%
                 filter(country == 'Brazil')

Infectados   <- Dados_Brasil$cum_confirm
Dias         <- seq(1,length(Infectados),by =1)
N            <- 211305000 # População do Brasil segundo IBGE em 2020.

tail(Dados_Brasil)

##          time country cum_confirm cum_heal cum_dead
## 28 2020-03-24  Brazil        1965      150       34
## 29 2020-03-25  Brazil        2271      150       47
## 30 2020-03-26  Brazil        2587      150       61
## 31 2020-03-27  Brazil        2988      150       77
## 32 2020-03-28  Brazil        3477      150       93
## 33 2020-03-29  Brazil        3904      150      114

Na parte acima do Script, carregamos os pacotes que serão utilizados, utilizamos a função load_nCov2019() para baixar os dados mais recentes sobre reports de COVID-19 nos países. Observamos as últimas observações da série brasileira através da função tail, bastante conhecida.

O Modelo SIR (Susceptible, Infectious and Recovered)

O modelo SIR é composto por três EDOs, que representam: (a) Pessoas suscetíveis a infecção; (b) Pessoas Infectadas; (c) Pessoas recuperadas (imunizadas). Matematicamente, podem ser expressas pelo seguinte sistema: \[\begin{array}{l} \partial S(t)/\partial t=-(\beta/N).S(t).I(t)\\ \partial I(t)/\partial t=(\beta/N).I(t).S(t)-\gamma.I(t)\\ \partial R(t)/\partial t=\gamma.I(t) \end{array}\] Onde \(S(t)\) são os suscetíveis, \(I(t)\) População infectada, \(R(t)\) é o número de Recuperados (imunizados) do Modelo.

Na sequência definimos no R o sistema através da seguinte função:

### Utilizando o modelo SIR

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

### Valores para as variáveis estado
Inicial <- c(S = N-Infectados[1], I = Infectados[1], R = 0)

Iremos então resolver o sistema de EDOs (integrá-las) e estimar os parâmetros com base nos dados observados de infectados acumulados no Brasil. Para essa estimação, iremos fazer uso do algoritmo simplex de otimização com restrições para valores máximos e mínimos dos parâmetros. Na primeira parte, desenhamos uma função que gera a Soma dos Quadrados dos Resíduos e essa função é otimizada para obtenção dos parâmetros estimados.

SQR <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- ode(y = Inicial, times = Dias, func = SIR, parms = parameters)
  fit <- out[,3]
  sum((Infectados - fit)^2)
}

Opt <- optim(c(0.5, 0.5), SQR,
             method = "L-BFGS-B",
             lower = c(0, 0), 
             upper = c(1, 1))

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

Na sequência, fazemos a projeção da curva estimada tendo como base os parâmetros obtidos anteriormente.

t <- 1:120 # time in days
fit <- data.frame(ode(y = Inicial, times = t,
                      func = SIR, parms = Opt_par))

Datas = seq(Dados_Brasil$time[1], 
            length.out = length(t), 
            by='1 days')

Previsao_IF <- xts(fit$I, order.by = Datas)
Previsao_S  <- xts(fit$S, order.by = Datas)
Previsao_R  <- xts(fit$R, order.by = Datas)
Observado   <- xts(Dados_Brasil$cum_confirm, order.by = Dados_Brasil$time)


# Plot Highcharter
frameWidget(hchart(Previsao_IF, type = "line", color = "red", name = "Série Prevista") %>%
      hc_title(text = "Contaminados por COVID-19 e Previsão (Modelo SIR Estimado)") %>%
      hc_add_series(Observado, type = "scatter", color = "black", name = "Casos Confirmados") %>%
      hc_subtitle(text = "Dados Diários - Pacote nCovid19") %>%
      hc_tooltip(pointFormat = '({point.x: %Y-%m-%d}) 
                                 {point.y:.0f} Casos') %>%
      hc_xAxis(plotBands = list(
      list(label = list(text = "Previsão"),
      color = 'lightgrey',
      from = datetime_to_timestamp(Datas[length(Dias)]+1),
      to = datetime_to_timestamp(Datas[length(Dias)]+length(t))))))

frameWidget(hchart(Previsao_S, type = "line", color = "magenta", name = "Previsão Suscetíveis") %>%
  hc_title(text = "Dinâmica SIR entre Variáveis Suscetíveis e Recuperados") %>%
  hc_add_series(Previsao_R, type = "line", color = "blue", name = "Previsão Recuperados") %>%
  hc_tooltip(pointFormat = '({point.x: %Y-%m-%d}) 
                                 {point.y:.0f} Casos'))

A figura anterior, relacionda a curva de infectados, aponta uma informação importante. O pico da curva se dá em 30/04. Neste momento, levando em consideração que temos 20% dos casos necessitando internação hospitalar, teríamos a necessidade de cerca de 1 milhão e 600 mil leitos no dia. Certamente, nesse número sofreremos da escassez de leitos e o sistema não dá vazão a sua demanda. Por conta disso, a paralização é um instrumento utilizado no sentido de deixar a curva mais plana. Somente assim, o sistema daria conta de atender a demanda que surge.