Introdução

O tema da violência urbana suscita debates na imprensa, redes sociais e casas legislativas. A par da relevância eleitoral, suscitada pela exploração da política do medo, publicações científicas sérias divulgam estudos elaborados com rigor metodológico.

O Atlas da Violência, produzido pelo Instituto de Pesquisa Econômica Aplicada (Ipea) e pelo Fórum Brasileiro de Segurança Pública (FBSP), fornece dados para compreender a evolução da violência no país, a partir de 1996.

Com base nesse estudo, iniciamos análise exploratória do conjunto de dados, que representa a série histórica da taxa de homicídios dos municípios mineiros, com população acima de 200 mil habitantes, através de técnicas gráficas.

Para tanto, utilizamos o RStudio, interface funcional e amigável para o R, linguagem para cálculo estatístico e visualização gráfica, largamente utilizada por estatísticos e cientistas de dados.

Pré-requisitos

Bibliotecas

Para adicionar funcionalidades extras, necessárias para baixar e ler páginas web, manipular dados e ler arquivos em diversos formatos, dentre outras utilidades, carregamos as seguintes bibilotecas:

library(tidyverse)
library(rvest)
library(ggthemes)
library(gridExtra)

Definição do diretório de trabalho

Definimos o diretório de trabalho, onde serão salvos os arquivos importados:

setwd("/home/gf/Scripts/Atlas")



Recursos

Taxa de homicídios por município

Baixamos manualmente o arquivo no formato csv, disponível em http://www.ipea.gov.br/atlasviolencia/dados-series/20, contendo a série histórica da taxa de violência de todos os municípios brasileiros, de 1996 a 2016.

Listagem de municípios por população

Para identificar os municípios com o perfil desejado, utilizamos uma listagem disponível no seguinte endereço: https://pt.wikipedia.org/wiki/Lista_de_munic%C3%ADpios_do_Brasil_por_popula%C3%A7%C3%A3o.

Tarefas de pré-processamento

Coleta de dados

Extração da série histórica

Após definirmos o diretório de trabalho, importamos o arquivo no formato csv, renomeamos os cabeçalhos das colunas e arredondamos a taxa:

homicidios <- read.csv2('taxa-homicidios-brasil-por-municipio-1996-01-15-2016-01-15.csv', header=TRUE, sep=';', stringsAsFactors = FALSE)
homicidios <- plyr::rename(homicidios, c("cod"="COD", "nome"="Município", "período"="Ano", "valor"="Taxa"))
homicidios$Taxa <- round(as.numeric(homicidios$Taxa))

Extração da tabela de municípios x população

Para lermos a página web, utilizamos uma função que retorna todas as informações do código html:

url <- "https://pt.wikipedia.org/wiki/Lista_de_munic%C3%ADpios_do_Brasil_por_popula%C3%A7%C3%A3o"
page <- read_html(url)

Consultamos a quantidade de tabelas disponíveis na página:

tbls <- html_nodes(page, "table")
length((tbls))
## [1] 1
head(tbls)
## {xml_nodeset (1)}
## [1] <table class="wikitable sortable" style="text-align: center;"><tbody>\n<t ...

Em seguida, extraimos o conteúdo da única tabela existente:

tab <- as.data.frame(html_table(page, fill = FALSE))

Inspecionamos o objeto e efetuamos as modificações desejadas; como por exemplo, simplificar os nomes das variáveis, para facilitar a manipulação:

str(tab) 
## 'data.frame':    5570 obs. of  5 variables:
##  $ Posição           : chr  "1º" "2º" "3º" "4º" ...
##  $ Código.IBGE       : int  3550308 3304557 5300108 2927408 2304400 3106200 1302603 4106902 2611606 5208707 ...
##  $ Município         : chr  "São Paulo" "Rio de Janeiro" "Brasília" "Salvador" ...
##  $ Unidade.federativa: chr  "São Paulo" "Rio de Janeiro" "Distrito Federal" "Bahia" ...
##  $ População         : int  12252023 6718903 3015268 2872347 2669342 2512070 2182763 1933105 1645727 1516113 ...
names(tab)[which(names(tab)=="Código.IBGE")] <- "COD"
names(tab)[which(names(tab)=="Unidade.federativa")] <- "UF"
str(tab)
## 'data.frame':    5570 obs. of  5 variables:
##  $ Posição  : chr  "1º" "2º" "3º" "4º" ...
##  $ COD      : int  3550308 3304557 5300108 2927408 2304400 3106200 1302603 4106902 2611606 5208707 ...
##  $ Município: chr  "São Paulo" "Rio de Janeiro" "Brasília" "Salvador" ...
##  $ UF       : chr  "São Paulo" "Rio de Janeiro" "Distrito Federal" "Bahia" ...
##  $ População: int  12252023 6718903 3015268 2872347 2669342 2512070 2182763 1933105 1645727 1516113 ...

