ANÁLISE DESCRITIVA ESPACIAL COM O PACOTE geobr (PARTE III)

1ª ESCOLA DE VERÃO

PROF. GUILHERME AUGUSTO VELOSO

UFF - 20, 22 E 24 DE MARÇO DE 2023

Uma forma usual de visualização de dados agregados por áreas é o uso de mapas coloridos segundo o padrão espacial do fenômeno; são os chamados mapas temáticos, nos quais as cores indicam diferentes classes ou faixas de valores da variável mapeada.

A maioria dos usuários limita-se a essas operações de visualização, tirando conclusões preliminares, de modo intuitivo. Contudo, é possível ir muito além dessa análise superficial.

Quando se visualiza um padrão espacial, é muito útil traduzi-lo em questões objetivas, como as apresentadas a seguir:

• O padrão observado é aleatório ou apresenta uma agregação definida?
• Essa distribuição pode ser associada a causas mensuráveis?
• Existem agrupamentos de áreas com padrões diferenciados na região de estudo?

Tais questões são objeto das técnicas de análise exploratória de dados agregados por área. Dessa forma, a terceira parte do curso tem por objetivo explorar os padrões de dependência espacial com a introdução do semivariograma e os índices de Moran global e local no contexto uni e bivariado.

Precisaremos dos pacotes:

• geobr: Permite o acesso aos shapefiles de diferentes unidades geográficas brasileiras.

• ggplot2: construção e customização de mapas coropléticos.

• RColorBrewer: Paleta de Cores.

• openxlsx: leitura de arquivos no formato .xlsx.

• ggiraph: permite a construção de mapas dinâmicos.

• sp: classes e métodos para dados espaciais.

• gstat: construção do semivariograma.

• spdep: construção das matrizes de vizinhança

• rgeoda: cálculo do LISA

Segue o código para a instalação dos pacotes:

if(!require(geobr)){ install.packages("geobr"); require(geobr)}
if(!require(ggplot2)){ install.packages("ggplot2"); require(ggplot2)}
if(!require(RColorBrewer)){ install.packages("RColorBrewer"); require(RColorBrewer)}
if(!require(openxlsx)){ install.packages("openxlsx"); require(openxlsx)}
if(!require(ggiraph)){ install.packages("ggiraph"); require(ggiraph)}
if(!require(sp)){ install.packages("sp"); require(sp)}
if(!require(gstat)){ install.packages("gstat"); require(gstat)}
if(!require(spdep)){ install.packages("spdep"); require(spdep)}
if(!require(rgeoda)){ install.packages("rgeoda"); require(rgeoda)}

1 Semivariograma

A função semivariância é utilizada na estatística espacial para expressar a variabilidade numa distância definida previamente. Considere a investigação de um fenômeno de interesse e que \(z(x)\) seja o valor do fenômeno de interesse em uma área e \(z(x+h)\) o valor do fenômeno de interesse em uma outra área, cujos centroides estão separados por uma distância \(h\). A semivariância experimental é definida pela seguinte expressão:

em que \(N(h)\) é o número de pares experimentais.

O semivariograma constitui-se no gráfico das semivariâncias situadas a intervalos regulares e sua representação geral é dada pela seguinte figura

• Alcance: distância dentro da qual as amostras apresentam-se correlacionadas espacialmente.

• Patamar: é o valor do semivariograma correspondente a seu alcance. Deste ponto em diante, considera-se que não existe mais dependência espacial entre as amostras, porque o semivariograma se torna invariante com a distância.

• Efeito Pepita: é o valor que revela a semivariância para a menor distância entre as amostras

Os comandos a seguir calculam o semivariograma para os dados de Minas Gerais:

## Leitura dos dados 
Dados_MG <- read.xlsx("DadosMG.xlsx")  

## Leitura do Shape dos municípios Mineiros
MG<-read_municipality(code_muni=31,year=2018,simplified=T)

## Merge do arquivo contendo as taxas com o shapefile
Dados_MG <- merge(MG,Dados_MG,by="code_muni",all.x=T)

