Análise espacial do crescimento populacional de Pernambuco

Introdução

Esta análise aborda o crescimento populacional em Pernambuco entre os censos de 2010 e 2022. O objetivo é demonstrar a habilidade de trabalhar com dados espaciais, processá-los e apresentá-los de maneira eficiente para fins de análise.

Primeiros passos

Para começar a análise, vamos instalar e carregar os pacotes necessários. A função abaixo verifica se os pacotes estão instalados e faz o download dos que ainda não estão. Em seguida, carrega-os. Outra alternativa seria usar a função “require”.

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

# Instalar pacotes se não estiverem instalados
packages <- c("sf", "spdep", "tmap", "dplyr", "readr", "data.table", "openxlsx", "tmap", "classInt", "spdep", "spatialreg", "rsconnect")
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)

# Carregar pacotes
library(sf)
library(spdep)
library(tmap)
library(dplyr)
library(readr)
library(data.table)
library(openxlsx)
library(tmap)
library(classInt)
library(spdep)
library(spatialreg)
library(rsconnect)

Carregar dados

Utilizei os dados do censo geográfico de 2010 e os dados preliminares do censo de 2022, ambos disponíveis em https://www.ibge.gov.br/estatisticas/downloads-estatisticas.html. É possível que o IBGE disponibilize esses dados via API. Porém, optei por salvar os arquivos em CSV para simplificar a análise pois o objetivo deste relatório é demonstrar a capacidade de tratamento dos dados.

Além disso, utilizamos a malha geográfica divulgada pelo IBGE junto com o último censo de 2022. Este arquivo JSON contém as coordenadas dos polígonos dos setores geográficos.

# Caminhos para os arquivos de dados
shapefile_path <- "C:\\Users\\mallu\\Documents\\estudos\\PE_Malha_Preliminar_2022.json"
censo_2010_path <- "C:\\Users\\mallu\\Documents\\estudos\\censo_2010.csv"
censo_2022_path <- "C:\\Users\\mallu\\Documents\\estudos\\censo_2022.csv"

# Carregar malha
pernambuco_shapefile <- st_read(shapefile_path)

## Reading layer `PE_Malha_Preliminar_2022' from data source 
##   `C:\Users\mallu\Documents\estudos\PE_Malha_Preliminar_2022.json' 
##   using driver `GeoJSON'
## Simple feature collection with 17574 features and 22 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -41.35834 ymin: -9.482897 xmax: -32.37777 ymax: -3.804622
## Geodetic CRS:  WGS 84

print("Shapefile carregado com sucesso!")

## [1] "Shapefile carregado com sucesso!"

Os dados em CSV possuem algumas limitações de codificação. Por exemplo, os nomes dos municípios podem apresentar erros de ortografia devido a problemas de leitura de caracteres especiais. A função abaixo tenta ler os arquivos CSV de diferentes maneiras para evitar esses problemas. Além disso, corrige erros ao identificá-los.

# Função para tentar ler CSV com diferentes codificações e delimitadores
try_read_csv <- function(file_path, delim = ";") {
  encodings <- c("UTF-8", "latin1", "ASCII")
  for (encoding in encodings) {
    try({
      data <- fread(file_path, encoding = encoding, sep = delim, fill = TRUE)
      print(paste("Carregado com sucesso usando a codificação:", encoding))
      return(data)
    }, silent = TRUE)
  }
  stop("Falha ao carregar o arquivo com as codificações testadas.")
}

# Função para corrigir codificação de caracteres especiais
correct_encoding <- function(data) {
  data <- data.frame(lapply(data, function(x) {
    if (is.character(x)) {
      return(iconv(x, from = "latin1", to = "UTF-8", sub = ""))
    }
    return(x)
  }), stringsAsFactors = FALSE)
  return(data)
}

# Função específica para corrigir nomes de colunas
correct_column_names <- function(data) {
  colnames(data) <- iconv(colnames(data), from = "latin1", to = "UTF-8", sub = "")
  return(data)
}

# Carregar dados do censo de 2010 e corrigir codificação
try_read_csv_correct_encoding <- function(file_path, delim = ";") {
  data <- try_read_csv(file_path, delim)
  # Corrigir nomes das colunas
  data <- correct_column_names(data)
  # Corrigir codificação de caracteres especiais
  data <- correct_encoding(data)
  return(data)
}

# Carregar dados do censo de 2010 com correção de codificação
censo_2010 <- try_read_csv_correct_encoding(censo_2010_path)

