#instale apenas caso você ainda não tenha o pacote ----
# DICA: só é preciso instalar cada pacote uma única vez
{
library(microdatasus)
library(ggplot2)
library(tidyverse)
library(pdftools)
library(stringr)
library(ggpubr)
library(ggridges)
}O sarampo e a Vacinação no Brasil
Avaliação do impacto das políticas públicas de imunização e a incidência de casos e óbitos de Sarampo no Brasil
O Sarampo e a vacinação no Brasil
Carregando pacotes
Dados da análise
Baixar dados organizados
#conferir diretório de trabalho ----
getwd()
# caso não seja o mesmo dir onde está o seu código, defina corretamente com setwd()
# criar diretórios para os dados e para os resultados ----
dir.create(path = "arquivos")
dir.create(path = "results")
# baixar dados de entrada ----
download.file(
"https://raw.githubusercontent.com/heronoh/bioinfo/main/aulas/r/arquivos/sarampo--todos_dados.xlsx",
destfile = "arquivos/sarampo--todos_dados.xlsx",
mode = "wb",
method = "libcurl"
)Tabela de siglas
# ler dados de entrada ----
## siglas ----
sigla_reg <- readxl::read_xlsx(path = "arquivos/sarampo--todos_dados.xlsx",
sheet = "Siglas",
col_names = T,
.name_repair = "unique") %>%
tibble::tibble()
sigla_reg %>%
print(n=100)# A tibble: 27 × 3
Região `Unidade da Federação` Sigla
<chr> <chr> <chr>
1 "Região\\nCentro-Oeste" Distrito Federal DF
2 "Região\\nCentro-Oeste" Goiás GO
3 "Região\\nCentro-Oeste" Mato Grosso MT
4 "Região\\nCentro-Oeste" Mato Grosso do Sul MS
5 "Região\\nNordeste" Alagoas AL
6 "Região\\nNordeste" Bahia BA
7 "Região\\nNordeste" Ceará CE
8 "Região\\nNordeste" Maranhão MA
9 "Região\\nNordeste" Paraíba PB
10 "Região\\nNordeste" Pernambuco PE
11 "Região\\nNordeste" Piauí PI
12 "Região\\nNordeste" Rio Grande do Norte RN
13 "Região\\nNordeste" Sergipe SE
14 "Região\\nNorte" Acre AC
15 "Região\\nNorte" Amapá AM
16 "Região\\nNorte" Amazonas AP
17 "Região\\nNorte" Pará PA
18 "Região\\nNorte" Rondônia RO
19 "Região\\nNorte" Roraima RR
20 "Região\\nNorte" Tocantins TO
21 "Região\\nSudeste" Espírito Santo ES
22 "Região\\nSudeste" Minas Gerais MG
23 "Região\\nSudeste" Rio de Janeiro RJ
24 "Região\\nSudeste" São Paulo SP
25 "Região\\nSul" Paraná PR
26 "Região\\nSul" Rio Grande do Sul RS
27 "Região\\nSul" Santa Catarina SC Habitantes por estado e ano
## habitantes ----
habitantes_br <- readxl::read_xlsx(path = "arquivos/sarampo--todos_dados.xlsx",
sheet = "Habitantes",
col_names = T,
.name_repair = "unique") %>%
tibble::tibble() %>%
rename_with(~ paste0("Ano ",
sub("^(\\d{4})$", "\\1", .x)),
matches("^\\d{4}$")) %>%
tidyr::pivot_longer(cols = dplyr::starts_with("Ano "),
names_to = "Ano",
values_to = "Value") %>%
dplyr::mutate(Ano = str_remove_all(Ano, pattern = "Ano ")) %>%
dplyr::left_join(y = BiocGenerics::unique(sigla_reg),
by = "Unidade da Federação") %>%
dplyr::mutate("Categoria" = "Habitantes")
habitantes_br# A tibble: 621 × 6
`Unidade da Federação` Ano Value Região Sigla Categoria
<chr> <chr> <dbl> <chr> <chr> <chr>
1 Distrito Federal 2001 2097447 "Região\\nCentro-Oeste" DF Habitantes
2 Distrito Federal 2002 2145839 "Região\\nCentro-Oeste" DF Habitantes
3 Distrito Federal 2003 2189789 "Região\\nCentro-Oeste" DF Habitantes
4 Distrito Federal 2004 2282049 "Região\\nCentro-Oeste" DF Habitantes
5 Distrito Federal 2005 2333108 "Região\\nCentro-Oeste" DF Habitantes
6 Distrito Federal 2006 2383784 "Região\\nCentro-Oeste" DF Habitantes
7 Distrito Federal 2008 2557159 "Região\\nCentro-Oeste" DF Habitantes
8 Distrito Federal 2009 2606885 "Região\\nCentro-Oeste" DF Habitantes
9 Distrito Federal 2010 2570160 "Região\\nCentro-Oeste" DF Habitantes
10 Distrito Federal 2011 2609998 "Região\\nCentro-Oeste" DF Habitantes
# ℹ 611 more rowsDoses vacinais por estado e ano
## doses aplicadas ----
doses_sarampo <- readxl::read_xlsx(path = "arquivos/sarampo--todos_dados.xlsx",
sheet = "Doses") %>%
rename_with(~ str_remove(.x, "\\.0$"), matches("^\\d{4}\\.0$")) %>%
dplyr::mutate(
dplyr::across(
dplyr::matches("^\\d{4}$"),
~ as.numeric(stringr::str_replace_all(as.character(.), ",", "."))
)
) %>%
tidyr::pivot_longer(cols = dplyr::starts_with(c("1","2")),
names_to = "Ano",
values_to = "Value") %>%
dplyr::left_join(y = BiocGenerics::unique(sigla_reg),
by = "Unidade da Federação") %>%
dplyr::mutate("Categoria" = "num doses")
doses_sarampo# A tibble: 783 × 6
`Unidade da Federação` Ano Value Região Sigla Categoria
<chr> <chr> <dbl> <chr> <chr> <chr>
1 Acre 1994 247271 "Região\\nNorte" AC num doses
2 Acre 1995 532068 "Região\\nNorte" AC num doses
3 Acre 1996 381238 "Região\\nNorte" AC num doses
4 Acre 1997 750997 "Região\\nNorte" AC num doses
5 Acre 1998 613288 "Região\\nNorte" AC num doses
6 Acre 1999 2095997 "Região\\nNorte" AC num doses
7 Acre 2000 897103 "Região\\nNorte" AC num doses
8 Acre 2001 803376 "Região\\nNorte" AC num doses
9 Acre 2002 730728 "Região\\nNorte" AC num doses
10 Acre 2003 721201 "Região\\nNorte" AC num doses
# ℹ 773 more rows1ª Dose - Cobertura vacinal por estado e ano
## cobertura 1a dose ----
cov_1a_dose <- readxl::read_xlsx(path = "arquivos/sarampo--todos_dados.xlsx",
sheet = "cobertura 1a dose",
col_names = T,
.name_repair = "unique"
) %>%
rename_with(~ str_remove(.x, "\\.0$"), matches("^\\d{4}\\.0$")) %>%
mutate(
across(
matches("^\\d{4}$"),
~ as.numeric(str_replace_all(as.character(.), ",", "."))
)
) %>%
tidyr::pivot_longer(cols = dplyr::starts_with(c("1","2")),
names_to = "Ano",
values_to = "Value") %>%
dplyr::left_join(y = BiocGenerics::unique(sigla_reg),
by = "Unidade da Federação") %>%
dplyr::mutate("Categoria" = "cov 1a dose")
cov_1a_dose# A tibble: 621 × 6
`Unidade da Federação` Ano Value Região Sigla Categoria
<chr> <chr> <dbl> <chr> <chr> <chr>
1 Acre 2003 99.4 "Região\\nNorte" AC cov 1a dose
2 Acre 2004 123. "Região\\nNorte" AC cov 1a dose
3 Acre 2005 95.1 "Região\\nNorte" AC cov 1a dose
4 Acre 2006 90.8 "Região\\nNorte" AC cov 1a dose
5 Acre 2007 109. "Região\\nNorte" AC cov 1a dose
6 Acre 2008 94.2 "Região\\nNorte" AC cov 1a dose
7 Acre 2009 103. "Região\\nNorte" AC cov 1a dose
8 Acre 2010 96.9 "Região\\nNorte" AC cov 1a dose
9 Acre 2011 105. "Região\\nNorte" AC cov 1a dose
10 Acre 2012 90.3 "Região\\nNorte" AC cov 1a dose
# ℹ 611 more rows2ª Dose - Cobertura vacinal por estado e ano
## cobertura 2a dose ----
cov_2a_dose <- readxl::read_xlsx(path = "arquivos/sarampo--todos_dados.xlsx",
sheet = "cobertura 2a dose",
col_names = T,
.name_repair = "unique"
) %>%
rename_with(~ str_remove(.x, "\\.0$"), matches("^\\d{4}\\.0$")) %>%
mutate(
across(
matches("^\\d{4}$"),
~ as.numeric(str_replace_all(as.character(.), ",", "."))
)
) %>%
tidyr::pivot_longer(cols = dplyr::starts_with(c("1","2")),
names_to = "Ano",
values_to = "Value") %>%
dplyr::left_join(y = BiocGenerics::unique(sigla_reg),
by = "Unidade da Federação") %>%
dplyr::mutate("Categoria" = "cov 2a dose")
cov_2a_dose# A tibble: 351 × 6
`Unidade da Federação` Ano Value Região Sigla Categoria
<chr> <chr> <dbl> <chr> <chr> <chr>
1 Acre 2013 27.1 "Região\\nNorte" AC cov 2a dose
2 Acre 2014 61.6 "Região\\nNorte" AC cov 2a dose
3 Acre 2015 51.7 "Região\\nNorte" AC cov 2a dose
4 Acre 2016 64.2 "Região\\nNorte" AC cov 2a dose
5 Acre 2017 57 "Região\\nNorte" AC cov 2a dose
6 Acre 2018 71.9 "Região\\nNorte" AC cov 2a dose
7 Acre 2019 78.6 "Região\\nNorte" AC cov 2a dose
8 Acre 2020 41.2 "Região\\nNorte" AC cov 2a dose
9 Acre 2021 26.0 "Região\\nNorte" AC cov 2a dose
10 Acre 2022 37.3 "Região\\nNorte" AC cov 2a dose
# ℹ 341 more rowsÓbitos por estado e ano
## óbitos ----
obitos_sarampo <- readxl::read_xlsx(path = "arquivos/sarampo--todos_dados.xlsx",
sheet = "Óbitos",
col_names = T,
.name_repair = "unique"
) %>%
rename_with(~ paste0("", sub("^(\\d{4})$", "\\1", .x)), matches("^\\d{4}$")) %>%
tibble::tibble() %>%
dplyr::mutate(
dplyr::across(
dplyr::matches("^\\d{4}$"),
~ as.numeric(stringr::str_replace_all(as.character(.), ",", "."))
)
) %>%
tidyr::pivot_longer(cols = dplyr::starts_with(c("1","2")),
names_to = "Ano",
values_to = "Value") %>%
dplyr::mutate(Ano = str_remove_all(Ano, pattern = "Ano ")) %>%
dplyr::left_join(y = sigla_reg,
by = "Unidade da Federação") %>%
dplyr::mutate("Categoria" = "Óbitos")
obitos_sarampo# A tibble: 972 × 6
`Unidade da Federação` Ano Value Região Sigla Categoria
<chr> <chr> <dbl> <chr> <chr> <chr>
1 Rondônia 1990 6 "Região\\nNorte" RO Óbitos
2 Rondônia 1991 1 "Região\\nNorte" RO Óbitos
3 Rondônia 1992 0 "Região\\nNorte" RO Óbitos
4 Rondônia 1993 0 "Região\\nNorte" RO Óbitos
5 Rondônia 1994 0 "Região\\nNorte" RO Óbitos
6 Rondônia 1995 0 "Região\\nNorte" RO Óbitos
7 Rondônia 1996 0 "Região\\nNorte" RO Óbitos
8 Rondônia 1997 0 "Região\\nNorte" RO Óbitos
9 Rondônia 1998 0 "Região\\nNorte" RO Óbitos
10 Rondônia 1999 0 "Região\\nNorte" RO Óbitos
# ℹ 962 more rowscasos por estado e ano
## casos ----
casos_sarampo <- readxl::read_xlsx(path = "arquivos/sarampo--todos_dados.xlsx",
sheet = "Casos",
col_names = T,
.name_repair = "unique"
) %>%
rename_with(~ paste0("", sub("^(\\d{4})$", "\\1", .x)), matches("^\\d{4}$")) %>%
tibble::tibble() %>%
dplyr::mutate(
dplyr::across(
dplyr::matches("^\\d{4}$"),
~ as.numeric(stringr::str_replace_all(as.character(.), ",", "."))
)
) %>%
tidyr::pivot_longer(cols = dplyr::starts_with(c("1","2")),
names_to = "Ano",
values_to = "Value") %>%
dplyr::mutate(Ano = str_remove_all(Ano, pattern = "Ano ")) %>%
dplyr::left_join(y = sigla_reg,
by = "Unidade da Federação") %>%
dplyr::mutate("Categoria" = "casos")
casos_sarampo# A tibble: 972 × 6
`Unidade da Federação` Ano Value Região Sigla Categoria
<chr> <chr> <dbl> <chr> <chr> <chr>
1 Rondônia 1990 165 "Região\\nNorte" RO casos
2 Rondônia 1991 211 "Região\\nNorte" RO casos
3 Rondônia 1992 48 "Região\\nNorte" RO casos
4 Rondônia 1993 71 "Região\\nNorte" RO casos
5 Rondônia 1994 35 "Região\\nNorte" RO casos
6 Rondônia 1995 23 "Região\\nNorte" RO casos
7 Rondônia 1996 7 "Região\\nNorte" RO casos
8 Rondônia 1997 17 "Região\\nNorte" RO casos
9 Rondônia 1998 22 "Região\\nNorte" RO casos
10 Rondônia 1999 1 "Região\\nNorte" RO casos
# ℹ 962 more rowsCombinar tabelas
# juntar tabelas e renomear categorias ----
dados_sarampo <- bind_rows(cov_1a_dose,
cov_2a_dose,
doses_sarampo,
obitos_sarampo,
casos_sarampo,
habitantes_br) %>%
dplyr::mutate(Ano = as.numeric(Ano)) %>%
dplyr::mutate(Região = stringr::str_replace_all(Região,
pattern = "\\\\n",
replacement ="\\\n"))
dados_sarampo %>%
dplyr::select(-c(Ano,Value,Categoria)) %>%
BiocGenerics::unique() %>%
print(n=100)# A tibble: 27 × 3
`Unidade da Federação` Região Sigla
<chr> <chr> <chr>
1 Acre "Região\nNorte" AC
2 Alagoas "Região\nNordeste" AL
3 Amapá "Região\nNorte" AM
4 Amazonas "Região\nNorte" AP
5 Bahia "Região\nNordeste" BA
6 Ceará "Região\nNordeste" CE
7 Distrito Federal "Região\nCentro-Oeste" DF
8 Espírito Santo "Região\nSudeste" ES
9 Goiás "Região\nCentro-Oeste" GO
10 Maranhão "Região\nNordeste" MA
11 Mato Grosso "Região\nCentro-Oeste" MT
12 Mato Grosso do Sul "Região\nCentro-Oeste" MS
13 Minas Gerais "Região\nSudeste" MG
14 Pará "Região\nNorte" PA
15 Paraíba "Região\nNordeste" PB
16 Paraná "Região\nSul" PR
17 Pernambuco "Região\nNordeste" PE
18 Piauí "Região\nNordeste" PI
19 Rio de Janeiro "Região\nSudeste" RJ
20 Rio Grande do Norte "Região\nNordeste" RN
21 Rio Grande do Sul "Região\nSul" RS
22 Rondônia "Região\nNorte" RO
23 Roraima "Região\nNorte" RR
24 Santa Catarina "Região\nSul" SC
25 São Paulo "Região\nSudeste" SP
26 Sergipe "Região\nNordeste" SE
27 Tocantins "Região\nNorte" TO dados_sarampo %>%
# dplyr::select(-c(Ano,Value,Categoria)) %>%
# BiocGenerics::unique() %>%
dplyr::filter(!is.na(Value)) %>%
dplyr::group_by(Categoria) %>%
dplyr::summarise( "Num observações" = length(Categoria),
"Ano de iníco" = min(Ano),
"Ano de fim" = max(Ano)) # A tibble: 6 × 4
Categoria `Num observações` `Ano de iníco` `Ano de fim`
<chr> <int> <dbl> <dbl>
1 Habitantes 621 2001 2025
2 casos 972 1990 2025
3 cov 1a dose 621 2003 2025
4 cov 2a dose 351 2013 2025
5 num doses 783 1994 2022
6 Óbitos 972 1990 2025Visualizar cobertura temporal
plot_dados_panorama <- dados_sarampo %>%
dplyr::mutate(Categoria = dplyr::case_when(
Categoria %in% c("cov 2a dose") ~ "Cobertura\nda2ª dose",
Categoria %in% c("cov 1a dose") ~ "Cobertura\nda1ª dose",
Categoria %in% c("num doses") ~ "Nº de\ndoses",
Categoria %in% c("Óbitos") ~ "Nº de óbitos",
Categoria %in% c("casos") ~ "Nº de casos",
Categoria %in% c("Habitantes") ~ "Nº de habitantes")) %>%
# dplyr::select(-c(Ano,Value,Categoria)) %>%
# BiocGenerics::unique() %>%
dplyr::filter(!is.na(Value)) %>%
dplyr::group_by(Categoria, `Unidade da Federação`) %>%
ggplot(aes(y = `Unidade da Federação`,
x = Categoria,
fill = Categoria)) +
geom_tile(width = 0.75,
height = 0.75) +
facet_grid(cols = vars(Ano),
rows = vars(`Região`),
scales = "free_y",
space = "free_y"
# axes = "all_x"
) +
scale_fill_manual(values = viridis::turbo(n=8)[2:7]) +
# scale_x_continuous(breaks = seq(1990,2025,1)) +
theme(strip.text.y = element_text(angle = 0),
axis.text.y = element_text(size = 10),
axis.text.x = element_blank(),
axis.title.x = element_blank(),
axis.ticks.x = element_blank()) +
# coord_equal() +
guides(fill = guide_legend(nrow = 1)) +
theme(legend.position = "bottom")
plot_dados_panorama
ggsave(plot = plot_dados_panorama,
filename = "results/plot_dados_panorama.png",
device = "png",
units = "cm",
dpi = 600,
height = 14,
width = 52,
)Transformar dado par ao formato amplo (wide)
# tabela wide ----
dados_sarampo_wide <- dados_sarampo %>%
tidyr::pivot_wider(id_cols = c("Sigla", "Unidade da Federação", "Ano", "Região"),
names_from = "Categoria",
values_from = "Value") %>%
dplyr::mutate("Casos por Milhão" = casos/Habitantes*1000000,
"Óbitos por Milhão" = `Óbitos`/Habitantes*1000000)
dados_sarampo_wide %>% write_csv(file = "arquivos/dados_sarampo_wide.csv",
col_names = T)Análises gráficas
Configurações
#Gráficos ----
## Configurações ----
#Função para visualizar decimais apenas para valores menores que 1
label_condicional <- function(x) {
sapply(x, function(val) if (val < 1) formatC(val, format = "f", digits = 2) else formatC(val, format = "f", digits = 0))
}Gráfico
Gráfico
Gráfico
Casos e óbitos após anos de vacinação
Casos após N anos de vacinação
tbl_lag_pearson_casos <- tibble::tibble()
tbl_lag_spearman_casos <- tibble::tibble()
for (LAG in 0:4) {
## cobertura vacinal 2 anos antes dos casos ----
dados_sarampo_wide_N_anos <- dados_sarampo_wide %>%
tidyr::pivot_longer(cols = c("cov 1a dose", "cov 2a dose", "num doses", "Óbitos", "casos", "Habitantes", "Casos por Milhão","Óbitos por Milhão"),
names_to = "Categoria",
values_to = "Value") %>%
dplyr::mutate(Ano = dplyr::case_when(Categoria %in% c("cov 1a dose", "cov 2a dose") ~ Ano - LAG,
TRUE ~ Ano)) %>%
tidyr::pivot_wider(id_cols = c("Sigla", "Unidade da Federação", "Ano", "Região"),
names_from = "Categoria",
values_from = "Value")
# correlação com o teste de spearman ----
message(paste0("Realizando o teste de Spearman para ", LAG, " anos"))
spearman_test <- cor.test(
log10(dados_sarampo_wide_N_anos$`Casos por Milhão`[dados_sarampo_wide_N_anos$casos > 0]),
dados_sarampo_wide_N_anos$`cov 1a dose`[dados_sarampo_wide_N_anos$casos > 0],
method = "spearman"
)
spearman_reult <- tibble(
"Categoria" = "Casos",
"Teste" = "Spearman",
"Vacinação N anos antes" = LAG,
"Valor (S)" = spearman_test$statistic,
"Significância (Valor p)" = spearman_test$p.value,
"Correlação (rho)" = spearman_test$estimate,
"Graus de liberdade" = spearman_test$parameter
)
tbl_lag_spearman_casos <- bind_rows(tbl_lag_spearman_casos,spearman_reult)
# correlação com o teste de pearson ----
message(paste0("Realizando o teste de Pearson para ", LAG, " anos"))
pearson_test <- cor.test(
log10(dados_sarampo_wide_N_anos$`Casos por Milhão`[dados_sarampo_wide_N_anos$casos > 0]),
dados_sarampo_wide_N_anos$`cov 1a dose`[dados_sarampo_wide_N_anos$casos > 0],
method = "pearson"
)
pearson_reult <- tibble(
"Categoria" = "Casos",
"Teste" = "Pearson",
"Vacinação N anos antes" = LAG,
"Valor (t)" = pearson_test$statistic,
"Significância (Valor p)" = pearson_test$p.value,
"Correlação (cor)" = pearson_test$estimate,
"Graus de liberdade" = pearson_test$parameter
)
tbl_lag_pearson_casos <- bind_rows(tbl_lag_pearson_casos,pearson_reult)
}
tbl_lag_spearman_casos# A tibble: 5 × 6
Categoria Teste `Vacinação N anos antes` `Valor (S)` Significância (Valor…¹
<chr> <chr> <int> <dbl> <dbl>
1 Casos Spearman 0 232458 0.0361
2 Casos Spearman 1 251266 0.0000129
3 Casos Spearman 2 242614 0.000000201
4 Casos Spearman 3 219382 0.000316
5 Casos Spearman 4 155750. 0.385
# ℹ abbreviated name: ¹`Significância (Valor p)`
# ℹ 1 more variable: `Correlação (rho)` <dbl>tbl_lag_pearson_casos# A tibble: 5 × 7
Categoria Teste `Vacinação N anos antes` `Valor (t)` Significância (Valor …¹
<chr> <chr> <int> <dbl> <dbl>
1 Casos Pearson 0 -1.81 0.0739
2 Casos Pearson 1 -4.34 0.0000339
3 Casos Pearson 2 -5.78 0.0000000915
4 Casos Pearson 3 -3.73 0.000321
5 Casos Pearson 4 -1.35 0.181
# ℹ abbreviated name: ¹`Significância (Valor p)`
# ℹ 2 more variables: `Correlação (cor)` <dbl>, `Graus de liberdade` <int>Óbitos após N anos de vacinação
tbl_lag_pearson_obitos <- tibble::tibble()
tbl_lag_spearman_obitos <- tibble::tibble()
for (LAG in 0:4) {
## cobertura vacinal 2 anos antes dos casos ----
dados_sarampo_wide_N_anos <- dados_sarampo_wide %>%
tidyr::pivot_longer(cols = c("cov 1a dose", "cov 2a dose", "num doses", "Óbitos", "casos", "Habitantes", "Casos por Milhão","Óbitos por Milhão"),
names_to = "Categoria",
values_to = "Value") %>%
dplyr::mutate(Ano = dplyr::case_when(Categoria %in% c("cov 1a dose", "cov 2a dose") ~ Ano - LAG,
TRUE ~ Ano)) %>%
tidyr::pivot_wider(id_cols = c("Sigla", "Unidade da Federação", "Ano", "Região"),
names_from = "Categoria",
values_from = "Value")
# correlação com o teste de spearman ----
message(paste0("Realizando o teste de Spearman para ", LAG, " anos"))
spearman_test <- cor.test(
log10(dados_sarampo_wide_N_anos$`Óbitos por Milhão`[dados_sarampo_wide_N_anos$Óbitos > 0]),
dados_sarampo_wide_N_anos$`cov 1a dose`[dados_sarampo_wide_N_anos$Óbitos > 0],
method = "spearman"
)
spearman_reult <- tibble(
"Categoria" = "Óbitos",
"Teste" = "Spearman",
"Vacinação N anos antes" = LAG,
"Valor (S)" = spearman_test$statistic,
"Graus de liberdade" = spearman_test$parameter,
"Significância (Valor p)" = spearman_test$p.value,
"Correlação (rho)" = spearman_test$estimate,
"Graus de liberdade" = spearman_test$parameter
)
tbl_lag_spearman_obitos <- bind_rows(tbl_lag_spearman_obitos,spearman_reult)
# correlação com o teste de pearson ----
message(paste0("Realizando o teste de Pearson para ", LAG, " anos"))
pearson_test <- cor.test(
log10(dados_sarampo_wide_N_anos$`Óbitos por Milhão`[dados_sarampo_wide_N_anos$Óbitos > 0]),
dados_sarampo_wide_N_anos$`cov 1a dose`[dados_sarampo_wide_N_anos$Óbitos > 0],
method = "pearson"
)
pearson_reult <- tibble(
"Categoria" = "Óbitos",
"Teste" = "Pearson",
"Vacinação N anos antes" = LAG,
"Valor (t)" = pearson_test$statistic,
"Significância (Valor p)" = pearson_test$p.value,
"Correlação (cor)" = pearson_test$estimate,
"Graus de liberdade" = pearson_test$parameter
)
tbl_lag_pearson_obitos <- bind_rows(tbl_lag_pearson_obitos,pearson_reult)
}
tbl_lag_spearman_obitos# A tibble: 5 × 6
Categoria Teste `Vacinação N anos antes` `Valor (S)` Significância (Valor…¹
<chr> <chr> <int> <dbl> <dbl>
1 Óbitos Spearman 0 608 0.240
2 Óbitos Spearman 1 502 0.727
3 Óbitos Spearman 2 634 0.165
4 Óbitos Spearman 3 656 0.116
5 Óbitos Spearman 4 736 0.0213
# ℹ abbreviated name: ¹`Significância (Valor p)`
# ℹ 1 more variable: `Correlação (rho)` <dbl>tbl_lag_pearson_obitos# A tibble: 5 × 7
Categoria Teste `Vacinação N anos antes` `Valor (t)` Significância (Valor …¹
<chr> <chr> <int> <dbl> <dbl>
1 Óbitos Pearson 0 -1.20 0.252
2 Óbitos Pearson 1 -0.943 0.364
3 Óbitos Pearson 2 -1.82 0.0936
4 Óbitos Pearson 3 -1.86 0.0875
5 Óbitos Pearson 4 -2.01 0.0677
# ℹ abbreviated name: ¹`Significância (Valor p)`
# ℹ 2 more variables: `Correlação (cor)` <dbl>, `Graus de liberdade` <int>Gráfico
tibble::tibble(
"Estatística" = c(
tbl_lag_spearman_casos$Categoria,
tbl_lag_pearson_casos$Categoria,
tbl_lag_spearman_obitos$Categoria,
tbl_lag_pearson_obitos$Categoria),
"Anos antes" = c(
tbl_lag_spearman_casos$`Vacinação N anos antes`,
tbl_lag_pearson_casos$`Vacinação N anos antes`,
tbl_lag_spearman_obitos$`Vacinação N anos antes`,
tbl_lag_pearson_obitos$`Vacinação N anos antes`),
"Teste" = c(
tbl_lag_spearman_casos$Teste,
tbl_lag_pearson_casos$Teste,
tbl_lag_spearman_obitos$Teste,
tbl_lag_pearson_obitos$Teste),
"Correlação" = c(
tbl_lag_spearman_casos$`Correlação (rho)`,
tbl_lag_pearson_casos$`Correlação (cor)`,
tbl_lag_spearman_obitos$`Correlação (rho)`,
tbl_lag_pearson_obitos$`Correlação (cor)`),
"Significância" = c(
tbl_lag_spearman_casos$`Significância (Valor p)`,
tbl_lag_pearson_casos$`Significância (Valor p)`,
tbl_lag_spearman_obitos$`Significância (Valor p)`,
tbl_lag_pearson_obitos$`Significância (Valor p)`)) %>%
ggplot2::ggplot(aes(x = `Anos antes`,
y = `Correlação`,
group = Teste,
fill = log10(`Significância`))) +
geom_bar(stat = "identity",
position = "dodge") +
facet_grid(rows = vars(`Estatística`),
cols = vars(Teste)) +
scale_fill_gradientn(colors = rev(viridis::magma(n=100)[1:90]))
# Casos após 2 anos de vacinação
dose1_x_casos_por_milhao_plot <- dados_sarampo_wide %>%
tidyr::pivot_longer(cols = c("cov 1a dose", "cov 2a dose", "num doses", "Óbitos", "casos", "Habitantes", "Casos por Milhão","Óbitos por Milhão"),
names_to = "Categoria",
values_to = "Value") %>%
dplyr::mutate(Ano = dplyr::case_when(Categoria %in% c("cov 1a dose", "cov 2a dose") ~ Ano-2,
TRUE ~ Ano)) %>%
tidyr::pivot_wider(id_cols = c("Sigla", "Unidade da Federação", "Ano", "Região"),
names_from = "Categoria",
values_from = "Value") %>%
# dplyr::filter(!`casos` %in% c(0)) %>%
dplyr::filter(!is.na(Habitantes)) %>%
ggplot(aes(x = `cov 1a dose` ,
y = `Casos por Milhão`)) +
geom_point(aes(col = `Unidade da Federação`,
shape = Região)) +
geom_smooth(method = "lm")+
stat_regline_equation(label.x=120, label.y=3.5)+
stat_cor(label.x=120, label.y=3.75) +
labs(
title = "Correlação entre casos de Sarampo e cobertura vacinal",
x = "Cobertura da 1ª dose vacinal",
y = "Casos por Milhão de habitantes",
subtitle = "Correlação realizada entre o número de casos de um dado ano e a cobertura vacinal de 2 anos anteriores"
) +
scale_colour_manual(values = viridis::turbo(n=27)) +
scale_y_continuous(trans = "log10",
breaks = c(0.1,0.25,0.5,0.75,1,2.5,5,10,25,50,100,250,500,1000,2500, 5000,10000),
labels = label_condicional) +
theme_bw() +
theme(
# legend.position = "bottom",
legend.location = "plot",
legend.direction = "vertical",
legend.text = element_text(size = 6),
legend.title = element_text(size = 8, face = "bold"))
ggsave(plot = dose1_x_casos_por_milhao_plot,
filename = "results/plot_corr_1aDose_Casos2Anos.png",
device = "png",
units = "cm",
dpi = 600,
height = 16,
width = 24,
)
Casos e cobertura por estado e ano
tile_cov_casos_plot <- ggplot() +
# geom_tile(data = dplyr::filter(dados_sarampo_wide, `Casos por Milhão` != 0),
geom_tile(data = dados_sarampo_wide %>% dplyr::mutate(`Casos por Milhão` = case_when(`Casos por Milhão` == 0 ~ NA,
TRUE ~ (`Casos por Milhão`))) %>%
filter(Ano >= 2000),
aes(x = Ano,
y = `Unidade da Federação`,
fill = `Casos por Milhão`,
group = `Unidade da Federação`),
width = 0.5,
height = 0.5) +
scale_fill_gradientn(colours = RColorBrewer::brewer.pal(name = "Reds", n = 9 )[3:9],na.value = "#ffffff") +
geom_tile(data = dplyr::filter(dados_sarampo_wide, Ano >= 2000),
# geom_tile(data = dplyr::filter(dados_sarampo_wide, `cov 1a dose` <= 90),
# geom_tile(data = dados_sarampo_wide,
aes(x = Ano,
y = `Unidade da Federação`,
col = `cov 1a dose`,
group = `Unidade da Federação`),
alpha = 0.5,
fill = NA,
linewidth = .75,
width = 0.75,
height = 0.75) +
# scale_colour_gradientn(colours = rev(viridis::turbo(n=10))[1:6],
# scale_colour_gradientn(colours = c("darkred",rev(viridis::turbo(n=8))[1:6],"lightblue","darkblue"),
scale_colour_gradientn(colours = c("darkred","red","orange","yellow","#00ff00","darkgreen","#003300"),
breaks = c(50,75,90,100,125,150),
values = scales::rescale(x = c(50,60,75,90,100,125,150),
to = c(0,1)),
na.value = "#dddddd",
name = "Cobertura\nda 1ª dose") +
geom_point(data = dplyr::filter(dados_sarampo_wide, `Óbitos por Milhão` > 0),
aes(x = Ano,
y = `Unidade da Federação`,
size = `Óbitos por Milhão`,
group = `Unidade da Federação`),
col = "#dd2222",
alpha = 0.5) +
facet_grid(rows = vars(`Região`),
axes = "all_x",
space = "free",
scales = "free",
drop = T) +
theme_bw() +
theme(strip.text.y = element_text(angle = 0),
axis.text.x = element_text(size = 6)) +
scale_x_continuous(breaks = seq(2000,2025,1),expand = c(0.01,0.01)) +
geom_vline(xintercept = 2019.5,linetype = "dashed") +
labs(
title = "Casos e óbitos por Sarampo e cobertura vacinal da 1ª dose")
ggsave(plot = tile_cov_casos_plot,
filename = "results/tile_cov_casos_plot.png",
device = "png",
units = "cm",
dpi = 600,
height = 24,
width = 32,
)