A seguir apresentamos o código para obter a deconvolução para Ouro Branco. O objetivo é estimar a curva de incidência a partir da curva de notificações, considerando uma amostra de indivíduos que reportou ambas as datas. Neste script a cidade selecionada foi Ouro Branco. Para rodar para outras, basta fazer as devidas adaptações na leitura dos dados.
# Por sorte a salvamos antes
data <- read.csv("INFLUD-15-02-2021.csv", sep = ";")
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(fitdistrplus)
## Carregando pacotes exigidos: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## Carregando pacotes exigidos: survival
library(utils)
library(httr)
library(ggplot2)
library(EpiEstim)
library(incidental)
##### OB
# Filtrando dados apenas de OB
OB <- filter(data, ID_MUNICIP == "OURO BRANCO")
# Armazenando data de notificação
OB_dados <- as.data.frame(OB$ID_REGIONA)
OB_dados[, "DataNot"] <- as.data.frame(OB$DT_NOTIFIC)
# Armazenando data de ínico dos sintomas
OB_dados[, "DataInicio"] <- as.data.frame(OB$DT_SIN_PRI)
# Passando dados para o formato de data
OB_dados$DataNot <- as.Date(OB_dados$DataNot, "%d/%m/%y")
OB_dados$DataInicio <- as.Date(OB_dados$DataInicio, "%d/%m/%y")
# Criando coluna com atraso nos dias de notificação
OB_dados[, "OB_Atraso"] <- difftime(OB_dados$DataNot, OB_dados$DataInicio, units = "days")
# Armazenado dados (em valor absoluto) do atraso em um vetor
OB_Atraso <- abs(as.vector(OB_dados$OB_Atraso))
# Retirando valores negativos devido ao erro de ano na notificação dos dados
OB_Atraso2 <- OB_Atraso>0
OB_Atraso <- OB_Atraso[OB_Atraso2]
# Retirando valores maiores que 100 (prováveis de erros de digitação)
OB_Atraso3 <- OB_Atraso<100
OB_Atraso <- OB_Atraso[OB_Atraso3]
# Observando algumas estatísticas descritivas para auxiliar a tomada de decisão (Cullen and Frey graph):
OB_g2 <- descdist(OB_Atraso, discrete=FALSE, boot=500)
# Ajustando às distribuições escolhidas:
fit_w_OB <- fitdist(OB_Atraso, "weibull")
fit_g_OB <- fitdist(OB_Atraso, "gamma")
fit_ln_OB <- fitdist(OB_Atraso, "lnorm")
## Graficos no ggplot para artigo
library(ggplot2)
p1 <- ggplot(data.frame(OB_Atraso), aes(x=OB_Atraso)) +
geom_histogram(binwidth = 1, bins= 55, boundary = 0.5, closed = "left", aes(y=..density..),
color="black", fill="grey90") + theme_bw()
p1 +
stat_function(fun = dweibull,
args = list(shape = fit_w_OB$estimate[1], scale = fit_w_OB$estimate[2]),
aes(colour = "Weibull", linetype = "Weibull"), size=1) +
stat_function(fun = dgamma,
args = list(rate = fit_g_OB$estimate[2], shape = fit_g_OB$estimate[1]),
aes(colour = "Gamma", linetype = "Gamma"), size=1) +
stat_function(fun = dlnorm, args = list(meanlog = fit_ln_OB$estimate[1], sdlog = fit_ln_OB$estimate[2]),
aes(colour = "Lognormal", linetype = "Lognormal"), size=1) +
xlab("Observações de atraso - OB") + ylab("Densidade") +
scale_color_manual(name = "Distribuição",
values = c("Weibull" = "green3",
"Gamma" = "red",
"Lognormal" = "blue")) +
scale_linetype_manual(name = "Distribuição",
values = c("Weibull" = "solid",
"Gamma" = "dashed",
"Lognormal" = "dotted")) +
theme(legend.position=c(0.7,0.7)) + xlim(0,40)
## Warning: Removed 1 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
OB_Atraso_df <- data.frame(OB_Atraso)
A partir da distribuição do atraso realiza-se a deconvolução para estimar a curva de indivíduos com início de manifestação dos sintomas.
###*******************************************************************************###
# Lendo planilha com os dados
dados_OB <- read.csv("DadosCOVID_OuroBranco.csv")
# Organizando dados
Infectados_OB <- dados_OB$Infectados
Mortos_OB <- dados_OB$Mortos
Recuperados_OB <- dados_OB$Recuperados
Dia <- 1:(length(Infectados_OB))
# Doentes <- Infectados - Recuperados - Mortos
Doentes_OB <- Infectados_OB - Mortos_OB - Recuperados_OB
### Deconvolution_OB
OB_model <- fit_incidence(reported = Doentes_OB, #round(yhat.sp1[,1]),
delay_dist = hist(OB_Atraso,
breaks = 55,
plot = F)$density,
regularization_order = 1,
dof_grid = 150)
### Convolucao
### dados de incidencia estimados via deconvolucao
a <- OB_model$Ihat
### probabilidades amostrais (depois fazer com weibull)
B <- hist(OB_Atraso,
breaks = 55,
plot = F)$density
b <- rep(0, length(a))
b[1:length(B)] <- B
### convolucao para obter dados reportados
yt <- convolve(OB_model$Ihat,rev(b), type = "o")
B2 <- B[1:40] ### probabilidades de atraso
plot(OB_model$Chat, pch = 21, type="n", xlab = "Tempo", ylab = "Doentes", ylim = c(0,130))
grid()
for(i in 1:100){
Doentes_OB_rep <- rep(0, length(Doentes_OB))
for(t in 1: length(Doentes_OB)){
atrasos <- sample(1:40, OB_model$Ihat[t], replace = T, prob = B2)
count_atraso <- tabulate(atrasos) # hist(atrasos, breaks = 60, plot = F)$counts
l_atr <- length(count_atraso)
Doentes_OB_rep[t:(t+l_atr-1)] <- Doentes_OB_rep[t:(t+l_atr-1)] + count_atraso
}
lines(Doentes_OB_rep, col = "gray85")
}
points(Doentes_OB, pch = 3, col = "red")
points(OB_model$Ihat, type = "l", col = "black", lty = 1, lwd = 2)
points(OB_model$Chat, type = "l", col = "blue", lty = 2, lwd = 2)
legend("topleft", legend = c("Incidencia estimada", "Convolução", "Dados reportados"),
col = c("black", "blue", "red"), lty = c(1,2,NA), pch = c(NA,NA,3), cex = 0.8)
Após obtenção da curva de indivíduos com manifestação dos sintomas, faz-se a estimativa de Rt corrigido.
# Estimando Rt para dados brutos
res_OB <- estimate_R(Doentes_OB,
method="parametric_si",
config = make_config(
list(mean_si = 4.8,
std_si = 2.3)))
## Default config will estimate R on weekly sliding windows.
## To change this change the t_start and t_end arguments.
R0_ic_dados_OB <- res_OB$R
colnames(R0_ic_dados_OB) <- c("t_start", "t_end", "mean_R", "std_R", "q0025", "q005", "q025", "median", "q075", "q095", "q0975")
### estimando Rt para dados corrigidos
res_OB2 <- estimate_R(OB_model$Ihat,
method="parametric_si",
config = make_config(
list(mean_si = 4.8,
std_si = 2.3)))
## Default config will estimate R on weekly sliding windows.
## To change this change the t_start and t_end arguments.
R0_ic_dados_OB_deconv <- res_OB2$R
colnames(R0_ic_dados_OB_deconv) <- c("t_start", "t_end", "mean_R", "std_R", "q0025", "q005", "q025", "median", "q075", "q095", "q0975")
line_types <- c("Rt corrigido" = "solid", "Rt bruto" = "dashed")
cols <- c("Rt corrigido" = "gray12", "Rt bruto" = "gray50")
p6 <- ggplot() +
geom_line(data = R0_ic_dados_OB_deconv, aes(x=t_start, y = mean_R, linetype = "Rt corrigido", colour = "Rt corrigido"), size=0.75) +
geom_line(data = R0_ic_dados_OB, aes(x=t_start, y = mean_R, linetype = "Rt bruto", colour = "Rt bruto"), size=0.75) +
geom_hline(yintercept = 1, color = "tomato2") +
scale_x_continuous(breaks = c(0, 30, 60, 90, 120, 150),
label = c("18/05","03/07","18/08","30/09","16/11","30/12")) +
ylab("Rt") +
xlab("Dias") +
scale_colour_manual(name=NULL,values=cols,
guide = guide_legend(override.aes=aes(fill=NA))) +
scale_linetype_manual(name=NULL,values=line_types) +
ggtitle("Ouro Branco") + theme_bw()
p6 + theme(legend.position=c(0.83,0.9), legend.key.size = unit(0.4, 'cm'), legend.key.width = unit(0.85, 'cm'))