## Obtenção dos centróides
sf_cent_MG <- st_centroid(Dados_MG)
## Extraindo a taxa 
taxa_MG<-sf_cent_MG$taxa
## Mudando o formato do dado espacial para aplicar a função variogram
sf_cent <- sf:::as_Spatial(sf_cent_MG$geometry)
## Construção do Variograma
Var_MG <- gstat::variogram(taxa ~ 1, sf_cent_MG)
plot(Var_MG,type="o")

## "Experimentando" modelos teóricos vgm("Exp" ou "Sph" ou "Gau")
plot(Var_MG,gstat::fit.variogram(Var_MG, gstat::vgm("Exp")))

1.1 Atividade

Construa o semivariograma para os dados do estado do Rio de Janeiro.

2 Índice Moran Global

Um aspecto fundamental da análise exploratória espacial é a caracterização da dependência espacial, mostrando como os valores estão correlacionados no espaço.

Uma vez determinada a estrutura espacial de análise, expressa pela matriz \(W\), qualquer medida particular de autocorrelação pode ser concebida mediante a definição de um modo de se mensurar a diferença entre valores do atributo associado às localizações (áreas).

A medida mais usada para este fim é o Índice de Moran Global (I), que é a versão espacial do coeficiente de correlação linear de Pearson e que representa o coeficiente de correlação para o relacionamento entre os valores de uma variável espacial (atributo) e o valor médio desta variável, ponderado pela estrutura de vizinhança.

O Índice Moran Global I é dado por:

\(\Large I=\frac{\sum_{i=1}^n\sum_{j=1}^nw_{ij}(z_i-\bar{z})(z_j-\bar{z})}{\sum_{i=1}^n(z_i-\bar{z})^2}\)

Na equação acima:

• \(n\) é o número de áreas;
• \(z_i\) o valor do atributo considerado na área \(i\);
• \(\bar{z}\) é o valor médio do atributo na região de estudo ;
• \(w_{ij}\) os elementos da matriz normalizada de proximidade espacial.

Observa-se na equação do índice I:

1- O numerador da fração é um termo de covariância. O cálculo do produto das diferenças entre o valor do atributo, em duas áreas, e a média global determina o quanto essas diferenças covariam. Se ambos os valores estão do mesmo lado da média (acima ou abaixo), o produto é positivo.

2- Os termos de covariância são multiplicados por \(w_{ij}\), o que faz com que as parcelas de covariância sejam ponderadas de acordo com a intensidade de sua interação espacial.

3- O denominador da equação de I representa a divisão pela variância.

O resultado da equação de I é análogo ao que se obtém para um coeficiente de correlação não espacial. Se houver correlação positiva dos dados, então a maioria dos polígonos vizinhos terá valores do mesmo lado da média e o índice será positivo, ou seja, \(I>0\), indicando correlação espacial direta.

Se os dados se correlacionam negativamente, então a maioria dos polígonos vizinhos terá valores de atributos em lados opostos da média e o índice será negativo, ou seja, \(I<0\), caso este seja de correlação espacial inversa.

Por fim, \(I=0\) indica ausência de correlação espacial.

2.1 Significância estatística

Uma vez calculado o índice global de Moran, é importante estabelecer sua validade estatística. Em outras palavras, será que os valores medidos representam correlação espacial significativa? Para estimar a significância do índice, será preciso associar a este uma distribuição estatística, sendo mais usual relacionar a estatística de teste à distribuição normal.

Outra possibilidade, sem pressupostos em relação à distribuição, e a abordagem mais comum é um teste de pseudo-significância. Neste caso,

• São geradas diferentes permutações dos valores de atributos associados às regiões;
• Cada permutação produz um novo arranjo espacial, onde os valores estão redistribuídos entre as áreas;
• Pode-se construir uma distribuição empírica de I e apenas um dos arranjos corresponde à situação observada;
• Se o valor do índice I medido originalmente corresponder a um “extremo” da distribuição simulada, então trata-se de valor com significância estatística.

Os comandos a seguir calculam o valor do índice global de Moran e realiza o teste de hipóteses para o índice I considerando os dados de Minas Gerais:

