📘 Memorial Descritivo – Jornada Analítica Espacial com Modelos Globais e Locais

Este notebook apresenta uma jornada completa de análise espacial com foco em padrões socioeconômicos simulados para Porto Alegre (RS). A partir de uma grade hexagonal urbana, dados de frequência de eventos e renda são atribuídos de forma controlada e contextual, permitindo explorar relações geográficas entre variáveis.

🔎 Objetivos

  • Simular variáveis socioeconômicas com variação espacial realista.
  • Aplicar regressões espaciais globais (SAR, SEM, SAC) e locais (GWR).
  • Identificar clusters significativos com estatísticas LISA univariadas e bivariadas.
  • Comparar modelos globais e locais para entender variações espaciais nos efeitos.
  • Interpretar resultados com visualizações e tabelas estilizadas.

🧩 Etapas Analíticas

  1. Construção da malha hexagonal urbana sobre o município de Porto Alegre com projeção UTM.
  2. Simulação de variáveis espaciais com hotspots controlados próximos a coordenadas reais, induzindo padrões espaciais na renda e na frequência.
  3. Aplicação de modelos SAR, SEM e SAC, estimando coeficientes globais e preparando estrutura de vizinhança para regressão espacial.
  4. Execução da regressão geograficamente ponderada (GWR), revelando coeficientes locais com variação espacial, ideais para identificar heterogeneidade territorial.
  5. Identificação dos extremos da GWR, com tabelas dos 25 hexágonos mais positivos e negativos.
  6. Classificação dos hotspots GWR com base nos quantis dos coeficientes locais, gerando mapas interpretativos.
  7. Comparação direta com o modelo SAR, destacando áreas onde GWR supera ou diverge fortemente do coeficiente global.
  8. Análise estatística da distribuição dos coeficientes GWR, incluindo medidas de assimetria e dispersão.
  9. Elaboração de tabelas interpretativas, incluindo implicações estatísticas e espaciais dos achados.

🗺️ Principais Insights

  • Os coeficientes locais variam de negativo a positivo, sugerindo forte heterogeneidade espacial.
  • O SAR estima um único coeficiente médio que subestima áreas com forte associação local.
  • A distribuição GWR apresenta cauda longa à direita, indicando que o efeito da renda é muito mais intenso em alguns bairros.
  • Hexágonos com diferenças marcantes entre GWR e SAR são potencialmente relevantes para políticas públicas regionalizadas.

🧠 Aplicabilidade

Este notebook serve como instrumento didático e técnico para pesquisadores, estudantes, analistas urbanos e profissionais da área de geoinformação. A abordagem permite explorar efeitos espaciais com profundidade, simulando dados com propósito interpretativo, e oferecendo recursos para replicação em outros contextos territoriais.

Bibliotecas

library(dplyr)
library(gt)
library(sf)
library(spdep)
library(geobr)
library(ggplot2)
library(mapview)
library(RColorBrewer)
library(sfdep)
library(spatialreg)
library(biscale)
library(purrr)
library(GWmodel)
library(moments)
🧭 Resumo da Jornada Analítica Espacial
📌 Etapa 🎯 Objetivo 🛠️ Técnica Utilizada 💡 Insight Gerado
🧱 Construção da grade hexagonal Representar o espaço urbano de Porto Alegre st_make_grid() + st_intersects() Base geográfica para modelagem espacial
🎲 Simulação de dados com hotspots Criar variáveis com padrões espacialmente realistas Funções de distância + exp(-d/λ) Hotspots e variação territorial simulada
📈 Modelos SAR/SEM/SAC Estimar efeito global da renda sobre frequência lagsarlm(), errorsarlm(), spatialreg Modelo global pode mascarar nuances locais
📡 Regressão GWR Estimar coeficiente local da renda em cada hexágono bw.gwr() + gwr.basic() GWR revela heterogeneidade espacial
🔍 Visualização e interpretação Destacar extremos e discrepâncias entre modelos ggplot2, mapview, gt Áreas onde SAR subestima o efeito da renda
🧠 Implicações estatísticas Avaliar impacto das variações espaciais na análise Skewness, dispersão, comparação GWR − SAR GWR detecta extremos que SAR suaviza

1. Carrega Porto Alegre (IBGE: 4314902)

poa <- read_municipality(code_muni = 4314902, year = 2020)

2. Reprojeção para UTM zona 22S

poa_utm <- st_transform(poa, crs = 31982)

3. Geração da grade hexagonal (~1km)

hex_grid <- st_make_grid(poa_utm, cellsize = 1000, square = FALSE)
hex_sf <- st_sf(id = 1:length(hex_grid), geometry = hex_grid)
hex_sf <- hex_sf[st_intersects(hex_sf, poa_utm, sparse = FALSE), ]

4. Seus 6 pontos de hotspot em latitude/longitude

coords_df <- data.frame(
  id = 1:6,
  lon = c(-51.2853, -51.1869, -51.1757, -51.1576, -51.2541, -51.1847),
  lat = c(-30.0144, -30.0250, -30.0258, -30.0214, -30.1007, -30.0455)
)

5. Converte para sf e reprojeta para UTM

pts_sf <- st_as_sf(coords_df, coords = c("lon","lat"), crs = 4326)
pts_sf <- st_transform(pts_sf, crs = 31982)

6. Calcula distância do centro dos hexágonos até os pontos

hex_centroid <- st_centroid(hex_sf)
dist_matrix <- st_distance(hex_centroid, pts_sf)
min_dist <- apply(dist_matrix, 1, min)  # menor distância por hexágono

7. Simula frequência com decaimento exponencial + ruído aleatório

