Let’s practice analyzing some point patterns using simulated data!

You’ll need to load the spatstat package for this exercise:

require(spatstat)

Complete Spatial Randomness

First, let’s create an observation window, then simulate some patterns. For the examples we’ll use a square, 1-unit by 1-unit window:

ppp_owin = owin(xrange = c(0, 1), yrange = c(0, 1))

Next, let’s choose a value for the point intensity, lanbda, and simulate some CSR points:

# Set a random seed so we can re-create our results
set.seed(123456789)

lambda_sim = 100
ppp_sim_1 = 
  rpoispp(
    lambda = lambda_sim,
    win = ppp_owin)

plot(ppp_sim_1)

Prove to yourself that CSR is Poisson-distributed

We can create some evenly-spaced quadrats. Let’s try a 5 by 5 grid:

quad_25 = quadratcount(ppp_sim_1, nx = 5, ny = 5)
hist(
  quad_25,
  main = "Histogram of quadrat counts",
  xlab = "count"
)

Now let’s plot some random Poisson-distributed data using the same parameters.

We know that the true lambda for our simulated point pattern is 100 points per unit. We also know that we simulated a homogeneous Poisson process, so all the quadrats have the same value of lambda - we just have to adjust for their area.

Since we divided the area into 25 quadrats, we need to simulated 25 Poisson-distributed numbers, and the value of lambda will simply be our original point intensity divided by the number of quadrats: \(100 / 25 = 4\).

We can use the rpois() function to generate random Poisson-distributed numbers.

set.seed(1)
pois_25 = rpois(n = 25, lambda = lambda_sim / 25)
hist(
  pois_25,
  main = "Histogram of random\nPoisson-distributed numbers",
  xlab = "count")

Compare Histograms

Do the distributions look similar?

Let’s plot them one on top of the other:

par(mfrow = c(2, 1))

hist(
  quad_25,
  main = "Histogram of quadrat counts",
  xlab = "count"
)
hist(
  pois_25,
  main = "Histogram of random\nPoisson-distributed numbers",
  xlab = "count")

They both appear right-skewed (like a Poisson should), but a qualitative test isn’t enough.

Compare mean and variance

We know that the Poisson distribution has one parameter that is equal to its mean and variance. A very simple check for whether numbers are Poisson-distributed is to compare the mean and variances:

mean(quad_25)
## [1] 4.2
var(c(quad_25))
## [1] 3.916667
mean(pois_25)
## [1] 4.28
var(pois_25)
## [1] 4.626667


Note that I had to surround the quadratcount object with the c() function in order to get it to work with the var() function

The means and variances are all near 4, which is what we expected, however this isn’t a formal test.

Compare the empirical CDFs

We could also compare the empirical cumulative distribution functions:

ecdf_pts = ecdf(quad_25)
ecdf_rnd = ecdf(pois_25)
plot(
  ecdf_pts(seq(0, 12)), type = "l",
  main = "Empirical CDFs",
  ylab = "Pr(x < i)", xlab = "i")
lines(ecdf_rnd(seq(0, 12)), lty = 2)

They overlap surprisingly well, but this still isn’t a formal test!

Test for CSR

Quadrat test

To conduct a formal test for CSR, we can use the quadrat.test() function, which by default performs a Chi-square test:

quadrat.test(quad_25)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 22.381, df = 24, p-value = 0.8869
## alternative hypothesis: two.sided
## 
## Quadrats: 5 by 5 grid of tiles


Note that the null hypothesis for the quadrat.test() is CSR.

Have a Go!

Try to repeat the above exercise using different values of lambda and different numbers of quadrats.

For example, if I use a very high intensity, lambda = 1000, and a larger number of quadrats (144) I get histograms that look like this:

The ECDFs look like:

My formal test for CSR looks good:

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 143.79, df = 143, p-value = 0.9315
## alternative hypothesis: two.sided
## 
## Quadrats: 12 by 12 grid of tiles
env = envelope(ppp_sim_1, fun = Lest, 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.
plot(env, . - r ~ r)