###########################################################
## PGM - Revista Dilemas - Dados de Encarceramento no BR ##
## Última atualização: 31/01/2018                        ##
###########################################################
rm(list=ls(all=TRUE)) # Limpar memória do R 
library(foreign); library(descr); library(car);
library(ggplot2); library(arm);

# Dados sobre encarceramento no Brasil - 2014 
setwd('C:/Users/Vinicius/Desktop/2018 Israel Bachini/BD') # Desktop

# --------------------------------------------------- #
# Lendo e padronizando as bases de dados - BD-NATASHA #
# --------------------------------------------------- #
# BD-Encarceramento.csv é uma atualização do BD encarceramento_corrigir_dados.csv
bd = read.table("BD-Encarceramento.csv", header= TRUE, dec=".", sep=";", nrows=27)
head(bd)
summary(bd)

# Variável dependente - Taxa de Presos por 100 mil hab.
y = bd$TAXAENCARCERAMENTO
summary(y); 
op=par(mfrow=c(2,2)) 
  hist(y); qqnorm(y); qqline(y);
  hist(log(y)); qqnorm(log(y)); qqline(log(y));
par(op)
shapiro.test(y)      # Teste de normalidade para y
shapiro.test(log(y)) # Teste de normalidade para log(y)

# Covariáveis (p = 17)
# BRAN.100 = Taxa de brancos
bd$BRAN.100
# EVA = Taxa de evangélicos
bd$EVA=bd$EVA/100
# FRONT = Dummy se o estado faz fronteira com outros países
bd$FRONT=as.factor(bd$FRONT)
# FUNDINC.100 = Taxa de população com ensino fundamental incompleto
bd$FUNDINC=bd$FUNDINC/100
# Índice GINI
bd$GINI=bd$GINI/1000
# H.100 = Taxa de homens
# HOMI.100 = Taxa de homicídios
bd$HOMI.100 = bd$HOMI.100/1000
# IDH
bd$IDHM=bd$IDHM/1000
# MIGESTICID = Taxa de migrantes
bd$MIGEST
# PIBBRASIL = Contriuição do estado para o PIB nacional
bd$PIBBRASIL
# TAXDES = Taxa de desemprego
bd$TAXDES
# TAXPOL = Taxa de policiais
bd$TAXPOL
# TAXJOV = Taxa de jovens
bd$TAXAENCJOVENS = bd$TAXAENCJOVENS/1000
# TXNEP = Taxa de negros ou pardos
bd$TXNEP #??? Muito altto
# TXPUPURB = Taxa de população urbana
bd$TXPOPURB
# VDILMA = Voto em Dilma
bd$VDILMA = as.factor(bd$VDILMA)
summary(bd)

# ----------------------------- #
# Análise Exploratória de Dados #
# ----------------------------- #
# Correlação entre as variáveis
corr = as.data.frame(cor(cbind(
    log(bd$TAXAENCARCERAMENTO), bd$IDHM, bd$GINI, bd$HOMI.100, bd$FUNDINC.100,
        bd$TXJOV, bd$TXNEP, bd$EVA, bd$TXPOPURB, bd$TAXDES, bd$TAXPOL, bd$H.100, 
        bd$PIBBRASIL,bd$DESARM, bd$VDILMA, bd$FRONT,bd$BRAN.100, bd$MIGESTCI)))
colnames(corr)=rownames(corr)=c('logtx','idh','gini','%homi','%fund',
        '%jov','%nep','%eva','%urb','%desemp','%policia','%H',
        'pibbr','desar','vdilma','front','%brancos','%miges')
fix(corr)
cc = (abs(corr)>.7)

# FIGURA 1
# Fig 1 a
op=par(mfrow=c(1,1))
  y = bd$TAXAENCARCERAMENTO
  ord = order(bd$TAXAENCARCERAMENTO, decreasing = T)
  par(mar=c(10,5,4,5)+.1)
  x=barplot(y[ord],  main='INFOPEN 2012', xlab='', cex.lab=1.2,
            ylab='Taxa de encarceramento (por 100 mil hab.)', xaxt='n')
  labs = bd$ESTADO[ord]
  text(cex=1, x=x+.5, y=-3, labs, xpd=TRUE, srt=45, pos=2)
par(op)  
  
