Phần mềm GIS

Một số phần mềm GIS phổ biến cho phép tương tác bằng các thao tác chuột để phát triển bản đồ và phân tích không gian. Những công cụ này đi kèm với những lợi thế như không cần phải học code và dễ dàng lựa chọn và đặt các biểu tượng và tính năng trên bản đồ theo cách thủ công. Dưới đây là hai phần mềm phổ biến:

ArcGIS - Một phần mềm GIS thương mại do công ty ESRI phát triển, rất phổ biến nhưng khá đắt

QGIS - Một phần mềm GIS mã nguồn mở nhưng làm được hầu hết mọi thứ mà ArcGIS có thể làm được. Bạn có thể tải QGIS tại đây

Sử dụng R để thao tác GIS thoạt đầu có vẻ đáng sợ hơn bởi vì thay bằng các “thao tác chuột”, nó có “giao diện dòng lệnh” (bạn phải viết code để có được kết quả mong muốn). Tuy nhiên, đây là một lợi thế lớn nếu bạn cần tạo bản đồ lặp đi lặp lại hoặc tạo các phân tích có thể tái lập được.

Dữ liệu không gian

Hai dạng dữ liệu không gian chính được sử dụng trong GIS là dữ liệu vectơ và raster:

Vector Data - Định dạng phổ biến nhất của dữ liệu không gian được sử dụng trong GIS, dữ liệu vectơ bao gồm các đặc điểm hình học của vertices and paths. Dữ liệu không gian dạng vectơ có thể được chia thành ba loại nhỏ được sử dụng rộng rãi:

• Điểm (Points) - Một điểm bao gồm một cặp tọa độ (x, y) đại diện cho một vị trí cụ thể trong một hệ tọa độ. Điểm là dạng dữ liệu không gian cơ bản nhất và có thể được sử dụng để biểu thị một trường hợp (vd: nhà bệnh nhân) hoặc một vị trí (vd: bệnh viện) trên bản đồ.

• Đường (Lines) - Một đường bao gồm hai điểm kết nối với nhau. Đường có độ dài và có thể được sử dụng để biểu thị những thứ như con đường hoặc sông.

• Đa giác (Polygons) - A polygon is composed of at least three line segments connected by points. Polygon features have a length (i.e. the perimeter of the area) as well as an area measurement. Polygons may be used to note an area (i.e. a village) or a structure (i.e. the actual area of a hospital).

Raster Data - Một đa giác bao gồm ít nhất ba đoạn thẳng được nối với nhau bằng các điểm. Các đối tượng đa giác có chiều dài (vd: chu vi của khu vực) cũng như số đo diện tích. Đa giác có thể được sử dụng để biểu diễn một khu vực (vd: một ngôi làng) hoặc một cấu trúc (vd: diện tích thực tế của một bệnh viện).

Trực quan hóa dữ liệu không gian

Để thể hiện trực quan dữ liệu không gian trên bản đồ, phần mềm GIS yêu cầu bạn cung cấp đầy đủ thông tin về vị trí của các đối tượng địa lý khác nhau, trong mối quan hệ của chúng với nhau. Nếu bạn đang sử dụng dữ liệu vectơ, điều này sẽ đúng cho hầu hết các trường hợp sử dụng, thông tin này thường sẽ được lưu trữ trong một shapefile:

Shapefiles - Shapefile là một định dạng dữ liệu phổ biến để lưu trữ dữ liệu không gian “vectơ” bao gồm hoặc đường, điểm hoặc đa giác. Một shapefile thực chất là một tập hợp của ít nhất ba tệp - .shp, .shx và .dbf. Tất cả các tệp thành phần phụ này phải nằm trong cùng một thư mục để chúng có thể đọc được. Các tệp liên quan này có thể được nén vào một thư mục ZIP để gửi qua email hoặc tải xuống từ một trang web.

Shapefile sẽ chứa thông tin về bản thân các đối tượng địa lý cũng như vị trí định vị chúng trên bề mặt Trái đất. Điều này rất quan trọng bởi vì trong khi Trái đất là một quả địa cầu, các bản đồ thường là hai chiều; các lựa chọn về cách “làm phẳng” dữ liệu không gian có thể có tác động lớn đến giao diện và cách giải thích các kết quả..

Hệ trục tọa độ tham chiếu (Coordinate Reference Systems - CRS) - CRS là một hệ thống dựa trên tọa độ được sử dụng để xác định vị trí các đối tượng địa lý trên bề mặt Trái đất. Nó có một số thành phần chính:

• Hệ tọa độ - Có nhiều hệ tọa độ khác nhau, vì vậy hãy đảm bảo rằng bạn biết hệ tọa độ của mình là từ hệ nào. Các kinh độ/vĩ độ là phổ biến nhẩt, nhưng bạn cũng có thể gặp hệ tọa độ UTM.

• Đơn vị - Các đơn vị dành cho hệ tọa độ của bạn (ví dụ: độ thập phân, mét)

• Dữ liệu - Một phiên bản được mô hình hóa cụ thể của Trái đất. Chúng đã được sửa đổi trong nhiều năm, vì vậy hãy đảm bảo rằng các lớp bản đồ của bạn đang sử dụng cùng một dữ liệu..

• Phép chiếu - Tham chiếu đến phương trình toán học được sử dụng để chiếu trái đất hình tròn lên một bề mặt phẳng (bản đồ).

Hãy nhớ rằng bạn có thể tóm tắt dữ liệu không gian mà không cần sử dụng các công cụ lập bản đồ được hiển thị bên dưới. Đôi khi một bảng đơn giản phân chia theo địa lý (ví dụ: huyện, quốc gia, v.v.) là tất cả những gì cần thiết!

Các loại bản đồ để trực quan hóa dữ liệu

Bản đồ Choropleth - một loại bản đồ trong đó màu sắc, sự đổ bóng, hoặc các họa tiết được sử dụng để thể hiện các vùng địa lý liên quan đến giá trị của một thuộc tính. Ví dụ: giá trị lớn hơn có thể được biểu thị bằng màu tối hơn giá trị nhỏ hơn. Loại bản đồ này đặc biệt hữu ích khi trực quan hóa một biến số và xem cách nó thay đổi trên các vùng hoặc khu vực địa chính trị xác định.

Bản đồ nhiệt mật độ trường hợp - một loại bản đồ trong đó màu sắc được sử dụng để thể hiện cường độ của một giá trị, tuy nhiên, nó không sử dụng các vùng hoặc ranh giới địa chính trị xác định để nhóm dữ liệu. Loại bản đồ này thường được sử dụng để hiển thị ‘điểm nóng’ hoặc các khu vực có mật độ hoặc tập trung điểm cao.

