### Aula 5 - Construcao de matrizes de vizinhanca, Moran local e Taxas Bayesianas Empíricas no R

# limpando objetos

#rm(list=ls())

# setando diretorio de trabalho

getwd()
#setwd("U:/Chico2/aula4")

########################################################

# carregando imagem

#load("aula4.Rdata")

#######################################################

# Instalar o pacote "maptools"

# carregar os pacotes

library(maptools)
library (spdep)
library(sf)


#######################################################

## 1 - Matriz de vizinhanca de um shape de poligonos - viz queen e rook

### olhar o html do police

police.map <- st_read("./bd_aula5/police.shp") # lendo mapa

plot(st_geometry(police.map))
names(police.map)
str(police.map) ## ou olhar no 'Environment'
head(police.map) ## ou olhar no 'Environment'

#### 1 - Diferentes tipos de matrizes de vizinhanca

## 1.1 criando diferentes matrizes de contig

plot(st_geometry(police.map))

# queen
pol.gal.q <- poly2nb(as(police.map, "Spatial"),queen=TRUE) # matriz de viz cont tipo queen
print(pol.gal.q)
str(pol.gal.q) ## ou olhar no 'Environment'

# rook
pol.gal.r <- poly2nb(as(police.map, "Spatial"),queen=F) # matriz de viz cont tipo rook
print(pol.gal.r)
str(pol.gal.r)

# diferenca entre as duas matrizes

pol.gal.dif <- diffnb(pol.gal.q,pol.gal.r)
print(pol.gal.dif)
str(pol.gal.dif)

# plotando as matrizes de viz e identif as diferencas entre elas

# obtendo o mapa de centroides

pol.centr <- st_geometry(st_centroid(police.map))
plot(st_geometry(police.map))
plot(st_geometry(pol.centr),pch=3, col='red', add=TRUE)
str(pol.centr)

## plotando somente os centroides
plot(st_geometry(pol.centr),pch=3, col='red')

# obtendo a matriz com as coord dos centroides

pol.coord <- st_coordinates(pol.centr)
class(pol.coord)

# plotando matriz de contig queen

plot(st_geometry(police.map))
plot(pol.gal.q, pol.coord, add=TRUE)

# plotanto matriz de cont rook

plot(st_geometry(police.map))
plot(pol.gal.r, pol.coord, add=TRUE)

# rodar os comandos seguintes juntos para ver a dif entre as duas matrizes
plot(st_geometry(police.map, border = "grey"))
plot(pol.gal.r, pol.coord, add = TRUE)
plot(pol.gal.q, pol.coord, add = TRUE)
plot(pol.gal.dif, pol.coord, add = TRUE, col = "red")

# matriz de pesos

pol.gal.q.pes <- nb2listw(pol.gal.q)
print(pol.gal.q.pes)
str(pol.gal.q.pes)

### 1.2 - Vizinhanca baseada na distancia

# Escolha dos k vizinhos mais proximos
# na grande maioria dos casos - matrizes assimetricas

# 'knearneigh' produz uma forma intermediaria convertida
# em um objeto 'nb' pelo comando 'knn2nb'
# 'knearneigh' pode tb ser utilizada com um argumento
# 'longlat' para lidar com coord geograficas.

### para um shape de pologonos - usar as coordenadas
# como ja foi feito acima
# coords.pol<-coordinates(police_map)

# 1 vizinho

# usando os dois comando 'knearneigh' e 'knn2nb' em separado
pol1 <- knearneigh(pol.coord, k = 1)
pol.1viz <- knn2nb(pol1)
print(pol.1viz)
str(pol.1viz)

# matriz de pesos
pol.1vz.peso <- nb2listw(pol.1viz)
str(pol.1vz.peso)

# usando os dois comando em conj
pol.1viz <- knn2nb(knearneigh(pol.coord, k = 1))
print(pol.1viz)
str(pol.1viz)

# 2 vizinhos
pol.2viz <- knn2nb(knearneigh(pol.coord, k = 2))
print(pol.2viz)
str(pol.2viz)

# matriz de pesos
pol.2vz.peso <- nb2listw(pol.2viz)
str(pol.2vz.peso)

# 4 vizinhos
pol.4viz <- knn2nb(knearneigh(pol.coord, k = 4))
print(pol.4viz)
str(pol.4viz)

# matriz de pesos
pol.4vz.peso <- nb2listw(pol.4viz)
str(pol.4vz.peso)

## 1.3 - Matriz de vizinhanca de um shape de pontos - k vizinhos

# lendo mapa de pontos

ret.diag <- st_read("./bd_aula5/ret_diag_79.shp") 
class(ret.diag)
plot(st_geometry(ret.diag))
names(ret.diag)

# construindo matriz de vizinhan?a de 6 viz mais pr?ximos

ret.diag.6v <- knearneigh(as(ret.diag, "Spatial"),k=6)
print(ret.diag.6v)
ret.diag.6v.nb <- knn2nb(ret.diag.6v)
str(ret.diag.6v.nb)

# matriz de pesos

ret.diag.6v.pes <- nb2listw(ret.diag.6v.nb) # convert lista de viz em pesos espaciais
str(ret.diag.6v.pes)

moran.test(ret.diag$RET_DIAG,ret.diag.6v.pes, zero.policy=NULL)

###############################################

### 3 Moran local

# lista de comandos do pacote 'spdep'

??localmoran

#localmoran(x, listw, zero.policy=NULL, na.action=na.fail,
#           alternative = "greater", p.adjust.method="none", mlvar=TRUE,
#           spChk=NULL)


pol.moran.local <- localmoran(police.map$CRIME,listw = pol.gal.q.pes,alternative="two.sided")
class(pol.moran.local)
pol.moran.local
head(pol.moran.local)