set.seed(123)
hex_sf$frequencia <- round(120 * exp(-min_dist / 1500)) +
  sample(1:30, nrow(hex_sf), replace = TRUE)

8. Visualização com ggplot

ggplot(hex_sf) +
  geom_sf(aes(fill = frequencia), color = "gray40") +
  scale_fill_viridis_c(option = "C") +
  labs(
    title = "Frequência Simulada com Hotspots",
    subtitle = "Valores mais altos próximos às coordenadas fornecidas",
    fill = "Frequência"
  ) +
  theme_minimal()

#st_crs(hex_sf)
hex_sf <- st_transform(hex_sf, crs = 4326)
summary(st_is_valid(hex_sf))
##    Mode    TRUE 
## logical     678
#hex_sf <- st_make_valid(hex_sf)

names(hex_sf)[names(hex_sf) == "frequencia"] <- "frequencia"

mapview(hex_sf, zcol = "frequencia")
plot(hex_sf["frequencia"], main = "Hexágonos com frequencia")

Vizinhança por toque (queen contiguity)

nb <- poly2nb(as(hex_sf, "Spatial"), queen = TRUE)
listw <- nb2listw(nb, style = "W", zero.policy = TRUE)

moran.test(hex_sf$frequencia, listw, zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  hex_sf$frequencia  
## weights: listw  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 29.138, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.6810788515     -0.0014792899      0.0005487168
lisa <- localmoran(hex_sf$frequencia, listw, zero.policy = TRUE)

Adiciona os resultados ao objeto sf

hex_sf$lisa_I <- lisa[,1]
hex_sf$lisa_p <- lisa[,5]

media_freq <- mean(hex_sf$frequencia)

hex_sf$cluster <- ifelse(
  hex_sf$frequencia > media_freq & hex_sf$lisa_I > 0, "Alta-Alta",
  ifelse(hex_sf$frequencia < media_freq & hex_sf$lisa_I > 0, "Baixa-Baixa",
         ifelse(hex_sf$frequencia > media_freq & hex_sf$lisa_I < 0, "Alta-Baixa",
                ifelse(hex_sf$frequencia < media_freq & hex_sf$lisa_I < 0, "Baixa-Alta", "Não Significativo")))
)
ggplot(hex_sf) +
  geom_sf(aes(fill = cluster), color = "gray30") +
  scale_fill_brewer(palette = "Set1") +
  theme_minimal() +
  labs(title = "Mapa LISA com Hexágonos",
       fill = "Cluster Espacial")

LISA Bivariado — Explorando Associações Entre Duas Variáveis

📌 Exemplo:

Mapa de hexágonos com:
frequencia de eventos
renda média da população
Permitindo investigar se áreas com alta frequencia de eventos estão próximas de áreas com alta ou baixa renda.
Supõe duas variáveis no hex_sf: frequencia e renda

1. Reutilizar os centros de alta renda (pts_sf já existe com 6 pontos)

set.seed(123)
hex_centroid <- st_centroid(hex_sf)

2. Calcular distância dos hexágonos até os centros

pts_sf <- st_transform(pts_sf, st_crs(hex_centroid))
dist_matrix_renda <- st_distance(hex_centroid, pts_sf)
min_dist_renda <- apply(dist_matrix_renda, 1, min)

3. Gerar renda maior próxima dos centros

hex_sf$renda <- round(1500 * exp(-min_dist_renda / 1500)) +
  sample(100:300, nrow(hex_sf), replace = TRUE)

hex_sf$renda = hex_sf$renda * hex_sf$frequencia
cor(hex_sf$renda, hex_sf$frequencia)
## [1] 0.8921704

4. Simular frequência como função da renda + ruído local

Cria a matriz de pesos espaciais como antes

nb <- poly2nb(as(hex_sf, "Spatial"), queen = TRUE)

Verifica número de vizinhos por polígono

num_vizinhos <- card(nb)

Filtra apenas os hexágonos com pelo menos 1 vizinho

hex_sf <- hex_sf[num_vizinhos > 0, ]
nb <- subset(nb, num_vizinhos > 0)
listw <- nb2listw(nb, style = "W", zero.policy = TRUE)

Calcula o LISA bivariado: frequencia vs renda

# Usando sfdep::local_moran_bv()
lisa_biv <- local_moran_bv(
  x = hex_sf$frequencia,
  y = hex_sf$renda,
  nb = nb,
  wt = listw$weights,
  nsim = 499
)

lisa_biv agora é um data.frame com $Ib e $p_sim

hex_sf$lisa_biv_I <- lisa_biv$Ib
hex_sf$lisa_biv_p <- lisa_biv$p_sim

Regressão Espacial — Modelando Influência de Fatores com Dependência Espacial

Quando há suspeita de que a frequencia é influenciada por variáveis como renda, densidade populacional, etc., mas com dependência espacial, é melhor usar um modelo de regressão espacial.

📊 Tabela Comparativa de Modelos de Regressão Espacial
🔍 Modelo 📐 Lag na Y 📶 Erro Espacial 🛠️ Função 📘 Nome Técnico
SAR ✅ ❌ lagsarlm() Lag Espacial
SEM ❌ ✅ errorsarlm() Erro Espacial
SAC ✅ ✅ sacsarlm() Combinado (SARAR)

1. Modelo de regressão espacial tipo SAR (Spatial Autoregressive)

(modelo_sar <- lagsarlm(frequencia ~ renda, data = hex_sf, listw = listw, zero.policy = TRUE)) |> summary()
## 
## Call:lagsarlm(formula = frequencia ~ renda, data = hex_sf, listw = listw, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -24.10216  -6.17984   0.49074   6.14984  19.77829 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 1.2011e+01 7.4242e-01  16.178 < 2.2e-16
## renda       6.4499e-04 2.7728e-05  23.262 < 2.2e-16
## 
## Rho: 0.22042, LR test value: 38.205, p-value: 6.3704e-10
## Asymptotic standard error: 0.037917
##     z-value: 5.8132, p-value: 6.1275e-09
## Wald statistic: 33.794, p-value: 6.1275e-09
## 
## Log likelihood: -2346.834 for lag model
## ML residual variance (sigma squared): 59.477, (sigma: 7.7121)
## Number of observations: 677 
## Number of parameters estimated: 4 
## AIC: 4701.7, (AIC for lm: 4737.9)
## LM test for residual autocorrelation
## test value: 0.013213, p-value: 0.90849

2. Modelo SEM (Spatial Error Model)

O modelo SEM considera que os erros da regressão estão autocorrelacionados espacialmente, ou seja, há dependência espacial nos resíduos.

summary(
  modelo_sem <- errorsarlm(frequencia ~ renda, data = hex_sf, listw = listw, zero.policy = TRUE)
)
## 
## Call:errorsarlm(formula = frequencia ~ renda, data = hex_sf, listw = listw, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -31.72974  -6.10205   0.56253   6.12617  20.34556 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 1.6049e+01 4.6491e-01  34.520 < 2.2e-16
## renda       7.6492e-04 1.8932e-05  40.403 < 2.2e-16
## 
## Lambda: 0.29846, LR test value: 29.174, p-value: 6.6144e-08
## Asymptotic standard error: 0.053936
##     z-value: 5.5336, p-value: 3.1371e-08
## Wald statistic: 30.621, p-value: 3.1371e-08
## 
## Log likelihood: -2351.349 for error model
## ML residual variance (sigma squared): 59.769, (sigma: 7.731)
## Number of observations: 677 
## Number of parameters estimated: 4 
## AIC: 4710.7, (AIC for lm: 4737.9)

3. Modelo SAC (Spatial Autoregressive Combined)

O modelo SAC (também chamado SARAR) combina:
Lag espacial na variável dependente (efeito de vizinhos sobre o valor observado)
Erro espacial autocorrelacionado (como no SEM)

summary(
modelo_sac <- sacsarlm(frequencia ~ renda, data = hex_sf, listw = listw, zero.policy = TRUE)
)
## 
## Call:sacsarlm(formula = frequencia ~ renda, data = hex_sf, listw = listw, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -23.54233  -6.22402   0.43718   6.13449  19.70200 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 1.1790e+01 9.3980e-01  12.545 < 2.2e-16
## renda       6.3684e-04 3.4785e-05  18.308 < 2.2e-16
## 
## Rho: 0.23318
## Asymptotic standard error: 0.050556
##     z-value: 4.6122, p-value: 3.9842e-06
## Lambda: -0.025304
## Asymptotic standard error: 0.084459
##     z-value: -0.2996, p-value: 0.76448
## 
## LR test value: 38.239, p-value: 4.9717e-09
## 
## Log likelihood: -2346.816 for sac model
## ML residual variance (sigma squared): 59.397, (sigma: 7.7069)
## Number of observations: 677 
## Number of parameters estimated: 5 
## AIC: 4703.6, (AIC for lm: 4737.9)

Etapas para Criar Mapas Bivariados com Quadrantes Cruzados

1. Classificar as variáveis em categorias

Transformamos frequencia e renda em classes (ex: baixa, média, alta) usando quantis:
Classificar em 3 categorias usando quantis

summary(hex_sf)
##        id                  geometry     frequencia         lisa_I         
##  Min.   :  21.0   POLYGON      :677   Min.   :  1.00   Min.   :-0.424589  
##  1st Qu.: 415.0   epsg:4326    :  0   1st Qu.: 12.00   1st Qu.:-0.007362  
##  Median : 630.0   +proj=long...:  0   Median : 21.00   Median : 0.099809  
##  Mean   : 621.5                       Mean   : 23.78   Mean   : 0.682085  
##  3rd Qu.: 820.0                       3rd Qu.: 29.00   3rd Qu.: 0.323980  
##  Max.   :1281.0                       Max.   :121.00   Max.   :19.088614  
##      lisa_p         cluster              renda          lisa_biv_I       
##  Min.   :0.0000   Length:677         Min.   :   103   Min.   :-0.181536  
##  1st Qu.:0.1663   Class :character   1st Qu.:  2688   1st Qu.:-0.004036  
##  Median :0.3833   Mode  :character   Median :  4484   Median : 0.097162  
##  Mean   :0.4129                      Mean   : 10132   Mean   : 0.706905  
##  3rd Qu.:0.6377                      3rd Qu.:  7830   3rd Qu.: 0.276650  
##  Max.   :0.9985                      Max.   :177000   Max.   :23.141492  
##    lisa_biv_p    
##  Min.   :0.0000  
##  1st Qu.:0.0420  
##  Median :0.1400  
##  Mean   :0.1845  
##  3rd Qu.:0.3160  
##  Max.   :0.5000
quantile(hex_sf$renda)
##     0%    25%    50%    75%   100% 
##    103   2688   4484   7830 177000
hex_sf <- hex_sf %>%
  mutate(
    freq_cat = ntile(frequencia, 3),
    renda_cat = ntile(renda, 3),
    bi_class = paste0(freq_cat, "-", renda_cat)
  )
quantile(hex_sf$renda_cat)
##   0%  25%  50%  75% 100% 
##    1    1    2    3    3

2. Criar paleta de cores bivariada

Utilizando o pacote biscale para aplicar uma paleta que mistura duas escalas:

# Classificar com biscale
hex_sf <- bi_class(hex_sf, x = frequencia, y = renda, style = "quantile", dim = 3)

3. Plotar o mapa com ggplot2

ggplot(hex_sf) +
  geom_sf(aes(fill = bi_class), color = "gray30") +
  bi_scale_fill(pal = "DkBlue", dim = 3) +
  labs(title = "Mapa Bivariado: frequencia vs Renda",
       fill = "frequencia × Renda") +
  theme_minimal()

3. Plotar o mapa com mapview

mapView(hex_sf, zcol="bi_class")
📌 Interpretação dos Quadrantes Bivariados
🔀 Quadrante 🧭 Significado 🔣 Ícone
1-1 Baixa frequencia e baixa renda ⬇️⬇️
1-3 Baixa frequencia e alta renda (outlier) ⬇️⬆️
3-1 Alta frequencia e baixa renda (outlier) ⬆️⬇️
3-3 Alta frequencia e alta renda ⬆️⬆️

hexágonos com associação significativa

significativos = hex_sf %>%
  filter(lisa_p < 0.05) %>%
  select(id, frequencia, renda, lisa_I, cluster, lisa_p)
significativos
## Simple feature collection with 80 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -51.27259 ymin: -30.21987 xmax: -51.08063 ymax: -29.99062
## Geodetic CRS:  WGS 84
## First 10 features:
##     id frequencia  renda     lisa_I   cluster       lisa_p
## 1  187         74  57646  9.4042464 Alta-Alta 8.161872e-09
## 2  209         49  37436  3.8722981 Alta-Alta 5.527719e-08
## 3  210         88  97768 14.4342524 Alta-Alta 3.719880e-12
## 4  215         36  20736  0.6070451 Alta-Alta 4.951020e-02
## 5  230         50  30050  1.9173439 Alta-Alta 9.583912e-03
## 6  231        109 126767 13.4904515 Alta-Alta 2.606603e-12
## 7  253         53  38107  3.6559479 Alta-Alta 5.454328e-08
## 8  254         98 123480 13.5153840 Alta-Alta 2.744103e-13
## 9  255         41  26117  1.4350265 Alta-Alta 1.077834e-02
## 10 274         40  20000  0.9429464 Alta-Alta 1.167916e-02
##                          geometry
## 1  POLYGON ((-51.2674 -30.1103...
## 2  POLYGON ((-51.26223 -30.118...
## 3  POLYGON ((-51.26219 -30.102...
## 4  POLYGON ((-51.26199 -30.024...
## 5  POLYGON ((-51.25706 -30.125...
## 6  POLYGON ((-51.25702 -30.110...
## 7  POLYGON ((-51.25185 -30.118...
## 8  POLYGON ((-51.25181 -30.102...
## 9  POLYGON ((-51.25178 -30.086...
## 10 POLYGON ((-51.24668 -30.125...
mapview(significativos, zcol="cluster")

hexágonos com associação bivariada significativa

significativos_bivariado = hex_sf %>%
  filter(lisa_biv_p < 0.05) %>%
  select(id, frequencia, renda, lisa_biv_I,lisa_biv_p, bi_class)
mapview(significativos_bivariado, zcol="bi_class")
🧭 Como Interpretar Cada Hexágono
📦 Coluna 🧬 Significado 🔍 Interpretação
frequencia Valor da variável principal (ex: número de eventos) Quanto maior o valor, maior a concentração de eventos naquele hexágono
renda Valor da variável secundária (ex: renda média) Usada para comparação cruzada com frequencia (ex: áreas ricas vs. pobres)
lisa_I Índice de Moran Local (autocorrelação univariada) Se positivo e significativo, indica agrupamento com vizinhos similares
lisa_p Significância estatística da autocorrelação local (p < 0.05 é significativo) Permite identificar clusters significativos ou padrões aleatórios
cluster Classificação do padrão espacial: Alta-Alta, Baixa-Baixa, outliers ou neutro Mostra se o hexágono pertence a um cluster ou é um outlier espacial
lisa_biv_I Índice de autocorrelação espacial entre duas variáveis Se positivo e significativo, há correlação local entre frequencia e renda
lisa_biv_p Significância do índice bivariado entre frequencia e renda Se < 0.05, o padrão entre frequencia e renda é estatisticamente relevante

Subset

Selecionar Hexágonos com Cluster Alta-Alta e P-Valor Significativo (LISA uni e bivariado)
Filtro: cluster Alta-Alta com p < 0.05 tanto para lisa univariado quanto bivariado

hex_aa_sig <- hex_sf %>%
  filter(
    cluster == "Alta-Alta",
    lisa_p < 0.05,
    lisa_biv_p < 0.05
  )

Visualiza os hexágonos selecionados

mapview(hex_aa_sig, zcol = "frequencia")

Pegando os vizinhos do hexágono Alta-Alta
Pegar os índices das linhas (não os IDs) dos hexágonos Alta-Alta

linhas_alta <- which(hex_sf$id %in% hex_aa_sig$id)

Coletar todos vizinhos corretamente

vizinhos_idx <- linhas_alta %>%
  map(~ nb[[.x]]) %>%
  unlist() %>%
  unique()

Remover os próprios hexágonos Alta-Alta

vizinhos_idx <- setdiff(vizinhos_idx, linhas_alta)

Extrair os vizinhos da camada

hex_vizinhos <- hex_sf[vizinhos_idx, ]

Vizinhos

mapview(hex_vizinhos, zcol = "frequencia")

Hot-spots com vizinhos

mapview(hex_vizinhos, zcol = "frequencia") + 
  mapview(hex_aa_sig, col.regions = "red")

GWR

SAR/SEM/SAC usam a estrutura de vizinhança (matriz W) para modelar o efeito espacial.Ambos tentam incorporar a dimensão espacial na modelagem estatística. GWR usa pesos geográficos (kernels) que variam conforme a posição no espaço.
SAR/SEM/SAC são “modelos espaciais globais”, enquanto GWR é um “modelo espacial local” — ideal para quando os efeitos variam geograficamente.

🌐 Diferenças e Relações entre GWR e SAR/SEM/SAC
📊 Característica 📦 SAR / SEM / SAC 🗺️ GWR
🔗 Tipo de Dependência Espacial Global, via matriz de vizinhança (listw) Local, via kernel de distância
📐 Estimação dos Coeficientes Coeficiente único por variável (modelo global) Coeficientes que variam por localidade
🧠 Complexidade Computacional Moderada (paramétrico) Alta (ajuste local em cada polígono)
📦 Pacotes Comuns `spdep`, `spatialreg` `GWmodel`, `spgwr`, `mgwr`
🎯 Tipo de Modelo Modelos espaciais globais Modelos espaciais locais
🗺️ Quando Utilizar Quando há padrão global de autocorrelação Quando há heterogeneidade espacial dos efeitos

Extrair coordenadas dos centroides

#devtools::install_github("GWmodel-Lab/GWmodel3")
#library(GWmodel3)
hex_coords <- st_coordinates(st_centroid(hex_sf))

Adicionar ao data.frame

hex_df <- hex_sf %>%
  st_drop_geometry() %>%
  mutate(X = hex_coords[,1], Y = hex_coords[,2])

Fórmula: frequencia explicada por renda

formula_gwr <- frequencia ~ renda

Bandwidth ideal (pode demorar um pouco)

bw.gwr() primeiro encontra o tamanho ideal de vizinhança

bw <- bw.gwr(formula_gwr, data = hex_sf, adaptive = TRUE)
## Adaptive bandwidth: 426 CV score: 34329.31 
## Adaptive bandwidth: 271 CV score: 29805.6 
## Adaptive bandwidth: 175 CV score: 26737.94 
## Adaptive bandwidth: 115 CV score: 24125.55 
## Adaptive bandwidth: 79 CV score: 21650.34 
## Adaptive bandwidth: 55 CV score: 19595.08 
## Adaptive bandwidth: 42 CV score: 18736.77 
## Adaptive bandwidth: 32 CV score: 18061.83 
## Adaptive bandwidth: 28 CV score: 18042.92 
## Adaptive bandwidth: 23 CV score: 18107.42 
## Adaptive bandwidth: 28 CV score: 18042.92

Aplicar GWR

Com esse bw, se alimenta a gwr.basic() para estimar os coeficientes locais

Analogia simples
bw.gwr() é como escolher o alcance do radar 🛰️ (quantos hexágonos em volta cada ponto considera).
gwr.basic() é como fazer a leitura personalizada para cada local 🗺️ (calculando os efeitos da variável em cada região).

modelo_gwr <- gwr.basic(formula_gwr, data = hex_sf, bw = bw, adaptive = TRUE)
summary(modelo_gwr)
##               Length Class  Mode
## GW.arguments  11     -none- list
## GW.diagnostic  8     -none- list
## lm            14     lm     list
## SDF           13     sf     list
## timings        5     -none- list
## this.call      5     -none- call
## Ftests         0     -none- list
tabela_gwr_funcoes <- tibble::tibble(
  Função = c("bw.gwr()", "gwr.basic()"),
  Finalidade = c(
    "Selecionar a largura de banda ótima (quantos vizinhos considerar)",
    "Aplicar a regressão GWR com os dados e a largura de banda definida"
  ),
  Entrada = c(
    "Fórmula, dados, coordenadas; opções de kernel e adaptatividade",
    "Fórmula, dados, coordenadas, bandwidth; opções de kernel"
  ),
  Retorno = c(
    "Valor numérico de bandwidth ideal",
    "Objeto com coeficientes locais, estatísticas e predições"
  ),
  Analogia = c(
    "Ajusta o alcance do radar espacial 🛰️",
    "Faz a leitura local personalizada em cada hexágono 📡"
  )
)

# Monta a tabela com gt
tabela_gwr_funcoes %>%
  gt() %>%
  tab_header(
    title = "🧠 Comparação entre `bw.gwr()` e `gwr.basic()` no GWmodel"
  ) %>%
  cols_label(
    Função = "🛠️ Função",
    Finalidade = "🎯 Finalidade",
    Entrada = "📥 Entradas",
    Retorno = "📤 Retorno",
    Analogia = "💡 Analogia"
  )
🧠 Comparação entre `bw.gwr()` e `gwr.basic()` no GWmodel
🛠️ Função 🎯 Finalidade 📥 Entradas 📤 Retorno 💡 Analogia
bw.gwr() Selecionar a largura de banda ótima (quantos vizinhos considerar) Fórmula, dados, coordenadas; opções de kernel e adaptatividade Valor numérico de bandwidth ideal Ajusta o alcance do radar espacial 🛰️
gwr.basic() Aplicar a regressão GWR com os dados e a largura de banda definida Fórmula, dados, coordenadas, bandwidth; opções de kernel Objeto com coeficientes locais, estatísticas e predições Faz a leitura local personalizada em cada hexágono 📡

Coeficiente da renda por hexágono

hex_sf$coef_renda <- modelo_gwr$SDF$renda

Mapa dos coeficientes locais

Valores positivos indicam que a renda está associada a maior frequência de eventos naquele local.
Valores negativos indicam que a renda está associada a menor frequência — ou seja, pode haver um padrão inverso.

ggplot(hex_sf) +
  geom_sf(aes(fill = coef_renda), color = "gray30") +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Coeficiente Local da Renda (GWR)",
       fill = "Efeito da Renda") +
  theme_minimal()

