8 Prática - Regressão Espacial

Docentes

Nota

Essa prática foi adaptada do roteiro criado por Luis Felipe Bortolatto da Cunha para edições passadas dessa disciplina.

Download do arquivo .qmd

Sugestão: abra o projeto no RStudio (que está usando ao longo do curso), e faça o download deste arquivo (em formato .qmd) com o código abaixo.

download.file(
  url = "https://raw.githubusercontent.com/beatrizmilz/ESHT011-21-analise-dados-planejamento-territorial/refs/heads/main/praticas/08_regressa-espacial.qmd",
  destfile = "08_regressa-espacial.qmd",
)

8.1 Introdução

Este roteiro tem como objetivo auxiliar na execução de modelos de regressão espacial (Spatial Lag, Spatial Error e GWR) com o software R.

Para executar os modelos de regressão espacial vamos utilizar dados demográficos e de consumo de água de 2010, para uma amostra de 4.417 municípios, extraídos do Censo Demográfico (IBGE) e do Sistema Nacional de Informações sobre Saneamento (SNIS), conforme análise apresentada por Carmo et al., 2013.

A base de dados está disponível para download no endereço abaixo:

https://1drv.ms/u/s!AjettDH-3Gbni9kP8cPtA04EBQQ7xQ?e=4dyWz5

8.2 Carregando os pacotes

Os pacotes necessários para essa prática são:

tidyverse: para manipulação e visualização de dados;
sf: para trabalhar com dados geoespaciais vetoriais (shapefile ou geopackage)
spdep: para calcular a autorrelação espacial e a matriz de pesos espaciais
spatialreg: para calcular os modelos de regressão espacial global
spgrw: para calcular o modelo de regressão espacial local

Podemos instalar os pacotes:

install.packages("tidyverse")
install.packages("sf")
install.packages("spdep")
install.packages("spatialreg")
install.packages("spgwr")

Aviso

Atenção: dependendo da versão de R que você está usando, você pode precisar instalar o pacote Rcpp antes dos pacotes de análise espacial. Se você encontrar um erro na instalação dos pacotes acima, instale o pacote citado e tente novamente. Se o problema persistir, por favor entre em contato com a equipe responsável pela disciplina.

install.packages("Rcpp")

E então carregá-los:

library(tidyverse)
library(sf)
library(spdep)
library(spatialreg)
library(spgwr)

8.3 Importando os dados

Para facilitar o download, você pode usar o código abaixo, que cria uma pasta chamada “dados” e baixa o arquivo geopackage para essa pasta. O arquivo será salvo com o nome agua_rede_sf.gpkg:

# Criando a pasta para armazenar os dados
fs::dir_create("dados")

# Download do arquivo geopackage
download.file("https://github.com/beatrizmilz/ESHT011-21-analise-dados-planejamento-territorial/raw/refs/heads/main/praticas/dados/agua_rede_sf.gpkg", 
              destfile = "dados/agua_rede_sf.gpkg", 
              mode = "wb")

A base de dados está em formato geopackage (.gpkg) e pode ser importada com a função read_sf().

agua_rede_sf <- read_sf("dados/agua_rede_sf.gpkg")

8.4 Explorando os dados

Os arquivos geoespaciais que trabalharemos são chamados de simple features (sf) no R. A sua estrutura é semelhante à de uma tabela, mas trazendo informações do Sistema de Referência de Coordenadas (CRS) e uma variável adicional: a geometria (geom ou geometry), com as coordenadas que permitem mapear os pontos, linhas ou polígonos.

Executando o nome do objeto é possível obter algumas dessas informações.

agua_rede_sf