# ordenando a matriz 

head(pol.moran.local)
pol.moran.local[order(pol.moran.local[,"Pr(z != E(Ii))"]),] # 3 areas com p < 0.05

# juntando o objeto 'sf' com a matriz

police.moran.local <- cbind(police.map,pol.moran.local)
class(police.moran.local)
str(police.moran.local)

# criando a  coluna 'moran.local' e dizendo q seus valores sao nulos


police.moran.local$moran.local <- 0

# classificando o moran local segundo as categorias 1 (alto-alto); 2 (baixo-baixo); 3 (baixo-alto) e 4 (alto-baixo)

str(police.moran.local)

# m?dia da vari?vel CRIME

mean(police.moran.local$CRIME) # 194.5

### poder?amos usar tb os valores da coluna Z.Ii q rerpresenta a m?dia dos vizinhos

# atribuindo os valores corretos

head(police.moran.local) ## para copiar nome da coluna 'Pr.z....E.Ii..'
police.moran.local$moran.local <- ifelse(police.moran.local$Pr.z....E.Ii..>=0.05 , 0,
                         ifelse(police.moran.local$CRIME>194.5 & police.moran.local$Ii > 0, 1,
                          ifelse (police.moran.local$CRIME<194.5 & police.moran.local$Ii>0,2,
                          ifelse(police.moran.local$CRIME<194.5 & police.moran.local$Ii<0,3,4))))

# mostrando o resultado
str(police.moran.local)
police.moran.local$moran.local

print(police.moran.local[c(16,22,25,26,28)],n=82)

# mapeando

plot(police.moran.local["moran.local"])

# exportando o shape

st_write(police.moran.local,"./bd_aula5/police.moran.local.shp")

### Moran plot

# moran.plot(x, listw, zero.policy=NULL, spChk=NULL, labels=NULL,
#           xlab=NULL, ylab=NULL, quiet=NULL, ...)

moran.plot(police.map$CRIME,listw = pol.gal.q.pes)

#####################################

### 4 - Taxa bayesiana empirica 

## 4.1 Global

# lendo shape

mama <- st_read("./bd_aula5/munic_esp_ob_pop_2010.shp") # lendo mapa

plot(st_geometry(mama))
names(mama)
str(mama)
head(mama)

# taxa bruta

mama$tx_br <- (mama$obt/mama$p_tot)*100000
head(mama)

# taxa bayesiana global

# EBest(n, x, family="poisson")
# n: casos
# x: popula??o

tx.gl <- EBest(mama$obt, mama$p_tot, family="poisson")
class(tx.gl)
str(tx.gl)

# juntando a taxa bayesiana global com o shape

mama <- cbind(mama,tx.gl$estmm)

str(mama)

# taxa bayesiana por 100000 hab

mama$tx.gl.100 <- (mama$tx.gl.estmm)*100000

print(mama[c(6,9)],n=100)

# todas as taxas foram em direcao a media - excesso de zero!!

plot(mama["tx_br"])
plot(mama["tx.gl.100"])

### 4.2 Taxa bayesiana local

# EBlocal(ri, ni, nb, zero.policy = NULL, spChk = NULL, geoda=FALSE)

# ri: casos
# ni: populacao
# nb: matriz de vizinhanca
# zero.policy: default NULL, use global option value; if TRUE assign zero to the lagged value
# of zones without neighbours, if FALSE assign NA
# geoda=T: para obter mesmo resultado q GeoDa

str(mama)

# matriz de vizinhanca - contig

# queen
mama.gal.q <- poly2nb(as(mama, "Spatial"),queen=TRUE) # matriz de viz cont tipo queen
print(mama.gal.q)
str(mama.gal.q)

# matriz de pesos

mama.gal.q.pes <- nb2listw(mama.gal.q) # erro! municipios sem vizinhos: ilhas

??nb2listw # olhar o help para ver os argumentos da fun??o - precisa do 'zero.policy=T', para permitir viz nula

mama.gal.q.pes <- nb2listw(mama.gal.q,zero.policy=T)
print(mama.gal.q.pes) # nao roda!
str(mama.gal.q.pes)

# outras possibilidades: editar a matriz, usar o snap (nao muito interessante, pq sao ilhas) 
# ou usar matriz de dist ou viz mais prox

??EBlocal

tx.lc <- EBlocal(mama$obt, mama$p_tot, mama.gal.q) # uso matriz de viz e n?o de pesos
str(tx.lc)

# juntando a taxa bayesiana local com o shape

mama <- cbind(mama,tx.lc$est)

str(mama)

# taxa bayesiana local por 100000 hab

mama$tx.lc.100 <- (mama$tx.lc.est)*100000
str(mama)

print(mama[c(6,11)],n=100)

plot(mama["tx_br"])
plot(mama["tx.lc.100"])

### 5 - Usando outro pacote leitura do shape - rgdal

library(rgdal)

# lendo shape

mama.rgdal <- readOGR('./bd_aula5', 'munic_esp_ob_pop_2010')
plot(mama.rgdal)
slotNames(mama.rgdal)
str(mama.rgdal@data)
mama.rgdal@data
mama.rgdal@bbox

# matriz de vizinhanca

#library(spdep)

mama.r.gal.q<-poly2nb (mama.rgdal, queen=TRUE)
print(mama.r.gal.q)
str(mama.r.gal.q)

# matriz de pesos

mama.r.gal.q.pes <- nb2listw(mama.r.gal.q,zero.policy = TRUE)
print(mama.r.gal.q.pes)
str(mama.r.gal.q.pes)

############################################################

## salvando imagem
save.image("aula5.Rdata")

############################################################