Bản đồ mật độ điểm - một loại bản đồ sử dụng các điểm để biểu thị các giá trị thuộc tính trong dữ liệu của bạn. Loại bản đồ này được sử dụng tốt nhất để trực quan hóa phân tán dữ liệu của bạn và scan một cách trực quan các cụm.

Bản đồ ký hiệu tỷ lệ (bản đồ ký hiệu chia độ) - một loại bản đồ tương tự như bản đồ choropleth, nhưng thay vì sử dụng màu sắc để biểu thị giá trị của một thuộc tính, nó sử dụng một ký hiệu (thường là một vòng tròn) liên quan đến giá trị. Ví dụ, một giá trị lớn hơn có thể được biểu thị bằng một ký hiệu lớn hơn một giá trị nhỏ hơn. Loại bản đồ này được sử dụng tốt nhất khi bạn muốn trực quan hóa kích thước hoặc số lượng dữ liệu của mình trên các vùng địa lý.

Bạn cũng có thể kết hợp một số loại trực quan hóa khác nhau để hiển thị các trường hợp địa lý phức tạp. Ví dụ, các trường hợp (điểm) trong bản đồ dưới đây được tô màu theo cơ sở y tế gần nhất của họ (xem phần chú thích). Các vòng tròn lớn hiển thị cơ sở y tế gần nhất trong khu vực ở một bán kính xác định, và các điểm màu đỏ tươi là những cơ sở y tế không nằm trong bất kỳ bán kính nào:

Lưu ý: Trọng tâm chính trong chương GIS này được dựa trên tình huống đáp ứng các vụ dịch tại thực địa. Do đó, nội dung của chương sẽ bao gồm các thao tác, hình ảnh hóa và phân tích dữ liệu không gian cơ bản.

Gọi package

Đoạn code này hiển thị việc gọi các package cần thiết cho các phân tích. Trong cuốn sách này, chúng tôi nhấn mạnh việc sử dụng hàm p_load() từ package pacman, giúp cài đặt các package cần thiết và gọi chúng ra để sử dụng. Bạn cũng có thể gọi các packages đã cài đặt với hàm library() của base R. Xem thêm chương R cơ bản để có thêm thông tin về các packages trong R.

pacman::p_load(
  rio,           # to import data
  here,          # to locate files
  tidyverse,     # to clean, handle, and plot the data (includes ggplot2 package)
  sf,            # to manage spatial data using a Simple Feature format
  tmap,          # to produce simple maps, works for both interactive and static maps
  janitor,       # to clean column names
  OpenStreetMap, # to add OSM basemap in ggplot map
  spdep          # spatial statistics
  ) 

Bạn có thể xem tổng quan về tất cả các package trong R xử lý dữ liệu không gian tại CRAN “Spatial Task View”.

Dữ liệu trường hợp mẫu

Với mục đích minh họa, chúng ta sẽ làm việc với một mẫu ngẫu nhiên gồm 1000 trường hợp từ bộ dữ liệu một vụ dịch Ebola mô phỏng có tên linelist (về mặt tính toán, việc làm việc với ít trường hợp hơn sẽ dễ hiển thị hơn trong sổ tay này). Để tiện the dõi, bấm để tải dữ liệu linelist “đã được làm sạch” (dưới dạng tệp .rds).

Vì chúng ta đang lấy một mẫu ngẫu nhiên của các trường hợp, nên kết quả của bạn có thể hơi khác so với những gì được minh họa ở đây khi bạn tự chạy code của mình.

Nhập dữ liệu vào bằng hàm import() từ package rio (nó xử lý nhiều loại tệp như .xlsx, .csv, .rds - xem thêm chương Nhập xuất dữ liệu để biết thêm chi tiết).

# import clean case linelist
linelist <- import("linelist_cleaned.rds")  

Tiếp theo, chúng ta chọn một mẫu ngẫu nhiên gồm 1000 hàng bằng hàm sample() từ base R.

# generate 1000 random row numbers, from the number of rows in linelist
sample_rows <- sample(nrow(linelist), 1000)

# subset linelist to keep only the sample rows, and all columns
linelist <- linelist[sample_rows,]

Bây giờ chúng ta muốn chuyển đổi bộ dữ liệu linelist đang là một dataframe, thành một đối tượng “sf” (spatial features - các đặc tính không gian). Do bộ dữ liệu linelist có hai cột “lon” và “lat” đại diện cho kinh độ và vĩ độ của nơi cư trú của từng trường hợp, nên việc chuyển đổi này khá dễ dàng.

Chúng ta sử dụng package sf (spatial features - các đặc tính không gian) và hàm st_as_sf() của nó để tạo đối tượng mới có tên linelist_sf. Đối tượng mới này về cơ bản giống với linelist, nhưng các cột lon và lat đã được chỉ định là cột tọa độ, và một hệ thống tham chiếu tọa độ (CRS) đã được gán khi các điểm được hiển thị. 4326 xác định tọa độ của chúng ta dựa trên Hệ thống trắc địa thế giới 1984 (WGS84) - là tiêu chuẩn cho tọa độ GPS.

# Create sf object
linelist_sf <- linelist %>%
     sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)

Dưới đây là cách bộ dữ liệu linelist gốc được hiển thị. Trong phần minh họa này, chúng ta sẽ chỉ sử dụng cột date_onset và geometry (được xây dựng từ các trường kinh độ và vĩ độ ở trên và là cột cuối cùng trong khung dữ liệu).

DT::datatable(head(linelist_sf, 10), rownames = FALSE, options = list(pageLength = 5, scrollX=T), class = 'white-space: nowrap' )

Shapefiles phân cấp hành chính

Sierra Leone: shapefiles phân cấp hành chính

Trước đó, chúng tôi đã tải xuống tất cả các phân cấp hành chính của Sierra Leone từ website của Humanitarian Data Exchange (HDX). Ngoài ra, bạn có thể tải xuống các dữ liệu này và tất cả các dữ liệu mẫu khác cho sổ tay này qua gói R của chúng tôi, như được giải thích trong trang Tải sách và dữ liệu.

Bây giờ chúng ta sẽ thực hiện những bước sau để lưu shapefile hành chính cấp 3 vào R:

1. Nhập shapefile
   
2. Làm sạch tên cột
   
3. Lọc các hàng để chỉ giữ lại các khu vực quan tâm

