library(spdep)
library(maptools)
library(RColorBrewer)

base=read.shape("neigh-parana.shp")

dadospolys=Map2poly(base)


dadosnb=poly2nb(dadospolys)

dadoscents=get.Pcent(base)

plot.polylist(dadospolys,border="blue")
points(dadoscents,col="red")

plot.polylist(dadospolys,border="green")
plot(dadosnb,dadoscents,add=TRUE,col="blue")

nome=base$att.data$NOME_MUNIC
codigo=base$att.data$ID_IBGE
nascidos=base$att.data$NASCIDOS
obitos=base$att.data$OBITOS

dados=data.frame(nome,codigo,nascidos,obitos)

res=probmap(obitos,nascidos)

#Fazendo análise exploratória da variável raw*1000


hist(res$raw*1000)

#variável raw
brks=c(-Inf,4,13,20,28,Inf)
nclr=length(brks)-1
plotclr=brewer.pal(nclr,"Greens")
cols=plotclr
plot(dadospolys,col=cols[findInterval(res$raw*1000,brks)],forcefill=FALSE,
 ylab='Longitude',xlab='Latitude',cex.axis=0.8)
text(-51.5,-22,'Taxa de mort. menores de 1 ano - 2002')
legend(c(-51,-48),c(-27.5,-26.5),fill=cols,legend=leglabs(brks),
cex=0.6,ncol=2,bty="n")

#outra forma de utilizar as cores
#brks=c(-Inf,4,13,20,28,Inf)
#cols=c("#FDDBC7","#EF8A62","#D1E5F0", "#67A9CF","#2166AC")
#plot(dadospolys,col=cols[findInterval(res$raw*1000,brks)],forcefill=FALSE,
# ylab='Longitude',xlab='Latitude',cex.axis=0.8)
#text(-51.5,-22,'Taxa de mort. menores de 1 ano - 2002')
#legend(c(-51,-48),c(-27.5,-26.5),fill=cols,legend=leglabs(brks),
#cex=0.6,ncol=2,bty="n")

#variável relRisk
brks=c(-Inf,23,76,118,163,Inf)
nclr=length(brks)-1
plotclr=brewer.pal(nclr,"Greens")
cols=plotclr
plot(dadospolys,col=cols[findInterval(res$relRisk,brks)],forcefill=FALSE,
 ylab='Longitude',xlab='Latitude',cex.axis=0.8)
text(-51.5,-22,'Risco relativo')
legend(c(-51,-48),c(-27.5,-26.5),fill=cols,legend=leglabs(brks),
cex=0.6,ncol=2,bty="n")

#variável pmap
brks=c(0,0.05,0.1,0.90,0.95,1)
nclr=length(brks)-1
plotclr=brewer.pal(nclr,"Oranges")
cols=plotclr
plot(dadospolys,col=cols[findInterval(res$pmap,brks)],forcefill=FALSE,
 ylab='Longitude',xlab='Latitude',cex.axis=0.8)
text(-51.5,-22,"Distribuicao de Poisson - pmap")
legend(c(-51,-48),c(-27.5,-26.5),fill=cols,legend=leglabs(brks),
cex=0.6,ncol=2,bty="n")

#regressão de Poisson
res=glm(obitos~offset(log(nascidos)),data=dados,family="quasipoisson")
stdres=rstandard(res)
brks=c(-Inf,-2,-1,1,2,+Inf)
nclr=length(brks)-1
plotclr=brewer.pal(nclr,"Greens")
cols=plotclr
plot(dadospolys,col=cols[findInterval(stdres,brks)],forcefill=FALSE,
ylab='Longitude',xlab='Latitude',cex.axis=0.8)
text(-51.5,-22,'Regressao de Poisson')
legend(c(-51,-48),c(-27.5,-26.5),fill=cols,legend=leglabs(brks),
cex=0.6,ncol=2,bty="n")


par(mfrow=c(1,2))

#Bayes empírico
res=EBest(obitos,nascidos)
brks=c(-Inf,4,13,20,28,Inf)
nclr=length(brks)-1
plotclr=brewer.pal(nclr,"Oranges")
cols=plotclr
plot(dadospolys,col=cols[findInterval(res$estmm*1000,brks)],forcefill=FALSE,
 ylab='Longitude',xlab='Latitude',cex.axis=0.8)
text(-51.5,-22,'Bayes empirico')
legend(c(-51,-48),c(-27.5,-26.5),fill=cols,legend=leglabs(brks),
cex=0.6,ncol=2,bty="n")


#Bayes empírico local
res=EBlocal(obitos,nascidos,dadosnb,zero.policy=TRUE)
brks=c(-Inf,4,13,20,28,Inf)
nclr=length(brks)-1
plotclr=brewer.pal(nclr,"Oranges")
cols=plotclr
sub=is.finite(res$est)
plot(dadospolys,forcefill=FALSE,ylab='Longitude',xlab='Latitude',cex.axis=0.8)
legend(c(-51,-48),c(-27.5,-26.5),fill=cols,legend=leglabs(brks),
cex=0.6,ncol=2,bty="n")
plot(subset(dadospolys,sub),col=cols[findInterval(res$est[sub]*1000,brks)],
   add=TRUE, forcefill=FALSE)
points(dadoscents[!sub,],pch=8)
text(-51.5,-22,'Bayes empirico local')