---
title: "Spatial Analysis: Endemism, Fragmentation, and SAR"
author: "Gilles Colling"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Spatial Analysis: Endemism, Fragmentation, and SAR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  dev = "svglite",
  fig.ext = "svg"
)
```

```{r transparent, include = FALSE}
transparent <- if (requireNamespace("ggplot2", quietly = TRUE)) {
  ggplot2::theme(
    panel.background = ggplot2::element_rect(fill = "transparent", colour = NA),
    plot.background  = ggplot2::element_rect(fill = "transparent", colour = NA)
  )
} else {
  NULL
}
```

## Introduction

Conservation planning rarely asks how many species live in a region. It asks
which fragments of that region carry species found nowhere else, how much
richness is lost when continuous habitat is broken into patches, and how
survey effort distorts the apparent richness of poorly sampled areas. These
questions sit beside the standard species accumulation curve, and spacc
collects the tools that answer them under one set of objects.

This vignette covers the model family aimed at spatial structure and
conservation value: endemism-area curves, the species-fragmented area
relationship (SFAR), the sampling-effort corrected species-area relationship
(SESARS), per-site accumulation metrics for spatial prioritisation, Moran
eigenvector maps (MEM) for decomposing spatial autocorrelation, the wavefront
expansion diagnostic, and grid-based spatial subsampling. Each function takes
either a site-by-species matrix with coordinates or a fitted `spacc` object,
and each returns a typed S3 object with `print()`, `summary()`, and `plot()`
methods.

The worked examples use a single simulated study region with known ground truth.
Range size varies across species on purpose: a subset of narrow-ranged
endemics share the region with widespread generalists. That contrast drives
every result below, from the shape of the endemism curve to the sign of the
fragmentation exponent. Working from a simulation with controlled inputs makes
it possible to read each output against the structure that produced it, rather
than against an unknown truth.

The seven functions divide into three groups by what they consume. Endemism
curves, per-site metrics, and the wavefront diagnostic take the raw
site-by-species matrix and coordinates, and run their own accumulation
internally. SFAR and SESARS take a fitted `spacc` object and a covariate (patch
labels or effort), since they reshape an existing curve rather than building a
new one. The MEM tools take coordinates to build a spatial basis, then accept
any per-site response for partitioning. Knowing which group a function belongs
to is the quickest way to see what it needs and what it returns.

## Models and theory

### Endemism-area concept

The species-area relationship counts how richness grows with area. The
endemism-area relationship asks a sharper question: of the species present in
an accumulated set of sites, how many occur *only* within that set. Let
$S(A)$ be the richness accumulated over area $A$, and let a species be endemic
to the accumulated region when every one of its occurrences falls inside it.
At each accumulation step $k$ over the visit order, a species with cumulative
occurrence count $a_k$ and global occurrence count $t$ is endemic when
$a_k = t$ and $a_k > 0$. The endemic count $E(A)$ rises as narrow-ranged
species are captured and approaches total richness only when every site is
included, since at full extent every species is trivially contained.

The diagnostic value lies at intermediate area. A site set that holds many
endemic species is irreplaceable: those species disappear from the protected
network if the set is dropped. Endemism concentrated in small areas signals a
setting where reserve placement matters more than reserve size.

The shape of $E(A)$ carries information beyond its endpoints. A curve that
climbs slowly and only catches up to richness near full extent describes a pool
of widespread species with diffuse ranges, where almost no area is irreplaceable
until nearly everything is included. A curve that tracks richness closely from
the start describes a pool of narrow endemics, where even small areas already
hold their full complement of locally restricted species. The endemism
proportion $E(A)/S(A)$ makes this readable on a single axis: it is the fraction
of the accumulated species pool that is local to the area sampled, and its early
behaviour separates landscapes built from widespread species from those built
from endemics.

### Species-fragmented area relationship (SFAR)

Habitat loss and habitat fragmentation are distinct processes. Removing area
reduces richness through the classic power law. Splitting the remaining area
into more, smaller pieces reduces it further, through edge effects, isolation,
and minimum viable area. The SFAR (Hanski et al. 2013) separates the two with
an additive log-linear model,

$$\log(S) = c + z \cdot \log(A) + f \cdot \log(n),$$

where $S$ is species richness, $A$ is total habitat area, $n$ is the number of
fragments, $z$ is the ordinary area exponent of the SAR, $f$ is the
fragmentation exponent, and $c$ is the intercept on the log scale (so $e^{c}$
is the richness expected at unit area and a single fragment). A negative $f$
means that, holding area fixed, more fragments support fewer species. spacc
reports $f$ with the sign convention that a positive value indicates
fragmentation reduces richness beyond the area effect, so the model is written
$S = e^{c} \cdot A^{z} \cdot n^{-f}$ in the printed output.

### Sampling-effort corrected SAR (SESARS)

Atlas and citizen-science data confound area with effort. A large region that
is sparsely surveyed can look species-poor for the wrong reason. SESARS
(Dennstadt et al. 2019) adds an effort term to the SAR,

$$\log(S) = c + z \cdot \log(A) + e \cdot \log(E),$$

where $E$ is sampling effort (visits, hours, trap-nights), $e$ is the effort
exponent, and $A$, $z$, $c$, $S$ carry their SAR meanings. A positive $e$
quantifies how much of the apparent richness gradient tracks survey intensity
rather than area. With effort in the model, the area exponent $z$ estimates the
area effect net of unequal sampling.

### Moran eigenvector maps

Spatial autocorrelation is structure worth describing in its own right. Moran
eigenvector maps build an orthogonal basis of spatial patterns directly from
coordinates. The construction truncates the pairwise
distance matrix at a threshold (by default the largest edge of the minimum
spanning tree, the smallest distance that keeps all sites connected),
double-centres the squared truncated matrix, and extracts its eigenvectors.
Each eigenvector is a map: broad-scale patterns carry large eigenvalues and
strong positive Moran's I, fine-scale patterns carry small eigenvalues. Moran's
I for eigenvector $v$ with connectivity matrix $W$ is

$$I = \frac{n}{\sum_{ij} w_{ij}}
      \frac{\sum_{ij} w_{ij}\,(v_i - \bar v)(v_j - \bar v)}
           {\sum_i (v_i - \bar v)^2},$$

with $n$ sites and weights $w_{ij} = 1$ when sites $i$ and $j$ lie within the
threshold. Regressing a per-site response on a selected set of these
eigenvectors partitions its variance into a spatial component and a residual,
and reveals at which scales the structure lives.

Two methods share the construction. PCNM retains every eigenvector with a
positive eigenvalue, including those describing negative autocorrelation
(checkerboard-like patterns where neighbours differ). dbMEM keeps only the
eigenvectors with positive Moran's I, the patterns where neighbours resemble
each other, which is the usual target for modelling ecological gradients. The
threshold matters as much as the method: a small threshold connects only the
nearest sites and yields many fine-scale vectors, while a large threshold
connects distant sites and emphasises broad gradients. The minimum-spanning-tree
default is the smallest threshold that leaves no site isolated, a principled
starting point that adapts to irregular sampling.

## Simulating range-size variation

The simulated region holds 80 sites scattered over a 100 by 100 unit square and
40 species. Each species occupies a Gaussian neighbourhood around a random
centre, so occurrence probability falls off with squared distance from that
centre. The decay constant `spread` sets range size: a large value confines a
species to a tight cluster near its centre, a small value lets it spill across
the whole extent.

```{r simulate}
library(spacc)