Por fim, extraimos o subconjunto representado pelos municípios que atendam, simultaneamente, os critérios definidos, isto é, pertençam ao Estado de Minas Gerais e tenham população igual ou superior a 200.000 habitantes. No total, treze municípios se encaixam no perfil desejado.

pop <- tab %>%
  filter(UF == "Minas Gerais" & População >= 200000) %>%
  select(COD, Município, População)

pop
##        COD            Município População
## 1  3106200       Belo Horizonte   2512070
## 2  3170206           Uberlândia    691305
## 3  3118601             Contagem    663855
## 4  3136702         Juiz de Fora    568873
## 5  3106705                Betim    439340
## 6  3143302        Montes Claros    409341
## 7  3154606   Ribeirão das Neves    334858
## 8  3170107              Uberaba    333783
## 9  3127701 Governador Valadares    279885
## 10 3131307             Ipatinga    263410
## 11 3167202          Sete Lagoas    239639
## 12 3122306          Divinópolis    238230
## 13 3157807          Santa Luzia    219134

Conjunto de dados

Construimos o conjunto de dados a ser analisado, a partir da junção de duas tabelas: a primeira, denominada “homicidios”, contendo a série histórica da taxa de homicidios de todos os municípios brasileiros, de 1996 a 2016; a segunda, denominada “pop”, contendo a relação dos municípios mineiros com mais de 200 mil habitantes. Ao mesclar as duas tabelas, adotando como chave o código do município no IBGE, extraímos um subconjunto que é a representação tabular do recorte amostral, contendo vinte e uma observações, por município, uma para cada registro das variáveis mapeadas:

tm <- pop %>%
  select(COD, População) %>%
  left_join(homicidios, by = "COD")

Análise exploratória de dados

Ao todo temos então 273 observações de 4 variáveis:

str(tm)
## 'data.frame':    273 obs. of  5 variables:
##  $ COD      : int  3106200 3106200 3106200 3106200 3106200 3106200 3106200 3106200 3106200 3106200 ...
##  $ População: int  2512070 2512070 2512070 2512070 2512070 2512070 2512070 2512070 2512070 2512070 ...
##  $ Município: chr  "Belo Horizonte" "Belo Horizonte" "Belo Horizonte" "Belo Horizonte" ...
##  $ Ano      : int  1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 ...
##  $ Taxa     : num  16 18 22 23 29 30 34 48 52 45 ...



O quadro abaixo apresenta um resumo das estatísticas descritivas para a variável numérica que descreve a taxa de homicídios:

summary(tm$Taxa)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00   13.00   23.00   27.32   41.00   86.00



Podemos utilizar o diagrama de caule e folha apara apresentar, de forma concisa e numérica, o conjunto de dados referente às taxas de homicídio:

stem(tm$Taxa)
## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##   0 | 111111222222333344444
##   0 | 55566666777778888888999999999999
##   1 | 0000011112222233334444444
##   1 | 555566666666666677777788888888999999
##   2 | 000000011111112222233333333444444
##   2 | 5555556666788899999
##   3 | 001111233444444
##   3 | 56666666677777788899
##   4 | 00011111122223333444444
##   4 | 566778888999999
##   5 | 001112233
##   5 | 5668889
##   6 | 0023344
##   6 | 55
##   7 | 0002
##   7 | 556
##   8 | 0
##   8 | 6



Embora o diagrama de caule e folha permita visualizar a distribuição dos valores, uma boa representação da distribuição dos dados pode ser obtida com o histograma:

histtm <- ggplot(tm) + geom_histogram(aes(x=Taxa), binwidth=5, fill="gray")
histtm



Percebemos que se trata de uma distribuição positivamente assimétrica; isto é, os dados não se distribuem normalmente ao redor da média. Nesse caso, a mediana é a medida mais apropriada para representar os dados.
Como se comportam os dados agrupados por município?