Etapas para Comparação: GWR vs Modelos Globais

  1. 📦 Resgatar o coeficiente global dos modelos, basta extrair os coeficientes da variável renda:
#coef_sar <- summary(modelo_sar)$coefficients["renda", "Estimate"]
#coef_sem <- summary(modelo_sem)$coefficients["renda", "Estimate"]
#coef_sac <- summary(modelo_sac)$coefficients["renda", "Estimate"]

coef_sar <- summary(modelo_sar)$coefficients["renda"]
coef_sem <- summary(modelo_sem)$coefficients["renda"]
coef_sac <- summary(modelo_sac)$coefficients["renda"]

Resumo dos coeficientes locais (mínimo, máximo, média)

📡 Coeficientes Locais da Renda (GWR)
Estatística GWR
Min 0.0004745822
Média 0.002668102
Max 0.005093496

25 maiores coeficientes

top25_gwr <- hex_sf %>%
  arrange(desc(coef_renda)) %>%
  slice_head(n = 25) %>%
  st_drop_geometry() %>%
  select(ID_top = id, Renda_top = coef_renda)

25 menores coeficientes

bottom25_gwr <- hex_sf %>%
  arrange(coef_renda) %>%
  slice_head(n = 25) %>%
  st_drop_geometry() %>%
  select(ID_bottom = id, Renda_bottom = coef_renda)