set.seed(42)
n_sites <- 80
n_species <- 40

coords <- data.frame(
  x = runif(n_sites, 0, 100),
  y = runif(n_sites, 0, 100)
)

# Species with varying range sizes (some endemic, some widespread)
species <- matrix(0L, n_sites, n_species)
for (sp in seq_len(n_species)) {
  cx <- runif(1, 10, 90)
  cy <- runif(1, 10, 90)
  # First 10 species are narrow-ranged (endemics)
  spread <- if (sp <= 10) 0.005 else 0.001
  prob <- exp(-spread * ((coords$x - cx)^2 + (coords$y - cy)^2))
  species[, sp] <- rbinom(n_sites, 1, prob)
}
colnames(species) <- paste0("sp", seq_len(n_species))
```

The first 10 species use `spread = 0.005`, which gives a Gaussian half-width of
roughly 12 units, about an eighth of the map edge. These are the
endemics: tight clusters that vanish if their corner of the map is dropped.
The remaining 30 use `spread = 0.001`, a half-width near 26 units, so they
appear across much of the area. The split is the ground truth that every
downstream analysis recovers.

```{r range-sizes}
range_size <- colSums(species)
tapply(range_size, c(rep("endemic", 10), rep("widespread", 30)), mean)
```

Mean occupancy is lower for the narrow-ranged group, as designed. The contrast
is moderate rather than extreme, which keeps the simulated curves realistic:
endemics are not single-site rarities, and widespread species do not saturate
every site.

## Endemism-area curves

`spaccEndemism()` tracks endemic richness as the sampled area expands. It runs
a nearest-neighbour accumulation from each of `n_seeds` random starts, and at
each step counts species whose every occurrence falls inside the accumulated
sites. The print method reports the final mean endemism, which equals total
richness once all sites are in.

```{r endemism}
end <- spaccEndemism(species, coords, n_seeds = 20, progress = FALSE, seed = 1)
end
```

The `summary()` method returns the curve as a data frame, with mean richness,
mean endemism, a percentile confidence band, and the endemism proportion at
each accumulated extent. The proportion is the share of accumulated species
that are local to the area sampled so far.

```{r endemism-summary}
es <- summary(end)
head(es[, c("sites", "mean_richness", "mean_endemism", "endemism_proportion")], 4)
```

```{r plot-endemism, eval = requireNamespace("ggplot2", quietly = TRUE), fig.cap = "Total richness and endemic richness as area expands."}
plot(end) + transparent
```

At small extents the endemism proportion is high: most species seen so far are
confined to the few sites visited, because a handful of sites can only contain
narrow-ranged occupants. As more sites enter, widespread species accumulate
and the proportion falls, then both curves converge as the remaining endemics
are absorbed. The early peak in proportion marks the spatial scale at which
irreplaceability is concentrated. Sites contributing to that peak are the ones
a reserve network cannot afford to lose, because the species they hold occur
nowhere else in the dataset.

Setting `map = TRUE` runs the accumulation from every site and records the
endemic count reachable from each starting point, which feeds the map plot and
the `as_sf()` export covered alongside the per-site metrics below.

## Species-fragmented area relationship (SFAR)

SFAR needs two inputs: a fitted `spacc` object and a vector assigning each site
to a habitat fragment. Patches here come from k-means clustering of the
coordinates, a stand-in for mapped habitat patches. `kmeans()` is base R, so no
guard is required.

```{r sfar}
sac <- spacc(species, coords, n_seeds = 20, progress = FALSE, seed = 1)