## [1] "Carregado com sucesso usando a codificação: UTF-8"

censo_2022 <- try_read_csv(censo_2022_path)

## [1] "Carregado com sucesso usando a codificação: UTF-8"

# Verificar estrutura dos dados
str(pernambuco_shapefile)

## Classes 'sf' and 'data.frame':   17574 obs. of  23 variables:
##  $ CD_SETOR  : chr  "260005405000001P" "260005405000002P" "260005405000003P" "260005405000004P" ...
##  $ AREA_KM2  : num  0.0682 0.0528 0.0663 0.084 0.0575 ...
##  $ CD_REGIAO : chr  "2" "2" "2" "2" ...
##  $ NM_REGIAO : chr  "Nordeste" "Nordeste" "Nordeste" "Nordeste" ...
##  $ CD_UF     : chr  "26" "26" "26" "26" ...
##  $ NM_UF     : chr  "Pernambuco" "Pernambuco" "Pernambuco" "Pernambuco" ...
##  $ CD_MUN    : chr  "2600054" "2600054" "2600054" "2600054" ...
##  $ NM_MUN    : chr  "Abreu e Lima" "Abreu e Lima" "Abreu e Lima" "Abreu e Lima" ...
##  $ CD_DIST   : chr  "260005405" "260005405" "260005405" "260005405" ...
##  $ NM_DIST   : chr  "Abreu e Lima" "Abreu e Lima" "Abreu e Lima" "Abreu e Lima" ...
##  $ CD_SUBDIST: chr  "26000540500" "26000540500" "26000540500" "26000540500" ...
##  $ NM_SUBDIST: chr  NA NA NA NA ...
##  $ CD_MICRO  : chr  "26017" "26017" "26017" "26017" ...
##  $ NM_MICRO  : chr  "Recife" "Recife" "Recife" "Recife" ...
##  $ CD_MESO   : chr  "2605" "2605" "2605" "2605" ...
##  $ NM_MESO   : chr  "Metropolitana de Recife" "Metropolitana de Recife" "Metropolitana de Recife" "Metropolitana de Recife" ...
##  $ CD_RGI    : chr  "260001" "260001" "260001" "260001" ...
##  $ NM_RGI    : chr  "Recife" "Recife" "Recife" "Recife" ...
##  $ CD_RGINT  : chr  "2601" "2601" "2601" "2601" ...
##  $ NM_RGINT  : chr  "Recife" "Recife" "Recife" "Recife" ...
##  $ CD_CONCURB: chr  "2611606" "2611606" "2611606" "2611606" ...
##  $ NM_CONCURB: chr  "Recife/PE" "Recife/PE" "Recife/PE" "Recife/PE" ...
##  $ geometry  :sfc_MULTIPOLYGON of length 17574; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:37, 1:2] -34.9 -34.9 -34.9 -34.9 -34.9 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:22] "CD_SETOR" "AREA_KM2" "CD_REGIAO" "NM_REGIAO" ...

str(censo_2010)

