Spatial Stratified Heterogeneity Test

This vignette explains how to do Spatial Stratified Heterogeneity Test in spEcula package.

Schematic diagram of Spatial Stratified Heterogeneity Test

Load package and pre-processing data.

library(sf)
library(tidyverse)
library(spEcula)

See layers in NTDs.gpkg:

ntdspath = system.file("extdata", "NTDs.gpkg",package = 'spEcula')
st_layers(ntdspath)
## Driver: GPKG 
## Available layers:
##   layer_name geometry_type features fields crs_name
## 1    disease       Polygon      189      2  unknown
## 2  watershed       Polygon        9      2  unknown
## 3  elevation Multi Polygon        7      2  unknown
## 4   soiltype Multi Polygon        6      2  unknown

In NTDs.gpkg, disease is a dependent variable, which is a continuous numerical variable, while others are independent and discrete variables.

Now we need to put these layers together:

watershed = read_sf(ntdspath,layer = 'watershed')
elevation = read_sf(ntdspath,layer = 'elevation')
soiltype = read_sf(ntdspath,layer = 'soiltype')
disease = read_sf(ntdspath,layer = 'disease')

Plot them together:

library(cowplot)

f1 = ggplot(data = disease) +
  geom_sf(aes(fill = incidence),lwd = .1,color = 'grey') +
  viridis::scale_fill_viridis(option="mako", direction = -1) +
  theme_bw() +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    legend.position = 'inside',
    legend.position.inside = c(.1,.25),
    legend.background = element_rect(fill = 'transparent',color = NA)
  )
f2 = ggplot(data = watershed) +
  geom_sf(aes(fill = watershed),lwd = .1,color = 'grey') +
  tidyterra::scale_fill_whitebox_c() +
  coord_sf(crs = NULL) +
  theme_bw() +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    legend.position = 'inside',
    legend.position.inside = c(.1,.25),
    legend.background = element_rect(fill = 'transparent',color = NA)
  )
f3 = ggplot(data = elevation) +
  geom_sf(aes(fill = elevation),lwd = .1,color = 'grey') +
  tidyterra::scale_fill_hypso_c() +
  theme_bw() +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    legend.position = 'inside',
    legend.position.inside = c(.1,.25),
    legend.background = element_rect(fill = 'transparent',color = NA)
  )
f4 = ggplot(data = soiltype) +
  geom_sf(aes(fill = soiltype),lwd = .1,color = 'grey') +
  tidyterra::scale_fill_wiki_c() +
  theme_bw() +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    legend.position = 'inside',
    legend.position.inside = c(.1,.25),
    legend.background = element_rect(fill = 'transparent',color = NA)
  )

plot_grid(f1,f2,f3,f4, nrow = 2,label_fontfamily = 'serif',
          labels = paste0('(',letters[1:4],')'),
          label_fontface = 'plain',label_size = 10,
          hjust = -1.5,align = 'hv')  -> p
p

Attribute spatial join

NTDs = disease %>%
  st_centroid() %>%
  st_join(watershed[,"watershed"]) %>%
  st_join(elevation[,"elevation"]) %>%
  st_join(soiltype[,"soiltype"])

Check whether has NA in NTDs:

NTDs %>%
  dplyr::filter(if_any(everything(),~is.na(.x)))
## Simple feature collection with 4 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 301567.5 ymin: 3989433 xmax: 318763.3 ymax: 3991906
## Projected CRS: +proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## # A tibble: 4 × 6
##   SP_ID incidence               geom watershed elevation soiltype
## * <chr>     <dbl>        <POINT [m]>     <int>     <int>    <int>
## 1 141        6.48 (318763.3 3991847)         9        NA        3
## 2 165        6.53 (316574.1 3989433)         4        NA        2
## 3 166        6.43 (311439.1 3990674)         4        NA        2
## 4 188        6.26 (301567.5 3991906)        NA         2        3

NTDs %>%
  dplyr::filter(if_all(everything(),~!is.na(.x))) -> NTDs

Factor detector

NTDs = st_drop_geometry(NTDs)
fd = ssh.test(incidence ~ watershed + elevation + soiltype,
              data = NTDs,type = 'factor')
fd
## Spatial Stratified Heterogeneity Test 
##  
##           Factor detector

variable	Q-statistic	P-value
watershed	0.6378	0.0001288
elevation	0.6067	0.04338
soiltype	0.3857	0.3721

Interaction detector

id = ssh.test(incidence ~ watershed + elevation + soiltype,
              data = NTDs,type = 'interaction')
id
## Spatial Stratified Heterogeneity Test 
##  
##          Interaction detector

Interactive variable	Interaction
watershed ∩ elevation	Enhance, bi-
watershed ∩ soiltype	Enhance, bi-
elevation ∩ soiltype	Enhance, bi-

Risk detector

rd = ssh.test(incidence ~ watershed + elevation + soiltype,
              data = NTDs,type = 'risk')