## Transformando objeto sf em sp
Dados_MG_sp <- sf:::as_Spatial(Dados_MG$geometry)
## Calculando a contiguidade queen
Queen_MG<-spdep::poly2nb(Dados_MG_sp, queen=T)
## 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_MG,listw,alternative="two.sided")
GlobalMoran

## 
##  Moran I test under randomisation
## 
## data:  taxa_MG  
## 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.0004487458

2.2 Atividade

Reproduza o exemplo anterior para os dados dos municípios do Rio de Janeiro.

2.3 Diagrama de Espalhamento de Moran

O diagrama de espalhamento de Moran é uma maneira adicional de visualizar a dependência espacial, sendo uma ferramenta visual útil para análise exploratória. Ele permite avaliar o quão semelhante é um valor observado com suas observações vizinhas.

Seu eixo horizontal é baseado nos valores do fenômeno de intersse e também é conhecido como eixo de resposta. O eixo Y vertical é baseado na média ponderada da observação correspondente no eixo X horizontal. O gráfico de dispersão de Moran fornece uma representação visual das associações espaciais na vizinhança em torno de cada observação.

No R esse diagrama é construído com a seguinte função:

## Plot do diagrama de espalhamento de Moran
moran.plot(taxa_MG,listw,label=Dados_MG$name_muni)

Para interpretar o gráfico, observe os quatro quadrantes:

• Quadrante 1: O quadrante superior direito contém os casos onde o valor do atributo de cada polígono e o valor médio do atributo nos polígonos vizinhos são maiores que a média global (ambos representados por alto): alto-alto (AA);

• Quadrante 2: O quadrante superior esquerdo contém os casos onde o *valor do atributo de cada polígono está abaixo da média global (representado por baixo), enquanto o valor médio do atributo em polígonos vizinhos está acima daquela média (representado por alto): baixo-alto (BA).

• Quadrante 3: O quadrante inferior esquerdo contém os casos onde o valor do atributo de cada polígono e o valor médio do atributo nos polígonos vizinhos são menores que a média global (ambos representados por baixo): baixo-baixo (BB);

• Quadrante 4: O quadrante inferior direito contém os casos onde o valor do atributo de cada polígono está acima da média global (representado por alto), enquanto o valor médio do atributo em polígonos vizinhos está abaixo daquela média (representado por baixo): alto-baixo (AB);

A maior concentração de pontos nos quadrantes 1 e 3 é uma indicação de autocorrelação espacial positiva \(I> 0\). Se os quadrantes 2 e 4 tivessem maior quantidade de pontos, a indicação seria de autocorrelação espacial negativa \(I < 0\). Caso os pontos estivessem igualmente distribuídos pelos quatro quadrantes, haveria uma indicação de ausência de autocorrelação \(I\approx0\).

Com base na teoria estatística de regressão linear, é correto afirmar que o índice \(I\) equivale ao coeficiente de regressão linear da reta de regressão sendo o coeficiente de correlação para o relacionamento entre os valores de atributos e os valores médios do atributo dos vizinhos.

2.4 Atividade

Reproduza o diagrama de espalhamento de Moran para os dados dos municípios do Rio de Janeiro.

3 Índice Moran Local (LISA)

Os indicadores globais de autocorrelação espacial, como o índice de Moran, fornecem um único valor como medida da associação espacial para todo o conjunto de dados, o que é útil na caracterização da região de estudo como um todo. Quando lidamos com grande número de áreas, é muito provável que ocorram diferentes regimes de associação espacial, onde a dependência espacial é ainda mais pronunciada.

Assim, muitas vezes é desejável examinar padrões em maior detalhe. Para tanto, é preciso utilizar indicadores de associação espacial que possam ser associados às diferentes localizações de uma variável distribuída espacialmente.

Os indicadores locais produzem um valor específico para cada área, permitindo assim a identificação de agrupamentos. O índice local de Moran pode ser expresso para cada área \(i\) a partir dos valores normalizados \(z_i^*\) do atributo \(z_i\) como:

\(\Large I_i=\frac{z_i^*\sum_{j=1}^nw_{ij}z_j^*}{\sum_{j=1}^n(z_j^*)^2},\)

em que \(z_i^*=\frac{z_i-\bar{z}}{\sigma_z}\) é a variável padronizada.

