---
title: "5. Reducing sampling bias"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{5. Reducing sampling bias}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,
  fig.width = 7,
  fig.height = 4,
  fig.align='center'
)
```

## Introduction

For most studies relying on primary biodiversity data, sampling bias can be a problem because it creates clusters of points that reflect accessible areas and the preferences of the researchers who collected the data, rather than the ecological preferences of the species.

Two common approaches to mitigate sampling bias involve thinning occurrence records (i.e., removing records that are close to each other). Thinning can be performed in either geographic or environmental space. The `RuHere` package allows the use of both approaches and also enables combining the results of both methods.

As an example, let’s use the records of *Araucaria angustifolia* available in the package, after removing records flagged as potentially problematic:

```{r import data}
# Load packages
library(RuHere)
library(terra)
library(mapview)

# Import occurrence data
data("occ_flagged", package = "RuHere")
# Remove flagged records
occ <- remove_flagged(occ = occ_flagged)
# Plot records
ggmap_here(occ = occ)
```
<br>

## Heatmap for occurrence data

Before starting the thinning process, let’s create a heatmap based on a kernel density estimation of the records:

```{r}
# Generate heatmap
heatmap <- spatial_kde(occ = occ, resolution = 0.2, buffer_extent = 50,
                       radius = 2, zero_as_NA = TRUE)
```
<br>

We can use the function `ggmap_here()` to plot the occurrences and the heatmap (or use `map_here()` if you prefer an interactive version).

```{r}
ggmap_here(occ = occ, size_points = 0.5, heatmap = heatmap)
```
<br>

We can observe a notable hotspot with a cluster of records around Curitiba city (capital of Paraná, southern Brazil), which indicates strong sampling bias. Let’s evaluate the effect of thinning these records in geographic space.

## Thinning in geographic space

The `thin_geo()` function flags occurrence records for thinning by keeping only one record per species within a radius of `d` kilometers. This function is similar to `thin()` from the [spThin package](https://nsojournals.onlinelibrary.wiley.com/doi/10.1111/ecog.01132), but with an important difference: it allows specifying a priority order for retaining records.

When a thinning distance is provided (e.g., 10 km), the function identifies clusters of records within that distance. Within each cluster, it retains the record with the highest priority according to the column defined in prioritary_column (for example, using the most recent record when `prioritary_column = "year"`), and flags all remaining nearby records for removal. If `prioritary_column` is `NULL`, the priority follows the original row order of the input occ data frame.

Let’s thin the records using a 10 km radius and keep the most recent record as the priority:

```{r}
# Thin records using a 10 km distance threshold
occ_thin <- thin_geo(occ = occ, d = 10, prioritary_column = "year")
sum(!occ_thin$thin_geo_flag)  # Number of records flagged for removal
```
<br>

The function flagged 1,860 records for removal. Let’s visualize the flagged records using `map_here()` and create a heatmap with the remaining records.

```{r}
# Remove flagged records
occ_thin_geo <- remove_flagged(occ = occ_thin)
# Create heatmap
heatmap_thin_geo <- spatial_kde(occ = occ_thin_geo, resolution = 0.2, 
                                buffer_extent = 50, radius = 2, 
                                zero_as_NA = TRUE)
# Plot
ggmap_here(occ_thin, size_points = 0.5, heatmap = heatmap_thin_geo)
```
<br>

The thinned dataset (in green) produces a more spatially uniform heatmap and reduces the strong sampling bias around Curitiba. We explore different distance thresholds in the next section.

<br>

## Selecting the best distance to thin records

A key question is determining the optimal thinning distance, since we rarely have sufficient biological justification for choosing a specific value. To address this issue, we adapted the approach of [Velazco et al. (2020)](https://onlinelibrary.wiley.com/doi/10.1111/ddi.13215), which computes spatial autocorrelation (Moran’s I) for datasets generated using different thinning distances and selects the distance that yields the lowest average spatial autocorrelation.

We extended this procedure by adjusting the selection rules to avoid choosing datasets with too few records or unrealistically low Moran’s I values. See `help(flag_geo_moran)` for full details of the selection procedure.

As an example, let’s test the effect of thinning using distances of 1, 3, 5, 7, 10, 15, 20, and 30 km. For computing spatial correlation, we need a raster containing the environmental variables. Here, we again specify a priority for retaining records (the most recent ones).

```{r}
# Load example of raster variables
data("worldclim", package = "RuHere")
# Unwrap Packed raster
r <- terra::unwrap(worldclim)
# Select thinned occurrences
occ_geo_moran <- flag_geo_moran(occ = occ, 
                                d = c(1, 3, 5, 7, 10, 15, 20, 30), 
                                prioritary_column = "year",
                                env_layers = r)
