28 Introdução ao GIS
28.1 Visão Geral
Os aspectos especiais dos seus dados podem trazer vários insights à situação de um surto, e pode responder perguntas como:
- Onde estão os hotspots da doença?
- Como os hotspots mudaram ao longo do tempo?
- Como é o acesso às unidades de saúdo? Melhorias são necessárias?
O foco dessa página sobre GIS é atender às necessidades dos epidemiologistas na resposta aos surtos. Iremos explorar métodos básicos de visualização de dados espaciais utilizando os pacotes tmap e ggplot2. Também vamos caminhar por métodos básicos de gerenciamento e consulta de dados espaciais com o pacote sf. Por fim, vamos dar uma olhada em conceitos de estatística espacial tais como relações espaciais, autocorrelação espacial e regressão espacial utilizando o pacote spdep.
28.2 Terminologia
Abaixo introduziremos a terminologia principal. Para uma introdução mais profunda ao GIS e à análise espacial, sugerimos que você revise algum dos tutoriais maiores ou cursos listados na seção de Referências.
Geographic Information System (GIS) (Sistema de Informação Geográfica) - Um GIS é um framework ou ambiente para reunir, gerenciar, analisar e visualizar dados espaciais.
Softwares GIS
Alguns Softwares GIS mais conhecidos permitem interação do tipo “apontar-e-clicar” (point-and-click) para o desenvolvimento de mapas e análises espaciais. Essas ferramentas vêm com as vantagens de não precisar aprender código e com a facilidade de selecionar e posicionar manualmente ícones e outras características em um mapa. Aqui temos dois dos mais conhecidos:
ArcGIS - Um software GIS comercial desenvolvido pela empresa ESRI, que é muito famoso mas bastante caro.
QGIS - Um software GIS gratuito e open-source que pode fazer praticamente qualquer coisa que o ArcGIS também pode. Você pode baixar o QGIS aqui
Utilizar o R como um GIS pode parecer mais intimidador no início, pois em vez do “apontar-e-clicar”, ele tem uma interface de linha de comando (você precisa escrever código para chegar no resultado esperado). No entanto, essa é uma grande vantagem se você precisa produzir mapas repetitivamente ou criar análises que sejam reprodutíveis.
Dados espaciais
As duas formas primordiais de dados espaciais utilizados no GIS são dados vetoriais e rasterizados.
Dados Vetoriais - O formato de dados espaciais mais comum utilizado no GIS. Dados vetoriais compreendem as propriedades geométricas de vértices e trajetos (paths). Os dados vetoriais espaciais também podem se subdividir em três tipos amplamente utilizados:
Pontos - Um ponto consiste em um par ordenado (x,y) representando um local específico em um sistema de coordenadas. Pontos são a forma mais básica de dados espaciais e podem ser utilizados para representar um caso (casa de um paciente) ou um local (hospital) em um mapa.
Linhas - Uma linha é composta de dois pontos conectados. Linhas possuem um comprimento, e podem ser utilizadas para representar coisas como rodovias ou rios.
Polígonos - Um polígono é composto de pelo menos três segmentos de linha conectados por pontos. As propriedades de um Polígono são comprimento (o perímetro da área) bem como a própria medida da área. Polígonos podem ser utilizado para demarcar uma área (uma vila) ou uma estrutura (a área de um hospita ).
Dados Rasterizados - Um formato alternativo para dados espaciais, dados rasterizados são uma matriz de células (por exemplo, pixels) em que cada célula contém informações como altura, temperatura, inclinação, cobertura florestal, etc. Essas geralmente são fotos aéreas, imagens de satélite, etc. Os “Rasters” também podem ser utilizados como “mapas base” sob dados vetoriais.
Visualizando dados espaciais
Para representar visualmente os dados espaciais em um mapa, os softwares GIS exigem que você forneça informação suficiente sobre onde posicionar cada característica em relação umas às outras. Se você estiver utilizando dados vetoriais, o que vai acontecer na maioria dos casos, essas informações normalmente estarão armazenadas em um “shapefile”.
Shapefiles - Um shapefile é um formato de dado comum para armazenar dados espaciais “vetoriais” que consistem de linhas, pontos ou polígonos. Um único shapefile na verdade é um conjunto de pelo menos três arquivos - .shp, .shx, e .dbf. Todos esses sub-componentes devem estar presente em um mesmo diretório (pasta) para que o shapefile possa ser lido. Esses arquivos associados podem ser comprimidos em uma pasta ZIP para ser enviado por email ou baixado de um site.
O shapefile vai conter informação sobre as propriedades propriamente ditas, bem como onde encontrá-las na superfície da Terra. Isso é importante pois apesar da Terra ser um globo, os mapas são tipicamente bidimensionais; escolhas sobre como “achatar” os dados espaciais podem ter grande impactos no visual e na interpretação do mapa produzido.
Coordinate Reference Systems (CRS) (Sistemas de Referências de Coordenadas) - Um SRC é um sistema baseado em coordenadas utilizado para localizar propriedades geográficas na superfície da Terra. Ele possui alguns componentes principais:
Sistema de Coordenadas - Existem muitos tipos diferentes de sistemas de coordenadas, então certifique-se de qual dos tipos suas coordenadas fazem parte. Graus de latitude/longitude são comuns, mas você pode também ver coordenadas UTM.
Unidades - Saiba qual são as unidades para o seu sistema de coordenadas (exemplo: graus decimais, metros)
Datum - Uma versão modelada específica da Terra. Essa versão vem sendo revisada ao longo dos anos, então certifique-se de que as camadas do seu mapa utilizam o mesmo datum.
Projeção - Uma referência à equação matemática que foi utilizada para projetar a verdadeiramente redonda Terra em uma superfície plana (mapa).
Lembre-se que você pode resumir dados espaciais sem utilizar as ferramentas de mapas mostradas abaixo. Às vezes, uma simples tabela por dados geográficos (distritos, países, etc.) é todo o necessário!
28.3 Iniciando com GIS
Existem alguns itens importantes que você precisará ter ou pensar a respeito para fazer um mapa. Esses incluem:
- Uma base de dados – pode ser em formato de dados espaciais (tais como shapefiles, como mostrado acima) ou pode não ser em formato espacial (apenas um csv).
- Se sua base de dados não estiver em um formato espacial você também vai precisar de uma base de dados de referência. Dados de referência consistem nas representações espaciais de dados e seus atributos relacionados, que deverão incluir material contendo o local e informações sobre o endereço de propriedades específicas.
- Se você estiver trabalhando com fronteiras geográficas pré-definidas (por examplo, regiões administrativas), shapefiles de referências geralmente estão disponíveis gratuitamente para download nos sites das agências governamentais ou organizações de compartilhamento de dados. Se estiver em dúvida, um bom jeito de começar é buscar no Google por “[região] shapefile”
- Se você tiver informação sobre endereços, mas não tiver latitude e longitude, talvez você precisa utilizar um motor de geocodificação (geocoding engine) para coletar os dados de referências espaciais para seus regristros.
- Uma ideia sobre como você quer apresentar a informação em suas bases para seu público alvo. Existem muitos tipos diferentes de mapas, e é importante pensar sobre qual se encaixa melhor à suas necessidades.
Tipos de mapas para visualizar seus dados
Mapa coroplético (Choropleth map) - um tipo de mapa temático onde cores, sombreamento ou padrões são utilizados para representar regiões geográficas de acordo com seus valores de um atributo. Por exemplo, um valor maior pode ser indicado por uma cor mais escura que um valor menor. Esse tipo de mapa é particularmente útil quando estamos visualizando uma variável e como seu valor muda através de regiões definidas ou áreas geopolíticas.
Mapa de calor de densidade de casos (Case density heatmap) - um tipo de mapa temático em que as cores são utilizada para representar a intensidade de um valor, no entanto, ele não utiliza regiões definidas ou fronteiras geopolíticas para agrupar os dados. Esse tipo de mapa é tipicamente utilizado para mostrar ‘hot spots’ ou áreas com uma alta densidade ou concentração de pontos.
Mapa de densidade de pontos (Dot density map) - um tipo de mapa temático que usa pontos para reprensentar valores de atributos nos seus dados. Esse tipo de mapa é melhor utilizado para visualizar o espalhamento dos dados e permitir um “escaneamento” visual de agrupamentos (clusters).
Mapa de símbolos proporcionais (mapa de símbolos graduados) (Proportional symbols map (graduated symbols map)) - um mapa temático similar ao mapa de coroplético mas, em vez de utilizar cores para indicar o valor de um atributo, ele utiliza um símbolo (geralmente um círculo) de acordo com o valor. Por exemplo, um valor maior pode ser indicado por um símbolo maior que um valor menor. Esse tipo de mapa é melhor utilizado quando se quer visualizar o tamanho ou quantidade dos seus dados através de regiões geográficas.
Você também pode combinar vários tipos diferentes de visualização para mostrar padrões geográficos complexos. Por exemplo, os casos (pontos) no mapa abaixo estão coloridos de acordo com a unidade de saúde mais próxima (veja a legenda). Os círculos maiores mostram áreas de alcance das unidades de saúde a partir de um certo raio e os pontos vermelho-vivo mostram os casos fora do alcance de qualquer uma das unidades de saúde mostradas:
Nota: O foco principal dessa página sobre GIS é baseado no contexo de uma resposta de campo a um surto. Por isso, os conteúdos da página irão cobrir o básico de manipulação, visualização e análise de dados espaciais.
28.4 Preparação
Carregue os pacotes R
O código abaixo realiza o carregamento dos pacotes necessários para a análise dos dados. Neste manual, enfatizamos o uso da função p_load()
, do pacman, que instala os pacotes, caso não estejam instalados, e os carrega no R para utilização. Também é possível carregar pacotes instalados utilizando a função library()
, do R base. Para mais informações sobre os pacotes do R, veja a página Introdução ao R.
::p_load(
pacman# para importar os dados
rio, # para localizar arquivos
here, # para limpar, manipular e gerar gráficos com os dados (inclui o pacote ggplot2)
tidyverse, # para gerenciar dados espaciais utilizando o formato Simple Feature
sf, # para produzir mapas simples, funciona tanto para mapas interativos quanto estáticos
stars, # para limpar nomes das colunas
janitor, #dependencia tmap
terra,
tmap,# para adicionar o mapa-base do OpenStreeMap ao mapa do ggplot
OpenStreetMap,
tmaptools,
spdep
)
Você pode ter uma visão geral de todos os pacotes do R que lidam com dados espaciais em CRAN “Spatial Task View”.
Amostra de dados dos casos
Para fins de demonstração, iremos trabalhar com uma amostra aleatória de 1000 casos do dataframe da linelist
de uma epidemia simulada de Ebola (computacionalmente, trabalhar com menos casos é mais fácil para mostrar nesse handbook). Se você quiser acompanhar a análise, clique aqui para baixar a linelist “limpa” (como um arquivo .rds).
Uma vez que estamos tomando uma amostra aleatória dos casos, seus resultados podem parecer ligeiramente diferentes do que está demonstrado aqui quando você rodar os códigos por conta própria.
Importe os dados com a função import()
do pacote rio (ela lida com vários tipos de arquivo como .xlsx, .csv, .rds - veja a página Importar and exportar para mais detalhes).
# importa a linelist de casos limpa
<- rio::import(here::here("data", "case_linelists", "linelist_cleaned.rds")) linelist
Next we select a random sample of 1000 rows using sample()
from base R.
# gera 1000 números aleatórios de linhas, a partir do número de linhas da linelist
<- sample(nrow(linelist), 1000)
sample_rows
# gera um subconjunto da linelist para manter apenas as linhas da amostra, e todas as colunas
<- linelist[sample_rows,] linelist
Queremos converter essa linelist
, que possui a classe dataframe, para um objeto da classe “sf” (spatial features). Dado que a linelist possui duas colunas “lon” e “lat” representando a longitude e latitude da casa de cada caso, isso será fácil.
Utilizamos o pacote sf (spatial features) e sua função st_as_sf()
para criar um novo objeto que chamaremos de linelist_sf
. Esse novo objeto parece essencialmente a mesma coisa que a linelist, mas as colunas lon
e lat
foram designadas como colunas de coordenadas, e um sistema de referência de coordenadas (SRC) foi atribuído para quando os pontos forem mostrados. O número 4326 identifica nossas coordenadas de acordo com o World Geodetic System 1984 (WGS84) - que é o padrão para coordenadas de GPS.
# Cria um objeto do tipo sf
<- linelist %>%
linelist_sf ::st_as_sf(coords = c("lon", "lat"), crs = 4326) sf
O dataframe da linelist
original tem esta forma. Nesta demonstração, iremos utilizar apenas a coluna date_onset
e geometry
(que foi construída a partir dos campos de longitude e latitude acima e são as últimas colunas no dataframe).
::datatable(head(linelist_sf, 10), rownames = FALSE, options = list(pageLength = 5, scrollX=T), class = 'white-space: nowrap' ) DT
Shapefiles das fronteiras administrativas
Serra Leoa: Shapefiles das fronteiras administrativas
De antemão, já baixamos todas as fronteiras administrativas de Serra Leoa do site da Humanitarian Data Exchange (HDX) que você encontra aqui. Caso prefira, você pode baixar esse e todos os outros dados de exemplos para esse manual através de nosso pacote de R, conforme explicado na página Download dos dados no manual.
Agora faremos o seguinte para salvar no R o shapefile do nível administrativo 3
- Importar o shapefile
- Limpar o nome das colunas
- Filtar linhas para manter apenas as áreas de interesse
Para importar um shapefile usamos a função read_sf()
do pacote sf. À ela é fornecido o caminho do arquivo através da função here()
. - em nosso caso o arquivo está em nosso projeto R, nas subpastas “data”, “gis”, e “shp”, com o nome de arquivo “sle_adm3.shp” (veja as páginas Importar e exportar e Projetos R para mais informações). Você precisará fornecer seu próprio caminho de arquivo.
Em seguida usamos a função clean_names()
do pacote janitor para padronizar os nomes das colunas do shapefile. Também utilizamos filter()
para manter apenas as linhas com a área administrativa 2 (admin2name) chamada “Western Area Urban” ou “Western Area Rural”.
# nível ADM3 limpo
<- sle_adm3_raw %>%
sle_adm3 ::clean_names() %>% # padronizar o nome das colunas
janitorfilter(admin2name %in% c("Western Area Urban", "Western Area Rural")) # filtrar para manter certas áreas
Abaixo, podemos ver como fica o shapefile depois da importação e limpeza. Role para a direita para ver como há colunas com nível administrativo 0 (país), nível administrativo 1, nível administrativo 2, e finalmente nível administrativo 3. Cada nível tem um nome composto de letras e um identificador único chamado “pcode”. O pcode expande com cada nível administrativo que aumenta ex: SL (Serra Leoa) -> SL04 (Oeste) -> SL0410 (Area Rural do Oeste) -> SL040101 (Koya Rural).
Dados populacionais
Serra Leoa: População por ADM3
Novamente, esses dados podem ser baixadados do site da HDX (link aqui) ou pelo pacote epirhandbook explicado nessa página. Usamos import()
para carregar o arquivo .csv. Também passamos o arquivo importado para clean_names()
para padronizar a sintaxe dos nomes das colunas.
# População por nível ADM3
<- import(here("data", "gis", "population", "sle_admpop_adm3_2020.csv")) %>%
sle_adm3_pop ::clean_names() janitor
Essa é a aparência do arquivo de população. Role para a direita para ver como cada jurisdição tem colunas com a população masculina (male
), feminina (female
), e a população total
, e as subdivisões nas colunas columns por grupo de idade.
Unidades de Saúde
Serra Leoa: Dados das Unidades de Saúde do OpenStreetMap
Novamente, baixamos os locais das unidades de saúde do HDX aqui ou por meio das instruções na página manual de Download e dados.
Nós importamos o shapefile dos pontos das unidades com read_sf()
, e novamente limpamos os nomes das colunas, e depois, filtramos para manter apenas os pontos marcados com “hospital”, clínica (“clinic”), ou doutores (“doctors”).
# shapefile das unidades de saúde do OSM (open street map)
<- sf::read_sf(here("data", "gis", "shp", "sle_hf.shp")) %>%
sle_hf ::clean_names() %>%
janitorfilter(amenity %in% c("hospital", "clinic", "doctors"))
Aqui está o dataframe resultante - role para a direita para ver o nome da unidade e as coordenadas de geometria (geometry
).
28.5 Criando gráficos das coordenadas
A forma mais fácil de criar um gráfico de coordenadas X-Y (longitude/latitude, pontos), nessa situação dos casos, é desenhenhá-los como pontos diretamente do objeto linelist_sf
que criamos na seção de preparação.
O pacote tmap oferece algumas capacidades simples de fazer mapas tanto de forma estática ( modo “plot”) quanto interativa ( modo de visualização - “view”) com apenas algumas linhas de código. A sintaxe do tmap é parecida com a do ggplot2, de forma que os comandos são somados uns aos outros com o operador +
. Veja mais detalhes nessa vignette.
- Defina o modo do tmap. Neste caso iremos utilizar o modo “plot”, que produz saídas estáticas.
tmap_mode("plot") # escolha entre "view" ou "plot"
Abaixo, os pontos são projetados sozinhos. A função tm_shape()
é disponibilizada com os objetos linelist_sf
. Em seguida, adicionamos pontos via tm_dots()
, especificando o tamanho e a cor. Pelo fato de linelist_sf
ser um objeto do tipo sf, nós já temos as duas colunas designadas que contêm as coordenadas de latitude e longitude e o sistema de referência de coordenadas (CRS):
Sozinhos os pontos não nos dizem muita coisa. Então devemos também projetar as fronteiras administrativas:
Novamente, utilizaremos tm_shape()
(veja a documentação) mas em vez de fornecer o shapefile com os pontos dos casos, forneceremos o shapefile das fronteiras administrativas (polygons).
Com o argumento bbox =
(bbox significa “bounding box”) podemos especificar as coordenadas das fronteiras. Primeiro mostramos o mapa sem bbox
, e depois com ele.
# Apenas as fronteiras administrativas (polígonos)
tm_shape(sle_adm3) + # shapefile das fronteiras administrativas
tm_polygons(col = "#F7F7F7") + # exibir os polígonos em cinza claro
tm_borders(col = "#000000", # exibir as bordas com cor e peso de linha
lwd = 2) +
tm_text("admin3name") # coluna de texto para mostrar para cada polígono
# Mesmo que acima, mas com zoom à partir do bounding box
tm_shape(sle_adm3,
bbox = c(-13.3, 8.43, # vértice (do bounding box)
-13.2, 8.5)) + # vértice (do bounding box)
tm_polygons(col = "#F7F7F7") +
tm_borders(col = "#000000", lwd = 2) +
tm_text("admin3name")
E agora ambos, pontos e polígonos juntos:
# Todos juntos
tm_shape(sle_adm3, bbox = c(-13.3, 8.43, -13.2, 8.5)) + #
tm_polygons(col = "#F7F7F7") +
tm_borders(col = "#000000", lwd = 2) +
tm_text("admin3name")+
tm_shape(linelist_sf) +
tm_dots(size=0.08, col='blue', alpha = 0.5) +
tm_layout(title = "Distribution of Ebola cases") # Dá título ao gráfico
Para ler uma boa comparação das diferentes opções de criação de mapas no R, leia esse blog post.
28.6 “Joins” Espaciais
Você deve estar familiarizado com o conceito de fazer join dos dados de uma base em outra. Vários métodos são discutidos na página Juntando dados desse manual. Um “join espacial” tem um propósito similar, mas lidando com os relacionamentos espaciais dos dados. Em vez de utilizar valores comuns às colunas para parear corretamente as observações, você pode utilizar as relações espaciais deles, tais como o fato de uma característica estar contida em outra, ou de um ponto ser o vizinho mais próximo de outro, ou contido na área de influência (buffer) de um certo raio de outro, etc.
O pacote sf oferece vários métodos para fazer joins espaciais. Veja mais na documentação sobre o métodos st_join e outros tipos de join espaciais nessa referência.
Pontos no polígono
Relaciona unidades administrativas aos casos
Aqui temos um dilema interessante: a linelist dos casos não contém nenhuma informação sobre as unidades administrativas dos casos. Apesar do ideal ser coletar esse tipo de informação nas fases iniciais da coleta de dados, nós também podemos relacionar as unidades administrativas aos casos individuais baseado em seus relacionamentos espaciais (ex: pontos que intersectam com um polígono).
Abaixo, vamos fazer a intersecção espacial dos locais de nossos casos (pontos) com as fronteiras ADM3 (polígonos):
- Comece com a linelist (pontos)
- Faça o join espacial com as fronteiras, defina o tipo do join em “st_intersects”
- Utilize
select()
para manter apenas algumas das novas colunas das fronteiras administrativas
<- linelist_sf %>% # faça join do arquivo com as fronteiras administrativas com o da linelist, baseado na intersecção espacial
linelist_adm ::st_join(sle_adm3, join = st_intersects) sf
Todas as colunas de sle_adms
foram adicionadas à linelist! Agora cada caso possui colunas detalhando os níveis administativos onde eles se encontram. Nesse exemplo, queremos manter apenas duas das novas colunas (admin level 3), então usamos select()
nos nomes das colunas antigas e apenas nas duas novas que temos interesse:
<- linelist_sf %>% # faça join do arquivo com as fronteiras administrativas com o da linelist, baseado na intersecção espacial
linelist_adm ::st_join(sle_adm3, join = st_intersects) %>%
sf# Mantém os nomes das colunas antigas e as duas novas de interessa das áreas administrativas
select(names(linelist_sf), admin3name, admin3pcod)
Você pode ver abaixo, apenas para fins de conferência, os primeiros dez casos e as jurisdições de seus níveis administrativos 3 (ADM3) que foram anexadas baseados em onde o ponto intersectou espacialmente com as formas dos polígonos.
# Agora você verá os nomes AD3 anexados a cada caso.
%>% select(case_id, admin3name, admin3pcod) linelist_adm
Simple feature collection with 1000 features and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -13.27074 ymin: 8.447887 xmax: -13.20522 ymax: 8.491748
Geodetic CRS: WGS 84
First 10 features:
case_id admin3name admin3pcod geometry
3693 9c3676 East III SL040205 POINT (-13.2103 8.466647)
3521 5fbeea West III SL040208 POINT (-13.26599 8.473993)
4363 8e0af5 West II SL040207 POINT (-13.2459 8.468553)
3159 a509e6 Central I SL040201 POINT (-13.23163 8.478547)
3544 f5bc7b West III SL040208 POINT (-13.25092 8.45933)
1472 9fcf61 West III SL040208 POINT (-13.26314 8.465306)
341 8968ee West II SL040207 POINT (-13.24839 8.465598)
3802 bc6562 Mountain Rural SL040102 POINT (-13.21879 8.463824)
4230 a5b60c West II SL040207 POINT (-13.2487 8.470481)
128 6a007d West II SL040207 POINT (-13.23488 8.467028)
Agora podemos descrever nossos casos por unidade administrativa - algo que não era possível antes do join espacial!
# Cria um novo dataframe contendo as contagens de casos por unidade administrativa
<- linelist_adm %>% # inicia com a lineliste contendo as novas colunas administrativas
case_adm3 as_tibble() %>% # converte para tibble para melhor exibição
filter(!is.na(admin3pcod))%>%
group_by(admin3pcod, admin3name) %>% # agrupa por unidade administrativa, por nome e pcode
summarise(cases = n()) %>% # resume e conta as linhas
arrange(desc(cases)) # organiza em ordem descentente
case_adm3
# A tibble: 9 × 3
# Groups: admin3pcod [9]
admin3pcod admin3name cases
<chr> <chr> <int>
1 SL040102 Mountain Rural 279
2 SL040208 West III 233
3 SL040207 West II 174
4 SL040204 East II 120
5 SL040201 Central I 58
6 SL040203 East I 48
7 SL040206 West I 36
8 SL040205 East III 25
9 SL040202 Central II 24
Também podemos criar um gráfico de barras com as contagens dos casos por unidade administrativa.
Nesse exemplo, iniciamos o ggplot()
com linelist_adm
, de forma que podemos aplicar as funções de fatores como fct_infreq()
que ordena as barras por frequência (veja a página de Fatores para dicas).
ggplot(
data = linelist_adm, # inicia com linelist contendo informações das unidades administrativas
mapping = aes(
x = fct_rev(fct_infreq(admin3name))))+ # eixo x são unidades administrativas, ordenadas por frequência (reversa)
geom_bar()+ # cria as barras, altura é o número de linhas
coord_flip()+ # inverte os eixos X e Y para facilitar a leitura das unidades administrativas
theme_classic()+ # simplifica a cor de fundo do gráfico
labs( # títulos e rótulos
x = "Admin level 3",
y = "Number of cases",
title = "Number of cases, by adminstative unit",
caption = "As determined by a spatial join, from 1000 randomly sampled cases from linelist"
)
Vizinho mais próximo
Achando a unidade de saúde / área de influência mais próxima
Pode ser útil saber onde as unidades de saúde estão localizadas em relação aos hotspots das doenças.
Podemos utilizar o método de join st_nearest_feature da função st_join()
(do pacote sf) para visualizar as unidades de saude mais próximas aos casos individuais.
- Iniciamos com o shapefile da linelist:
linelist_sf
- Fazemos o join espacial com
sle_hf
, que possui os locais das unidades de saúde e clínicas (pontos)
# Unidade de saúde mais próxima a cada caso
<- linelist_sf %>% # inicia com o shapefile da linelist
linelist_sf_hf st_join(sle_hf, join = st_nearest_feature) %>% # dados da clínica mais próxima com join nos dados dos casos
select(case_id, osm_id, name, amenity) %>% # mantém as colunas de interesse, incluindo id, nome, tipo, e geometria da unidade de saúde
rename("nearest_clinic" = "name") # renomeia para uma maior clareza
Podemos ver abaixo (primeiras 50 linhas) que cada caso agora possui dados sobre as clínicas/hospitais mais próximos
Podemos ver que “Den Clinic” é a unidade de saúde mais próxima de cerca de 30% dos casos.
# Conta os casos por unidade de saúde
<- linelist_sf_hf %>% # inicia com a linelist incluindo os dados das clínicas mais próximas
hf_catchment as.data.frame() %>% # converte o shapefile em dataframe
count(nearest_clinic, # conta linhas por nome (da clínica)
name = "case_n") %>% # define a nova coluna de contagens como "case_n"
arrange(desc(case_n)) # organiza em ordem descrescente
# imprime no console hf_catchment
nearest_clinic case_n
1 Den Clinic 355
2 Shriners Hospitals for Children 326
3 GINER HALL COMMUNITY HOSPITAL 173
4 panasonic 59
5 ARAB EGYPT CLINIC 29
6 Princess Christian Maternity Hospital 29
7 MABELL HEALTH CENTER 17
8 <NA> 12
Para visualizar os resultados, podemos utilizar tmap - dessa vez, no modo interativo para facilitar a visualização
tmap_mode("view") # define o modo do tmap para interativo
# gera o gráfico dos pontos dos casos e clínicas
tm_shape(linelist_sf_hf) + # gera o gráfico dos casos
tm_dots(size=0.08, # colore os casos pelas clínicas mais próximas
col='nearest_clinic') +
tm_shape(sle_hf) + # adiciona as clínicas ao gráfico como grandes pontos pretos
tm_dots(size=0.3, col='black', alpha = 0.4) +
tm_text("name") + # sobrepõe com o nome
tm_view(set.view = c(-13.2284, 8.4699, 13), # ajusta o zoom (centra as coordenadas, e define o zoom)
set.zoom.limits = c(13,14))+
tm_layout(title = "Cases, colored by nearest clinic")
Áreas de Influência (Buffers)
Também podemos explorar quantos casos estão localizados em uma distância a pé de 2.5km (~30 min) da unidade de saúde mais próxima.
Nota: Para cálculos de distância mais precisos, é melhor re-projetar seu objeto sf ao sistema de projeção do local respectivo, tal como UTM (Terra projetada sobre uma superfície plana). Nesse exemplo, para fins de simplicidade vamos nos ater ao sistema de coordenadas do Sistema Geodésico Mundial (WGS84) (Terra representada em uma superfície esférica / redonda, portanto, as unidades estão em graus decimais). Utilizaremos uma conversão geral de: 1 grau decimal = ~111km.
Veja mais informações sobre projeções de mapas e sistemas de coordenadas nesse artigo da esri. Esse blog fala sobre diferentes tipos de projeções de mapas e como cada uma delas pode ser escolhida dependendo da área de interesse e o contexto do seu mapa / análise.
Primeiro, crie uma área de influência circular com um raio de ~2.5km ao redor de cada unidade de saúde. Isso é feito com a função st_buffer()
do pacote tmap. Pelo fato das unidades do mapa estarem em lat/long em graus decimais, é assim que “0.02” será interpretado. Se o sistema de coordenadas do seu mapa for em metros, o número deve ser fornecido em metros.
<- sle_hf %>%
sle_hf_2k st_buffer(dist=0.02) # graus decimais representando aprox. 2.5km
Abaixo criamos o gráfico das zonas de influência propriamente ditas, com:
tmap_mode("plot")
# Cria as zonas de influência circulares
tm_shape(sle_hf_2k) +
tm_borders(col = "black", lwd = 2)+
tm_shape(sle_hf) + # adiciona as clínicas ao gráfico como grandes pontos vermelhos
tm_dots(size=0.3, col='black')
Segundo, intersectamos essas áreas de influência com os casos (pontos) utilizando st_join()
e com o tipo de join st_intersects. Ou seja, join dos dados das áreas de influência com os pontos que eles intersectam.
# Intersecta os casos com as áreas de influência
<- linelist_sf_hf %>%
linelist_sf_hf_2k st_join(sle_hf_2k, join = st_intersects, left = TRUE) %>%
filter(osm_id.x==osm_id.y | is.na(osm_id.y)) %>%
select(case_id, osm_id.x, nearest_clinic, amenity.x, osm_id.y)
Agora podemos contar os resultados: nrow(linelist_sf_hf_2k[is.na(linelist_sf_hf_2k$osm_id.y), ])
dos 1000 casos não intersectam com nenhuma área de influência (esses valores estão faltando), e então residem a mais de 30 min a pé da da unidade de saúde mais próxima.
# Casos que não intesectaram com nenhuma das áreas de influência das unidades de saúde
%>%
linelist_sf_hf_2k filter(is.na(osm_id.y)) %>%
nrow()
[1] 1000
Podemos visualizar os resultados de forma que os casos que não intersectam apareçam em vermelho.
tmap_mode("view")
# Primeiro mostra os casos como pontos
tm_shape(linelist_sf_hf) +
tm_dots(size=0.08, col='nearest_clinic') +
# adiciona as clínicas como grandes pontos pretos
tm_shape(sle_hf) +
tm_dots(size=0.3, col='black')+
# Depois sobrepõe as áreas de influência das unidades de saúde como polilinhas
tm_shape(sle_hf_2k) +
tm_borders(col = "black", lwd = 2) +
# Destaca os casos que não são parte de nenhuma área de influência das unidades de saúde
# como pontos vermelhos
tm_shape(linelist_sf_hf_2k %>% filter(is.na(osm_id.y))) +
tm_dots(size=0.1, col='red') +
tm_view(set.view = c(-13.2284,8.4699, 13), set.zoom.limits = c(13,14))+
# adiciona título
tm_layout(title = "Cases by clinic catchment area")
Outros joins espaciais
Valores alternativos para o argumento join
incluem (a partir da documentação)
- st_contains_properly
- st_contains
- st_covered_by
- st_covers
- st_crosses
- st_disjoint
- st_equals_exact
- st_equals
- st_is_within_distance
- st_nearest_feature
- st_overlaps
- st_touches
- st_within
28.7 Mapas coropléticos
Mapas coropléticos podem ser úteis para visualizar seus dados por áreas pré-definidas, geralmente unidades administrativas ou de saúde. Nas respostas aos surtos isso pode ajudar a direcionar a alocação de recursos para áreas específicas com altas taxas de incidência, por exemplo.
Agora que já temos os nomes das unidades administrativas associadas a todos os casos (veja a seção sobre joins espaciais, acima), podemos começar a mapear as contagens dos casos por área (mapas coropléticos)
Uma vez que também temos dados de população por ADM3, podemos adicionar essas informações à tabela case_adm3 criada anteriormente.
Vamos começar com o dataframe case_adm3
criado no passo anterior, que é uma tabela resumo com cada unidade administrativa e seus números de casos.
- Os dados populacionais
sle_adm3_pop
são obtidos utilizandoleft_join()
do pacote dplyr com base nos valores comuns ao longo das colunasadmin3pcod
no dataframecase_adm3
, e colunaadm_pcode
no dataframesle_adm3_pop
. Veja a página Juntando dados).
select()
é aplicado ao novo dataframe, para manter apenas as colunas úteis -total
é a população total- Casos por 10,000 é calculado como uma nova coluna com
mutate()
# Adiciona dados populacionais e calcula casos por 10mil
<- case_adm3 %>%
case_adm3 left_join(sle_adm3_pop, # adiciona colunas da base 'pop'
by = c("admin3pcod" = "adm3_pcode")) %>% # faz join baseado nos valores comum ao longo dessas duas colunas
select(names(case_adm3), total) %>% # mantém apenas as colunas importantes, incluindo população total
mutate(case_10kpop = round(cases/total * 10000, 3)) # cria uma nova coluna com taxa de casos por 10000, arredondando para 3 casas decimais
# imprime no console para visualização case_adm3
# A tibble: 9 × 5
# Groups: admin3pcod [9]
admin3pcod admin3name cases total case_10kpop
<chr> <chr> <int> <int> <dbl>
1 SL040102 Mountain Rural 279 33993 82.1
2 SL040208 West III 233 210252 11.1
3 SL040207 West II 174 145109 12.0
4 SL040204 East II 120 99821 12.0
5 SL040201 Central I 58 69683 8.32
6 SL040203 East I 48 68284 7.03
7 SL040206 West I 36 60186 5.98
8 SL040205 East III 25 500134 0.5
9 SL040202 Central II 24 23874 10.1
Faz join dessa tabela com os polígonos do shapefile de ADM3 para poder mapeá-los
<- case_adm3 %>% # inicia com casos e taxas por unidade administrativa
case_adm3_sf left_join(sle_adm3, by="admin3pcod") %>% # faz join dos dados do shapefile por coluna comum
select(objectid, admin3pcod, # mantém apenas algumas colunas de interesse
admin3name = admin3name.x, # limpa o nome de uma coluna
admin2name, admin1name,
cases, total, case_10kpop,%>% # mantém a geometria de forma que os polígonos possam ser adicionados ao gráfico
geometry) st_as_sf() # converte para shapefile
Mapeando os resultados
# modo do tmap
tmap_mode("plot") # visualiza o mapa estático
# adiciona os polígonos
tm_shape(case_adm3_sf) +
tm_polygons("cases") + # colore pela coluna de número de casos
tm_text("admin3name") # nome exibido
Também podemos mapear as taxas de incidência
# Casos por 10mil
tmap_mode("plot") # modo de visualização estática
# gráfico
tm_shape(case_adm3_sf) + # adiciona polígonos ao mapa
tm_polygons("case_10kpop", # colore pelas colunas contendo a taxa de casos
breaks=c(0, 10, 50, 100), # define os pontos de quebra para as cores
palette = "Purples" # utiliza uma paleta de cor roxa
+
) tm_text("admin3name") # exibe o texto
28.8 Criando mapas com ggplot2
Se você já está familiarizado em usar o ggplot2, você pode usar esse pacote para criar mapas estáticos de seus dados. A função geom_sf()
vai desenhar diferentes objetos baseados em quais características (pontos, linhas ou polígonos) estão em seus dados. Por exemplo, você pode usar geom_sf()
em um ggplot()
utilizando dados sf
com geometria de polígonos para criar um mapa coroplético.
Para ilustrar como isso funciona, podemos iniciar com os polígonos do shapefile de ADM3 que utilizamos anteriormente. Lembre-se que esses são regiões de Níveis Administrativos 3 em Serra Leoa:
sle_adm3
Simple feature collection with 12 features and 19 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -13.29894 ymin: 8.094272 xmax: -12.91333 ymax: 8.499809
Geodetic CRS: WGS 84
# A tibble: 12 × 20
objectid admin3name admin3pcod admin3ref_n admin2name admin2pcod admin1name
* <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
1 155 Koya Rural SL040101 Koya Rural Western A… SL0401 Western
2 156 Mountain Ru… SL040102 Mountain R… Western A… SL0401 Western
3 157 Waterloo Ru… SL040103 Waterloo R… Western A… SL0401 Western
4 158 York Rural SL040104 York Rural Western A… SL0401 Western
5 159 Central I SL040201 Central I Western A… SL0402 Western
6 160 East I SL040203 East I Western A… SL0402 Western
7 161 East II SL040204 East II Western A… SL0402 Western
8 162 Central II SL040202 Central II Western A… SL0402 Western
9 163 West III SL040208 West III Western A… SL0402 Western
10 164 West I SL040206 West I Western A… SL0402 Western
11 165 West II SL040207 West II Western A… SL0402 Western
12 167 East III SL040205 East III Western A… SL0402 Western
# ℹ 13 more variables: admin1pcod <chr>, admin0name <chr>, admin0pcod <chr>,
# date <date>, valid_on <date>, valid_to <date>, shape_leng <dbl>,
# shape_area <dbl>, rowcacode0 <chr>, rowcacode1 <chr>, rowcacode2 <chr>,
# rowcacode3 <chr>, geometry <MULTIPOLYGON [°]>
Podemos utilizar a função left_join()
do pacote dplyr para adicionar dados que queremos mapear ao objeto shapefile. Nesse caso, vamos utilizar o dataframe case_adm3
que criamos anterioremente para resumir a contagem de casos por regiões administrativas; no entanto, podemos utilizar essa mesma abordagem para mapear qualquer dado armazenado no dataframe.
<- sle_adm3 %>%
sle_adm3_dat inner_join(case_adm3, by = "admin3pcod") # inner join = mantém apenas se existir em ambos os objetos
select(sle_adm3_dat, admin3name.x, cases) # imprime as variáveis selecionadas no console
Simple feature collection with 9 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -13.29894 ymin: 8.384533 xmax: -13.12612 ymax: 8.499809
Geodetic CRS: WGS 84
# A tibble: 9 × 3
admin3name.x cases geometry
<chr> <int> <MULTIPOLYGON [°]>
1 Mountain Rural 279 (((-13.21496 8.474341, -13.21479 8.474289, -13.21465 8.4…
2 Central I 58 (((-13.22646 8.489716, -13.22648 8.48955, -13.22644 8.48…
3 East I 48 (((-13.2129 8.494033, -13.21076 8.494026, -13.21013 8.49…
4 East II 120 (((-13.22653 8.491883, -13.22647 8.491853, -13.22642 8.4…
5 Central II 24 (((-13.23154 8.491768, -13.23141 8.491566, -13.23144 8.4…
6 West III 233 (((-13.28529 8.497354, -13.28456 8.496497, -13.28403 8.4…
7 West I 36 (((-13.24677 8.493453, -13.24669 8.493285, -13.2464 8.49…
8 West II 174 (((-13.25698 8.485518, -13.25685 8.485501, -13.25668 8.4…
9 East III 25 (((-13.20465 8.485758, -13.20461 8.485698, -13.20449 8.4…
Para fazer um gráfico de colunas com contagem de casos por região, utilizando ggplot2, podemos então utilizar geom_col()
dessa forma:
ggplot(data=sle_adm3_dat) +
geom_col(aes(x=fct_reorder(admin3name.x, cases, .desc=T), # reordena o eixo x 'cases' de forma decrescente
y=cases)) + # eixo y é o número de casos por região
theme_bw() +
labs( # define o texto da figura
title="Number of cases, by administrative unit",
x="Admin level 3",
y="Number of cases"
+
) guides(x=guide_axis(angle=45)) # angula os rótulos do eixo x para 45 graus para caber melhor
Se quisermos utilizar o ggplot2 para fazer um mapa coroplético da contagem dos casos, podemos utilizar uma sintaxe parecida para chamar a função geom_sf()
:
ggplot(data=sle_adm3_dat) +
geom_sf(aes(fill=cases)) # define a cor de preenchimento (fill) para variar de acordo com a variável de contagem dos casos
Podemos customizar a aparência de nosso mapa utilizando a gramática que é consistente no ggplot2, por exemplo:
ggplot(data=sle_adm3_dat) +
geom_sf(aes(fill=cases)) +
scale_fill_continuous(high="#54278f", low="#f2f0f7") + # muda a cor do gradiente
theme_bw() +
labs(title = "Number of cases, by administrative unit", # define o texto da figura
subtitle = "Admin level 3"
)
Para usuários de R que estão confortáveis em trabalhar com ggplot2, a função geom_sf()
oferece uma implementação simples e direta que é adequada para visualizações básicas de mapas. Para aprender mais, leia a vinheta do geom_sf() ou o livro do ggplot2.
28.9 Mapas Base (Basemaps)
OpenStreetMap
Descrevemos abaixo como gerar um mapa base (basemap) para um mapa do ggplot2 utilizando funcionalidades do OpenStreetMap. Métodos alternativos incluem utilizar o ggmap que requer um cadastro gratuito via Google (detalhes).
OpenStreetMap é um projeto colaborativo para criar um mapa do mundo editável e gratuito. Os dados de geolocalização fundamentais (ex: localidades de cidades, estradas, características naturais, aeroportos, escolas, hospitais, etc) são considerados as saídas (outputs) principais do projeto.
Primeiro carregamos o pacote OpenStreetMap, de onde pegaremos nosso mapa base.
Depois, criamos o objeto map
, o qual definimos utilizando a função openmap()
do pacote OpenStreetMap (documentação). Passaremos os seguintes argumentos para função:
upperLeft
elowerRight
dois pares de coordenadas especificando os limites do tile do mapa base- Nesse caso inserimos o max e min das linhas da linelist, para que o mapa responda dinamicamente aos dados
zoom =
(se não for preenchido será determinado automaticamente)
type =
qual tipo de mapa base - listamos várias possibilidades aqui e o código está utilizando a primeira delas ([1]
) “osm”
mergeTiles =
optamos por TRUE para que os tiles base sejam mesclados em um só
# carrega o pacote
::p_load(OpenStreetMap)
pacman
# Ajusta o mapa base pelo intervalo das coordenadas de lat/long. Define o tipo do tile
<- openmap(
map upperLeft = c(max(linelist$lat, na.rm=T), max(linelist$lon, na.rm=T)), # limites do tile do mapa base
lowerRight = c(min(linelist$lat, na.rm=T), min(linelist$lon, na.rm=T)),
zoom = NULL,
type = c("osm", "stamen-toner", "stamen-terrain", "stamen-watercolor", "esri","esri-topo")[1])
Se executarmos esse mapa base nesse momento, utilizando autoplot.OpenStreetMap()
do pacote OpenStreetMap, você verá que as unidades nos eixos não são coordenadas de latitude/longitude. Ele está utilizando um sistema de coordenadas diferente. Para mostrar corretamente as residências dos casos (que estão armazenados em lat/long), isso precisa ser modificado.
Assim, nossa intenção é converter o mapa para latitude/longitude com a função openproj()
do pacote OpenStreetMap. Nós passamos para a função o mapa base map
e também o Sistema de Referência de Coordenadas (SRC) que queremos. Fazemos isso passando a string do tipo “proj.4” para a projeção WGS 1984, mas você pode passar o SRC de outras formas também. (veja esta página para compreender melhor do que se trata uma string do tipo proj.4)
# Projeção WGS84
<- openproj(map, projection = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") map_latlon
Agora, quando executamos o mapa, vemos que ao longo dos eixos estão as coordenadas de latitude e longitude. O sistema de coordenadas foi convertido. Agora todos os nossos casos serão mostrados corretamente se sobrepostos ao mapa!
# Gera o mapa. É necessário usar "autoplot" para funcionar com o ggplot
autoplot.OpenStreetMap(map_latlon)
28.10 Mapas de calor de densidade com contornos
Abaixo descrevemos como gerar um mapa de calor de densidade com contornos para os dados dos casos, sobre um mapa base, a partir de uma linelist (uma linha por caso).
- Cria o tile do mapa base a partir do OpenStreetMap, como descrito acima.
- Adiciona os casos da
linelist
utilizando as colunas de latitude e longitude - Converte os pontos para um mapa de calor de densidade com
stat_density_2d()
do ggplot2,
Quando temos um mapa base com coordenadas lat/long, podemos sobrepor nossos casos utilizando as coordenadas lat/long de suas residências.
Utilizando a função autoplot.OpenStreetMap()
para criar o mapa base, as funções do ggplot2 podem facilmente sobrepor camadas sobre ele, como mostrado com geom_point()
abaixo:
# Gera o mapa. É necessário usar "autoplot" para funcionar com o ggplot
autoplot.OpenStreetMap(map_latlon)+ # inicia com o mapa base
geom_point( # adiciona os pontos xy a partir das colunas lo e lat da linelist
data = linelist,
aes(x = lon, y = lat),
size = 1,
alpha = 0.5,
show.legend = FALSE) + # remove a legenda completamente
labs(x = "Longitude", # títulos e rótulos
y = "Latitude",
title = "Cumulative cases")
O mapa acima pode ser difícil de interpretar, principalmente pelos pontos se soprepondo. Então, você pode, ao invés, gerar um mapa de densidade 2d utilizando a função stat_density_2d()
do ggplot2. Você ainda estará utilizando as coordenadas lat/lon da linelist, mas uma estimativa de densidade por kernel 2D será calculada e os resultados serão mostrados como linhas de contorno - como em um mapa topográfico. Leia a documentação aqui.
# inicia com um mapa base
autoplot.OpenStreetMap(map_latlon)+
# adiciona o mapa de densidade
::stat_density_2d(
ggplot2data = linelist,
aes(
x = lon,
y = lat,
fill = ..level..,
alpha = ..level..),
bins = 10,
geom = "polygon",
contour_var = "count",
show.legend = F) +
# especifica a escala de cores
scale_fill_gradient(low = "black", high = "red")+
# rótulos
labs(x = "Longitude",
y = "Latitude",
title = "Distribution of cumulative cases")
Mapa de calor de séries temporais
O mapa de calor de densidade acima mostra os casos cumulativos. Podemos investigar o surto ao longo do tempo e espaço dividindo o mapa de calor em facetas baseadas no mês de início dos sintomas (symptom onset), derivado da linelist.
Começamos com a linelist
, criando uma nova coluna com o Ano (Year) e Mês (Month) do início. A função format()
do R base modifica a forma como a data é exibida. Nesse caso, queremos o formato “YYYY-MM”.
# Extrai o mês de início
<- linelist %>%
linelist mutate(date_onset_ym = format(date_onset, "%Y-%m"))
# Investiga os valores
table(linelist$date_onset_ym, useNA = "always")
2014-05 2014-06 2014-07 2014-08 2014-09 2014-10 2014-11 2014-12 2015-01 2015-02
10 19 42 97 180 181 127 104 64 50
2015-03 2015-04 <NA>
45 35 46
Agora, nós simplesmente adicionamos as facetas via ggplot2 ao mapa de calor de densidade. Aplicamos a função facet_wrap()
, utilizando a nova coluna como linhas. Nós definimos o número de colunas das facetas para 4, para fins de clareza.
# pacotes
::p_load(OpenStreetMap, tidyverse)
pacman
# inicia com o mapa base
autoplot.OpenStreetMap(map_latlon)+
# adiciona o mapa de densidade
::stat_density_2d(
ggplot2data = linelist,
aes(
x = lon,
y = lat,
fill = ..level..,
alpha = ..level..),
bins = 10,
geom = "polygon",
contour_var = "count",
show.legend = F) +
# especifica a escala de cores
scale_fill_gradient(low = "black", high = "red")+
# rótulos
labs(x = "Longitude",
y = "Latitude",
title = "Distribution of cumulative cases over time")+
# gera as facetas do gráfico por mês-ano de início
facet_wrap(~ date_onset_ym, ncol = 4)
28.11 Estatística espacial
Grande parte de nossa discussão até aqui estava focada na visualização de dados espaciais. Em alguns casos, você também pode estar interessado em estatística espacial para quantificar as relações espaciais entre os atributotos de seu dados. Essa seção vai passar uma visão global sobre alguns conceitos chave da estatística espacial, e sugerir alguns recursos que podem explorados caso você deseje fazer uma análise espacial mais aprofundada.
Relações espaciais
Antes de calcularmos qualquer estatística espacial, precisamos especificar as relações entre as características de nossos dados. Existem várias formas de conceituar relações espaciais, mas uma forma simples e um modelo comumente aplicado é utilizar o conceito de adjacência - especificamente, o fato de esperarmos uma relação geográfica entre áreas que compartilham uma fronteira ou são vizinhos uns dos outros.
Podemos quantificar relações de adjacências entre polígonos de regiões administrativas da base que viemos utilizando (sle_adm3
) com o pacote spdep. Iremos especificar contiguidade do tipo queen (rainha), o que significa que as regiões serão vizinhas se compartilharem pelo menos um ponto ao longo de suas fronteiras. O método alternativo seria contiguidade do tipo rook (torre), que requer que as regiões compartilhem uma aresta - em nosso caso, com polígonos irregulares, a diferença é trivial, mas em alguns casos a escolha entre queen e rook pode fazer diferença. (Nota do tradutor: Os termos queen (rainha) e rook (torre) correspondem à peças do xadrez e o tipo de vizinhança que representam descritas acima estão relacionadas aos movimentos que fazem no referido jogo: na diagonal para rainha, ou seja apenas com um ponto de contato com a casa vizinha, e para frente ou para os lados com a torre, ou seja com uma aresta -a lateral do quadrado- completa de contato com a casa vizinha).
<- spdep::poly2nb(sle_adm3_dat, queen=T) # cria vizinhos
sle_nb <- spdep::nb2mat(sle_nb) # cria uma matriz resumindo os relacionamentos entre vizinhos
sle_adjmat <- spdep::nb2listw(sle_nb) # cria um um objeto listw (lista de pesos (weights) ) -- iremos precisar disso posteriormente
sle_listw
sle_nb
Neighbour list object:
Number of regions: 9
Number of nonzero links: 30
Percentage nonzero weights: 37.03704
Average number of links: 3.333333
round(sle_adjmat, digits = 2)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
1 0.00 0.20 0.00 0.20 0.00 0.2 0.00 0.20 0.20
2 0.25 0.00 0.00 0.25 0.25 0.0 0.00 0.25 0.00
3 0.00 0.00 0.00 0.50 0.00 0.0 0.00 0.00 0.50
4 0.25 0.25 0.25 0.00 0.00 0.0 0.00 0.00 0.25
5 0.00 0.33 0.00 0.00 0.00 0.0 0.33 0.33 0.00
6 0.50 0.00 0.00 0.00 0.00 0.0 0.00 0.50 0.00
7 0.00 0.00 0.00 0.00 0.50 0.0 0.00 0.50 0.00
8 0.20 0.20 0.00 0.00 0.20 0.2 0.20 0.00 0.00
9 0.33 0.00 0.33 0.33 0.00 0.0 0.00 0.00 0.00
attr(,"call")
spdep::nb2mat(neighbours = sle_nb)
A matriz exibida acima mostra as relações entre as 9 regiões de nossa base sle_adm3
. Um score de 0 indica duas regiões que não são vizinhas, enquanto quaisquer outros valores diferentes de zero indicam uma relação de vizinhança. Os valores na matriz estão escalados de forma que a linha de cada região tenha um peso total de 1.
Uma melhor forma de visualizar essas relações de vizinhança é gerando um gráfico para elas:
plot(sle_adm3_dat$geometry) + # cria gráficos com as fronteiras das regiões
::plot.nb(sle_nb,as(sle_adm3_dat, 'Spatial'), col='grey', add=T) # adiciona relações de vizinhança spdep
Utilizamos uma abordagem de adjacência para identificar polígonos vizinhos; os vizinhos identificados também podem ser chamados de vizinhos baseados em contiguidade (contiguity-based neighbors). Mas, essa é apenas um das formas de escolher quais regiões são esperadas de possuírem uma relação geográfica. As abordagens alternativas mais comuns para identificar relações geográficas geram vizinhos baseados em distância (distance-based neighbors); basicamente, eles são:
- K-vizinhos mais próximos (K-nearest neighbors) - Baseados nas distâncias entre centróides (centros geograficamente-pesados dos polígonos de cada região), seleciona as n regiões mais próximas como vizinhas. Um limite de proximidade de máxima-distância também pode ser especificado. Em spdep, você pode usar
knearneigh()
(veja documentação). - Vizinhos por limite de distância - Seleciona todos os vizinhos contidos em um limite de distância. Em spdep, essa relação de vizinhos pode ser identificada utilizando a função
dnearneigh()
(veja documentação).
Autocorrelação espacial
A conhecida primeira lei da geografia de Tobler diz que “todas as coisas estão relacionadas com todas as outras, mas coisas próximas estão mais relacionadas do que coisas distantes.” Na epidemiologia, isso geralmente significa que o risco de um certo desfecho de saúde em uma dada região é mais similar ao de suas regiões vizinhas que daquelas regiões mais longínquas. Esse conceito foi formalizado como autocorrelação espacial - a propridedade estatística em que características geográficas com valores similares são agregados no espaço. Medidas estatísticas de autocorrelação espacial podem ser utilizadas para quantificar a extensão da agregação espacial em seus dados, localizar onde a agregação ocorre, e identificar padrões compartilhados da autocorrelação espacial entre diferentes variáveis dos seus dados. Essa seção te dá uma visão global de algumas medidas de autocorrelação espacial mais comuns e como calculá-las no R.
I de Moran - Esta é uma estatística resumo global da correlação entre o valor de uma varíavel em uma região e os valores da mesma variável em regiões vizinhas. O I de Moran geralmente varia entre -1 e 1. O valor 0 indica que não há nenhum padrão de correlação espacial, enquanto valores próximos a 1 ou -1 indicam uma autocorrelação espacial mais forte (valores similares próximos de si) ou dispersão espacial (valores dissimilares próximoas de si), respectivamente.
Por exemplo, podemos calcular o I de Moran para quantificar a autocorrelação espacial nos casos de Ebola que mapeamos mais cedo (lembre-se, esse é um subset dos casos do dataframe da linelist
da epidemia simulada). O pacote spdep tem uma função moran.test
, que pode fazer esse cálculo para nós:
<-spdep::moran.test(sle_adm3_dat$cases, # vetor numérico com a variável de interesse
moran_i listw=sle_listw) # objeto listw resumindo as relações entre vizinhos
# imprime os resultados do teste do I de Moran moran_i
Moran I test under randomisation
data: sle_adm3_dat$cases
weights: sle_listw
Moran I statistic standard deviate = 1.6653, p-value = 0.04793
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.22479823 -0.12500000 0.04412088
A saída da função moran.test()
nos mostra um valor do I de Moran de round(moran_i$estimate[1],2)
. Isso indica a presença de autocorrelação espacial em nossos dados - especificamente, que regiões com números parecidos de casos de Ebola estão potencialmente próximos entre si. O p-valor fornecido por moran.test()
é gerado por comparação à expectativa, sob a hipótese nula, de não haver autocorrelação espacial, e pode ser utilizado caso você precise mostrar os resultados de um teste de hipóteses formal.
I de Moran Local - Podemos decompor o I de Moran (global) calculado acima para identificar uma autocorrelação espacial localizada; ou seja, para identificar agrupamentos específicos em nossos dados. Essa estatística, que as vezes é chamada de Indicador Local de Associação Espacial Local Indicator of Spatial Association (LISA), resume a extensão da autocorrelação espacial sobre cada região individual. Ela pode ser útil para achar regiões “quentes” ou “frias” no mapa.
Para mostrar um exemplo, podemos calcular e mapear I’s de Moran locais para as contagens de casos de Ebola utilizados acima, com a função local_moran()
do pacote spdep:
# calcula o I de Moran local
<- spdep::localmoran(
local_moran $cases, # variável de interesse
sle_adm3_datlistw=sle_listw # objeto listw com pesos para os vizinhos
)
# faz join dos resultados aos dados do sf (shapefile)
<- cbind(sle_adm3_dat, local_moran)
sle_adm3_dat
# gera o mapa
ggplot(data=sle_adm3_dat) +
geom_sf(aes(fill=Ii)) +
theme_bw() +
scale_fill_gradient2(low="#2c7bb6", mid="#ffffbf", high="#d7191c",
name="Local Moran's I") +
labs(title="Local Moran's I statistic for Ebola cases",
subtitle="Admin level 3 regions, Sierra Leone")
Gi (ou G) de Getis-Ord - Essa é outra estatística comumente utilizada para análises de hotspots; em grande parte, a popularidade dessa estatística está relacionada à sua utilização na ferramenta de análise de hotspots no ArcGIS. Isso é baseado na premissa de que tipicamente, a diferença entre os valores de uma variável entre regiões vizinhas deve seguir uma distribuição normal. Ela usa a abordagem do z-score para identificar regiões que possuem valores de uma dada variável que sejam significativamente mais altas (hot spot) ou significativamente mais baixas (cold spot), comparados aos seus vizinhos.
Podemos calcular e mapear a estatística Gi* utilizando a função localG()
do spdep:
# Faz uma análise G
<- spdep::localG(
getis_ord $cases,
sle_adm3_dat
sle_listw
)
# faz join dos resultados aos dados do sf
$getis_ord <- as.numeric(getis_ord)
sle_adm3_dat
# gera o mapa
ggplot(data=sle_adm3_dat) +
geom_sf(aes(fill=getis_ord)) +
theme_bw() +
scale_fill_gradient2(low="#2c7bb6", mid="#ffffbf", high="#d7191c",
name="Gi*") +
labs(title="Getis-Ord Gi* statistic for Ebola cases",
subtitle="Admin level 3 regions, Sierra Leone")
Como você pode ver, o mapa de Gi* de Getis-Ord tem um visual ligeiramente diferente do mapa de I de Moran Local produzido anteriormente. Isso demonstra que o método utilizado para calcular essas duas estatísticas são ligeiramente diferentes; aquele que você vai utilizar depende do seu caso de uso específico e sua pergunta experimental de interesse.
Test L de Lee - Esse é um teste estatístico para correlação espacial bivariada. Ele permite que você teste se um padrão espacial para uma dada variável x é similar ao padrão espacial de outra variável y, cuja hipótese é que seja espacialmente relacionada com x.
Para dar um exemplo, vamos testar se um padrão espacial dos casos de Ebola da epidemia simulada é correlacionado com o padrão espacial da população. Para iniciar, precisamos de uma variável population
em nossa base sle_adm3
. Podemos utilizar a variável total
do dataframe sle_adm3_pop
que carregamos anteriormente.
<- sle_adm3_dat %>%
sle_adm3_dat rename(population = total) # renomeia 'total' para 'population'
Podemos visualizar rapidamente os padrões espaciais das duas variáveis lado a lado, para ver se eles se parecem:
tmap_mode("plot")
<- tm_shape(sle_adm3_dat) + tm_polygons("cases") + tm_layout(main.title="Cases")
cases_map <- tm_shape(sle_adm3_dat) + tm_polygons("population") + tm_layout(main.title="Population")
pop_map
tmap_arrange(cases_map, pop_map, ncol=2) # organiza em facetas 2x1
Visualmente, os padrões parecem dissimilares. Podemos utilizar a função lee.test()
do pacote spdep para testar estatisticamente se o padrão de autocorrelação espacial nas duas variáveis é relacionado. A estística L será próxima de 0 se não houver correlação entre os padrões, e próxima de 1 se tiver uma forte correlação positiva (ou seja, os padrões são similares), e próximo de -1 se houver uma forma correlação negativa (ou seja, os padrões são inversos).
<- spdep::lee.test(
lee_test x=sle_adm3_dat$cases, # variável 1 para comparar
y=sle_adm3_dat$population, # variável 2 para comparar
listw=sle_listw # objeto listw com peso dos vizinhos
)
lee_test
Lee's L statistic randomisation
data: sle_adm3_dat$cases , sle_adm3_dat$population
weights: sle_listw
Lee's L statistic standard deviate = -0.86785, p-value = 0.8073
alternative hypothesis: greater
sample estimates:
Lee's L statistic Expectation Variance
-0.13219674 -0.03580398 0.01233658
A saída acima mostra que a estatística L de Lee para as nossas duas variáveis é round(lee_test$estimate[1],2)
, o que indica uma fraca correlação negativa. Isso confirma nossa avaliação visual de que o padrão dos casos e população não estão relacionados entre si, e traz evidência de que o padrão espacial dos casos não é estritamente um resultado da densidade populacional em áreas de alto risco.
A estatística L de Lee pode ser útil para fazer esses tipos de inferência sobre a relação entre variáveis espacialmente distribuídas; no entanto, para descrever a natureza do relacionamento entre duas variáveis com mais detalhes, ou ajustar para fatores de confusão, técnicas de regressão espacial serão necessárias. Essas são descritas brevemente na seção seguinte.
Regressão espacial
Talvez você queira fazer inferências estatísticas a respeito dos relacionamentos entre variáveis nos seus dados espaciais. Nesses casos, pode ser útil considerar técnicas de regressão espacial - ou seja, abordagens para regressão que consideram explicitamente a organização espacial das unidades dos seus dados. Algumas razões que podem fazer você considerar modelos de regressão espacial, em vez de modelos de regressão padrão tais como GLMs, incluem:
- Modelos padrão de regressão assumem que os resíduos são independentes uns dos outros. Na presença de fortes autocorrelações espaciais, os resíduos de um modelo padrão de regressão têm a potencialidade de estarem espacialmente autocorrelacionados também, portanto, violando essa premissa. Isso pode levar a problemas com a interpretação dos resultados do modelo, o que para o caso, um modelo espacial seria mais adequado.
- Modelos de regressão também assumem, tipicamente, que o efeito de uma variável x é constante sobre todas as observações. No caso de heterogeneidade espacial, os efeitos que queremos estimar podem variar sobre o espaço, e podemos estar interessados em quantificar essas diferenças. Nesse caso, modelos de regressão espacial oferecem mais flexibilidade para estimar e interpretar os efeitos.
Os detalhes das abordagens de regressão espacial estão além do escopo desse manual. Essa seção vai, ao invés, passar uma visão global dos modelos de regressão espacial mais comuns e seus usos, e sugerir referências que podem ser úteis caso você queira explorar essa área mais a fundo.
Modelos de Erro Espacial - Esses modelos assumem que os termos de erro ao longo de unidades espaciais são correlacionados, sendo esse o caso, os dados violariam as premissas de um modelo OLS padrão. Modelos de Erro Espacial também podem ser chamados de modelos autorregressivos simultâneos (simultaneous autoregressive (SAR) models). Eles podem ser ajustados utilizando a função errorsarlm()
do pacote spatialreg (funções de regressão espacial que costumavam fazer parte do pacote spdep).
Modelos de Defasagem Espacial (spacial lag models) - Esses modelos assumem que a variável dependente para uma região i é influenciada não apenas pelo valor das variáveis independentes em i, mas também pelos valores daquelas variáveis em regiões vizinhas a i. Assim como os modelos de erro espacial, modelos de defasagem espacial também são descritos algumas vezes como modelos autorregressivos simultâneos (simultaneous autoregressive (SAR) models). Eles podem ser ajustados utilizando a função lagsarlm()
do pacote spatialreg.
O pacote spdep contém vários testes diagnóstico úteis para decidir entre os modelos OLS padrão, de defasagem espacial ou de erro espacial. Esses testes, chamados diagnósticos de Multiplicador de (Lagrange Multiplier diagnostics), podem ser utilizados para identificar o tipo de dependência espacial em seus dados e escolher qual modelo é mais apropriado. A função lm.LMtests()
pode ser utilizada para calcular todos os testes de Multiplicador de Lagrange. Anselin (1988) também demonstra uma ferramenta de fluxogramas para decidir qual modelo de regressão espacial a ser utilizado baseado nos resultados dos testes de Multiplicador de Lagrange:
Modelos Hierárquicos Bayesianos - Abordagens Bayesianas são comumente utilizadas para algumas aplicação em análise espacial, principalmente para mapeamento de doenças. Elas são preferidas em ocasiões onde os dados dos casos estão esparçamente distribuídas (por exemplo, no caso de um desfecho raro) ou com muito “ruído” estatístico, de forma que elas podem ser utilizadas para gerar estimativas “suavizadas” do risco de doenças ao considerar processos espaciais subjacentes que estavam latentes. Isso pode aumentar a qualidade das estimativas. Elas também permitem ao investigador a pre-especificação (via escolha de uma priori) de complexos padrões de correlação espaciais que podem existir nos dados, o que pode considerar variações espaço-dependente e -independente tanto nas variáveis independentes quanto nas dependentes. No R, Modelos Hierarquicos Bayesianos podem ser ajustados utilizando o pacote CARbayes (veja documentação) ou R-INLA (veja website e livro). R também pode ser utilizado para chamar softwares externos que fazem estimativas Bayesianas, tais como JAGS ou WinBUGS.
28.12 Recursos
- Pacote Simple Features vignette
- Pacote tmap vignette
- ggmap: Visualização Espacial com ggplot2
- Introdução à produção de mapas com R, visão geral de diferentes pacotes
- Dados Espaciais no R (EarthLab course)
- Análise de Dados Espaciais Aplicada com R livro
- SpatialEpiApp - um Shiny app que pode ser baixado como um pacote R, permitindo que você insira seus próprios dados e faça a criação de mapas, análises de agrupamento (clusters) e estatísticas espaciais.
- Introdução à Economoetria Espacial com R workshop