rd
## Spatial Stratified Heterogeneity Test 
##  
##              Risk detector             
## 
## --------------------------------------
## Variable elevation:
## 
## |1   |2   |3   |4   |5   |6   |7   |
## |:---|:---|:---|:---|:---|:---|:---|
## |NA  |Yes |Yes |Yes |Yes |Yes |Yes |
## |Yes |NA  |No  |Yes |Yes |Yes |Yes |
## |Yes |No  |NA  |Yes |Yes |Yes |Yes |
## |Yes |Yes |Yes |NA  |Yes |Yes |Yes |
## |Yes |Yes |Yes |Yes |NA  |Yes |Yes |
## |Yes |Yes |Yes |Yes |Yes |NA  |Yes |
## |Yes |Yes |Yes |Yes |Yes |Yes |NA  |
## --------------------------------------
## Variable soiltype:
## 
## |1   |2   |3   |4   |5   |
## |:---|:---|:---|:---|:---|
## |NA  |Yes |Yes |Yes |No  |
## |Yes |NA  |No  |Yes |Yes |
## |Yes |No  |NA  |Yes |Yes |
## |Yes |Yes |Yes |NA  |Yes |
## |No  |Yes |Yes |Yes |NA  |
## --------------------------------------
## Variable watershed:
## 
## |1   |2   |3   |4   |5   |6   |7   |8   |9   |
## |:---|:---|:---|:---|:---|:---|:---|:---|:---|
## |NA  |Yes |Yes |Yes |Yes |Yes |Yes |Yes |Yes |
## |Yes |NA  |Yes |No  |Yes |Yes |Yes |Yes |Yes |
## |Yes |Yes |NA  |Yes |Yes |Yes |No  |No  |No  |
## |Yes |No  |Yes |NA  |Yes |Yes |Yes |Yes |Yes |
## |Yes |Yes |Yes |Yes |NA  |No  |Yes |Yes |Yes |
## |Yes |Yes |Yes |Yes |No  |NA  |Yes |Yes |Yes |
## |Yes |Yes |No  |Yes |Yes |Yes |NA  |No  |No  |
## |Yes |Yes |No  |Yes |Yes |Yes |No  |NA  |Yes |
## |Yes |Yes |No  |Yes |Yes |Yes |No  |Yes |NA  |

You can change the significant interval by assign alpha argument,the default value of alpha argument is 0.95.

rd99 = ssh.test(incidence ~ watershed + elevation + soiltype,
              data = NTDs,type = 'risk',alpha = 0.99)
rd99
## Spatial Stratified Heterogeneity Test 
##  
##              Risk detector             
## 
## --------------------------------------
## Variable elevation:
## 
## |1   |2   |3   |4   |5   |6   |7   |
## |:---|:---|:---|:---|:---|:---|:---|
## |NA  |Yes |Yes |Yes |Yes |Yes |Yes |
## |Yes |NA  |No  |Yes |Yes |Yes |Yes |
## |Yes |No  |NA  |Yes |Yes |Yes |Yes |
## |Yes |Yes |Yes |NA  |Yes |Yes |Yes |
## |Yes |Yes |Yes |Yes |NA  |Yes |Yes |
## |Yes |Yes |Yes |Yes |Yes |NA  |Yes |
## |Yes |Yes |Yes |Yes |Yes |Yes |NA  |
## --------------------------------------
## Variable soiltype:
## 
## |1   |2   |3   |4   |5   |
## |:---|:---|:---|:---|:---|
## |NA  |Yes |Yes |Yes |No  |
## |Yes |NA  |No  |Yes |Yes |
## |Yes |No  |NA  |Yes |Yes |
## |Yes |Yes |Yes |NA  |Yes |
## |No  |Yes |Yes |Yes |NA  |
## --------------------------------------
## Variable watershed:
## 
## |1   |2   |3   |4   |5   |6   |7   |8   |9   |
## |:---|:---|:---|:---|:---|:---|:---|:---|:---|
## |NA  |Yes |Yes |Yes |Yes |Yes |Yes |Yes |Yes |
## |Yes |NA  |Yes |No  |Yes |Yes |Yes |Yes |Yes |
## |Yes |Yes |NA  |Yes |Yes |Yes |No  |No  |No  |
## |Yes |No  |Yes |NA  |Yes |Yes |Yes |Yes |Yes |
## |Yes |Yes |Yes |Yes |NA  |No  |Yes |Yes |Yes |
## |Yes |Yes |Yes |Yes |No  |NA  |Yes |Yes |Yes |
## |Yes |Yes |No  |Yes |Yes |Yes |NA  |No  |No  |
## |Yes |Yes |No  |Yes |Yes |Yes |No  |NA  |Yes |
## |Yes |Yes |No  |Yes |Yes |Yes |No  |Yes |NA  |

Ecological detector

ed = ssh.test(incidence ~ watershed + elevation + soiltype,
              data = NTDs,type = 'ecological')
ed
## Spatial Stratified Heterogeneity Test 
##  
##           ecological detector          
## 
## --------------------------------------
## 
## 
## |          |elevation |soiltype |
## |:---------|:---------|:--------|
## |watershed |No        |No       |
## |elevation |NA        |No       |

You can also change the significant interval by assign alpha argument,the default value of alpha argument is 0.95.

ed99 = ssh.test(incidence ~ watershed + elevation + soiltype,
              data = NTDs,type = 'ecological',alpha = 0.99)
ed99
## Spatial Stratified Heterogeneity Test 
##  
##           ecological detector          
## 
## --------------------------------------
## 
## 
## |          |elevation |soiltype |
## |:---------|:---------|:--------|
## |watershed |No        |No       |
## |elevation |NA        |No       |

Wenbo Lv

Load package and pre-processing data.

Factor detector

Interaction detector

Risk detector

Ecological detector