Nowcasting de óbitos
a partir dos boletins epidemiológicos do Ministério da Saúde
Observatório Covid-19 BR
08/06/2020 13:55:55
Boletins do Ministério da Saúde
A partir dos boletins epidemiológicos do ministério da saúde compilamos os óbitos por SRAG com confirmação de Covid-19, cada boletim traz um gráfico de barras que conta os óbitos por data de óbito, ou seja, a data da ocorrência do óbito. É importante deixar claro que não há somete este modo de informar óbitos, o ministério mesmo já informou em alguns de seus boletins óbitos via a data de notificação do óbito. Todo o objetivo do nowcasting é realizar uma correçãopara se chegar à um número mais próximo da realidade, sem que se tenha que esperar para isso a chegada de dados no sistema. O nowcasting corrige esses atrasos que são inerentes do sistema de notificação.
uol<-read_csv("SRAGs-tabela-last-updated_revised3.csv")## Parsed with column specification:
## cols(
## .default = col_double(),
## Data = col_character()
## )## See spec(...) for full column specifications.uol<-as.data.frame(uol)
uol$Data<-as.Date(uol$Data, format = "%d/%m/%Y")
uol_melted<-reshape::melt(uol, id.vars = "Data")
p.uol <-
ggplot(uol_melted, aes(x = Data, y = value, col = variable)) +
geom_point(shape = 1)+
geom_line()+
geom_vline(xintercept = as.Date("2020-03-17", format = "%Y-%m-%d"), colour = "indianred3", size = 0.45, linetype = "dashed")+
scale_color_viridis_d(name = "Data Boletim", option = "viridis", direction = 1)+
labs(x = "Data", y = "Número de Óbitos") +
theme_bw() +
theme(legend.position = "right")+
theme(axis.text= element_text(size=14),
axis.title = element_text(size=14))
p.uol## Warning: Removed 1275 rows containing missing values (geom_point).## Warning: Removed 1275 row(s) containing missing values (geom_path).
O efeito dos atrasos nos boletins é o deslocamento do pico para a direita, mais perto do presente, e aumento do número de óbitos até pra datas do passado. A linha pontilhada vermelha marca o dia de anúncio do primeiro óbito, veja-se que há óbitos ocorridos antes dessa data, os boletins mais recentes realizaram essa correção. Outro efeito é no fim da série, em que há uma queda abrupta, mais perto do presente. Como se demora algum tempo para que o óbito já ocorrido seja confirmado como devido à Covid-19, os dias mais pertos do presente sempre terão menos óbitos que a atualidade. Em resumo, os óbitos não constam instantaneamente no sistema de notificação.
Montando as matrizes de datas
Com a tabela de boletins do MS, montamos a seguinte matriz:
uol <- read.csv("SRAGs-tabela-last-updated_revised3.csv") ##load da CSV Utilizada
uol2<-uol #variavel auxiliar
uol2[is.na(uol2)] <- 0 ##preenchendo com 0 onde não há dados ###
uol2.1<-uol2 ###variavel auxiliar
###loop pra contar os valores do boletins mais antigo e contabilizar esse valor nos boletins mais recentes###
###para manter a data de divulgação como a data do boletim em que a morte é contada pela primeira vez #######
for(i in 3:length(uol2)){
for(k in 1:dim(uol2)[1]){
if(uol2.1[k,i]>uol2.1[k,(i-1)] & uol2.1[k,i]>0){
uol2.1[k,i]=uol2.1[k,(i-1)]
}
}
}
uol2<-uol2.1 ###devolvendo o resultado na variável utilizada ###
## matriz para guardar os resultados
uol3 <- matrix( nrow=nrow(uol2), ncol=ncol(uol2)-1)
## Coloca a ultima coluna dos dados na ultima coluna da matriz
uol3[,ncol(uol3)] <- uol2[,ncol(uol2)]
## Loop que vai preenchendo a nova matriz com a difreença da coluna mais a direita pra seguinte
for(i in ncol(uol2):3){
uol3[,(i-2)] <- uol2[,(i-1)] - uol2[,i]
}
## Conferindo
all(apply(uol3, 1, sum) == uol2[,2]) ## ok!## [1] TRUE# print(head(uol3))
head(uol3)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 2 0 0 1 0 0 0 0 0 0 0 0
## [5,] 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 1 0 0 0 0 1 0 1 0 0 0 1 1 0 1 0 0
## [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 1 0 0 2
## [3,] 0 0 1 0 0 0 0 0 0 0 0 0 0 0 2
## [4,] 0 0 0 0 0 0 1 0 0 0 0 0 0 0 4
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0 1 0 6
## [6,] 0 1 1 0 0 0 1 0 0 0 0 1 1 0 6As colunas mais a esquerda são os boletins mais recentes. As linhas são as data em que ocorreram óbitos, mais a abaixo mais perto do presente.
## Junta com coluna de datas dos boletins e nomeia as colunas
uol3.df <- as.data.frame(uol3)
names(uol3.df) <- names(uol2)[-1]
uol4 <- cbind(data=uol2[,1], uol3.df)
## Monta data.frame de datas de obito e de regitro
datas.boletins <- as.character(as.Date(names(uol3.df), "X%d.%m.%y"))
datas.obito <- as.character(as.Date(as.character(uol4$data), "%d/%m/%y"))
## Vetores para guardar as datas de evento e
onset.dates <- record.dates <- c()
## Loop sobre as linhas e então colunas da matriz
for(i in 1:length(datas.obito)){
for(j in 1:length(datas.boletins)){
if(uol3[i,j]>0){
onset.dates <- c(onset.dates, rep(datas.obito[i], uol3[i,j]))
record.dates <- c(record.dates, rep(datas.boletins[j], uol3[i,j]))
}
}
}
## Conferindo o total de óbitos ##
length(record.dates)## [1] 26888length(onset.dates)## [1] 26888sum(uol2[,2])## [1] 26888Vamos percorrendo a matriz de modo que os óbitos para cada data são atualizados conforme se percorre os boletins, colunas. Com isso criamos vetores de datas, teremos que ter dois vetores de tamanho igual a soma da segunda (teste no fim do chunk) coluna da matriz original, boletim mais recente. Assim teremos para cada óbito um par de datas, a data de óbito e a data de reporte do óbito, que aqui tomamos com a data do boletim. Como não temos uma data em que realmente o óbito foi notificado no sistema sivep-gripe, tomamos a data de notificação como as datas dos boletins. Conforme novos óbitos são notificados para um mesmo dia, somente somamos a esse dias esse novo óbito e lhe damos a data de notificação conforme a data do boletim que ele começou a constar.
## Monta o data.frame
uol_df <- data.frame(uol_death_date = as.Date(onset.dates), uol_report_date = as.Date(record.dates))
## Transformando em datas que entram no nowcasting##
uol_df = uol_df %>%
mutate(Death_date = as.Date(uol_df$uol_death_date, format = "%d/%m/%Y"),
Report_date = as.Date(uol_df$uol_report_date, format = "%d/%m/%Y")) %>%
as.data.frame()
### somando os óbitos ###
uol_df2 = uol_df %>%
group_by(Death_date)%>%
dplyr::summarise(N=n())%>%
as.data.frame()
print(head(uol_df))## uol_death_date uol_report_date Death_date Report_date
## 1 2020-03-15 2020-04-02 2020-03-15 2020-04-02
## 2 2020-03-16 2020-04-04 2020-03-16 2020-04-04
## 3 2020-03-16 2020-03-29 2020-03-16 2020-03-29
## 4 2020-03-16 2020-03-29 2020-03-16 2020-03-29
## 5 2020-03-17 2020-05-12 2020-03-17 2020-05-12
## 6 2020-03-17 2020-04-15 2020-03-17 2020-04-15print(head(uol_df2))## Death_date N
## 1 2020-03-15 1
## 2 2020-03-16 3
## 3 2020-03-17 4
## 4 2020-03-18 8
## 5 2020-03-19 9
## 6 2020-03-20 17Temos dois dataframes, um com o total de óbitos em cada data, e outro com o par de datas para cada óbito ocorrido.
Nowcasting
A função NobBS, recebe o par de datas, utilizamos como unidade de tempo 1 dia. Os specs são o default, eles dão os priors que serão ajustados aos atrasos e computarão uma distribuição para eles. O atributo now recebe a última data do vetor de datas de ocorência de eventos. Os atributos onset_date e report_date, recebem respectivamente a coluna das datas de ocorrência do evento, e a data de notificação do evento.
###Nowcasting ###
nowcasting<- NobBS(data = uol_df,
now = max(uol_df$Death_date),
onset_date = "Death_date",
report_date = "Report_date",
units = "1 day",
specs = list(nAdapt = 3000, nBurnin = 3000, nThin = 1, nSamp = 10000)
)## Computing a nowcast for 2020-06-03## Warning in NobBS(data = uol_df, now = max(uol_df$Death_date), onset_date = "Death_date", : Warning!
## The line list has zero case reports for one or more possible onset dates at one or more delays.
## Proceeding under the assumption that the true number of cases at the associated delay(s) and
## week(s) is zero.## Warning in jags.model(file = ifelse(specs[["dist"]] == "Poisson", JAGSmodPois, : Adaptation
## incomplete## NOTE: Stopping adaptationbetas<-beta.summary(nowcasting) #### função em funcoes.RPlots
Plots distribuição de atrasos e Nowcasting dia a dia
################################################################################
## Plots: objetos ggplot2
################################################################################
## Tempos de atraso ##
p.betas <-
ggplot(betas, aes(atraso, mean)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.25) +
xlab("Dias após dia do óbito") +
ylab("Probabilidade de notificação") +
theme_bw() +
theme(legend.position="none")+
ggtitle("Atraso de notificação de óbito")
## N de casos previstos e seus ICS ##
p.prev.ic <- ggplot(nowcasting$estimates, aes(x = onset_date, y = estimate)) +
geom_line(data = uol_df2, aes(x = Death_date, y = N),
col = "skyblue", lwd = 1.2) +
geom_line(aes(x = onset_date, y = n.reported), col = "blue", lwd = 1.2)+
geom_line(col = "indianred3") +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha =0.15, fill = "indianred3") +
xlab("Dia do Óbito") +
ylab("Nº de Óbitos") +
theme_bw() +
theme(legend.position = "none") +
ggtitle("Nowcasting de óbitos")
ggarrange(p.betas, p.prev.ic)## Warning: Removed 5 row(s) containing missing values (geom_path).
Gráfico da esquerda, é a distrubuição de probabilidades de atraso, quanto mais alto maior a probabilidade de sofrer um atraso de até “X” dias. Gráfico da direita é o nowcasting para a série de óbitos. Linha vermelha e seus tons, são o nowcasting por dia e o intervalo de confiança respectivo. Linha azul é a série temporal de óbitos que o ministério divulga.
Plot do Nowcasting para números acumulados
uol_final<-nowcasting$estimates
uol_final<-uol_final[, c("estimate", "lower", "upper")]
uol_final2<-as.data.frame(cbind(uol_final,
"Data" = as.Date(uol_df2$Death_date, "%Y-%m-%d"),
"Boletim ultimo" = uol[,2]))
uol_final3<-apply(t(uol_final2[,-4]), 1, cumsum)
colnames(uol_final3)<-c("estimate Cumsum", "lower Cumsum", "upper Cumsum", "Boletim Cumsum")
uol_final4<-as.data.frame(cbind("Data" = uol_final2$Data, uol_final2[,-4], uol_final3))
p.prev.ic.cumsum <- ggplot(uol_final4, aes(x = Data, y = `estimate Cumsum`)) +
geom_line(data = uol_final4, aes(x = Data, y = `Boletim Cumsum`, color="Notificados"), lwd = 1.5) +
geom_line(aes(col = "Estimado")) +
geom_ribbon(aes(ymin =`lower Cumsum`, ymax = `upper Cumsum`), fill="red", alpha =0.15) +
xlab("Dia do Óbito") +
ylab("Nº de Óbitos Acumulados") +
theme_bw() +
theme(legend.position = c(0.2,0.8), legend.title= element_blank()) +
scale_colour_manual(values = c("red", "blue"), aesthetics = c("colour", "fill")) +
ggtitle("Nowcasting de óbitos de COVID-19 cumulativo")
p.prev.ic.cumsum
Gráfico com o nowcasting porém em números cumulativos.
Referências
Bastos, L.S., Economou, T., Gomes, M.F., Villela, D.A., Coelho, F.C., Cruz, O.G., Stoner, O., Bailey, T. and Codeço, C.T. (2019). A modelling approach for correcting reporting delays in disease surveillance data. Statistics in medicine, 38(22), pp.4363-4377.
Sarah F. McGough, Michael A. Johansson, Marc Lipsitch, Nicolas A. Menzies(2019). Nowcasting by Bayesian Smoothing: A flexible, generalizable model for real-time epidemic tracking. bioRxiv 663823; doi: https://doi.org/10.1101/663823
Sarah McGough, Nicolas Menzies, Marc Lipsitch and Michael Johansson (2020). NobBS: Nowcasting by Bayesian Smoothing. R package version 0.1.0. https://CRAN.R-project.org/package=NobBS