Combina as duas tabelas lado a lado

gwr_extremos <- bind_cols(top25_gwr, bottom25_gwr)
📊 Extremos dos Coeficientes Locais da Renda (GWR)
🔺 Top 25 – ID 🔺 Coef. Local Renda 🔻 Bottom 25 – ID 🔻 Coef. Local Renda
242 0.00509 544 0.00047
286 0.00504 567 0.00048
264 0.00503 588 0.00048
330 0.00496 566 0.00049
308 0.00492 611 0.00049
241 0.00492 543 0.00049
352 0.00487 632 0.00049
646 0.00481 522 0.00050
220 0.00481 655 0.00051
577 0.00480 610 0.00051
198 0.00479 523 0.00052
285 0.00478 676 0.00052
511 0.00476 587 0.00052
667 0.00475 500 0.00053
533 0.00473 589 0.00053
329 0.00471 654 0.00055
555 0.00467 565 0.00055
690 0.00466 499 0.00055
800 0.00463 699 0.00055
840 0.00461 521 0.00055
600 0.00455 545 0.00055
373 0.00452 633 0.00055
863 0.00451 478 0.00055
1281 0.00450 677 0.00056
1259 0.00450 298 0.00057

Quantis

quantile(hex_sf$coef_renda)
##           0%          25%          50%          75%         100% 
## 0.0004745822 0.0015889106 0.0029120698 0.0036203441 0.0050934956

