O que acontece num município depende do vizinho?

Antes de tudo: que lugar é esse

MATOPIBA é a junção das siglas de quatro estados: MAranhão, TOcantins, PIauí e BAhia. Não é região do IBGE nem unidade administrativa. É um recorte criado pelo Decreto 8.447, de 2015, para orientar política agrícola federal numa faixa específica de cerrado.

A região é conhecida como a última fronteira agrícola do cerrado. Até os anos 1980 era pecuária extensiva e cerrado nativo; o solo é ácido e pobre, e a lavoura só se tornou viável com correção por calcário e cultivares adaptadas ao clima quente. A partir dos anos 1990 a soja avançou sobre os chapadões, primeiro no oeste da Bahia, depois no sul do Maranhão, no sudoeste do Piauí e no leste do Tocantins.

Hoje são 5.390.709 hectares de soja plantada nos 997 municípios, concentrados de um jeito bem desigual: 763 deles plantam zero, e um só, São Desidério, na Bahia, planta 493.639 hectares. Essa desigualdade é o assunto da aula.

Como é a tabela que você vai usar

É esta a tabela que você abre no passo 1. Uma linha por município, e a coluna soja_ha é a variável da aula inteira. Estas são as primeiras linhas, na ordem em que o QGIS abre, por código do IBGE:

E estas são as cinco maiores áreas plantadas do recorte, que vão aparecer em todos os resultados da aula:

Quatro colunas de identificação e uma de medida. Repare em CD_MUN: são sete dígitos de texto, não um número, e é essa coluna que liga a malha à tabela do SIDRA. Repare também em Aguiarnópolis, com zero: é a situação de três quartos do recorte, e ela domina quase tudo que vem depois.

Na aula 1 você viu a primeira lei de Tobler: tudo se relaciona com tudo, mas coisas próximas se relacionam mais. Ela descreve o que se vê em quase todo mapa temático brasileiro. O que ela não diz é por quê, e existem duas razões possíveis que produzem manchas parecidas e exigem tratamentos opostos.

A matriz de vizinhança

Para dizer que um valor depende do vizinho, é preciso primeiro dizer quem é vizinho. Isso vira uma matriz quadrada W, de tamanho n por n, em que a célula w_ij vale alguma coisa se j é vizinho de i e zero se não é. A diagonal é sempre zero: ninguém é vizinho de si mesmo.

Existem várias maneiras de preencher essa matriz, e elas não concordam entre si.

Depois de decidir quem é vizinho, falta decidir o peso. A prática dominante é a padronização por linha: cada linha da matriz é dividida pela sua soma, de modo que os pesos de cada município somem 1. Um município com 12 vizinhos dá peso 1/12 a cada um; um com 2 vizinhos dá 1/2. O efeito é que a matriz passa a calcular médias em vez de somas, e nenhum município influencia mais que outro só por ter muita divisa.

A média móvel espacial

Com W padronizada por linha, multiplicar a matriz pelo vetor de valores produz, para cada município, a média dos valores dos seus vizinhos. Esse novo vetor se chama defasagem espacial, ou spatial lag, e é a peça central de tudo o que vem depois.

lag(i) = soma de w(i,j) · x(j)   -- com pesos somando 1, é a média dos vizinhos de i

É a mesma ideia da média móvel de uma série temporal, com uma diferença: no tempo, o passado é único e vem antes. No espaço não existe antes, e cada unidade tem vários vizinhos, todos simultâneos. Por isso a estatística espacial precisou de aparato próprio em vez de reaproveitar a de séries.

No QGIS não existe botão de matriz de pesos, e isso é uma vantagem pedagógica. Você vai construir a defasagem com uma junção espacial por resumo, que faz exatamente a conta acima, e vai ver o número aparecer coluna por coluna.

O índice de Moran

Patrick Moran propôs em 1950 um índice que mede se um valor alto tende a aparecer perto de outro valor alto. A forma que interessa aqui, com W padronizada por linha, é uma razão entre duas somas.

I = soma de (x(i) − x̄)·(lag(i) − x̄)
    ────────────────────────────────
        soma de (x(i) − x̄)²

numerador   : quanto o desvio de cada município acompanha o desvio dos vizinhos
denominador : a variação total da variável
resultado   : uma correlação entre o que você tem e o que os seus vizinhos têm

O índice varia aproximadamente entre −1 e 1. Perto de 1, valores parecidos se agrupam. Perto de −1, valores opostos se alternam, como um tabuleiro de xadrez. Perto de zero, a distribuição no espaço é indiferente ao valor.

Essa conta é a que o GeoDa faz no passo 4, sobre a soja do MATOPIBA. O resultado é I = 0,3256, e é esse número que as próximas duas seções vão interrogar: primeiro se ele poderia ter saído por acaso, depois o que ele esconde.

