covid-19 e modelo sir Archives - Análise Macro

Nas últimas semanas, o mundo foi assolado pela maior pandemia das últimas décadas, com impactos severos sobre os mercados globais e sobre os sistemas de saúde de grande parte dos países. No Brasil, o primeiro caso foi registrado no último dia 26 de fevereiro e a curva de crescimento dos infectados segue uma trajetória exponencial. De forma a compreender o comportamento dos casos no país, nós vamos implementar o famoso modelo SIR proposto por Kermack e McKendrick (1927) aos dados de covid-19 no Brasil. Meu objetivo aqui é basicamente construir o gráfico abaixo com as previsões de casos confirmados a partir de um modelo SIR.

Para estimar o modelo SIR, estou basicamente replicando o código do blog Learning Machines, com algumas mudanças tópicas para visualização e previsão dos casos infectados nos próximos dias em nosso país. Os dados para o Covid-19 foram obtidos através do pacote nCov2019. É possível, diga-se, que existam dados mais atualizados no repositório da JHU ou na plataforma do Ministério da Saúde. Os demais pacotes utilizados são carregados a seguir.

require(nCov2019)
require(dplyr)
require(ggplot2)
require(scales)
require(gridExtra)
require(deSolve)

Os dados são obtidos a seguir.

## Obtendo os dados
data = load_nCov2019()
data_global = data["global"] #extract global data
data_br = filter(data_global, country == 'Brazil')

A seguir, visualizamos os dados.

Infected = data_br$cum_confirm
Day <- 1:(length(Infected))
N <- 211289547 # População do Brasil

df = tibble(Day, Infected, lnInfected = log(Infected))

g1 = ggplot(df, aes(x=Day, y=Infected))+
geom_line()+
geom_point()+
labs(x='Dias', y='Casos')

g2 = ggplot(df, aes(x=Day, y=lnInfected))+
geom_line()+
geom_point()+
geom_smooth(method='lm', se=FALSE)+
labs(x='Dias', y='log(Casos)')

grid.arrange(g1, g2,
ncol=2, nrow=1,
top='Total de Infectados por Covid-19 no Brasil')

O modelo SIR é um dos modelos epidemiológicos mais simples, e muitos modelos são derivações do mesmo. O modelo consiste em três compartimentos: S para o número de pessoas suscetíveis à doença, I para o número de infectados e R para o número de indivíduos recuperados (ou imunes).

Para modelar a dinâmica de uma epidemia, precisamos de três equações diferenciais, uma para a mudança em cada compartimento, com o parâmetro controlando a transição entre S e I e controlando a transição entre I e R:

(1)  

A seguir implementamos o código do blog Learning Machines para ajustar o modelo aos nossos dados:

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

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

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

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

Com o modelo estimado, nós podemos ver o ajuste do mesmo abaixo.

df2 = tibble(time = data_br$time,
observado = data_br$cum_confirm,
modelo = round(fit$I[1:length(Infected)],0),
dias = Day)

ggplot(df2, aes(x=time))+
geom_point(aes(y=observado, colour='Casos Confirmados'),
stat='identity')+
geom_line(aes(y=modelo, colour='Modelo'))+
scale_colour_manual('',
values=c('Casos Confirmados'='black',
'Modelo'='red'))+
scale_x_date(breaks = date_breaks("1 day"),
labels = date_format("%d/%b"))+
theme(axis.text.x=element_text(angle=45, hjust=1),
plot.title = element_text(size=12, face='bold'),
legend.position = c(.3,.6))+
labs(x='', y='Casos Confirmados',
title='Modelando casos de Covid-19 no Brasil',
subtitle='Modelo SIR aplicado ao Brasil',
caption='analisemacro.com.br')

Na sequência nós visualizamos os casos previstos para os próximos dias.

dates = seq(data_br$time[1], length.out = length(t), by='1 days')
df3 = tibble(time = dates,
observado = c(data_br$cum_confirm,
rep(NA,
length(fit$I)-length(data_br$cum_confirm))),
modelo = round(fit$I,0))

ggplot(df3, aes(x=time))+
annotate("rect", fill = "orange", alpha = 0.3,
xmin = as.Date('2020-03-24'),
xmax = as.Date('2020-03-29'),
ymin = -Inf, ymax = Inf)+
annotate('text', x=as.Date('2020-03-26'), y=1000,
label='Previsão',
colour='black', size=4.5)+
geom_point(aes(y=observado, colour='Casos Confirmados'),
stat='identity', size=3)+
geom_line(aes(y=modelo, colour='Modelo'),
linetype='dashed', size=.8)+
scale_colour_manual('',
values=c('Casos Confirmados'='black',
'Modelo'='red'))+
scale_x_date(breaks = date_breaks("1 day"),
labels = date_format("%d/%b"))+
theme(axis.text.x=element_text(angle=45, hjust=1),
plot.title = element_text(size=12, face='bold'),
legend.position = c(.3,.6))+
labs(x='', y='Casos Confirmados',
title='Projetando casos de Covid-19 no Brasil',
subtitle='Modelo SIR aplicado ao Brasil',
caption='Fonte: analisemacro.com.br')

Mantida a trajetória de crescimento, teremos quase 9 mil casos no próximo dia 29/03. Obviamente que essa tendência não leva em consideração as medidas de distanciamento social que visam reduzir a taxa de reprodução, isto é, quantas pessoas saudáveis uma pessoa infectada contamina, em média. Ela é obtida pela razão entre os parâmetros e . Para os dados disponibilizados até aqui, a taxa de reprodução está em 1,79 no Brasil, abaixo de 2 que tem sido a taxa obtida em outros países.

________________________

Kermack, William Ogilvy, and Anderson G. McKendrick. 1927. “A Contribution to the Mathematical Theory of Epidemics.” Proceedings of the Royal Society of London, Series A 115, no. 772: 700–721.

(*) Aprenda R nos nossos Cursos Aplicados de R.

(**) Um pdf com os códigos estará disponível amanhã no Clube do Código.

Assine nossa Newsletter e receba nosso conteúdo exclusivo