## 'data.frame':    12379 obs. of  33 variables:
##  $ Cod_setor          :integer64 260005405000001 260005405000002 260005405000003 260005405000004 260005405000005 260005405000006 260005405000007 260005405000008 ... 
##  $ Cod_Grandes.Regiões: int  2 2 2 2 2 2 2 2 2 2 ...
##  $ Nome_Grande_Regiao : chr  "Região Nordeste" "Região Nordeste" "Região Nordeste" "Região Nordeste" ...
##  $ Cod_UF             : int  26 26 26 26 26 26 26 26 26 26 ...
##  $ Nome_da_UF         : chr  "Pernambuco" "Pernambuco" "Pernambuco" "Pernambuco" ...
##  $ Cod_meso           : int  2605 2605 2605 2605 2605 2605 2605 2605 2605 2605 ...
##  $ Nome_da_meso       : chr  "Metropolitana de Recife" "Metropolitana de Recife" "Metropolitana de Recife" "Metropolitana de Recife" ...
##  $ Cod_micro          : int  26017 26017 26017 26017 26017 26017 26017 26017 26017 26017 ...
##  $ Nome_da_micro      : chr  "Recife" "Recife" "Recife" "Recife" ...
##  $ Cod_RM             : int  11 11 11 11 11 11 11 11 11 11 ...
##  $ Nome_da_RM         : chr  "RM Recife" "RM Recife" "RM Recife" "RM Recife" ...
##  $ Cod_municipio      : int  2600054 2600054 2600054 2600054 2600054 2600054 2600054 2600054 2600054 2600054 ...
##  $ Nome_do_municipio  : chr  "ABREU E LIMA" "ABREU E LIMA" "ABREU E LIMA" "ABREU E LIMA" ...
##  $ Cod_distrito       : int  260005405 260005405 260005405 260005405 260005405 260005405 260005405 260005405 260005405 260005405 ...
##  $ Nome_do_distrito   : chr  "ABREU E LIMA" "ABREU E LIMA" "ABREU E LIMA" "ABREU E LIMA" ...
##  $ Cod_subdistrito    :integer64 26000540500 26000540500 26000540500 26000540500 26000540500 26000540500 26000540500 26000540500 ... 
##  $ Nome_do_subdistrito: chr  "ABREU E LIMA" "ABREU E LIMA" "ABREU E LIMA" "ABREU E LIMA" ...
##  $ Cod_bairro         :integer64 2600054006 2600054006 2600054006 2600054006 2600054002 2600054006 2600054008 2600054001 ... 
##  $ Nome_do_bairro     : chr  "Centro" "Centro" "Centro" "Centro" ...
##  $ Situacao_setor     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ V001               : int  211 223 140 318 258 301 184 457 305 354 ...
##  $ V002               : int  724 750 417 971 890 1012 623 1556 976 1116 ...
##  $ V003               : chr  "3,43" "3,36" "2,98" "3,05" ...
##  $ V004               : chr  "2,24" "2,11" "2,12" "2,07" ...
##  $ V005               : chr  "452,8" "632,62" "973,5" "637,76" ...
##  $ V006               : chr  "369490,04" "260188,28" "1573094,48" "396666,92" ...
##  $ V007               : chr  "764,32" "727,19" "1117,13" "765,31" ...
##  $ V008               : chr  "385725,2" "230159,08" "1645327,64" "378314,62" ...
##  $ V009               : chr  "322,3" "386,56" "623,45" "382,99" ...
##  $ V010               : chr  "266389,59" "229024,17" "958614,38" "295420,17" ...
##  $ V011               : chr  "693,75" "686,42" "1005,48" "658,24" ...
##  $ V012               : chr  "315881,23" "200776,62" "1162765,85" "326616,9" ...
##  $ V33                : logi  NA NA NA NA NA NA ...

str(censo_2022)

## Classes 'data.table' and 'data.frame':   17574 obs. of  29 variables:
##  $ CD_SETOR  : chr  "260005405000001P" "260005405000002P" "260005405000003P" "260005405000004P" ...
##  $ AREA_KM2  : num  0.0682 0.0528 0.0663 0.084 0.0575 ...
##  $ CD_REGIAO : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ NM_REGIAO : chr  "Nordeste" "Nordeste" "Nordeste" "Nordeste" ...
##  $ CD_UF     : int  26 26 26 26 26 26 26 26 26 26 ...
##  $ NM_UF     : chr  "Pernambuco" "Pernambuco" "Pernambuco" "Pernambuco" ...
##  $ CD_MUN    : int  2600054 2600054 2600054 2600054 2600054 2600054 2600054 2600054 2600054 2600054 ...
##  $ NM_MUN    : chr  "Abreu e Lima" "Abreu e Lima" "Abreu e Lima" "Abreu e Lima" ...
##  $ CD_DIST   : int  260005405 260005405 260005405 260005405 260005405 260005405 260005405 260005405 260005405 260005405 ...
##  $ NM_DIST   : chr  "Abreu e Lima" "Abreu e Lima" "Abreu e Lima" "Abreu e Lima" ...
##  $ CD_SUBDIST:integer64 26000540500 26000540500 26000540500 26000540500 26000540500 26000540500 26000540500 26000540500 ... 
##  $ NM_SUBDIST: chr  "Subdistrito sem denominação" "Subdistrito sem denominação" "Subdistrito sem denominação" "Subdistrito sem denominação" ...
##  $ CD_MICRO  : int  26017 26017 26017 26017 26017 26017 26017 26017 26017 26017 ...
##  $ NM_MICRO  : chr  "Recife" "Recife" "Recife" "Recife" ...
##  $ CD_MESO   : int  2605 2605 2605 2605 2605 2605 2605 2605 2605 2605 ...
##  $ NM_MESO   : chr  "Metropolitana de Recife" "Metropolitana de Recife" "Metropolitana de Recife" "Metropolitana de Recife" ...
##  $ CD_RGI    : int  260001 260001 260001 260001 260001 260001 260001 260001 260001 260001 ...
##  $ NM_RGI    : chr  "Recife" "Recife" "Recife" "Recife" ...
##  $ CD_RGINT  : int  2601 2601 2601 2601 2601 2601 2601 2601 2601 2601 ...
##  $ NM_RGINT  : chr  "Recife" "Recife" "Recife" "Recife" ...
##  $ CD_CONCURB: int  2611606 2611606 2611606 2611606 2611606 2611606 2611606 2611606 2611606 2611606 ...
##  $ NM_CONCURB: chr  "Recife/PE" "Recife/PE" "Recife/PE" "Recife/PE" ...
##  $ v0001     : int  627 589 278 902 828 734 649 778 922 865 ...
##  $ v0002     : int  256 240 157 409 339 325 255 348 400 372 ...
##  $ v0003     : int  256 240 157 409 339 325 255 347 399 372 ...
##  $ v0004     : int  0 0 0 0 0 0 0 1 1 0 ...
##  $ v0005     : num  2.85 2.82 2.28 2.64 2.79 ...
##  $ v0006     : num  12.73 3.83 6.56 5.85 1.68 ...
##  $ v0007     : int  220 209 122 342 297 274 230 295 342 313 ...
##  - attr(*, ".internal.selfref")=<externalptr>

