Objetivo: explorar padrões de autocorrelação espacial no desmatamento
dos municípios da Amazônia Legal usando R — do carregamento dos dados
georreferenciados até o mapa LISA de clusters locais.

Dados: PRODES/INPE 2022 · Geometrias: geodata-br (GitHub)
Pacotes principais: sf, spdep, ggplot2, dplyr

1 Configuração

1.1 Pacotes necessários

# Instalar se necessário:
# install.packages(c("sf","spdep","ggplot2","dplyr","tidyr",
#                    "RColorBrewer","patchwork","scales"))

library(sf)
library(spdep)
library(ggplot2)
library(dplyr)
library(tidyr)
library(RColorBrewer)
library(patchwork)
library(scales)

2 Bloco 1 — Dados espaciais no R

2.1 Download das geometrias

O GeoJSON com todos os 5.564 municípios brasileiros está hospedado no repositório público geodata-br no GitHub. O script verifica se o arquivo já existe antes de baixar novamente — boa prática para não sobrecarregar a conexão em sala de aula.

url_municipios <- paste0(
  "https://raw.githubusercontent.com/",
  "tbrugz/geodata-br/master/geojson/geojs-100-mun.json"
)

dir.create("dados", showWarnings = FALSE)
destino <- file.path("dados", "municipios_br.geojson")

if (!file.exists(destino)) {
  message("Baixando geometrias dos municípios...")
  download.file(url_municipios, destfile = destino, mode = "wb", quiet = TRUE)
} else {
  message("Arquivo já existe, pulando download.")
}

mun_br <- st_read(destino, quiet = TRUE)

2.2 Explorando o objeto sf

Um objeto sf é um data.frame com uma coluna especial geometry que guarda as formas geométricas. Todas as funções do tidyverse funcionam normalmente sobre ele.

class(mun_br)          # "sf" e "data.frame"
[1] "sf"         "data.frame"
dim(mun_br)            # linhas = municípios, colunas = atributos + geometry
[1] 5564    4
names(mun_br)          # nomes das colunas
[1] "id"          "name"        "description" "geometry"   
# Sistema de Referência de Coordenadas
st_crs(mun_br)$input   # WGS 84 (EPSG:4326)
[1] "WGS 84"
# Extensão espacial
st_bbox(mun_br)
      xmin       ymin       xmax       ymax 
-73.990944 -33.751583 -32.392190   5.272156

Nota sobre CRS: O referencial oficial para o Brasil é o SIRGAS 2000 (EPSG:4674). O WGS 84 (EPSG:4326) usado aqui é praticamente equivalente para análises em escala municipal.

2.4 Dados PRODES (desmatamento)

O script tenta baixar o CSV consolidado do repositório TerraBrasilis/INPE. Se a URL não estiver disponível, gera dados sintéticos realistas com base nos valores publicados no Relatório PRODES 2022 — apenas para fins didáticos.

url_prodes <- paste0(
  "https://raw.githubusercontent.com/",
  "terrabrasilis/prodes-data-plots/master/",
  "data/prodes_annual_rates_amazon_municipalities.csv"
)

destino_prodes <- file.path("dados", "prodes_municipios.csv")

if (!file.exists(destino_prodes)) {
  resultado <- tryCatch({
    download.file(url_prodes, destfile = destino_prodes, mode = "wb", quiet = TRUE)
    "ok"
  }, error = function(e) "erro")

  if (resultado == "erro") {
    # ------------------------------------------------------------------
    # FALLBACK: dados sintéticos baseados no PRODES 2022 (INPE)
    # Valores por UF derivados do relatório oficial (km²/ano)
    # ------------------------------------------------------------------
    set.seed(2022)

    uf_total <- c(
      "11" = 1465, "12" = 459,  "13" = 2763, "14" = 543,
      "15" = 3287, "16" = 122,  "17" = 179,  "21" = 812, "51" = 1788
    )

    prodes_list <- lapply(names(uf_total), function(uf) {
      muns_uf <- amz |> filter(uf_code == uf) |> st_drop_geometry()
      n  <- nrow(muns_uf)
      mu <- uf_total[uf] / n
      vals <- rlnorm(n, meanlog = log(max(mu, 0.5)), sdlog = 1.2)
      data.frame(cod_ibge = muns_uf$cod_ibge, desmat_km2 = round(vals, 2))
    })

    prodes_df <- do.call(rbind, prodes_list)
    write.csv(prodes_df, destino_prodes, row.names = FALSE)
  }
}