```
<br>

The results for each tested distance are returned in the `imoran` data.frame. It includes Moran’s I for each variable and summary statistics across variables (mean, median, minimum, and maximum), along with the number of retained records (`n_filtered`) and the proportion of records flagged (`prop_lost`).

```{r}
occ_geo_moran$imoran
```
<br>

The "best" distance that reduces the spatial autocorrelation without discarding too many records was 15km. Using this threshold, 2,040 records were flagged. 

```{r}
# Best distance selected
occ_geo_moran$distance
# Number of flagged records using this distance to thin
sum(!occ_geo_moran$occ$thin_geo_flag)
```
<br>

Visual inspection shows an even more uniform heatmap:

```{r}
# Remove flagged records
occ_thin_geo_moran <- remove_flagged(occ = occ_geo_moran$occ)
# Create heatmap
heatmap_thin_geo_moran <- spatial_kde(occ = occ_thin_geo_moran, 
                                      resolution = 0.2, 
                                      buffer_extent = 50, radius = 2,
                                      zero_as_NA = TRUE)
ggmap_here(occ = occ_geo_moran$occ, size_points = 0.5, 
           heatmap = heatmap_thin_geo_moran)
```
<br>

A potential issue when filtering records in geographic space is that two nearby records may actually occur in distinct environmental conditions, especially in highly heterogeneous regions. This can lead to the loss of unique information about the species’ niche and environmental tolerances. To address this, we can instead apply thinning in environmental space, as explored in the next sections.
<br>

## Thinning in environmental space

Thinning in environmental space removes records with similar environmental conditions, representing redundant ecological information. This is achieved by building a multidimensional grid in environmental space, dividing each variable into `n_bins` equally sized intervals. Each record is assigned to a unique environmental block (a combination of bins), and records within the same block (i.e., environmentally similar) are flagged for removal.

To illustrate how the environmental grid works, let’s use `get_env_bins()` with 10 bins:

```{r}
# Get bins
b <- get_env_bins(occ = occ, env_layers = r, n_bins = 10)
head(b$data)
```
<br>

The function returns the environmental block IDs for each record. We can visualize the grid for any two variables:

```{r}
# Plot
plot_env_bins(b, x_var = "bio_1", y_var = "bio_12",
              xlab = "Temperature", ylab = "Precipitation")
```
<br>

We can see that several records fall into the same block, meaning they are environmentally similar and therefore redundant. Let’s flag these redundant records using `thin_env()`:

```{r}
# Flag records that are close to each other in the enviromnetal space
occ_thin_env <- thin_env(occ = occ, env_layers = r, n_bins = 10, 
                         prioritary_column = "year")
# Number of flagged (redundant) records
sum(!occ_thin_env$thin_env_flag) #Number of flagged records
```
<br>

The function flagged 2,227 records. Let's visualize these and create a heatmap of the remaining data.

```{r}
# Remove flagged records
occ_thinned_env <- remove_flagged(occ = occ_thin_env)
# Create heatmap
heatmap_thin_env <- spatial_kde(occ = occ_thinned_env, resolution = 0.2, 
                                buffer_extent = 50, radius = 2, 
                                zero_as_NA = TRUE)