Để nhập một shapefile, chúng ta sử dụng hàm read_sf() từ package sf. Nó được cung cấp đường dẫn tệp thông qua hàm here(). - trong trường hợp này tệp nằm bên trong Dự án R của bạn tại thư mục “data”, rồi tới “gis”, và thư mục con “shp”, với tệp có tên “sle_adm3.shp” (xem chương Nhập xuất dữ liệu và Dự án R để biết thêm chi tiết). Bạn sẽ cần cung cấp đường dẫn tệp của riêng mình.

Tiếp theo, chúng ta sử dụng hàm clean_names() từ package janitor để chuẩn hóa tên các cột của shapefile. Chúng ta cũng sử dụng hàm filter() để chỉ giữ lại các hàng có tên admin2 là “Western Area Urban” hoặc “Western Area Rural”.

# ADM3 level clean
sle_adm3 <- sle_adm3_raw %>% 
  janitor::clean_names() %>% # standardize column names
  filter(admin2name %in% c("Western Area Urban", "Western Area Rural")) # filter to keep certain areas

Dưới đây, bạn có thể thấy shapefile trông như thế nào sau khi nhập và làm sạch. Cuộn sang phải to see how there are columns hành chính cấp 0 (quốc gia), hành chính cấp 1, hành chính cấp 2, và cuối cùng là hành chính cấp 3. Mỗi cấp độ có một tên và một mã định danh duy nhất “pcode”. Pcode mở rộng với mỗi sự tăng lên của cấp hành chính, vd. SL (Sierra Leone) -> SL04 (Western) -> SL0410 (Western Area Rural) -> SL040101 (Koya Rural).

Dữ liệu dân số

Sierra Leone: Dân số theo ADM3

Những dữ liệu này có thể tải xuống từ HDX (link tại đây) hoặc thông qua R package epirhandbook như đã được giải thích trong chương này. Chúng ta cũng sử dụng hàm import() để nạp tệp .csv. Chúng tôi cũng chuyển file được nhập tới hàm clean_names() để chuẩn hóa cú pháp tên cột.

# Population by ADM3
sle_adm3_pop <- import(here("data", "gis", "population", "sle_admpop_adm3_2020.csv")) %>%
  janitor::clean_names()

Bộ dữ liệu dân số trông sẽ như bên dưới. Cuộn sang bên phải để xem các cột: dân số nam (male population), dân số nữ (female populaton), tổng dân số (total population), và dân số theo từng nhóm tuổi.

Các cơ sở y tế

Sierra Leone: Dữ liệu cơ sở y tế từ OpenStreetMap

Một lần nữa, bạn có thể tải xuống thông tin vị trí của các cơ sở y tế từ HDX tại đây hoặc thông qua hướng dẫn trong chương Tải sách và dữ liệu.

Chúng ta nhập shapefile tọa độ điểm các cơ sở với hàm read_sf(), làm sạch lại tên cột, sau đó lọc để chỉ giữ lại các điểm được gắn thẻ là “hospital”, “clinic”, hoặc “doctors”.

# OSM health facility shapefile
sle_hf <- sf::read_sf(here("data", "gis", "shp", "sle_hf.shp")) %>% 
  janitor::clean_names() %>%
  filter(amenity %in% c("hospital", "clinic", "doctors"))

Dưới đây là dataframe kết quả - Cuộn phải để xem tên cơ sở và tọa độ geometry.

Điểm trong đa giác

Gán đơn vị hành chính cho các trường hợp theo không gian

Sau đây là một câu hỏi hóc búa thú vị: linelist không chứa bất kỳ thông tin nào về các đơn vị hành chính của các trường hợp. Mặc dù lý tưởng nhất là thu thập thông tin như vậy trong giai đoạn thu thập dữ liệu ban đầu, chúng ta cũng có thể gán các đơn vị hành chính cho các trường hợp riêng lẻ dựa trên mối quan hệ không gian của chúng (tức là điểm giao với một đa giác)..

Dưới đây, chúng ta sẽ giao các vị trí (điểm) theo không gian với địa giới hành chính cấp 3 (đa giác):

1. Bắt đầu với linelist (các điểm)
   
2. Nối theo không gian tới các địa giới, thiết lập kiểu nối là “st_intersects”
   
3. Sử dụng hàm select() để chỉ giữ lại một số cột địa giới hành chính mới

linelist_adm <- linelist_sf %>%
  
  # join the administrative boundary file to the linelist, based on spatial intersection
  sf::st_join(sle_adm3, join = st_intersects)

Tất cả các cột từ sle_adms đã được thêm vào linelist! Mỗi trường hợp bây giờ có các cột thể hiện chi tiết cấp hành chính mà nó nằm trong đó. Trong ví dụ này, chúng ta chỉ muốn giữ lại hai trong số các cột mới (địa giới hành chính cấp 3), vì vậy chúng dùng hàm select() để chọn tên các cột cũ và chỉ hai cột bổ sung quan tâm:

linelist_adm <- linelist_sf %>%
  
  # join the administrative boundary file to the linelist, based on spatial intersection
  sf::st_join(sle_adm3, join = st_intersects) %>% 
  
  # Keep the old column names and two new admin ones of interest
  select(names(linelist_sf), admin3name, admin3pcod)

Dưới đây, chỉ với mục đích hiển thị, bạn có thể thấy mười trường hợp đầu tiên và địa giới hành chính cấp 3 của chúng đã được gắn với nhau, dựa trên vị trí giao nhau trong không gian giữa các điểm với đa giác.

# Now you will see the ADM3 names attached to each case
linelist_adm %>% select(case_id, admin3name, admin3pcod)

## Simple feature collection with 1000 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -13.27052 ymin: 8.448185 xmax: -13.20545 ymax: 8.490177
## Geodetic CRS:  WGS 84
## First 10 features:
##      case_id     admin3name admin3pcod                   geometry
## 1814  1a1418 Mountain Rural   SL040102 POINT (-13.24142 8.455619)
## 2693  6f38aa Mountain Rural   SL040102    POINT (-13.22 8.450453)
## 1875  cd9a93        East II   SL040204 POINT (-13.21458 8.483898)
## 4450  3f9246       East III   SL040205 POINT (-13.20893 8.480541)
## 4160  191866       West III   SL040208 POINT (-13.25953 8.456642)
## 3281  cfeed3 Mountain Rural   SL040102 POINT (-13.21144 8.460869)
## 2836  237167      Central I   SL040201 POINT (-13.23522 8.485997)
## 4379  c104fd Mountain Rural   SL040102 POINT (-13.21395 8.461462)
## 2460  11945b       West III   SL040208 POINT (-13.25969 8.457996)
## 3706  b69c7c        West II   SL040207 POINT (-13.23563 8.471182)