set.seed(123)
patches <- kmeans(coords, centers = 5)$cluster

sfar_fit <- sfar(sac, patches)
sfar_fit
```

The fitted coefficients are stored in `$coef` as a named vector: `log_c` is the
log-scale intercept, `z` the area exponent, and `f` the fragmentation exponent
with the sign convention that positive means fragmentation costs richness.

```{r sfar-coef}
round(sfar_fit$coef, 3)
```

```{r plot-sfar, eval = requireNamespace("ggplot2", quietly = TRUE), fig.cap = "SFAR fit: observed richness and the area-plus-fragmentation prediction."}
plot(sfar_fit) + transparent
```

The area exponent $z$ takes a value in the usual SAR range, confirming that
richness climbs with the number of accumulated sites. The fragmentation
exponent $f$ measures the residual penalty once area is accounted for. A
negative printed $f$ (reported as "no additional fragmentation penalty") means
the patch count adds nothing beyond area in this dataset, which is the
expected result for a simulation where patches are an arbitrary partition of a
continuous surface rather than isolated habitat islands. In a genuinely
fragmented region, where each patch loses edge-sensitive and area-demanding species, $f$
turns positive and the SFAR line falls below the pure-area expectation. The
interpretation hinges on that gap: SFAR is read as the distance between the
fitted curve and what area alone would predict.

The mechanism behind a positive $f$ is the reason the term is added at all.
Theory predicts that fragmentation matters most when total habitat is already
low, because at that point the remaining area is divided into pieces too small
to sustain populations that need large territories or interior conditions away
from edges. A single block of a given area then holds more species than the same
area split into many fragments. The SFAR captures this as a multiplicative
penalty $n^{-f}$ that grows with fragment count, so doubling the number of
fragments at fixed area reduces predicted richness by a factor $2^{-f}$.
Estimating $f$ therefore requires patches that are genuinely isolated and that
vary in number across the dataset, conditions the practical guidance below makes
explicit.

## Sampling-effort corrected SAR (SESARS)

SESARS corrects the SAR for uneven survey intensity. Effort enters as a
per-site vector with one entry per site. Here effort is drawn from a Poisson
distribution, mimicking a survey where some sites received many visits and
others few.

```{r sesars}
set.seed(7)
effort <- rpois(n_sites, 10) + 1

