Spatial Point Patterns

Recommended reading:

chapter 7 in Applied Spatial Data Analysis with R, second edition ; Springer (2013).
Statistical Analysis and Modelling of Spatial Point Patterns Dr. Janine Illian, Prof. Antti Penttinen, Dr. Helga Stoyan, Dietrich Stoyan
Spatial Point Patterns: methodology and applications with R. Badelley, Rubak, Turner. (The spatstat book) http://spatstat.org/

Spatial Point Pattern analysis is concerned with

describing the distribution of points (a pattern) in space;
makeing inference about the process that (could have) generated this pattern

Usually a 2-D Euclidian (metric) space is assumed:

library(sf)
## Linking to GEOS 3.7.1, GDAL 2.4.2, PROJ 5.2.0
demo(nc, ask = FALSE, echo = FALSE) # load nc
## Reading layer `nc.gpkg' from data source `/home/edzer/R/x86_64-pc-linux-gnu-library/3.6/sf/gpkg/nc.gpkg' using driver `GPKG'
## Simple feature collection with 100 features and 14 fields
## Attribute-geometry relationship: 0 constant, 8 aggregate, 6 identity
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## epsg (SRID):    4267
## proj4string:    +proj=longlat +datum=NAD27 +no_defs
pts = st_centroid(nc)
## Warning in st_centroid.sf(nc): st_centroid assumes attributes are constant
## over geometries of x
## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
## of_largest_polygon): st_centroid does not give correct centroids for
## longitude/latitude data
library(maptools) # contains the Spatial -> spatstat routines
## Loading required package: sp
## Checking rgeos availability: TRUE
library(spatstat)
## Loading required package: spatstat.data
## 
## Attaching package: 'spatstat.data'
## The following object is masked _by_ '.GlobalEnv':
## 
##     redwoodfull.extra
## Loading required package: nlme
## Loading required package: rpart
## 
## spatstat 1.61-0       (nickname: 'Puppy zoomies') 
## For an introduction to spatstat, type 'beginner'
## 
## Note: spatstat version 1.61-0 is out of date by more than 11 weeks; a newer version should be available.
as.ppp(as(pts, "Spatial"))
## Error in as.ppp.SpatialPointsDataFrame(as(pts, "Spatial")): Only projected coordinates may be converted to spatstat class objects

Over a give window \(W\), we observe a number of \(n\) events (points) \(S_i\), \(i=1,…n\).

We distinguish:

a point pattern: the observed locations of events in their window
a point process: the stochastic process that generates point pattern as its outcome

We can describe:

first order properties: concerning the density of the events (is the density of the point pattern constant over the window?)
second order properties: concerning the interaction between events (do events attract each other, repulse each other, or act indifferently with respect to each other?)

The baseline process is completely spatial random (CSR):

the density of the process is constant, and
points occur completely independently, meaning there is no interaction.

Testing for constant density can be done with the quadrat counts test:

our point pattern has \(N\) events
cut the area in \(n\) boxes of equal size
under CSR, the Expected number of points per box is then \(E_i = N/n\)
count the number of Observed events per box, \(O_i\)
under the null hypothesis of CSR, \(G = \sum_{i=1}^{n}\frac{(O_i - E_i)^2}{E_i}\) now follows a \(\chi^2\) distribution with \(n-1\) degrees of freedom

In case you are unfamiliar with what \(\chi^2\) is, an approximate permutation test involves

generate a large (\(>1000\)) number of independent CSR patterns of size \(N\) over the target region (spatstat::runifpoint), and consider them observed,
compute \(G\) for each of them, and take the 95-percentile of these \(G\) values
compare the test statistic from the real observed pattern: is it larger than this percentile, we can reject the CSR hypothesis with 95% confidence

Other tests involve e.g. distances between points.

Caveats:

how to determine \(n\)?
can we get equally sized boxes in real cases?
as usual: not significant does not mean that the process is CSR, it can also indicate lack of power (small sample)

Alternatives to CSR:

attracting or repulsing processes, constant density (e.g. Strauss)
no interaction process, varying density (log Gaussian Cox process)

Software:

spatial: spatstat (JSS papers :1, 2)
spatio-temporal: CRAN: stpp JSS paper; lgcp JSS papers: 1, 2

Typical spatio-temporal questions involve whether the density (if varying) varies over time in a way that is independent from the variation over space. In case of an epidemics, movement of the epidemic would imply interaction (FMD examples in the stpp package/paper).

Functions:

K-function (Ripley's K): the average number of other points found within the distance \(r\) from the typical point.
L-function: similar to K, but CSR is a line rather than a parabolic curve
g-function: (pair-correlation function); recommended by Illian et al., similar to L and K but \(g® = 1\) for CSR cases

Inference:

inference often follows a simulation approach: point patterns are simulated given a particular process, and the correspondence of the (set of) simulations to the observed processed is considered.

library(spatstat)
data(japanesepines)
summary(japanesepines)
## Planar point pattern:  65 points
## Average intensity 65 points per square unit (one unit = 5.7 metres)
## 
## Coordinates are given to 2 decimal places
## i.e. rounded to the nearest multiple of 0.01 units (one unit = 5.7 metres)
## 
## Window: rectangle = [0, 1] x [0, 1] units
## Window area = 1 square unit
## Unit of length: 5.7 metres

library(maptools)
spjpines <- as(japanesepines, "SpatialPoints")
summary(spjpines)
## Object of class SpatialPoints
## Coordinates:
##      min max
## [1,]   0 5.7
## [2,]   0 5.7
## Is projected: NA 
## proj4string : [NA]
## Number of points: 65

spjpines1 <- elide(spjpines, scale=TRUE, unitsq=TRUE)
summary(spjpines1)
## Object of class SpatialPoints
## Coordinates:
##      min max
## [1,]   0   1
## [2,]   0   1
## Is projected: NA 
## proj4string : [NA]
## Number of points: 65

pppjap <- as(spjpines1, "ppp")
summary(pppjap)
## Planar point pattern:  65 points
## Average intensity 65 points per square unit
## 
## Coordinates are given to 2 decimal places
## i.e. rounded to the nearest multiple of 0.01 units
## 
## Window: rectangle = [0, 1] x [0, 1] units
## Window area = 1 square unit

data(redwoodfull)
spred <- as(redwoodfull, "SpatialPoints")
data(cells)
spcells <- as(cells, "SpatialPoints")
dpp<-data.frame(rbind(coordinates(spjpines1), coordinates(spred), 
   coordinates(spcells)))
njap<-nrow(coordinates(spjpines1))
nred<-nrow(coordinates(spred))
ncells<-nrow(coordinates(spcells))
dpp<-cbind(dpp,c(rep("JAPANESE",njap), rep("REDWOOD", nred), rep("CELLS", ncells))) 
names(dpp)<-c("x", "y", "DATASET")


library(lattice)
## 
## Attaching package: 'lattice'
## The following object is masked from 'package:spatstat':
## 
##     panel.histogram
print(xyplot(y~x|DATASET, data=dpp, pch=19, aspect=1))

plot of chunk unnamed-chunk-3


envjap <- envelope(as(spjpines1, "ppp"), fun=Gest, r=r, nrank=2, nsim=99)
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
envred <- envelope(as(spred, "ppp"), fun=Gest, r=r, nrank=2, nsim=99)
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
envcells <- envelope(as(spcells, "ppp"), fun=Gest, r=r, nrank=2, nsim=99)
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
Gresults <- rbind(envjap, envred, envcells) 
Gresults <- cbind(Gresults, 
y=rep(c("JAPANESE", "REDWOOD", "CELLS"), each=length(r)))
## # CHANGED DATASET TO y RSB
## #    save(Gresults, envjap, envred, envcells, file="sppaGestEnv.RData")


###################################################
### code chunk number 21: sppa.Rnw:562-563
###################################################
#load("Gresults.RData")

print(xyplot(obs~theo|y , data=Gresults, type="l", 
xlab = "theoretical", ylab = "observed", # EJP
panel=function(x, y, subscripts) {
   lpolygon(c(x, rev(x)), 
   c(Gresults$lo[subscripts], rev(Gresults$hi[subscripts])),
   border="gray", col="gray"
)
llines(x, y, col="black", lwd=2)
}))

plot of chunk unnamed-chunk-3

density:

plot(density(as.ppp(spred)), axes=TRUE)
points(as.ppp(spred), col="green2", pch=19)

plot of chunk unnamed-chunk-4

K-functions, homogeneous, inhomogeneous

opar <- par(mfrow=c(1,2))
plot(envelope(as.ppp(spred), Kest, verbose=FALSE), main="Homogeneous")
plot(envelope(as.ppp(spred), Kinhom, verbose=FALSE), main="Inhomogeneous")

plot of chunk unnamed-chunk-5

plot(density(as.ppp(spred)), axes=TRUE)
points(as.ppp(spred), col="green2", pch=19)

plot of chunk unnamed-chunk-6

plot(density(as.ppp(spred), adjust = 0.5), axes=TRUE,  main = "0.5 x the default bandwidth")
points(as.ppp(spred), col="green2", pch=19)

plot of chunk unnamed-chunk-6

opar <- par(mfrow=c(1,2))
plot(envelope(as.ppp(spred), Kest, verbose=FALSE), main="Homogeneous")
plot(envelope(as.ppp(spred), Kinhom, adjust = 0.5, verbose=FALSE), main="Inhomogeneous")

plot of chunk unnamed-chunk-7