ggmap_here(occ_thin_env, size_points = 0.5, heatmap = heatmap_thin_env)
```
<br>

Thinning in environmental space produces a spatial pattern that differs from the dataset filtered exclusively in geographic space.

Similar to the thinning process in geographic space, we can test different numbers of bins (see next section).

## Selecting the best number of environmental bins

In geographic thinning, the key parameter is distance. In environmental thinning, it is the number of bins. More bins result in finer partitions of environmental space, reducing the chances of records falling into the same block. As with geographic thinning, we can test multiple bin values and select the one that reduces spatial autocorrelation without discarding many records.

Here, we test 5, 10, 20, 30, 40, 50, 60, 70, and 80 bins:

```{r}
# Select thinned occurrences
occ_env_moran <- flag_env_moran(occ = occ, 
                                n_bins = c(5, 10, 20, 30, 40, 50, 60, 70, 80), 
                                prioritary_column = "year",
                                env_layers = r)
```

<br>

Again, results are returned in an imoran data frame. It includes Moran’s I for each variable and summary statistics across variables (mean, median, minimum, and maximum), along with the number of retained records (`n_filtered`) and the proportion of records flagged (`prop_lost`).

```{r}
occ_env_moran$imoran
```
<br>

The "best" number of bins that reduces the spatial autocorrelation without discarding too many records was 70. With this threshold, 1,659 records were flagged.

```{r}
# Best distance selected
occ_env_moran$n_bins
# Number of flagged records using this distance to thin
sum(!occ_env_moran$occ$thin_env_flag)
```
<br>

Let's check the distribution of these records and the heatmap generated with the unflagged records:

```{r}
# Remove flagged records
occ_thin_env_moran <- remove_flagged(occ = occ_env_moran$occ)
# Create heatmap
heatmap_thin_env_moran <- spatial_kde(occ = occ_thin_env_moran, 
                                      resolution = 0.2, 
                                      buffer_extent = 50, radius = 2,
                                      zero_as_NA = TRUE)
ggmap_here(occ = occ_env_moran$occ, size_points = 0.5, 
           heatmap = heatmap_thin_env_moran)
```
<br>

## Consensus between environmental and geographic thinning

Thinning in environmental space can suffer from the opposite issue of geographic thinning: environmentally similar records may be geographically far apart, potentially removing important information about the species’ geographic range.

To address this, the `flag_consensus()` function can be used to flag records only when they are redundant in both geographic and environmental space.

```{r}
# Flag occurrences by thinning in geographic space
occ_geo_moran <- flag_geo_moran(occ = occ, 
                                d = c(1, 3, 5, 7, 10, 15, 20, 30), 
                                prioritary_column = "year",
                                env_layers = r)

# Flag occurrences by thinning in environmental space
occ_env_moran <- flag_env_moran(occ = occ_geo_moran$occ, 
                                n_bins = c(5, 10, 20, 30, 40, 50, 60, 70, 80), 
                                prioritary_column = "year",
                                env_layers = r)
# Get consensus
occ_consensus <- flag_consensus(occ = occ_env_moran$occ,
                                flags = c("thin_geo", "thin_env"),
                                consensus_rule = "any_true",
                                flag_name = "thin_geo_env_flag")
# Remove flagged
occ_consensus_filtered <- remove_flagged(occ = occ_consensus, flags = NULL,
                                         additional_flags = c("thin_geo_env" = "thin_geo_env_flag"))

# Create heatmap
heatmap_consensus_filtered <- spatial_kde(occ = occ_consensus_filtered,
                                          resolution = 0.2,
                                          buffer_extent = 50,
                                          radius = 2, zero_as_NA = TRUE)

```

<br>

Let’s visualize which records were flagged by geographic thinning, environmental thinning, or both:

```{r}
ggmap_here(occ = occ_consensus, 
          flags = c("thin_geo", "thin_env"), 
          additional_flags = "thin_geo_env_flag", 
          names_additional_flags = "Thinned geo-env", 
          col_additional_flags = "blue",
          size_points = 0.5,
          heatmap = heatmap_consensus_filtered)
```
<br>

Red points indicate records thinned in geographic space; yellow points indicate thinning in environmental space; and blue points indicate records thinned in both.

Note that we retained the records flagged only in the geographic thinning, those flagged only in the environmental thinning, and those not flagged by either method.
