Introdução à análise estatística espacial: da visualização à krigagem

A crescente disponibilidade de dados georreferenciados tem transformado a maneira como investigamos fenômenos em áreas tão diversas como saúde pública, ciências ambientais, agricultura e planejamento urbano. No cerne dessa transformação está a Análise Estatística Espacial, um conjunto de técnicas que permite não apenas visualizar a distribuição geográfica de dados, mas também identificar padrões, detectar dependências espaciais e construir modelos preditivos que consideram explicitamente a localização no espaço.

Este minicurso tem como objetivo servir como um guia introdutório e prático para essa construção, conduzindo o leitor desde os passos fundamentais de construção e visualização de mapas até a aplicação de uma das técnicas de modelagem espacial amplamente utilizada que é a Krigagem. A premissa dessa técnica é a de que a localização não é detalhe, mas uma variável crucial que, quando incorporada à análise, revela o que métodos tradicionais são incapazes de capturar. Padrões de agrupaentos, tendências regionais e a autocorrelação espacial (onde observações próximas tendem a ser mais similares do que observações distantes) são evidências que só se tornam aparentes quando adotamos uma lógica de modelagem estatística espacial.

A escolha do ambiente R como ferramenta para este fim não é aleatória. Sua vasta coleção de pacotes dedicados à análise espacial, como sf para manipulação de dados vetoriais, tmap e ggplot2 para a criação de mapas estáticos e interativos, spdep para análise de autocorrelação e gstat para interpolação e krigagem.Isso,o torna uma plataforma completa e coerente. Este minicurso buscará, portanto, conectar teoria e prática, apresentando os conceitos estatísticos de forma acessível e imediatamente aplicando-os por meio de exemplos reprodutíveis em R.

A estrutura do minicurso reflete um fluxo de análise espacial robusto. Iniciaremos com a visualização, explorando técnicas para criar mapas temáticos que comunicam efetivamente a distribuição espacial de uma variável. Em seguida, avançaremos para a análise exploratória de dados espaciais. E, teremos a etapa de modelagem, com um foco especial na Krigagem. Esta técnica, que vai além de simples interpolações, fornece um modelo estatístico ótimo para predição espacial, fornecendo não apenas uma estimativa para locais não amostrados, mas também uma medida do erro associado a essa estimativa.

Ao final deste minicurso, espera-se que o leitor, seja estudante, pesquisador ou profissional, adquira as bases necessárias para planejar, executar e interpretar suas próprias análises espaciais e se sinta espirado a estudar materiais mais avançados. Mais do que dominar funções em um software, o objetivo é incentivar uma compreensão crítica sobre o papel do espaço na geração de dados, capacitando-o a contribuir de forma mais específica e fundamentada para a sua área de atuação, com evidências robustas para a tomada de decisão.

1.1 Comunicação usando gráficos

  A eficácia da comunicação visual por meio de mapas depende diretamente de sua relevância para o conteúdo. Assim, gráficos ou mapas com informações irrelevantes, mesmo que sejam interessantes, podem distrair o interesse do público e dificultar a compreensão da mensagem central. Diante disso, a utilização desses para a visualização de um tema deve ser feita de forma estratégica, visando comunicar uma informação com clareza e impacto para a comunidade ou o público em geral.

  A análise espacial e o geoprocessamento na vigilância em saúde são ferramentas que integram dados geográficos para mapear, analisar e gerenciar riscos, surtos de doenças e recursos de saúde. Essas técnicas permitem visualizar padrões espaciais e temporais de problemas de saúde, entender a relação entre fatores geográficos e o ambiente com a saúde da população e aprimorar o planejamento de políticas públicas e estratégias de vigilância.

Figura 2. Análise da Situação de Saúde (ASIS): Diagnóstico.
Fonte: Magalhães e Carrijo, 2024- Análise Espacial e Geoprocessamento na Vigilância em Saúde (PROFEPI-MSOPAS).


  O Sistema de Informações Geográficas (SIG) atua como uma plataforma central para integrar, gerenciar e analisar diferentes tipos de dados espaciais, tudo dentro de um contexto específico de uma temática patogênica — como o estudo da disseminação de uma doença. No centro do processo, o SIG é representado por um conjunto de camadas sobrepostas, no qual cada uma corresponde a um tipo diferente de informação geográfica ( uma camada para rios, outra para estradas, uma para casos da doença, outra para densidade populacional). O foco está na capacidade de sobrepor e analisar essas camadas simultaneamente.

Figura 3. Sistema de Informações Geográficas (SIG).
Fonte: Magalhães e Carrijo, 2024- Análise Espacial e Geoprocessamento na Vigilância em Saúde.


  A imagem seguinte exibe o Painel de Controle da COVID-19 (COVID-19 Dashboard), desenvolvido pela Universidade Johns Hopkins. Esta ferramenta foi uma referência mundial para o acompanhamento da pandemia, consolidando dados de diversas fontes globais e exemplificando a aplicação da análise espacial na temática de vigilância epidemiológica.

Figura 4. Painel de Controle da COVID-19 (COVID-19 Dashboard).
Fonte: Centro de Ciência e Engenharia de Sistemas (CSSE) da Universidade Johns Hopkins (JHU).


  O processo de integração e linkage (vinculação) permite que os dados provenientes de diferentes sistemas de informação criem uma única e poderosa Base de Dados Analítica. O objetivo é combinar fontes distintas para permitir análises mais ricas e complexas, que seriam impossíveis de realizar utilizando apenas uma das bases de dados originais. Em especial, existem 4 fontes de dados (input) destaques para gerar a base de Dados Analítica.

  • SIM (Sistema de Informação de Mortalidade): este registra os dados de todas as declarações de óbito do Brasil.

  • IPEA (Instituto de Pesquisa Econômica Aplicada): fornece dados e estudos sobre o contexto socioeconômico com seus indicadores (renda, escolaridade, saneamento, entre outros).

  • TABnet (Dados Populacionais): fornece os denominadores essenciais para o cálculo de taxas e coeficientes (por exemplo, para calcular a taxa de mortalidade, é preciso saber o número de óbitos do SIM e o total da população do TABnet).

  • IBGE ShapeFile: fornece os arquivos geográficos (ShapeFiles) que contêm os mapas digitais do Brasil.

Figura 5. Vinculação da base de dados analítica.
Fonte: Lizzi et. al. Modelagem bayesiana espaço-temporal da taxa de mortalidade feminina no estado Paraná - Brasil - 2023.


  Com essa base de dados unificada, um pesquisador pode, por exemplo, responder a perguntas complexas como: “A taxa de mortalidade por doenças cardíacas (dados do SIM e TABnet) é maior em municípios (mapa do IBGE) com piores indicadores de renda e escolaridade (dados do IPEA)?”

  Outro exemplo que ressalta a importância da criação de mapeamentos para a vigilância em saúde está na série de mapas a seguir. Tais ferramentas permitem identificar rapidamente tanto as tendências temporais quanto as geográficas relacionadas a mortalidade bruta de mulheres no Paraná, a qual apresentou uma tendência de elevação entre 2009 e 2019. Nota-se que ela não se distribui de forma homogênea, concentrando-se de maneira mais crítica em certas regiões do estado, senda essa informação crucial para direcionar investigações mais aprofundadas e o planejamento de políticas públicas de saúde focadas nessas áreas de maior vulnerabilidade.

Figura 6. Recorte temporal da mortalidade bruta das mulheres no Paraná entre 2009-2011 e 2017-2019.
Fonte: Lizzi et. al. Modelagem bayesiana espaço-temporal da taxa de mortalidade feminina no estado Paraná - Brasil - 2023.


  Em paralelo, o mapa detalhado seguinte ilustra os deslocamentos populacionais dentro e para fora do estado do Rio Grande do Sul especificamente quando há a utilização de serviços de saúde de alta complexidade. Complementado por uma legenda abrangente, o mapa oferece insights cruciais sobre a organização da rede de saúde e os padrões de acesso da população.


Figura 7. Deslocamentos para utilização dos serviços de saúde no Rio Grande do Sul.
Fonte: SCHROEDER et. al. Deslocamentos para utilização dos serviços de saúde no Rio Grande do Sul. 2020.


  No R, os mapas são altamente extensíveis e oferecem um vasto leque de personalização, indo muito além das visualizações geográficas básicas. Isso se deve à flexibilidade da plataforma, que permite duas abordagens principais: a primeira é o uso do “Base R”, na qual funções nativas como plot() podem ser usadas para renderizar dados espaciais, sendo uma opção viável tanto para gráficos rápidos e simples ou para exemplos complexos.


2. Mapas de redes (Networks Maps)

  O programa em R a seguir foi criado para visualizar uma rede de dados sobre um mapa, na qual um arquivo nodes.txt, contendo as localizações geográficas (nós), serve de base para que sejam criadas conexões aleatórias (arestas) entre tais localidades, atribuindo a cada uma um peso e uma categoria. Utilizando a biblioteca ggplot2, o script gera um mapa que mostra os nós como pontos e as conexões como linhas. A visualização é enriquecida de duas formas: o tamanho de cada ponto representa o número de conexões que ele possui (calculado com a biblioteca igraph), e a cor de cada linha representa o tipo de categoria da conexão, permitindo uma análise visual da estrutura da rede. O arquivo se encontra comentado para facilitar o entendimento geral sobre o algoritmo.

