Em ciência, é importante ter o seu trabalho visto e revisto pelos pares. Então, achei de muito bom tom que os autores do modelo do Imperial College sobre infecção por Covid-19, que tem sido citado por um monte de gente, terem disponibilizados os códigos de R do exercício. Os códigos podem ser acessados aqui e o artigo pode ser lido aqui.
Na terça-feira, publiquei o nosso comentário de conjuntura semanal com um modelo SIR ajustado aos dados do Brasil. Para quem quiser replicar, todos os códigos de R foram disponibilizados, de modo que o post é totalmente reprodutível. Já aqui, quero apresentar as curvas geradas pelo modelo, considerando um horizonte de 90 dias.
O modelo divide a população em compartimentos, como ilustra a figura a seguir.
As pessoas são suscetíveis à doença, depois infectadas e, por fim, recuperadas. Ao deixar a doença tomar seu curso, sem fazer nada a respeito, o total de infectados atingiria seu pico no início de maio, se reduzindo ao longo das semanas seguintes. O número do pico fica pouco abaixo de 24 milhões de pessoas, com todas as consequências conhecidas sobre o sistema de saúde.
Como o país está tomando medidas de distanciamento social, o pico de infectados deve ocorrer antes disso, em meados de abril e com um número menor.
(*) Aprenda a rodar modelos como esse ao aprender R em nossos Cursos Aplicados de R.
Ao longo dos últimos dias, tenho publicado nesse espaço alguns posts e exercícios sobre a pandemia do coronavírus. Ontem, a propósito, publiquei o comentário de conjuntura dessa semana com um modelo SIR ajustado aos dados da doença no Brasil. Como é possível observar nesse exercício, o Brasil está no início da transmissão, com um crescimento exponencial dos casos confirmados.
Se olharmos, contudo, para os dados da China, primeiro país exposto à pandemia, a curva de casos confirmados parece seguir um formato logístico. Isso, a propósito, está em linha com a tese de "imunidade de grupo", ou seja, quanto mais pessoas vão sendo expostas ao vírus, mais pessoas ganham imunidade e a contaminação passa a desacelerar.
As curvas em formato de sino que têm sido divulgadas por aí, nesse aspecto, derivam justamente do modelo SIR, onde as pessoas são "compartimentadas" nos grupos de suscetíveis, infectados e recuperados. Ou seja, as pessoas saem de um para outro grupo, daí o formato da curva.
O formato logístico, por seu turno, não significa que devemos simplesmente abandonar as medidas de distanciamento social. Isso porque, quanto mais pessoas forem expostas ao vírus, mais casos graves serão registrados, o que tende a congestionar o sistema de saúde, como temos visto na Itália.
A pergunta de um trilhão de reais, portanto, é onde é o "limite superior" da curva logística.
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 blogLearning 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.
## 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)
}
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.