Interpretação

Se os coeficientes da GWR tiverem ampla variação (ex: de negativo a positivo), isso revela heterogeneidade espacial — ou seja, o efeito da renda não é uniforme em Porto Alegre.
Os modelos SAR/SEM/SAC, por serem globais, estimam um valor médio que pode mascarar efeitos locais fortes ou inversos.

🗺️ Etapas para Visualização dos “Hotspots GWR”
1. Criar classificação baseada no coeficiente GWR

# hex_sf <- hex_sf %>%
#   mutate(
#     efeito_renda_gwr = case_when(
#       coef_renda > quantile(coef_renda, 0.9, na.rm = TRUE) ~ "Efeito Forte Positivo",
#       coef_renda < quantile(coef_renda, 0.1, na.rm = TRUE) ~ "Efeito Fraco Positivo",
#       TRUE ~ "Efeito Moderado"
#     )
#   )
hex_sf <- hex_sf %>%
  mutate(
    efeito_renda_gwr = case_when(
      coef_renda > quantile(coef_renda, 0.85, na.rm = TRUE) ~ "Efeito Forte Positivo",
      coef_renda < quantile(coef_renda, 0.15, na.rm = TRUE) ~ "Efeito Fraco Positivo",
      TRUE ~ "Efeito Moderado"
    )
  )
  1. Mapear os efeitos locais da renda