# Fig 1 b
op=par(mfrow=c(1,1))
  boxplot(y,  main='INFOPEN - 2012', lwd=2, cex.lab=1.2,
          xlab='Estados brasileiros',ylab='Taxa de encarceramento (por 100 mil hab.)')  
par(op)
    
# Dados de 2012 e de 2014 -------------------------- #
# Dados 2014: Meu trabalho com a Fernanda Mencarelli #
bd2 = read.csv2("EstadosBR.csv", header = T) 
summary(bd2)
bd2$SIGLA = c('AC', 'AL', 'AP', 'AM', 'BH', 'CE', 'DF', 'ES', 'GO', 
              'MA', 'MT', 'MS', 'MG', 'PA', 'PB', 'BR', 'PE', 'PI',
              'RJ', 'RN', 'RS', 'RO', 'RR', 'SC', 'SP', 'SE', 'TO')

op=par(mfrow=c(1,2))
  hist(bd2$TAXA, main='INFOPEN 2014', xlab='taxa de encarceramento',
       ylab='probabilidade',probability = T)
  lines(density(bd2$TAXA))
  boxplot(bd2$TAXA, xlab='INFOPEN 2014')
par(op)  

op=par(mfrow=c(2,2))
  qqnorm(bd2$TAXA, main='Taxa original')
  qqline(bd2$TAXA)
  qqnorm(log(bd2$TAXA), main='Escala do log')
  qqline(log(bd2$TAXA))
  qqnorm(scale(bd2$TAXA), main='Scaled')
  qqline(scale(bd2$TAXA))
  qqnorm(sqrt(bd2$TAXA), main='Escala da raiz')
  qqline(sqrt(bd2$TAXA))
par(op)
shapiro.test(bd2$TAXA) # Teste de normalidade para log(bd2$TAXA)


# http://www.ibge.gov.br/apps/populacao/projecao/
bd2$Pop2004 = c(646548, 3049431, 577786, 3170740, 14075570, 8042368, 2278824,
                3397224, 5534201, 6135099, 2763537, 2273874, 19037702, 6904392,
                3612078, 10128262, 8484308, 3007333, 15374696, 3062933, 10628806,
                1515151, 392392, 5801932, 39824526, 1950985, 1285028)
summary(bd2$Pop2004)


# Site do Ministério da Justiça - Pop Prisional referente a dez 2014.
# http://www.justica.gov.br/seus-direitos/politica-penal/transparencia-institucional/estatisticas-prisional/relatorios-estatisticos-sinteticos
bd2$Presos2004 = c(2548, 2541, 1574, 3012, 7144, 10116, 7299, 5221, 6226,
                   2964, 7891, 6289, 7221, 6076, 6118, 10817, 15817, 1785,
                   23054, 2243, 22621, 4124, 972,  9570, 120601, 2142, 933)
bd2$tx2004 = (bd2$Presos2004/bd2$Pop2004)*100000
summary(bd2$tx2004)

ord = order(bd2$tx2004, decreasing = T)
op =par(mar=c(10,4,4,5)+.1)
  x=barplot(bd2$tx2004[ord],  ylab='taxa de encarceramento',
          main='INFOPEN 2014', xaxt='n')
  names=bd2$ESTADO[ord]
  text(cex=1, x=x+.05, y=-1.3, labs, xpd=TRUE, srt=45, pos=2)
par(op)

# PIB per capita Estados brasileiros - 2002-2014.
# http://saladeimprensa.ibge.gov.br/noticias?view=noticia&id=1&busca=1&idnoticia=3315
aux0 =c('Rondônia', 5147.41, 19462.61,
       'Acre',     4876.17, 17034.15,
       'Amazonas', 7353.15, 22373.36,
       'Roraima',  6736.70, 19608.40,
       'Pará',     4043.64, 15430.53,
       'Amapá',    5977.03, 17845.34,
       'Tocantins',4344.12, 17495.94,
       'Maranhão', 2718.05, 11216.37,
       'Piauí',    2440.70, 11808.08,
       'Ceará',    3712.24, 14255.05,
       'Rio Grande do Norte', 4709.83, 15849.33,
       'Paraíba',  3627.98, 13422.42,
       'Pernambuco',4426.56,16722.05,
       'Alagoas',  3962.88, 12335.44,
       'Sergipe',  5529.80, 16882.71,
       'Bahia',    4388.28, 14803.95,
       'Minas Gerais', 6703.46, 24917.12,
       'Espírito Santo', 8348.80, 33148.56,
       'Rio de Janeiro', 12414.77, 40767.26,
       'São Paulo', 13443.91, 42197.87,
       'Paraná', 8927.46, 31410.74,
       'Santa Catarina', 9745.87, 36055.90,
       'Rio Grande do Sul', 9423.79, 31927.16,
       'Mato Grosso do Sul', 7599.05, 30137.58,
       'Mato Grosso', 7265.37, 31396.81,
       'Goiás', 7307.95, 25296.60,
       'Distrito Federal', 24721.18, 69216.80)