Inferência por permutação

Ter o número não basta. 0,3256 é muito ou é pouco? Sem uma régua, não dá para dizer: talvez qualquer distribuição de soja por 997 municípios produza algo parecido, e o que parece padrão seja só o jeito como os números caíram.

A régua tradicional seria uma tabela estatística, que exige supor que o dado segue uma curva normal. A área plantada de soja não segue: três quartos dos municípios estão em zero e um deles planta quase meio milhão de hectares. A saída é construir a régua com os próprios dados, embaralhando.

O procedimento é direto. Mantenha a geometria e a matriz W paradas. Sorteie uma nova atribuição dos 997 valores de soja aos 997 municípios, ao acaso. Calcule I. Repita 999 vezes. Você obtém uma distribuição de referência construída com os seus próprios dados, sem hipótese sobre a forma da distribuição.

Vale ver isso acontecendo. Na animação abaixo o exemplo é menor, 25 unidades em vez de 997, para caber na tela, mas o procedimento é exatamente o mesmo.

O diagrama de espalhamento de Moran

Luc Anselin propôs em 1996 uma leitura gráfica do índice que é mais informativa que o número. No eixo horizontal, o valor padronizado de cada município. No eixo vertical, a média padronizada dos seus vizinhos. Cada município vira um ponto, e a inclinação da reta ajustada é exatamente o índice de Moran.

Os quatro quadrantes têm nome e significado.

LISA e o problema das múltiplas comparações

O índice global responde uma pergunta só: existe padrão. Não diz onde. Anselin resolveu isso decompondo o índice em uma parcela por unidade, e chamou o resultado de indicador local de associação espacial, LISA na sigla em inglês.

Cada município ganha o seu próprio I local, o seu quadrante e o seu pseudo p-valor, calculado por permutação condicional: os valores dos outros 996 municípios são embaralhados, o do município em questão fica parado, e se mede quantas vezes o I local sorteado supera o observado.

Na prática, o mapa LISA é usado como ferramenta exploratória, não como bateria de testes confirmatórios. A leitura correta é: aqui há indício de agrupamento, vá olhar. A leitura incorreta é tratar cada município significativo como uma descoberta estatística independente.

Processo real ou artefato da unidade de agregação

Falta o problema mais desconfortável, e ele vem da aula 4. O índice de Moran é calculado sobre unidades administrativas que ninguém desenhou pensando em soja. Município é uma unidade política, de tamanho arbitrário, criada por lei estadual em datas diferentes. Se a unidade muda, o índice muda.

Há um segundo efeito, do mesmo tipo. A área plantada em hectares depende do tamanho do município: um município grande tende a ter mais soja simplesmente por ser grande. Trocando os hectares pela fração da área municipal ocupada por soja, o índice desta aula sobe de 0,3256 para 0,4457. Os dois são positivos e significativos, e medem coisas diferentes: o primeiro fala de onde está a produção, o segundo de onde a paisagem virou lavoura. No topo do segundo está Luís Eduardo Magalhães, com 53,2% do município plantado.

Seis passos no GeoDa. Nenhum deles tem conta para fazer na mão: o trabalho é escolher o que pedir e entender o que voltou.

A segunda atividade prática supervisionada é lançada hoje e usa exatamente este roteiro em outro território. Cada aluno recebe um recorte do país sorteado pelo trimestre de nascimento, com a variável já preparada, e a pergunta de maior peso é sobre os outliers espaciais: por que um município destoa da vizinhança. O enunciado está em APS 2 · autocorrelação espacial em quatro recortes do Brasil.

O QGIS não tem autocorrelação espacial no núcleo. Dá para montar a conta na mão, e quem quiser ver de onde sai cada pedaço do índice encontra aqui o caminho completo. Não é necessário para a aula nem para a prova.

("soja_ha" - mean("soja_ha")) * ("soja_ha_mean" - mean("soja_ha"))

("soja_ha" - mean("soja_ha")) ^ 2

sum("dxdlag") / sum("dx2")   -- o mesmo valor em todas as 997 linhas

CASE
  WHEN "soja_ha" >= mean("soja_ha") AND "soja_ha_mean" >= mean("soja_ha") THEN 'Alto-Alto'
  WHEN "soja_ha" <  mean("soja_ha") AND "soja_ha_mean" <  mean("soja_ha") THEN 'Baixo-Baixo'
  WHEN "soja_ha" >= mean("soja_ha") AND "soja_ha_mean" <  mean("soja_ha") THEN 'Alto-Baixo'
  ELSE 'Baixo-Alto'
END