#----------------------------------------------#
# MAPA MORTALIDADE FEMININA EM PE (2009 -2013) #
#----------------------------------------------#


# instalar e carregar pacotes necessarios para a realizacao das analises
install.packages("maptools", "plotly", "rgdal", "readxl", "dplyr", "stringr")

library(maptools)
library(plotly)
library(rgdal)
library(readxl)
library(dplyr)
library(stringr)

# carregar arquivo shapefile de mapa de pernambuco
geo_pernambuco <- readShapeSpatial("pernambuco.shp")

# carregar banco de dados de violencia contra a mulher no brasil 
violencia_mulher <- read_excel("MapaViolencia2015_mulheres (1).xlsx")

# ordenar por nome da cidade em PE
violencia_mulher_pe <- subset(violencia_mulher, UF == "PE")
violencia_mulher_pe <- with(violencia_mulher_pe,  violencia_mulher_pe[order(municipio) , ])

geo_pernambuco <- with(geo_pernambuco,  geo_pernambuco[order(geo_pernambuco$NOME), ])

# juntar bases
mapaviolencia_mulher_pe <- data.frame(geo_pernambuco, violencia_mulher_pe)

# fixar banco mapa da mapaviolencia_mulher_pe na area de trabalho do R
attach(mapaviolencia_mulher_pe)

# calcular taxa de homicidio feminino em PE entre 2009 e 2013 por 10000 hab
mapaviolencia_mulher_pe$taxa_homicidio_pe <- ((hominicio_mulher_2009 + hominicio_mulher_2010 +
                                                 hominicio_mulher_2011 + hominicio_mulher_2012 + 
                                                 hominicio_mulher_2013) / mapaviolencia_mulher_pe$pop_media) * 10000



#---GRAFICO COMPARATIVO DA TAXA DE HOMICIDIO DE MULHERES EM PE e BRASIL (2009 - 2013)---#

# retirar missing cases (casos faltantes) do banco
violencia_mulher2 <- violencia_mulher[!is.na(violencia_mulher$hominicio_mulher_2009),]
violencia_mulher2 <- violencia_mulher[!is.na(violencia_mulher$hominicio_mulher_2010),]
violencia_mulher2 <- violencia_mulher[!is.na(violencia_mulher$hominicio_mulher_2011),]
violencia_mulher2 <- violencia_mulher[!is.na(violencia_mulher$hominicio_mulher_2012),]
violencia_mulher2 <- violencia_mulher[!is.na(violencia_mulher$hominicio_mulher_2013),]

# criar variavel com medias de homicidio do brasil e pernambuco
homicidio_mulher_pe_bra <- c( 
  
  # calcular media brasileira em cada ano
  sum(violencia_mulher2$hominicio_mulher_2009) / length(violencia_mulher2$hominicio_mulher_2009),
  sum(violencia_mulher2$hominicio_mulher_2010) / length(violencia_mulher2$hominicio_mulher_2010),
  sum(violencia_mulher2$hominicio_mulher_2011) / length(violencia_mulher2$hominicio_mulher_2011),
  sum(violencia_mulher2$hominicio_mulher_2012) / length(violencia_mulher2$hominicio_mulher_2012),
  sum(violencia_mulher2$hominicio_mulher_2013) / length(violencia_mulher2$hominicio_mulher_2013),
  
  # calcular media pernambucana em cada ano
  sum(violencia_mulher_pe$hominicio_mulher_2009) / length(violencia_mulher_pe$hominicio_mulher_2009),
  sum(violencia_mulher_pe$hominicio_mulher_2010) / length(violencia_mulher_pe$hominicio_mulher_2010),
  sum(violencia_mulher_pe$hominicio_mulher_2011) / length(violencia_mulher_pe$hominicio_mulher_2011),
  sum(violencia_mulher_pe$hominicio_mulher_2012) / length(violencia_mulher_pe$hominicio_mulher_2012),
  sum(violencia_mulher_pe$hominicio_mulher_2013) / length(violencia_mulher_pe$hominicio_mulher_2013)
)

# criar variavel com o nome do local de cada caso
local <- c("Brasil","Brasil","Brasil","Brasil","Brasil",
           "Pernambuco","Pernambuco","Pernambuco","Pernambuco","Pernambuco")

# criar variavel com os anos de cada caso
anos2 <- c(2009, 2010, 2011, 2012, 2013, 
           2009, 2010, 2011, 2012, 2013)

# arredondar com 3 casas decimais os dados de homicidios
homicidio_mulher_pe_bra <- round(homicidio_mulher_pe_bra, 3)

# criar banco de dados com dados da media brasileira e pernambucana (2009 - 2013)
violencia_mulher_pe_bra_anos <- data.frame(homicidio_mulher_pe_bra, local, anos2)


## GRAFICO ##

# criar grafico 
plot2 <- plot_ly(data = violencia_mulher_pe_bra_anos,
                 x = anos2,
                 y = homicidio_mulher_pe_bra,
                 type = "scatter",
                 mode = "line",
                 color = local,
                 colors = "Set1")
plot2

# lista criada para inserir valores dos homicidios no grafico
a <- list()
for (i in seq_len(nrow(violencia_mulher_pe_bra_anos))) {
  m <- violencia_mulher_pe_bra_anos[i, ]
  a[[i]] <- list(
    x = m$anos2, 
    y = m$homicidio_mulher_pe_bra,
    text = c(m$homicidio_mulher_pe_bra),
    xref = "x", 
    yref = "y",
    showarrow = TRUE, 
    arrowhead = 0,
    ax = 5, 
    ay = -30)
}

