[GIS] Test for spatial randomness with Ripley’s K or L function in R

rspatial statistics

I have a spatial point pattern that I would like to test against complete spatial randomness.
I managed to plot a Ripley's K function for my data and a function for a random distribution with a confidence envelope. From this I can infer whether or not my data deviates from spatial randomness.

However, I wonder whether it is possible to statistically test these two functions against one another (and thus get a p-value that indicates the deviation from a random distribution)?

Best Answer

Getting a p value is not that easy. Marcon and Lang (Testing randomness of spatial point patterns with the Ripley statistic, http://arxiv.org/abs/1006.1567) have demonstrated that (K1, K2..Kn) where K1, K2..Kn are the values of the Ripley s K function at different distances is a Gaussian vector. They have then computed its mean and covarianace matrix in the square and used a chi 2 test for spatial randomness. You can use the same idea by computing empirically (with Monte-Carlo simulations) the mean/covariance matrix in your domain. Regards Thibault Lagache

Related Solutions

[GIS] Simulation envelopes and significance levels

The +1 is a convention. It's all about converting the ranks to percentiles. Consider 99 iterations. The rank will run from 1 through 99 (in whole steps). You can convert the rank to a percentile by dividing by 99 and multiplying by 100. That would produce percentiles from 100/99 = 1.01% to 99*100/99 = 100%. That lacks a desirable symmetry: you're saying the lowest value is 1.01% from the bottom but the highest value is right smack at the top of the range. To restore that symmetry, you should shift the ranks down by 0.5: the lowest value gets assigned a percentile of (1-0.5)*100/99 = 0.50% and the highest value gets a percentile of (99-0.5)*100/99 = 99.50% = 100 - 0.50%. Everything becomes nice and symmetric. You can avoid adjusting the ranks simply by adding 1 to the denominator. Now the percentiles in the example range from 1*100/100 = 1% up to 99% in even, symmetrical, 1% steps. (I illustrate this, and discuss it more generally, at http://www.quantdec.com/envstats/notes/class_02/characterizing_distributions.htm : look towards the bottom of that page.)
The purpose of the simulation is to compute what is likely to happen by chance alone (the "null distribution"). For instance, we can ask whether Joe Paterno's football bowl record at Penn State of 24 wins in 37 attempts is just the result of blind luck, no better (or worse) than if each game had been decided by a fair coin flip. To do so, we might actually flip a coin 37 times, count the heads (to represent wins), and repeat this for a total of 99 iterations. We would find that very few of those 99 attempts achieved 24 or more heads--most likely 5 or fewer. That's pretty good evidence that over time his teams were better than the bowl competition, indicating they tended to be under-rated. (99 is too small a number to use here, really: I actually used 100,000, in which I observed 4972 occasions of 24 or more heads.)

The one aspect of this example that you control is the weight of evidence: is a 5/100 chance enough to convince you the result was not due to randomness (or luck)? Depending on the circumstance, some people need stronger evidence, others weaker. That's the role of m. When you draw the envelopes at the most extreme out of (say) 99 iterations, then you are estimating the lowest 1% and the highest 1%. You would guess that random fluctuations would place the curve between these envelopes 100 - 1 - 1 = 98% of the time. This corresponds (very roughly, because 99 iterations is still too little) to a "two-sided p-value of 0.02." If you don't need such strong evidence, you might choose (say) m = 3. Now the lower envelope represents the bottom 3/99 = 3.03% of the chance distribution and the upper envelope represents the top 3.03% of that distribution. Your two-sided p-value is around 6%. (Because 99 is so small, the true p-value could be as large as 15% or so, so beware! Do many more Monte-Carlo iterations if you possibly can.)
In some sense you're trying to depict a distribution of envelopes (curves). That's a complicated thing. One way is to choose some percentiles in advance, such as 1% (and therefore, by symmetry, 99%), 5% (and 95%), 10% (and 90%), 25% (and 75%). Out of 99 simulated curves these would correspond roughly to the lowest (and highest), the fifth lowest (and fifth highest), etc. (Comparing curves that can intersect each other in this way is problematic, though: when they cross, which is more extreme than the other?) Plotting those selected curves at least gives one a visual idea of the spread that is produced by chance mechanisms alone.

I hope you get the sense that this approach of generating a small bunch of curves (under the null hypothesis) and drawing a selected few is exploratory, approximate, and somewhat crude, because it is all of those. But it's far better than just guessing whether or not your data could have arisen by chance.

ArcGIS Desktop – How to Implement Bivariate Ripley’s K Function

After much searching in the back corners of ESRI documentation, I've concluded that there is no reasonable way run a bivariate Ripley's K function in Arcpy/ArcGIS. However, I have found a solution using R:

# Calculates an estimate of the cross-type L-function for a multitype point pattern.
library(maptools)
library(spatstat)
library(sp)

# Subset certain areas within a points shapefile.  In this case, features are grouped by gap number
gap = 1

# Read the shapefile
sdata = readShapePoints("C:/temp/GapPoints.shp")  #Read the shapefile
data = sdata[sdata$SITE_ID == gap,]  # segregate only those points in the given cluster

# Get the convex hull of the study area measurements
gapdata = readShapePoints("C:/temp/GapAreaPoints_merged.shp")  #Read the shapefile that is used to estimate the study area boundary
data2 = gapdata[gapdata$FinalGap == gap,]  # segregate only those points in the given cluster
whole = coordinates(data2) # get just the coords, excluding other data
win = convexhull.xy(whole) # Convex hull is used to get the study area boundary
plot(win)

# Converting to PPP
points = coordinates(data) # get just the coords, excluding other data
ppp = as.ppp(points, win) # Convert the points into the spatstat format
ppp = setmarks(ppp, data$SPECIES) # Set the marks to species type YB or EH
summary(ppp) # General info about the created ppp object
plot(ppp) # Visually check the points and bounding area

# Plot the cross type L function
# Note that the red and green lines show the effects of different edge corrections
plot(Lcross(ppp,"EH","YB"))

# Use the Lcross function to test the spatial relationship between YB and EH
L <- envelope(ppp, Lcross, nsim = 999, i = "EH", j = "YB")
plot(L)