7  Processos Pontuais Espaço-Temporais

No Capítulo 5, exploramos a configuração de eventos cujas localizações espaciais \(s_i \in W\) eram tratadas como estáticas, representando uma realização pontual em um recorte temporal fixo. Entretanto, tal como discutido na transição da geoestatística, a maioria dos processos pontuais de interesse científico desde a propagação de epidemias e a ocorrência de crimes até a sismicidade possui uma dinâmica intrínseca. Onde um evento ocorre é frequentemente indissociável de quando ele ocorre.

Nesta seção, avançamos para o estudo dos Processos Pontuais Espaço-Temporais (STPP), onde o objeto de análise é o par \((s_i, t_i)\), representando a localização e o instante de um evento. De acordo com Jonatan A. González et al. (2016) e Schoenberg, Brillinger, e Guttorp (2002), essa transição exige que abandonemos a visão bidimensional para adotar um domínio misto \(W \times T\), onde \(W \subset \mathbb{R}^2\) é a janela espacial e \(T \subset \mathbb{R}\) é o intervalo temporal. É fundamental ressaltar, como adverte Diggle (2013), que a geometria deste domínio não é puramente euclidiana em \(\mathbb{R}^3\): o tempo possui uma ordem causal e unidades distintas das espaciais, o que impede que o tratemos simplesmente como uma terceira dimensão geométrica.

7.0.1 Definição e Representação

Um processo pontual espaço-temporal \(N\) é definido como uma medida aleatória sobre um subconjunto do espaço-tempo. Para qualquer região \(A \subseteq W\) e intervalo \(B \subseteq T\), \(N(A \times B)\) representa o número de pontos que caem no cilindro espácio-temporal definido por estas fronteiras (Jonatan A. González et al. 2016; Schoenberg, Brillinger, e Guttorp 2002).

A caracterização fundamental de um STPP é feita através da sua intensidade de primeira ordem \(\lambda(s, t)\). Enquanto na vertente puramente espacial a intensidade media pontos por área, aqui ela define a taxa esperada de eventos por unidade de volume espácio-temporal (Sun 2014):

\[\lambda(s, t) = \lim_{|ds \times dt| \to 0} \frac{\mathbb{E}[N(ds \times dt)]}{|ds||dt|} \tag{7.1}\]

Onde:

  • \(\mathbb{E}[N(ds \times dt)]\) é o número esperado de pontos em uma região infinitesimal ao redor de \((s, t)\).

  • \(|ds|\) representa a área da vizinhança espacial infinitesimal.

  • \(|dt|\) representa a duração do intervalo temporal infinitesimal.

Seguindo a taxonomia proposta por Jonatan A. González et al. (2016), os processos pontuais espaço-temporais podem ser classificados em três grandes famílias, dependendo da natureza do fenômeno:

  • Eventos Instantâneos (Snapshots): Eventos que ocorrem em um ponto fixo no tempo e no espaço e não possuem duração (ex: epicentros de terremotos, disparos criminais). Este é o foco principal da modelagem via intensidades condicionais (Ogata 1998).

  • Processos de Nascimento e Morte: Onde os pontos surgem, permanecem por um tempo e depois são removidos (ex: árvores em uma floresta ou indivíduos em um estudo epidemiológico).

  • Trajetórias: Onde os objetos se movem continuamente, gerando caminhos \(Y_i(t) \in W\) (ex: monitoramento de animais via GPS).

Para representar a vizinhança de um evento no espaço-tempo, não utilizamos esferas 3D, mas sim a vizinhança cilíndrica \(B[(s, t), r, h]\), definida por Jonatan A. González et al. (2016) como:

\[B[(s, t), r, h] = \{ (s', t') : \|s - s'\| \le r, |t - t'| \le h \}\]

Esta métrica respeita a separação física entre a distância euclidiana espacial \(r\) e o lag temporal \(h\).

Uma das hipóteses mais cruciais na modelagem inicial de STPP é a separabilidade de primeira ordem. Um processo é dito separável se a sua intensidade puder ser decomposta de forma multiplicativa (Jonatan A. González et al. 2024; Sun 2014):

\[\lambda(s, t) = \lambda_1(s) \lambda_2(t)\]

Nesta formulação, \(\lambda_1(s)\) captura o padrão espacial estável (ex: a densidade populacional de uma cidade) e \(\lambda_2(t)\) captura a tendência temporal global (ex: sazonalidade ou ciclos semanais). Contudo, Medialdea Villanueva (2025) e Siino, Adelfio, e Mateu (2018) enfatizam que, embora matematicamente conveniente, a separabilidade falha em capturar interações dinâmicas, como a propagação de um incêndio florestal ou o contágio de uma doença, onde o padrão espacial muda agressivamente com o tempo.

7.1 Análise Exploratória de Processos Pontuais Espaço-Temporais (ST-ESDA)

A exploração visual de Processos Pontuais Espaço-Temporais (STPP) é uma etapa crítica que precede a modelagem formal, permitindo diagnosticar a estrutura de dependência e a presença de interações dinâmicas. Diferente da vertente geoestatística, onde se analisa a continuidade de um valor, na ST-ESDA de pontos o foco recai sobre a densidade de ocorrências e a configuração geométrica dos eventos no domínio \(W \times T\).

O Cubo Espaço-Tempo

A representação fundamental de um STPP é o Cubo Espaço-Tempo (3D). Nesta construção, as coordenadas geográficas \((x, y)\) formam a base, e o eixo vertical \(Z\) representa o tempo \(t\). De acordo com Diggle (2005) e Schoenberg, Brillinger, e Guttorp (2002), essa visualização permite identificar imediatamente nuvens de pontos que sinalizam agrupamentos espaço-temporais.