sesars_fit <- sesars(sac, effort)
sesars_fit
```

The coefficients live in `$coef`: `log_c` (or `c` under the additive model), the
area exponent `z`, and the effort exponent (named `w` in the object, the $e$ of
the equation above). The R-squared in `$r_squared` reports how much of the
richness variation the joint area-and-effort model captures.

```{r sesars-coef}
round(sesars_fit$coef, 3)
sesars_fit$r_squared
```

```{r plot-sesars, eval = requireNamespace("ggplot2", quietly = TRUE), fig.cap = "SESARS fit with the area-plus-effort prediction overlaid on observed richness."}
plot(sesars_fit) + transparent
```

A positive effort exponent means apparent richness rises with survey intensity
independent of area, the bias SESARS exists to expose. When effort and area are
correlated in field data, ignoring effort inflates the area exponent, because
the SAR absorbs the effort signal into $z$. Fitting both terms returns an area
exponent net of effort, which is the quantity comparable across regions
surveyed at different intensities.

The effort enters as a cumulative quantity aligned with the accumulation order,
so the model relates accumulated richness to accumulated area and accumulated
effort step by step. In this simulation effort is independent of the species
data, so its exponent reflects only the Poisson draw rather than a real
survey-intensity gradient. The instructive case is the field dataset where
under-surveyed regions are also species-poor in the raw counts: there the naive
SAR reads the low richness of sparsely visited sites as an area effect, and only
the effort term separates true poverty from missing detections. SESARS is most
useful precisely when effort and richness move together, since that is the
configuration in which the uncorrected SAR is most misleading.

## Per-site metrics for spatial prioritisation

`spaccMetrics()` runs an accumulation curve starting from each site and reduces
each curve to summary statistics, giving one row per site. Three metrics are
requested here: `slope_10` (the initial slope over the first 10 percent of
sites), `half_richness` (the number of sites needed to reach half the total
species pool), and `auc` (the area under the accumulation curve, normalised by
length).

```{r metrics}
met <- spaccMetrics(species, coords,
                    metrics = c("slope_10", "half_richness", "auc"),
                    progress = FALSE)
