1 Problematização

Em análise de dados espaciais, uma matriz W reflete as características do arranjo espacial, ou seja, o padrão de conectividade de certo fenômeno no espaço. Nesse sentido, o grau de interação característico do fenômeno tende a reduzir com a distância geográfica. O problema surge quando a mudança no arranjo é significativa o suficientemente ao ponto de comprometer a estabilidade paramétrica. Uma estatística muito popular entre estudos regionais é o coeficiente I de Moran. Sob condições de normalidade, o I de Moran mensura o grau de autocorrelação global subjacente a dado fenômeno. Assim, recorre-se ao método de simulação de Monte Carlo para testar a estabilidade do coeficiente I de Moran, extraído da taxa de variação da produtividade da lavoura temporária na Amazônia Legal, considerando diferentes arranjos espaciais. Para este fim, três matrizes W são avaliadas:1

2 Dados e metodologia

Como demonstrado adiante, o banco de dados é uma junção de arquivo shapefile extraído da Malha Municipal do IBGE com dados sobre o valor da produção e área colhida da lavoura temporária extraídos Censo Agropecuário 2017 do IBE.

A estatística I de Moran é inferida a partir dos resíduos de um modelo de regressão linear. Trata-se da distribuição da taxa de variação da produtividade agrícola \((\ln Y_i)\) condicionada à área colhida \((\ln A_i)\) . \[ \ln Y_i = \alpha + \beta \ln A_i + u_i. \] Em que, \(\alpha\) é o intercepto, \(\beta\) é o coeficiente angular e \(u_i\) são distúrbios.

Após a estimação, os resíduos do modelo de regressão fornecem as informações necessárias para estatística I de Moran.
\[ I = \frac{\sum_i \sum_j w_{ij} \hat{u}_i \hat{u}_j}{\hat{\sigma}^2 \sum_i \sum_j w_{ij}} \] Em que, \(\hat{\sigma}^2\) é a variância estimada e \(w_{ij} \in W\) são pesos espaciais que estabelecem a interação entre \(i\) e seu vizinho \(j\).

A hipótese nula do teste é ausência de autocorrelação espacial global.

3 Desenvolvimento

3.1 Carregando os pacotes

library(ggplot2)
library(sf)
library(spdep)
library(dplyr)

Os pacotes sf e ggplot2 funcionam conjuntamente na leitura e visualização dos metados. O pacote spdep fornece as funções necessárias à construção das matrizes W. A função moran.mc é usada na inferência da estatística I de Moran. A função left_join do pacote dplyr integra o banco de dados.

3.2 Organizando o banco de dados

O primeiro passo é ler o arquivo shapefile com a distribuição geográfica dos 772 municípios da Amazônia Legal.

bsf <- read_sf(dsn = "Mun_Amazonia_Legal_2022.shp", layer = "Mun_Amazonia_Legal_2022")
ggplot() +
  geom_sf(data=bsf)+
  theme_minimal()

O segundo passo é intergrar os dados do Censo Agropecuário 2017 ao arquivo shapefile

## Lendo dados do Censo Agropecuário 2017
b0 <- read.delim("PROD_LAV_TEMP.csv", header = TRUE, sep = ",", fill = TRUE)
b0$CD_MUN <- factor(b0$CODMUN)
## Juntando as bases
bsf <- left_join(bsf,b0,by="CD_MUN")
## Estimando a produtividade da lavoura temporária
resml <- lm(log(bsf$VPROD) ~ log(bsf$AREA))
## Extraíndo os resíduos do modelo de regressão
bsf$PRM_u <- resml$residuals

Estimando o modelo de regressão para a taxa de produtividade agrícola

summary(resml)
## 
## Call:
## lm(formula = log(bsf$VPROD) ~ log(bsf$AREA))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0067 -0.3802 -0.0603  0.3708  2.0710 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.11567    0.08172   25.89   <2e-16 ***
## log(bsf$AREA)  0.90423    0.01029   87.84   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6117 on 770 degrees of freedom
## Multiple R-squared:  0.9093, Adjusted R-squared:  0.9092 
## F-statistic:  7717 on 1 and 770 DF,  p-value: < 2.2e-16

3.3 Distribuição espacial da produtividade agrícola

A dispersão dos resíduos do modelo de regressão oferece as primeiras evidências acerca da distribuição espacial da produtividade agrícola na Amazômia Legal. Notadamente, produtores com desempenho acima da média estão localizados é áreas mais isoladas, entre municípios do Amazonas, Acre e Roraima. Produtores com desempenho abaixo da média estão localizados no Maranhão, Sudeste do Pará e Tocantins, principalmente.

