Estimativa de R0 para MG

A seguir apresentamos o código para a estimativa de \(R_0\) para o estado de Minas Gerais.

# Bibliotecas
library(utils)
library(httr)
library(ggplot2)
library(bbmle)
## Carregando pacotes exigidos: stats4
#download dos dados
GET("https://raw.githubusercontent.com/wcota/covid19br/master/cases-brazil-states.csv", authenticate(":", ":", type="ntlm"), write_disk(tf <- tempfile(fileext = ".csv")))
## Response [https://raw.githubusercontent.com/wcota/covid19br/master/cases-brazil-states.csv]
##   Date: 2022-04-26 14:36
##   Status: 200
##   Content-Type: text/plain; charset=utf-8
##   Size: 3.31 MB
## <ON DISK>  D:\Usuários\Temp\RtmpiW2xR3\file46d821fce1d.csv
#leitura de dados
data <- read.csv(tf)

# Mudar para o estado que deseja modelar
estado <- c("MG")

# Filtrando para o estado selecionado
data_UF <- subset(data, state == estado)

# Ordenando por data
data_UF <- data_UF[order(as.Date(data_UF$date,format="%Y/%m/%d")),]

# Removendo dias sem infectados
data_UF <- data_UF[data_UF$totalCases != 0,]

# Encontrando R0
Dia_log <- 5:16

fit <- lm(log(data_UF$totalCases[5:16])~Dia_log)
summary(fit)
## 
## Call:
## lm(formula = log(data_UF$totalCases[5:16]) ~ Dia_log)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.269383 -0.070103  0.000161  0.121034  0.232379 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.53075    0.15995   -9.57 2.37e-06 ***
## Dia_log      0.39830    0.01447   27.52 9.29e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.173 on 10 degrees of freedom
## Multiple R-squared:  0.987,  Adjusted R-squared:  0.9857 
## F-statistic: 757.6 on 1 and 10 DF,  p-value: 9.287e-11
#
r <- fit$coef["Dia_log"]

# V, intervalo serial, conforme Nishiura et al. (2020), V = 4.8 
V <- 4.8

# R0 de MG
R0_MG <- V*r+1
R0_MG
##  Dia_log 
## 2.911855
# Limites p/ R0 de MG
V*confint(fit)[2,1] + 1
## [1] 2.757085
V*confint(fit)[2,2] + 1
## [1] 3.066624
# Plotando
data_UF3 <- data.frame(Day = 1:100,
                       totalCases = data_UF$totalCases[1:100])
data_UF3_fit <- data.frame(Day = Dia_log, 
                           totalCases = exp(predict(fit, newdata = data.frame(Dia_log))),
                           L = exp(predict(fit, newdata = data.frame(Dia_log), interval = "confidence")[,2]),
                           U = exp(predict(fit, newdata = data.frame(Dia_log), interval = "confidence")[,3]))

ggplot(data_UF3) + 
  geom_line(data = data_UF3_fit, aes(x=Day, y=totalCases)) +
  geom_ribbon(data = data_UF3_fit, aes(x = Day, ymin=L, ymax=U), alpha=0.2, 
              col = "aquamarine3", fill = "aquamarine3") + 
  ylab("log Incidência acumulada - MG") + xlab("Dia") + 
  geom_point(aes(x=Day, y=totalCases), col = "aquamarine4", size = 2, 
             shape = 21, fill = "aquamarine3") + 
  scale_y_continuous(trans='log10') + 
  theme_bw()