Da equação de \(I_i\), percebe-se que resultados positivos são obtidos onde ocorrem concentrações de valores baixos ou de valores altos do atributo, enquanto resultados negativos decorrem da proximidade entre valores baixos e altos na mesma área. Assim, o Índice de Moran Local dá uma indicação da homogeneidade e da diversidade dos dados.

3.1 Significância Estatística

Como fazer para determinar a significância estatística do resultado obtido para cada valor de \(I_i\)? O procedimento é similar ao adotado para o Índice de Moran Global, utilizando o método de permutações:

• Para cada área, uma vez calculado e fixado o respectivo Índice de Moran Local, permuta-se, aleatoriamente, o valor de atributo das demais áreas;
• Em cada permutação, calcula-se o valor correspondente da estatística local até que se obtenha uma pseudodistribuição, para a qual seja possível calcular os parâmetros de significância de \(I_i\)
• Um p-valor, ou seja, um valor de pseudossignificância, pode ser determinado a partir da posição do valor da estatística local \(I_i\), relativamente ao conjunto de valores obtidos para todas as permutações realizadas.

Os comandos a seguir auxiliam na construção e teste dos indicadores de Moran locais:

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

É recomendável gerar um mapa temático, indicando as regiões com correlação local significativamente diferente do resto dos dados. Essas áreas podem ser expressas como agrupamentos de estacionariedade, com dinâmica espacial própria e, por essa razão, merecem estudo mais detalhado. Os comandos à seguir constroem o mapa

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

Dinamicamente:

## 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")
  x <- girafe(ggobj = gg)
print(x)

3.2 Atividade

Faça o mapa com o LISA para o estado do Rio de Janeiro

4 Índice Brasileiro de Privação

É uma medida que informa níveis de privação material ou, de um modo mais geral, níveis de posição socioeconômica, em diferentes áreas geográficas do Brasil: setores censitários, municípios, estados, macrorregiões e Brasil como um todo.

O IBP foi calculado utilizando informações de renda, escolaridade e condições do domicílio. Com o IBP, é possível monitorar e avaliar as condições de privação sobre a saúde da população.

O IBP foi desenvolvido por pesquisadores do Centro de Integração de Dados e Conhecimentos para Saúde (Cidacs/Fiocruz Bahia) e da Universidade de Glasgow-Escócia.

Ao oferecer um nível particular de visualização de áreas de maior e menor privação, o IBP apoia gestores, servidores, profissionais da saúde e pesquisadores na tomada de decisões relativa ao planejamento e avaliação de políticas públicas, nos âmbitos da saúde, da educação, da habitação, saneamento básico, entre outras.

Os códigos à seguir representam o IBP para os municípios de Minas Gerais:

## Função utilizada para categorizar os pontos de corte
age.cat <- function(x, breaks,
                    sep = " to ", above.char = "+",min,max) {
  labs <- c(paste(c(min,breaks)[-(length(breaks)+1)],
                  breaks-0.01,
                  sep = sep),
            paste(breaks[length(breaks)], max, sep = " to "))
  cut(x, breaks = c(min,breaks,Inf),
      right = FALSE, labels = labs)
}

A construção do mapa coroplético considerando percentis como pontos de corte é feito a partir do seguinte código

## Leitura dos dados 
IBP_MG <- read.xlsx("IBP_MG.xlsx")  