Bây giờ chúng ta có thể mô tả các trường hợp theo đơn vị hành chính - điều mà chúng ta đã không thể làm trước khi nối theo không gian!

# Make new dataframe containing counts of cases by administrative unit
case_adm3 <- linelist_adm %>%          # begin with linelist with new admin cols
  as_tibble() %>%                      # convert to tibble for better display
  group_by(admin3pcod, admin3name) %>% # group by admin unit, both by name and pcode 
  summarise(cases = n()) %>%           # summarize and count rows
  arrange(desc(cases))                 # arrange in descending order

case_adm3

## # A tibble: 10 x 3
## # Groups:   admin3pcod [10]
##    admin3pcod admin3name     cases
##    <chr>      <chr>          <int>
##  1 SL040102   Mountain Rural   275
##  2 SL040208   West III         224
##  3 SL040207   West II          185
##  4 SL040204   East II          104
##  5 SL040201   Central I         65
##  6 SL040203   East I            51
##  7 SL040206   West I            43
##  8 SL040202   Central II        30
##  9 SL040205   East III          19
## 10 <NA>       <NA>               4

Chúng ta cũng có thể tạo một biểu đồ cột về số lượng trường hợp theo đơn vị hành chính.

Trong ví dụ này, chúng ta bắt đầu hàm ggplot() với bộ dữ liệu linelist_adm, trong đó chúng ta có thể áp dụng các hàm làm việc với factor như fct_infreq() để sắp xếp các cột theo tần suất (xem chương Factors để biết thêm các mẹo).

ggplot(
    data = linelist_adm,                       # begin with linelist containing admin unit info
    mapping = aes(
      x = fct_rev(fct_infreq(admin3name))))+ # x-axis is admin units, ordered by frequency (reversed)
  geom_bar()+                                # create bars, height is number of rows
  coord_flip()+                              # flip X and Y axes for easier reading of adm units
  theme_classic()+                           # simplify background
  labs(                                      # titles and labels
    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"
  )

Hàng xóm gần nhất

Tìm cơ sở y tế gần nhất / khu vực cung cấp dịch vụ y tế

Có thể sẽ hữu ích nếu biết các cơ sở y tế nằm ở đâu trong mối liên quan đến các điểm nóng về dịch bệnh.

Chúng ta có thể sử dụng phương pháp nối st_nearest_feature trong hàm st_join() (sf package) để trực quan hóa cơ sở y tế gần nhất với các trường hợp bệnh nhân.

1. Chúng ta bắt đầu với shapefile linelist có tên linelist_sf
   
2. Chúng ta nối theo không gian với đối tượng sle_hf, chứa thông tin về vị trí của các cơ sở y tế và phòng khám (các điểm)

# Closest health facility to each case
linelist_sf_hf <- linelist_sf %>%                  # begin with linelist shapefile  
  st_join(sle_hf, join = st_nearest_feature) %>%   # data from nearest clinic joined to case data 
  select(case_id, osm_id, name, amenity) %>%       # keep columns of interest, including id, name, type, and geometry of healthcare facility
  rename("nearest_clinic" = "name")                # re-name for clarity

Chúng ta có thể thấy bên dưới (50 hàng đầu tiên) rằng mỗi trường hợp hiện đã có dữ liệu về phòng khám/bệnh viện gần nhất

Chúng ta có thể thấy rằng “Den Clinic” là cơ sở y tế gần nhất với khoảng ~30% các trường hợp.

# Count cases by health facility
hf_catchment <- linelist_sf_hf %>%   # begin with linelist including nearest clinic data
  as.data.frame() %>%                # convert from shapefile to dataframe
  count(nearest_clinic,              # count rows by "name" (of clinic)
        name = "case_n") %>%         # assign new counts column as "case_n"
  arrange(desc(case_n))              # arrange in descending order

hf_catchment                         # print to console

##                          nearest_clinic case_n
## 1                            Den Clinic    383
## 2       Shriners Hospitals for Children    330
## 3         GINER HALL COMMUNITY HOSPITAL    173
## 4                             panasonic     54
## 5 Princess Christian Maternity Hospital     23
## 6                     ARAB EGYPT CLINIC     20
## 7                  MABELL HEALTH CENTER      9
## 8                                  <NA>      8

Để trực quan hóa kết quả, chúng ta có thể sử dụng tmap - lúc này interactive mode sẽ xem dễ dàng hơn

tmap_mode("view")   # set tmap mode to interactive  

# plot the cases and clinic points 
tm_shape(linelist_sf_hf) +            # plot cases
  tm_dots(size=0.08,                  # cases colored by nearest clinic
          col='nearest_clinic') +    
tm_shape(sle_hf) +                    # plot clinic facilities in large black dots
  tm_dots(size=0.3, col='black', alpha = 0.4) +      
  tm_text("name") +                   # overlay with name of facility
tm_view(set.view = c(-13.2284, 8.4699, 13), # adjust zoom (center coords, zoom)
        set.zoom.limits = c(13,14))+
tm_layout(title = "Cases, colored by nearest clinic")

Vùng đệm

Chúng ta cũng có thể tìm hiểu xem có bao nhiêu trường hợp nằm trong khoảng cách đi bộ 2,5 km (~30 phút) từ cơ sở y tế gần nhất.

Chú ý: Để tính toán khoảng cách chính xác hơn, tốt hơn nên chiếu lại đối tượng sf của bạn lên hệ thống chiếu bản đồ địa phương tương ứng chẳng hạn như UTM (Trái đất được chiếu lên bề mặt phẳng). Trong ví dụ này, để đơn giản hơn, chúng ta sẽ dựa vào Hệ tọa độ địa lý của Hệ thống trắc địa thế giới (WGS84) (Trái đất được biểu diễn trong một bề mặt hình cầu/tròn, do đó các đơn vị được tính bằng độ thập phân). Chúng ta sẽ sử dụng quy đổi chung là: 1 độ thập phân = ~111km.

Xem thêm thông tin về phép chiếu bản đồ và hệ tọa độ tại bài viết này esri article. Blog này nói về các loại phép chiếu bản đồ khác nhau và cách người ta có thể chọn phép chiếu phù hợp tùy thuộc vào khu vực quan tâm và bối cảnh của bản đồ/phân tích của bạn.