### Carregamento de bibliotecas
### Descrição das bibliotecas
# assertthat: Pacote para realizar verificações e validações de dados.
# dplyr e purrr: Ferramentas essenciais para manipulação de dados e programação funcional.
# igraph: A principal biblioteca para análise e manipulação de grafos (redes).
# ggplot2, ggraph e ggmap: Bibliotecas poderosas para visualização de dados, redes e mapas.
# maps: Pacote para obter os dados de polígonos do mapa-múndi.

library(assertthat)
library(dplyr)
library(purrr)
library(igraph)
library(ggplot2)
library(ggraph)
library(ggmap)
library(maps) # <-- PACOTE ADICIONADO

#################################################################################
### Leitura do banco de dados
# Carrega os dados dos nós (vértices) a partir de um arquivo de texto.
nodes <- read.table("nodes.txt", header=T)

# Renomeia as colunas para um padrão consistente em minúsculas.
colnames(nodes) <- c('id', 'lon', 'lat', 'name')

attach(nodes) # 'attach' não é recomendado, mas mantido.
# head(nodes) 

##########################################################################################
### Configuração de parâmetros aleatórios
# O arquivo 'nodes.txt' contém os "países" (nós) do gráfico.
# Agora, vamos criar conexões aleatórias (arestas/edges) entre eles.

set.seed(123) # Define a semente aleatória para reprodutibilidade.

N_EDGES_PER_NODE_MIN <- 1
N_EDGES_PER_NODE_MAX <- 4
N_CATEGORIES <- 4

###############################################################################################
### Geração das arestas (edges)
# Itera sobre cada 'id' de nó para criar conexões aleatórias.

edges <- map_dfr(nodes$id, function(id) {
  # 1. Define um número aleatório 'n' de conexões
  n <- floor(runif(1, N_EDGES_PER_NODE_MIN, N_EDGES_PER_NODE_MAX + 1))
  
  # 2. Sorteia 'n' IDs de destino
  to <- sample(1:max(nodes$id), n, replace = FALSE)
  
  # 3. Remove "auto-loops" (conexões do nó com ele mesmo)
  to <- to[to != id]
  
  # 4. Sorteia categorias e pesos
  categories <- sample(1:N_CATEGORIES, length(to), replace = TRUE)
  weights <- runif(length(to))
  
  # 5. Retorna um data.frame (usando tibble, que é o moderno)
  # (map_dfr irá juntar todos os data.frames no final)
  tibble(from = id, to = to, weight = weights, category = categories) # <-- FUNÇÃO ATUALIZADA
})

# Converte 'category' para fator (ideal para legendas no ggplot2)
edges <- edges %>% mutate(category = as.factor(category))

###############################################################################################
### Criação do objeto grafo (igraph)
# Constrói um grafo não-direcionado com os nós e arestas.
g <- graph_from_data_frame(edges, directed = FALSE, vertices = nodes)

###############################################################################################
### Preparação de dados para o plot (ggplot2)
# Adiciona coordenadas (x, y, xend, yend) à tabela de arestas para o ggplot2.

edges_for_plot <- edges %>%
  # Junta as coordenadas (lon, lat) do nó de ORIGEM (from)
  inner_join(nodes %>% select(id, lon, lat), by = c('from' = 'id')) %>%
  rename(x = lon, y = lat) %>%
  
  # Junta as coordenadas (lon, lat) do nó de DESTINO (to)
  inner_join(nodes %>% select(id, lon, lat), by = c('to' = 'id')) %>%
  rename(xend = lon, yend = lat)

# Verificação de segurança: garante que nenhuma aresta foi perdida.
assert_that(nrow(edges_for_plot) == nrow(edges))

###############################################################################################
### Adição de pesos aos nós (vértices)
# Calcula o "grau" (número de conexões) de cada nó.
nodes$weight = degree(g)

###############################################################################################
### Preparação do Mapa-Múndi (NOVA ETAPA)
# Carrega os dados de polígonos do mapa-múndi do pacote 'maps'
world_map <- map_data("world")

###############################################################################################
### Visualização: Criação do mapa com ggplot2

mapa_plot <- ggplot() +
  
  # 1. Desenha o mapa-múndi como a camada de fundo (NOVA CAMADA)
  # 'group = group' é essencial para o geom_polygon conectar os pontos corretamente.
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "#f0f0f0", color = "dark grey", linewidth = 0.2) +
  
  # 2. Desenha as conexões (arestas) sobre o mapa
  geom_segment(data = edges_for_plot, 
               aes(x = x, y = y, xend = xend, yend = yend, color = category),
               alpha = 0.5) +
  
  # 3. Desenha os países (nós) sobre o mapa
  geom_point(data = nodes, 
             aes(x = lon, y = lat, size = weight), 
             color = "black", fill = "skyblue", shape = 21, alpha = 0.8) +
  
  # 4. Adiciona os rótulos de texto (nomes dos países)
  geom_text(data = nodes, aes(x = lon, y = lat, label = name), 
            vjust = -1.5, size = 3, check_overlap = TRUE) +
  
  # 5. Corrige a projeção do mapa (IMPORTANTE)
  # Isso garante que a proporção lon/lat esteja correta, evitando distorções.
  coord_quickmap() +
  
  # 6. Usa tema limpo (sem eixos ou grades)
  theme_void() +
  
  # 7. Define os títulos e legendas
  labs(title = "Mapa de Conexoes Aleatorias",
       subtitle = "O tamanho do no representa seu numero de conexoes (grau)",
       size = "Num. de Conexoes",
       color = "Categoria")

# Exibe o mapa na aba 'Plots'
# Esta é a única saída que será mostrada no documento final.
print(mapa_plot)


3. Mapas do Brasil usando o pacote GeoBR

  Este script demonstra as capacidades do pacote geobr para análise espacial em R. O código mostra como baixar, com uma única função, diversos tipos de dados geoespaciais oficiais do Brasil, incluindo limites de estados, municípios, e setores censitários. Utilizando os pacotes sf e ggplot2, o script primeiro renderiza mapas simples dessas geometrias. Em seguida, ele avança para técnicas de visualização mais complexas, demonstrando como: 1) criar um mapa coroplético (temático) ao juntar dados externos (como expectativa de vida) à geometria dos estados; 2) gerar um mapa de processo pontual (ou “bubble map”) ao sobrepor dados simulados a um mapa base; e 3) integrar-se a outros pacotes, como o censobr, para processar microdados do censo e mapear indicadores sociais, como a cobertura de rede de esgoto, em nível municipal.

################################
### instalação
# Instala o pacote 'geobr', que dá acesso direto aos dados
# geoespaciais oficiais do Brasil (IBGE, IPEA, etc.).
#install.packages("geobr")


#################################
### Carregando pacotes necessários
library(geobr)   # Para baixar os dados de mapas (shapefiles)
library(sf)      # O pacote principal para manipular dados espaciais (ler, escrever, juntar)
library(dplyr)   # Para manipulação de dados (joins, filtros, sumarizações)
library(ggplot2) # Para visualização de dados e criação dos mapas


#################################
### bancos de dados disponíveis
# A função list_geobr() mostra uma tabela com TODOS os datasets
# que o pacote pode baixar, incluindo os anos disponíveis e a cobertura geográfica.
datasets <- list_geobr()
head(datasets) # Mostra os primeiros 6 datasets disponíveis

#################################
# Exemplo 1: O estado do Sergipe
# read_state() baixa os polígonos (shapefiles) dos estados.
# Ao especificar 'code_state="SE"', baixamos apenas o contorno de Sergipe.
state <- read_state(
  code_state="SE",
  year=2018,
  showProgress = FALSE
)

# Plotando o mapa:
ggplot() +  
  # geom_sf() é a função do ggplot2 para desenhar dados espaciais (objetos 'sf').
  geom_sf(data = state, color=NA, fill = 'lightgray') + # color=NA remove a borda
  theme_void() # theme_void() remove todos os eixos, títulos e fundo cinza.

#################################
# Exemplo 2: O município de São Paulo
# read_municipality() baixa os polígonos de municípios.
# Usamos o código IBGE de 7 dígitos para um município específico.
muni <- read_municipality(
  code_muni = 3550308,  # Este é o código IBGE para São Paulo
  year=2010, 
  showProgress = FALSE
)

ggplot() +  
  geom_sf(data = muni, color=NA, fill = '#1ba185') +
  theme_void()

################################# 
# Exemplo 3: Todos os municípios de Minas Gerais
# Um truque útil: ao passar a sigla de um estado (ex: "MG") para a função
# read_municipality(), o geobr baixa TODOS os municípios daquele estado.
muni <- read_municipality(code_muni = "MG", 
                          year = 2007,
                          showProgress = FALSE)

ggplot() +  
  # Desta vez, definimos 'color="black"' para desenhar as bordas entre os municípios.
  geom_sf(data = muni, color="black", fill = 'gray60', size = 0.1) + # size = 0.1 linhas finas
  theme_void()