# Identificar problemas ao ler os dados (se necessário)
if (exists("problems")) {
  problems_2010 <- problems(censo_2010)
  problems_2022 <- problems(censo_2022)
  
  # Mostrar problemas encontrados
  print(problems_2010)
  print(problems_2022)
} else {
  print("Função problems não está disponível.")
}

## # A tibble: 0 × 4
## # ℹ 4 variables: row <int>, col <int>, expected <chr>, actual <chr>
## # A tibble: 0 × 4
## # ℹ 4 variables: row <int>, col <int>, expected <chr>, actual <chr>

O IBGE divulga os dados para os setores geográficos, que possuem uma dimensão menor que os municípios. Utilizamos o código dos municípios para agregar esses dados, assumindo que a população dos municípios é a soma das populações dos setores que os compõem. Além disso, fiz alguns ajustes nos nomes das colunas para facilitar a análise, mantendo apenas as colunas de interesse, como dados sobre a população residente, o código e o nome dos municípios, entre outros. Vale ressaltar que os dados de população residente estão codificados como v0001 no censo de 2022 e V002 no censo de 2010. Renomeei esses valores para “censo_2022” e “censo_2010”, respectivamente.

# Agregar dados por município para o censo de 2010
censo_2010_agg <- censo_2010 %>%
  group_by(Cod_municipio) %>%
  summarise(censo_2010 = sum(V002, na.rm = TRUE))

# Agregar dados por município para o censo de 2022
censo_2022_agg <- censo_2022 %>%
  group_by(CD_MUN, NM_MUN, NM_MESO) %>%
  summarise(censo_2022 = sum(v0001, na.rm = TRUE))

# Renomear colunas para facilitar o merge
censo_2022_agg <- censo_2022_agg %>%
  rename(Cod_municipio = CD_MUN)

Em seguida, vamos combinar os dados em um único dataframe e calcular a variação percentual da população residente entre os censos.

# Combinar dados do censo de 2010 e 2022
censo_combined <- merge(censo_2010_agg, censo_2022_agg, by = c("Cod_municipio"))

# Calcular variação percentual da população
censo_combined<- censo_combined %>%
  mutate(variacao_percentual = ((censo_2022 - censo_2010) / censo_2010) * 100)

# Verificar combinação dos dados
print("Dados combinados:")

## [1] "Dados combinados:"

print(head(censo_combined))

##   Cod_municipio censo_2010                NM_MUN                    NM_MESO
## 1       2600054      93892          Abreu e Lima    Metropolitana de Recife
## 2       2600104      35018 Afogados da Ingazeira        Sertão Pernambucano
## 3       2600203      17558               Afrânio São Francisco Pernambucano
## 4       2600302      22645             Agrestina       Agreste Pernambucano
## 5       2600401      32427            Água Preta          Mata Pernambucana
## 6       2600500      39973           Águas Belas       Agreste Pernambucano
##   censo_2022 variacao_percentual
## 1      98462            4.867294
## 2      40241           14.915186
## 3      18674            6.356077
## 4      23779            5.007728
## 5      26461          -18.398248
## 6      41548            3.940160

