In antiphon/Kdirectional: Analysis of anisotropy in point patterns using second order statistics

Originally written around 2014. Updated 10.12.2022.

Introduction

The package Kdirectional is a collection of tools for analysing the isotropy of 2- and 3-dimensional point patterns. This document shows how to use some of these tools by analysing an synthetic example pattern.

In this vignette we use the following packages:

library(spatstat)
library(Kdirectional)

library(spatstat)
library(Kdirectional)

Example pattern

We simulate Strauss process and transform it. Assume the transformation is

comp <- 0.7
D <- diag( c(comp, 1/comp))
angle <- pi/8
R <- rotationMatrix3(az=angle)[-3,-3] # internal, not rgl one
C <- R%*%D

that is, we compress the x-axis by factor 0.7 and stretch the y-axis to keep volume 1. Then we simulate the original process in a window that will be square after the transformation

set.seed(1)
x0 <- rStrauss(500, gamma = 0.05, R = 0.03, W = affine(square(), solve(C)))

and transform it to the observed pattern

x <- affine(x0, C)

Check things went ok (define a little helper for plotting):

plotp <- function(...) plot(..., cex=.3)
par(mfrow=c(1,2), mar=c(1,1,1,1))
plotp(x0)
plotp(x)

The pattern on the left is what we want to estimate.

Nearest neighbour vectors

Check for starters the distribution of nearest neighbours vectors.

nna <- nnangle(coords(x))
ang <- nna[,1]-pi
par(mfrow=c(1,2))
angles.flower(ang, ci = T, k=24, bty = "n")
# Compare to spatstat
rose(ang, unit = "rad", nclass = 24+1, main = "")

The blue rings give approximate 95% pointwise confidence intervals. Hard to see anything definite, perhaps some peakyness at the 2:30'oclock axis. Nearest neighbours are not very informative, so let's analyse the full set of pairwise vectors.

Fry-points and ellipsoids

As in the paper [1], the Fry-points tell us interesting things about the isotropy.

f <- fry_ellipsoids(x, nvec=1:20, nangles = 30, eps=0.1, r_adjust = .2)

r_adjust=.2 reduces the computations to pairs with length within 20% window width, eps adds a bit of angular smoothing (in radians), and nangles chooses the amount of directions to look at.

Plot the estimated ellipsoids:

plot(f, zoom=.3, used_points = FALSE)

Oval shape is visible, some ambiguity of rotation. Check the ellipticity with a contrast $axis_x-axis_y=0$:

ci <- confint(f, nsim=500)
print(summary(ci)[1:5,], digits=2)
plot(ci, ylim = c(-1,1)*.1)

The first contour is poorly estimated as there are some short range pairs, kind of like Poisson noise. Others seem to hint of compression (the p-values and confidence intervals are biased and conservative).

Let's compute the average rotation, where we scale each ellipsoid to have det = 1 by setting just_rotation=TRUE.

contours <- f$ellipsoids[-1]
mean_e <- mean_ellipsoids(contours, just_rotation = TRUE)
summary(mean_e)

Rotation seems to be non-zero (or non $\pm\pi$). Ignore the error variance, its not correct.

c(true=angle, est=mean_e$rot_angle + pi)

Backrotate the pattern:

Rhat <- rotationMatrix3(az=pi)[-3,-3]%*%mean_e$rot 
xr <- affine(x, solve(Rhat))

Next we estimate the compression by minimizing the difference of K-functions along the main axes. We need to provide a grid to optimise over:

gamma_vec <- seq(0.5, 1, by=0.05)
grid <- cbind(gamma_vec, 1/gamma_vec)

Then compute the anisotropy profile:

r <- seq(0, 0.2, length=50)
eps <- pi/10
ani <- anisotropy_profile(xr, grid = grid, r=r, eps=eps, power=2)
ani0 <- anisotropy_profile(x, grid = grid, r=r, eps=eps, power=2)

We set the sector-K angle with eps and use integral of squared differences (power=2) as the distance.

Plot the profiles, scaled to 0-1 for plotting (it's the default):

plot(ani, type="b", scale=TRUE)
lines(ani0, type="b", col=2)

( s <- summary(ani) )

Estimated compression

c(true=comp, est=s$estimate[1])

Backtransform:

Dhat <- s$estimate
xhat <- affine(xr, solve(Dhat))

Recap of the steps:

par(mfrow=c(2,2), mar=c(1,1,1,1))
plotp(x)
plotp(xr)
plotp(xhat)
plotp(x0)

Finally, let's illustrate the estimated transformation. Use simulation and estimation to get an ellipsoid object:

Chat <- Rhat %*% Dhat
e0 <- ellipsoid_OLS( rellipsoid(1000, c(1,1), R = C) )
e <- ellipsoid_OLS( rellipsoid(1000, c(1,1), R = Chat) )

plot(e0, add=FALSE)
plot(e, col=3)

References

[1] T Rajala, A Särkkä, C Redenbach, M Sormani (2015): Estimating geometric anisotropy in spatial point patterns, Spatial Statistics

antiphon/Kdirectional documentation built on Feb. 13, 2023, 6:26 a.m.