prodes_df <- read.csv(destino_prodes) |>
  mutate(cod_ibge = as.character(cod_ibge))   # garantir character para join

summary(prodes_df$desmat_km2)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
   0.060    2.505    7.640   29.791   26.835 1424.880

2.5 Junção e mapa exploratório

amz <- amz |>
  left_join(prodes_df, by = "cod_ibge") |>
  mutate(
    desmat_km2 = tidyr::replace_na(desmat_km2, 0),
    log_desmat = log1p(desmat_km2)          # log(1+x): trata assimetria
  )

message("Municípios sem dado PRODES: ", sum(is.na(amz$desmat_km2)))
tema_mapa <- theme_minimal(base_size = 11) +
  theme(
    panel.grid       = element_blank(),
    axis.text        = element_blank(),
    axis.ticks       = element_blank(),
    legend.position  = "bottom",
    legend.key.width = unit(2, "cm"),
    plot.title       = element_text(face = "bold", size = 13),
    plot.subtitle    = element_text(color = "grey40", size = 10)
  )

p_mapa_bruto <- ggplot(amz) +
  geom_sf(aes(fill = desmat_km2 + 1), color = NA) +
  scale_fill_gradientn(
    colours = brewer.pal(9, "YlOrRd"),
    trans   = "log10",
    name    = "km² (escala log)",
    labels  = label_comma(accuracy = 1),
    breaks  = c(1, 10, 100, 500)
  ) +
  labs(
    title    = "Desmatamento acumulado — Amazônia Legal",
    subtitle = "PRODES/INPE 2022  |  municípios",
    caption  = "Fonte: INPE/TerraBrasilis"
  ) +
  tema_mapa

p_mapa_bruto
Desmatamento acumulado nos municípios da Amazônia Legal. Escala logarítmica necessária dado o forte assimetria da distribuição (mediana: 7,6 km²; máximo: 1.425 km²).

Desmatamento acumulado nos municípios da Amazônia Legal. Escala logarítmica necessária dado o forte assimetria da distribuição (mediana: 7,6 km²; máximo: 1.425 km²).

Observação: a escala logarítmica é essencial — sem ela, os poucos municípios com desmatamento extremo (Codajás/AM: 1.425 km²) comprimem toda a variação dos demais para uma única cor clara.

3 Bloco 2 — Matriz de vizinhança

A matriz de pesos espaciais \(W\) codifica quem é vizinho de quem. É o ingrediente fundamental do I de Moran e do LISA.

3.1 Critérios de contiguidade

Critério Definição Vizinhos médios
Rainha Compartilha ponto ou aresta 5,54
Torre Apenas aresta compartilhada 5,38 
nb_queen <- poly2nb(amz, queen = TRUE)   # critério Rainha
nb_rook  <- poly2nb(amz, queen = FALSE)  # critério Torre

summary(nb_queen)
Neighbour list object:
Number of regions: 807 
Number of nonzero links: 4470 
Percentage nonzero weights: 0.6863734 
Average number of links: 5.539033 
2 disjoint connected subgraphs
Link number distribution:

  1   2   3   4   5   6   7   8   9  10  11  12  13 
  3  29  65 156 185 158  94  48  37  18   8   4   2 
3 least connected regions:
605 638 674 with 1 link
2 most connected regions:
508 794 with 13 links