ggplot(hex_sf) +
  geom_sf(aes(fill = efeito_renda_gwr), color = "gray50") +
  scale_fill_manual(values = c(
    "Efeito Forte Positivo" = "#2166AC",  # azul escuro
    "Efeito Fraco Positivo" = "#B2182B",  # vermelho escuro
    "Efeito Moderado" = "gray80"
  )) +
  labs(
    title = "Mapa GWR: Efeito Local da Renda sobre Frequência",
    fill = "Força do Efeito"
  ) +
  theme_minimal()

mapView(hex_sf, zcol="efeito_renda_gwr")

Por que um hexágono pode ter alta renda e alta frequência, mas efeito negativo da renda:

Porque a associação estatística entre renda e frequência depende da distribuição local, e o GWR capta como a variação da renda afeta a variação da frequência.
Exemplo hipotético:
- Um hexágono tem renda alta e frequência alta → entra no quadrante 3-3.
- Mas nas redondezas, hexágonos vizinhos com renda ainda mais alta não apresentam aumentos proporcionais na frequência.
- Nesse cenário, a GWR pode detectar um coeficiente negativo, ou seja: em áreas com mais renda, a frequência tende a cair.
Isso não contradiz os valores altos absolutos, mas indica um padrão inverso na associação local.

hex_sf$coef_renda <- modelo_gwr$SDF$renda
# Gera a variável de sinal (positivo, negativo, nulo)
hex_sf <- hex_sf %>%
  mutate(
    sinal_gwr = case_when(
      coef_renda > 0 ~ "positivo",
      coef_renda < 0 ~ "negativo",
      TRUE ~ "nulo"
    ),
    classe_cruzada = paste0(bi_class, "_", sinal_gwr)
  )

