Introdução

Este trabalho visa aplicar conceitos de geoestatística na análise espacial de variáveis sociodemográficas integradas com dados geográficos do Brasil. Os dados utilizados são provenientes de fontes públicas, via função geobr::read_*().

Objetivos

Importar e visualizar mapas geográficos com geobr;
Integrar variáveis sociodemográficas;
Aplicar métodos de análise espacial: autocorrelação (Moran), krigagem, regressão espacial;
Visualizar padrões com ggplot2, tmap e leaflet.

Carregamento de Pacotes

# 📊 Visualização Gráfica
library(ggplot2)          # Gráficos em camadas
library(tmap)             # Mapas temáticos
library(leaflet)          # Mapas interativos
library(leaflet.extras2)  # Extensões do leaflet
library(leafem)           # Funcionalidades extra para leaflet
library(ggspatial)        # Escalas e bússolas para ggplot2
library(vioplot)          # Violin plots
library(RColorBrewer)     # Paletas predefinidas
library(colorspace)       # Manipulação/criação de cores

# 🌍 Dados Espaciais e Geoprocessamento
library(sf)               # Simple Features
library(sp)               # Estrutura espacial clássica
library(geobr)            # Mapas do Brasil
library(terra)
library(raster)           # Dados raster
library(lattice)          # Gráficos multivariados
library(sidrar)
library(dplyr)
library(rnaturalearth)
library(devtools)

# 📈 Estatística Espacial e Modelagem
library(spatstat)         # Análise de padrões espaciais 
library(spdep)            # Autocorrelação espacial
library(spatialreg)       # Regressão espacial
library(spgwr)            # GWR (Geographically Weighted Regression)
library(gstat)            # Krigagem e geoestatística
library(automap)          # Krigagem automatizada
library(tmaptools)
#-------------------------------------------------------------------#

Construção do mapa: Brasil

De forma ilustrativa, montaremos o mapa do Brasil para incorporar os indicadores posteriormente.

brasil.ufs <- read_state(code_state = "all", year = 2019, showProgress = FALSE)

ggplot(brasil.ufs) + geom_sf()

ggplot(brasil.ufs) + geom_sf(aes(fill = name_region)) + labs(fill = "Região")

Integração com Dados Socioeconômicos da população (Censo 2010)

Aqui geramos um mapa colorimétrico, que tem como objetivo descrever os indicadores populacionais em cada Unidade de Federação.

# Data frame com população por UF - Censo 2010 (limpo e direto)
pop_censo2010 <- data.frame(
  sigla = c("AC", "AL", "AP", "AM", "BA", "CE", "DF", "ES", "GO", "MA", "MT", "MS", 
            "MG", "PA", "PB", "PR", "PE", "PI", "RJ", "RN", "RS", "RO", "RR", "SC", 
            "SP", "SE", "TO"),
  populacao = c(733559, 3245979, 668689, 3516620, 14930634, 8778575, 2570160,
                3516542, 6009975, 6850943, 3184412, 2470140, 19597331, 7603791,
                3766528, 10284552, 8796448, 3119580, 15989929, 3168027,
                10787393, 1749770, 450479, 6248436, 41262199, 2065380, 1383441),
  alfab = c(86.9, 75.4, 90.0, 88.0, 82.2, 84.8, 97.4, 93.8, 93.5, 82.3, 93.7, 93.7,
                    91.4, 85.0, 83.3, 95.2, 89.6, 82.8, 97.3, 88.8, 96.5, 93.3, 96.3, 96.9,
                    97.3, 85.0, 89.6)
)

pop_censo2010["alfab_abs"] <- pop_censo2010["alfab"]*pop_censo2010["populacao"]/100

# Juntar com o mapa
brasil.ufs <- brasil.ufs %>%
  select(abbrev_state, geom) %>%
  left_join(pop_censo2010, by = c("abbrev_state" = "sigla"))

# Verifique se as colunas estão corretas
names(brasil.ufs)

## [1] "abbrev_state" "populacao"    "alfab"        "alfab_abs"    "geom"

# Calcular os centróides dos estados
centroides <- st_centroid(brasil.ufs)