Sobre o aviso “2 disjoint connected subgraphs”: indica municípios insulares (ex.: arquipélagos fluviais do PA e AP) sem contiguidade direta com o continente. Como não há municípios com card = 0, o aviso é informativo e não afeta os cálculos.

# Verificar e corrigir municípios sem vizinhos (card == 0)
n_ilhas <- sum(card(nb_queen) == 0)
message("Municípios sem vizinho: ", n_ilhas)

if (n_ilhas > 0) {
  coords   <- st_coordinates(st_centroid(st_geometry(amz)))
  nb_knn1  <- knearneigh(coords, k = 1) |> knn2nb()
  nb_queen <- union.nb(nb_queen, nb_knn1)
}

3.2 Lista de pesos (listw)

# style = "W": linha-normalizada (cada linha soma 1)
# Interpretação: média ponderada dos vizinhos
lw_queen <- nb2listw(nb_queen, style = "W", zero.policy = TRUE)

# Verificação: pesos do primeiro município somam 1
cat("Pesos do município 1:", round(lw_queen$weights[[1]], 4))
Pesos do município 1: 0.1429 0.1429 0.1429 0.1429 0.1429 0.1429 0.1429
cat("\nSoma:", sum(lw_queen$weights[[1]]))

Soma: 1

3.3 Visualização da rede

coords <- st_coordinates(st_centroid(st_geometry(amz)))

plot(st_geometry(amz),
     border = "grey70", col = "grey95", lwd = 0.2,
     main   = "Rede de contiguidade — critério Rainha")

plot(nb_queen, coords,
     add    = TRUE,
     col    = adjustcolor("steelblue", alpha.f = 0.4),
     lwd    = 0.5,
     points = FALSE)
Grafo de contiguidade (critério Rainha). Cada aresta conecta os centróides de dois municípios vizinhos. A densidade é maior no interior do AM e PA (municípios menores).

Grafo de contiguidade (critério Rainha). Cada aresta conecta os centróides de dois municípios vizinhos. A densidade é maior no interior do AM e PA (municípios menores).

4 Bloco 3 — I de Moran Global

4.1 Conceito

O I de Moran mede a autocorrelação espacial global:

\[I = \frac{n}{S_0} \cdot \frac{\mathbf{z}^\top W \mathbf{z}}{\mathbf{z}^\top \mathbf{z}}\]

onde \(\mathbf{z}\) é a variável padronizada, \(W\) a matriz de pesos e \(S_0 = \sum_i \sum_j w_{ij}\).

Interpretação:

  • \(I > E(I)\) → autocorrelação positiva (clustering de valores similares)
  • \(I \approx E(I)\) → distribuição aleatória no espaço
  • \(I < E(I)\) → autocorrelação negativa (padrão em tabuleiro de xadrez)

O valor esperado sob \(H_0\) é \(E(I) = -1/(n-1) \approx 0\).

Atenção: \(I\) não está restrito ao intervalo \([-1, 1]\). Os limites teóricos dependem da estrutura da matriz \(W\).

4.2 Teste analítico

moran_test <- moran.test(
  x           = amz$desmat_km2,
  listw       = lw_queen,
  alternative = "greater",     # H1: clustering positivo
  zero.policy = TRUE
)

print(moran_test)

    Moran I test under randomisation

data:  amz$desmat_km2  
weights: lw_queen    

Moran I statistic standard deviate = 13.578, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.2667728314     -0.0012406948      0.0003895961 
I_obs   <- moran_test$estimate["Moran I statistic"]
I_esp   <- moran_test$estimate["Expectation"]
I_var   <- moran_test$estimate["Variance"]
z_score <- (I_obs - I_esp) / sqrt(I_var)
p_valor <- moran_test$p.value

cat("I observado :", round(I_obs,   4), "\n")
I observado : 0.2668 
cat("I esperado  :", round(I_esp,   6), "\n")
I esperado  : -0.001241 
cat("z-score     :", round(z_score, 3), "\n")
z-score     : 13.578 
cat("p-valor     :", format(p_valor, scientific = TRUE), "\n")
p-valor     : 2.688484e-42