O último tratamento necessário foi combinar os dados populacionais com a malha geográfica do arquivo JSON. Isso nos permitirá construir mapas e realizar análises de estatística espacial. Adicionei algumas partes do código para validação, a fim de garantir que todos os municípios foram cruzados corretamente com a malha e que não há valores faltantes (“NAs”). Os limites dos municípios são dados pela união dos limites externos dos setores que os compõem.

# Combinar dados agregados com a malha geográfica
pernambuco_shapefile <- pernambuco_shapefile %>%
  mutate(Cod_municipio = as.character(CD_MUN)) %>%
  group_by(Cod_municipio) %>%
  summarise(geometry = st_union(geometry))

# Unir com os dados agregados
municipios_shapefile <- merge(pernambuco_shapefile, censo_combined, by = "Cod_municipio")

# Remover linhas com valores NA na coluna variacao_percentual
municipios_shapefile <- municipios_shapefile[!is.na(municipios_shapefile$variacao_percentual), ]

# Verificar a quantidade de valores nulos restantes
sum(is.na(municipios_shapefile$variacao_percentual))

## [1] 0

Mapas e análises visuais

A seguir, apresento os mapas populacionais de Pernambuco. Mantive a mesma escala para os mapas de 2010 e 2022, permitindo uma comparação visual direta entre os dois períodos. Para essa visualização, utilizei o pacote “tmap”.

tmap_mode("view")

# Obter os intervalos kmeans do mapa de 2010
breaks_2010 <- classIntervals(municipios_shapefile$censo_2010, n = 7, style = "kmeans")$brks

# Plotar o mapa para o Censo 2010
tm_shape(municipios_shapefile) +
  tm_polygons("censo_2010", style = "kmeans", n = 7, palette = "Blues", title = "Censo 2010",
              popup.vars = c("Nome do Município" = "NM_MUN", "População 2010" = "censo_2010")) +
  tm_layout(title = "População do Censo 2010 em Pernambuco",
            legend.outside = TRUE, frame = FALSE, inner.margins = c(0, 0, 0, 0),
            outer.margins = c(0, 0, 0, 0), asp = 0.0)

# Usar os mesmos intervalos para o mapa de 2022
tm_shape(municipios_shapefile) +
  tm_polygons("censo_2022", breaks = breaks_2010, palette = "Blues", title = "Censo 2022",
              popup.vars = c("Nome do Município" = "NM_MUN", "População 2022" = "censo_2022")) +
  tm_layout(title = "População do Censo 2022 em Pernambuco",
            legend.outside = TRUE, frame = FALSE, inner.margins = c(0, 0, 0, 0),
            outer.margins = c(0, 0, 0, 0), asp = 0.0)

A seguir, apresento o crescimento populacional, em pontos percentuais, entre 2010 e 2022 no estado de Pernambuco.

tm_shape(municipios_shapefile) +
  tm_polygons("variacao_percentual", style = "jenks", n = 7, palette = "Blues", title = "Variação Percentual",
              popup.vars = c("Nome do Município" = "NM_MUN", "Variação Percentual" = "variacao_percentual")) +
  tm_layout(title = "Variação Percentual da População em Pernambuco (2010-2022)",
            legend.outside = TRUE, frame = FALSE, inner.margins = c(0, 0, 0, 0),
            outer.margins = c(0, 0, 0, 0), asp = 0.0)

Visualmente, parece que existem algumas relações espaciais importantes. Os municípios com maior crescimento populacional parecem ter vizinhos com crescimento intermediário. Além disso, há uma formação de concentrações populacionais em torno de cidades populosas em ambos os censos, como no caso da região metropolitana do Recife. Vamos selecionar os dados de crescimento populacional para uma análise rápida utilizando estatística espacial.

Matriz de pesos espaciais

Primeiro, tratamos os conjuntos de vizinhos vazios para garantir que cada município seja considerado vizinho de si mesmo, essencial para incluir municípios isolados como Fernando de Noronha na análise. Em seguida, criamos a matriz de pesos espaciais usando nb2listw com o estilo “W”, que normaliza os pesos para que a soma dos pesos de cada município seja 1. A verificação da matriz de pesos espaciais mostra que há dois componentes desconexos: o componente 1, que contém a maioria dos municípios de Pernambuco, e o componente 2, que contém apenas Fernando de Noronha, como esperado e já tratado em nossa pré-análise.

# Função para tratar conjuntos de vizinhos vazios
nb_self <- include.self(poly2nb(municipios_shapefile))

# Criar a matriz de pesos espaciais
lw <- nb2listw(nb_self, style="W")