# Mapa
ggplot() +
  # Polígonos dos estados, preenchidos com a população
  geom_sf(data = brasil.ufs, aes(fill = populacao / 1e6), color = "black") +
  
  # Pontos nos centróides, com tamanho proporcional à alfabetização
  geom_sf(data = centroides, aes(size = alfab_abs / 1e6), color = "red", alpha = 0.7) +
  
  scale_fill_viridis_c(
    name = "População (Mi)",
    direction = -1,
    na.value = "grey80"
  ) +
  scale_size_continuous(
    name = "Alfabetização (Mi)",
    range = c(2, 8)
  ) +
  xlab("Longitude") +
  ylab("Latitude") +
  annotation_scale(location = "bl", width_hint = 0.25) +
  annotation_north_arrow(location = "br", which_north = "true") +
  theme_minimal()

A seguir, podemos gerar um mapa baseado em centróides calculados a partir dos polígonos que delimitam os Estados brasileiros.

# Gerar o gráfico
coord_pontos <- brasil.ufs |> st_centroid()

ggplot(brasil.ufs) +
  geom_sf(fill = "gray95", color = "gray30") +
  geom_sf(data = coord_pontos, aes(size = populacao/1e6),
          color = "blue", alpha = 0.65, show.legend = "point") +
  scale_size_continuous(name = "População por UF (Milhões)", range = c(1, 12)) +
  xlab("Longitude") + ylab("Latitude") +
  coord_sf(crs = st_crs(brasil.ufs)) +
  annotation_scale(location = "bl", width_hint = 0.25) +
  annotation_north_arrow(location = "br", which_north = "true", 
                         style = north_arrow_fancy_orienteering()) +
  theme_minimal()

E para mapas interativos, ou seja, gastar uma onda desnecessária (no meu caso):

leaflet(coord_pontos) |>
  addProviderTiles(providers$OpenStreetMap, group = "Mapa Político") |>
  addProviderTiles(providers$Esri.WorldImagery, group = "Satélite") |>
  addCircleMarkers(
    label = ~ paste0(abbrev_state, ": ", round(populacao / 1e7, 2), " Mi"), 
    labelOptions = labelOptions(textsize = "13px"), 
    radius = ~ populacao / 1e6,  # Escala ajustada (em milhões)
    fillOpacity = 0.5,
    color = "blue",
    stroke = FALSE,
    group = "Dados"
  ) |>
  addLayersControl(
    baseGroups = c("Mapa Político", "Satélite"),
    overlayGroups = "Dados",
    options = layersControlOptions(collapsed = FALSE)
  )

Análise Estatística Espacial

I) Formatação dos dados

Agora, iremos ler um shape file contendo a malha das Unidades de Federação brasileiras.

Primeiramente, geramos o contorno da malha e projetamos os dados:

br <- read_sf("C:/Users/Usuario/OneDrive/Doutorado/BR_UF_2024.shp", crs = 31981)

br_joined <- br |>
  mutate(idsetor = as.numeric(CD_REGIA)) |>
  inner_join(pop_censo2010, by = c("SIGLA_UF" = "sigla"))

# Unir todos os estados em um único polígono (o "contorno do Brasil")
br_union <- st_union(br_joined)

# Visualizar o contorno
plot(st_geometry(br_union), main = "Contorno do Brasil")

Aqui, deveremos tratar as classes até gerarmos a classe “owin”, pertinente para modelagem espacial. Além disso, os centróides devem ser gerados.

br.owin <- as.owin(st_geometry(br_union))

centroides <- st_centroid(st_geometry(br_joined))

# Transformando em os centróides em formato sp (spatstat)
centroides.sp <- as.data.frame(as_Spatial(centroides))
names(centroides.sp) <- c("X", "Y")

plot(centroides.sp)

Agora, unificando os centróides ao contorno:

centroides.ppp <- ppp(centroides.sp$X, centroides.sp$Y, br.owin)

plot(centroides.ppp, pch = 20, cex = 0.9)

II) Cálculo de taxas e mapa de densidade (kernel)

Vamos começar gerando ruídos através de uma distribuição uniformemente variada. O intervalo [-15%,+15%] resguarda as características dos dados por UF. Isso evitará viéses na análise espacial.