histtm <- ggplot(tm) + geom_histogram(aes(x=Taxa), binwidth=15, fill="gray")
histtm <- histtm + facet_wrap(~Município, scales="free")
histtm

Taxa Média

Definida a medida de tendência central apropriada para o conjunto de dados, calculamos a taxa média da série histórica, por município, com base na mediana. Após classificar o resultado em ordem decrescente, giramos o eixo x em noventa graus, no sentido horário, para que a informação visual seja extraída de modo mais eficiente.

hpm <- aggregate(Taxa ~ Município, data = tm, FUN = median, na.rm = TRUE)

gshpm <- ggplot(hpm) + geom_bar(aes(x=reorder(Município, +Taxa), y= Taxa),
                                stat="identity",
                                fill="gray") 
  
gshpm <- gshpm + coord_flip() +
        theme(axis.text.y = element_text(size=rel (0.8))) +
        xlab("") +
        ylab("Taxa média")

gshpm

Curvas de densidade

Adicionalmente, podemos construir um gráfico de densidade, que é uma variação do histograma, para estudar a distribuição das variáveis, utilizando a escala contínua, para o conjunto de dados:

dtm <- ggplot(tm) + geom_density(aes(x=Taxa)) +
scale_x_continuous()
dtm



Observamos que as formas de distribuição dos dados variam entre si, quando visualizamos os subconjuntos, separados por município, como já demonstrado anteriormente, quando utilizamos a mesma técnica para visualização de múltiplos histogramas.

dtma <- ggplot(tm) + geom_density(aes(x=Taxa)) +
scale_x_continuous()
dtma <- dtma + facet_wrap(~Município, scales="free")
dtma



Se sobrepusermos um modelo teórico, representado pela distribuição normal, à distribuição do conjunto de dados em análise, utilizando a função de densidade para o eixo vertical, e acrescentarmos a média, ilustrada por uma linha tracejada, visualizamos mais facilmente a assimetria, consoante o recorte intencional de dados, obtido segundo o método de seleção não-probabilístico.

gg <- ggplot(tm, aes(x=tm$Taxa)) 
gg <- gg + geom_histogram(binwidth=5, colour="black", 
                          aes(y=..density.., fill=..count..))

gg <- gg + scale_fill_gradient("Count", low="#DCDCDC", high="#7C7C7C")
gg <- gg + stat_function(fun=dnorm,
                         color="red",
                         args=list(mean=mean(tm$Taxa), 
                                   sd=sd(tm$Taxa)))
gg <- gg + xlab("Taxa de homicídios/100 mil habitantes")

gg + geom_vline(xintercept = mean(tm$Taxa, na.rm = T), size = 1, colour = "#FF3721",
                linetype = "dashed")

Diagramas de caixa

Para uma visão condensada, utilizamos o diagrama de caixa, que permite identificar valores extremos, além de exibir graficamente as cinco medidas fundamentais: a mínima, o primeiro quartil, a mediana, o segundo quartil, a máxima. O intervalo interquatílico, que é uma medida da variabilidade, representa 50% dos dados. Logo, metade das observações está situada na faixa de valores entre a mínima e a mediana. Como a mediana está próxima do primeiro quartil, confirmamos mais uma vez que os dados são positivamente assimétricos. As caudas, representadas pelas linhas verticais, indicam a direção da assimetria.

boxtm <- ggplot(tm, aes(x = factor(0), y = Taxa)) + geom_boxplot() + xlab("") +
scale_x_discrete(breaks = NULL)

out <- ggplot_build(boxtm)[["data"]][[1]][["outliers"]]
names(out) <- levels(factor(tm$x))
tidyout <- purrr::map_df(out, tibble::as_tibble, .id = "Taxa")

boxtm <- boxtm + geom_text(data = tidyout, aes(Taxa, value, label = value), 
              hjust = -.3)

boxtm



Embora seja unidimensional, isto é, dependa de apenas uma variável para ser construído, variáveis categóricas podem ser utilizadas como critério de agrupamento. Desse modo, desenhamos diagramas de caixa, separadamente, para todos os níveis da variável.

boxtm <- ggplot(tm, aes(Município, Taxa)) +
geom_boxplot() + xlab("") 

boxtm <- boxtm + theme(axis.text.x = element_text(angle = 45, hjust = 1))
boxtm

Valores extremos - Juiz de Fora