# Verificar a matriz de pesos espaciais
print(lw)

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 185 
## Number of nonzero links: 1135 
## Percentage nonzero weights: 3.316289 
## Average number of links: 6.135135 
## 2 disjoint connected subgraphs
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1       S2
## W 185 34225 185 65.62513 755.4484

# Identificar os componentes conectados
components <- n.comp.nb(nb_self)

# Listar os nomes dos municípios em cada componente desconexo
for (i in 1:components$nc) {
  cat("Componente", i, "contém os seguintes municípios:\n")
  print(municipios_shapefile$NM_MUN[which(components$comp.id == i)])
}

## Componente 1 contém os seguintes municípios:
##   [1] "Abreu e Lima"              "Afogados da Ingazeira"    
##   [3] "Afrânio"                   "Agrestina"                
##   [5] "Água Preta"                "Águas Belas"              
##   [7] "Alagoinha"                 "Aliança"                  
##   [9] "Altinho"                   "Amaraji"                  
##  [11] "Angelim"                   "Araçoiaba"                
##  [13] "Araripina"                 "Arcoverde"                
##  [15] "Barra de Guabiraba"        "Barreiros"                
##  [17] "Belém de Maria"            "Belém do São Francisco"   
##  [19] "Belo Jardim"               "Betânia"                  
##  [21] "Bezerros"                  "Bodocó"                   
##  [23] "Bom Conselho"              "Bom Jardim"               
##  [25] "Bonito"                    "Brejão"                   
##  [27] "Brejinho"                  "Brejo da Madre de Deus"   
##  [29] "Buenos Aires"              "Buíque"                   
##  [31] "Cabo de Santo Agostinho"   "Cabrobó"                  
##  [33] "Cachoeirinha"              "Caetés"                   
##  [35] "Calçado"                   "Calumbi"                  
##  [37] "Camaragibe"                "Camocim de São Félix"     
##  [39] "Camutanga"                 "Canhotinho"               
##  [41] "Capoeiras"                 "Carnaíba"                 
##  [43] "Carnaubeira da Penha"      "Carpina"                  
##  [45] "Caruaru"                   "Casinhas"                 
##  [47] "Catende"                   "Cedro"                    
##  [49] "Chã de Alegria"            "Chã Grande"               
##  [51] "Condado"                   "Correntes"                
##  [53] "Cortês"                    "Cumaru"                   
##  [55] "Cupira"                    "Custódia"                 
##  [57] "Dormentes"                 "Escada"                   
##  [59] "Exu"                       "Feira Nova"               
##  [61] "Ferreiros"                 "Flores"                   
##  [63] "Floresta"                  "Frei Miguelinho"          
##  [65] "Gameleira"                 "Garanhuns"                
##  [67] "Glória do Goitá"           "Goiana"                   
##  [69] "Granito"                   "Gravatá"                  
##  [71] "Iati"                      "Ibimirim"                 
##  [73] "Ibirajuba"                 "Igarassu"                 
##  [75] "Iguaracy"                  "Inajá"                    
##  [77] "Ingazeira"                 "Ipojuca"                  
##  [79] "Ipubi"                     "Itacuruba"                
##  [81] "Itaíba"                    "Ilha de Itamaracá"        
##  [83] "Itambé"                    "Itapetim"                 
##  [85] "Itapissuma"                "Itaquitinga"              
##  [87] "Jaboatão dos Guararapes"   "Jaqueira"                 
##  [89] "Jataúba"                   "Jatobá"                   
##  [91] "João Alfredo"              "Joaquim Nabuco"           
##  [93] "Jucati"                    "Jupi"                     
##  [95] "Jurema"                    "Lagoa do Carro"           
##  [97] "Lagoa de Itaenga"          "Lagoa do Ouro"            
##  [99] "Lagoa dos Gatos"           "Lagoa Grande"             
## [101] "Lajedo"                    "Limoeiro"                 
## [103] "Macaparana"                "Machados"                 
## [105] "Manari"                    "Maraial"                  
## [107] "Mirandiba"                 "Moreno"                   
## [109] "Nazaré da Mata"            "Olinda"                   
## [111] "Orobó"                     "Orocó"                    
## [113] "Ouricuri"                  "Palmares"                 
## [115] "Palmeirina"                "Panelas"                  
## [117] "Paranatama"                "Parnamirim"               
## [119] "Passira"                   "Paudalho"                 
## [121] "Paulista"                  "Pedra"                    
## [123] "Pesqueira"                 "Petrolândia"              
## [125] "Petrolina"                 "Poção"                    
## [127] "Pombos"                    "Primavera"                
## [129] "Quipapá"                   "Quixaba"                  
## [131] "Recife"                    "Riacho das Almas"         
## [133] "Ribeirão"                  "Rio Formoso"              
## [135] "Sairé"                     "Salgadinho"               
## [137] "Salgueiro"                 "Saloá"                    
## [139] "Sanharó"                   "Santa Cruz"               
## [141] "Santa Cruz da Baixa Verde" "Santa Cruz do Capibaribe" 
## [143] "Santa Filomena"            "Santa Maria da Boa Vista" 
## [145] "Santa Maria do Cambucá"    "Santa Terezinha"          
## [147] "São Benedito do Sul"       "São Bento do Una"         
## [149] "São Caitano"               "São João"                 
## [151] "São Joaquim do Monte"      "São José da Coroa Grande" 
## [153] "São José do Belmonte"      "São José do Egito"        
## [155] "São Lourenço da Mata"      "São Vicente Férrer"       
## [157] "Serra Talhada"             "Serrita"                  
## [159] "Sertânia"                  "Sirinhaém"                
## [161] "Moreilândia"               "Solidão"                  
## [163] "Surubim"                   "Tabira"                   
## [165] "Tacaimbó"                  "Tacaratu"                 
## [167] "Tamandaré"                 "Taquaritinga do Norte"    
## [169] "Terezinha"                 "Terra Nova"               
## [171] "Timbaúba"                  "Toritama"                 
## [173] "Tracunhaém"                "Trindade"                 
## [175] "Triunfo"                   "Tupanatinga"              
## [177] "Tuparetama"                "Venturosa"                
## [179] "Verdejante"                "Vertente do Lério"        
## [181] "Vertentes"                 "Vicência"                 
## [183] "Vitória de Santo Antão"    "Xexéu"                    
## Componente 2 contém os seguintes municípios:
## [1] "Fernando de Noronha"