# Adicionar ruído uniforme entre -15 e +15 da população
set.seed(123)
ruido_1 <- runif(n = nrow(br_joined), min = -15, max = 15)
ruido_2 <- runif(n = nrow(br_joined), min = -0.05, max = 0.05)

br_joined$AREA_KM2 <- abs(br_joined$AREA_KM2 * (1 + ruido_1))
br_joined$alfab <- abs(br_joined$alfab * (1 + ruido_2))
br_joined$populacao <- abs(br_joined$populacao * (1 + ruido_1))

Seguindo, calcularemos a densidade populacional e Nº de alfabetizados por UF, para direcionar o kernel, que preencherá de forma continua o contorno.

# Coordenadas territoriais
br_joined_proj <- st_transform(br_joined, crs = 5880)

# Agora sim, calcule a área corretamente em m² e converta para km²
br_joined$area_km2 <- br_joined$AREA_KM2 


# Densidade populacional
br_joined$tx1 <- br_joined$populacao / br_joined$area_km2

# Calcular alfabetizacao absoluta
br_joined$tx2 <- (br_joined$alfab_abs) * br_joined$populacao

kernel.tx1 <- density(centroides.ppp, 1000, weights = br_joined$tx1, scalekernel = TRUE)
kernel.tx2 <- density(centroides.ppp, 1000, weights = br_joined$tx2, scalekernel = TRUE)

plot(kernel.tx1)

plot(kernel.tx2)

III) Matriz das vizinhanças

Para verificarmos autocorrelação espacial, precisamos construir a matriz das vizinhanças.

viz <- poly2nb(br_joined)

br_sp <- as(br_joined, "Spatial")  # convertendo em formato sp
coord <- coordinates(br_sp)  # coordenadas dos centroides

class(br_sp)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

viz

## Neighbour list object:
## Number of regions: 27 
## Number of nonzero links: 102 
## Percentage nonzero weights: 13.99177 
## Average number of links: 3.777778

Devemos verificar as malhas para geração do gráfico.

viz.sf <- as(nb2lines(viz, coords = coord), "sf")
viz.sf <- st_set_crs(viz.sf, st_crs(br_joined))

mapa.viz <- ggplot(br_joined) + geom_sf(fill = "lightpink", color = "white") + geom_sf(data = viz.sf) +
    theme_minimal() + ggtitle("Vizinhança por \n conectividade") + ylab("Latitude") +
    xlab("Longitude")
mapa.viz

IV) Correlação dos pesos vizinhos

O teste de Moran avalia se existe autocorrelação espacial entre valores de uma variável geográfica — ou seja, se lugares próximos tendem a ter valores parecidos (ou não).

a) Densidade Populacional

pesos.viz <- nb2listw(viz)
moran.test(br_joined$tx1, pesos.viz)

## 
##  Moran I test under randomisation
## 
## data:  br_joined$tx1  
## weights: pesos.viz    
## 
## Moran I statistic standard deviate = 1.1498, p-value = 0.1251
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.10055114       -0.03846154        0.01461601

Para variável “densidade populacional” a hipótese nula é rejeitada a 95% de confiança, sem evidência sobre presença de autocorrelação espacial.

b) Número de alfabetizados

moran.test(br_joined$tx2, pesos.viz)

## 
##  Moran I test under randomisation
## 
## data:  br_joined$tx2  
## weights: pesos.viz    
## 
## Moran I statistic standard deviate = 2.5435, p-value = 0.005487
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.22938465       -0.03846154        0.01108923

Para variável “Número de alfabetizados”, p-value é limiar: há evidências para rejeitarmos H0, o que é um ótimo sinal para modelagem espacial.

V) Correlograma

Pelo correlograma, podemos ver a similaridade com os vizinhos próximos, indicados pelo “Lag”. Espera-se que os vizinhos viscinais apresentem as maiores autocorrelações e com o aumento de Lag essa autocorrelação diminua.

correl1 <- sp.correlogram(viz, 
                          br_joined$tx1, 
                          order = 6, 
                          method = "I", 
                          zero.policy = TRUE)
correl2 <- sp.correlogram(viz, 
                          br_joined$tx2, 
                          order = 6, 
                          method = "I",
                          zero.policy = TRUE)