# Conta combinações
tabela_cruzada <- hex_sf %>%
  count(bi_class, sinal_gwr) %>%
  rename(Classe_Bivariada = bi_class, Sinal_GWR = sinal_gwr, Hexágonos = n)
tabela_cruzada$geometry = NULL
# Exibe com gt
tabela_cruzada[,1:3] %>%
  gt() %>%
  tab_header(title = "📌 Cruzamento entre Classe Bivariada e Sinal do Efeito da Renda (GWR)") %>%
  cols_label(
    Classe_Bivariada = "🔀 Classe Bivariada (frequência × renda)",
    Sinal_GWR = "📈 Sinal do Coef. GWR",
    Hexágonos = "🔢 Nº de Hexágonos"
  ) %>%
  fmt_number(columns = Hexágonos, decimals = 0)
📌 Cruzamento entre Classe Bivariada e Sinal do Efeito da Renda (GWR)
🔀 Classe Bivariada (frequência × renda) 📈 Sinal do Coef. GWR 🔢 Nº de Hexágonos
1-1 positivo 196
1-2 positivo 44
2-1 positivo 29
2-2 positivo 144
2-3 positivo 46
3-1 positivo 1
3-2 positivo 37
3-3 positivo 180

Criar uma nova variável que combine o quadrante bivariado com o sinal do coeficiente GWR. Assim fica fácil ver onde existe “associação alta-alta com efeito inverso”.

hex_sf %>%
  filter(bi_class == "3-3") %>%
  select(id, frequencia, renda, coef_renda) %>%
  arrange(coef_renda)
## Simple feature collection with 180 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -51.29837 ymin: -30.25121 xmax: -51.04428 ymax: -29.95918
## Geodetic CRS:  WGS 84
## First 10 features:
##     id frequencia  renda   coef_renda                       geometry
## 1  544        121 172183 0.0004745822 POLYGON ((-51.18422 -30.032...
## 2  567        100 100700 0.0004793492 POLYGON ((-51.17902 -30.024...
## 3  588        106 150626 0.0004815490 POLYGON ((-51.17385 -30.032...
## 4  566         76  74404 0.0004881768 POLYGON ((-51.17905 -30.040...
## 5  611         69  66585 0.0004929488 POLYGON ((-51.16865 -30.024...
## 6  543        120 177000 0.0004933899 POLYGON ((-51.18425 -30.047...
## 7  632         85  86445 0.0004938524 POLYGON ((-51.16348 -30.032...
## 8  522         68  66096 0.0004993409 POLYGON ((-51.18942 -30.040...
## 9  655        118 171926 0.0005084106 POLYGON ((-51.15828 -30.024...
## 10 610         59  50563 0.0005092856 POLYGON ((-51.16868 -30.040...

Criar uma nova variável que combine o quadrante bivariado com o sinal do coeficiente GWR. Assim fica fácil ver onde existe “alta-alta com efeito inverso”:

hex_sf <- hex_sf %>%
  mutate(
    sinal_gwr = case_when(
      coef_renda > 0 ~ "positivo",
      coef_renda < 0 ~ "negativo",
      TRUE ~ "nulo"
    ),
    classe_cruzada = paste0(bi_class, "_", sinal_gwr)
  )

ggplot(hex_sf) +
  geom_sf(aes(fill = classe_cruzada), color = "gray30") +
  scale_fill_manual(
    values = c(
      "3-3_negativo" = "#b2182b",
      "3-3_positivo" = "#2166ac",
      "1-1_positivo" = "#67a9cf",
      "1-1_negativo" = "#ef8a62"
      # adicione outras combinações se quiser
    )
  ) +
  labs(
    title = "🌐 Mapa de Conflito: Onde a Associação é paralela e o Efeito é forte",
    fill = "Classe Cruzada"
  ) +
  theme_minimal()

mapView(hex_sf, zcol="classe_cruzada")
  1. Contrastar com o modelo SAR
ggplot(hex_sf, aes(x = coef_renda)) +
  geom_histogram(bins = 30, fill = "#67a9cf", color = "white") +
  geom_vline(xintercept = coef_sar, linetype = "dashed", color = "red", size = 1.2) +
  labs(
    title = "Distribuição dos Coeficientes Locais (GWR)",
    x = "Coeficiente GWR da Renda",
    y = "Número de Hexágonos"
  ) +
  theme_minimal()

A linha vermelha mostra o valor global estimado pelo SAR.
Se a distribuição GWR for bem dispersa em torno dele, significa que o modelo global “apaga” nuances locais!

O que significa a linha vermelha deslocada para a esquerda? A linha vermelha representa o coeficiente global da renda estimado via SAR. Ela estando mais à esquerda sugere que o modelo SAR estimou um efeito da renda mais fraco ou menos positivo, comparado à maioria dos coeficientes locais.
Isso pode indicar que:
Há zonas com efeito positivo intenso da renda, que puxam a cauda direita da distribuição GWR.
O efeito da renda não é uniforme — algumas áreas possuem associação muito maior entre renda e frequência.
O que significa a cauda longa à direita na distribuição GWR?
Uma cauda direita mais pronunciada sugere que existem hexágonos onde o coeficiente local de renda é especialmente elevado, ou seja, a renda está # muito fortemente associada à frequência naquele local.
Isso pode estar ligado a zonas de alta renda com densidade alta de eventos — ou padrões sociais/geográficos que amplificam o impacto da renda ali.

Identificar os hexágonos na cauda direita para sabe onde estão os trechos com maior efeito da renda.

hex_sf %>%
  filter(coef_renda > quantile(coef_renda, 0.95, na.rm = TRUE))
