### 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") ############################################################