summary(met)
```

The per-site values are stored in `met$metrics`, a data frame carrying the
coordinates and a column for each metric, which makes it directly mappable and
joinable.

```{r metrics-df}
head(met$metrics[, c("site_id", "x", "y", "slope_10", "half_richness", "auc")], 4)
```

```{r plot-metrics, eval = requireNamespace("ggplot2", quietly = TRUE), fig.cap = "Spatial heatmap of the initial accumulation slope per site."}
plot(met, metric = "slope_10", type = "heatmap") + transparent
```

A steep initial slope marks a site whose neighbourhood adds species fast: many
distinct species are packed nearby, the spatial signature of a diversity
hotspot. A low `half_richness` flags a site from which half the species pool is
reachable within a short accumulation, another hotspot indicator. Reading the
heatmap, the clusters of high slope coincide with the corners where endemics
were seeded, since those tight ranges concentrate species turnover. These maps
turn an abstract curve into a spatial prioritisation layer: rank sites by
slope or inverse half-richness and the top entries are candidate priority
areas.

The three metrics capture different stages of the curve and need not agree. The
slope reads the steepness at the origin, sensitive to immediate neighbours.
Half-richness reads the middle, where the curve crosses the halfway mark and
turnover slows. The AUC integrates the whole curve and rewards sites that
accumulate quickly and reach a high plateau. A site can have a steep slope but
a high half-richness if its rich neighbourhood is small and the curve then
stalls. Comparing the maps of the three metrics shows whether the hotspots are
local peaks or extended rich zones, a distinction that matters when a reserve
must enclose a contiguous area rather than a single point.

The metrics object exports to sf for use in a GIS workflow or for spatial joins
against habitat layers. The conversion is guarded because sf is a Suggests
dependency.

```{r metrics-sf, eval = requireNamespace("sf", quietly = TRUE)}
met_sf <- as_sf(met, crs = 32631)
class(met_sf)
```

## Moran eigenvector maps and spatial partitioning

`spatialEigenvectors()` builds the MEM basis from coordinates alone. The print
method reports the truncation threshold, the number of eigenvectors retained,
and the range of Moran's I across them.

```{r mem}
mem <- spatialEigenvectors(coords)
mem
```

The `summary()` method tabulates each eigenvector's eigenvalue, Moran's I, and
share of total variance. The leading vectors carry the broad-scale, strongly
autocorrelated patterns; later vectors describe progressively finer structure.

```{r mem-summary}
head(summary(mem), 5)
```

```{r plot-mem, eval = requireNamespace("ggplot2", quietly = TRUE), fig.cap = "First spatial eigenvectors mapped over the landscape."}
plot(mem, type = "map", n_vectors = 4) + transparent
```

Mapping the first few eigenvectors shows the scales explicitly: MEM1 splits the
region into broad halves, later vectors carve it into finer tiles. To learn
how much of a diversity pattern is spatially structured, `spatialPartition()`
regresses a per-site response on the MEMs with forward selection, then reports
the spatial and non-spatial shares of variance. The response here is the
per-site initial slope from the metrics object.

```{r partition}
slope <- met$metrics$slope_10
part <- spatialPartition(slope, mem)
part
```

```{r plot-partition, eval = requireNamespace("ggplot2", quietly = TRUE), fig.cap = "Variance in per-site slope split into spatial and non-spatial components."}
plot(part) + transparent
```

The spatial R-squared is the fraction of slope variation explained by the
selected eigenvectors, and the selected MEMs name the scales that carry it. A
large spatial share with mostly low-numbered MEMs means the hotspots are
organised at broad scales; a share spread across high-numbered MEMs means the
structure is patchy and fine-grained. Partitioning this way separates the
question "is the pattern spatial" from "is the pattern explained by a
covariate", and the two analyses can be combined to ask how much spatial
structure survives after accounting for an environmental predictor.

Forward selection guards against overfitting. It adds eigenvectors one at a
time, keeping each only when it improves the model fit by enough to justify the
extra parameter, and stops when no remaining vector helps. The result is a
parsimonious set rather than the full basis, so the reported spatial R-squared
reflects structure the data actually support. Passing `forward = FALSE` includes
every eigenvector and inflates the fit, which is appropriate only when the goal
is to remove all spatial signal rather than to describe it. The selected MEMs
are worth inspecting individually: their maps name the spatial scales at which
the response is organised, turning a single R-squared into a scale-resolved
account of where the pattern lives.

## Wavefront expansion

`wavefront()` accumulates species inside an expanding radius from seed points,
a model of spread from an introduction front. As a diagnostic it answers a
different question than the nearest-neighbour curve: how richness grows with
geographic reach rather than with site count, which exposes how clumped the
species pool is around typical starting points. The method originates in
invasion biology, where it describes the species an organism encounters as its
range front sweeps outward, but it serves equally as a spatial diagnostic for
any dataset where distance from a point is the relevant axis.

```{r wavefront}
wf <- wavefront(species, coords, n_seeds = 20, n_steps = 40,
                progress = FALSE, seed = 1)