################################# 
# Exemplo 4: Todos os setores censitários do RJ
# Setores censitários são as menores unidades geográficas do Censo IBGE.
# Note que isso pode ser um download bem grande e demorado.
cntr <- read_census_tract(
  code_tract = "RJ", 
  year = 2010,
  showProgress = FALSE
)

ggplot() +  
  geom_sf(data = cntr, color="black", fill = 'gray60', size = 0.05) + # size bem pequeno
  theme_void()

################################# 
# Exemplo 5: Regiões Intermediárias (substitutas das Mesorregiões)
# Baixa as divisões de regiões intermediárias do IBGE.
# Como não especificamos um estado, ele baixa para o Brasil inteiro.
inter <- read_intermediate_region(
  year = 2017,
  showProgress = FALSE
)

ggplot() +  
  geom_sf(data = inter, color="yellow", fill = 'blue') +
  theme_void()

################################# 
# Exemplo 6: Todos os estados brasileiros
# Ao não passar nenhum argumento para 'code_state' em read_state(),
# o geobr baixa o shapefile do Brasil com todos os estados.
states <- read_state(
  year = 2019, 
  showProgress = FALSE
)

ggplot() +  
  geom_sf(data = states, color="yellow", fill = 'blue') +
  theme_void()

# Criando um tema customizado para remover todos os eixos e bordas.
# Isso é uma alternativa ao theme_void() e pode ser reutilizado.
no_axis <- theme(axis.title=element_blank(),
                 axis.text=element_blank(),
                 axis.ticks=element_blank())

# Gráfico customizado do Brasil por estados
ggplot() +
  # Usamos cores customizadas (fill e color) e uma borda fina (size=.15)
  geom_sf(data=states, fill="#2D3E50", color="#FEBF57", size=.15, show.legend = FALSE) +
  labs(subtitle="Estados do Brasil, 2019", size=8) +
  theme_minimal() + # theme_minimal() é um tema limpo, mas que mantém o fundo branco
  no_axis # Adicionamos nossa regra 'no_axis' para remover o que restou.

# Exemplo 7: Mapa customizado dos municípios do Rio
all_muni <- read_municipality(
  code_muni = "RJ", 
  year= 2010,
  showProgress = FALSE
)

# Plotando com o mesmo estilo do mapa do Brasil
ggplot() +
  geom_sf(data=all_muni, fill="#2D3E50", color="#FEBF57", size=.15, show.legend = FALSE) +
  labs(subtitle="Municípios do Rio de Janeiro, 2010", size=10) +
  theme_minimal() +
  no_axis

#################################  
# MAPA TEMÁTICO (COROPLÉTICO)
# O objetivo é juntar dados (ex: expectativa de vida) ao nosso mapa (shapefile).

# 1. Ler um data.frame comum (não-espacial) com os dados.
# Este CSV de exemplo vem dentro do pacote 'geobr'.
df <- utils::read.csv(system.file("extdata/br_states_lifexpect2017.csv", package = "geobr"), encoding = "UTF-8")

# 2. Preparar as "chaves" para a junção (join).
# Para juntar 'states' (mapa) com 'df' (dados), precisamos de uma coluna em comum.
# Vamos usar o nome do estado, mas eles precisam estar idênticos.
states$name_state <- tolower(states$name_state) # "Rio de Janeiro" -> "rio de janeiro"
df$uf <- tolower(df$uf)                         # "Rio de Janeiro" -> "rio de janeiro"

# 3. Unir os dados
# Usamos left_join() do dplyr para adicionar as colunas do 'df' ao 'states'.
# O objeto 'states' continua sendo um objeto 'sf', mas agora com mais colunas.
states <- dplyr::left_join(states, df, by = c("name_state" = "uf"))


################################# 
#### Desenhar mapa temático (Coroplético)
# A grande mudança é no 'aes()' (estética).
ggplot() +
  # Mapeamos a coluna 'ESPVIDA2017' (que acabamos de juntar) ao preenchimento 'fill'.
  # O ggplot irá automaticamente colorir cada estado com base no valor dessa coluna.
  geom_sf(data=states, aes(fill=ESPVIDA2017), color= NA, size=.15) +
  labs(subtitle="Expectativa de vida ao nascer - Brasil, 2017", size=8) +
  # scale_fill_distiller() aplica uma paleta de cores (ex: "Reds")
  # e nos permite configurar a legenda (nome, limites).
  scale_fill_distiller(palette = "Reds", name="Expectativa de vida", direction = 1, limits = c(70,80)) +
  theme_minimal() +
  no_axis

################################################################################
############ MAPA DE PROCESSO PONTUAL (Bubble Map)
# O objetivo é plotar pontos (ex: casos de uma doença, lojas) sobre um mapa base.

# Carregar pacotes (já carregados, mas é boa prática listá-los)
library(geobr)
library(ggplot2)
library(dplyr)

# 1. Obter o shapefile do mapa base (polígono)
pr_map <- read_state(code_state = "PR", year = 2020)

# 2. Simular dados de pontos.
# Note que este é um data.frame comum, NÃO um objeto 'sf'.
# O importante é ter colunas de longitude (x) e latitude (y).
set.seed(6789) # Para reprodutibilidade
n_points <- 50
data_points <- data.frame(
  lon = runif(n_points, min = -54.6, max = -48),  # Intervalo de longitude do Paraná
  lat = runif(n_points, min = -26.7, max = -22.5), # Intervalo de latitude do Paraná
  value = runif(n_points, min = 10, max = 100) # Um valor para mapear ao tamanho
)

# 3. Plotar em camadas
ggplot() +
  # Camada 1: O mapa base do Paraná (polígono)
  geom_sf(data = pr_map, fill = "white", color = "black", size = 0.3) +
  # Camada 2: Os pontos, plotados por cima do mapa.
  # O ggplot usa as colunas 'lon' e 'lat' para posicionar os pontos.
  geom_point(data = data_points, 
             # Mapeamos a coluna 'value' ao TAMANHO 'size' do ponto.
             aes(x = lon, y = lat, size = value), 
             color = "blue", alpha = 0.7) + 
  # Controla a escala de tamanho (o tamanho mínimo e máximo das bolhas)
  scale_size_continuous(range = c(2, 8), name = "Valor Simulado") +
  labs(title = "Gráfico de Processo Pontual - Paraná",
       x = "Longitude", y = "Latitude") +
  theme_minimal() # Usamos theme_minimal() desta vez para manter os eixos (lon/lat)

################################################################################
###### Usando dados do Censo (Pacote 'censobr')
## Objetivo: Mapear a proporção de domicílios conectados a
## uma rede de esgoto nos municípios brasileiros.

### Carregar os pacotes
library(censobr) # Pacote para baixar microdados do Censo IBGE
library(arrow)   # Pacote para manipular dados grandes em disco (exigido pelo censobr)

# 1. Baixar os microdados de domicílios do censo 2010.
# Isso NÃO carrega os dados na memória RAM. Ele cria um 'Arrow dataset',
# que aponta para os arquivos baixados no disco (muito eficiente).
hs <- read_households(year = 2010, 
                      showProgress = FALSE)

# 2. Processar os dados (usando dplyr + arrow)
# Este pipeline de 'dplyr' é executado "preguiçosamente" (lazy) no disco.
esg <- hs |> 
  collect() |> # 'collect()' força os dados do disco para a memória RAM.
  group_by(code_muni) |>            # (a) Agrupar por município
  
  # (b) Sumarizar os dados
  # V0010 = Peso do domicílio (quantas pessoas ele representa)
  # V0207 = Tipo de esgotamento (Se '1' = Rede geral)
  summarize(
    # Soma o peso dos domicílios que estão conectados à rede
    rede = sum(V0010[which(V0207=='1')]), 
    # Soma o peso total de todos os domicílios
    total = sum(V0010)
  ) |> 
  
  # (c) Calcular a proporção
  mutate(cobertura = rede / total) |> 
  
  # (d) Coletar os resultados finais
  # O 'esg' agora é um data.frame pequeno na memória, com
  # uma linha por município e sua % de cobertura.
  collect() 


# 3. Geometria dos municípios
# Baixar o shapefile de TODOS os municípios do Brasil em 2010
muni_sf <- geobr::read_municipality(year = 2010,
                                    showProgress = FALSE)


# 4. Unir os dados
# Juntamos nossos dados de esgoto ('esg') ao nosso mapa ('muni_sf')
esg_sf <- left_join(muni_sf, esg, by = 'code_muni')


# 5. Desenhar o mapa temático (Coroplético)
ggplot() +
  geom_sf(data = esg_sf, aes(fill = cobertura), color=NA) + # Mapeia a cobertura ao 'fill'
  labs(title = "Proporção de domicílios conectados a redes de esgoto (2010)") +
  # Usamos a paleta "Greys" e 'direction = 1' (escuro = maior cobertura)
  scale_fill_distiller(palette = "Greys", direction = 1, 
                       name='Cobertura de Esgoto', 
                       # scales::percent formata a legenda (0.2 -> 20%)
                       labels = scales::percent) +
  theme_void()


