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