wf
```

```{r plot-wavefront, eval = requireNamespace("ggplot2", quietly = TRUE), fig.cap = "Species captured against expanding radius from seed points."}
plot(wf) + transparent
```

The curve rises steeply at small radii when seeds land in species-dense
neighbourhoods, then flattens as the front reaches sparse ground. Plotting
against `xaxis = "sites"` instead recovers a site-based accumulation, so the
two axes give the geographic and the sampling view of the same expansion. A
front that saturates early relative to its maximum radius indicates a pool
concentrated near the seeds, the spatial signature the endemism curve also
detects.

The radius step `dr` defaults to the maximum pairwise distance divided by
`n_steps`, so the front sweeps the whole region in the requested number of
increments. The contrast between the radius axis and the site axis is the
diagnostic content. A curve that climbs fast on the radius axis but slow on the
site axis means each unit of distance brings in many sites, the signature of a
dense, evenly covered network. The opposite pattern, slow on radius and fast on
sites, marks a network with large gaps the front must cross before reaching the
next cluster. Read this way, the wavefront measures how the spatial layout of
the sampling, not just the species, shapes the pace of discovery.

## Spatial subsampling

Sites clustered in space inflate the apparent accumulation rate and bias
macroecological fits, because nearby sites share species. `subsample()` thins
the network to reduce that autocorrelation. Grid thinning overlays a mesh of a
chosen cell size and keeps one random site per occupied cell.

```{r subsample}
keep <- subsample(coords, method = "grid", cell_size = 20, seed = 1)
length(keep)

sac_full <- spacc(species, coords, n_seeds = 20, progress = FALSE, seed = 1)
sac_thin <- spacc(species[keep, ], coords[keep, ], n_seeds = 20,
                  progress = FALSE, seed = 1)
