Top 10 municípios da Bahia em casos de COVID-19
Cid Edson Póvoas
2020-04-16
Resultados
library(raster)
library(sf)
library(ggplot2)
library(DT)
library(ggrepel)
library(geobr)
library(leaflet)
library(leaflet.extras)
library(RColorBrewer)
#remotes::install_github("liibre/coronabr", force=T)
library(coronabr)
library(magrittr)
library(tidyverse)
muni <- read_municipality(code_muni= "BA", year=2017, showProgress = F)
sede <- read_municipal_seat(year= 2010, showProgress = F)
sede$longitude <- as.numeric(gsub(".*?([-]*[0-9]+[.][0-9]+).*", "\\1", sede$geom))
sede$latitude <- as.numeric(gsub(".* ([-]*[0-9]+[.][0-9]+).*", "\\1", sede$geom))
sede <- subset(sede, abbrev_state=="BA")
baz <- as.data.frame(sede[,c(2,4,9:10)])
baz <- baz[-5]
mapba <- merge(muni,baz)
options(DT.options = list(scrollY="100vh"))
dad <-get_corona_br()
ba <- subset(dad, state=="BA" & place_type == "city" & is_last == "True")
ba$city <- iconv(ba$city, "UTF-8")
ba <- subset(ba, city != "Importados/Indefinidos")
bahia.all <- ba %>% mutate(ranking = dense_rank(desc(confirmed)))
k <- 10
top.citys <- bahia.all %>% filter(ranking <= k) %>%
arrange(ranking) %>% pull(city) %>% as.character()
top.citys %<>% c('Outras')
data.bahia <- bahia.all %>% filter(!is.na(city)) %>%
mutate(city=ifelse(ranking <= k, as.character(city), 'Outras'))
top10 <- subset(data.bahia, city != "Outras")
top10$deaths[top10$deaths == 0] <- NA
max.date <- top10$date[1]
idx <- order(top10$confirmed, decreasing = T)
levels <- top10$city[idx]
top10$city <- factor(top10$city, levels=levels, ordered=TRUE)
top.10.long <- top10 %>% select(city,
confirmed,
confirmed_per_100k_inhabitants,
death_rate,
deaths) %>% mutate(death_rate = death_rate * 100) %>%
gather(key = type, value = count,-city)
top.10.long %<>% mutate(type=recode_factor(type,
confirmed='Casos Confirmados',
confirmed_per_100k_inhabitants='Por 100 mil habitantes',
deaths='Óbitos',
death_rate='Letalidade (%)'))
top.10.long %>% ggplot(aes(x=city, y=count, fill=city, group=city, label=count)) +
geom_bar(stat='identity') +
#geom_text(aes(label=count, y=count), size=21, vjust=-1) +
xlab('') + ylab('') +
geom_text(aes(label= round(count,1)),
vjust = -0.2,
angle=0,
color="black",
position = position_dodge(1),
size =2.5)+
labs(title=paste0('Top 10 municípios da Bahia em casos de COVID-19 - ', max.date),
caption = "Cid Póvoas / Fonte: Secretaria de Saúde da Bahia") +
scale_fill_discrete(name='', labels=aes(count)) +
theme_minimal()+
theme(legend.title=element_blank(),
legend.position='none',
plot.title=element_text(size=11),
axis.text=element_text(size=9),
axis.text.x=element_text(angle=45, hjust=1)) +
facet_wrap(~type, ncol=2, scales='free_y')
BAm <- as(mapba, "sf")
names(ba)[9] <- "code_muni"
bacorona = merge(BAm, ba, by = "code_muni")
bacorona <- st_transform(bacorona, 4326)
darkcols <- brewer.pal(8, "YlOrRd")
darkcols <- darkcols[2:8]
bins <- c(0, 2, 5, 15, 20, 30, 50, 100, 200)
pal_fun <- colorBin(darkcols, domain = bacorona$confirmed_per_100k_inhabitants, bins = bins)
bins2 <- c(1, 5, 10 ,20, 50, 100, 200, 500)
pal_fun2 <- colorBin(darkcols, domain = bacorona$confirmed, bins = bins2)
p_popup <- paste0("Casos/100mil hab. em ", bacorona$city, "<strong>: </strong>" , round(bacorona$confirmed_per_100k_inhabitants,1))
p_popup2 <- paste0("Número de casos em ", bacorona$city, "<strong>: </strong>" , round(bacorona$confirmed,0))Mapa de concentração de casos por 100 mil habitantes
Heatmap
leaflet(bacorona) %>%
addHeatmap(lng = ~as.numeric(bacorona$longitude), lat = ~as.numeric(bacorona$latitude), intensity = ~bacorona$confirmed_per_100k_inhabitants,
blur = 20, max = 10, radius = 15) %>% addTiles() Mapa de número de casos por 100 mil habitantes
leaflet(bacorona) %>%
addPolygons(
stroke = FALSE,
fillColor = ~ pal_fun(confirmed_per_100k_inhabitants),
fillOpacity = 5,
smoothFactor = 0.1,
popup = p_popup
) %>%
addTiles() %>%
addLegend("bottomright",
pal = pal_fun,
values = ~ confirmed_per_100k_inhabitants,
title = 'Casos/100mil hab.')Mapa de número de casos
leaflet(bacorona) %>%
addPolygons(
stroke = FALSE,
fillColor = ~ pal_fun2(confirmed),
fillOpacity = 5,
smoothFactor = 0.1,
popup = p_popup2
) %>%
addTiles() %>%
addLegend("bottomright",
pal = pal_fun2,
values = ~ confirmed,
title = 'Número de casos')Tabela
ba[c(-1:-2, -4, -7,-9)] %>% mutate(confirmed_per_100k_inhabitants =
confirmed_per_100k_inhabitants %>%
round(1),
death_rate = death_rate*100 %>% round(1)) %>%
datatable(
.,
rownames = FALSE,
filter = "top",
extensions = c('Buttons', 'Scroller'),
options = list(
dom = 'Blfrtip',
bPaginate = TRUE,
buttons = 'excel',
deferRender = T,
scroller = TRUE
),
colnames = c('Cidade' = 1, 'Confirmados' = 2, 'Por 100k' = 5, 'Letalidade (%)' = 6, 'Óbitos' = 3, 'População' = 4 )
)ba2 <- subset(dad, state == "BA" & place_type == "city")
ba2$city <- iconv(ba2$city, "UTF-8")
ba2 <- subset(ba2, city != "Importados/Indefinidos")
ba2 %<>% group_by(city) %>% mutate(cum_cases = cumsum(confirmed))
ba2.long <- ba2 %>%
select(c(city, date, confirmed_per_100k_inhabitants)) %>%
gather(key=type, value=count, -c(city, date))
ba2.long## # A tibble: 1,100 x 4
## # Groups: city [84]
## city date type count
## <chr> <date> <chr> <dbl>
## 1 Abaíra 2020-04-15 confirmed_per_100k_inhabitants 11.4
## 2 Adustina 2020-04-15 confirmed_per_100k_inhabitants 5.87
## 3 Alagoinhas 2020-04-15 confirmed_per_100k_inhabitants 3.30
## 4 Araci 2020-04-15 confirmed_per_100k_inhabitants 1.84
## 5 Aurelino Leal 2020-04-15 confirmed_per_100k_inhabitants 8.67
## 6 Barra 2020-04-15 confirmed_per_100k_inhabitants 3.73
## 7 Barra do Choça 2020-04-15 confirmed_per_100k_inhabitants 3.16
## 8 Barra do Rocha 2020-04-15 confirmed_per_100k_inhabitants 52.5
## 9 Barreiras 2020-04-15 confirmed_per_100k_inhabitants 0.643
## 10 Belmonte 2020-04-15 confirmed_per_100k_inhabitants 8.57
## # ... with 1,090 more rowsba2.long %<>% mutate(type=recode_factor(type,
confirmed_per_100k_inhabitants='Casos confirmados popr 100 mil habitantes',
deaths='Óbitos'))
gg_color_hue <- function(n) {
hues = seq(15, 375, length = n + 1)
hcl(h = hues, l = 65, c = 100)[1:n]
}
cgroup_cols = gg_color_hue(10)
ba2.long$city <- as.character(ba2.long$city)
focus_cn <- c("Salvador", "Ilhéus", "Itabuna", "Uruçuca", "Ipiaú")
ba2.long %>% filter(count > 1) %>%
mutate(days_elapsed = (date) - min(date),
end_label = ifelse(date == max(date), city, NA),
end_label = case_when(city %in% focus_cn ~ end_label,
TRUE ~ NA_character_),
cgroup = case_when(city %in% focus_cn ~ city,
F ~ "ZZOTHER")) %>%
ggplot(mapping = aes(x = days_elapsed, y = count,
color = cgroup, label = end_label,
group = city)) +
geom_line(size = 0.7) +
geom_text_repel(nudge_x = 0.75,
segment.color = NA) +
guides(color = FALSE) +
scale_color_manual(values = cgroup_cols) +
scale_y_continuous(labels = scales::comma_format(accuracy = 1),
breaks = 2^seq(1, 19, 1),
trans = "log2") +
labs(x = "Dias desde o primeiro caso confirmado",
y = "Número acumulado de casos relatados (log2)",
title = "Número acumulado de casos relatados de COVID-19 por 100 mil habitantes",
subtitle = paste("Secretaria de Saúde do Estado", format(max(ba2.long$date), "%A, %B %e, %Y")),
caption = "Cid Póvoas @cidedson / Data: Secretaria de Saúde da Bahia") +
theme_classic()
SIR Model
ios <- ba2 %>% filter(city == "Ilhéus") %>%
select(c(city, date, confirmed)) %>% arrange(date) %>%
gather(key=type, value=count, -c(city, date))
Infectados <- c(ios$count)
inicio <- min(ios$date)
dia <- 1:(length(Infectados))
Day <- inicio+dia
N <- 162327 # população de Ilhéus
old <- par(mfrow = c(1, 2))
plot(Day, Infectados, type = "b")
plot(Day, Infectados, log = "y")
abline(lm(log10(Infectados) ~ Day))
title("Confirmed Cases 2019-nCoV Ilhéus",
outer = TRUE,
line = -2)
par(mfrow = c(1, 2))
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))
})
}
library(deSolve)
init <- c(S = N - Infectados[1],
I = Infectados[1],
R = 0)
RSS <- function(parameters) {
names(parameters) <- c("beta", "gamma")
out <- ode(
y = init,
times = dia,
func = SIR,
parms = parameters
)
fit <- out[, 3]
sum((Infectados - fit) ^ 2)
}
Opt <-
optim(
c(0.5, 0.5),
RSS,
method = "L-BFGS-B",
lower = c(0, 0),
upper = c(1, 1)
)
#Opt$message
Opt_par <- setNames(Opt$par, c("beta", "gamma"))
#Opt_par
t <- 1:120
fit <-
data.frame(ode(
y = init,
times = t,
func = SIR,
parms = Opt_par
))
col <- 1:3
matplot(
fit$time,
fit[, 2:4],
type = "l",
xlab = "Day",
ylab = "Number of subjects",
lwd = 2,
lty = 1,
col = col
)
matplot(
fit$time,
fit[, 2:4],
type = "l",
xlab = "Day",
ylab = "Number of subjects",
lwd = 2,
lty = 1,
col = col,
log = "y"
)
points(dia, (Infectados))
legend(
"bottomright",
c("Susceptibles", "Infecteds", "Recovereds"),
lty = 1,
lwd = 2,
col = col,
inset = 0.05
)
title("SIR Model 2019-nCoV Ilhéus",
outer = TRUE,
line = -2)
par(old)
library(benford.analysis)
R0 <- setNames(Opt_par["beta"] / Opt_par["gamma"], "R0")
print(paste("R0 = ",round(as.numeric(R0),1)))## [1] "R0 = 1.5"print(paste("Número de infectados = ",(round(fit[fit$I == max(fit$I), "I", drop = FALSE]$I,0))))## [1] "Número de infectados = 9741"print(paste("Data do pico = ",(as.Date(inicio+as.numeric(rownames(fit[fit$I == max(fit$I), "I", drop = FALSE]))) %>% format('%d %B %Y') )))## [1] "Data do pico = 20 maio 2020"print(paste("Letalidade em Ilhéus = ",paste(round((2/max(ios$count)*100),2), "%")))## [1] "Letalidade em Ilhéus = 4.08 %"print(paste("Estimativa de mortos com atual letalidade = ",(round((max(fit$I) * (2/max(ios$count))),0))))## [1] "Estimativa de mortos com atual letalidade = 398"