Surtos localizados aparecem como colunas densas, enquanto eventos puramente aleatórios (Poisson) assemelham-se a ruído uniformemente distribuído no volume. No R, o pacote stpp (Gabriel, Rowlingson, e Diggle 2013) permite a manipulação desses objetos e a visualização via stan(), que utiliza a biblioteca rgl para interatividade tridimensional.

Dinâmica Temporal e Curvas de Acumulação

Para diagnosticar a tendência temporal global, a ferramenta primária é a curva cumulativa de eventos \(N(t)\). Conforme discutido por (rasmussen_temporal?) e observado nos dados de epidemias de Diggle (2005), um processo de Poisson homogêneo no tempo resulta em uma reta. Desvios para formatos em “S” (sigmoides) indicam processos de aceleração, típicos de mecanismos auto-excitáveis (Hawkes), onde um evento gera réplicas próximas, como terremotos (Ogata 1998) ou contágio criminal (Mohler et al. 2011).

Projeções de Hovmöller e Snapshots de Intensidade

Adaptando o rigor da geoestatística espaço-temporal (ver Capítulo 6), utilizamos os diagramas de Hovmöller para detectar a propagação geográfica. Ao colapsar uma das dimensões espaciais (ex: longitude) e plotar a outra (latitude) contra o tempo, identificamos se a densidade de eventos está se deslocando (ex: uma frente epidêmica ou o espalhamento de um incêndio) (Jonatan A. González et al. 2016).

Complementarmente, Jonatan A. González e Moraga (2023) sugerem a técnica de Snapshots de Intensidade. Em vez de uma única estimativa de kernel espacial, o tempo é fatiado em janelas \([t_k, t_{k+1}]\), e para cada janela calcula-se a intensidade \(\hat{\lambda}(s)\) (Jonatan A. González et al. 2016):

\[\hat{\lambda}(s, t) = \sum_{i=1}^{n} \kappa_\epsilon(s - s_i) \kappa_\delta(t - t_i) \tag{7.2}\]

Onde \(\kappa_\epsilon\) e \(\kappa_\delta\) são os kernels espacial e temporal, respectivamente. Essa técnica permite observar a dissipação ou o surgimento de novos focos de risco.

7.1.1 Pré-processamento e Harmonização de Dados

Antes da análise, os dados brutos exigem um tratamento para garantir que as premissas de um processo pontual ordenado (orderly process) sejam respeitadas (Diggle 2005):

  • Tratamento de Duplicatas: Processos pontuais não permitem múltiplos eventos no mesmo local e tempo exatos. Duplicatas espaciais (ex: vários crimes no mesmo edifício) podem exigir um pequeno jittering (perturbação nas coordenadas) para viabilizar estimadores de segunda ordem (Jonatan A. González et al. 2016).

  • Harmonização Temporal: Conversão de fusos horários e unidades. Gabriel, Rowlingson, e Diggle (2013) recomenda a conversão de datas para formatos numéricos contínuos para facilitar o cálculo de lags temporais.

  • Binning e Cadência: A escolha do passo temporal (ex: diário, semanal) e da malha espacial é crucial. Um passo muito fino gera excesso de zeros (esparsidade), dificultando a convergência de modelos como o LGCP (Diggle 2013). Por outro lado, um passo muito largo mascara interações de curto alcance, como o efeito near-repeat em crimes (Mohler et al. 2011).

Código 
pacman::p_load(stpp, viridis, splancs, ggplot2, ragg)

# Utilizando o conjunto de dados da febre aftosa (fmd)
data(fmd)
data(northcumbria)
fmd_st <- as.3dpoints(fmd[,1], fmd[,2], fmd[,3])

# Plot do Cubo Espaço-Tempo (Projeção estática)
# O tempo é o eixo vertical
plot(fmd_st) 

title(main = "Cubo Espaço-Tempo: Febre Aftosa")

# Histograma Espaço-Tempo (Binning)
# Útil para visualizar a densidade em células discretas
st_hist <- function(pts) {
  hist(pts[,3], breaks = "FD", main = "Contagens Temporais (Bursts)", 
       xlab = "Tempo (Dias)", col = "steelblue", border = "white")
}
st_hist(fmd_st)

animation(fmd_st, runtime = 15, cex = 0.7)
Código 
pacman::p_load(scatterplot3d)

s3d <- scatterplot3d(fmd_st[,1], fmd_st[,2], fmd_st[,3],
                    pch = 16, color = viridis(100)[as.numeric(cut(fmd_st[,3], 100))],
                    type = "h",                    
                    main = "Cubo Espaço-Tempo: Febre Aftosa",
                    xlab = "Longitude (X)",
                    ylab = "Latitude (Y)",
                    zlab = "Tempo (Dias)",
                    angle = 45,                    
                    grid = TRUE, box = FALSE)

s3d$plane3d(0, 0, 0, lty = "dotted")

7.2 Intensidade espaço-temporal \(\lambda(s,t)\)

A análise de processos pontuais espaço-temporais (STPP) inicia-se, pela caracterização das suas propriedades de primeira ordem. Enquanto no caso puramente espacial ( Capítulo 5) focamos na densidade de pontos por unidade de área, no domínio \(W \times T\) a atenção volta-se para a taxa de ocorrência de eventos em um volume infinitesimal de espaço-tempo (Jonatan A. González et al. 2016).

A função intensidade \(\lambda(s, t)\) representa o número esperado de eventos por unidade de área e unidade de tempo em uma vizinhança infinitesimal do ponto \((s, t) \in W \times T\). Formalmente, ela é definida como o limite da esperança matemática da medida de contagem \(N\) sobre um volume infinitesimal dividido pela magnitude desse volume (Sun 2014; Schoenberg, Brillinger, e Guttorp 2002):

\[\lambda(s, t) = \lim_{|ds \times dt| \to 0} \frac{\mathbb{E}[N(ds \times dt)]}{|ds||dt|} \tag{7.3}\]