4.3 Teste de permutação (Monte Carlo)

A inferência analítica assume normalidade. Como o desmatamento segue uma distribuição assimétrica, o teste de permutação é mais apropriado: embaralha os valores entre os municípios 999 vezes, recalcula \(I\) e reporta a proporção de \(I_{sim} \geq I_{obs}\).

set.seed(42)
moran_mc <- moran.mc(
  x           = amz$desmat_km2,
  listw       = lw_queen,
  nsim        = 999,
  alternative = "greater",
  zero.policy = TRUE
)

print(moran_mc)

    Monte-Carlo simulation of Moran I

data:  amz$desmat_km2 
weights: lw_queen  
number of simulations + 1: 1000 

statistic = 0.26677, observed rank = 1000, p-value = 0.001
alternative hypothesis: greater
plot(moran_mc,
     main = "Distribuição de referência — permutação de Moran (999 sim.)",
     xlab = "I de Moran simulado",
     col  = "steelblue")
abline(v = moran_mc$statistic, col = "red", lwd = 2, lty = 2)
Distribuição de referência do I de Moran sob a hipótese nula de aleatoriedade espacial (999 permutações). A linha vermelha tracejada indica o I observado (0,267), bem além do intervalo simulado.

Distribuição de referência do I de Moran sob a hipótese nula de aleatoriedade espacial (999 permutações). A linha vermelha tracejada indica o I observado (0,267), bem além do intervalo simulado.

4.4 Diagrama de dispersão de Moran

O diagrama de Moran plota o valor padronizado de cada município (\(z_i\)) no eixo x contra o seu lag espacial (\(Wz_i\) = média ponderada dos vizinhos) no eixo y. A inclinação da reta de regressão é numericamente igual ao I.

mp <- moran.plot(
  x           = amz$desmat_km2,
  listw       = lw_queen,
  zero.policy = TRUE,
  main        = "Diagrama de dispersão de Moran — desmatamento",
  xlab        = "Desmatamento padronizado (z)",
  ylab        = "Lag espacial padronizado (Wz)",
  labels      = FALSE
)
Diagrama de dispersão de Moran para o desmatamento. Pontos fora do intervalo são os municípios 'influentes' — aqueles com maior leverage no cálculo do I global.

Diagrama de dispersão de Moran para o desmatamento. Pontos fora do intervalo são os municípios ‘influentes’ — aqueles com maior leverage no cálculo do I global. 

# Municípios com maior leverage (hat > mean)
influentes <- amz[mp$is_inf, ] |>
  select(nome_municipio, uf_sigla, desmat_km2) |>
  st_drop_geometry() |>
  arrange(desc(desmat_km2))

head(influentes, 10)

4.5 Moran com variável log-transformada

amz <- amz |>
  mutate(log_desmat = log1p(desmat_km2))   # log(1 + x)

moran_log <- moran.test(
  x           = amz$log_desmat,
  listw       = lw_queen,
  alternative = "greater",
  zero.policy = TRUE
)

cat("I de Moran (log_desmat):", round(moran_log$estimate[1], 4), "\n")
I de Moran (log_desmat): 0.4387 
cat("p-valor               :", format(moran_log$p.value, scientific = TRUE), "\n")
p-valor               : 5.259296e-90

Por que o I aumentou? Os outliers extremos (Codajás, Caapiranga) inflam a variância total e diluem o sinal espacial. A transformação log reduz esse efeito, tornando o padrão de clustering mais nítido (\(I = 0{,}439\)).

5 Bloco 4 — LISA (I de Moran Local)

5.1 Conceito

O LISA decompõe o I global em contribuições locais para cada município:

\[I_i = z_i \sum_{j} w_{ij}\, z_j\]