Đầu tiên, tạo một vùng đệm hình tròn với bán kính ~2.5km xung quanh mỗi cơ sở y tế. Điều này được thực hiện với hàm st_buffer() trong package tmap. Bởi vì đơn vị của bản đồ là kinh/vĩ độ thập phân, đó là cách “0.02” được diễn giải. Nếu hệ tọa độ bản đồ của bạn tính bằng mét, thì số đó phải được cung cấp bằng mét.

sle_hf_2k <- sle_hf %>%
  st_buffer(dist=0.02)       # decimal degrees translating to approximately 2.5km 

Sau đây chúng ta vẽ biểu đồ của chính các vùng đệm, với:

tmap_mode("plot")
# Create circular buffers
tm_shape(sle_hf_2k) +
  tm_borders(col = "black", lwd = 2)+
tm_shape(sle_hf) +                    # plot clinic facilities in large red dots
  tm_dots(size=0.3, col='black')      

Thứ hai, chúng ta giao các vùng đệm này với các trường hợp (điểm) bằng cách sử dụng hàm st_join() và kiểu nối là st_intersects. Tức là, dữ liệu từ vùng đệm được nối với các điểm mà chúng giao với nhau.

# Intersect the cases with the buffers
linelist_sf_hf_2k <- linelist_sf_hf %>%
  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)

Bây giờ chúng ta có thể đếm kết quả: có nrow(linelist_sf_hf_2k[is.na(linelist_sf_hf_2k$osm_id.y),]) trường hợp trong số 1000 trường hợp không giao nhau với bất kỳ vùng đệm nào (giá trị đó bị thiếu), và do đó họ sống ở nơi cách nhiều hơn 30 phút đi bộ tới cơ sở y tế gần nhất.

# Cases which did not get intersected with any of the health facility buffers
linelist_sf_hf_2k %>% 
  filter(is.na(osm_id.y)) %>%
  nrow()

## [1] 203

Chúng ta có thể trực quan hóa kết quả sao cho các trường hợp không giao nhau với bất kỳ vùng đệm nào sẽ xuất hiện với màu đỏ.

tmap_mode("view")

# First display the cases in points
tm_shape(linelist_sf_hf) +
  tm_dots(size=0.08, col='nearest_clinic') +

# plot clinic facilities in large black dots
tm_shape(sle_hf) +                    
  tm_dots(size=0.3, col='black')+   

# Then overlay the health facility buffers in polylines
tm_shape(sle_hf_2k) +
  tm_borders(col = "black", lwd = 2) +

# Highlight cases that are not part of any health facility buffers
# in red dots  
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))+

# add title  
tm_layout(title = "Cases by clinic catchment area")

Các hàm nối theo không gian khác

Các giá trị tùy chọn cho đối số join bao gồm (lấy từ tài liệu này)

• 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

OpenStreetMap

Dưới đây, chúng tôi mô tả cách lấy được bản đồ cơ sở cho bản đồ ggplot2 bằng cách sử dụng các tính năng của OpenStreetMap. Các phương pháp thay thế bao gồm sử dụng ggmap, yêu cầu bạn đăng ký miễn phí với Google (chi tiết).

OpenStreetMap là một dự án hợp tác nhằm tạo ra một bản đồ thế giới có thể chỉnh sửa miễn phí. Dữ liệu vị trí địa lý nền tảng ví dụ: vị trí của các thành phố, đường xá, đặc điểm tự nhiên, sân bay, trường học, bệnh viện, đường xá, v.v.) được coi là đầu ra chính của dự án.

Đầu tiên, chúng ta gọi package OpenStreetMap để lấy bản đồ cơ sở.

Sau đó chúng ta tạo đối tượng map, được xác định bằng cách sử dụng hàm openmap() từ package OpenStreetMap (tài liệu). Chúng ta cung cấp những thông tin sau:

• upperLeft và lowerRight Hai cặp tọa độ xác định giới hạn của ô bản đồ cơ sở

  • Trong trường hợp này, chúng tôi đã đưa giá trị tối đa và tối thiểu từ các hàng trong linelist, vì vậy bản đồ sẽ tương tác động với dữ liệu
    

• zoom = (nếu null nó sẽ được xác định tự động)
  

• type = loại bản đồ cơ sở nào - chúng tôi đã liệt kê một số khả năng ở đây và code hiện đang sử dụng cái đầu tiên ([1]) “osm”
  

• mergeTiles = chúng tôi đã chọn TRUE để tất cả các lớp nền được hợp nhất thành một

# load package
pacman::p_load(OpenStreetMap)

# Fit basemap by range of lat/long coordinates. Choose tile type
map <- OpenStreetMap::openmap(
  upperLeft = c(max(linelist$lat, na.rm=T), max(linelist$lon, na.rm=T)),   # limits of basemap tile
  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])

Nếu chúng ta vẽ bản đồ cơ sở này ngay bây giờ, sử dụng hàm autoplot.OpenStreetMap() từ package OpenStreetMap, bạn sẽ thấy rằng các đơn vị trên các trục không phải là kinh/vĩ độ. Nó đang sử dụng một hệ tọa độ khác. Để hiển thị chính xác nơi cứ trú các trường hợp (được lưu trữ theo vĩ độ/kinh độ), điều này phải được thay đổi.

autoplot.OpenStreetMap(map)

Vì vậy, chúng ta muốn chuyển đổi bản đồ thành vĩ độ/kinh độ với hàm openproj() từ package OpenStreetMap. Chúng ta cung cấp bản đồ cơ sở map và cũng cung cấp Hệ thống tham chiếu tọa độ (CRS) mà chúng ta muốn. Chúng ta thực hiện điều này bằng cách cung cấp chuỗi ký tự “proj.4” cho phép chiếu WGS 1984, nhưng bạn cũng có thể cung cấp CRS theo những cách khác. (xem trang này để hiểu rõ hơn chuỗi proj.4 là gì)