Apesar da baixa variabilidade, indicada pela forma achatada do diagrama, e valor da mediana situado na parte inferior do eixo y, Juiz de Fora contém a maior quantidade de outiliers. Para investigar esta anomalia, primeiro selecionamos o subconjunto de dados para o município; depois, calculamos as estatísticas básicas:

tmjf <- tm %>%
  filter(Município == "Juiz de Fora") %>%
  select(Ano, Taxa)

summary(tmjf$Taxa)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    6.00    8.00    9.00   12.14   12.00   27.00



Em seguida, visualizamos o diagrama de caixa isoladamente:

boxtmjf <- ggplot(tmjf, aes(x = factor(0), y = Taxa)) + geom_boxplot() + xlab("") +
scale_x_discrete(breaks = NULL)

out <- ggplot_build(boxtmjf)[["data"]][[1]][["outliers"]]
names(out) <- levels(factor(tmjf$x))
tidyout <- purrr::map_df(out, tibble::as_tibble, .id = "Taxa")

boxtmjf <- boxtmjf + geom_text(data = tidyout, aes(Taxa, value, label = value), 
              hjust = -.3)

boxtmjf

Observem a ausência da cauda superior: 75% dos dados estão situados entre a mínima e o primeiro quartil. Do mesmo modo, a distribuição dos dados é positivamente assimétrica. No entanto, cabe investigar por que um quarto dos dados são considerados outliers. Ou seja, os dados situados entre o terceiro quartil e a máxima se distanciam muito da maioria, ou por erro de mensuração ou por ocorrência de valores atípicos.

Dados discrepantes

Os dados do IPEA estão defasados dois anos em relação ao ano da publicação. Portanto, após a importação do arquivo de dados da série, do Atlas 2017, adicionei à tabela o valor obtido para o ano de 2016, para o município de Juiz de Fora, calculado com base em informações divulgadas na imprensa: http://g1.globo.com/mg/zona-da-mata/noticia/2017/01/delegado-analisa-numero-de-homicidios-em-juiz-de-fora.html.

O Atlas 2018 incluiu os indicadores de 2016; no entanto, foram identificadas divergências entre as duas bases de dados, sem justificativa aparente. Questionado, o IPEA não se manifestou a respeito.

Atlas 2017

Nesse ponto, vamos direto aos dados.

h <- read.csv2('taxa-homicidios-brasil-por-municipio-1996-2015.csv', header=FALSE, sep=';', stringsAsFactors = FALSE)
h <- plyr::rename(h, c("V1"="Município", "V2"="Ano", "V3"="Taxa"))

hjf <- h %>%
  filter(Município == "Juiz de Fora") %>%
  select(Ano, Taxa)

hjf <- hjf %>% add_row(Ano = "2016", Taxa = "25.194")

hjf$Taxa <- round(as.numeric(hjf$Taxa))

hjf
##     Ano Taxa
## 1  1996   16
## 2  1997   13
## 3  1998   21
## 4  1999   17
## 5  2000   18
## 6  2001   14
## 7  2002   16
## 8  2003   17
## 9  2004   21
## 10 2005   11
## 11 2006   12
## 12 2007   16
## 13 2008   16
## 14 2009   16
## 15 2010   20
## 16 2011   24
## 17 2012   28
## 18 2013   34
## 19 2014   46
## 20 2015   22
## 21 2016   25
summary(hjf$Taxa)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   11.00   16.00   17.00   20.14   22.00   46.00



Colocamos os dois diagramas de caixa lado a lado, para facilitar a comparação.

boxhjf <- ggplot(hjf, aes(x = factor(0), y = Taxa)) + geom_boxplot() + xlab("") +
scale_x_discrete(breaks = NULL)

out <- ggplot_build(boxhjf)[["data"]][[1]][["outliers"]]
names(out) <- levels(factor(hjf$x))
tidyout <- purrr::map_df(out, tibble::as_tibble, .id = "Taxa")

boxhjf <- boxhjf + geom_text(data = tidyout, aes(Taxa, value, label = value), 
              hjust = -.3)


boxhjf <- boxhjf + xlab("Atlas 2017")
boxtmjf <- boxtmjf + xlab("Atlas 2018")
grid.arrange(boxhjf, boxtmjf, ncol=2)

Referências:
Kabacoff, R. (2011). R in action: Data analysis and graphics with R.
Zumel, Nina and Mount, John. (2014). Practical Data Science with R.