Cada município é então classificado em um quadrante com base no sinal de \(z_i\) (valor próprio) e de \(\sum_j w_{ij} z_j\) (lag espacial = média dos vizinhos):

Quadrante \(z_i\) lag Interpretação
HH \(> 0\) \(> 0\) Alto cercado de alto — hotspot
LL \(< 0\) \(< 0\) Baixo cercado de baixo — coldspot
HL \(> 0\) \(< 0\) Outlier: alto cercado de baixo
LH \(< 0\) \(> 0\) Outlier: baixo cercado de alto

5.2 Cálculo com permutação

Nota de versão: a partir do spdep ≥ 1.3, as permutações locais usam localmoran_perm() separada de localmoran(). Misturar as duas é o erro mais comum nesta versão do pacote.

set.seed(42)

# Passo 1: I local com p-valor assintótico
lisa <- localmoran(
  x           = amz$log_desmat,
  listw       = lw_queen,
  alternative = "two.sided",
  zero.policy = TRUE
)

# Passo 2: p-valor por permutação (999 simulações)
lisa_perm <- localmoran_perm(
  x           = amz$log_desmat,
  listw       = lw_queen,
  nsim        = 999,
  alternative = "two.sided",
  zero.policy = TRUE
)

# Identificar coluna do p-valor simulado (à prova de versão do spdep)
col_p_sim <- grep("Sim", colnames(lisa_perm), value = TRUE)[1]
cat("Coluna usada para p-valor:", col_p_sim, "\n")
Coluna usada para p-valor: Pr(z != E(Ii)) Sim 
amz <- amz |>
  mutate(
    lisa_I = lisa[, "Ii"],
    lisa_Z = lisa[, "Z.Ii"],
    lisa_p = lisa_perm[, col_p_sim]
  )

5.3 Classificação dos quadrantes

z     <- scale(amz$log_desmat)[, 1]
z_lag <- lag.listw(lw_queen, z)

amz <- amz |>
  mutate(
    z_padrao  = z,
    z_lag     = z_lag,
    quadrante = case_when(
      lisa_p >= 0.05               ~ "Não sig.",
      z_padrao > 0 & z_lag > 0    ~ "HH",
      z_padrao < 0 & z_lag < 0    ~ "LL",
      z_padrao > 0 & z_lag < 0    ~ "HL",
      z_padrao < 0 & z_lag > 0    ~ "LH",
      TRUE                         ~ "Não sig."
    ),
    quadrante = factor(quadrante,
                       levels = c("HH","LL","HL","LH","Não sig."))
  )

table(amz$quadrante)

      HH       LL       HL       LH Não sig. 
      79      102       15       20      591

5.4 Mapa LISA cluster map

cores_lisa <- c(
  "HH"       = "#d7191c",
  "LL"       = "#2c7bb6",
  "HL"       = "#fdae61",
  "LH"       = "#abd9e9",
  "Não sig." = "#ebebeb"
)

p_lisa <- ggplot(amz) +
  geom_sf(aes(fill = quadrante), color = NA) +
  scale_fill_manual(
    values = cores_lisa,
    name   = "Cluster LISA\n(p < 0,05)",
    drop   = FALSE
  ) +
  labs(
    title    = "Mapa LISA — Desmatamento (log)",
    subtitle = "I de Moran local · 999 permutações · Amazônia Legal",
    caption  = "Fonte: PRODES/INPE 2022  •  Matriz de contiguidade Rainha"
  ) +
  tema_mapa +
  theme(legend.position = "right")

p_lisa
Mapa LISA de clusters de desmatamento. Vermelho (HH) = municípios com alto desmatamento cercados por vizinhos também com alto desmatamento — corresponde ao 'arco do desmatamento'. Azul (LL) = florestas remotas com baixa pressão antrópica.

Mapa LISA de clusters de desmatamento. Vermelho (HH) = municípios com alto desmatamento cercados por vizinhos também com alto desmatamento — corresponde ao ‘arco do desmatamento’. Azul (LL) = florestas remotas com baixa pressão antrópica.

