7  Basic statistics for spatial analysis

This section aims at providing some basic statistical tools to study the spatial distribution of epidemiological data.

7.1 Import and visualize epidemiological data

In this section, we load data that reference the cases of an imaginary disease throughout Cambodia. Each point correspond to the geolocalisation of a case.

library(sf)

#Import Cambodia country border
country = st_read("data_cambodia/cambodia.gpkg", layer = "country", quiet = TRUE)
#Import provincial administrative border of Cambodia
education = st_read("data_cambodia/cambodia.gpkg", layer = "education", quiet = TRUE)
#Import district administrative border of Cambodia
district = st_read("data_cambodia/cambodia.gpkg", layer = "district", quiet = TRUE)

# Import locations of cases from an imaginary disease
cases = st_read("data_cambodia/cambodia.gpkg", layer = "cases", quiet = TRUE)
cases = subset(cases, Disease == "W fever")

The first step of any statistical analysis always consists on visualizing the data to check they were correctly loaded and to observe general pattern of the cases.

# View the cases object
head(cases)
Simple feature collection with 6 features and 2 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: 255891 ymin: 1179092 xmax: 506647.4 ymax: 1467441
Projected CRS: WGS 84 / UTM zone 48N
  id Disease                           geom
1  0 W fever MULTIPOINT ((280036.2 12841...
2  1 W fever MULTIPOINT ((451859.5 11790...
3  2 W fever  MULTIPOINT ((255891 1467441))
4  5 W fever MULTIPOINT ((506647.4 12322...
5  6 W fever  MULTIPOINT ((440668 1197958))
6  7 W fever MULTIPOINT ((481594.5 12714...
# Map the cases
library(mapsf)

mf_map(x = district, border = "white")
mf_map(x = country,lwd = 2, col = NA, add = TRUE)
mf_map(x = cases, lwd = .5, col = "#990000", pch = 20, add = TRUE)

In epidemiology, the true meaning of point is very questionable. If it usually gives the location of an observation, its not clear if this observation represents an event of interest (e.g. illness, death, …) or a person at risk (e.g. a participant that may or may not experience the disease). Considering a ratio of event compared to a population at risk is often more informative than just considering cases. Administrative divisions of countries appears as great areal units for cases aggreagation since they make available data on population count and structures. In this study, we will use district as the areal unit of the study.

# Aggregate cases over districts
district$cases <- lengths(st_intersects(district, cases))

The incidence (\(\frac{cases}{population}\)) is commonly use to represent cases distribution related to population density but other indicators exists. As example, the standardized incidence ratios (SIRs) represents the deviation of observed and expected number of cases and is expressed as \(SIR = \frac{Y_i}{E_i}\) with \(Y_i\), the observed number of cases and \(E_i\), the expected number of cases. In this study, we computed the expected number of cases in each district by assuming infections are homogeneously distributed across Cambodia, i.e. the incidence is the same in each district.

# Compute incidence in each district (per 100 000 population)
district$incidence = district$cases/district$T_POP * 100000

# Compute the global risk
rate = sum(district$cases)/sum(district$T_POP)

# Compute expected number of cases 
district$expected = district$T_POP * rate

# Compute SIR
district$SIR = district$cases / district$expected
par(mfrow = c(1, 3))
# Plot number of cases using proportional symbol 
mf_map(x = district) 
mf_map(
  x = district, 
  var = "cases",
  val_max = 50,
  type = "prop",
  col = "#990000", 
  leg_title = "Cases")
mf_layout(title = "Number of cases of W Fever")

# Plot incidence 
mf_map(x = district,
       var = "incidence",
       type = "choro",
       pal = "Reds 3",
       leg_title = "Incidence \n(per 100 000)")
mf_layout(title = "Incidence of W Fever")

# Plot SIRs

# create breaks and associated color palette
break_SIR = c(0, exp(mf_get_breaks(log(district$SIR), nbreaks = 8, breaks = "pretty")))
col_pal = c("#273871", "#3267AD", "#6496C8", "#9BBFDD", "#CDE3F0", "#FFCEBC", "#FF967E", "#F64D41", "#B90E36")

mf_map(x = district,
       var = "SIR",
       type = "choro",
       breaks = break_SIR, 
       pal = col_pal, 
       cex = 2,
       leg_title = "SIR")
mf_layout(title = "Standardized Incidence Ratio of W Fever in Cambodia")

These maps illustrates the spatial heterogenity of the cases. The incidence shows how the disease vary from one district to another while the SIR highlight districts that have :

  • higher risk than average (SIR > 1) when standardized for population

  • lower risk than average (SIR < 1) when standardized for population

  • average risk (SIR ~ 1) when standardized for population

In this example, we standardized the cases distribution for population count. This simple standardization assume that the risk of contracting the disease is similar for each person. Howerver, this case does not apply for all disease and for all observed events (e.g. the number of childhood illness and death outcomes are usually related to the age pyramid) and you should keep in mind that other standardization can be performed based on variables known to have an effect but that you don’t want to analyze (e.g. sex ratio, occupations, age pyramid).

7.2 Basics statistics

The problem is usually expressed by defining two hypothesis : the null hypothesis (H0), i.e. an a priori hypothesis of the studied phenomenon (e.g. the situation is a random) and the alternative hypothesis (HA), e.g. the situation is not random. The main principle is to measure how likely the observed situation belong to the ensemble of situation that are possible under the H0 hypothesis.

The statistical analysis performed relies on the type of data.

7.2.1 Spatial autocorrelation (Moran’s I test)

A popular test for spatial autocorrelation is the Moran’s test.

Moran’s I test tells us whether nearby units tend to exhibit similar rates. It ranges from -1 to +1, whith a value of -1 denoting that units with low rates are located near other units with high rates, while a Moran’s I value of +1 indicates a concentration of spatial units exhibiting similar rates.

We will compute the Moran’s statistics using spdep and Dcluster packages. spdep package provides a collection of functions to analyze spatial correlations of polygons and works with sp objects. Dcluster package provides a set of functions for the detection of spatial clusters of disease using count data.

# Plot the incidence histogramm
hist(log(district$incidence))

7.3 Cluster analysis

In epidemiology, the definition of a cluster

7.3.1 Population-based clusters (kulldorf statistic)

Kulldorff ’s spatial scan statistic identifies the most likely disease clusters maximizing the likelihood that disease cases are located within a set of concentric circles that are moved across the study area.

7.3.2 Expectation-based cluster

In many case, population is not specific enough to

7.3.3 To go further …