# Projection WGS84
map_latlon <- openproj(map, projection = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

Bây giờ khi chúng ta tạo biểu đồ, chúng ta thấy rằng dọc theo các trục là tọa độ vĩ độ và kinh độ. Hệ tọa độ đã được chuyển đổi. Bây giờ các trường hợp của chúng ta sẽ được biểu diễn chính xác nếu được phủ lên!

# Plot map. Must use "autoplot" in order to work with ggplot
autoplot.OpenStreetMap(map_latlon)

Xem các hướng dẫn tại đây và tại đây để biết thêm chi tiết.

Bản đồ nhiệt chuỗi thời gian

Bản đồ nhiệt mật độ ở trên cho thấy các trường hợp tích lũy. Chúng ta có thể khảo sát các vụ dịch theo thời gian và không gian bằng cách faceting bản đồ nhiệt dựa trên tháng khởi phát triệu chứng, như được rút ra từ linelist.

Chúng ta bắt đầu với linelist, tạo một cột mới với Năm và Tháng khởi phát. Hàm format() trong base R thay đổi cách hiển thị ngày. Trong trường hợp này, chúng ta muốn “YYYY-MM”..

# Extract month of onset
linelist <- linelist %>% 
  mutate(date_onset_ym = format(date_onset, "%Y-%m"))

# Examine the values 
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 2015-03 2015-04    <NA> 
##      12      20      27      75     185     191     122      98      70      63      57      33      47

Bây giờ, chúng ta chỉ cần đơn giản facet thông qua ggplot2 vào bản đồ nhiệt mật độ. facet_wrap() được áp dụng, sử dụng cột mới làm hàng. Chúng ta thiết lập số lượng cột được facet là 4 cho rõ ràng.

# packages
pacman::p_load(OpenStreetMap, tidyverse)

# begin with the basemap
autoplot.OpenStreetMap(map_latlon)+
  
  # add the density plot
  ggplot2::stat_density_2d(
        data = linelist,
        aes(
          x = lon,
          y = lat,
          fill = ..level..,
          alpha = ..level..),
        bins = 10,
        geom = "polygon",
        contour_var = "count",
        show.legend = F) +                          
  
  # specify color scale
  scale_fill_gradient(low = "black", high = "red")+
  
  # labels 
  labs(x = "Longitude",
       y = "Latitude",
       title = "Distribution of cumulative cases over time")+
  
  # facet the plot by month-year of onset
  facet_wrap(~ date_onset_ym, ncol = 4)               

Mối quan hệ không gian

Trước khi có thể tính toán bất kỳ số liệu thống kê không gian nào, chúng ta cần xác định mối quan hệ giữa các đối tượng địa lý trong dữ liệu của mình. Có nhiều cách để khái niệm hóa các mối quan hệ không gian, nhưng một mô hình đơn giản và thường được áp dụng để sử dụng là mô hình adjancy - liền kề - cụ thể là chúng ta mong đợi mối quan hệ địa lý giữa các khu vực có chung biên giới hoặc “láng giềng” với nhau..

Chúng ta có thể định lượng các mối quan hệ liền kề giữa các vùng địa giới hành chính trong dữ liệu sle_adm3 mà chúng ta đang sử dụng với package spdep. Chúng ta sẽ cụ thể sự tiếp giáp queen, có nghĩa là các khu vực sẽ là hàng xóm của nhau nếu chúng có chung ít nhất một điểm dọc theo biên giới của chúng. Phương pháp thay thế có thể là sự tiếp giáp rook, đòi hỏi các khu vực phải chia sẻ một cạnh - trong trường hợp của chúng ta, với các đa giác không đều, sự phân biệt là không đáng kể, nhưng trong một số trường hợp, sự lựa chọn giữa queen và rook có thể có ảnh hưởng.

sle_nb <- spdep::poly2nb(sle_adm3_dat, queen=T) # create neighbors 
sle_adjmat <- spdep::nb2mat(sle_nb)    # create matrix summarizing neighbor relationships
sle_listw <- spdep::nb2listw(sle_nb)   # create listw (list of weights) object -- we will need this later

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)

Ma trận được in ở trên hiển thị mối quan hệ giữa 9 vùng trong dữ liệu sle_adm3 của chúng ta. Điểm 0 cho biết hai vùng không phải là láng giềng, trong khi bất kỳ giá trị nào khác 0 cho biết mối quan hệ láng giềng. Các giá trị trong ma trận được chia tỷ lệ để mỗi vùng có tổng trọng số hàng là 1.

Một cách tốt hơn để trực quan hóa những mối quan hệ láng giềng này là vẽ biểu đồ chúng:

plot(sle_adm3_dat$geometry) +                                           # plot region boundaries
  spdep::plot.nb(sle_nb,as(sle_adm3_dat, 'Spatial'), col='grey', add=T) # add neighbor relationships

Chúng ta đã sử dụng phương pháp tiếp cận adjacency để xác định các đa giác hàng xóm; những người hàng xóm mà chúng ta đã xác định đôi khi cũng được gọi là những người hàng xóm dựa trên sự tiếp giáp. Nhưng đây chỉ là một cách để lựa chọn các khu vực dự kiến có mối quan hệ địa lý. Các cách tiếp cận thay thế phổ biến nhất để xác định các mối quan hệ địa lý tạo ra những người hàng xóm dựa trên khoảng cách; ngắn gọn thì chúng là:

• K-nearest neighbors - Dựa trên khoảng cách giữa các tâm (trung tâm có trọng số địa lý của mỗi vùng đa giác), hãy chọn n vùng gần nhất làm vùng lân cận. Ngưỡng khoảng cách tối đa cũng có thể được chỉ định. Trong spdep, bạn có thể sử dụng hàm knearneigh() (xem tài liệu).

• Distance threshold neighbors - Chọn tất cả hàng xóm trong một ngưỡng khoảng cách. Trong spdep, các mối quan hệ láng giềng này có thể được xác định bằng cách sử dụng hàm dnearneigh() (xem tài liệu).

Tự tương quan không gian

Định luật địa lý đầu tiên được trích dẫn của Tobler tuyên bố rằng “mọi thứ đều liên quan đến mọi thứ khác, nhưng những thứ ở gần có liên quan hơn những thứ ở xa”. Trong dịch tễ học, điều này thường có nghĩa là nguy cơ về một kết quả sức khỏe cụ thể ở một vùng nhất định tương tự với các vùng lân cận hơn là những vùng xa. Khái niệm này đã được chính thức hóa dưới dạng tự tương quan không gian - thuộc tính thống kê mà các đối tượng địa lý có giá trị tương tự được nhóm lại với nhau trong không gian. Các phép đo thống kê về tự tương quan không gian có thể được sử dụng để định lượng mức độ phân cụm không gian trong dữ liệu của bạn, xác định vị trí xảy ra phân cụm, và xác định các xu hướng chung về tự tương quan không gian giữa các biến khác nhau trong dữ liệu của bạn. Phần này giới thiệu tổng quan về một số thước đo phổ biến của tự tương quan không gian và cách tính chúng trong R.