Simple feature collection with 5567 features and 25 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -73.99045 ymin: -33.75208 xmax: -28.83609 ymax: 5.271841
Geodetic CRS:  SIRGAS 2000
# A tibble: 5,567 × 26
   ID_IBGE DOMICIL  REDE PROPREDE ID_SNIS NOME_MUN UF    REGIAO    PIB RENDAPITA
     <dbl>   <int> <int>    <dbl>   <int> <chr>    <chr> <chr>   <dbl>     <dbl>
 1 1100015    7270  1774    0.244  110001 Alta Fl… RO    Norte  1.67e5      468.
 2 1100023   27166  9576    0.352  110002 Ariquem… RO    Norte  6.13e5      673.
 3 1100031    1975   474    0.24   110003 Cabixi   RO    Norte  4.17e4      447.
 4 1100049   24215 18428    0.761  110004 Cacoal   RO    Norte  5.83e5      719.
 5 1100056    5348  1816    0.340  110005 Cerejei… RO    Norte  8.69e4      553.
 6 1100064    5959  3617    0.607  110006 Colorad… RO    Norte  1.09e5      508.
 7 1100072    2649   344    0.130  110007 Corumbi… RO    Norte  5.84e4      402.
 8 1100080    3696   609    0.165  110008 Costa M… RO    Norte  4.97e4      353.
 9 1100098    8683  2385    0.275  110009 Espig<e… RO    Norte  1.53e5      572.
10 1100106   10684  4780    0.447  110010 Guajar<… RO    Norte  2.35e5      476.
# ℹ 5,557 more rows
# ℹ 16 more variables: GINI <dbl>, IDH <dbl>, IDH_CLASS <chr>, GE012 <int>,
#   AG001 <int>, AG020 <dbl>, AG022 <int>, CONSUMO1 <dbl>, CONSUMO2 <dbl>,
#   .fitted <dbl>, .resid <dbl>, .std.resid <dbl>, .hat <dbl>, .sigma <dbl>,
#   .cooksd <dbl>, geom <MULTIPOLYGON [°]>

A função names() exibe os nomes das variáveis.

names(agua_rede_sf)

 [1] "ID_IBGE"    "DOMICIL"    "REDE"       "PROPREDE"   "ID_SNIS"   
 [6] "NOME_MUN"   "UF"         "REGIAO"     "PIB"        "RENDAPITA" 
[11] "GINI"       "IDH"        "IDH_CLASS"  "GE012"      "AG001"     
[16] "AG020"      "AG022"      "CONSUMO1"   "CONSUMO2"   ".fitted"   
[21] ".resid"     ".std.resid" ".hat"       ".sigma"     ".cooksd"   
[26] "geom"

Os nomes das variáveis estão codificados, apresentando correspondência às variáveis descritas abaixo:

Código	Descrição
ID_IBGE	Código IBGE (7 dígitos)
DOMICIL	Quantidade de Domicílios
REDE	Quantidade de Domicílios com Acesso à Rede Geral de Água
PROPREDE	Proporção de Domicílios com com Acesso à Rede Geral de Água (REDE/DOMICIL)
ID_SNIS	Código IBGE (6 dígitos)
NOME_MUN	Nome do Município
UF	Unidade da Federação
REGIAO	Região do País
PIB	Produto Interno Bruto 2010
RENDAPITA	Renda per Capita 2010
GINI	Índice GINI 2010
IDH	Índice de Desenvolvimento Humano 2010
IDH_CLASS	Classificação do Índice de Desenvolvimento Humano 2010: Muito Alto >= 0,9; Alto >= 0,8; Médio >= 0,5; Baixo < 0,5.
GE012	População Total Residente no Município
AG001	População Total Atendida com Abastecimento de Água
AG020	Volume Micromedido nas Economias Residenciais Ativas de Agua - 1.000 m3/ano
AG022	Quantidade de Economias Residenciais Ativas Micromedidas
CONSUMO1	Consumo de Água per capita - População Total - m3/ano (AG020/GE012)
CONSUMO2	Consumo de Água per capita - População Atendida - m3/ano (AG020/AG001)

Além delas, estão presentes outras variáveis que trazem os resultados do modelo de regressão linear (.fitted, .resid, .std.resid, .hat, .sigma, .cooksd), conforme apresentado no Roteiro 5.

ATENÇÃO: Algumas considerações sobre os dados antes de rodar o modelo:

