Local Moran and Local Getis-Ord Maps of OD Deaths by Town of Residence using the {rgeoda} package in R

Using the {rgeoda} package to create Local Moran and Local Getis-Ord Maps

Here are the steps for making a cluster map and a hot spots map using the {rgeoda} package:

1. Set working directory.

setwd("~/Documents/GIS Shapefiles/OD Drug Deaths Resident Town 2022")

2. Load libraries.

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

#install.packages("rgeoda")
library(rgeoda)

## Loading required package: digest

library(digest)
library(sf)

## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

3. Load spatial data from working directory.

OD_Drug_Deaths <- st_read("2022_OD_Drug_Deaths_by_Resident_Town.shp")

4. Plot OD Deaths (i.e.,‘COUNT’) variable.

plot(OD_Drug_Deaths["COUNT"], main = "Connecticut OD Deaths by Town of Residence, 2022")

5. Create a ‘Queen’ contiguity weights using the sf object ‘OD_Drug_Deaths’.

queen_w <- queen_weights(OD_Drug_Deaths)
summary(queen_w)

##                      name              value
## 1 number of observations:                169
## 2          is symmetric:                TRUE
## 3               sparsity: 0.0320717061727531
## 4        # min neighbors:                  1
## 5        # max neighbors:                 10
## 6       # mean neighbors:   5.42011834319527
## 7     # median neighbors:                  6
## 8           has isolates:              FALSE

Spatial Autocorrelations

Local Moran (Cluster Analysis)

The Local Moran statistic is a method to identify local clusters and local spatial outliers.

6. Examine the local Moran statistics of variable “COUNT” (Number of OD Drugs Deaths by Town of Residence):

lisa <- local_moran(queen_w,  OD_Drug_Deaths['COUNT'])

Create Local Moran Map

7. With the LISA results, we can make a local moran cluster map.

lisa_colors <- lisa_colors(lisa)
lisa_labels <- lisa_labels(lisa)
lisa_clusters <- lisa_clusters(lisa)

plot(st_geometry(OD_Drug_Deaths), 
     col=sapply(lisa_clusters, function(x){return(lisa_colors[[x+1]])}), 
     border = "#333333", lwd=0.2)
title(main = "Local Moran Map of OD Deaths by Town of Residence, 2022")
legend('topleft', legend = lisa_labels, fill = lisa_colors, border = "#eeeeee", cex=0.7)

8. Create a significance map that is associated with the local Moran map.

lisa_p <- lisa_pvalues(lisa)
p_labels <- c("Not significant", "p <= 0.05", "p <= 0.01", "p <= 0.001")
p_colors <- c("#eeeeee", "#84f576", "#53c53c", "#348124")
plot(st_geometry(OD_Drug_Deaths), 
     col=sapply(lisa_p, function(x){
       if (x <= 0.001) return(p_colors[4])
       else if (x <= 0.01) return(p_colors[3])
       else if (x <= 0.05) return (p_colors[2])
       else return(p_colors[1])
     }), 
     border = "#333333", lwd=0.2)
title(main = "Local Moran Map of OD Deaths by Town of Residence, 2022")
legend('topleft', legend = p_labels, fill = p_colors, border = "#eeeeee", cex=0.7)

Local Getis-Ord Statistics (Hot Spot Analysis)

A value larger than the mean suggests a high-high cluster or hot spot and a value smaller than the mean indicates a low-low cluster or cold spot.

9. We can call the function local_g() with the created Queen weights and the variable “COUNT” as input parameters.

localg_COUNT <- local_g(queen_w,  OD_Drug_Deaths['COUNT'])

10. With the local_g results, we can make a Local Getis-Ord hot spots map.

lisa_colors_g <- lisa_colors(localg_COUNT)
lisa_labels_g <- lisa_labels(localg_COUNT)
lisa_clusters_g <- lisa_clusters(localg_COUNT)

plot(st_geometry(OD_Drug_Deaths), 
     col=sapply(lisa_clusters_g, function(x){return(lisa_colors_g[[x+1]])}), 
     border = "#333333", lwd=0.2)
title(main = "Local Getis-Ord Map of OD Deaths by Town of Residence, 2022")
legend('topleft', legend = lisa_labels_g, fill = lisa_colors_g, border = "#eeeeee", cex=0.7)

11. Create a significance map that is associated with the Local Getis-Ord hot spots map.

lisa_p_g <- lisa_pvalues(localg_COUNT)
p_labels_g <- c("Not significant", "p <= 0.05", "p <= 0.01", "p <= 0.001")
p_colors_g <- c("#eeeeee", "#84f576", "#53c53c", "#348124")
plot(st_geometry(OD_Drug_Deaths), 
     col=sapply(lisa_p_g, function(x){
       if (x <= 0.001) return(p_colors_g[4])
       else if (x <= 0.01) return(p_colors_g[3])
       else if (x <= 0.05) return (p_colors_g[2])
       else return(p_colors_g[1])
     }), 
     border = "#333333", lwd=0.2)
title(main = "Local Getis-Ord Map of OD Deaths by Town of Residence, 2022")
legend('topleft', legend = p_labels_g, fill = p_colors_g, border = "#eeeeee", cex=0.7)

A.M.D.G.