plot(correl1)

plot(correl2)

VII) LISA

Em cada UF do mapa do Brasil, podemos observar como o p-value dos alfabetizados se distribui no território nacional, como uma forma extremamente intrigante e poderosa de avaliar visualmente a autocorrelação.

Muitas interpretações diretas podem ser feitas com esse tipo de mapa.

br_joined$pval <- localmoran(br_joined$tx2, pesos.viz)[, 5]

tm_shape(br_joined) + tm_polygons(col = "pval", title = "p-valores", breaks = c(0,
    0.01, 0.05, 0.1, 1), style = "fixed", palette = "-Oranges", border.col = "grey")

library(scales)

resI <- localmoran.sad(lm(br_joined$tx2 ~ 1), 1:length(viz), viz, style = "W")
summary(resI)[1:10, ]

br_joined$MoranLocal <- summary(resI)[, 1]
ggplot(br_joined) + geom_sf(aes(fill = MoranLocal), color = "black") + scale_fill_gradientn(colours = c("blue",
    "white", "red"), values = rescale(c(min(br_joined$MoranLocal), 0, max(br_joined$MoranLocal))),
    guide = "colorbar") + ggtitle("Moran local") + theme_void()

Mapas com classificação de clusteres são úteis para identificar verdadeiros positivos e falsos positivos. Por exemplo, “Alto-Alto” e “Baixo-Baixo” são casos de UF positivos e negativos verdadeiros. Casos como “Alto-Baixo” ou “Baixo-Alto” devem ser analisados com cautela: Nem sempre serão falso-positivos.

library(spdep)

# Convertendo uma estrutura de vizinhança (objeto nb) em uma matriz de pesos
# espaciais padronizada
w <- nb2listw(viz, style = "W")

# Matriz de pesos espaciais padronizada (row-standardized)
resI <- localmoran(br_joined$tx2, nb2listw(viz, style = "W"))

# Armazena estatísticas no objeto espacial
br_joined$Ii <- resI[, 1]  # valor de Moran Local
br_joined$pvalue <- resI[, 5]  # p-valor

# Valor padronizado da variável tx
br_joined$z_tx2 <- scale(br_joined$tx2)[, 1]

# Média dos vizinhos (lagged value)
lag_tx2 <- lag.listw(w, br_joined$z_tx2)

# Classificação dos tipos de associação local dos clusters
br_joined$cluster_type <- "Nao significativo"

br_joined$cluster_type[br_joined$pvalue <= 0.05 & br_joined$z_tx2 > 0 & lag_tx2 > 0] <- "Alto-Alto"
br_joined$cluster_type[br_joined$pvalue <= 0.05 & br_joined$z_tx2 < 0 & lag_tx2 < 0] <- "Baixo-Baixo"
br_joined$cluster_type[br_joined$pvalue <= 0.05 & br_joined$z_tx2 > 0 & lag_tx2 < 0] <- "Alto-Baixo"
br_joined$cluster_type[br_joined$pvalue <= 0.05 & br_joined$z_tx2 < 0 & lag_tx2 > 0] <- "Baixo-Alto"

# Modo estático (ideal para relatórios impressos ou PDF)
tmap_mode("plot")

# Tema visual limpo com fundo branco
tmap_style("white")

tm_shape(br_joined) +
  tm_polygons(
    fill = "cluster_type",
    fill.scale = tm_scale_nominal(
      values = c(
        "Alto-Alto"         = "#E60000",  # vermelho forte
        "Baixo-Baixo"       = "#0033CC",  # azul escuro
        "Baixo-Alto"        = "#9999FF",  # azul claro
        "Alto-Baixo"        = "#FF9999",  # rosa claro
        "Nao significativo" = "#FFFFFF"   # branco
      )
    ),
    fill.legend = tm_legend(
      title = "Cluster Local (LISA)",
      orientation = "portrait"
    ),
    border.col = "gray40"
  ) +
  tm_title("Mapa LISA - Clusters Locais de pessoas alfabetizadas") +
  tm_layout(
    legend.outside = TRUE,
    frame = FALSE,
    bg.color = "white"
  )

VIII) Modelos de Regressão Espacial

a) MRLS