Onde:

  • \(N(ds \times dt)\) é o número de eventos em uma região espacial \(ds\) ao redor de \(s\) e intervalo de tempo \(dt\) ao redor de \(t\).

  • \(|ds|\) e \(|dt|\) denotam, respectivamente, a medida de Lebesgue (área e duração) de tais intervalos.

A intensidade captura a propensão de um evento ocorrer. Um processo é dito homogêneo se \(\lambda(s, t) = \lambda\) (constante em todo o domínio), o que implica em uma distribuição uniforme de eventos sob a hipótese de independência (Processo de Poisson). No entanto, na prática, a maioria dos fenômenos exibe inhomogeneidade, onde \(\lambda(s, t)\) varia sistematicamente em função da localização e do tempo (Diggle 2013).

7.2.1 Estimadores Não Paramétricos

Quando não se dispõe de conhecimento prévio sobre a forma funcional da intensidade, recorre-se a estimadores não paramétricos, sendo o método de Kernel o mais difundido.

Kernel Espaço-Temporal

O estimador de kernel para STPP é uma generalização da densidade de kernel bidimensional. A intensidade em um ponto arbitrário \((s, t)\) é estimada como a soma ponderada das distâncias espaciais e temporais de todos os \(n\) eventos observados (Jonatan A. González e Moraga 2023; Jonatan A. González et al. 2016):

\[\hat{\lambda}(s, t) = \sum_{i=1}^{n} \frac{\kappa_{\epsilon}(s - s_i)}{e_s(s_i)} \frac{\kappa_{\delta}(t - t_i)}{e_t(t_i)} \tag{7.4}\]

Onde:

  • \(\kappa_{\epsilon}\) e \(\kappa_{\delta}\) são as funções de densidade kernel (usualmente Gaussianas ou de Quartic) com largura de banda (bandwidth) espacial \(\epsilon\) e temporal \(\delta\).

  • \(e_s(s_i)\) e \(e_t(t_i)\) são fatores de correção de borda (ex: correção de Jones-Diggle) que compensam a perda de massa do kernel nas fronteiras de \(W \times T\) (ver Capítulo 5).

A escolha de kernels pode ser separável (como o produto acima) ou anisotrópica, onde as dimensões são tratadas conjuntamente através de uma métrica de distância ajustada, embora o produto de kernels seja preferido pela interpretabilidade e facilidade computacional (Gabriel, Rowlingson, e Diggle 2013).

Escolha da Largura de Banda (Bandwidth)

O desempenho do estimador de kernel depende de \(\epsilon\) e \(\delta\).

  • Heurísticas: Regras baseadas na distribuição dos dados, como a regra de Scott ou a de Silverman (Jonatan A. González et al. 2016).

  • Validação Cruzada (Cross-validation): Minimiza critérios como o erro quadrático médio (MSE) ou maximiza a verossimilhança de deixar-um-fora (leave-one-out likelihood). Como discutido por Jonatan A. González e Moraga (2023), este último é robusto mas computacionalmente intensivo em grandes conjuntos de dados.

Estimadores Separáveis

Sob a hipótese de que o padrão espacial é invariante ao longo do tempo (apenas sua magnitude global muda), o estimador pode ser simplificado como o produto das intensidades marginais (Jonatan A. González et al. 2016):

\[\hat{\lambda}(s, t) = \frac{\hat{\lambda}_{space}(s) \hat{\lambda}_{time}(t)}{n} \tag{7.5}\]

Aqui, \(\hat{\lambda}_{space}(s)\) é a intensidade espacial acumulada (ignorando o tempo) e \(\hat{\lambda}_{time}(t)\) é a intensidade temporal global.

Testes e Diagnóstico de Separabilidade

A separabilidade é uma premissa que se ela for rejeitada, modelos simples não serão suficientes para explicar os dados (Medialdea Villanueva 2025).

A separabilidade implica que \(\lambda(s, t) = \lambda_1(s)\lambda_2(t)\). Fisicamente, significa que a localização de um surto de doença ou crime não se desloca no mapa ao longo do tempo.

O teste baseado na estatística \(\chi^2\) (qui-quadrado), onde o domínio \(W \times T\) é dividido em células (bins). Compara-se o número de pontos observados em cada bin com o esperado sob o produto das margens (Jonatan A. González e Moraga 2023).

Para validar a significância, simula-se \(m\) vezes um processo de Poisson separável com as intensidades marginais estimadas. O p-valor é calculado comparando a estatística dos dados reais com a distribuição gerada pelas simulações.

7.3 Modelos de Poisson não homogêneo

Para explicar a variação da intensidade por meio de fatores externos, utilizamos modelos paramétricos conhecidos como regressões de intensidade (Diggle 2005).

Assume-se que \(\lambda(s, t)\) segue uma forma log-linear associada a um vetor de covariáveis \(Z(s, t)\) (Jonatan A. González et al. 2024):

\[\log \lambda(s, t) = \beta_0 + \sum_{k=1}^{p} \beta_k Z_k(s, t) + \text{offset}(s) \tag{7.6}\]

Onde:

  • \(\beta_0\): Intercepto representando a intensidade de fundo.

  • \(Z_k(s, t)\): Covariáveis que podem variar no espaço e tempo (ex: temperatura diária, proximidade de infraestrutura).

  • Offsets: Termos conhecidos (como log-densidade populacional) que ajustam a intensidade base sem serem estimados (Diggle 2005).

  • Sazonalidade: Modelada via regressão harmônica, usando termos de seno e cosseno para capturar ciclos periódicos (Medialdea Villanueva 2025):

\(\beta_{sin} \sin(2\pi t / \omega) + \beta_{cos} \cos(2\pi t / \omega)\).

Exemplo: Considere o caso dos incêndios florestais no Nepal (Medialdea Villanueva 2025). A intensidade não é apenas espacialmente heterogênea devido à altitude, mas temporalmente dinâmica devido ao clima.