O modelo não vai funcionar em uma base de dados com valores faltantes (NA).
É necessário trazer em uma coluna as coordenadas do centróide dos polígonos (longitude/latitude).
Os modelos de regressão espacial são complexos e exigem um tempo de processamento muito maior que os modelos lineares! Nos nossos testes, o tempo de processamento variou entre 10 a 30 minutos. Portanto, recomendamos rodar o modelo apenas em uma amostra dos dados (ex: UF == "SP").

As funções apresentadas nas aulas anteriores funcionam em dados geoespaciais de forma semelhante à tabelas. Vamos usá-las para resolver a primeira e segunda consideração.

Para filtrar apenas as observações completas, ou seja, sem valores faltantes (NA), podemos aplicar a função drop_na(), levando em consideração as variáveis que serão utilizadas (ID, Variável Y e Variável X), para remover as linhas que contém NA nas variáveis indicadas.

agua_rede_sf_sem_na <- agua_rede_sf |>
  drop_na(ID_IBGE, CONSUMO1, RENDAPITA)

Dos 5567 municípios presentes na base de dados, apenas 4417 tinham informações completas do ID_IBGE, CONSUMO1 e RENDAPITA. O modelo será aplicado apenas nessas observações.

Amostragem: o código abaixo cria uma amostra apenas com os dados do estado de São Paulo.

agua_rede_sf_sp <- agua_rede_sf_sem_na |>
  filter(UF == "SP")

Para executar o modelo de regressão espacial local (GWR), é necessário ter variáveis que descrevem a longitude e latitude do centróide dos polígonos.

A função st_centroid() calcula o centróide de um polígono, porém ela só funciona se as geometrias forem consideradas “válidas”. Podemos fazer essa checagem com a função st_is_valid():

# Checando se a geometria é valida
geometrias_validas <- st_is_valid(agua_rede_sf_sp)

# Somando quantas geometrias são inválidas
sum(geometrias_validas == FALSE)

[1] 0

Agora podemos calcular o centróide:

agua_rede_sf_sp_centroide <- agua_rede_sf_sp |>
  mutate(LON = st_coordinates(st_centroid(agua_rede_sf_sp))[,1],
         LAT = st_coordinates(st_centroid(agua_rede_sf_sp))[,2])

Warning: There were 2 warnings in `stopifnot()`.
The first warning was:
ℹ In argument: `LON = st_coordinates(st_centroid(agua_rede_sf_sp))[, 1]`.
Caused by warning:
! st_centroid assumes attributes are constant over geometries
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.

Agora temos as colunas LAT e LON, mantendo a coluna geom:

agua_rede_sf_sp_centroide |> 
  select(LAT, LON, geom)

