# Optimally thin or subset sample points

I have several hundred points that represent samples in a landscape: As you can see, the points are unevenly spread and form large, dense clusters. To account for spatial autocorrelation, I want to subset or thin these points so that the points that remain are a minimum distance from each other. At the moment I've set up a pairwise distance matrix for them and am randomly removing points that are below the distance threshold until the only ones remaining meet my minimum distance criteria. However, this does not seem optimal from two different perspectives: first, different potential random subsets using this method can return different numbers of samples. Ideally, I'd retain as many points as possible. Second, different subsets may be distributed differently across the landscape. Ideally, I want a point subset with a distribution that more uniformly represents the landscape. Are there any recommended or general algorithms for optimally subsetting points (at least more optimal than my approach)?

There is a function `pp.subsample` in the spatialEco package that will subset your data using a spatial intensity function. This will, in effect, thin your data to a defined n, using the expected density as sample weights.

However, in looking at the spatial distribution of your data, unless you subsample to a handful of observations, this will not mitigate the spatial autocorrelation (pseudoreplication) issues in this particular data. You have a notable sample distribution problem.

Here is a worked example, from the functions help, that illustrates how this works.

First, lets add the required packages and some data to work with.

``````require(spatialEco)
require(spatstat)
require(sp)

data(bei)
trees <- as(bei, 'SpatialPoints')
plot(trees, pch=20, col='black', cex=0.50)
box()
title("Observed population")
`````` Now we can set a 10% subsample target, which is fairly arbitrary and one should define n based on their data and analysis goals. We then run the `pp.subsample` function to use the intensity estimate as sample weights in pulling the subsample.

``````n = round(length(trees) * 0.10, digits=0)
trees.wrs <- pp.subsample(trees, n=n, window='hull')
``````

Here are the points plotted on the estimated intensity function. Please note that the sigma parameter is quite important. I am defaulting to the first-order Scott bandwidth estimate. This smooths the point process quite a bit. If you want more of a second-order process, either specify your own sigma value, or use a second-order sigma estimate such as diggle or stoyan. Plotting the results illustrates the expected spatially weighted sample distribution and a notable thinning of the data.

``````plot(trees, pch=20, col='black', cex=0.50)
plot(trees.wrs, pch=20, col='red', cex=1.25, add=TRUE)
box()
title('10% subsample (n=360) of population (n=3604)')
legend('bottomright', legend=c('Original sample', 'Subsample'),
col=c('black','red'),pch=c(20,20),bg="white")
`````` • Thank you for the detailed answer. It's certainly better than the approach I had when I wrote the question, and better than anything else I was able to find online. I ended up figuring out my own algorithm that has the specific properties I'm looking for, and will also post it as an answer when I publish my work using it. – anjama Dec 20 '17 at 17:46