ggplot() +
  geom_sf(data = bsf, aes(fill=PRM_u), color="gray50", size = .2) +
  scale_fill_distiller(palette = "RdBu", direction = -1, na.value = "grey60")+
  labs(title = "Dispersão produtiva da lavoura temporária 2017", x=NULL, y=NULL, caption = NULL, fill = "Resíduos") +
  theme_bw()

3.4 Matrizes W

Observe que a matriz de contiguidade queen wq é produzida a partir de uma lista (poly2nb) com os vizinhos de primeira ordem, extraída do polígono bsf. Com base nos centróides determinados pelas coordenadas geográficas, a matriz wknn é identificada com k=2 vizinhos mais próximos. A matriz wdc é identificada tomando a distância de corte equivalente a 2.57 graus (dmax). Sumariamente, 0.74% dos pesos da matriz wq são diferentes de zero. Esse percentual é de 0.26% na matriz wknn e de 8.9% em wdc.

## Contiguidade de primeira ordem normalizada na linha
list.nb <- poly2nb(bsf)
wq <- nb2listw(list.nb, style = "W")
## Matriz k-vizinhos mais próximos 
coords <- st_coordinates(st_centroid(bsf))
knn <- knn2nb(knearneigh(coords, k=2)) # Objeto "nb"
wknn <- nb2listw(knn,zero.policy = TRUE)
## Matriz baseada na distância de corte.
dmax <- max(unlist(nbdists(knn, coords)))
dnn <- dnearneigh(coords, 0, dmax, longlat = NULL)
wdc <- nb2listw(dnn)

3.5 Estatística I de Moran

A função moran.test retorna o coeficiente I de Moral global, sua expectância e variância, bem como oferece um teste de significância baseado na normalidade. A função moran.test retorna apenas o coeficiente I de Moran, com significância avaliada por meio de simulação de Monte Carlo – é usual proceder 99 permutações, que corresponde a um p-value igual a 0.01.

Embora com valores ligeiramente diferentes, os testes mostram que as três matrizes W geram coeficientes positivos e significantes, a 1% de probabilidade de erro. A matriz k=2 vizinhos gera o coeficiente com maior valor absoluto, oferecendo controle eficiente para a autocorrelação espacial. Não obstante, a matriz de contiguidade de primeira ordem, também ofercere controle eficiente, levando-se em conta a magnitude de sua variância. A matriz W baseada na distância de corte oferece gera um I de Moram com menor valor absoluto.

Por fim, mesmo relaxando a hipótese de normalidade, o teste de robustez baseado em simulação de Monte Carlo confirma a consistência do coeficiente I de Moral.

n <- 99
## Teste I de Moran com matriz W contiguidade de primeira ordem
moran.test(x = bsf$PRM_u, listw = wq)
## 
##  Moran I test under randomisation
## 
## data:  bsf$PRM_u  
## weights: wq    
## 
## Moran I statistic standard deviate = 23.308, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.5115296460     -0.0012970169      0.0004841097
moran.mc(x = bsf$PRM_u, listw = wq, nsim = n)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  bsf$PRM_u 
## weights: wq  
## number of simulations + 1: 100 
## 
## statistic = 0.51153, observed rank = 100, p-value = 0.01
## alternative hypothesis: greater
## Teste I de Moran com matriz W k-vizinhos mais próximos 
moran.test(x = bsf$PRM_u, listw = wknn)
## 
##  Moran I test under randomisation
## 
## data:  bsf$PRM_u  
## weights: wknn    
## 
## Moran I statistic standard deviate = 17.721, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.581904696      -0.001297017       0.001083052
moran.mc(x = bsf$PRM_u, listw = wknn, nsim = n)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  bsf$PRM_u 
## weights: wknn  
## number of simulations + 1: 100 
## 
## statistic = 0.5819, observed rank = 100, p-value = 0.01
## alternative hypothesis: greater
## Teste I de Moran com matriz W distância de corte
moran.test(x = bsf$PRM_u, listw = wdc)
## 
##  Moran I test under randomisation
## 
## data:  bsf$PRM_u  
## weights: wdc    
## 
## Moran I statistic standard deviate = 46.098, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      3.893242e-01     -1.297017e-03      7.180283e-05
moran.mc(x = bsf$PRM_u, listw = wdc, nsim = n)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  bsf$PRM_u 
## weights: wdc  
## number of simulations + 1: 100 
## 
## statistic = 0.38932, observed rank = 100, p-value = 0.01
## alternative hypothesis: greater

  1. KOPCZEWSKA, K. Applied spatial statistics and econometrics: data analysis in R. London and New York: Routledge, 2021.↩︎