aux =c('Rondônia', 'Acre', 'Amazonas', 'Roraima', 'Pará',
        'Amapá',  'Tocantins', 'Maranhão', 'Piauí', 'Ceará',
        'Rio Grande do Norte', 'Paraíba', 'Pernambuco', 'Alagoas',
        'Sergipe', 'Bahia', 'Minas Gerais', 'Espírito Santo',
        'Rio de Janeiro', 'São Paulo', 'Paraná', 'Santa Catarina', 
        'Rio Grande do Sul', 'Mato Grosso do Sul', 'Mato Grosso', 
        'Goiás', 'Distrito Federal')
aux1 =matrix(c(5147.41, 19462.61,
        4876.17, 17034.15,
        7353.15, 22373.36,
        6736.70, 19608.40,
        4043.64, 15430.53,
        5977.03, 17845.34,
        4344.12, 17495.94,
        2718.05, 11216.37,
        2440.70, 11808.08,
        3712.24, 14255.05,
        4709.83, 15849.33,
        3627.98, 13422.42,
        4426.56,16722.05,
        3962.88, 12335.44,
        5529.80, 16882.71,
        4388.28, 14803.95,
        6703.46, 24917.12,
        8348.80, 33148.56,
        12414.77, 40767.26,
        13443.91, 42197.87,
        8927.46, 31410.74,
        9745.87, 36055.90,
        9423.79, 31927.16,
        7599.05, 30137.58,
        7265.37, 31396.81,
        7307.95, 25296.60,
        24721.18, 69216.80), 2,27)
aux1=t(aux1)
pib = data.frame('Estado'=aux,'PIB2002'=aux1[,1],'PIB2014'=aux1[,2])
od = order(pib$Estado)
pib = pib[od,]
pib$aumento = log(pib$PIB2014/pib$PIB2002)

# ----------------------------------------------
# Taxa de aumento do PIB per capita de 2002-2014
bd2$PIBtxAUMENTO = pib$aumento
par(mfrow=c(1,3))
  boxplot(bd2$PIBtxAUMENTO, main='Variação do PIB per capita (2002-2014)')
  hist(bd2$PIBtxAUMENTO, nclass=15, main='Variação do PIB per capita (2002-2014)')
  barplot(bd2$PIBtxAUMENTO,names.arg = bd2$SIGLA, main='Variação do PIB per capita (2002-2014)')
par(op)

# -------------------------
# Taxa de mudança 2014-2004
bd2$r = log(bd2$TAXA/bd2$tx2004)
summary(bd2$r); sd(bd2$r)
op = par(mfrow=c(1,2))
  qqnorm(bd2$r);  qqline(bd2$r);
  qqnorm(log(bd2$r));  qqline(log(bd2$r));
par(op)
shapiro.test(bd2$r)      # Teste de normalidade para r.
shapiro.test(log(bd2$r)) # Teste de normalidade para log(r).


ord = order(bd2$r, decreasing = T)
op = par(mar=c(10,4,4,5)+.1, mfrow=c(1,1))
  x=barplot(bd2$r[ord],  ylab=paste('r=',expression(log(Tx2014/Tx2004))),
          cex.lab=1.2, main='Taxa de aumento do encarceramento',xaxt='n')
  labs = bd2$ESTADO[ord]
  text(cex=1, x=x, y=-0.01, labs, xpd=TRUE, srt=45, pos=2)
par(op)


# ----------------------------------#
# Comparando as taxas temporalmente #
# ----------------------------------#
boxplot(bd2$tx2004, bd$TAXAENCARCERAMENTO, bd2$TAXA, names=c('2005','2012','2014'))
cor(cbind(bd2$tx2004,bd$TAXAENCARCERAMENTO,bd2$TAXA))