4. Introdução teórica: Modelos gaussianos para dados geoestatísticos

  Os modelos Gaussianos são amplamente utilizados na prática como modelos para dados geoestatísticos, sendo aplicados como modelos empíricos convenientes que podem capturar uma ampla gama de comportamento espacial de acordo com a especificação de sua estrutura de correlação. Historicamente, uma razão muito boa para se concentrar nos modelos Gaussianos era que eles são unicamente tratáveis como modelos para dados dependentes. Com o uso crescente de métodos computacionalmente intensivos e, em particular, de métodos de inferência baseados em simulação, a tratabilidade analítica dos modelos Gaussianos está se tornando uma razão menos convincente para usá-los. No entanto, ainda é conveniente trabalhar dentro de uma classe de modelo padrão em aplicações de rotina. O escopo da classe de modelos Gaussianos pode ser estendido usando uma transformação da variável de resposta original e, com essa flexibilidade extra, o modelo frequentemente fornece um bom ajuste empírico aos dados.

  Além disso, no contexto específico da geoestatística, a suposição Gaussiana é a contrapartida baseada em modelo de alguns métodos de predição geoestatística amplamente utilizados, incluindo a krigagem simples, ordinária e universal.


4.1 Funções de covariância e o variograma

  Um processo espacial Gaussiano, \(\{S(x) : x \in \mathbb{R}^2\}\), é um processo estocástico com a propriedade de que para qualquer coleção de localizações \(x_1, \dots, x_n\) com cada \(x_i \in \mathbb{R}^2\), a distribuição conjunta de \(S = \{S(x_1), \dots, S(x_n)\}\) é Gaussiana multivariada. Qualquer processo desse tipo é completamente especificado por sua função de média, \(\mu(x) = E[S(x)]\) e, sua função de covariância, \(\gamma(x, x') = \text{Cov}\{S(x), S(x')\}\).

  Em qualquer processo desse tipo, considere um conjunto arbitrário de localizações \(x_1, \dots, x_n\), defina \(S = \{S(x_1), \dots, S(x_n)\}\), escreva \(\mu_S\) para o vetor de \(n\) elementos com elementos \(\mu(x_i)\) e \(G\) para a matriz \(n \times n\) com elementos \(G_{ij} = \gamma(x_i, x_j)\). Então, \(S\) segue uma distribuição Gaussiana multivariada com vetor de média \(\mu_S\) e matriz de covariância \(G\). Escrevemos isso como \(S \sim \text{MVN}(\mu_S, G)\).

  Agora, seja \(T = \sum_{i=1}^n a_i S(x_i)\). Então \(T\) é Gaussiano univariado com média \(\mu_T = \sum_{i=1}^n a_i \mu(x_i)\) e variância