Precisaremos normalizar as nossas variáveis, pois as escalas são discrepantes. Iremos explorar a relação do número de alfabetizados com a população por UF.

br_joined$tx2_z    <- scale(br_joined$tx2)
br_joined$tx1_z    <- scale(br_joined$tx1)

br_joined$populacao[is.na(br_joined$populacao)] <- 0
br_joined$pop_z <- scale(br_joined$populacao)

pop_z <- br_joined$pop_z

br.lm <- lm(br_joined$tx2_z ~ br_joined$tx1_z, data = br_joined)

summary(br.lm)

## 
## Call:
## lm(formula = br_joined$tx2_z ~ br_joined$tx1_z, data = br_joined)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4237 -0.2918 -0.2721 -0.0409  3.9374 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)
## (Intercept)     1.709e-16  1.887e-01   0.000    1.000
## br_joined$tx1_z 2.753e-01  1.923e-01   1.432    0.165
## 
## Residual standard error: 0.9804 on 25 degrees of freedom
## Multiple R-squared:  0.0758, Adjusted R-squared:  0.03883 
## F-statistic:  2.05 on 1 and 25 DF,  p-value: 0.1646

br.lm$lmresid <- residuals(br.lm)
moran.test(br.lm$lmresid, pesos.viz)

## 
##  Moran I test under randomisation
## 
## data:  br.lm$lmresid  
## weights: pesos.viz    
## 
## Moran I statistic standard deviate = 1.124, p-value = 0.1305
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.08342923       -0.03846154        0.01176023

b) CAR

br.car <- errorsarlm(br_joined$tx2_z ~ br_joined$tx1_z, data = br_joined, 
                     listw = pesos.viz)

br.car$carresid <- residuals(br.car)

moran.test(br.car$carresid, pesos.viz)

## 
##  Moran I test under randomisation
## 
## data:  br.car$carresid  
## weights: pesos.viz    
## 
## Moran I statistic standard deviate = 0.50023, p-value = 0.3085
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.01421713       -0.03846154        0.01108977

summary(br.car)

## 
## Call:errorsarlm(formula = br_joined$tx2_z ~ br_joined$tx1_z, data = br_joined, 
##     listw = pesos.viz)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.142436 -0.312780 -0.248792 -0.030881  3.945157 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)     -0.043246   0.234351 -0.1845   0.8536
## br_joined$tx1_z  0.151314   0.185182  0.8171   0.4139
## 
## Lambda: 0.23921, LR test value: 0.57344, p-value: 0.4489
## Asymptotic standard error: 0.22155
##     z-value: 1.0797, p-value: 0.28026
## Wald statistic: 1.1658, p-value: 0.28026
## 
## Log likelihood: -36.45098 for error model
## ML residual variance (sigma squared): 0.85701, (sigma: 0.92575)
## Number of observations: 27 
## Number of parameters estimated: 4 
## AIC: 80.902, (AIC for lm: 79.475)

c) GWR

library(spgwr)

# Estimando a largura de banda “ideal” para o kernel
GWRbanda <- gwr.sel(br_joined$tx2_z ~ br_joined$tx1_z, data = br_joined, 
                    coords = cbind(centroides.sp$X, centroides.sp$Y),
    adapt = T)

## Adaptive q: 0.381966 CV score: 35.10261 
## Adaptive q: 0.618034 CV score: 34.69812 
## Adaptive q: 0.763932 CV score: 34.5467 
## Adaptive q: 0.854102 CV score: 34.52642 
## Adaptive q: 0.8416852 CV score: 34.53176 
## Adaptive q: 0.9098301 CV score: 34.42356 
## Adaptive q: 0.9442719 CV score: 34.35732 
## Adaptive q: 0.9655581 CV score: 34.33143 
## Adaptive q: 0.9787138 CV score: 34.32989 
## Adaptive q: 0.9739747 CV score: 34.33053 
## Adaptive q: 0.9868444 CV score: 34.32864 
## Adaptive q: 0.9918694 CV score: 34.32777 
## Adaptive q: 0.994975 CV score: 34.32721 
## Adaptive q: 0.9968944 CV score: 34.32685 
## Adaptive q: 0.9980806 CV score: 34.32663 
## Adaptive q: 0.9988138 CV score: 34.32649 
## Adaptive q: 0.9992669 CV score: 34.3264 
## Adaptive q: 0.9995469 CV score: 34.32635 
## Adaptive q: 0.99972 CV score: 34.32631 
## Adaptive q: 0.9998269 CV score: 34.32629 
## Adaptive q: 0.999893 CV score: 34.32628 
## Adaptive q: 0.9999339 CV score: 34.32627 
## Adaptive q: 0.9999339 CV score: 34.32627

