library(sf)
#Import Cambodia country border
= st_read("data_cambodia/cambodia.gpkg", layer = "country", quiet = TRUE)
country #Import provincial administrative border of Cambodia
= st_read("data_cambodia/cambodia.gpkg", layer = "education", quiet = TRUE)
education #Import district administrative border of Cambodia
= st_read("data_cambodia/cambodia.gpkg", layer = "district", quiet = TRUE)
district
# Import locations of cases from an imaginary disease
= st_read("data_cambodia/cambodia.gpkg", layer = "cases", quiet = TRUE)
cases = subset(cases, Disease == "W fever") cases
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.
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
$cases <- lengths(st_intersects(district, cases)) district
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)
$incidence = district$cases/district$T_POP * 100000
district
# Compute the global risk
= sum(district$cases)/sum(district$T_POP)
rate
# Compute expected number of cases
$expected = district$T_POP * rate
district
# Compute SIR
$SIR = district$cases / district$expected district
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
= c(0, exp(mf_get_breaks(log(district$SIR), nbreaks = 8, breaks = "pretty")))
break_SIR = c("#273871", "#3267AD", "#6496C8", "#9BBFDD", "#CDE3F0", "#FFCEBC", "#FF967E", "#F64D41", "#B90E36")
col_pal
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