combined <- c(full = sac_full, thinned = sac_thin)
```

```{r plot-subsample, eval = requireNamespace("ggplot2", quietly = TRUE), fig.cap = "Accumulation curve before and after grid thinning."}
plot(combined) + transparent
```

The thinned curve sits above the full curve at equal site counts: with
redundant near-neighbours removed, each retained site adds more new species, so
the per-site accumulation rate climbs. The two curves converge at their right
ends toward total richness, since thinning drops sites, not species. For SAR
and SFAR fits, the thinned set gives a less biased area exponent because it
breaks the spatial redundancy that otherwise flattens the curve. Thinning gives
up sample size to gain independence, and the cell size sets the balance directly.

Grid thinning is one of three methods the function offers. Random subsampling
takes a simple sample of a target size and ignores geography, useful only for
matching sample sizes across datasets. Distance thinning removes sites greedily
until no two retained points fall closer than `min_dist`, dropping the member of
each close pair that has more near neighbours, which gives finer control than a
grid when the sampling is clustered irregularly. The grid method is the fastest
and the most predictable: one site per occupied cell, with the cell size read
directly off the autocorrelation scale. Whichever method is chosen, the return
value is a vector of row indices, so the same indices subset both the species
matrix and the coordinates and keep them aligned by position.

## Practical guidance

The spatial model family rewards attention to sample size and to whether the
question matches the method. The rules below use the conventions of this
package and the cited literature.

**Rules of thumb.**

- SFAR needs at least 4 fragments to estimate $f$ separately from $z$, and
  closer to 8 to 10 for a stable exponent. With fewer patches the fragmentation
  term and the area term are nearly collinear, and the split between them is
  unidentifiable.
- SESARS needs effort to span at least one order of magnitude (a 10-fold range
  from least to most surveyed site). When effort is near-constant, $\log(E)$ has
  no variance to fit and the effort exponent collapses toward zero with a wide
  standard error.
- For MEM, retain on the order of $n/2$ candidate eigenvectors at most for a
  network of $n$ sites, and let forward selection prune from there. Selecting
  more than half the eigenvectors risks fitting noise, since each MEM consumes a
  degree of freedom.
- Endemism curves stabilise around 30 to 50 sites; below 20 sites the early
  endemism proportion is dominated by sampling order and the percentile band is
  too wide to read.
- For grid subsampling, set `cell_size` near the spatial scale of
  autocorrelation. A cell smaller than the nearest-neighbour spacing keeps
  nearly every site and removes no redundancy; a cell several times the spacing
  discards most of the data.

**Minimum sample sizes.**

| Method          | Minimum input            | Comfortable input        |
|-----------------|--------------------------|--------------------------|
| Endemism curve  | 20 sites                 | 40 or more sites         |
| SFAR            | 4 fragments              | 8 to 10 fragments        |
| SESARS          | effort spanning 10x      | effort spanning 30x+     |
| MEM partition   | 20 sites                 | 50 or more sites         |
| Per-site metrics| 30 sites                 | 50 or more sites         |
| Wavefront       | 30 sites                 | 50 or more sites         |

**Choosing a method.**

| Question                                        | Method        |
|-------------------------------------------------|---------------|
| Which sites hold species found nowhere else     | spaccEndemism |
| Does fragmentation cost richness beyond area    | sfar          |
| Is richness biased by uneven survey effort      | sesars        |
| Where are the accumulation hotspots             | spaccMetrics  |
| At what scale is the diversity pattern organised| spatialEigenvectors + spatialPartition |
| How does richness grow with geographic reach    | wavefront     |
| How to break spatial autocorrelation            | subsample     |

**When not to use each.**

- Do not fit SFAR when patches are an arbitrary partition of continuous
  habitat, as the k-means patches here are. SFAR assumes fragments are real,
  isolated habitat units; on a continuous surface $f$ carries no biological
  meaning and the example above returns no penalty for exactly this reason.
- Do not use SESARS when effort is essentially uniform across sites. The model
  reduces to an ordinary SAR with an extra unidentified parameter, and the
  effort exponent will be noise.
- Do not read endemism proportions from fewer than 20 sites. With a short visit
  order the early steps are dominated by which site happened to be visited
  first, not by genuine spatial restriction.
- Do not feed MEM more eigenvectors than roughly half the site count. The basis
  then has enough degrees of freedom to fit any per-site response, and the
  spatial R-squared loses meaning.

## References

- Hanski, I., Zurita, G.A., Bellocq, M.I. & Rybicki, J. (2013).
  Species-fragmented area relationship. Proceedings of the National
  Academy of Sciences, 110, 12715-12720.
- Rybicki, J. & Hanski, I. (2013). Species-area relationships and extinctions
  caused by habitat loss and fragmentation. Ecology Letters, 16, 27-38.
- Dennstadt, F., Horak, J. & Martin, M.D. (2019). Predictive sampling effort
  and species-area relationship models for estimating richness in fragmented
  landscapes. Diversity and Distributions, 26, 1112-1123.
- Borcard, D. & Legendre, P. (2002). All-scale spatial analysis of ecological
  data by means of principal coordinates of neighbour matrices. Ecological
  Modelling, 153, 51-68.
- Dray, S., Legendre, P. & Peres-Neto, P.R. (2006). Spatial modelling: a
  comprehensive framework for principal coordinate analysis of neighbour
  matrices (PCNM). Ecological Modelling, 196, 483-493.
- Kier, G., Kreft, H., Lee, T.M., et al. (2009). A global assessment of
  endemism and species richness across island and mainland regions.
  Proceedings of the National Academy of Sciences, 106, 9322-9327.
- Aiello-Lammens, M.E., Boria, R.A., Radosavljevic, A., et al. (2015).
  spThin: an R package for spatial thinning. Ecography, 38, 541-545.