br.gwr = gwr(br_joined$tx2_z ~ br_joined$tx1_z, data = br_joined, 
             coords = cbind(centroides.sp$X, centroides.sp$Y),
    adapt = GWRbanda, hatmatrix = TRUE, se.fit = TRUE)

br.gwr

## Call:
## gwr(formula = br_joined$tx2_z ~ br_joined$tx1_z, data = br_joined, 
##     coords = cbind(centroides.sp$X, centroides.sp$Y), adapt = GWRbanda, 
##     hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss 
## Adaptive quantile: 0.9999339 (about 26 of 27 data points)
## Summary of GWR coefficient estimates at data points:
##                       Min.    1st Qu.     Median    3rd Qu.       Max. Global
## X.Intercept.    -0.0265309 -0.0057276  0.0040564  0.0294579  0.0406189 0.0000
## br_joined.tx1_z  0.2452803  0.2537648  0.2573461  0.2680764  0.2801634 0.2753
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 2.434741 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 24.56526 
## Sigma (residual: 2traceS - traceS'S): 0.9802575 
## Effective number of parameters (model: traceS): 2.227642 
## Effective degrees of freedom (model: traceS): 24.77236 
## Sigma (model: traceS): 0.9761514 
## Sigma (ML): 0.9350158 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 80.64802 
## AIC (GWR p. 96, eq. 4.22): 75.22196 
## Residual sum of squares: 23.60487 
## Quasi-global R2: 0.09212027

results <- as.data.frame(br.gwr$SDF)
head(results)

results$residuos <- br_joined$tx2_z - results$pred
moran.test(results$residuos, pesos.viz)

## 
##  Moran I test under randomisation
## 
## data:  results$residuos  
## weights: pesos.viz    
## 
## Moran I statistic standard deviate = 1.0813, p-value = 0.1398
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.07832296       -0.03846154        0.01166575

br_joined$coef.tx1 <- results$br_joined.tx1_z
br_joined$localR2 <- results$localR2
br_joined$pred.gwr <- results$pred

hist(results$localR2)
abline(v = median(results$localR2), col = "red")

Por fim, faremos um mapa com a distribuição dos coeficientes e de R² local:

map.gwr <- ggplot(br_joined) +
  geom_sf(aes(fill = coef.tx1), color = "black", size = 0.1) +
  scale_fill_gradient2(
    low = "blue",      # para coeficientes negativos
    mid = "white",     # para valores próximos de zero
    high = "pink",      # para coeficientes positivos
    midpoint = 0,      # centraliza a escala no zero
    name = "coef.tx1"
  ) +
  ggtitle("Coef. locais da variável densidade populacional (km²)") +
  theme_void()

map.gwr

Aqui, observamos que os coeficientes significativos no ajuste do modelo GWR (apresentando baixo Quasi-R² ~ 9%), pertinentes ao mapa anterior. Basicamentem, o modelo não é capaz de mapear a alfabetização em função da densidade populacional.

# Calculando a estatística de wald
br_joined$wald.teste <- abs(results$br_joined.tx1_z/results$br_joined.tx1_z_se)

# Dicotomizando a estatística de wald
br_joined$wald.teste <- ifelse(br_joined$wald.teste < 2, 0, 1)

map.wald <- ggplot(br_joined) + geom_sf(aes(fill = factor(wald.teste)), color = "black") +
    scale_fill_manual(values = c("white", "purple"), labels = c("< 2", ">= 2"), name = "Wald") +
    ggtitle("Coef. significativos (Teste de Wald)") + theme_void()


map.wald

Nenhum coef. foi significativo. Por fim, o mapa da distribuição dos coeficientes de determinação (R²).

