Já deve ter passado pelo seu radar os relatórios de mobilidade que a Google liberou recentemente. O site que hospeda o conjunto de relatórios de mobilidade, os chamados COVID-19 Community Mobility Reports, começa com a seguinte frase: “Veja como sua comunidade está se deslocando de forma diferente devido ao COVID19“.
Do que se trata esses relatórios de mobilidade?
A Google oferece um recurso de Location History nos seus serviços/sistemas que monitora a localização, e consequentemente o deslocamento, dos usuários. Esses dados podem ser acessados e o serviço desativado a qualquer momento pelos usuários. De acordo com eles, esse recurso precisa ser ativado voluntariamente, pois vem desativado por padrão. Com base nessas informações, eles observaram como e para onde esses indivíduos costumavam se deslocar em um período anterior ao surto do COVID19 e como e para onde estão se deslocando agora, durante o surto. Existe um viés claro aqui. Pessoas que não possuem celular ou tablet, ou que não ativaram esse recurso, estão fora da amostra deles e isso pode ter impacto nas conclusões do relatório. Ainda assim, vale a pena dar uma olhada.
Os relatórios são divididos por países e obviamente o Brasil chamou minha atenção (você pode acessar o relatório do Brasil clicando aqui). Eles não disponibilizam os dados em um formato que facilite a sua manipulação e análise, infelizmente. No entanto, O Vitor Baptista fez o favor de resolver esse impedimento para nós 🙂 Utilizando o Python, ele conseguiu extrair os dados do relatório (tanto brasileiro como de alguns outros países) e gerar um arquivo CSV (comma-separated value) que você pode baixar clicando aqui.
Hora de começar a análise. Pronto?
A primeira coisa que me veio em mente foi ver como essa variação nos deslocamentos tem relação com variáveis desses estados (densidade populacional, por exemplo) e na taxa de crescimento do número de casos de COVID19. Com isso em mente, baixei os dados mais atuais do dashboard oficial do Ministério da Saúde, também um CSV, assim como dados sobre os estados no site do IBGE. Coloquei todos esses arquivos em um diretório, entrei nele pelo R, carreguei algumas bibliotecas que irei precisar já agora no início da análise e li os arquivos.
setwd('/home/mribeirodantas/COVID mobility/') | |
library(readr) | |
library(dplyr) | |
library(lubridate) | |
# Source: https://covid.saude.gov.br/ | |
brazil <- read_csv2(file = 'brazil n_cases.csv') | |
# Source: https://github.com/vitorbaptista/google-covid19-mobility-reports | |
mobility <- read_delim(file = 'mobility_reports.csv', delim = ',') | |
# Source: https://www.ibge.gov.br/cidades-e-estados | |
ibge <- read_delim(file = 'IBGE.csv', ',', skip = 1) |
Dataset do Ministério da Saúde
O dataset que baixei do dashboard do Ministério da Saúde tem 7 variáveis: regiao, estado, data, casosNovos, casosAcumulados, obitosNovos e obitosAcumulados. A idéia é trabalhar estado a estado, ou seja, agrupar por estado, ordenar pela data, e calcular a taxa de crescimento de um dia para o outro. O código do que acabei de descrever segue abaixo:
# Função para calcular a taxa de crescimento diária | |
growth_rate <- function(x) { | |
(x/lag(x)-1)*100 | |
} | |
# Transformando o campo data para uma data válida | |
brazil$data <- dmy(brazil$data) | |
# Calculando a taxa de cresc diária do número de casos | |
brazil %>% | |
select(estado, data, casosAcumulados) %>% | |
arrange(estado, data) %>% | |
group_by(estado) %>% | |
mutate(growth_rate_cases = growth_rate(casosAcumulados)) -> df | |
# Calculando a taxa de cresc diária do número de óbitos | |
brazil %>% | |
select(estado, data, obitosAcumulados) %>% | |
arrange(estado, data) %>% | |
group_by(estado) %>% | |
mutate(growth_rate_obitos = growth_rate(obitosAcumulados)) -> df2 | |
# Juntando isso em um só dataframe | |
df <- merge(df, df2, by=c('estado', 'data')) |
Um outro ponto que precisamos tomar cuidado é o de comparar as mesmas semanas epidemiológicas. Como o número de casos costuma dobrar com uma frequência alarmante (dias, quando não poucos dias), comparar o cenário de um estado após um mês com um estado que acabou de identificar seus primeiros casos não seria adequado. O correto seria comparar o começo com o começo, e assim por diante. Com base nisso, nós iremos ignorar as linhas, por estado, anteriores ao primeiro caso confirmado.
df %>% | |
group_by(estado) %>% | |
slice(which(growth_rate_casos == Inf)[1]:n()) -> df |
Talvez você fique curioso com o Inf no código acima. Se você calcular a taxa de crescimento de 0 para qualquer outro número, em uma parte da sua equação você terá um número sobre zero e isso é igual a infinito. Logo, o primeiro momento que encontramos um Inf é porque ocorreu o primeiro caso. E removemos o que estava antes.
Agora que temos o nosso dataframe com as informações de taxa de crescimento, iremos calcular a média das taxas diárias (de novos casos e de novos óbito levando em consideração os valores acumulados) para cada estado e estamos prontos com relação a esse dataset. Após isso, apagamos as variáveis que não iremos mais utilizar. O filtro de is.finite é justamente para remover os Inf que sobraram.
df %>% | |
group_by(estado) %>% | |
filter(is.finite(growth_rate_casos)) %>% | |
filter(is.finite(growth_rate_obitos)) %>% | |
summarise_at(c('growth_rate_casos', 'growth_rate_obitos'), funs(mean)) -> casos_obitos | |
rm(df2, brasil, df, growth_rate) |
Dataset de Mobilidade do Google
Esse dataset contém as seguintes variáveis: country, region, updated_at, retail_and_recreation, grocery_and_pharmacy, parks, transit_stations, workplaces e residential. Os dados que eu baixei foram atualizados pela última vez em 29 de Março. Se você abriu o relatório, verá a definição dessas variáveis, que trago abaixo:
retail_and_recreation (Varejo e recreação): Refere-se ao deslocamento para lugares como restaurantes, cafés, shopping centers, parques temáticos, museus, bibliotecas e cinemas.
grocery_and_pharmacy (Mercearia e farmácias): Refere-se ao deslocamento para lugares como mercearia, mercados, armazéns de alimentos, lojas de alimentos especializados, drogarias e farmácias.
parks (Parques): Refere-se ao deslocamento para lugares como parques nacionais, praias públicas, marinas, parques para cães, praças e jardins públicos.
transit_stations (Estações de transporte): Refere-se ao deslocamento para lugares como centros de transporte público, como estações de metrô, ônibus e trem.
workplaces (Locais de trabalho): Refere-se ao deslocamento para o seu local de trabalho.
residential (Área residencial): Refere-se ao deslocamento para residências.
A primeira coisa que iremos fazer é remover a coluna de updated_at, que não nos importa, filtrar apenas para os dados do Brasil e corrigir o nome dos estados para o formato do arquivo do Ministério da Saúde, que são as siglas das unidades da federação (RN, SP, RJ, etc). Isso será importante na hora de mesclarmos os dados.
mobility %>% | |
filter(country == 'Brazil') %>% | |
select(-c('updated_at')) %>% | |
# Updating state names | |
mutate(region = case_when(region == 'State of Acre' ~ 'AC', | |
region == 'State of Amapá' ~'AP', | |
region == 'State of Espírito Santo' ~ 'ES', | |
region == 'State of Maranhão' ~ 'MA', | |
region == 'State of Mato Grosso do Sul' ~ 'MS', | |
region == 'State of Paraná' ~ 'PR', | |
region == 'State of Pará' ~ 'PA', | |
region == 'State of Piauí' ~ 'PI', | |
region == 'State of Rio Grande do Sul' ~ 'RS', | |
region == 'State of Rondônia' ~ 'RO', | |
region == 'State of Santa Catarina' ~ 'SC', | |
region == 'State of São Paulo' ~ 'SP', | |
region == 'Federal District' ~ 'DF', | |
region == 'State of Alagoas' ~ 'AL', | |
region == 'State of Amazonas' ~ 'AM', | |
region == 'State of Ceará' ~ 'CE', | |
region == 'State of Goiás' ~ 'GO', | |
region == 'State of Mato Grosso' ~ 'MT', | |
region == 'State of Minas Gerais' ~ 'MG', | |
region == 'State of Paraíba' ~ 'PB', | |
region == 'State of Pernambuco' ~ 'PE', | |
region == 'State of Rio Grande do Norte' ~ 'RN', | |
region == 'State of Rio de Janeiro' ~ 'RJ', | |
region == 'State of Roraima' ~ 'RR', | |
region == 'State of Sergipe' ~ 'SE', | |
region == 'State of Bahia' ~ 'BA', | |
region == 'State of Tocantins' ~ 'TO')) -> mobility |
Dataset do IBGE
O dataset do IBGE também precisa de um leve pré-processamento. Iremos alterar também o nome dos estados, ignorar linhas que não tem relação com a informação que queremos e colunas também. Esse dataset tem as seguintes variáveis: UF, Código, Gentílico, Governador [2019], Capital [2010], Área Territorial – km² [2018], População estimada – pessoas [2019], Densidade demográfica – hab/km² [2010], Matrículas no ensino fundamental – matrículas [2018], IDH [2010], Receitas realizadas… e várias outras colunas que não me interessam para essa análise. Mantive apenas população, densidade demográfica e nome do estado.
ibge %>% | |
select(`UF [-]`, `População estimada – pessoas [2019]`, | |
`Densidade demográfica – hab/km² [2010]`) %>% | |
mutate(`UF [-]` = case_when(`UF [-]` == 'Acre' ~ 'AC', | |
`UF [-]` == 'Alagoas' ~ 'AL', | |
`UF [-]` == 'Amapá' ~ 'AP', | |
`UF [-]` == 'Amazonas' ~ 'AM', | |
`UF [-]` == 'Bahia' ~ 'BA', | |
`UF [-]` == 'Ceará' ~ 'CE', | |
`UF [-]` == 'Distrito Federal' ~ 'DF', | |
`UF [-]` == 'Espírito Santo' ~ 'ES', | |
`UF [-]` == 'Goiás' ~ 'GO', | |
`UF [-]` == 'Maranhão' ~ 'MA', | |
`UF [-]` == 'Mato Grosso' ~ 'MT', | |
`UF [-]` == 'Mato Grosso do Sul' ~ 'MS', | |
`UF [-]` == 'Minas Gerais' ~ 'MG', | |
`UF [-]` == 'Pará' ~ 'PA', | |
`UF [-]` == 'Paraíba' ~ 'PB', | |
`UF [-]` == 'Paraná' ~ 'PR', | |
`UF [-]` == 'Pernambuco' ~ 'PE', | |
`UF [-]` == 'Piauí' ~ 'PI', | |
`UF [-]` == 'Rio de Janeiro' ~ 'RJ', | |
`UF [-]` == 'Rio Grande do Norte' ~ 'RN', | |
`UF [-]` == 'Rio Grande do Sul' ~ 'RS', | |
`UF [-]` == 'Rondônia' ~ 'RO', | |
`UF [-]` == 'Roraima' ~ 'RR', | |
`UF [-]` == 'Santa Catarina' ~ 'SC', | |
`UF [-]` == 'São Paulo' ~ 'SP', | |
`UF [-]` == 'Sergipe' ~ 'SE', | |
`UF [-]` == 'Tocantins' ~ 'TO')) -> ibge | |
colnames(ibge) <- c('estado', 'pop', 'dens_demo') | |
# Ignorar a última linha que não tem informação relacionada aos estados | |
ibge <- slice(ibge, 1:27) |
Juntei os três datasets em apenas um, com todas as informações que preciso, corrigi o tipo dos dados em cada coluna e apaguei as variáveis que não serão mais utilizadas.
df <- merge(x = casos_obitos, y = mobility, by.x = c('estado'), by.y = c('region')) | |
df <- merge(x = df, y = ibge, by.x = c('estado')) | |
rm(casos_obitos, mobility, ibge) | |
# Setting the correct column types | |
df[, c(2,3,5:12)] <- sapply(df[, c(2,3,5:12)], as.numeric) |
Visualização dos dados
Com base nesses dados, algo que eu gostaria de ver é a taxa média de crescimento diário do número de casos e de óbito com relação às mudanças nos deslocamentos que os indivíduos daquele estado tiveram. Por exemplo, o código abaixo irá plotar o gráfico abaixo dele. O tamanho dos pontos tem relação com a população dos estados e a cor com a densidade demográfica daquele estado.
ggplot(df, aes(x=growth_rate_casos, y=retail_and_recreation,size=pop,color=dens_demo)) + | |
geom_point(shape=19) + | |
geom_smooth() + | |
scale_color_gradient(name = "Densidade\ndemográfica", low="blue", high="red", | |
labels = c('100 hab/km²', | |
'200 hab/km²', | |
'300 hab/km²', | |
'400 hab/km²', | |
'500 hab/km²')) + | |
scale_size_continuous(name = "População", labels = c('1 milhão hab', | |
'10 milhões hab', | |
'20 milhões hab', | |
'30 milhões hab', | |
'40 milhões hab', | |
'50 milhões hab')) + | |
xlab('Taxa de crescimento médio de número de casos por dia desde o primeiro caso confirmado') + | |
ylab('Variação de mobilidade para fins de\n lazer e varejo') + | |
labs(caption = "Ribeiro-Dantas, MC. April, 2020. mribeirodantas.xyz/blog") + | |
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + | |
scale_x_continuous(labels = function(x) paste0(sprintf("%0.0f", x), '%')) + | |
ggtitle("Taxa média diária de casos novos e deslocamento") + | |
ggrepel::geom_label_repel(aes(label = df$estado), size=3) |
A maior variação ocorreu nesse tipo de mobilidade, para locais de varejo e lazer, o que é positivo do ponto de vista do lockdown (lazer), mas polêmico do ponto de vista econômico (varejo), ainda que o e-commerce esteja com tudo em alguns segmentos e países. O estado com a menor redução foi o Pará e ainda assim com 60% de redução, o que é uma redução substancial.
O código abaixo irá plotar o gráfico abaixo dele, dessa vez com relação a taxa média diária de número de óbitos novos.
ggplot(df, aes(x=growth_rate_obitos, y=workplaces,size=pop,color=dens_demo)) + | |
geom_point(shape=19) + | |
geom_smooth() + | |
scale_color_gradient(name = "Densidade\ndemográfica", low="blue", high="red", | |
labels = c('100 hab/km²', | |
'200 hab/km²', | |
'300 hab/km²', | |
'400 hab/km²', | |
'500 hab/km²')) + | |
scale_size_continuous(name = "População", labels = c('1 milhão hab', | |
'10 milhões hab', | |
'20 milhões hab', | |
'30 milhões hab', | |
'40 milhões hab', | |
'50 milhões hab')) + | |
xlab('Taxa de crescimento médio de número de óbitos por dia desde o primeiro caso confirmado') + | |
ylab('Variação de mobilidade para ambiente\n de trabalho') + | |
labs(caption = "Ribeiro-Dantas, MC. April, 2020. mribeirodantas.xyz/blog") + | |
ggtitle("Taxa média de óbitos novos e deslocamento") + | |
ggrepel::geom_label_repel(aes(label = df$estado), size=3) |
Eu pensei em não publicar a imagem acima, por receio de que alguém a utilizasse com objetivos de desinformação. No entanto, preferi deixá-la e explicar aqui por que ela passa uma informação equivocada e porque eu decidi não plotar as outras imagens com base na taxa média diária d óbitos novos. O Espírito Santo teve o primeiro óbito de COVID19 dia 2 de Abril, no dia 3 de Abril teve mais três óbitos (300% de crescimento) e mais 1 no dia 4 de Abril (25% de crescimento). Na média dessas duas taxas de crescimento, nós temos mais de ~162%, o que joga o Espírito Santo para o canto direito do gráfico acima, bastante diferente dos outros estados. Ele merece isso? O Ceará teve 3 óbitos em 26 de Março, e mais um óbito em 28 Março (33% de aumento). Média? 16.6%.
Alterando o nome da variável no código do primeiro gráfico, você poderá obter os outros gráficos. Seguem eles abaixo.
Esse é o único gráfico onde tivemos um aumento de mobilidade, em vez de diminuição, mas faz sentido. As pessoas estão ficando em casa, ou se deslocando no entorno para praticar exercícios físicos, como caminhada.
Como experimento, decidi utilizar o número total de casos confirmados em vez da taxa média diária de novos casos. Segue código e gráficos abaixo.
brasil <- read_csv2(file = 'brasil n_casos.csv') | |
brasil %>% | |
group_by(estado) %>% | |
summarise(total = sum(casosNovos)) -> brasil | |
df <- merge(x = df, y = brasil, by = c('estado')) | |
library(ggplot2) | |
# Scatter plots | |
# retail_and_recreation | |
ggplot(df, aes(x=growth_rate_obitos, y=retail_and_recreation,size=pop,color=dens_demo)) + | |
geom_point(shape=19) + | |
geom_smooth() + | |
scale_color_gradient(name = "Densidade\ndemográfica", low="blue", high="red", | |
labels = c('100 hab/km²', | |
'200 hab/km²', | |
'300 hab/km²', | |
'400 hab/km²', | |
'500 hab/km²')) + | |
scale_size_continuous(name = "População", labels = c('1 milhão hab', | |
'10 milhões hab', | |
'20 milhões hab', | |
'30 milhões hab', | |
'40 milhões hab', | |
'50 milhões hab')) + | |
xlab('Número total de casos confirmados em 5 de Abril') + | |
ylab('Variação de mobilidade para fins de\n lazer e varejo') + | |
labs(caption = "Ribeiro-Dantas, MC. April, 2020. mribeirodantas.xyz/blog") + | |
ggtitle("Total de casos confirmados e deslocamento") + | |
ggrepel::geom_label_repel(aes(label = df$estado), size=3) |
Como você pode ver, há uma disparidade muito grande no eixo x com São Paulo deslocado quando comparado aos outros estados. Por essa razão, decidi aplicar uma transformação logarítmica (log10) no eixo x, como você pode ver no código abaixo e gráficos a seguir.
ggplot(df, aes(x=total, y=retail_and_recreation,size=pop,color=dens_demo)) + | |
geom_point(shape=19) + | |
geom_smooth() + | |
scale_color_gradient(name = "Densidade\ndemográfica", low="blue", high="red", | |
labels = c('100 hab/km²', | |
'200 hab/km²', | |
'300 hab/km²', | |
'400 hab/km²', | |
'500 hab/km²')) + | |
scale_size_continuous(name = "População", labels = c('1 milhão hab', | |
'10 milhões hab', | |
'20 milhões hab', | |
'30 milhões hab', | |
'40 milhões hab', | |
'50 milhões hab')) + | |
xlab('Número total de casos confirmados em 5 de Abril') + | |
ylab('Variação de mobilidade para fins de\n lazer e varejo') + | |
labs(caption = "Ribeiro-Dantas, MC. April, 2020. mribeirodantas.xyz/blog") + | |
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + | |
scale_x_continuous(labels = function(x) paste0(sprintf("%0.0f", x), '%')) + | |
scale_x_continuous(trans='log10') + | |
ggtitle("Total de casos confirmados e deslocamento") + | |
ggrepel::geom_label_repel(aes(label = df$estado), size=3) |
Eu poderia seguir plotando os dados dos 6 tipos de deslocamento com a taxa de letalidade, dentre outras métricas, mas acredito que é melhor nos encaminharmos para os “finalmentes” desse post.
Ressalvas
Há uma vontade incontrolável dentro de todos nós de querer ser útil, principalmente em um momento como o atual, de pandemia. Nos sentimos afoitos, ansiosos para pular na primeira oportunidade que podemos para tirar conclusões de informações que nos são dadas. Embora a boa intenção seja louvável, precisamos ter muita responsabilidade. Algumas pessoas podem se sentir tentadas a avaliar a intervenção de quarentena com base nesses dados. Isso é uma péssima ideia. Os dados do Google comparam dois momentos no tempo, antes da intervenção (um ponto no tempo que na opinião de alguns pode ser visto como uma escolha arbitrária) e o dia 29 de Março de 2020. Ainda que tivéssemos essa informação ao longo de todo o período, por dia, ainda estaríamos estabelecendo apenas uma correlação. E, como todos nós sabemos, correlação não implica em causalidade.
De mesmo modo, algumas de minhas decisões podem parecer arbitrárias para terceiros. Vamos olhar para a métrica que eu utilizei, por exemplo, a taxa média diária de crescimento do número de casos ou de óbitos. Digamos que exista um estado A que ao longo dos 10 primeiros dias após o primeiro caso confirmado teve seus casos distribuídos acumuladamente da seguinte maneira: 1, 1, 1, 1, 1, 1, 1, 1, 8, 10, com um total de 10 pacientes confirmados. Isto é, 1 paciente confirmado no primeiro dia e de repente mais sete pacientes confirmados no nono dia e mais dois no décimo dia. Um outro estado B teve ao longo dos 10 primeiros dias após o primeiro caso confirmado os casos distribuídos da seguinte maneira: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10. Também 10 pacientes. Vamos calcular a taxa?
> A <- c(1,1,1,1,1,1,1,1,8,10) | |
> B <- c(1,2,3,4,5,6,7,8,9,10) | |
> growth_rate(A) | |
[1] NA 0 0 0 0 0 0 0 700 25 | |
> sum(growth_rate(A), na.rm=TRUE)/length(A) | |
[1] 72.5 | |
> | |
> growth_rate(B) | |
[1] NA 100.00000 50.00000 33.33333 25.00000 20.00000 16.66667 14.28571 12.50000 | |
[10] 11.11111 | |
> sum(growth_rate(B), na.rm=TRUE)/length(B) | |
[1] 28.28968 |
Como você pode ver, o estado A teve uma taxa diária média de novos casos mais de duas vezes maior que o estado B. Se esses dados refletem exatamente a realidade, eu ficaria um pouco mais tranquilo, mas infelizmente não temos como ter certeza. Talvez, o estado A fosse o estado B em um universo paralelo onde o laboratório que realiza os testes recebeu muitas amostras e acabou liberando o resultado de 7 pacientes de uma vez no nono dia. Ainda que eu imagine que os estados atualizem seus dados para confirmar o caso no dia que o paciente teve a amostra coletada, não há garantias de que isso esteja sendo feito dessa maneira, e por todos os estados.
Em um levantamento da Open Knowledge Brasil, que avalia conteúdo, formato e detalhamento do que é divulgado em portais de governos estaduais e federal, os autores concluíram que a transparência com relação ao COVID-19 é insuficiente em 90% dos estados brasileiros. Eu confesso que estou até impressionado com algumas ações do governo federal, dos governos estaduais e até da iniciativa privada, mas para fazermos análises responsáveis e profissionais, precisamos de mais.
Estou publicando esse post de modo a divulgar esses datasets, ensinar um pouco de R e análise de dados e estimular outros que por ventura possam ter mais/melhores dados e possam realizar uma análise mais profissional. Como dever de casa, para quem está atrás de um 😛 seria interessante ir atrás das publicações nos diários oficiais de quando cada estado instituiu o lockdown (proibindo o varejo, por exemplo) e fazer uma análise das variações com base nisso 😉
Em breve vou organizar os códigos e disponibilizar em um repositório do GitHub.
PS: O Rafael Gomes escreveu um código (e um post utilizando o código) para extrair mais do que o valor do dia 29 de Março, mas ao longo de todo o período que aparece nos gráficos do relatório. Serve de empurrãozinho para quem quer fazer algo mais elaborado.