A verificação da matriz de pesos espaciais revelou 185 regiões e 1135 links não nulos, com um percentual de pesos não nulos de 3,32% e uma média de 6,14 links por município.

Explorar autocorrelação espacial

moran_test <- moran.test(municipios_shapefile$variacao_percentual, lw)
moran_test

## 
##  Moran I test under randomisation
## 
## data:  municipios_shapefile$variacao_percentual  
## weights: lw    
## 
## Moran I statistic standard deviate = 10.881, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.460156427      -0.005434783       0.001830902

Executamos o teste de Moran para avaliar a autocorrelação espacial da variação percentual da população entre os municípios de Pernambuco. O teste utiliza a matriz de pesos espaciais previamente calculada (lw). O resultado mostra um valor estatístico de Moran I de 0.460, significativamente maior do que o esperado sob a hipótese nula de ausência de autocorrelação (valor esperado de -0.005). Com um desvio padrão de 10.881 e um p-valor muito pequeno (< 2.2e-16), rejeitamos a hipótese nula, indicando uma autocorrelação espacial positiva significativa. Isso sugere que municípios com crescimentos populacionais semelhantes tendem a estar geograficamente próximos, revelando padrões espaciais no crescimento populacional de Pernambuco.

Ajustar modelos espaciais

slm <- lagsarlm(variacao_percentual ~ 1, data = municipios_shapefile, listw = lw)
summary(slm)

## 
## Call:lagsarlm(formula = variacao_percentual ~ 1, data = municipios_shapefile, 
##     listw = lw)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -34.57724  -5.28365  -0.34518   4.26225  30.41043 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.47886    0.65738  0.7284   0.4663
## 
## Rho: 0.59505, LR test value: 39.371, p-value: 3.5045e-10
## Asymptotic standard error: 0.073551
##     z-value: 8.0903, p-value: 6.6613e-16
## Wald statistic: 65.453, p-value: 5.5511e-16
## 
## Log likelihood: -694.5115 for lag model
## ML residual variance (sigma squared): 78.548, (sigma: 8.8627)
## Number of observations: 185 
## Number of parameters estimated: 3 
## AIC: 1395, (AIC for lm: 1432.4)
## LM test for residual autocorrelation
## test value: -72.403, p-value: 1