Um modelo robusto para este fenômeno incluiria:

  • Componente Espacial Estática: Distância a estradas ou elevação.

  • Componente Temporal Dinâmica: Velocidade do vento, precipitação e temperatura média do dia.

A equação de predição resultante no modelo IPP (Inhomogeneous Poisson Process) seria:

\[\lambda(s, t) = \exp \{ \beta_0 + \beta_1 \text{Vento}(t) + \beta_2 \text{Temp}(t) + \beta_3 \text{Sazonalidade}(t) + f(s) \}\]

Onde \(f(s)\) pode ser um termo spline espacial (GAM) para capturar o efeito da geografia. Segundo Medialdea Villanueva (2025), a inclusão dessas covariáveis climáticas dinâmicas é o que permite transitar de um modelo meramente descritivo para um modelo preditivo de risco, capturando a interação que a hipótese de separabilidade ignoraria.

7.4 Dependência e interação espaço-tempo

Enquanto a intensidade de primeira ordem (\(\lambda\)) descreve a densidade média de eventos no domínio \(W \times T\), as propriedades de segunda ordem investigam a correlação entre pares de eventos. Elas permitem responder se a ocorrência de um fenômeno em uma coordenada \((s, t)\) aumenta ou diminui a probabilidade de outros eventos ocorrerem nas proximidades espaciais e temporais (Jonatan A. González et al. 2016; Diggle 2013).

Intensidade de Segunda Ordem e Correlação de Pares

O ponto de partida teórico é a intensidade de segunda ordem \(\lambda^{(2)}(\xi_i, \xi_j)\), que define a densidade conjunta de encontrar dois eventos nos locais \(\xi_i = (s_i, t_i)\) e \(\xi_j = (s_j, t_j)\). De acordo com Jonatan A. González et al. (2016), uma forma mais interpretável de descrever esta relação é a função de correlação de pares espaço-temporal (PCF\(_{st}\)), denotada por \(g(r, u)\):

\[g(r, u) = \frac{\lambda^{(2)}(s_i, s_j, t_i, t_j)}{\lambda(s_i, t_i)\lambda(s_j, t_j)} \tag{7.7}\]

Onde:

  • \(r = \|s_i - s_j\|\) é a distância euclidiana espacial.

  • \(u = |t_i - t_j|\) é o lag temporal.

  • \(\lambda(s, t)\) é a intensidade de primeira ordem em cada ponto.

Se \(g(r, u) = 1\), o processo é independente (Poisson), indicando que os eventos não interagem.

  • Se \(g(r, u) > 1\), existe atração ou agrupamento (clustering), comum em processos contagiosos (Mohler et al. 2011).

  • Se \(g(r, u) < 1\), existe repulsão ou inibição, sugerindo competição por espaço ou recursos (Raeisi, Bonneu, e Gabriel 2021).

Extensão da Função K de Ripley: \(K_{st}(r, u)\)

A função \(K\) espaço-temporal (\(K_{ST}\)) é também utilizada para detectar desvios da aleatoriedade em múltiplas escalas. Ela representa a esperança do número de eventos adicionais dentro de um cilindro de raio espacial \(r\) e janela temporal \(u\) ao redor de um evento típico, normalizada pela intensidade (Gabriel, Rowlingson, e Diggle 2013; Diggle 2005):