5.5 Mapa de significância local

amz <- amz |>
  mutate(
    sig_cat = cut(
      lisa_p,
      breaks = c(0, 0.01, 0.05, 0.10, 1),
      labels = c("p < 0,01", "p < 0,05", "p < 0,10", "Não sig."),
      right = TRUE, include.lowest = TRUE
    )
  )

p_sig <- ggplot(amz) +
  geom_sf(aes(fill = sig_cat), color = NA) +
  scale_fill_manual(
    values = c("p < 0,01" = "#b2182b",
               "p < 0,05" = "#ef8a62",
               "p < 0,10" = "#fddbc7",
               "Não sig." = "#ebebeb"),
    name = "Significância\n(permutação)"
  ) +
  labs(title    = "Significância local — I de Moran",
       subtitle = "Amazônia Legal") +
  tema_mapa +
  theme(legend.position = "right")

p_sig
Distribuição espacial dos p-valores locais (permutação). Municípios em vermelho escuro apresentam a evidência mais forte de não-aleatoriedade espacial.

Distribuição espacial dos p-valores locais (permutação). Municípios em vermelho escuro apresentam a evidência mais forte de não-aleatoriedade espacial.

5.6 Painel comparativo

p_mapa_bruto + p_lisa +
  plot_annotation(
    title   = "Desmatamento na Amazônia Legal — análise exploratória espacial",
    caption = "PRODES/INPE 2022  •  spdep + ggplot2",
    theme   = theme(plot.title = element_text(face = "bold", size = 14))
  )
Comparação entre o mapa de desmatamento bruto (esquerda) e o mapa LISA (direita). O mapa bruto mostra *onde* o desmatamento é alto; o LISA mostra *onde esse padrão é estatisticamente não-aleatório* e qual a relação com os municípios vizinhos.

Comparação entre o mapa de desmatamento bruto (esquerda) e o mapa LISA (direita). O mapa bruto mostra onde o desmatamento é alto; o LISA mostra onde esse padrão é estatisticamente não-aleatório e qual a relação com os municípios vizinhos.

6 Bloco 5 — Síntese e discussão

6.1 Resultados principais

I de Moran global para o desmatamento nos 807 municípios da Amazônia Legal.
Variável I de Moran p (analítico) p (permutação)
desmat_km2 (original) 0.2668 2.688484e-42 0,001
log_desmat (log-transf.) 0.4387 5.259296e-90 —
Distribuição dos clusters LISA (p < 0,05 · 999 permutações).
Quadrante Municípios Interpretacao %
HH 79 Alto cercado de alto — hotspot 9.8
LL 102 Baixo cercado de baixo — coldspot 12.6
HL 15 Outlier: alto cercado de baixo 1.9
LH 20 Outlier: baixo cercado de alto 2.5
Não sig. 591 Sem padrão espacial significativo 73.2

6.2 Pontos de atenção metodológicos

6.2.1 1. MAUP — Modifiable Areal Unit Problem

Os resultados são sensíveis à escala e ao recorte espacial. Clusters identificados em municípios podem se dissolver ou se intensificar ao analisar mesorregiões ou estados. A unidade municipal é uma escolha convencional, não neutra.

6.2.2 2. Escolha da matriz de pesos

Rainha, Torre, KNN-\(k\), distância — cada critério produz resultados distintos. A análise de sensibilidade (rodar com 2–3 especificações) é recomendada antes de publicar.

6.2.3 3. Problema da significância múltipla

Com 807 testes simultâneos, espera-se \(807 \times 0{,}05 \approx 40\) falsos positivos a \(p < 0{,}05\) por acaso. A correção FDR (Benjamini-Hochberg) é mais conservadora que Bonferroni e mais adequada para dados espacialmente correlacionados:

amz$lisa_p_fdr <- p.adjust(amz$lisa_p, method = "BH")
cat("HH com p < 0,05 (original):", sum(amz$quadrante == "HH"), "\n")
HH com p < 0,05 (original): 79 
cat("HH com p_FDR < 0,05       :",
    sum(amz$quadrante == "HH" & amz$lisa_p_fdr < 0.05), "\n")
HH com p_FDR < 0,05       : 23

6.2.4 4. Autocorrelação ≠ causalidade

Clusters HH indicam similaridade espacial, não mecanismo causal. Podem refletir:

  • Contágio: o desmatamento se espalha para municípios vizinhos (borda)
  • Homofilia: municípios com características similares ficam próximos
  • Variável de confusão: distância a portos/rodovias correlacionada

Para inferência causal, são necessários modelos de econometria espacial (SAR, SEM, SDM) com variáveis de controle adequadas.

6.2.5 5. Excesso de zeros

A transformação log1p() atenua a assimetria, mas municípios com desmatamento zero inflam o cluster LL. Em análises inferenciais, considere restringir a amostra a municípios com floresta nativa remanescente acima de um limiar mínimo.

7 Exercícios propostos

1. Repita a análise com o critério Torre (queen = FALSE). Os clusters LISA mudam? Em quais regiões a diferença é maior?

2. Experimente matrizes KNN com \(k = 3\), \(k = 5\) e \(k = 8\). Como o I de Moran global responde à mudança de \(k\)?

coords   <- st_coordinates(st_centroid(st_geometry(amz)))
nb_knn5  <- knearneigh(coords, k = 5) |> knn2nb()
lw_knn5  <- nb2listw(nb_knn5, style = "W")
moran.test(amz$log_desmat, lw_knn5)

3. Calcule o I de Moran apenas para os municípios do Pará. Lembre: a matriz de vizinhança precisa ser reconstruída sobre o subconjunto.

pa  <- amz |> filter(uf_sigla == "PA")
nb_pa <- poly2nb(pa, queen = TRUE)
lw_pa <- nb2listw(nb_pa, style = "W", zero.policy = TRUE)
moran.test(pa$log_desmat, lw_pa)

4. Compare os quadrantes LISA de log_desmat vs. desmat_km2 bruto. Quais municípios mudam de categoria? O que isso sugere sobre o impacto dos outliers?

5. Aplique a correção FDR e mapeie os clusters que sobrevivem:

amz$quadrante_fdr <- with(amz, case_when(
  lisa_p_fdr >= 0.05               ~ "Não sig.",
  z_padrao > 0 & z_lag > 0        ~ "HH",
  z_padrao < 0 & z_lag < 0        ~ "LL",
  z_padrao > 0 & z_lag < 0        ~ "HL",
  z_padrao < 0 & z_lag > 0        ~ "LH",
  TRUE                              ~ "Não sig."
))

8 Referências

Anselin, L. (1995). Local indicators of spatial association — LISA. Geographical Analysis, 27(2), 93–115. https://doi.org/10.1111/j.1538-4632.1995.tb00338.x

Bivand, R., Pebesma, E., & Gomez-Rubio, V. (2013). Applied Spatial Data Analysis with R (2a ed.). Springer.

Bivand, R., & Wong, D. (2018). Comparing implementations of global and local indicators of spatial association. TEST, 27(3), 716–748.

INPE. (2022). Projeto PRODES: Monitoramento do Desmatamento da Floresta Amazônica Brasileira por Satélite. Instituto Nacional de Pesquisas Espaciais. https://terrabrasilis.dpi.inpe.br

Pebesma, E. (2018). Simple features for R: Standardized support for spatial vector data. The R Journal, 10(1), 439–446.

Sessão R:
  R version: R version 4.5.2 (2025-10-31 ucrt) 
  sf     : 1.0.23 
  spdep  : 1.4.1 
  ggplot2: 4.0.1 
  dplyr  : 1.1.4

Documento produzido com R Markdown · DCAD/PPGOP – UFSM · 2025