# editar layout do grafico
layout(plot2,   annotations = a,
       xaxis = list(title = ""),
       yaxis = list(title = "Taxa de Homicídio de Mulheres p/ 10.000 hab"),
       title = "")

#---- calcular razao de chance de morte feminina em PE em relacao ao Brasil ----#

# visualizar banco de dados
violencia_mulher_pe_bra_anos

# criar objetos da media de mortalidade feminina em Pernambuco e Brasil
x1 <- subset(violencia_mulher_pe_bra_anos, local == "Brasil")
x2 <- subset(violencia_mulher_pe_bra_anos, local == "Pernambuco")

# calculo da razao
(sum(x2$homicidio_mulher_pe_bra) / length(x2$homicidio_mulher_pe_bra) ) /
  (sum(x1$homicidio_mulher_pe_bra) / length(x1$homicidio_mulher_pe_bra))

#---- MAPA HOMICIDIO FEMININO EM PE por municipio (media 2009 - 2013) ----#

# ver minimo e maximo da taxa de homicidio de mulheres 
# em PE e definir intervalos a serem usados no mapa
max(mapaviolencia_mulher_pe$taxa_homicidio_pe)
min(mapaviolencia_mulher_pe$taxa_homicidio_pe)

# definir intervalos a serem utilizados no mapa
intervalos <- c(7.5, 5, 2.5, 0)

# definir cortes e niveis
cortes <- cut(mapaviolencia_mulher_pe$taxa_homicidio_pe, intervalos, include.lowest = T)
niveis <- levels(cortes)

## MAPA ##

# definir cores dos niveis
levels(cortes) <- c("#ff0000","#b20000","#660000")

# plotar mapa
plot(geo_pernambuco, border = gray(0.9), lwd = .1,
     axes = F, las = 1, col = as.character(cortes))

# inserir legenda no mapa
legend("topright", niveis, fill = c("#ff0000","#b20000","#660000"), bty = "n", 
       title = "Homicidio de Mulheres p/ 10.000 hab", cex = 1)

##----MESOREGIOES----##

# visualizar niveis de MESOREGIAO
levels(mapaviolencia_mulher_pe$MESOREGIAO)

# criar objetos das mesoregioes
agreste <- subset(mapaviolencia_mulher_pe, MESOREGIAO == "AGRESTE PERNAMBUCANO")
zonadamata  <- subset(mapaviolencia_mulher_pe, MESOREGIAO ==  "MATA PERNAMBUCANA")
zonametro  <- subset(mapaviolencia_mulher_pe, MESOREGIAO == "METROPOLITANA DE RECIFE")
saofrancisco  <- subset(mapaviolencia_mulher_pe, MESOREGIAO == "SAO FRANCISCO PERNAMBUCANO")
sertao  <- subset(mapaviolencia_mulher_pe, MESOREGIAO == "SERTAO PERNAMBUCANO")

# cria variavel com media de cada mesoregiao
mapaviolencia_mulher_pe$media_mesoregiao[str_detect(MESOREGIAO, pattern = "AGRESTE PERNAMBUCANO")] <- mean(agreste$taxa_homicidio_pe)
mapaviolencia_mulher_pe$media_mesoregiao[str_detect(MESOREGIAO, pattern = "MATA PERNAMBUCANA")] <- mean(zonadamata$taxa_homicidio_pe)
mapaviolencia_mulher_pe$media_mesoregiao[str_detect(MESOREGIAO, pattern = "METROPOLITANA DE RECIFE")] <- mean(zonametro$taxa_homicidio_pe)
mapaviolencia_mulher_pe$media_mesoregiao[str_detect(MESOREGIAO, pattern = "SAO FRANCISCO PERNAMBUCANO")] <- mean(saofrancisco$taxa_homicidio_pe)
mapaviolencia_mulher_pe$media_mesoregiao[str_detect(MESOREGIAO, pattern = "SERTAO PERNAMBUCANO")] <- mean(sertao$taxa_homicidio_pe)

# arredondar, transformar em fator e definir niveis de MESOREGIAO
mapaviolencia_mulher_pe$media_mesoregiao2 <- round(mapaviolencia_mulher_pe$media_mesoregiao, 3)
mapaviolencia_mulher_pe$media_mesoregiao2 <- as.factor(mapaviolencia_mulher_pe$media_mesoregiao2)

niveis <- levels(mapaviolencia_mulher_pe$media_mesoregiao2)

## MAPA ##

# definir cores dos niveis
levels(mapaviolencia_mulher_pe$media_mesoregiao2) <-c("#ff0000","#cc0000","#7f0000",
                                                      "#4c0000", "#000000")

# plotar mapa
plot(geo_pernambuco, border = gray(0.9), lwd = .1, axes = F, las = 1, 
     col = as.character(mapaviolencia_mulher_pe$media_mesoregiao2))

# inserir legenda no mapa
legend("topright", niveis, fill = c("#ff0000","#cc0000","#7f0000","#4c0000", "#000000"),
       bty = "n", title = "Homicidio de Mulheres p/ 10.000 hab", cex = 1)


#----TESTE DE COMPARAÇÃO DE MEDIAS (ANOVA) POR MESOREGIAO----#

# realizar teste e atribuir resultados a um objeto
dados.anova <- aov(mapaviolencia_mulher_pe$taxa_homicidio_pe ~ 
                    mapaviolencia_mulher_pe$MESOREGIAO)

# visualizar resultados do teste
summary(dados.anova)