Chapter 5 Spatial Analysis

5.1 Processing

rm(list = ls())
library(tidyverse)
library(spdep)
library(tmap)
library(sf)
library(leaflet)

dat_ci_map <- readRDS("./_dat/dat_ci_map_v3.rds") %>%
  filter(sample.size>=10) %>%
  filter(rii !=Inf & rii>0)
dat_edu_map <- readRDS("./_dat/dat_edu_map_v2.rds") %>%
  filter(sample.size>=10) %>%
  filter(rii !=Inf & rii>0)
source("./mal_ineq_fun.R") 

breaks <- c(-Inf, -2.58, -1.96, -1.65, 1.65, 1.96, 2.58, Inf)
labels <- c("Cold spot: 99% confidence", "Cold spot: 95% confidence", "Cold spot: 90% confidence", "PSU Not significant","Hot spot: 90% confidence", "Hot spot: 95% confidence", "Hot spot: 99% confidence")

5.2 Wealth

dat_ci_sp <- as(dat_ci_map, "Spatial")

coords<-coordinates(dat_ci_sp)
IDs<-row.names(as.data.frame(coords))
Neigh_nb<-knn2nb(knearneigh(coords, k=1, longlat = TRUE), row.names=IDs)
dsts<-unlist(nbdists(Neigh_nb,coords))
max_1nn<-max(dsts)
Neigh_kd1<-dnearneigh(coords,d1=0, d2=max_1nn, row.names=IDs)

self<-nb2listw(Neigh_kd1, style="W", zero.policy = T)

w.getis <- dat_ci_map %>%
  st_set_geometry(NULL) %>%
  nest() %>%
  crossing(var1 = c("ci_val", "sii", "rii")) %>%
  mutate(data_sub = map2(.x = data, .y = var1, .f = ~dplyr::select(.x, .y, country, cluster_number, 
                                                                   prev, lat, long) %>% 
                           rename(output = 1)),
         LISA = map(.x = data_sub, .f = ~localG(.x$output, self)),
         clust_LISA = map(.x = LISA, .f = ~cut(.x, include.lowest = TRUE, 
                                               breaks = breaks, labels = labels))) %>%
  dplyr::select(data_sub, var1, LISA, clust_LISA) %>%
  unnest()

5.2.1 W-CI

dat_w_map_ci <- w.getis %>%
  filter(var1 =="ci_val") %>%
  st_as_sf(coords = c("long", "lat"), crs = 4326) %>%
  filter(!is.na(LISA)) %>%
  dplyr::mutate(lat = sf::st_coordinates(.)[,2],
                long = sf::st_coordinates(.)[,1])

pal <- colorFactor("RdBu", dat_w_map_ci$clust_LISA, reverse = T)

dat_w_map_ci %>%
  leaflet() %>%
  addTiles() %>%
  addCircles(lng = ~long, lat = ~lat, color = ~ pal(clust_LISA),
             popup = paste("Country", dat_w_map_ci$country, "<br>",
                           "Cluster:", dat_w_map_ci$cluster_number, "<br>",
                           "Malaria Prevalence:", dat_w_map_ci$prev, "<br>",
                           "CI:", dat_w_map_ci$output, "<br>",
                           "LISA:", dat_w_map_ci$LISA)) %>%
  addLegend("bottomright",
            pal = pal, values = ~clust_LISA,
            title = "LISA - Wealth CI") %>%
  addScaleBar(position = c("bottomleft"))

5.2.2 W-SII

dat_w_map_sii <- w.getis %>%
  filter(var1 =="sii") %>%
  st_as_sf(coords = c("long", "lat"), crs = 4326) %>%
  filter(!is.na(LISA)) %>%
  dplyr::mutate(lat = sf::st_coordinates(.)[,2],
                long = sf::st_coordinates(.)[,1])

pal <- colorFactor("RdBu", dat_w_map_sii$clust_LISA, reverse = T)

dat_w_map_sii %>%
  leaflet() %>%
  addTiles() %>%
  addCircles(lng = ~long, lat = ~lat, color = ~ pal(clust_LISA),
             popup = paste("Country", dat_w_map_sii$country, "<br>",
                           "Cluster:", dat_w_map_sii$cluster_number, "<br>",
                           "Malaria Prevalence:", dat_w_map_sii$prev, "<br>",
                           "CI:", dat_w_map_sii$output, "<br>",
                           "LISA:", dat_w_map_sii$LISA)) %>%
  addLegend("bottomright",
            pal = pal, values = ~clust_LISA,
            title = "LISA - Wealth SII") %>%
  addScaleBar(position = c("bottomleft"))

5.2.3 W-RII