# 2012-2014
a = lm(bd$TAXAENCARCERAMENTO ~ bd2$TAXA)
summary(a)
# Figura 2 b
op = par(mfrow=c(1,1))
  plot(bd$TAXAENCARCERAMENTO, bd2$TAXA, main='', xlab="2012", ylab="2014", cex.lab=1.2, lwd=3)
  text(x=bd$TAXAENCARCERAMENTO, y=bd2$TAXA+15, labels=bd2$SIGLA, col='red', cex=1)
  abline(a=0,b=1, lty=3, lwd=2)
  abline(a$coefficients, lty=3, lwd=2, col="blue")
op=par()
# 2004-2014
a = lm(bd2$TAXA ~ bd2$tx2004)
summary(a)
# Figura 2 a
op = par(mfrow=c(1,1))
  plot(bd2$tx2004, bd2$TAXA, main='', xlab="2004", ylab="2014", cex.lab=1.2, lwd=3)
  text(x=bd2$tx2004, y=bd2$TAXA+15, labels=bd2$SIGLA, col='red', cex=1)
  abline(a=0, b=1, lwd=2, lty=3)
  abline(a$coefficients, lty=3, lwd=2, col="blue")
op=par()
  

# ------------------------------------- #
# Modelos de Regressão Linear Múltipla  #
# ------------------------------------- #
# Todas
m1todas <- lm(log(TAXAENCARCERAMENTO)~ IDHM+GINI+HOMI.100+FUNDINC.100
              +TXJOV+TXNEP+EVA+TXPOPURB+TAXDES+TAXPOL+H.100 
              +PIBBRASIL+DESARM+VDILMA+FRONT+BRAN.100+MIGESTCID, data=bd)
summary(m1todas)
step(m1todas)
  
# Hipóteses da literatura/ mapa do encarceramento:
m5 <- lm(log(TAXAENCARCERAMENTO) ~  GINI + VDILMA + TXJOV + TXPOPURB, data=bd)
summary(m5)
plot(m5)
  

# ------------------- #
# IDEOLOGIA           #
# Evolução das taxas  #
# ------------------- #
# -------------------------------------------
# Tratamento de Ideologia - Braga et al. 2015
ideol = read.csv("BragaEtAl2015.csv")
bd2$CategoriaIdeologia = ideol$Categoria
bd2$ContIdeologia = ideol$Direita...Esquerda
summary(bd2$ContIdeologia); sd(bd2$ContIdeologia)
op=par(mfrow=c(1,2))
  plot(bd2$CategoriaIdeologia)
  plot(bd2$ContIdeologia, bd2$r)
  abline(lm(r~ContIdeologia,data=bd2), col='red')
op=par()  

Er = subset(bd2, as.numeric(ContIdeologia) < 0, select=c(r))
Hr = subset(bd2, as.numeric(ContIdeologia) ==0, select=c(r))
Dr = subset(bd2, as.numeric(ContIdeologia) >0 , select=c(r))
boxplot(c(Er,Hr,Dr), main='Aumento nas Taxas de Encarceramento',
        xlab='Posição política', ylab='r',
        names=c('Esquerda', 'Hibrido', 'Direita'))

# Modelo bayesiano gaussiano - Variáveis numéricas
mb1=bayesglm(r ~ bd2$ContIdeologia + PIBtxAUMENTO, data=bd2, family = gaussian())
summary(mb1)
mb2=bayesglm(r ~ bd2$ContIdeologia, data=bd2, family = gaussian())
summary(mb2)
# Modelo bayesiano - gamma
mb1g=bayesglm(r ~ bd2$ContIdeologia + PIBtxAUMENTO, data=bd2, family = Gamma(link='identity'))
summary(mb1g)
mb2g=bayesglm(r ~ bd2$ContIdeologia, data=bd2, family = Gamma(link='identity'))
summary(mb2g)


# ---------------------------
# Tarouco e Madeira 2013
tm = read.csv("TaroucoMadeira2013.csv")
bd2$ED = tm$Esquerda.direita
summary(bd2$ED); sd(bd2$ED)
bd2$LC = tm$Liberal.conservador
summary(bd2$LC); sd(bd2$LC)