## Pontos de corte
breaks <- round(quantile(IBP_MG$IBP,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
IBP_MG$Categoria <- age.cat(round(IBP_MG$IBP,2),breaks,min=min(round(IBP_MG$IBP,2)),max=max(round(IBP_MG$IBP,2)))
## Merge do arquivo contendo as taxas com o shapefile
Dados_MG <- merge(Dados_MG,IBP_MG,by="code_muni",all.x=T)

## Construção do mapa
Mapa_IBP_MG <- ggplot() +  geom_sf(data=Dados_MG,aes(fill=Categoria),colour="black") + 
            theme_bw() +
            scale_fill_manual(values = mycolors,name="IBP") +
            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_IBP_MG

4.1 Atividade

Reproduzir os dados de IBP para os municípios do Rio de Janeiro.

5 índice Moran Global Bivariado

O conceito de correlação espacial bivariada é complexo e muitas vezes mal interpretado. É normalmente considerada a correlação entre uma variável x e a média de outra variável y avaliada nos vizinhos. No entanto, isso não leva em conta a correlação inerente entre as duas variáveis. Sua fórmula é uma generalização da fórmula do índice de Moran Global e é dada por:

\(\Large I^B=\frac{\sum_{i=1}^n\sum_{j=1}^nw_{ij}x_i^*y_j^*}{\sum_{i=1}^n(x_i^*)^2},\)

em que \(x_i^*=\frac{x_i-\bar{x}}{\sigma_x}\) e \(y_i^*=\frac{y_i-\bar{y}}{\sigma_y}\).

No R, calcula-se o índice global de Moran bivariado com os seguintes comandos:

## Alternando formatos de matrizes de contiguidade
nb<- poly2nb(Dados_MG)
listw<- nb2listw(nb, style = "W")
## Variável X
x<-c(Dados_MG$taxa)
## Variável Y
y<-c(Dados_MG$IBP)
## Índice Moran Global Bivariado
moran_bv(x,y,nsim=999,listw)

## 
## DATA PERMUTATION
## 
## 
## Call:
## boot(data = xx, statistic = bvm_boot, R = nsim, sim = "permutation", 
##     listw = listw, parallel = parallel, ncpus = ncpus, cl = cl)
## 
## 
## Bootstrap Statistics :
##      original     bias    std. error
## t1* 0.1575831 -0.1578151  0.01580645

5.1 Atividade

Calcule o índice global de Moran bivariado para os municípios do estado do Rio de Janeiro.

6 Indice Moran Local Bivariado (LISA Bivariado)

Em essência, esse índice ele captura a relação entre o valor de uma variável no local i, e a média dos valores vizinhos para outra variável. Seu cálculo é dado por:

\(\Large I^B_i=x_i^*\sum_{j=1}^nw_{ij}y_j^*,\)

em que \(x_i^*=\frac{x_i-\bar{x}}{\sigma_x}\) e $y_i^*=

No R, tem-se os seguintes códigos para calcular os LISAs bivariados e testar suas significâncias

## Cálculo e teste dos LISAs bivariados.
bilisa <- local_bimoran(queen_w_MG, as.data.frame(cbind(x,y)))
## Colocando na shapefile as classificações 
Dados_MG$BILISA<-lisa_clusters(bilisa)
## Colocando na shape file os p-valores
Dados_MG$BILISA_PVALUE<-lisa_pvalues(bilisa)
## Substituindo números por caracteres
Dados_MG$BILISA<-gsub("0","Not Significant",Dados_MG$BILISA)
Dados_MG$BILISA<-gsub("1","High-High",Dados_MG$BILISA)
Dados_MG$BILISA<-gsub("2","Low-Low",Dados_MG$BILISA)
Dados_MG$BILISA<-gsub("3","Low-High",Dados_MG$BILISA)
Dados_MG$BILISA<-gsub("4","High-Low",Dados_MG$BILISA)

## Caleta de cores
mycolors <- c("red","pink","lightblue","blue","white")

## Mapa
ggplot() +  geom_sf(data=Dados_MG,aes(fill=BILISA),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")

Dinamicamente,

## O que será mostrado ao encostar em cada município
Dados_MG$tooltip <- c(paste0(Dados_MG$name_muni, "\n  ", Dados_MG$BILISA_PVALUE))
## Construção do mapa
gg <- ggplot(Dados_MG)  +
    geom_sf_interactive(aes(fill=BILISA, 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")
  x <- girafe(ggobj = gg)
print(x)

6.1 Atividade

Estude o comportamento do Lisa Bivariado para os municípios do estado do Rio de Janeiro.

7 Referências

Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) “Análise Espacial de Dados Geográficos”. Brasília, EMBRAPA, 2004 (ISBN: 85-7383-260-6).

https://cran.r-project.org/web/packages/geobr/index.html

https://www.healthdata.org/gbd

https://cidacs.bahia.fiocruz.br/ibp/indice/