Moran’s I - Đây là thống kê tóm tắt mang tính toàn cục (global) về mối tương quan giữa giá trị của một biến trong một vùng và giá trị của cùng một biến ở các vùng lân cận. Thống kê Moran’s I thường dao động từ -1 tới 1. Giá trị 0 cho thấy không có tương quan không gian, trong khi các giá trị gần 1 hoặc -1 cho thấy tự tương quan không gian (các giá trị tương tự gần nhau) hoặc phân tán không gian (các giá trị không tương tự ở gần nhau) mạnh hơn.

Ví dụ: chúng ta sẽ tính toán thống kê Moran’s I để định lượng tự tương quan theo không gian trong các trường hợp Ebola mà chúng ta đã mapping trước đó (hãy nhớ rằng đây là một tập con các trường hợp từ một vụ dịch mô phỏng trong dataframe linelist). Package spdep có hàm moran.test có thể giúp chúng ta thực hiện phép tính toán này:

moran_i <-spdep::moran.test(sle_adm3_dat$cases,    # numeric vector with variable of interest
                            listw=sle_listw)       # listw object summarizing neighbor relationships

moran_i                                            # print results of Moran's I test

## 
##  Moran I test under randomisation
## 
## data:  sle_adm3_dat$cases  
## weights: sle_listw    
## 
## Moran I statistic standard deviate = 1.6554, p-value = 0.04892
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.22298240       -0.12500000        0.04418648

Đầu ra từ hàm moran.test() cho chúng ta thấy kết quả thống kê Moran I statistic từ hàm round(moran_i$estimate[1],2). Điều này cho thấy sự hiện diện của tự tương quan không gian trong dữ liệu của chúng ta - cụ thể là những khu vực có số ca mắc Ebola tương tự nhau thì có khả năng ở gần nhau. Giá trị p cung cấp bởi kiểm định moran.test() được tạo ra bằng cách so sánh với kỳ vọng trong giả thuyết không là không có tự tương quan không gian, và có thể được sử dụng nếu bạn cần báo cáo kết quả của một kiểm định giả thuyết chính thức.

Local Moran’s I - Chúng ta có thể phân tách thống kê Moran’s I (toàn cục - global) được tính toán ở trên để xác định tự tương quan cục bộ (local); nghĩa là, để xác định các cụm cụ thể trong dữ liệu của chúng ta. Thống kê này, đôi khi còn được gọi là Local Indicator of Spatial Association (LISA), tóm tắt mức độ tự tương quan không gian xung quanh mỗi vùng riêng lẻ. Nó có thể hữu ích để tìm các điểm “nóng - hot” và “lạnh - cold” trên bản đồ.

Để hiển thị một ví dụ, chúng ta có thể tính toán và lập bản đồ Local Moran’s I cho các trường hợp Ebola ở trên với hàm local_moran() từ package spdep:

# calculate local Moran's I
local_moran <- spdep::localmoran(                  
  sle_adm3_dat$cases,                              # variable of interest
  listw=sle_listw                                  # listw object with neighbor weights
)

# join results to sf data
sle_adm3_dat<- cbind(sle_adm3_dat, local_moran)    

# plot map
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")

Getis-Ord Gi* - Đây là một thống kê khác thường được sử dụng để phân tích điểm nóng - hotspot analysis; mức độ phổ biến của thống kê này liên quan đến ứng dụng của nó trong công cụ Phân tích điểm nóng của ArcGIS. Nó dựa trên giả định rằng thông thường, sự khác biệt về giá trị của một biến giữa các vùng lân cận phải tuân theo phân phối chuẩn. Nó sử dụng cách tiếp cận z-score để xác định các vùng có giá trị cao hơn đáng kể (điểm nóng) hoặc thấp hơn đáng kể (điểm lạnh) của một biến cụ thể, so với các vùng lân cận của chúng..

Chúng ta có thể tính toán và lập bản đồ chỉ số Gi* sử dụng hàm localG() từ package spdep:

# Perform local G analysis
getis_ord <- spdep::localG(
  sle_adm3_dat$cases,
  sle_listw
)

# join results to sf data
sle_adm3_dat$getis_ord <- getis_ord

# plot map
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")

Như bạn có thể thấy, bản đồ được tạo ra bởi Getis-Ord Gi* trông hơi khác so với bản đồ của Local Moran’s I được tạo ra trước đó. Điều này phản ánh rằng phương pháp được sử dụng để tính toán hai số liệu thống kê này hơi khác nhau; bạn nên sử dụng cái nào tùy thuộc vào trường hợp sử dụng cụ thể của bạn và câu hỏi nghiên cứu quan tâm.

Lee’s L test - Đây là một kiểm định thống kê về mối tương quan không gian giữa hai biến. Nó cho phép bạn kiểm định xem xu hướng không gian của một biến x có tương tự với xu hướng không gian của một biến y khác hay không, với giả định là chúng có mối tương quan không gian với nhau

Để đưa ra một ví dụ, hãy kiểm định liệu rằng xu hướng không gian các trường hợp Ebola từ một vụ dịch mô phỏng có tương quan với xu hướng không gian của dân số hay không. Để bắt đầu, chúng ta cần có một biến population trong bộ dữ liệu sle_adm3 của mình. Chúng ta có thể sử dụng biến total từ bộ dữ liệu sle_adm3_pop mà chúng ta đã tải trước đó.

sle_adm3_dat <- sle_adm3_dat %>% 
  rename(population = total)                          # rename 'total' to 'population'

Chúng ta có thể trực quan hóa nhanh xu hướng không gian của hai biến cạnh nhau, để xem liệu chúng có tương tự nhau hay không:

tmap_mode("plot")

cases_map <- tm_shape(sle_adm3_dat) + tm_polygons("cases") + tm_layout(main.title="Cases")
pop_map <- tm_shape(sle_adm3_dat) + tm_polygons("population") + tm_layout(main.title="Population")

tmap_arrange(cases_map, pop_map, ncol=2)   # arrange into 2x1 facets

Một cách trực quan, các xu hướng có vẻ không giống nhau. Chúng ta có thể sử dụng hàm lee.test() trong package spdep để kiểm định liệu xu hướng tự tương quan không gian giữa hai biến số có liên quan có ý nghĩa thống kê với nhau hay không. Thống kê L sẽ gần tới 0 nếu không có mối tương quan giữa các xu hướng, gần bằng 1 nếu có mối tương quan dương mạnh (tức là các xu hướng tương tự nhau) và gần bằng -1 nếu có mối tương quan âm mạnh (tức là các xu hướng là ngược nhau).

