---
bibliography: references.bib
---

# Basic statistics for spatial analysis

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

## 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.

```{r load_cases, eval = TRUE, echo = TRUE, nm = TRUE, fig.width=8, class.output="code-out", warning=FALSE, message=FALSE}
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.

```{r cases_visualization, eval = TRUE, echo = TRUE, nm = TRUE, fig.width=8, class.output="code-out", warning=FALSE, message=FALSE}

# View the cases object
head(cases)

# 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 aggregation since they make available data on population count and structures. In this study, we will use district as the areal unit of the study.

```{r district_aggregate, eval = TRUE, echo = TRUE, nm = TRUE, fig.width=8, class.output="code-out", warning=FALSE, message=FALSE}
# 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.

```{r indicators, eval = TRUE, echo = TRUE, nm = TRUE, fig.width=8, fig.height=4, class.output="code-out", warning=FALSE, message=FALSE}

# 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
```

```{r inc_visualization, eval = TRUE, echo = TRUE, nm = TRUE, fig.width=8, fig.height=4, class.output="code-out", warning=FALSE, message=FALSE}
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")
```

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. However, assumption does not hold for all diseases and for all observed events since confounding effects can create nuisance into the interpretations (e.g. the number of childhood illness and death outcomes in a district 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).

## Cluster analysis

Since this W fever seems to have a heterogeneous distribution across Cambodia, it would be interesting to study where excess of cases appears, i.e. to identify clusters of the disease. The first question is to wonder if data are auto correlated or spatially independent, i.e. study if neighboring districts are likely to have similar incidence.

In statistics, problems are 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.

### Spatial autocorrelation (Moran's I test)

A popular test for spatial autocorrelation is the Moran's test. This test tells us whether nearby units tend to exhibit similar incidences. It ranges from -1 to +1. A value of -1 denote 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.

Here the statistics hypothesis are :

-   H0 : the distribution of cases is spatially independant, i.e. Moran's I value is 0.

-   H1: the distribution of cases is spatially autocorrelated, i.e. Moran's I value is different than 0.

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.

```{r MoransI, eval = TRUE, echo = TRUE, nm = TRUE, fig.width=8, class.output="code-out", warning=FALSE, message=FALSE}

library(spdep) # Functions for creating spatial weight, spatial analysis
library(DCluster)  # Package with functions for spatial cluster analysis)

qnb <- poly2nb(district)
q_listw <- nb2listw(qnb, style = 'W') # row-standardized weights

# Moran's I test
moranI.test(cases ~ offset(log(expected)), 
                  data = district,
                  model = 'poisson',
                  R = 499,
                  listw = q_listw,
                  n = 159,
                  S0 = Szero(q_listw))

```

### Spatial scan statistics

While Moran's indice focuses on finding correlation between neighboring polygons, the spatial scan statistic compare the incidence level of a given windows of observation with the incidence level outside of this windows.

The package `SpatialEpi`

### 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.

### Expectation-based cluster

In many case, population is not specific enough to

### To go further ...