\[ \sigma_T^2 = \sum_{i=1}^n \sum_{j=1}^n a_i a_j G_{ij} = a'Ga, \]

  Onde \(a = (a_1, \dots, a_n)\). Deve ser, portanto, o caso de que \(a'Ga \ge 0\). Essa condição, que deve ser válida para todas as escolhas de \(n\), \((x_1, \dots, x_n)\) e \((a_1, \dots, a_n)\), restringe \(G\) a ser uma matriz positiva definida, e a \(\gamma(\cdot, \cdot)\) correspondente a ser uma função positiva definida. Inversamente, qualquer função positiva definida \(\gamma(\cdot, \cdot)\) é uma função de covariância legítima para um processo espacial Gaussiano.

  Um processo espacial Gaussiano é estacionário se \(\mu(x) = \mu\), uma constante para todo \(x\), e \(\gamma(x, x') = \gamma(u)\), onde \(u = x - x'\) i.e., a covariância depende apenas da diferença vetorial entre \(x\) e \(x'\). Adicionalmente, um processo estacionário é isotrópico se \(\gamma(u) = \gamma(\|u\|)\), onde \(\|\cdot\|\) denota a distância Euclidiana i.e., a covariância entre os valores de \(S(x)\) em quaisquer duas localizações depende apenas da distância entre elas. Note que a variância de um processo estacionário é uma constante, \(\sigma^2 = \gamma(0)\). Nós então definimos a função de correlação como sendo \(\rho(u) = \gamma(u)/\sigma^2\). A função de correlação é simétrica em \(u\) i.e., \(\rho(-u) = \rho(u)\). Isso decorre do fato de que para qualquer \(u\), \(\text{Corr}\{S(x), S(x-u)\} = \text{Corr}\{S(x-u), S(u)\} = \text{Corr}\{S(x), S(x+u)\}\), a segunda igualdade seguindo da estacionariedade de \(S(x)\). Portanto, \(\rho(u) = \rho(-u)\). De agora em diante, usaremos \(u\) para significar tanto o vetor \(x - x'\) quanto o escalar \(\|x - x'\|\) de acordo com o contexto. Também usaremos o termo estacionário como uma abreviação para estacionário e isotrópico. Um processo para o qual \(S(x) - \mu(x)\) é estacionário é chamado de covariância estacionária. Processos desse tipo são muito amplamente utilizados na prática como modelos para dados geoestatísticos.

  Todos esses conceitos podem ser vistos com detalhes no livro texto “Model-based geostatistics”, neste livro é introduzido inicialmente o variograma empírico como uma ferramenta para análise exploratória de dados. Consideramos agora o variograma teórico como uma caracterização alternativa da dependência de segunda ordem em um processo estocástico espacial. O variograma de um processo estocástico espacial \(S(x)\) é a função

\[ V(x, x') = \frac{1}{2}\operatorname{Var}\{S(x) - S(x')\}. \]

  Note que \(V(x, x') = \frac{1}{2}[\text{Var}\{S(x)\} + \text{Var}\{S(x')\} - 2\text{Cov}\{S(x), S(x')\}]\). No caso estacionário, isso simplifica para \(V(u) = \sigma^2\{1 - \rho(u)\}\) o que, incidentalmente, explica por que o fator de um meio é convencionalmente incluído na definição do variograma. O variograma também é bem definido como uma função de \(u\) para uma classe limitada de processos não estacionários; um exemplo unidimensional é um passeio aleatório simples, para o qual \(V(u) = \alpha u\). Processos que são não estacionários, mas para os quais \(V(u)\) é bem definida são chamados de funções aleatórias intrínsecas (Matheron, 1973).

  No caso estacionário, o variograma é teoricamente equivalente à função de covariância, mas tem várias vantagens como ferramenta de análise de dados, especialmente quando as localizações dos dados formam um arranjo irregular. As condições para a validade teórica de uma classe especificada de variogramas são geralmente discutidas em termos da família correspondente de funções de covariância. Gneiting, Sasvári e Schlather (2001) apresentam resultados análogos em termos de variogramas.


4.2 Regularização

  A medição geoestatística pode afetar a escolha de um modelo para os dados. Quando o suporte para cada valor medido se estende por uma área, em vez de estar confinado a um único ponto, o sinal modelado \(S(x)\) deve estritamente ser representado como

\[ S(x) = \int w(r)S^*(x - r)dr, \]

  Onde \(S^*(\cdot)\) é um processo de sinal subjacente e não observado e, \(w(\cdot)\) é uma função de ponderação. Nesse caso, a forma de \(w(\cdot)\) restringe a forma admissível para a função de covariância de \(S(\cdot)\). Especificamente, se \(\gamma(\cdot)\) e \(\gamma^*(\cdot)\) são as funções de covariância de \(S(\cdot)\) e \(S^*(\cdot)\), respectivamente, segue de que

\[ \gamma(u) = \int \int w(r)w(s)\gamma^*(u + r - s)drds. \]

  Agora, faça uma mudança de variável na primeira expressão de \(s\) para \(t = r - s\), e defina

\[ W(t) = \int w(r)w(t - r)dr. \]

  Então, torna-se

\[ \gamma(u) = \int W(t)\gamma^*(u + t)dt. \]

  Funções de ponderação \(w(r)\) típicas seriam radialmente simétricas, de valor não negativo e funções não crescentes de \(\|r\|\); isso vale para o efeito da integração da câmara gama no Exemplo 1.3 do livro referenciado, onde \(w(r)\) não é conhecido explicitamente, mas é suavemente decrescente em \(\|r\|\). Em geral, o efeito de funções de ponderação desse tipo é fazer com que \(S(x)\) varie mais suavemente do que \(S^*(x)\), com um efeito similar em \(\gamma(u)\) em comparação com \(\gamma^*(u)\).

  Um resultado análogo vale para a relação entre os variogramas de \(S(\cdot)\) e \(S^*(\cdot)\). Usando a relação \(V(u) = \gamma(0) - \gamma(u)\), segue

\[ V(u) = \int W(t)\{V^*(t + u) - V^*(t)\}dt. \]

  Se a forma da função de ponderação \(w(\cdot)\) for conhecida, seria possível incorporá-la ao nosso modelo para os dados. Isso significaria especificar um modelo para a função de covariância de \(S^*(\cdot)\) e avaliar tal equação para derivar a função de covariância correspondente de \(S(\cdot)\). Note que isso permitiria que dados com diferentes suportes fossem combinados naturalmente. Uma estratégia pragmática, e a única disponível se \(w(\cdot)\) for desconhecido, é especificar diretamente um modelo apropriadamente suave para a função de covariância de \(S(\cdot)\).

  A questão da regularização também pode surgir em conexão com a predição, em vez da formulação do modelo. O problema de predição geoestatística simples é mapear o sinal espacial \(S(x)\), mas em algumas aplicações um alvo mais relevante para a predição pode ser um mapa de um sinal regularizado, na qual, a integral é sobre um disco com centro \(x\), i.e., \(T(x)\) é uma média espacial sobre o disco.

\[ T(x) = \int S(u)du, \]


4.3 Continuidade e diferenciabilidade de processos estocásticos

  A especificação da estrutura de covariância de um processo espacial \(S(x)\) afeta diretamente a suavidade das superfícies que o processo gera. Descritores matemáticos aceitos da suavidade de uma superfície são a sua continuidade e diferenciabilidade. No entanto, para superfícies geradas estocasticamente \(S(x)\) precisamos distinguir dois tipos de continuidade ou diferenciabilidade. No que se segue, consideraremos um espaço unidimensional \(x\), essencialmente por conveniência notacional.

  Primeiro, consideramos propriedades de média quadrática, definidas como se segue. Um processo estocástico \(S(x)\) é contínuo em média quadrática se \(E[\{S(x+h) - S(x)\}^2] \to 0\) quando \(h \to 0\). Além disso, \(S(x)\) é diferenciável em média quadrática, com derivada em média quadrática \(S'(x)\), se

\[ E\left[\left\{\frac{S(x + h) - S(x)}{h} - S'(x)\right\}^2\right] \to 0 \]

  quando \(h \to 0\). A diferenciabilidade em média quadrática de ordem superior é então definida sequencialmente da maneira óbvia; \(S(x)\) é duas vezes diferenciável em média quadrática se \(S'(x)\) for diferenciável em média quadrática, e assim por diante. Um resultado importante, descrito por exemplo em Bartlett (1955), é o seguinte.

  Teorema: Um processo estocástico estacionário com função de correlação \(\rho(u)\) é \(k\) vezes diferenciável em média quadrática se e somente se \(\rho(u)\) for \(2k\) vezes diferenciável em \(u = 0\).

  Para examinar a diferenciabilidade na origem de qualquer função de correlação particular \(\rho(u)\), precisamos considerar a forma estendida de \(\rho(u)\) na qual \(u\) pode assumir valores positivos ou argumentos negativos com \(\rho(-u) = \rho(u)\). Consequentemente, por exemplo, a função de correlação exponencial \(\rho(u) = \exp(-u/\phi)\) é contínua, mas não diferenciável na origem. Em contraste, a função de correlação Gaussiana, definida por \(\rho(u) = \exp\{-(u/\phi)^2\}\), é infinitamente diferenciável.

Figura 8. Realização de um processo estocástico de valores binários, contínuo em média quadrática.


  Uma segunda versão das propriedades de continuidade e diferenciabilidade diz respeito à continuidade de trajetória e diferenciabilidade de trajetória. Um processo \(S(x)\) é contínuo em trajetória, ou mais geralmente \(k\) vezes diferenciável em trajetória, se as suas realizações são funções contínuas ou \(k\) vezes diferenciáveis, respetivamente.

  Em geral, não precisa haver ligação entre as propriedades de média quadrática e as propriedades de trajetória dos processos estocásticos. Como um exemplo simples, podemos considerar um processo de valores binários \(S(x)\) no qual a reta real é particionada numa sequência de intervalos aleatórios, cujos comprimentos são realizações independentes de uma distribuição exponencial de média unitária, o valor de \(S(x)\) dentro de cada intervalo é zero com probabilidade \(p\), um caso contrário, e os valores de \(S(x)\) em intervalos sucessivos são determinados independentemente. A Figura 8 mostra uma realização com \(p = 0.5\). Claramente, este processo não é contínuo em trajetória. No entanto, a sua função de correlação é a exponencial, \(\rho(u) = \exp(-u)\), que é contínua em \(u = 0\), logo \(S(x)\) é contínuo em média quadrática.

  Kent (1989) apresenta uma discussão teórica rigorosa da continuidade de trajetória para processos estacionários, não necessariamente Gaussianos. Escreva \(\rho(u) = p_m(u) + r_m(u)\), onde \(p_m(u)\) é o polinômio de grau \(m\) dado pela expansão em série de Taylor de \(\rho(u)\) em torno de \(u = 0\). Então, uma condição suficiente para a existência de um processo estacionário bidimensional contínuo em trajetória com função de correlação \(\rho(\cdot)\) é que \(\rho(\cdot)\) seja duas vezes continuamente diferenciável e \(|r_2(u)| = O(u^2 / |\log u|^{3+\gamma})\) quando \(u \to 0\), para algum \(\gamma > 0\). Uma condição ligeiramente mais forte, que é mais fácil de verificar na prática, é que \(|r_2(u)| = O(u^{2+\epsilon})\) para algum \(\epsilon > 0\). Para processos Gaussianos estacionários em duas dimensões, uma condição suficiente para a continuidade de trajetória é que \(\rho(0) - \rho(u) = O(1 / |\log u|^{1+\epsilon})\), o que é apenas ligeiramente mais forte do que o requisito para a continuidade em média quadrática, nomeadamente que \(\rho(\cdot)\) seja contínua na origem.

  Isso justifica usar a diferenciabilidade em média quadrática como uma medida conveniente da suavidade de processos Gaussianos estacionários ao considerar sua adequação como modelos empíricos para fenômenos naturais.


4.4 Famílias de funções de covariância e suas propriedades

  Ser positiva definida é a condição necessária e suficiente para que uma família paramétrica de funções defina uma classe legítima de funções de covariância, mas esta não é uma condição fácil de verificar diretamente. Por esta razão, é útil ter disponível uma gama de famílias padrão que são conhecidas por serem positivas definidas, mas que, em outros aspetos, são suficientemente flexíveis para satisfazer as necessidades de aplicações a dados geoestatísticos. Nesta seção, apresentamos os detalhes de várias dessas famílias e delineamos as suas propriedades. A nossa preocupação aqui é com modelos para processos em duas dimensões espaciais. Todas as famílias de covariância que descrevemos são também válidas em uma ou três dimensões. Em geral, uma família de covariância válida em \(\mathbb{R}^d\) não permanece necessariamente válida em mais do que \(d\) dimensões espaciais, mas é automaticamente válida em dimensões menores que \(d\).

4.4.1 Família Materna

  A forma mais comum de comportamento empírico para uma estrutura de covariância estacionária é que a correlação entre \(S(x)\) e \(S(x')\) diminui à medida que a distância \(u = \|x - x'\|\) aumenta. É, portanto, natural procurar modelos cuja estrutura de correlação teórica se comporte dessa maneira. Além disso, podemos esperar que diferentes aplicações possam exibir diferentes graus de suavidade no processo espacial subjacente \(S(x)\).

  A família materna de funções de correlação, nomeada em homenagem a Matérn (1960), atende a ambos os requisitos. É uma família de dois parâmetros,

\[ \rho(u) = \{2^{\kappa-1}\Gamma(\kappa)\}^{-1} (u/\phi)^\kappa \mathcal{K}_\kappa(u/\phi), \]

  Na qual \(K_\kappa(\cdot)\) denota uma função de Bessel modificada de ordem \(\kappa\), \(\phi > 0\) é um parâmetro de escala com as dimensões de distância, e \(\kappa > 0\), chamado de ordem, é um parâmetro de forma que determina a suavidade analítica do processo subjacente \(S(x)\). Especificamente, \(S(x)\) é \(\lceil\kappa - 1\rceil\) vezes diferenciável em média quadrática, onde \(\lceil\kappa\rceil\) denota o menor inteiro maior ou igual a \(\kappa\).

  A figura a seguir mostra a função de correlação materna para cada \(\kappa = 0.5\), \(1.5\) e \(2.5\), correspondendo a processos \(S(x)\) que são contínuos em média quadrática, uma vez diferenciáveis e duas vezes diferenciáveis, respectivamente. No diagrama, os valores de \(\phi\) foram ajustados de modo a dar a todas as três funções o mesmo alcance prático, que definimos aqui como a distância \(u\) na qual a correlação é \(0.05\). Nessa figura, usamos \(u = 0.75\) como o valor do alcance prático. Para \(\kappa = 0.5\), a função de correlação de materna reduz-se à exponencial, \(\rho(u) = \exp(-u/\phi)\), enquanto que quando \(\kappa \to \infty\), \(\rho(u) \to \exp\{-(u/\phi)^2\}\), que também é chamada de função de correlação Gaussiana ou, de forma algo confusa no presente contexto, o modelo Gaussiano. Whittle (1954) propôs o caso especial da função de correlação materna com \(\kappa = 1\).

Figura 9. Funções de correlação materna com \(\kappa = 0.5\) (linha sólida), \(\kappa = 1.5\) (linha tracejada)
e \(\kappa = 2.5\) (linha pontilhada), e valores ajustados de \(\phi\) para alcances práticos equivalentes.


  Note que os parâmetros \(\phi\) e \(\kappa\) na equação são não-ortogonais, no seguinte sentido. Se a verdadeira estrutura de correlação é materna com parâmetros \(\phi\) e \(\kappa\), então a aproximação de melhor ajuste com ordem \(\kappa^* \ne \kappa\) também terá \(\phi^* \ne \phi\). Em outras palavras, parâmetros de escala correspondentes a diferentes ordens da correlação materna não são diretamente comparáveis. A relação entre o alcance prático e o parâmetro de escala \(\phi\) depende, portanto, do valor de \(\kappa\). Por exemplo, o alcance prático como definido acima é aproximadamente \(3\phi\), \(4.75\phi\) e \(5.92\phi\) para as funções maternas com \(\kappa = 0.5\), \(1.5\) e \(2.5\), respetivamente, e \(\sqrt{3}\phi\) para a função de correlação Gaussiana. Por esta razão, Handcock e Wallis (1994) sugerem uma reparametrização da equação de \(\kappa\) e \(\phi\) para um par quase ortogonal \(\kappa\) e \(\alpha = 2\phi\sqrt{\kappa}\).

  A Figura 10 mostra um traço unidimensional através de uma realização simulada de um processo espacial Gaussiano com cada uma das funções de correlação materna acima, usando a mesma semente aleatória para todas as três realizações. A crescente suavidade analítica do processo à medida que \(\kappa\) aumenta reflete-se na aparência visual das três realizações, mas a diferença mais notável é entre o caso não-diferenciável e o diferenciável, i.e., entre \(\kappa = 0.5\) por um lado e \(\kappa = 1.5\) ou \(\kappa = 2.5\) por outro.

  A Figura 11 mostra realizações bidimensionais simuladas de processos Gaussianos cujas funções de correlação são maternas com \(\kappa = 0.5\) e \(\kappa = 2.5\), novamente usando a mesma semente de número aleatório para tornar as realizações diretamente comparáveis. A diferença na suavidade entre os casos não-diferenciável e diferenciável é novamente visualmente impressionante.

Figura 10. Realizações unidimensionais de processos espaciais Gaussianos cujas funções de correlação são
maternas com \(\kappa = 0.5\) (linha sólida), \(\kappa = 1.5\) (linha tracejada) e \(\kappa = 2.5\) (linha pontilhada).


Figura 11. Simulações de processos Gaussianos com funções de correlação materna
com \(\kappa = 0.5\) e \(\phi = 0.25\) (esquerda) e \(\kappa = 2.5\) e \(\phi = 0.13\) (direita).


4.4.2 Família Exponencial Potencializada

  Esta família é definida pela função de correlação

\[ \rho(u) = \exp\left\{-\left(\frac{u}{\phi}\right)^\kappa\right\} \]

  Assim como a família materna, ela possui um parâmetro de escala \(\phi > 0\), um parâmetro de forma \(\kappa\), neste caso limitado por \(0 < \kappa \le 2\) e, gera funções de correlação que são monotonicamente decrescentes em \(u\). Também como na família materna, a relação entre o alcance prático e o parâmetro \(\phi\) dependerá do valor de \(\kappa\). No entanto, a família é menos flexível, no sentido de que o processo Gaussiano subjacente \(S(x)\) é contínuo em média quadrática e não diferenciável em média quadrática para todo \(0 < \kappa < 2\), mas infinitamente diferenciável em média quadrática quando \(\kappa = 2\), o valor legítimo máximo.

  A figura a seguir mostra a função de correlação exponencial potencializada para cada valor de \(\kappa = 0.7, 1\) e \(2\), e com valores de \(\phi\) ajustados para fornecer o mesmo alcance prático de 0.75.

Figura 12. Funções de correlação exponencial potencializada com \(\kappa = 0.7\) (linha tracejada), \(\kappa = 1\)
(linha contínua) e \(\kappa = 2\) (linha pontilhada) e valores de \(\phi\) ajustados de tal forma que o alcance prático é 0.75.


  A realização para o modelo exponencial potencializado com \(\kappa = 1\) é, portanto, a mesma que para o modelo materno com \(\kappa = 0.5\). Note que as realizações para \(\kappa = 0.7\) e \(\kappa = 1\), ambas correspondendo a processos contínuos em média quadrática, mas não diferenciáveis, parecem bastante similares em seu caráter.

  O caso extremo \(\kappa = 2\), que é equivalente ao caso limite de uma função de correlação materna quando \(\kappa \to \infty\), pode gerar uma estrutura de covariância muito mal condicionada. Um processo \(S(x)\) com esta função de correlação tem a propriedade teórica de que sua realização em um intervalo contínuo arbitrariamente pequeno determina a realização em toda a linha real. Para a maioria das aplicações, isso seria considerado irrealista.


5. Aplicação computacional: Modelos de krigagem espacial (geoestatística)

  Este script em R executa um fluxo de trabalho completo de interpolação geoestatística, conhecida como Krigagem. O processo inicia definindo uma área de estudo ao baixar o mapa do Paraná (com geobr) e, em seguida, simulando 50 pontos de dados amostrais com valores aleatórios dentro de seus limites. Utilizando o pacote geoR, o script modela a correlação espacial desses pontos (através de um semivariograma) para entender como a proximidade influencia os valores. O destaque do código é o uso do pacote sf para gerar uma grade de predição de alta resolução que se encaixa perfeitamente dentro dos contornos do estado, garantindo que a interpolação respeite os limites geográficos. Finalmente, a Krigagem é executada para estimar os valores em toda essa grade, e o ggplot2 é usado para visualizar o resultado como um mapa de calor contínuo, sobreposto pelos pontos de amostra originais.

################################
################################
### Carregando pacotes necessários
library(geobr)    # Para obter shapes do Brasil
library(sf)       # Para trabalhar com dados espaciais
library(dplyr)    # Para manipulação de dados
library(ggplot2)  # Para visualização
library(geoR)     # Para análise geoestatística

#########################
#########################
##### Processo Pontual ##
#########################

# PASSO 1: OBTER O MAPA BASE
# Buscamos o shapefile do Paraná para usar como área de estudo (pacoteGeobr)
pr_map <- read_state(code_state = "PR", year = 2020)

# Visualização rápida do mapa
ggplot() +
  geom_sf(data = pr_map, fill = "lightblue", color = "black", size = 0.5) +
  labs(title = "Mapa Base - Estado do Paraná",
       x = "Longitude", y = "Latitude")

# PASSO 2: SIMULAR DADOS ESPACIAIS
# Criamos pontos aleatórios com valores associados para simular dados reais
set.seed(6789) # Define semente para reproduzir os mesmos resultados

n_points <- 50
data_points <- data.frame(
  lon = runif(n_points, min = -54, max = -50),   # Longitude dentro do PR
  lat = runif(n_points, min = -26, max = -23),   # Latitude dentro do PR
  value = runif(n_points, min = 10, max = 100)   # Valores simulados
)

# PASSO 3: VISUALIZAR PROCESSO PONTUAL
# Mostra os pontos no mapa, com tamanho proporcional ao valor
ggplot() +
  geom_sf(data = pr_map, fill = "white", color = "black", size = 0.3) +
  geom_point(data = data_points, 
             aes(x = lon, y = lat, size = value), 
             color = "blue", alpha = 0.7) +
  scale_size_continuous(range = c(2, 8), name = "Valor") +
  labs(title = "Processo Pontual - Dados Simulados no Paraná",
       x = "Longitude", y = "Latitude") +
  theme_minimal()

################################################################################
### ANÁLISE GEOESTATÍSTICA - VARIOGRAMA E KRIGAGEM- Parte 01
################################################################################

# PASSO 4: CONVERTER PARA OBJETO GEOESTATÍSTICO
# O geoR precisa dos dados no formato 'geodata' para análise
geo_data <- as.geodata(data_points, coords.col = c("lon", "lat"), data.col = "value")

# PASSO 5: CALCULAR O SEMIVARIOGRAMA EXPERIMENTAL
# O variograma medida como a dissimilaridade entre pontos muda com a distância (processo estocástico espacial)

vario <- variog(geo_data, 
                max.dist = 1,    # Distância máxima para análise
                bin.cloud = TRUE) # Mostra nuvem de pontos

# Visualizar o variograma experimental
plot(vario, main = "Semivariograma Experimental")
# Pontos: semivariância calculada para cada distância
# Linha: tendência geral da dependência espacial


# PASSO 6: AJUSTAR MODELO TEÓRICO AO VARIOGRAMA
# Encontramos a curva matemática que melhor descreve o padrão espacial

vario_model <- variofit(vario, 
                        cov.model = "linear",  # Tipo de modelo (linear, exponencial, etc.)
                        nugget = 5,            # Efeito pepita (variação em distância zero)
                        ini.cov.pars = c(50, 0.5)) # Valores iniciais para ajuste

summary(vario_model) # Mostra parâmetros do modelo ajustado
lines(vario_model, col = "red") # Adiciona modelo ao gráfico

# ?variofit


# A função `variofit` é o coração da modelagem geoestatística - ela encontra a curva
# matemática que melhor descreve a estrutura espacial dos seus dados.

vario_model <- variofit(vario,               # Dados experimentais do variograma
                        cov.model = "linear",  # Modelo teórico (linear, exponencial, esférico, etc.)
                        nugget = 5,          # Efeito pepita (variação em distâncias muito curtas)
                        ini.cov.pars = c(50, 0.5)) # Parâmetros iniciais para o ajuste

summary(vario_model) # Mostra parâmetros do modelo ajustado
lines(vario_model, col = "red") # Adiciona modelo ao gráfico

####### PARÂMETROS DO VARIOFIT EXPLICADOS:

# • `cov.model = "linear"`: 
#   - Modelo de covariância linear - crescimento constante sem patamar definido
#   - Fórmula: γ(h) = nugget + psill × h
#   - Use quando a dependência espacial não satura no alcance estudado

# • `nugget = 5`:
#   - Efeito pepita fixo em 5 unidades
#   - Representa variabilidade em escala muito pequena (erro + micro-escala)
#   - Em modelos mais robustos, deixamos `fix.nugget = FALSE` para estimar automaticamente

# • `ini.cov.pars = c(50, 0.5)`:
#   - Valores iniciais para o algoritmo de otimização
#   - c(psill = 50, range = 0.5)
#   - Psill = patamar parcial (variabilidade estruturada espacialmente)
#   - Range = alcance (distância de influência espacial)

####### MODELOS ALTERNATIVOS RECOMENDADOS:

# Para uma análise mais robusta, considere usar:
#
# vario_model_robusto <- variofit(vario,
#     cov.model = "exponential",    # Modelo mais flexível com patamar definido
#     ini.cov.pars = c(50, 0.5),    # Valores iniciais (psill, range)
#     fix.nugget = FALSE,           # DEIXE estimar o nugget automaticamente!
#     nugget = 1,                   # Valor inicial para nugget
#     weights = "cressie",          # Pesos estatisticamente otimizados
#     minimisation.function = "optim") # Algoritmo robusto
#
# VANTAGENS DO MODELO ROBUSTO:
# 1. Exponential modela bem diversos tipos de dependência espacial
# 2. Nugget estimado é mais realista que nugget fixo
# 3. Pesos de Cressie dão menos peso a bins com poucos pares
# 4. Algoritmo 'optim' é menos sensível a valores iniciais



# COMPONENTES DO VARIOGRAMA:
# - Efeito Pepita (Nugget): Variabilidade em distância zero (erro + micro-escala)
# - Patamar (Sill): Variabilidade máxima quando a correlação espacial desaparece  
# - Alcance (Range): Distância onde atinge o patamar (limite da dependência espacial)

################################################################################
### INTERPOLAÇÃO POR KRIGAGEM

# PASSO 7: CRIAR GRADE REGULAR PARA INTERPOLAÇÃO
# Vamos criar pontos onde queremos estimar valores

# Grade retangular simples (para demonstração)
lon_range <- seq(min(data_points$lon), max(data_points$lon), length.out = 100)
lat_range <- seq(min(data_points$lat), max(data_points$lat), length.out = 100)
grid_rect <- expand.grid(lon = lon_range, lat = lat_range)

# PASSO 8: EXECUTAR KRIGAGEM NA GRADE RETANGULAR
# A krigagem usa o modelo do variograma para fazer interpolação ótima
krige_result_rect <- krige.conv(geo_data, 
                                loc = as.matrix(grid_rect), 
                                krige = krige.control(obj.model = vario_model))

grid_rect$value <- krige_result_rect$predict

# PASSO 9: VISUALIZAR MALHA DE INTERPOLAÇÃO INTERMEDIÁRIA
# Mostra a grade retangular com a interpolação
ggplot() +
  geom_tile(data = grid_rect, aes(x = lon, y = lat, fill = value), alpha = 0.7) +
  geom_sf(data = pr_map, fill = NA, color = "black", size = 0.5) +
  scale_fill_viridis_c(name = "Valor Interpolado") +
  geom_point(data = data_points, aes(x = lon, y = lat), color = "red", size = 1.5) +
  labs(title = "Malha de Interpolação - Grade Retangular",
       subtitle = "Pontos vermelhos: dados originais | Cores: valores interpolados",
       x = "Longitude", y = "Latitude") +
  theme_minimal()

################################################################################
### INTERPOLAÇÃO DENTRO DO POLÍGONO

# PASSO 10: CRIAR GRADE DENTRO DOS LIMITES DO ESTADO
# Método mais realista - só interpola dentro da área de interesse

pr_sf <- st_as_sf(pr_map) # Converter para formato sf

# Criar grade regular dentro do polígono
grid_polygon <- st_make_grid(pr_sf, 
                             cellsize = 0.05,  # Tamanho da célula em graus
                             what = "centers") %>% # Criar pontos no centro das células
  st_as_sf() %>% 
  st_filter(pr_sf)  # Manter apenas pontos dentro do polígono

# Extrair coordenadas
grid_coords <- st_coordinates(grid_polygon) %>% as.data.frame()
colnames(grid_coords) <- c("lon", "lat")

# PASSO 11: KRIGAGEM NA GRADE DO POLÍGONO
krige_result_poly <- krige.conv(geo_data, 
                                loc = as.matrix(grid_coords), 
                                krige = krige.control(obj.model = vario_model))

grid_coords$value <- krige_result_poly$predict

# PASSO 12: VISUALIZAÇÃO FINAL
ggplot() +
  geom_sf(data = pr_map, fill = "white", color = "black", size = 0.3) +
  geom_tile(data = grid_coords, aes(x = lon, y = lat, fill = value)) +
  scale_fill_viridis_c(name = "Interpolação") +
  geom_point(data = data_points, aes(x = lon, y = lat, size = value), 
             color = "red", alpha = 0.7) +
  scale_size_continuous(range = c(3, 8), name = "Valor Original") +
  labs(title = "Krigagem com Bolinhas Proporcionais - Paraná",
       subtitle = "Interpolação dentro do polígono com dados originais sobrepostos",
       x = "Longitude", y = "Latitude") +
  theme_minimal()

  Nesta etapa vamos mostrar algumas possibilidades de verificação dos pressupotos desses modelos para qualidade do ajustes.

################################################################################
### VERIFICAÇÃO DO MODELO GEOESTATÍSTICO
################################################################################

# PASSO 13: VALIDAÇÃO CRUZADA (CROSS-VALIDATION)
# A validação cruzada remove cada ponto um por um e estima seu valor usando os demais
################################################################################
### VERIFICAÇÃO DO MODELO GEOESTATÍSTICO - VERSÃO CORRIGIDA
################################################################################

# PASSO 13: VALIDAÇÃO CRUZADA (CROSS-VALIDATION) - CORRIGIDO
cat("=== VALIDAÇÃO CRUZADA DO MODELO ===\n")

# Executar validação cruzada
xvalid <- xvalid(geo_data, model = vario_model)

# Verificar a estrutura do objeto xvalid
cat("Estrutura do objeto xvalid:\n")
print(names(xvalid))
cat("Número de predições:", length(xvalid$predicted), "\n")

# PASSO 14: ANÁLISE DOS RESÍDUOS DA VALIDAÇÃO CRUZADA - CORRIGIDO

# CORREÇÃO: Calcular os resíduos manualmente pois xvalid não tem $residuals
validation_df <- data.frame(
  observed = geo_data$data,           # Valores observados reais
  predicted = xvalid$predicted,       # Valores preditos na validação
  # residuals = observed - predicted (calculamos manualmente)
  std_error = sqrt(xvalid$krige.var)  # Erro padrão (raiz da variância)
)

# CALCULAR RESÍDUOS MANUALMENTE
validation_df$residuals <- validation_df$observed - validation_df$predicted

# Verificar se as dimensões estão corretas
cat("\n=== VERIFICAÇÃO DE DIMENSÕES ===\n")
cat("Número de observações:", nrow(validation_df), "\n")
cat("Primeiras linhas do validation_df:\n")
print(head(validation_df))

# PASSO 15: CÁLCULO DAS MÉTRICAS DE QUALIDADE
MSE <- mean(validation_df$residuals^2)          # Mean Squared Error
RMSE <- sqrt(MSE)                               # Root Mean Squared Error
MAE <- mean(abs(validation_df$residuals))       # Mean Absolute Error
bias <- mean(validation_df$residuals)           # Viés (deve ser próximo de zero)
correlation <- cor(validation_df$observed, validation_df$predicted) # Correlação

cat("\n=== MÉTRICAS DE VALIDAÇÃO ===\n")
cat("RMSE:", round(RMSE, 3), "\n")
cat("MAE:", round(MAE, 3), "\n") 
cat("Viés (Bias):", round(bias, 3), "\n")
cat("Correlação (Obs vs Pred):", round(correlation, 3), "\n")

# PASSO 16: GRÁFICOS DIAGNÓSTICOS

# Gráfico 1: Valores Observados vs Preditos
plot_obs_vs_pred <- ggplot(validation_df, aes(x = observed, y = predicted)) +
  geom_point(alpha = 0.7, color = "blue", size = 2) +
  labs(title = "Validação Cruzada: Observados vs Preditos",
       subtitle = paste("Correlação =", round(correlation, 3), 
                        "| RMSE =", round(RMSE, 3)),
       x = "Valores Observados", 
       y = "Valores Preditos") +
  theme_minimal()

print(plot_obs_vs_pred)

# Gráfico 2: Resíduos vs Valores Preditos
plot_residuals <- ggplot(validation_df, aes(x = predicted, y = residuals)) +
  geom_point(alpha = 0.7, color = "purple", size = 2) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed", size = 1) +
  labs(title = "Análise de Resíduos",
       subtitle = paste("Viés =", round(bias, 3), "| Ideal: resíduos centrados em zero"),
       x = "Valores Preditos", 
       y = "Resíduos (Observado - Predito)") +
  theme_minimal()

print(plot_residuals)

# Gráfico 3: Histograma dos Resíduos
plot_hist_residuals <- ggplot(validation_df, aes(x = residuals)) +
  geom_histogram(aes(y = ..density..), bins = 12, fill = "lightblue", alpha = 0.7) +
  geom_density(color = "darkblue", size = 1) +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed", size = 1) +
  labs(title = "Distribuição dos Resíduos",
       subtitle = "Distribuição deve ser aproximadamente normal centrada em zero",
       x = "Resíduos", 
       y = "Densidade") +
  theme_minimal()

print(plot_hist_residuals)

# PASSO 17: ANÁLISE ESPACIAL DOS RESÍDUOS
spatial_residuals <- data.frame(
  lon = geo_data$coords[,1],
  lat = geo_data$coords[,2],
  residuals = validation_df$residuals,
  abs_residuals = abs(validation_df$residuals)
)

plot_spatial_residuals <- ggplot() +
  geom_sf(data = pr_map, fill = "white", color = "black", size = 0.3) +
  geom_point(data = spatial_residuals, 
             aes(x = lon, y = lat, color = residuals, size = abs_residuals), 
             alpha = 0.7) +
  scale_color_gradient2(low = "blue", mid = "white", high = "red", 
                        midpoint = 0, name = "Resíduos") +
  scale_size_continuous(range = c(2, 8), name = "|Resíduos|") +
  labs(title = "Distribuição Espacial dos Resíduos",
       subtitle = "Resíduos devem ser aleatórios no espaço (sem padrões)",
       x = "Longitude", y = "Latitude") +
  theme_minimal()

print(plot_spatial_residuals)

# PASSO 18: TESTES ESTATÍSTICOS
shapiro_test <- shapiro.test(validation_df$residuals)

cat("\n=== TESTES ESTATÍSTICOS ===\n")
cat("Teste de Normalidade de Shapiro-Wilk:\n")
cat("Estatística W:", round(shapiro_test$statistic, 4), "\n")
cat("Valor-p:", round(shapiro_test$p.value, 4), "\n")

if(shapiro_test$p.value > 0.05) {
  cat("Conclusão: Resíduos seguem distribuição normal (p > 0.05)\n")
} else {
  cat("Conclusão: Resíduos NÃO seguem distribuição normal (p <= 0.05)\n")
}

# PASSO 19: MAPA DE INCERTEZA
grid_coords$variance <- krige_result_poly$krige.var

plot_uncertainty <- ggplot() +
  geom_sf(data = pr_map, fill = "white", color = "black", size = 0.3) +
  geom_tile(data = grid_coords, aes(x = lon, y = lat, fill = variance)) +
  scale_fill_viridis_c(name = "Variância", option = "plasma") +
  geom_point(data = data_points, aes(x = lon, y = lat), color = "red", size = 1.5, alpha = 0.6) +
  labs(title = "Mapa de Incerteza da Krigagem",
       subtitle = "Áreas com maior variância têm maior incerteza nas predições",
       x = "Longitude", y = "Latitude") +
  theme_minimal()

print(plot_uncertainty)

# PASSO 20: INTERPRETAÇÃO FINAL
cat("\n=== INTERPRETAÇÃO DOS RESULTADOS ===\n")
cat("CRITÉRIOS DE AVALIAÇÃO:\n")
cat("1. RMSE e MAE: Quanto menor, melhor o modelo prediz\n")
cat("2. Viés (Bias): Idealmente entre -2 e 2 (próximo de zero)\n") 
cat("3. Correlação: Quanto mais próxima de 1, melhor\n")
cat("4. Resíduos: Devem ser aleatórios, sem padrões espaciais\n")
cat("5. Normalidade: Resíduos normais indicam modelo adequado\n")

# Avaliação prática do modelo
cat("\nAVALIAÇÃO DO SEU MODELO:\n")
if(abs(bias) < 10 && correlation > 0.3 && RMSE < sd(geo_data$data)) {
  cat("✅ MODELO CONSIDERADO ADEQUADO!\n")
  if(abs(bias) < 5) cat(" - Viés baixo (bom)\n") else cat(" - Viés moderado\n")
  if(correlation > 0.7) cat(" - Boa correlação\n") else cat(" - Correlação moderada\n")
  if(shapiro_test$p.value > 0.05) cat(" - Resíduos normais (bom)\n")
} else {
  cat("⚠️ MODELO PODE PRECISAR DE AJUSTES!\n")
  if(abs(bias) >= 10) cat(" - Viés muito alto\n")
  if(correlation <= 0.3) cat(" - Correlação muito baixa\n")
  if(RMSE >= sd(geo_data$data)) cat(" - Erro muito alto\n")
}

# Comparação com desvio padrão dos dados
cat("Desvio padrão dos dados originais:", round(sd(geo_data$data), 3), "\n")
cat("RMSE em relação ao desvio padrão:", round(RMSE/sd(geo_data$data), 3), "\n")

# PASSO 21: ANÁLISE DO VARIOGRAMA DOS RESÍDUOS
# Se o modelo for bom, os resíduos não devem ter estrutura espacial

# Converter resíduos para objeto geodata
residuals_geodata <- as.geodata(cbind(geo_data$coords, validation_df$residuals))

# Calcular variograma dos resíduos
vario_residuals <- variog(residuals_geodata, max.dist = 1)

plot(vario_residuals, main = "Variograma dos Resíduos",
     sub = "Se o modelo for bom, não deve mostrar estrutura espacial")
abline(h = var(validation_df$residuals), col = "red", lty = 2)
legend("topleft", legend = "Variância total dos resíduos", col = "red", lty = 2)

cat("\n=== ANÁLISE DO VARIOGRAMA DOS RESÍDUOS ===\n")
cat("Se o variograma dos resíduos for plano (sem estrutura), o modelo capturou\n")
cat("toda a dependência espacial dos dados. Se mostrar estrutura, o modelo\n") 
cat("pode estar mal especificado.\n")

  Chegamos ao fim deste minicurso, espero que tenha sido uma jornada enriquecedora do ponto de vista teórico e prático. Possibilitando ver a riqueza da análise estatística espacial aplicada.

Agradecimentos especiais

Agradecimentos especiais a: coordenação da EPBEST 2025 pelo convite e parceria; ao estudante monitor de apoio docente que digitalizou todo esse documento em linguagem Rmarkdown (Alex Gregório- estudante da UTFPR), ao Prof. Dr. Edson Martinez parceiro de pesquisa na área de estatística espacial, a UTFPR minha instituição de origem e ao Ministério da Saúde pela possibilidade de atuação técnica em cursos de formação.

7. Referencias bibliográficas

Livro DIGGLE, P. J.; RIBEIRO, P. J. Model-based geostatistics. New York: Springer, 2007. 228 p. (Springer Series in Statistics).

Documentos Online (Sites, Blogs, Repositórios) BRASIL. geobr: Loads Shapefiles of Official Brazilian Administrative Borders in an Easy Way. 2024. Disponível em: https://cran.r-project.org/web/packages/geobr/vignettes/intro_to_geobr.html. Acesso em: [DD m. aaaa]. (Substitua pela data de acesso).

INTERNAUT. Network Map. 2016. Disponível em: https://gist.github.com/internaut/a9a274c72181eaa7f5c3ab3a5f54b996. Acesso em: [DD m. aaaa].

LOVELACE, R.; NOWOSAD, J.; MÜNCHER, M. Spatial Data Science with applications in R. 2023. Disponível em: https://r.geocompx.org/adv-map. Acesso em: [DD m. aaaa].

MARTINEZ, Edson. Spatial Data Analysis Repository. 2024. Disponível em: https://github.com/edsonzmartinez/spatial/tree/main. Acesso em: [DD m. aaaa].

LIZZI, Elisangela. Spatial-Epi2024: Códigos para Análise Espacial. 2024. Disponível em: https://github.com/elisangelalizzi/Spatial-Epi2024. Acesso em: [DD m. aaaa].

Nota sobre Softwares e Pacotes/ “Ferramentas Computacionais”.

As análises e mapas foram produzidos utilizando o ambiente R (R CORE TEAM, 2024) com os pacotes geobr (BRASIL, 2024), entre outros, seguem referências correspondentes e complementares:

R CORE TEAM. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. 2024. Disponível em: https://www.R-project.org/.

BRASIL. geobr: Loads Shapefiles of Official Brazilian Administrative Borders in an Easy Way. R package version 1.8.0. 2024. Disponível em: https://CRAN.R-project.org/package=geobr.