lee_test <- spdep::lee.test(
  x=sle_adm3_dat$cases,          # variable 1 to compare
  y=sle_adm3_dat$population,     # variable 2 to compare
  listw=sle_listw                # listw object with neighbor weights
)

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 = -1.0254, p-value = 0.8474
## alternative hypothesis: greater
## sample estimates:
## Lee's L statistic       Expectation          Variance 
##       -0.15912924       -0.04520388        0.01234284

Kết quả từ thống kê Lee’s L cho hai biến là round(lee_test$estimate[1],2), chỉ ra rằng chúng có mối tương quan âm yếu. Điều này khẳng định đánh giá trực quan của chúng ta rằng xu hướng các trường hợp bệnh và dân số không liên quan đến nhau cũng như cung cấp bằng chứng rằng xu hướng không gian của các trường hợp không hoàn toàn là kết quả của mật độ dân số ở các khu vực có nguy cơ cao.

Thống kê Lee L có thể hữu ích để đưa ra các suy luận về mối quan hệ phân bố theo không gian giữa các biến; tuy nhiên, để mô tả bản chất của mối quan hệ giữa hai biến một cách chi tiết hơn, hoặc điều chỉnh cho nhiễu, sẽ cần đến các kỹ thuật hồi quy không gian. Chúng sẽ được mô tả ngắn gọn trong phần sau.

Hồi quy không gian

Bạn có thể muốn đưa ra các suy luận thống kê về mối quan hệ giữa các biến trong dữ liệu không gian của mình. Trong những trường hợp này, sẽ hữu ích khi xem xét các kỹ thuật hồi quy không gian - nghĩa là, các cách tiếp cận hồi quy xem xét rõ ràng tổ chức không gian của các đơn vị trong dữ liệu của bạn. Một số lý do mà bạn có thể cần xem xét các mô hình hồi quy không gian, thay vì các mô hình hồi quy tiêu chuẩn như GLM, bao gồm:

• Các mô hình hồi quy tiêu chuẩn giả định rằng các phần dư là độc lập với nhau. Khi có hiện tượng tự tương quan không gian mạnh, các phần dư của mô hình hồi quy chuẩn cũng có khả năng tự tương quan không gian, do đó vi phạm giả định này. Điều này có thể dẫn đến các vấn đề trong việc giải thích kết quả mô hình, trong trường hợp đó, một mô hình không gian sẽ được ưu tiên hơn.

• Các mô hình hồi quy cũng thường giả định rằng tác động của một biến x là không đổi đối với tất cả các quan sát. Trong trường hợp không đồng nhất về không gian - spatial heterogenity, các tác động mà chúng ta muốn ước tính có thể thay đổi theo không gian và chúng ta có thể quan tâm đến việc định lượng những khác biệt đó. Trong trường hợp này, các mô hình hồi quy không gian mang lại sự linh hoạt hơn cho việc ước lượng và giải thích các tác động.

Các chi tiết về phương pháp tiếp cận hồi quy không gian nằm ngoài phạm vi của sổ tay này. Thay vào đó, phần này sẽ cung cấp tổng quan về các mô hình hồi quy không gian phổ biến nhất và cách sử dụng của chúng, đồng thời giới thiệu cho bạn tài liệu tham khảo có thể sử dụng nếu bạn muốn khám phá thêm lĩnh vực này.

Các mô hình sai số không gian - Spatial error models - Các mô hình này giả định rằng sai số trên các đơn vị không gian có tương quan với nhau, trong trường hợp đó, dữ liệu sẽ vi phạm các giả định của mô hình OLS tiêu chuẩn. Các mô hình này đôi khi còn được biết đến với tên Các mô hình tự hồi quy đồng thời - simultaneous autoregressive (SAR) models. Mô hình được fit bằng cách sử dụng hàm errorsarlm() trong package spatialreg (các hàm hồi quy không gian là một phần của package spdep).

Các mô hình độ trễ không gian - Spatial lag models - Các mô hình này giả định rằng biến phụ thuộc của một vùng i không chỉ bị ảnh hưởng bởi giá trị của các biến độc lập trong i, mà còn bởi giá trị của các biến đó ở các vùng lân cận i. Tương tự như các mô hình sai số không gian, các mô hình độ trễ không gian cũng thường được biết đến với tên gọi Các mô hình tự hồi quy đồng thời. Mô hình được fit bằng cách sử dụng hàm lagsarlm() trong package spatialreg.

Package spdep chứa một số kiểm định hữu ích để quyết định lựa chọn giữa các mô hình OLS chuẩn, độ trễ không gian và sai số không gian. Các kiểm định này được gọi là chẩn đoán Lagrange Multiplier, có thể được sử dụng để xác định loại phụ thuộc không gian trong dữ liệu của bạn và chọn mô hình nào phù hợp nhất. Hàm lm.LMtests() có thể được sử dụng để tính toán tất cả các kiểm định Lagrange Multiplier. Anselin (1988) cũng cung cấp một sơ đồ khối hữu ích để quyết định sử dụng mô hình hồi quy không gian nào dựa trên kết quả của các kiểm định Lagrange Multiplier:

Các mô hình Bayesian phân tầng - Bayesian hierarchical models - Phương pháp tiếp cận Bayes thường được sử dụng cho một số ứng dụng trong phân tích không gian, phổ biến nhất là lập bản đồ dịch bệnh. Chúng được ưu tiên sử dụng trong các trường hợp dữ liệu trường hợp được phân phối thưa thớt (ví dụ trong trường hợp các outcome là hiếm gặp) hoặc “nhiễu” về mặt thống kê, vì chúng có thể được sử dụng để tạo ra các ước tính “mượt mà” về nguy cơ bệnh tật bằng cách tính đến quy trình không gian tiềm ẩn cơ bản. Điều này có thể giúp cải thiện chất lượng của các ước tính. Chúng cũng cho phép người điều tra cụ thể trước (thông qua việc lựa chọn prior) các xu hướng tương quan không gian phức tạp có thể tồn tại trong dữ liệu, mà có thể giải thích cho các sự biến động phụ thuộc hoặc không phụ thuộc không gian trong cả biến độc lập và phụ thuộc. Trong R, các mô hình Bayesian phân tầng có thể sử dụng thông qua package CARbayes (xem hưỡng dẫn) hoặc R-INLA (xem trang web và sách giáo khoa). R cũng có thể được sử dụng để gọi phần mềm bên ngoài thực hiện ước lượng Bayes, chẳng hạn như JAGS hoặc WinBUGS.