Simple feature collection with 510 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -53.10986 ymin: -25.31232 xmax: -44.16137 ymax: -19.77966
Geodetic CRS:  SIRGAS 2000
# A tibble: 510 × 3
     LAT   LON                                                              geom
   <dbl> <dbl>                                                <MULTIPOLYGON [°]>
 1 -21.6 -51.1 (((-51.09093 -21.47214, -51.08788 -21.46944, -51.08722 -21.46444…
 2 -21.3 -49.7 (((-49.69668 -21.3443, -49.69703 -21.34019, -49.69894 -21.33772,…
 3 -21.9 -46.7 (((-46.73069 -21.9447, -46.73192 -21.93454, -46.73487 -21.9303, …
 4 -22.5 -46.6 (((-46.635 -22.45781, -46.62814 -22.45148, -46.62519 -22.44678, …
 5 -22.9 -49.3 (((-49.28903 -22.77028, -49.28139 -22.77217, -49.27747 -22.77183…
 6 -22.6 -47.9 (((-47.86292 -22.60939, -47.86517 -22.61207, -47.86423 -22.61583…
 7 -22.6 -49.1 (((-49.0167 -22.38629, -49.00953 -22.38424, -48.99643 -22.38395,…
 8 -23.5 -47.9 (((-47.86567 -23.51306, -47.86243 -23.50943, -47.85984 -23.51047…
 9 -21.9 -51.4 (((-51.37102 -21.87702, -51.36865 -21.87703, -51.36898 -21.87891…
10 -20.5 -49.1 (((-49.14974 -20.57783, -49.1535 -20.57674, -49.15553 -20.57754,…
# ℹ 500 more rows

8.5 Modelos de regressão espacial global (Spatial Lag e Spatial Error)

Após a análise espacial dos resíduos - que permitiu observar que as observações não são independentes espacialmente, vamos incorportar a estrutura de dependência espacial no modelo.

Os modelos globais incluem no modelo de regressão um parâmetro para capturar a estrutura de autocorrelação espacial na área de estudo como um todo.

O Spatial Lag (SAR) atribue a autocorrelação espacial à variável resposta Y (lag), enquanto o Spatial Error (CAR) atribue a autocorrelação espacial ao erro.

8.5.1 Preparando os dados

O ponto de partida para o Spatial Lag e Spatial Error é o modelo de regressão linear (Roteiro 5), criado com a função lm().

modelo1 <- lm(formula = CONSUMO1 ~ RENDAPITA, 
              data = agua_rede_sf_sp_centroide, 
              na.action = na.exclude)

Também é necessário criar um arquivo com a lista de vizinhança, com a função poly2nb(). Foi usado o critério de vizinhança do tipo queen. Para usar outros critérios, lembre-se de consultar a documentação da função.

vizinhanca <- poly2nb(pl = agua_rede_sf_sp_centroide, 
             row.names = agua_rede_sf_sp_centroide$ID_IBGE)

Warning in poly2nb(pl = agua_rede_sf_sp_centroide, row.names = agua_rede_sf_sp_centroide$ID_IBGE): some observations have no neighbours;
if this seems unexpected, try increasing the snap argument.

Warning in poly2nb(pl = agua_rede_sf_sp_centroide, row.names = agua_rede_sf_sp_centroide$ID_IBGE): neighbour object has 5 sub-graphs;
if this sub-graph count seems unexpected, try increasing the snap argument.

Visualizando o sumário desse arquivo, é possível ter uma ideia da distribuição da vizinhança.

summary(vizinhanca)

Neighbour list object:
Number of regions: 510 
Number of nonzero links: 2290 
Percentage nonzero weights: 0.8804306 
Average number of links: 4.490196 
3 regions with no links:
3520400, 3543204, 3546702
5 disjoint connected subgraphs
Link number distribution:

  0   1   2   3   4   5   6   7   8   9  10  11  20 
  3  15  50 103 114  82  75  34  15   9   8   1   1 
15 least connected regions:
3500402 3500600 3501400 3503158 3504701 3504909 3508702 3523701 3527108 3534302 3534708 3539905 3546108 3546207 3551900 with 1 link
1 most connected region:
3550308 with 20 links

Agora é possível criar uma matriz de pesos espaciais, com as funções nb2mat() e mat2listw(). Em vários momentos é necessário definir o argumento adicional zero.policy = TRUE, que permite obter os resultados mesmo que uma ou mais observações não tenham nenhum vizinho. Não se esqueça de adicionar este parâmetro sempre que indicado neste roteiro.

wm <- nb2mat(neighbours = vizinhanca, style='B', zero.policy = TRUE)

rwm <- mat2listw(x = wm, style='W', zero.policy = TRUE)

Warning in mat2listw(x = wm, style = "W", zero.policy = TRUE): neighbour object
has 5 sub-graphs

Agora já é possível realizar o teste de autocorrelação espacial dos resíduos, com a função lm.morantest(), que levará como argumentos: o modelo (modelo1), a matriz de pesos espaciais (rwm) e a tipo de hipótese alternativa (“two-sided”, que representa o teste bicaudal).

lm.morantest(model = modelo1, 
             listw = rwm, 
             alternative = "two.sided",
             zero.policy = TRUE)


    Global Moran I for regression residuals

data:  
model: lm(formula = CONSUMO1 ~ RENDAPITA, data =
agua_rede_sf_sp_centroide, na.action = na.exclude)
weights: rwm

Moran I statistic standard deviate = 9.5448, p-value < 2.2e-16
alternative hypothesis: two.sided
sample estimates:
Observed Moran I      Expectation         Variance 
     0.297470015     -0.002570733      0.000988155

O resultado do teste é significativo (p-valor < 0,00000000000000022), portanto há dependência espacial nos resíduos e um modelo de regressão espacial global pode ser aplicado.

Mas antes disso serão aplicados quatro testes de hipótese para identificar qual é o modelo mais adequado (Spatial Lag ou Spatial Error):

LMerr: teste LM simples para dependência do erro
LMlag: teste LM simples para uma variável dependente espacialmente defasada
RLMerr: teste LM robusto para dependência do erro
RLMlag: teste LM robusto para uma variável dependente espacialmente defasada

lm.LMtests(model = modelo1,
          listw = rwm,
          test = c("LMerr","LMlag","RLMerr","RLMlag"),
          zero.policy = TRUE)

Please update scripts to use lm.RStests in place of lm.LMtests


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = CONSUMO1 ~ RENDAPITA, data =
agua_rede_sf_sp_centroide, na.action = na.exclude)
test weights: listw

RSerr = 89.586, df = 1, p-value < 2.2e-16


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = CONSUMO1 ~ RENDAPITA, data =
agua_rede_sf_sp_centroide, na.action = na.exclude)
test weights: listw

RSlag = 72.74, df = 1, p-value < 2.2e-16


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = CONSUMO1 ~ RENDAPITA, data =
agua_rede_sf_sp_centroide, na.action = na.exclude)
test weights: listw

adjRSerr = 17.354, df = 1, p-value = 3.102e-05


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = CONSUMO1 ~ RENDAPITA, data =
agua_rede_sf_sp_centroide, na.action = na.exclude)
test weights: listw

adjRSlag = 0.50775, df = 1, p-value = 0.4761

Como interpretar os resultados dos testes de hipótese dos multiplicadores de Lagrange?

Primeiro examinamos os testes simples (LMerr e LMlag). Se eles não forem significativos (p-valor > 0,05), a hipótese de modelagem do erro ou lag não se sustenta. Se apenas um deles for significativo (p-valor < 0,05): problema resolvido - mas é bastante comum que ambos sejam significativos. Se ambos forem significativos, é necessário conferir os testes robustos.

Olhando para os multiplicadores de Lagrange robustos (RLMlag e RLMerr), normalmente, apenas um deles será significativo, ou um terá uma ordem de magnitude mais significativo do que o outro (por exemplo, p < 0,000003 em comparação com p < 0,03). Nesse caso, a decisão é simples: estimar o modelo de regressão espacial correspondendo à estatística mais “robusta” (ou significativa). No caso raro de ambos serem altamente significativos, escolha o modelo com o maior valor para a estatística de teste. No entanto, nesta situação, é necessário algum cuidado, uma vez que podem existir outras fontes de erros de especificação. Outra ação a se tomar é alterar a especificação básica (ou seja, a parte não espacial) do modelo.

Também há casos raros em que nenhuma das estatísticas do teste robusto é significativa. Nesses casos, problemas de especificação mais sérios provavelmente estão presentes e devem ser resolvidos primeiro. Por outros erros de especificação, nos referimos a problemas com algumas das outras suposições da regressão linear (Roteiro 5).

Neste caso, os dois testes simples são significativos. Mas o teste de erro robusto é mais significativo que o teste de lag robusto, então o modelo mais apropriado seria o Spatial Error (CAR). De qualquer forma, vamos ensinar a aplicar os dois modelos (SAR e CAR) na sequência.

8.5.2 Spatial Autoregressive Model (Spatial Lag/SAR)

Para aplicar o modelo Spatial Lag usaremos a função lagsarlm(), especificando a fórmula, os dados e a matriz de pesos espaciais.

sar_modelo <- lagsarlm(formula = CONSUMO1 ~ RENDAPITA,
                       data = agua_rede_sf_sp_centroide,
                       listw = rwm,
                       zero.policy = TRUE)

Para ver os resultados do novo modelo, você pode acessar o sumário.

summary(sar_modelo)


Call:
lagsarlm(formula = CONSUMO1 ~ RENDAPITA, data = agua_rede_sf_sp_centroide, 
    listw = rwm, zero.policy = TRUE)

Residuals:
     Min       1Q   Median       3Q      Max 
-51.5374  -6.2923  -1.1463   3.5627 128.1935 

Type: lag 
Regions with no neighbours included:
 3520400 3543204 3546702 
Coefficients: (asymptotic standard errors) 
             Estimate Std. Error z value  Pr(>|z|)
(Intercept) 7.5769500  2.7147528   2.791  0.005254
RENDAPITA   0.0316423  0.0030092  10.515 < 2.2e-16

Rho: 0.35367, LR test value: 60.126, p-value: 8.8818e-15
Asymptotic standard error: 0.045679
    z-value: 7.7426, p-value: 9.77e-15
Wald statistic: 59.948, p-value: 9.77e-15

Log likelihood: -2061.04 for lag model
ML residual variance (sigma squared): 183.69, (sigma: 13.553)
Number of observations: 510 
Number of parameters estimated: 4 
AIC: 4130.1, (AIC for lm: 4188.2)
LM test for residual autocorrelation
test value: 1.9034, p-value: 0.1677

O modelo SAR possui um melhor desempenho que o modelo de regressão linear, indicado pelo AIC (4130), que é menor que o AIC do modelo de regressão linear (4188). Além disso, o teste de autocorrelação espacial dos resíduos não é significativo (p-valor = 0,16), indicando que após a modelagem os resíduos não são mais correlacionados espacialmente.

8.5.3 Conditional Autoregressive Model (Spatial Error/CAR)

Para aplicar o modelo Spatial Error usaremos a função errorsarlm(), especificando a fórmula, os dados e a matriz de pesos espaciais.

car_modelo <- errorsarlm(formula = CONSUMO1 ~ RENDAPITA,
                        data = agua_rede_sf_sp_centroide,
                        listw = rwm,
                        zero.policy = TRUE)

É possível visualizar os resultados do modelo usando o sumário.

summary(car_modelo)


Call:
errorsarlm(formula = CONSUMO1 ~ RENDAPITA, data = agua_rede_sf_sp_centroide, 
    listw = rwm, zero.policy = TRUE)

Residuals:
     Min       1Q   Median       3Q      Max 
-49.3458  -6.4307  -1.1948   3.6627 126.9162 

Type: error 
Regions with no neighbours included:
 3520400 3543204 3546702 
Coefficients: (asymptotic standard errors) 
              Estimate Std. Error z value  Pr(>|z|)
(Intercept) 22.5150675  2.4084788  9.3483 < 2.2e-16
RENDAPITA    0.0341102  0.0031534 10.8170 < 2.2e-16

Lambda: 0.42026, LR test value: 69.624, p-value: < 2.22e-16
Asymptotic standard error: 0.049799
    z-value: 8.4391, p-value: < 2.22e-16
Wald statistic: 71.218, p-value: < 2.22e-16

Log likelihood: -2056.291 for error model
ML residual variance (sigma squared): 177.77, (sigma: 13.333)
Number of observations: 510 
Number of parameters estimated: 4 
AIC: 4120.6, (AIC for lm: 4188.2)

O AIC do modelo CAR (4120) é menor que o AIC do modelo de regressão linear (4188), indicando um melhor resultado. Além disso, todos os valores estimados são significativos (p-valor < 0,05).

8.5.4 Exportando os resíduos do modelo

Para verificar a autocorrelação espacial no GeoDa, é necessário exportar os resíduos do modelo. O código abaixo indica como executar essa etapa para os dois modelos, criando as variáveis (colunas) sar_residuos e car_residuos.

agua_rede_sar_car <- agua_rede_sf_sp_centroide |>
  mutate(sar_residuos = residuals(sar_modelo),
         car_residuos = residuals(car_modelo))

Após criar essa nova base de dados, você pode exportá-la como geopackage ou shapefile com a função write_sf().

write_sf(agua_rede_sar_car, "dados/agua_rede_sar_car.gpkg")

write_sf(agua_rede_sar_car, "dados/agua_rede_sar_car.shp")

8.6 Modelo de regressão espacial local (GWR)

O modelo de regressão espacial local ou Geographically Weighted Regression (GWR) incorpora a estrutura de dependência espacial (verificada com o teste de autocorrelação espacial - Roteiro 6).

Nesse modelo, as variações espaciais são modeladas de forma contínua, com parâmetros variando no espaço. Ele ajusta um modelo de regressão para cada observação, ponderando todas as demais observações como função da distância a este ponto. Existe uma função (kernel) sobre cada ponto do espaço que determina todos os pontos da regressão local que é ponderada pela distância. Pontos mais próximos do ponto central tem maior peso. Assim como no kernel, a escolha da largura da banda é importante, pondendo ser fixa ou adaptável à densidade dos dados.

O processo de modelagem envolverá quatro etapas: (1) estimativa do kernel, (2) cômputo do modelo, (3) análise do sumário do modelo e (4) exportação dos resultados para análise espacial no QGIS.

É possível obter a melhor estimativa do kernel com a função gwr.sel(). Como argumentos, indicamos a fórmula, os dados, as coordenadas dos centróides e indicamos que a largura da banda será adaptativa (adapt = TRUE). O resultado será salvo no objeto gwr_kernel.

gwr_kernel <- gwr.sel(formula = CONSUMO1 ~ RENDAPITA, 
                        data = agua_rede_sf_sp_centroide,
                        coords = cbind(agua_rede_sf_sp_centroide$LON,
                                       agua_rede_sf_sp_centroide$LAT),
                        adapt = TRUE)

Adaptive q: 0.381966 CV score: 102015.2 
Adaptive q: 0.618034 CV score: 105439.5 
Adaptive q: 0.236068 CV score: 98555.9 
Adaptive q: 0.145898 CV score: 95869.4 
Adaptive q: 0.09016994 CV score: 93855.01 
Adaptive q: 0.05572809 CV score: 92954.56 
Adaptive q: 0.03444185 CV score: 92782.26 
Adaptive q: 0.03259011 CV score: 92883.41 
Adaptive q: 0.04359183 CV score: 92738.4 
Adaptive q: 0.04822747 CV score: 92801.64 
Adaptive q: 0.04080911 CV score: 92746.26 
Adaptive q: 0.04536248 CV score: 92721.9 
Adaptive q: 0.04645681 CV score: 92737.45 
Adaptive q: 0.04504466 CV score: 92718.97 
Adaptive q: 0.04448972 CV score: 92726.54 
Adaptive q: 0.04500397 CV score: 92719.53 
Adaptive q: 0.0451316 CV score: 92718.7 
Adaptive q: 0.04521979 CV score: 92719.92 
Adaptive q: 0.0451723 CV score: 92719.26 
Adaptive q: 0.04509091 CV score: 92718.33 
Adaptive q: 0.04509091 CV score: 92718.33

Após obter a melhor estimativa para o kernel, é possível fazer a modelagem com a função gwr(). Além da fórmula, dados e coordenadas, é necessário indicar o kernel (adapt = gwr_kernel) e dois argumentos adicionais para salvar os resultados completos.

gwr_modelo <- gwr(formula = CONSUMO1 ~ RENDAPITA,
                data = agua_rede_sf_sp_centroide,
                coords = cbind(agua_rede_sf_sp_centroide$LON,
                               agua_rede_sf_sp_centroide$LAT),
                adapt = gwr_kernel, 
                hatmatrix = TRUE,
                se.fit = TRUE)

É possível visualizar o sumário com alguns resultados executando como comando o nome do objeto onde o modelo foi salvo, gwr_modelo.

gwr_modelo

Call:
gwr(formula = CONSUMO1 ~ RENDAPITA, data = agua_rede_sf_sp_centroide, 
    coords = cbind(agua_rede_sf_sp_centroide$LON, agua_rede_sf_sp_centroide$LAT), 
    adapt = gwr_kernel, hatmatrix = TRUE, se.fit = TRUE)
Kernel function: gwr.Gauss 
Adaptive quantile: 0.04509091 (about 22 of 510 data points)
Summary of GWR coefficient estimates at data points:
                  Min.   1st Qu.    Median   3rd Qu.      Max.  Global
X.Intercept.  3.193279 16.736881 21.030497 26.436066 37.206325 20.5236
RENDAPITA     0.020480  0.029492  0.035011  0.041511  0.057511  0.0369
Number of data points: 510 
Effective number of parameters (residual: 2traceS - traceS'S): 32.04562 
Effective degrees of freedom (residual: 2traceS - traceS'S): 477.9544 
Sigma (residual: 2traceS - traceS'S): 13.39139 
Effective number of parameters (model: traceS): 22.38185 
Effective degrees of freedom (model: traceS): 487.6182 
Sigma (model: traceS): 13.25803 
Sigma (ML): 12.96385 
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 4109.837 
AIC (GWR p. 96, eq. 4.22): 4083.107 
Residual sum of squares: 85711.26 
Quasi-global R2: 0.3793284

Duas estatísticas muito importantes sobre o modelo são o AIC e o R². Neste caso o AIC do modelo de regressão espacial local (4083) é menor que o AIC do modelo de regressão linear (4188), indicando um melhor resultado. O R² “quasi-global” (0,37) é maior que o R² do modelo de regressão linear (0,21), também indicando um melhor resultado.

Também é importante verificar os parâmetros locais estimados e o coeficiente local de determinação (R²). Para visualizar esses parâmetros, vamos exportar uma nova base de dados (em formato geopackage ou shapefile) para abrir e visualizar no QGIS.

O código abaixo indica como executar essa última etapa.

agua_rede_gwr <- cbind(agua_rede_sf_sp_centroide,
                       as.matrix(as.data.frame(gwr_modelo$SDF))) |>
  mutate(t_beta1 = RENDAPITA.1 / RENDAPITA_se)

As novas variáveis são:

X.Intercept = valores locais de beta zero (intercepto)
RENDAPITA.1 = valores locais do coeficiente de renda per capita (X1)
X.Intercept_se = erro padrão de beta zero (intercepto)
RENDAPITA_se = erro padrão do coeficiente estimado para renda per capita (X1)
t_beta1 = estatística t local para o coeficiente de renda per capita
pred = valores previstos de Y
localR2 = valores locais de R2

Para interpretar os resultados do GWR, é importante mapear algumas estatísticas, tais como o R2 local, os coeficientes beta e estatísticas t.

Para obter a estatística t para cada coeficiente beta, deve-se dividir os valores dos coeficientes (por exemplo, RENDAPITA.1) por seus respectivos erros padrão (RENDAPITA_se). Se o módulo da estatística t for maior do que o valor do t crítico, pode-se rejeitar a hipótese nula de que beta é igual a zero (ou seja, beta é significativo).

A variável “pred” apresenta os parâmetros locais estimados para a variável Renda per capita, enquanto a variável “localR2” apresenta o coeficiente local de determinação (R²). Também é interessante aplicar um teste de autocorrelação espacial nos resíduos, para verificar se a dependência espacial foi resolvida pelo modelo.

Para exportar como geopackage:

write_sf(agua_rede_gwr, "dados/agua_rede_gwr.gpkg")

Para exportar como shapefile:

write_sf(agua_rede_gwr, "dados/agua_rede_gwr.shp")

8.7 Sugestões de materiais

Dados espaciais no R - material do curso “R intermediário e pesquisa reprodutível”