ggplot(br_joined) + geom_sf(aes(fill = localR2), color = "black", size = 0.1) + scale_fill_viridis_c(option = "viridis",
    , direction = -1, name = expression(R^2 ~ local)) + ggtitle("Mapa dos coeficientes de determinação regionalizados (R² local)") +
    theme_void()

Conclusão

a) Dados da população e da proporção de alfabetizados - CENSO IBGE 2010

Pela elaboração dos mapas, conseguimos identificar que a maior parte da população se concentra nas áreas mais urbanizadas com maior atividade industrial, como as grandes metrópoles. A proporção de alfabetizados não traz conclusões, e por isso deve ser transformada em um valor absoluto, como o número de pessoas alfabetizadas por UF.

Naturalmente, o número absoluto de alfabetizados acompanhará o número absoluto de pessoas por UF, já que o índice tem seu mínimo em 80% e máximo em 95%.

Porém, ao observarmos o mapa, nota-se estranheza sobre a densidade populacional e o número de pessoas alfabetizadas. Há maior concentração de pessoas em áreas mais urbanizadas não necessariamente extensas em área.

Isso levanta o questionamento: há alguma relação espacial sobre a alfabetização que possa ser explicada pelo número de pessoas por UF? Se sim, isso pode ser um indicador da desigualdade social?

Essa hipótese foi investigada por meio dos modelos RLS, CAR e GWR.

b) Teste de Moran para as variáveis

Com um p-valor = 0.1309, a variável densidade populacional (km²) rejeita H0, não apresentando autocorrelação espacial. Ou seja, outros fatores demográficos torna sua correlação não-linear com o território.

Já a variável Nº de alfabetizados apresenta um p-valor ~ 0.005 e estatística de Moran I ~ 0.22, apresentando correlação positiva razoável. O correlograma mostra o comportamento adequado, com corte em Lag 1.

(o.b.s.: O famoso Leg-press ou dia de perna).

Vale comentar, que os kerneis reforçam a ideia de espalhamento não-homogêneo pelo território, indicando falta de correlação espacial ou correlação defasada.

c) LISA

O Mapa Moran Local evidencia algo interessante. As correlações estatísticamente significativas se concentram no entorno dos maiores polos industriais do Brasil.

Os clusters Alto-Alto estão concentrados nas regiões Sudeste, refletindo o padrão histórico de investimento educacional, maior urbanização e infraestrutura.

Já o cluster Baixo-Alto, normalmente indicando outliers, presente no Centro-Oeste (Mato Grosso do Sul), pode sugerir valores locais inferiores à média regional, apontando possível desigualdade. O efeito de conurbação é visível e seu impacto social pode ser discutido nesse conglomerado.

d) MLSR, CAR e GWR

De maneira geral, a alfabetização não pode ser explicada pela densidade populacional.

Nenhum modelo de regressão espacial foi capaz de captar a relação espacial entre esses indicadores populacionais. O que nos gera alguns insights:

Primeiramente, o MRLS foi desastroso, com R² próximo de zero e um AIC próximo de 79. Os resíduos não apresentaram autocorrelação espacial, porém o modelo em si não ajusta, e o teste de Moran perde significado.
O modelo CAR ajustou tão mal quando o MRLS, com AIC equivalente. Testes LR e Wald mostraram um lambda não significativo.E o teste de coeficiente da variável regressora não rejeitou H0, sendo não significativa.
O modelo GWR também não ajustou, apresentando um Quasi-R² global ~ 0.09. Os coeficientes para cada UF foram próximos entre si, mostrando redundância do modelo, que não captou nenhuma natureza espacial na regressão. Totalmente inconclusivo.
Isso reflete que há uma natureza variacional complexa sobre a alfabetização. Uma ideia futura é utilizar mais indicadores de forma múltipla (dummies) ou não. Exemplos de indicadores: Renda per capta, índice de violência, índices sobre infra-estrutura etc. Principalmente os que refletirem as grandes metrópoles, onde vemos um padrão espacial nas vizinhanças.

Referências

IBGE (2010). Censo Demográfico.

Trabalho Final: Geoestatística Aplicada

Vinícius Nunes da Rocha

2025-06-23