dat_w_map_rii <- w.getis %>%
  filter(var1 =="rii") %>%
  st_as_sf(coords = c("long", "lat"), crs = 4326) %>%
  filter(!is.na(LISA)) %>%
  dplyr::mutate(lat = sf::st_coordinates(.)[,2],
                long = sf::st_coordinates(.)[,1])

pal <- colorFactor("RdBu", dat_w_map_rii$clust_LISA, reverse = T)

dat_w_map_rii %>%
  leaflet() %>%
  addTiles() %>%
  addCircles(lng = ~long, lat = ~lat, color = ~ pal(clust_LISA),
             popup = paste("Country", dat_w_map_rii$country, "<br>",
                           "Cluster:", dat_w_map_rii$cluster_number, "<br>",
                           "Malaria Prevalence:", dat_w_map_rii$prev, "<br>",
                           "CI:", dat_w_map_rii$output, "<br>",
                           "LISA:", dat_w_map_rii$LISA)) %>%
  addLegend("bottomright",
            pal = pal, values = ~clust_LISA,
            title = "LISA - Wealth RII") %>%
  addScaleBar(position = c("bottomleft"))

5.3 Education

dat_edu_sp <- as(dat_edu_map, "Spatial")

coords<-coordinates(dat_edu_sp)
IDs<-row.names(as.data.frame(coords))
Neigh_nb<-knn2nb(knearneigh(coords, k=1, longlat = TRUE), row.names=IDs)
dsts<-unlist(nbdists(Neigh_nb,coords))
max_1nn<-max(dsts)
Neigh_kd1<-dnearneigh(coords,d1=0, d2=max_1nn, row.names=IDs)

self<-nb2listw(Neigh_kd1, style="W", zero.policy = T)

w.getis_e <- dat_edu_map %>%
  st_set_geometry(NULL) %>%
  nest() %>%
  crossing(var1 = c("ci_val", "sii", "rii")) %>%
  mutate(data_sub = map2(.x = data, .y = var1, .f = ~dplyr::select(.x, .y, country, cluster_number, 
                                                                   prev, lat, long) %>% 
                           rename(output = 1)),
         LISA = map(.x = data_sub, .f = ~localG(.x$output, self)),
         clust_LISA = map(.x = LISA, .f = ~cut(.x, include.lowest = TRUE, 
                                               breaks = breaks, labels = labels))) %>%
  dplyr::select(data_sub, var1, LISA, clust_LISA) %>%
  unnest()

5.3.1 E-CI

dat_e_map_ci <- w.getis_e %>%
  filter(var1 =="ci_val") %>%
  st_as_sf(coords = c("long", "lat"), crs = 4326) %>%
  filter(!is.na(LISA)) %>%
  dplyr::mutate(lat = sf::st_coordinates(.)[,2],
                long = sf::st_coordinates(.)[,1])

pal <- colorFactor("RdBu", dat_e_map_ci$clust_LISA, reverse = T)

dat_e_map_ci %>%
  leaflet() %>%
  addTiles() %>%
  addCircles(lng = ~long, lat = ~lat, color = ~ pal(clust_LISA),
             popup = paste("Country", dat_e_map_ci$country, "<br>",
                           "Cluster:", dat_e_map_ci$cluster_number, "<br>",
                           "Malaria Prevalence:", dat_e_map_ci$prev, "<br>",
                           "CI:", dat_e_map_ci$output, "<br>",
                           "LISA:", dat_e_map_ci$LISA)) %>%
  addLegend("bottomright",
            pal = pal, values = ~clust_LISA,
            title = "LISA - Education CI") %>%
  addScaleBar(position = c("bottomleft"))

5.3.2 E-SII

dat_e_map_sii <- w.getis_e %>%
  filter(var1 =="sii") %>%
  st_as_sf(coords = c("long", "lat"), crs = 4326) %>%
  filter(!is.na(LISA)) %>%
  dplyr::mutate(lat = sf::st_coordinates(.)[,2],
                long = sf::st_coordinates(.)[,1])

pal <- colorFactor("RdBu", dat_e_map_sii$clust_LISA, reverse = T)

dat_e_map_sii %>%
  leaflet() %>%
  addTiles() %>%
  addCircles(lng = ~long, lat = ~lat, color = ~ pal(clust_LISA),
             popup = paste("Country", dat_e_map_sii$country, "<br>",
                           "Cluster:", dat_e_map_sii$cluster_number, "<br>",
                           "Malaria Prevalence:", dat_e_map_sii$prev, "<br>",
                           "CI:", dat_e_map_sii$output, "<br>",
                           "LISA:", dat_e_map_sii$LISA)) %>%
  addLegend("bottomright",
            pal = pal, values = ~clust_LISA,
            title = "LISA - Education SII") %>%
  addScaleBar(position = c("bottomleft"))