## Simple feature collection with 34 features and 15 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -51.27212 ymin: -30.21994 xmax: -51.00272 ymax: -29.92794
## Geodetic CRS:  WGS 84
## First 10 features:
##     id                       geometry frequencia      lisa_I     lisa_p
## 1  198 POLYGON ((-51.26694 -29.938...         11  0.70168933 0.09238787
## 2  220 POLYGON ((-51.26178 -29.946...         10  0.52802168 0.09577714
## 3  241 POLYGON ((-51.25662 -29.954...         26 -0.07613455 0.13996313
## 4  242 POLYGON ((-51.25658 -29.938...          2  0.80573464 0.18922025
## 5  264 POLYGON ((-51.25142 -29.946...          6  0.46253032 0.25736731
## 6  285 POLYGON ((-51.24626 -29.954...          9  0.28949091 0.39364851
## 7  286 POLYGON ((-51.24622 -29.938...         22  0.06819416 0.17230032
## 8  307 POLYGON ((-51.2411 -29.9618...         20  0.12022807 0.16625800
## 9  308 POLYGON ((-51.24106 -29.946...         25 -0.05198395 0.06809727
## 10 329 POLYGON ((-51.23589 -29.954...          3  0.66302039 0.16477123
##        cluster renda  lisa_biv_I lisa_biv_p freq_cat renda_cat bi_class
## 1  Baixa-Baixa  1199  0.30087518      0.026        1         1      1-1
## 2  Baixa-Baixa  2590  0.29701948      0.012        1         1      1-1
## 3   Alta-Baixa  5148 -0.04720336      0.012        2         2      2-2
## 4  Baixa-Baixa   510  0.48262458      0.020        1         1      1-1
## 5  Baixa-Baixa  1818  0.35547729      0.036        1         1      1-1
## 6  Baixa-Baixa  2322  0.28217087      0.060        1         1      1-1
## 7  Baixa-Baixa  3366  0.03857014      0.036        2         2      2-2
## 8  Baixa-Baixa  5320  0.07983730      0.022        2         2      2-2
## 9   Alta-Baixa  4150 -0.02681766      0.010        2         2      2-2
## 10 Baixa-Baixa   534  0.40774225      0.048        1         1      1-1
##     coef_renda      efeito_renda_gwr sinal_gwr classe_cruzada
## 1  0.004788730 Efeito Forte Positivo  positivo   1-1_positivo
## 2  0.004807184 Efeito Forte Positivo  positivo   1-1_positivo
## 3  0.004921306 Efeito Forte Positivo  positivo   2-2_positivo
## 4  0.005093496 Efeito Forte Positivo  positivo   1-1_positivo
## 5  0.005028641 Efeito Forte Positivo  positivo   1-1_positivo
## 6  0.004778180 Efeito Forte Positivo  positivo   1-1_positivo
## 7  0.005038281 Efeito Forte Positivo  positivo   2-2_positivo
## 8  0.004438433 Efeito Forte Positivo  positivo   2-2_positivo
## 9  0.004923139 Efeito Forte Positivo  positivo   2-2_positivo
## 10 0.004706058 Efeito Forte Positivo  positivo   1-1_positivo

Mapear a diferença GWR – SAR:

hex_sf <- hex_sf %>%
  mutate(dif_coef = coef_renda - coef_sar)

dif_coef para ver em quais regiões o SAR subestimou ou superestimou o efeito:

ggplot(hex_sf) +
  geom_sf(aes(fill = dif_coef), color = "gray20") +
  scale_fill_gradient2(
    low = "#2166AC", mid = "white", high = "#B2182B",
    midpoint = 0,
    name = "Diferença GWR − SAR"
  ) +
  labs(
    title = "🗺️ Diferença entre Coeficientes Locais (GWR) e Global (SAR)",
    subtitle = "Regiões onde GWR supera SAR indicam subestimação global",
    caption = paste("Valor global SAR:", round(coef_sar, 5))
  ) +
  theme_minimal()

hexágonos mais discrepantes

Essas regiões podem ser prioritárias em estudos mais refinados, pois mostram onde os modelos globais “apagam” nuances locais importantes.

🔍 Top 10 Hexágonos onde GWR supera SAR
ID Renda_GWR dif_coef
242 0.00509 0.00445
286 0.00504 0.00439
264 0.00503 0.00438
330 0.00496 0.00431
308 0.00492 0.00428
241 0.00492 0.00428
352 0.00487 0.00422
646 0.00481 0.00417
220 0.00481 0.00416
577 0.00480 0.00415

Análise estatística da assimetria

Valores de skewness > 0 confirmam a cauda longa à direita. Isso reforça que o SAR não captura bem os extremos.

skewness(hex_sf$coef_renda, na.rm = TRUE)
## [1] -0.300166
📊 Implicações Estatísticas e Espaciais da Distribuição GWR vs SAR
🔍 Observação 📐 Implicação Estatística 🗺️ Significado Espacial
Linha vermelha deslocada à esquerda do pico GWR SAR subestima o efeito positivo da renda em boa parte da cidade Hotspots locais são mais intensos do que o SAR detecta
Distribuição GWR com cauda longa à direita Existem locais com forte associação positiva entre renda e frequência Alguns bairros têm padrões muito distintos dos demais
Grande dispersão nos coeficientes GWR Alto desvio padrão indica variação espacial real Recomenda-se análise local (GWR) para capturar os extremos
Coeficientes locais positivos e negativos Sugere presença de heterogeneidade espacial Políticas públicas podem precisar de abordagens regionalizadas
SAR retorna valor médio único Modelo global ignora variações locais importantes SAR é útil para visão geral, mas GWR mostra nuance