\[K_{ST}(r, u) = 2\pi \int_{0}^{u} \int_{0}^{r} g(r', u') r' \, dr' \, du'\]

Para um processo de Poisson homogêneo, o valor teórico é \(K_{ST}(r, u) = \pi r^2 u\). Valores observados superiores a este benchmark indicam agregação espaço-temporal significativa (Jonatan A. González et al. 2016).

Estimadores Práticos e Correções de Borda

Para dados reais, o estimador \(\hat{K}_{ST}(r, u)\) exige o tratamento de eventos próximos às fronteiras do mapa e do período de estudo. Seguindo Gabriel, Rowlingson, e Diggle (2013) temos:

\[\hat{K}_{ST}(r, u) = \frac{1}{|W \times T|} \sum_{i=1}^{n} \sum_{j \neq i} \frac{\mathbb{I}(\|s_i - s_j\| \le r, |t_i - t_j| \le u)}{\hat{\lambda}(s_i, t_i)\hat{\lambda}(s_j, t_j) w_{ij} v_{ij}} \tag{7.8}\]

Onde:

  • \(\mathbb{I}(\cdot)\) é a função indicadora de pertinência ao cilindro.

  • \(w_{ij}\) é a correção de borda espacial (proporção da circunferência do círculo dentro de \(W\)).

  • \(v_{ij}\) é a correção de borda temporal (análoga unidimensional para o intervalo \(T\)).

  • \(\hat{\lambda}\) é a intensidade estimada (ver Seção 7.1).

Separabilidade e a Estatística \(D(r, u)\)

Uma questão central na análise de segunda ordem é se a interação é puramente espacial somada a uma interação temporal, ou se existe uma interação sinérgica entre as duas dimensões (Siino, Adelfio, e Mateu 2018).

Sob a hipótese de separabilidade de segunda ordem, a função \(K\) pode ser decomposta como o produto das funções marginais espacial (\(K_s\)) e temporal (\(K_t\)):

\[K_{ST}(r, u) = K_s(r) K_t(u)\]

Para diagnosticar desvios desta hipótese, Diggle (2013) e Jonatan A. González et al. (2016) propõem a estatística de diagnóstico:

\[D(r, u) = \hat{K}_{ST}(r, u) - \hat{K}_s(r) \hat{K}_t(u)\]

Se \(D(r, u) > 0\), há uma interação espaço-temporal positiva (atração específica no espaço-tempo), o que significa que o agrupamento não é explicado apenas pelas tendências isoladas de cada dimensão (Medialdea Villanueva 2025).

Testes de Hipótese e Envelopes de Monte Carlo

Para validar se o agrupamento observado é estatisticamente significativo, utilizamos o Teste de Independência Espaço-Temporal via simulação (Gabriel, Rowlingson, e Diggle 2013; Jonatan A. González e Moraga 2023).

  • Nulidade: A hipótese nula (\(H_0\)) assume que não há interação espaço-temporal (os tempos são independentes das localizações).

  • Permutação: Uma técnica comum, citada por Sun (2014), consiste em permutar aleatoriamente os selos temporais (\(t_i\)) entre as localizações fixas (\(s_i\)). Isso preserva a estrutura espacial e a tendência temporal, mas quebra qualquer ligação entre ambas.

  • Envelopes: Geram-se \(m\) (ex: 99 ou 999) padrões permutados e calcula-se \(\hat{K}_{ST}\) para cada um. O envelope define os limites de variação esperados sob \(H_0\). Se a curva dos dados reais romper o limite superior do envelope, rejeita-se a independência em favor do agrupamento (clustering).

7.5 Modelos paramétricos espaço-temporais

A modelagem paramétrica de STPP visa descrever o mecanismo gerador dos dados através de funções matemáticas cujos parâmetros possuem interpretação científica. Enquanto os estimadores de kernel (ver Seção 7.1) são puramente descritivos, os modelos paramétricos permitem testar hipóteses sobre causalidade, prever riscos futuros e quantificar a incerteza estatística (Diggle 2013; Jonatan A. González et al. 2016).

7.5.1 Processos de Poisson Espaço-Temporais não homogêneos

O Processo de Poisson não homogêneos (IPP) é a classe fundamental onde a variação da intensidade é tratada como puramente determinística. A premissa central é a independência condicional: dado que conhecemos a função intensidade \(\lambda(s, t)\), a ocorrência de um evento em \((s, t)\) não altera a probabilidade de ocorrência em qualquer outro ponto (Diggle 2005; Schoenberg, Brillinger, e Guttorp 2002).

A intensidade é modelada vinculando-a a covariáveis \(Z(s, t)\) através de uma função de ligação logarítmica para garantir a não-negatividade (Jonatan A. González e Moraga 2023):

\[\lambda(s, t) = \exp\left( \beta_0 + \sum_{k=1}^p \beta_k Z_k(s, t) \right) \tag{7.9}\]

A estimativa dos parâmetros \(\beta\) é realizada via Máxima Verossimilhança (MLE). A função de log-verossimilhança para um padrão observado \(\{(s_i, t_i)\}_{i=1}^n\) é (Sun 2014):

\[\ell(\beta) = \sum_{i=1}^n \log \lambda(s_i, t_i; \beta) - \int_T \int_W \lambda(s, t; \beta) \, ds \, dt \tag{7.10}\]

7.5.2 Processos de Cox Log-Gaussianos (ST-LGCP)

Os Processos de Cox Log-Gaussianos (LGCP) são processos duplamente estocásticos. Eles estendem o IPP ao assumir que a própria intensidade é uma superfície aleatória, capaz de capturar heterogeneidades não explicadas pelas covariáveis (Diggle 2013; Møller, Syversveen, e Waagepetersen 1998).

A intensidade estocástica \(\Lambda(s, t)\) é definida como (Diggle 2013; Medialdea Villanueva 2025):

\[\log \Lambda(s, t) = \mu(s, t) + Y(s, t)\]

Onde:

  • \(\mu(s, t)\) é a componente de tendência (similar ao Poisson).

  • \(Y(s, t)\) é um Campo Aleatório Gaussiano (GRF) com média zero e função de covariância \(C((s, t), (s', t'))\).

Devido à natureza latente de \(Y(s, t)\), a verossimilhança é intratável analiticamente. Existem duas rotas principais:

  • MCMC (Monte Carlo via Cadeias de Markov): Utiliza o algoritmo MALA (Metropolis-Adjusted Langevin Algorithm) para amostrar o campo latente, conforme implementado no pacote lgcp (Taylor et al. 2013).

  • Abordagem SPDE/INLA: Proposta por Simpson et al. (2015), esta técnica converte o campo contínuo em uma malha de elementos finitos (mesh), permitindo uma inferência Bayesiana muito mais rápida através de aproximações de Laplace integradas.

7.6 Processos Auto-Excitáveis (Hawkes Espaço-Temporal)

Os Processos de Hawkes são fundamentais para modelar fenômenos de auto-excitação, onde a ocorrência de um evento desencadeia outros temporária e localmente. Eles são amplamente aplicados em sismologia (modelo ETAS) e criminologia (Ogata 1998; Mohler et al. 2011).

Diferente dos anteriores, a intensidade aqui depende explicitamente do histórico \(\mathcal{H}_t\) de eventos passados (Bernabeu, Zhuang, e Mateu 2024):

\[\lambda(s, t | \mathcal{H}_t) = \mu(s, t) + \sum_{i: t_i < t} g(s - s_i, t - t_i)\]

Onde:

  • \(\mu(s, t)\): Intensidade de fundo (background), representando eventos espontâneos.

  • \(g(\cdot)\): Função de disparo (triggering kernel). Ela define como o risco decai no espaço e no tempo após um evento. Um exemplo comum é o produto de um kernel espacial Gaussiano e um decaimento temporal de Omori (lei de potência) (Ogata 1998).

Declustering Estocástico

Uma técnica de inferência crucial para estes modelos é o algoritmo EM de declustering estocástico. Segundo Bernabeu, Zhuang, e Mateu (2024), este método calcula a probabilidade de cada evento ser um evento de fundo ou ter sido gerado por um evento anterior.

Uma nota prática importante feita por Mohler et al. (2011) é o desafio de distinguir entre uma intensidade de fundo \(\mu\) muito alta e variável (tendência) de um processo de excitação forte \(g\) (interação). O rigor no diagnóstico de segunda ordem (ver Seção 7.4) é essencial para resolver essa ambiguidade.

7.7 Prática no R

  1. Preparação e Estruturação de Dados

Diferente dos objetos espaciais estáticos (ppp), os dados espaço-temporais requerem objetos que vinculem as coordenadas \((x, y)\) a um índice temporal \(t\).

  • Pacote stpp : Utiliza a função as.3dpoints() para criar uma matriz tridimensional ou stppp() para definir o padrão dentro de uma janela espacial (owin) e um intervalo temporal (Gabriel, Rowlingson, e Diggle 2013).
Código 
pacman::p_load(stpp, lgcp, stopp, spatstat)
data(fmd)

# Definir a janela espacial (owin) e o limite temporal (tlim)
janela <- owin(c(280000, 400000), c(480000, 580000))
fmd_obj <- stppp(list(data = fmd, tlim = c(1, 200), window = janela))


fmd_3d <- as.3dpoints(fmd[,1], fmd[,2], fmd[,3])
head(fmd_3d)

pacman::p_load(stpp, stopp, lgcp, spatstat, splancs)

data(fmd)
data(northcumbria)

janela_poly <- owin(poly = northcumbria)

fmd_obj <- stppp(list(data = fmd, tlim = c(1, 220), window = janela_poly))

fmd_3d <- as.3dpoints(fmd[,1], fmd[,2], fmd[,3])

df_stopp <- data.frame(x = fmd[,1], y = fmd[,2], t = fmd[,3])
fmd_stp <- stp(df_stopp)

print(class(fmd_obj)) # Deve ser "stppp"
print(class(fmd_stp)) # Deve ser "stp"
  • Pacote lgcp: Exige a função stppp(). Um tratamento crítico aqui é a função integerise(). Como os algoritmos de MCMC e FFT operam de forma mais eficiente em grades discretas, esta função converte tempos reais em inteiros (ex: dias ou horas), garantindo que a intensidade temporal \(\mu(t)\) seja escalonada corretamente (Taylor et al. 2013).
Código 
fmd_int <- integerise(fmd_obj)
print(fmd_int)
  • Pacote stopp: Introduz a classe stp(), que armazena os dados em um data.frame otimizado. Se os eventos ocorrerem em ruas, utiliza-se stlp() para vincular o padrão a uma rede linear (linnet) (dangelo2025stopp?).
Código 
pacman::p_load(stpp, stopp,lgcp, stopp, spatstat,spatstat.data)
data(chicago)

df_raw <- data.frame(
  x = as.numeric(chicago$data$x),
  y = as.numeric(chicago$data$y),
  t = as.numeric(chicago$data$tp)
)

chicago_stp <- stp(df_raw, L = domain(chicago))

print(chicago_stp)
plot(chicago_stp, main = "Crimes em Chicago (Formatado para stopp)")
Código 
L_fmd <- linnet(vertices = ppp(c(300000, 380000), c(500000, 570000), 
                             window = owin(c(290000, 390000), c(490000, 580000))),
                edges = matrix(c(1, 2), ncol = 2))

fmd_stlp <- stp(fmd, L = L_fmd)

plot(fmd_stlp)
  1. Análise Exploratória (ST-ESDA) e Visualização

O objetivo é identificar visualmente o agrupamento (clustering) antes de iniciar o ajuste de parâmetros.

  • Animação: A função animation(obj) do pacote stpp gera uma sequência de frames.
Código 
data(northcumbria) 

fmd_3d <- as.3dpoints(fmd[,1], fmd[,2], fmd[,3])

# runtime define o tempo total em segundos; cex controla o tamanho dos pontos
animation(fmd_3d, runtime = 10, cex = 0.5, s.region = northcumbria)
Código 
pacman::p_load(ggplot2, gganimate, sf)

border_sf <- st_as_sf(SpatialPolygons(list(Polygons(list(Polygon(northcumbria)), "Cumbria"))))
fmd_df <- as.data.frame(fmd)

p_anim <- ggplot() +
  geom_sf(data = border_sf, fill = "gray95", color = "black") +
  geom_point(data = fmd_df, aes(x = X, y = Y, color = ReportedDay), size = 2, alpha = 0.7) +
  scale_color_viridis_c(option = "plasma") +
  labs(title = "Expansão da Febre Aftosa - Dia: {frame_time}",
       subtitle = "Cúmbria, Reino Unido",
       x = "Easting", y = "Northing") +
  theme_minimal() +
  transition_time(ReportedDay) +
  shadow_mark(alpha = 0.3, size = 0.5) 

animate(p_anim, 
        nframes = 100, 
        device = "ragg_png",
        renderer = magick_renderer(), 
        width = 800, height = 600)

Observe se os pontos surgem em “rajadas” em locais específicos (indicação de auto-excitação ou campo latente) ou se aparecem de forma esparsa (Poisson).

  • Cubo Espaço-Tempo Interativo: A função stan() do stpp permite rotacionar o volume 3D.
Código 
# bgpoly desenha a fronteira da janela espacial no fundo do cubo
pacman::p_load(rgl, rpanel)
stan(fmd_3d, bgpoly = northcumbria)

Colunas verticais de pontos indicam persistência espacial ao longo do tempo (ex: um ponto de venda de drogas ou uma fonte de poluição fixa).

  • Intensidade Marginal: plot(stp_obj) no pacote stopp decompõe a saída em três painéis: o mapa espacial, o histograma temporal e a curva cumulativa \(N(t)\) (dangelo2025stopp?).
Código 
# O painel da direita (tcum = TRUE) é essencial para diagnosticar aceleração da taxa
plot(fmd_stp, tcum = TRUE)
  1. Estimativa de Intensidade de Primeira Ordem

Para entender a tendência de fundo \(\lambda(s, t)\), utilizamos estimadores de Kernel.

  • Função kernel2d() (via splancs) ou density.ppp(): Permitem criar “snapshots” temporais.
Código 
fmd_km <- as.3dpoints(fmd[,1]/1000, fmd[,2]/1000, fmd[,3])
regiao_km <- northcumbria/1000

# intensidade temporal marginal mu(t)
Mt <- density(fmd_km[, 3], n = 1000)
mut <- Mt$y[findInterval(fmd_km[, 3], Mt$x)] * nrow(fmd_km)

# Usamos kernel2d do pacote splancs para a componente espacial
pacman::p_load(splancs)
h_fixo <- 5 
Ms <- kernel2d(as.points(fmd_km[, 1:2]), regiao_km, h = h_fixo, nx = 100, ny = 100)

image(Ms, main = "Intensidade Espacial de Fundo (Kernel)")

O stpp oferece a função mse2d() para selecionar o \(\epsilon\) (espacial) que minimiza o erro quadrático médio. Uma largura de banda muito pequena produz um mapa pontilhado; uma muito grande mascara focos epidêmicos (Jonatan A. González e Moraga 2023).

Código 
fmd_pts <- as.points(fmd[,1], fmd[,2])
poly_pts <- northcumbria

busca_h <- mse2d(fmd_pts, poly_pts, nsmse = 30, range = 5000) 

h_otimo <- busca_h$h[which.min(busca_h$mse)]
print(paste("Largura de banda ótima:", round(h_otimo, 2)))
  1. Diagnóstico de Segunda Ordem (Interação)

Aqui determinamos se o espaço e o tempo interagem de forma não separável.

  • Função STIKhat() (stpp): Estima a função \(K\) espaço-temporal.
Código 
pacman::p_load(stpp, spatstat, splancs)

win <- owin(poly = list(x = regiao_km[,1], y = regiao_km[,2]))
fmd_ppp <- ppp(fmd_km[,1], fmd_km[,2], window = win)

# h_otimo = 3.68 km
mus <- density(fmd_ppp, sigma = h_otimo, at = "points")

# 3. Componente Temporal mut (Marginal)
Mt <- density(fmd_km[, 3], n = 1000)
mut <- Mt$y[findInterval(fmd_km[, 3], Mt$x)] * nrow(fmd_km)

lambda_estimada <- mus * (mut / area(win))


print(paste("Vetor lambda criado com", length(lambda_estimada), "observações."))
Código 
fmd_st_km <- as.3dpoints(fmd_km[,1], fmd_km[,2], fmd_km[,3])

# Lags para KM (u) e Dias (v)
u_km <- seq(0, 10, by = 1) 
v_dias <- seq(0, 15, by = 1)

stik <- STIKhat(xyt = fmd_st_km, 
                lambda = lambda_estimada, 
                s.region = regiao_km, 
                t.region = c(28, 200),
                dist = u_km, 
                times = v_dias, 
                infectious = TRUE)

plotK(stik, type = "persp", main = "Superfície da Função K-ST (Febre Aftosa)")
  • Função PCFhat() (stpp): Estima a correlação de pares \(g(r, u)\).

Se o gráfico 3D da PCF mostrar um pico em distâncias \(r\) e \(u\) curtas, o processo é auto-excitável. Se houver um “buraco” (valores próximos a zero) em \(r\) pequeno, o processo exibe inibição (ex: o incêndio consumiu o combustível e impede outro fogo imediato no mesmo local) (Raeisi, Bonneu, e Gabriel 2021).

  1. Ajuste de Modelos Paramétricos
  • Processos de Poisson (IPP)

Utiliza-se a função stppm() do pacote stopp. Exemplo: stppm(dados, formula = ~ x + y + t + covariavel).

Código 
modelo_ipp <- stppm(fmd_stp, formula = ~ x + y + t)

summary(modelo_ipp)
plot(modelo_ipp)

O pacote gera pontos “dummy” para aproximar o integral da verossimilhança.

Verifique o summary(). Coeficientes positivos para uma covariável (ex: temperatura) indicam que o aumento desta eleva a taxa de eventos (dangelo2025stopp?).

Código 
summary(modelo_ipp)

plot(modelo_ipp)

5.2 Processos de Cox Log-Gaussianos (LGCP)

A inferência em LGCP é computacionalmente intensiva. Para viabilizá-la, o pacote lgcp utiliza a Transformada Rápida de Fourier (FFT). Isso exige a definição de uma grade retangular onde o número de células em cada dimensão é, preferencialmente, uma potência de 2 (ex: \(128 \times 128\)). De acordo com Taylor et al. (2013), os parâmetros de covariância (\(\sigma^2\) para variância, \(\phi\) para escala espacial e \(\theta\) para decaimento temporal) devem ser estimados previamente, geralmente por métodos de momentos ou contraste mínimo.

Código 
# theta: escala temporal
parametros <- lgcppars(sigma = 1.6, phi = 1.9, theta = 1.4)

largura_celula <- 2 

# mala.length: iterações totais; burnin: descarte inicial; retain: afinamento (thinning)
mcmc_ctrl <- mcmcpars(mala.length = 100000, 
                      burnin = 20000, 
                      retain = 100, 
                      adaptivescheme = andrieuthomsh(inith = 1))

Execução da Predição Espácio-Temporal

A função lgcpPredict() realiza a inferência Bayesiana para o campo latente \(Y\). Ela utiliza o algoritmo MALA (Metropolis-Adjusted Langevin Algorithm), que usa o gradiente da densidade para guiar as simulações. Segundo o rigor de Taylor et al. (2013), é possível calcular “on-line” as probabilidades de excedência, que indicam o risco de a intensidade ultrapassar um limiar crítico \(k\) (ex: o dobro da média).

Código 
fmd_lgcp <- fmd_obj 

grid_obj <- lgcpGrid(fmd_lgcp, cellwidth = 2000) 

pars <- lgcppars(sigma = 1.5, phi = 3000, theta = 10) # phi em metros

mcmc_ctrl <- mcmcpars(mala.length = 5000, burnin = 1000, retain = 10,
                      adaptivescheme = andrieuthomsh(inith = 1))

res_lgcp <- lgcpPredict(xyt = fmd_lgcp, 
                        T = max(fmd[,3]), 
                        laglength = 5,
                        model.parameters = pars,
                        grid = grid_obj,
                        mcmc.control = mcmc_ctrl)

Diagnóstico de Convergência (MALA)

Um aspecto crítico do MALA é a escolha do tamanho do passo \(h\). O pacote lgcp implementa um esquema adaptativo de (andrieu2008tutorial?). O pesquisador deve monitorar a convergência através da função hvals(). Se o gráfico mostrar que \(h\) estabilizou em um valor constante após o burn-in, a cadeia atingiu a ergodicidade necessária (Taylor et al. 2013).

Código 
plot(hvals(res_lgcp), type = "l", 
     main = "Convergência do Parâmetro h (MALA)", 
     xlab = "Iteração", ylab = "h")

trace_y <- extract(res_lgcp, x = 10, y = 10, t = 1, s = -1)
plot(trace_y, type = "l", main = "Trace Plot do Campo Latente Y")

O resultado final de um LGCP é uma superfície de risco probabilístico. A saída mais valiosa para a tomada de decisão é o mapa de excedência. Conforme destacado por Diggle (2013), uma afirmação como “há 90% de probabilidade de o risco nesta coordenada ser 3 vezes maior que a média histórica” permite que gestores identifiquem surtos (outbreaks) reais, filtrando a variação aleatória comum de processos de Poisson.

Código 
plotExceed(res_lgcp, fun = "exceed", index = 3)

media_intensidade <- intens(res_lgcp)
plot(media_intensidade, main = "Intensidade Média Predita (Poisson Rate)")

  • Campo Latente (\(Y\)): Representa “fatores ambientais escondidos”. Se for alto e persistente, há um fator local (ex: água parada) favorecendo o evento.

  • Probabilidade de Excedência: Valores próximos a 1 (cor vermelha no mapa) indicam áreas de intervenção urgente.

  • H-values: Se o gráfico for errático ou não estabilizar, os resultados de predição são matematicamente inválidos e o mala.length deve ser aumentado (Taylor et al. 2013).

  1. Validação e Verificação

A etapa final e mais rigorosa da modelagem de Processos Pontuais Espaço-Temporais (STPP) é a validação. Um modelo bem ajustado deve ser capaz de descrever não apenas a densidade média de eventos (primeira ordem), mas também a estrutura de interação residual (segunda ordem). Conforme discutido por (dangelo2025stopp?) e Taylor et al. (2013), utilizamos diagnósticos baseados em resíduos e simulações para garantir que o modelo não “mentiu” sobre o mecanismo gerador dos dados.

  • Envelopes de Monte Carlo

Para avaliar se a estrutura de dependência do modelo capturou o agrupamento real, utilizamos os Envelopes de Monte Carlo. O procedimento consiste em simular o modelo ajustado \(m\) vezes (geralmente \(n=39, 99\) ou \(199\)) para gerar uma faixa de variação esperada (área cinza). Segundo Gabriel, Rowlingson, e Diggle (2013), se a curva da estatística observada (como a função \(K_{ST}\)) romper os limites deste envelope, o modelo é estatisticamente inadequado para representar a interação entre os pontos.

Código 
env_modelo <- envelope(fmd_obj, STIKhat, nsim = 39, 
                       simulate = expression(rstpp(lambda = modelo_ipp$l)),
                       verbose = FALSE)

plot(env_modelo, main = "Envelopes de Monte Carlo para o Modelo IPP")

Diagnóstico Global e Local (STOPP)

O pacote stopp fornece ferramentas avançadas para identificar falhas estruturais e pontos influentes. A função globaldiag() realiza uma verificação de bondade de ajuste comparando a função \(K\) espacial e temporal simultaneamente. Complementarmente, a função localdiag() é utilizada para identificar “outliers” espácio-temporais pontos individuais que possuem um valor de resíduo significativamente alto, indicando que a intensidade do modelo falhou localmente (dangelo2025stopp?).

Código 
diag_global <- globaldiag(modelo_ipp)
plot(diag_global)

# p = 0.9 indica que estamos buscando os 10% de pontos menos explicados pelo modelo
res_local <- localdiag(fmd_stp, modelo_ipp$l, p = 0.9)

plot(res_local)

Interpretação de Pontos de Influência

Pontos identificados como discrepantes sugerem que o modelo de Poisson (IPP) pode ser insuficiente, possivelmente devido à falta de uma covariável importante ou à existência de um processo de auto-excitação (Hawkes) não contabilizado. A função infl() permite visualizar a contribuição individual de cada ponto “outlier” para a estatística de segunda ordem. Conforme destacado por (dangelo2025stopp?), esta análise local é o que permite transitar de uma análise global simplista para uma compreensão detalhada de surtos inesperados no espaço-tempo.

Código 
infl(res_local)

  • Envelopes: O rompimento do envelope indica que o modelo de Poisson é “pobre” e talvez um modelo de Cox (LGCP) ou Hawkes seja necessário para explicar a atração entre os pontos (Taylor et al. 2013).

  • Globaldiag: Uma diferença sistemática crescente indica que o modelo falha em grandes escalas espaciais ou temporais.

  • Localdiag: Se os pontos discrepantes se agrupam em uma região específica, isso é evidência de um “ponto quente” (hotspot) persistente que exige investigação de campo ou novas variáveis ambientais (Mohler et al. 2011).

Importante

  1. Alguns pacotes tem apresentado alguns erros, sugiro ver sempre se houve atualização do pacote.

  2. Os livros Cressie e Wikle (2011) e Wikle, Zammit-Mangion, e Cressie (2019) tem a teoria e a prática respectivamente.