Executamos um modelo SAR (Spatial Autoregressive Model) usando a função lagsarlm para analisar a variação percentual da população nos municípios de Pernambuco. O modelo considera a dependência espacial através da matriz de pesos espaciais (lw). O coeficiente de autocorrelação espacial (Rho) é 0.595, com um valor z de 8.0903 e um p-valor extremamente pequeno (6.6613e-16), indicando uma forte autocorrelação espacial positiva. O teste de razão de verossimilhança (LR test) também é significativo (p-valor: 3.5045e-10), reforçando a importância da autocorrelação espacial no modelo. O AIC (Akaike Information Criterion) para o modelo SAR é 1395, menor que o AIC do modelo linear simples (1432.4), sugerindo que o modelo SAR é mais adequado para os dados. O teste LM para autocorrelação residual não é significativo, indicando que a autocorrelação espacial foi adequadamente capturada pelo modelo. Esses resultados mostram que a variação percentual da população entre os municípios está fortemente influenciada por fatores espaciais.

Outras possibilidades de análise

Outra possibilidade de análise é realizar agregações diferentes, como por mesorregiões. Neste caso, repetimos o processo de agregação para criar os limites das mesorregiões. Visualmente, percebemos um aumento populacional maior nas mesorregiões interioranas, especialmente na região do Vale do São Francisco, que na última década se tornou uma importante área produtora do estado. Para trabalhos futuros, análises estatísticas podem cruzar dados como renda, PIB per capita, educação, disponibilidade habitacional, entre outros, para tentar explicar esse movimento de forma mais fundamentada.Entretanto, ressalta-se que mais alguns tratamentos serão necessários, pois o crescimento não foi homogêneo nos municípios de cada região e outliers podem estar elevando significativamente a média, como Petrolina no Vale do São Francisco.

# Agregar os municípios nas mesorregiões e criar os limites das mesorregiões
mesorregioes_shapefile <- municipios_shapefile %>%
  group_by(NM_MESO) %>%
  summarise(geometry = st_union(geometry),
            censo_2010 = sum(censo_2010),
            censo_2022 = sum(censo_2022))

# Calcular variação percentual da população nas mesorregiões
mesorregioes_shapefile <- mesorregioes_shapefile %>%
  mutate(variacao_percentual = ((censo_2022 - censo_2010) / censo_2010) * 100)

# Verificar a junção final
print("Dados geográficos combinados por mesorregião:")

## [1] "Dados geográficos combinados por mesorregião:"

print(mesorregioes_shapefile)

## Simple feature collection with 5 features and 4 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: -41.35834 ymin: -9.482897 xmax: -32.37777 ymax: -3.804622
## Geodetic CRS:  WGS 84
## # A tibble: 5 × 5
##   NM_MESO                     geometry censo_2010 censo_2022 variacao_percentual
## * <chr>                 <GEOMETRY [°]>      <int>      <int>               <dbl>
## 1 Agreste P… POLYGON ((-35.59181 -8.4…    2209144    2339935                5.92
## 2 Mata Pern… MULTIPOLYGON (((-35.203 …    1303236    1259882               -3.33
## 3 Metropoli… MULTIPOLYGON (((-34.8624…    3675687    3730141                1.48
## 4 Sertão Pe… POLYGON ((-38.17867 -8.2…     993864    1051783                5.83
## 5 São Franc… POLYGON ((-38.46981 -8.8…     575949     677190               17.6

tm_shape(mesorregioes_shapefile) +
  tm_polygons("variacao_percentual", style = "jenks", n = 7, palette = "Blues", title = "Variação Percentual",
              popup.vars = c("Nome do Município" = "NM_MESO", "Variação Percentual" = "variacao_percentual")) +
  tm_layout(title = "Variação Percentual da População em Pernambuco (2010-2022)",
            legend.outside = TRUE, frame = FALSE, inner.margins = c(0, 0, 0, 0),
            outer.margins = c(0, 0, 0, 0), asp = 0.0)

Conclusão

Nossa análise revelou uma autocorrelação espacial significativa no crescimento populacional dos municípios de Pernambuco entre 2010 e 2022, indicando que municípios próximos tendem a apresentar variações populacionais semelhantes. O modelo SAR mostrou que a dependência espacial é um fator importante na distribuição do crescimento populacional. Além disso, observamos um crescimento mais acentuado nas mesorregiões interioranas, especialmente no Vale do São Francisco, uma área de crescente importância produtiva. Esses resultados sugerem que fatores espaciais e socioeconômicos específicos devem ser considerados para entender melhor os padrões de crescimento populacional no estado. Esta análise é breve e simplificada, pois o objetivo é demonstrar a capacidade de coletar, tratar e visualizar dados espaciais e gerar conclusões a partir deles.