Estatística espacial - GET/IME/UFF
Patrícia Lusié Velozo da Costa
2023-2
Dados de Área
Material baseado na aula do professor Guilherme A. Veloso (https://rpubs.com/GuilhermeV/1014494, https://rpubs.com/GuilhermeV/1017519 e https://rpubs.com/GuilhermeV/1018554)
1 Regiões brasileiras
O pacote geobr contém diferentes unidades geográficas.
Algumas unidades geográficas utilizadas pelo Instituto Brasileiro de Geografia e Estatística (IBGE):
- Estados: São as maiores unidades geográficas do Brasil.
- Mesorregiões: São divisões estaduais que agrupam municípios com características semelhantes.
- Regiões Imediatas: Subdivisões das mesorregiões que contêm municípios com características ainda mais semelhantes.
- Microrregiões: Subdivisões das mesorregiões que representam áreas geográficas menores e mais homogêneas em termos socioeconômicos.
- Municípios: As menores unidades administrativas do Brasil.
As “regiões intermediárias” são áreas que abrangem grupos de municípios com características geográficas, econômicas e sociais semelhantes. Elas são definidas como uma unidade intermediária entre as mesorregiões e os municípios individuais. Essa subdivisão é útil para fins de análise e planejamento regional, pois permite uma visão mais detalhada das disparidades e características das áreas dentro de uma mesorregião. Para entender as “regiões intermediárias” específicas dentro do contexto do Brasil e utilizá-las em pesquisas, estudos ou análises, você pode acessar os dados fornecidos pelo IBGE, que fornece informações detalhadas sobre as regiões intermediárias, incluindo características demográficas, econômicas e sociais.
A codificação de cidades e estados pode ser consultada no site IBGE. Já a codificação de regiões intermediárias e imediatas no site WIKIPEDIA.
As unidades geográficas brasileiras são codificadas em números:
- Rio de Janeiro (cidade): 3304557
- Rio de Janeiro (região geográfica imediata): 330001
- Rio de Janeiro (região geográfica intermediária): 3301
- Rio de Janeiro (estado): 33
#Pacote que permite o acesso aos shapefiles de diferentes unidades geográficas brasileiras.
if(!require(geobr)){ install.packages("geobr"); require(geobr)}
#Construção e customização de gráficos
if(!require(ggplot2)){ install.packages("ggplot2"); require(ggplot2)}
#gridExtra: colocar diferentes graficos numa mesma tela
if(!require(gridExtra)){ install.packages("gridExtra"); require(gridExtra)}
#paleta de cores
if(!require(RColorBrewer)){ install.packages("RColorBrewer"); require(RColorBrewer)}
#manipulação e análise de dados espaciais
if(!require(sf)){ install.packages("sf"); require(sf)}
#classes e métodos para dados espaciais
if(!require(sp)){ install.packages("sp"); require(sp)}
#construção das matrizes de vizinhança
if(!require(spdep)){ install.packages("spdep"); require(spdep)}
#manipulação e análise de dados
if(!require(tidyverse)){ install.packages("tidyverse"); require(tidyverse)}
#LISA
if(!require(rgeoda)){ install.packages("rgeoda"); require(rgeoda)}
#Definindo a região de interesse como sendo o país
Brasil<-read_country(year=2018,simplified = T, showProgress = F)
ggplot()+geom_sf(data=Brasil,fill="#2D3E50", color="orange")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#Definindo a região de interesse como sendo o país delimitando por regiões
Regiao<-read_region(year=2018,simplified = T, showProgress = F)
ggplot()+geom_sf(data=Regiao,fill="#2D3E50", color="orange")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#Definindo a região de interesse como sendo o país e o delimitando pelos estados
Estado<-read_state(code_state="all",year=2018,simplified = T, showProgress = F)
ggplot()+geom_sf(data=Estado,fill="#2D3E50", color="orange")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#Definindo a região de interesse como sendo o estado do RJ
EstadoRJ<-read_state(code_state=33,year=2018,simplified = T, showProgress = F)
ggplot()+geom_sf(data=EstadoRJ,fill="#2D3E50", color="orange")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#Definindo a região de interesse como sendo o país e o delimitando pelos municípios
Municipio<- read_municipality(code_muni = "all",year=2018,simplified = T, showProgress = F)
ggplot()+geom_sf(data=Municipio,fill="#2D3E50", color="orange")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#Definindo a região de interesse como sendo o estado do RJ e o delimitando pelos municípios
MunicipioRJ<- read_municipality(code_muni = 33,year=2018,simplified = T, showProgress = F)
ggplot()+geom_sf(data=MunicipioRJ,fill="#2D3E50", color="orange")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#Definindo a região de interesse como sendo o município do Rio
RioDeJaneiro<- read_municipality(code_muni = 3304557,year=2018,simplified = T, showProgress = F)
ggplot()+geom_sf(data=RioDeJaneiro,fill="#2D3E50", color="orange")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#Definindo a região de interesse como sendo o país e o delimitando pelas regiões imediatas
Immediate<- read_immediate_region(code_immediate = "all",year=2020,simplified = T, showProgress = F)
ggplot()+geom_sf(data=Immediate,fill="#2D3E50", color="orange")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#Definindo a região de interesse como sendo o país e o delimitando pelas regiões intermediárias
Intermediate<- read_intermediate_region(code_intermediate = "all",year=2020,simplified = T, showProgress = F)
ggplot()+geom_sf(data=Intermediate,fill="#2D3E50", color="orange")+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
#Definindo a região de interesse como sendo o estado do RJ e o delimitando pelas regiões intermediárias
IntermediateRJ<- read_intermediate_region(code_intermediate = 33,year=2020,simplified = T, showProgress = F)
ggplot()+geom_sf(data=IntermediateRJ,fill=alpha("grey", 0.7), color="orange")+
geom_sf_text(data = IntermediateRJ, aes(label = name_intermediate),color="black",size=3)+
theme_bw()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
1.1 Exercício:
1 - Criar o shapefile dos municípios de São Paulo e representar em um mapa.
2 - Criar o shapefile do município de Niterói e representar em um mapa.
2 Mapa coroplético
Baixe os dados em: https://www.dropbox.com/scl/fi/j4f7lspc29u5rc2bvq2lg/DadosMG.xlsx?rlkey=7jsl3pf8ao2x11v5mm5clnvc4&dl=0 .
if(!require(openxlsx)){ install.packages("openxlsx"); require(openxlsx)}
## Leitura dos dados
Dados_MG <- read.xlsx("DadosMG.xlsx")
#Transformando os dados em uma variável não espacial
#Dados_MG$taxa <- rnorm(nrow(Dados_MG), mean(Dados_MG$taxa), sd(Dados_MG$taxa))
## Leitura do Shape dos municípios Mineiros
MG<-read_municipality(code_muni=31,year=2018,simplified=T)
## Função utilizada para categorizar os pontos de corte
age.cat <- function(x, breaks,
sep = "-", above.char = "+",min,max) {
labs <- c(paste(c(min,breaks)[-(length(breaks)+1)],
breaks-0.01,
sep = sep),
paste(breaks[length(breaks)], max, sep = "-"))
cut(x, breaks = c(min,breaks,Inf),
right = FALSE, labels = labs)
}
## Pontos de corte
breaks <- round(quantile(Dados_MG$taxa,probs = c(0.05,0.20,0.40,0.60,0.80,0.95)),2)
## Escolha da paleta de cores
mycolors <- rev(brewer.pal(7, "RdBu"))
## Categorizando os pontos de corte
Dados_MG$Categoria <- age.cat(round(Dados_MG$taxa,2),breaks,min=min(round(Dados_MG$taxa,2)),max=max(round(Dados_MG$taxa,2)))
## Merge do arquivo contendo as taxas com o shapefile
Dados_MG <- merge(MG,Dados_MG,by="code_muni",all.x=T)
## Construção do mapa
Mapa_MG <- ggplot() + geom_sf(data=Dados_MG,aes(fill=Categoria),colour="black") +
theme_bw() +
scale_fill_manual(values = mycolors,name="Taxa de mortalidade padronizada \npor idade / 100 mil habitantes 2009-2011") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "bottom")
Mapa_MG
## Categorizando os pontos de corte
breaks2 <- round(seq(min(Dados_MG$taxa), max(Dados_MG$taxa), l=8),2)[2:7]
Dados_MG$Categoria2 <- age.cat(round(Dados_MG$taxa,2),breaks2,min=min(round(Dados_MG$taxa,2)),max=max(round(Dados_MG$taxa,2)))
Mapa_MG2 <- ggplot() + geom_sf(data=Dados_MG,aes(fill=Categoria2),colour="black") +
theme_bw() +
scale_fill_manual(values = mycolors,name="Taxa de mortalidade padronizada \npor idade / 100 mil habitantes 2009-2011") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "bottom")
grid.arrange(Mapa_MG,Mapa_MG2, nrow = 2)
#Os comandos a seguir permitem a construção do mapa dinâmico:
if(!require(ggiraph)){ install.packages("ggiraph"); require(ggiraph)}
## O que será mostrado ao encostar em cada município
Dados_MG$tooltip <- c(paste0(Dados_MG$name_muni, "\n ", round(Dados_MG$taxa,2)))
## Construção do mapa
gg <- ggplot(Dados_MG) + geom_sf(aes(fill=Categoria)) +
geom_sf_interactive(aes(fill=Categoria, tooltip = tooltip, data_id = code_muni),color="black",size=0.1)+ theme_bw() +
scale_fill_manual(values = mycolors,name="Taxa de Mortalidade padronizada \npor idade / 100 mil habitantes 2009-2011") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "bottom")
girafe(ggobj = gg)3 Centróides
## Obtenção dos centroides
sf_cent_MG <- st_centroid(Dados_MG)
## Mapa dos centroides
ggplot() +
geom_sf(data = Dados_MG, fill = 'white') +
geom_sf(data = sf_cent_MG, color = 'red') + theme_bw() +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank())
4 Matrizes de Vizinhança
## Mudando de arquivo tipo SF para SP
Dados_MG_sp <- sf:::as_Spatial(Dados_MG$geometry)
## Calculando os vizinhos pela contiguidade Queen
Queen_MG<-spdep::poly2nb(Dados_MG_sp, queen=T)
## Construindo a matriz de vizinhança Queen
WQueen_MG<-matrix(0,nrow(Dados_MG),nrow(Dados_MG))
for(i in 1:nrow(Dados_MG))
WQueen_MG[i,Queen_MG[[i]]]=1
## Calculando os vizinhos pela contiguidade Rook
Rook_MG<-spdep::poly2nb(Dados_MG_sp, queen=F)
## Construindo a matriz de vizinhança Rook
WRook_MG<-matrix(0,nrow(Dados_MG),nrow(Dados_MG))
for(i in 1:nrow(Dados_MG))
WRook_MG[i,Rook_MG[[i]]]=1
## Plotando as duas estruturas de vizinhança lado a lado
par(mfrow=c(1,2))
sp::plot(Dados_MG_sp, border="lightgrey",main="QUEEN")
sp::plot(Queen_MG, coordinates(Dados_MG_sp), pch = 19, cex = 0.6, add = TRUE, col= "blue")
sp::plot(Dados_MG_sp, border="lightgrey",main="ROOK")
sp::plot(Rook_MG, coordinates(Dados_MG_sp), pch = 19, cex = 0.6, add = TRUE, col= "blue")
5 Médias móveis
## Criando a matriz de vizinhança Queen normalizada
NQueen_MG<-WQueen_MG/apply(WQueen_MG,1,sum)
## Extraindo o vetor de taxas
taxa<-Dados_MG$taxa
## Médias Móveis
Movel_MG<-NQueen_MG%*%taxa
## Inserindo as médias móveis no shapefile
Dados_MG$Movel<-Movel_MG
breaks <- round(quantile(c(Dados_MG$Movel),probs = c(0.05,0.20,0.40,0.60,0.80,0.95)),2)
Dados_MG$Categoria2 <- age.cat(round(Dados_MG$Movel,2),breaks,min=min(round(Dados_MG$Movel,2)),max=max(round(Dados_MG$Movel,2)))
Mapa_Movel_MG<-ggplot() + geom_sf(data=Dados_MG,aes(fill=Categoria2),colour="black") + theme_bw() +
scale_fill_manual(values = mycolors,name="Estimador de média móvel \n da taxa de Mortalidade padronizada \npor idade / 100 mil habitantes 2009-2011") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "bottom")
grid.arrange(Mapa_MG,Mapa_Movel_MG, nrow = 2)
6 Índice de Moran
6.1 Global
## Transformando os pesos em formato nb
listw<- nb2listw(Queen_MG, style = "W")
## Cálculo e teste do Indice Global de Moran
GlobalMoran<-moran.test(taxa,listw,alternative="two.sided")
GlobalMoran##
## Moran I test under randomisation
##
## data: taxa
## weights: listw
##
## Moran I statistic standard deviate = 11.556, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.2436303976 -0.0011737089 0.00044874586.2 Diagrama de Espalhamento de Moran
## Plot do diagrama de espalhamento de Moran
par(mfrow=c(1,1))
moran.plot(taxa,listw,label=Dados_MG$name_muni)
6.3 Índice Moran Local (LISA)
## Contiguidade Queen
queen_w_MG <- queen_weights(Dados_MG)
## Construção e teste de hipóteses do LISA
lisa <- rgeoda::local_moran(queen_w_MG, as.data.frame(Dados_MG$taxa))
## Colocar no shapefile a classificação dos clusters 0, 1, 2, 3 ou 4.
Dados_MG$LISA<-lisa_clusters(lisa)
## Colocar no shapefile os p-valor
Dados_MG$LISA_PVALUE<-lisa_pvalues(lisa)
## Substituindo os clusters numéricos por caracteres
Dados_MG$LISA<-gsub("0","Not Significant",Dados_MG$LISA)
Dados_MG$LISA<-gsub("1","High-High",Dados_MG$LISA)
Dados_MG$LISA<-gsub("2","Low-Low",Dados_MG$LISA)
Dados_MG$LISA<-gsub("3","Low-High",Dados_MG$LISA)
Dados_MG$LISA<-gsub("4","High-Low",Dados_MG$LISA)
## Cores
mycolors <- c("red","pink","lightblue","blue","white")
## Mapa
ggplot() + geom_sf(data=Dados_MG,aes(fill=LISA),color="black") +
theme_bw() +
scale_fill_manual(values = mycolors) +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "bottom")
#mapa dinâmico
## O que será mostrado ao encostar em cada município
Dados_MG$tooltip <- c(paste0(Dados_MG$name_muni, "\n ", Dados_MG$LISA_PVALUE))
## Construção do mapa
gg <- ggplot(Dados_MG) +
geom_sf_interactive(aes(fill=LISA, tooltip = tooltip, data_id = code_muni),color="black",size=0.1)+ theme_bw() +
scale_fill_manual(values = mycolors,name="Taxa de Mortalidade padronizada \npor idade / 100 mil habitantes 2009-2011") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "bottom")
girafe(ggobj = gg)