op=par(mfrow=c(2,2))
  barplot(bd2$ED, main='Esquerda-direita')
  barplot(bd2$LC, main='Liberal-conservador')
  plot(bd2$ED, bd2$r, xlab='esquerda-direita',ylab='r')
  abline(lm(r~ED,data=bd2), col='red')
  plot(bd2$LC, bd2$r, xlab='liberal-conservador',ylab='r')
  abline(lm(r~LC,data=bd2), col='red')
par(op)  


# FIGURA 3
od = order(bd2$r, decreasing = T)
op=par(mar=c(10,6,4,5)+.1, mfrow=c(1,1))
  x = barplot(bd2$r[od], main='Taxas de aumento', ylim=c(0,2.5),
              ylab='r  (2004-2014)', cex.lab=1.2, xaxt='n')
  labs = bd2$ESTADO[ord]
  text(cex=1, x=x+.5, y=-.05, labs, xpd=TRUE, srt=45, pos=2)
  text(x=x, y=bd2$r[od]+.075, labels=bd2$ED[od], col='red', cex=1.2)
  text(x=x, y=bd2$r[od]-.075, labels=bd2$ContIdeologia[od], col='blue', cex=1.5)
  points(x, bd2$PIBtxAUMENTO[od], pch=16)
  axis(side=4)
  mtext(side=4, line=3, 'Taxa de aumento do PIB per capita (2002-2014)', cex=1.2)  
  #legend('topright', legend=c('Tarouco e Madeira 2013', 'Braga et al. 2015'), col=c('red','blue'))
par(op)


# Modelo bayesiano gaussiano - Variáveis numéricas
mb1=bayesglm(r ~ ED+LC+PIBtxAUMENTO+ED:LC, data=bd2, family = gaussian())
summary(mb1)
mb2=bayesglm(r ~ ED+LC+PIBtxAUMENTO, data=bd2, family = gaussian())
summary(mb2)
mb3=bayesglm(r ~ ED+LC+ED:LC, data=bd2, family = gaussian())
summary(mb3)
mb4=bayesglm(r ~ ED, data=bd2, family = gaussian()) # Melhor modelo pelo AIC
summary(mb4) 
# Modelo bayesiano - gamma
mb1g=bayesglm(r ~ ED+LC+PIBtxAUMENTO+ED:LC, data=bd2, family = Gamma(link='identity'))
summary(mb1g)
mb2g=bayesglm(r ~ ED+LC+PIBtxAUMENTO, data=bd2, family = Gamma(link='identity'))
summary(mb2g)
mb3g=bayesglm(r ~ ED+LC+ED:LC, data=bd2, family = Gamma(link='identity'))
summary(mb3g)
mb4g=bayesglm(r ~ ED, data=bd2, family = Gamma(link='identity')) # Melhor modelo pelo AIC
summary(mb4g) 



# Esquerda, centro, direita
e = subset(bd2, as.numeric(bd2$ED) < 0, select=c(r))
c = subset(bd2, as.numeric(bd2$ED) == 0, select=c(r))
d = subset(bd2, as.numeric(bd2$ED) > 0, select=c(r))
# Liberal, centro, conservador
lib = subset(bd2, as.numeric(bd2$LC) < 0, select=c(r))
neu = subset(bd2, as.numeric(bd2$LC) == 0, select=c(r))
con = subset(bd2, as.numeric(bd2$LC) > 0, select=c(r))
op=par(mfrow=c(1,2))
  boxplot(c(e,c,d), main='', xlab='Posição política', ylab='r', ylim=c(0,2),
        names =c('Esquerda', 'Centro', 'Direita'))
  boxplot(c(lib,neu,con), main='', xlab='Posição política', ylab='r',  ylim=c(0,2),
          names=c('Liberal', 'Neutro','Conservador'))
op=par()  


# Modelo bayesiano - Variáveis categóricas
mbc=bayesglm(r ~ as.factor(ED), data=bd2, family = gaussian())
summary(mbc)
# Modelo clássico 
mcc = lm(r ~ as.factor(ED), data=bd2)
summary(mcc)


# Comparação entre as posições políticas
id = as.data.frame(cbind(bd2$ContIdeologia,bd2$ED, bd2$LC))
colnames(id) =  c('ideologia', 'ed','lc')
cor(id)

summary(id)