Setup

Packages

Install any packages you might be missing:

library(spatstat)
library(splancs)

Data bramblecanes, finpines, and spruces

# tree data

data(finpines)
#get info
print(finpines)
## Marked planar point pattern: 126 points
## Mark variables: diameter, height 
## window: rectangle = [-5, 5] x [-8, 2] metres
intensity(finpines)
## [1] 1.26
summary(finpines)
## Marked planar point pattern:  126 points
## Average intensity 1.26 points per square metre
## 
## Coordinates are given to 13 decimal places
## 
## Mark variables: diameter, height
## Summary:
##     diameter         height     
##  Min.   :0.000   Min.   :0.800  
##  1st Qu.:1.000   1st Qu.:1.825  
##  Median :2.000   Median :2.850  
##  Mean   :2.532   Mean   :2.828  
##  3rd Qu.:3.000   3rd Qu.:3.600  
##  Max.   :7.000   Max.   :5.400  
## 
## Window: rectangle = [-5, 5] x [-8, 2] metres
## Window area = 100 square metres
## Unit of length: 1 metre
#creates point pattern data format from the x and y coordinates
finpines.pts<-as.points(finpines$x,finpines$y)
plot(finpines.pts, pch=19, xlab="X", ylab="Y",main="Finish Pines Locations")

Exploring patterns

Plot and count Finpines data in quadrats

#divide into quadrants and count number of points in the quadrant
Q <- quadratcount(finpines, nx = 4, ny = 3)
Q
##                x
## y               [-5,-2.5) [-2.5,0) [0,2.5) [2.5,5]
##   [-1.33,2]            10        8      11      21
##   [-4.67,-1.33)         8        7       9       3
##   [-8,-4.67)           12       11      22       4
summary(Q)
## Number of cases in table: 126 
## Number of factors: 2 
## Test for independence of all factors:
##  Chisq = 20.215, df = 6, p-value = 0.002535
Q_dat<-data.frame(Q)
var(Q_dat$Freq)
## [1] 33.72727
mean(Q_dat$Freq)
## [1] 10.5
plot(finpines.pts,pch=19, main="Finpines", xlab="X",ylab="Y")
plot(Q, add = TRUE, cex = 2)

Compare Finpines to uniformly distributed points and regular points

# Finpines data
sp_point <- matrix(NA, nrow=length(finpines$x),ncol=2)
sp_point[,1] <- finpines$x
sp_point[,2] <- finpines$y
colnames(sp_point) <- c("x","y")
plot(x=sp_point[,1],y=sp_point[,2],main="Finpines Data", xlab="X",ylab="Y",cex=.5)

intensity(finpines)
## [1] 1.26
summary(finpines)
## Marked planar point pattern:  126 points
## Average intensity 1.26 points per square metre
## 
## Coordinates are given to 13 decimal places
## 
## Mark variables: diameter, height
## Summary:
##     diameter         height     
##  Min.   :0.000   Min.   :0.800  
##  1st Qu.:1.000   1st Qu.:1.825  
##  Median :2.000   Median :2.850  
##  Mean   :2.532   Mean   :2.828  
##  3rd Qu.:3.000   3rd Qu.:3.600  
##  Max.   :7.000   Max.   :5.400  
## 
## Window: rectangle = [-5, 5] x [-8, 2] metres
## Window area = 100 square metres
## Unit of length: 1 metre
# Random points
u.x <- runif(n=nrow(sp_point), min=bbox(sp_point)[1,1], max=bbox(sp_point)[1,2])
u.y <- runif(n=nrow(sp_point), min=bbox(sp_point)[2,1], max=bbox(sp_point)[2,2])

# Regular points
r.x <- seq(from=min(sp_point[,1]),to=max(sp_point[,1]),length=sqrt(nrow(sp_point)))
r.y <- seq(from=min(sp_point[,2]),to=max(sp_point[,2]),length=sqrt(nrow(sp_point)))
r.x <- jitter(rep(r.x,length(r.x)),.001)
r.y <- jitter(rep(r.y,each=length(r.y)),.001)

# Plot Finpines, uniform, and regular points
par(mfrow=c(1,3),mar=c(4,4,1.5,0.5))
plot(x=sp_point[,1],y=sp_point[,2],main="Finpines Data", xlab="X",ylab="Y",cex=.5,pch=19)
plot(x=u.x,y=u.y,main="Random Points", xlab="X",ylab="Y",cex=.5,pch=19)
plot(x=r.x,y=r.y,main="Regular Points", xlab="X",ylab="Y",cex=.5,pch=19)

Intro to homogeneous Poisson processes

Uniform point process

# CSR from random uniform point process
## Generates a CSR point process conditionally on a fixed number N of points

## Let N = 100
par(mfrow=c(2,2),mar=c(0.4,0.4,1,0.4))
pp.100.unif<-runifpoint(100,win=owin(c(0,1),c(0,1)))
pp.100.unif
## Planar point pattern: 100 points
## window: rectangle = [0, 1] x [0, 1] units
intensity(pp.100.unif)
## [1] 100
plot(pp.100.unif, pch=19, cex=0.6,main="")
title('100 uniform points')


# repeat a few times to see different realizations of the process

CSR from homogeneous poisson process (HPP)

# HPP with intensity 100 in the unit square
pp.100.unit <- rpoispp(100,win=owin(c(0,1),c(0,1)))
pp.100.unit
## Planar point pattern: 109 points
## window: rectangle = [0, 1] x [0, 1] units
intensity(pp.100.unit) 
## [1] 109
plot(pp.100.unit,main="")
title('HPP intensity = 100, unit square')

# repeat a few times to see different realizations of the process (including different numbers of points)

# HPP with intensity 1 in a 10 x 10 square
pp.1.10 <- rpoispp(1, win=owin(c(0,10),c(0,10)))

par(mfrow=c(1,2))
plot(pp.100.unit,main="")
title('intensity = 100, unit square')
plot(pp.1.10,main="")
title('intensity = 1, 10 x 10 square')

# HPP with intensity 200 in the unit square
par(mfrow=c(1,1))
pp.200.unit <- rpoispp(200)
plot(pp.200.unit)

Simulate pi by counting points from a HPP within a circle inscribed in a square

# HPP with intensity 100 in the unit square
side.square <-1
pp.100.unit <- rpoispp(10000,win=owin(c(0,side.square),c(0,side.square)))
points.tot <- npoints(pp.100.unit) 
plot(pp.100.unit,main="")

# inscribe a circle within the square
radius.circle <-  side.square/2
win.circle <- disc(radius.circle, c(radius.circle,radius.circle))

pp.incircle <- pp.100.unit[win.circle]
points.incircle <- npoints(pp.incircle)
summary(pp.incircle)
## Planar point pattern:  7943 points
## Average intensity 10117.4 points per square unit
## 
## Coordinates are given to 8 decimal places
## 
## Window: polygonal boundary
## single connected closed polygon with 128 vertices
## enclosing rectangle: [0, 1] x [0, 1] units
## Window area = 0.785083 square units
## Fraction of frame area: 0.785
plot(pp.incircle)

# calculate pi using ratio of areas and number of points
piEst <- ((side.square^2)*points.incircle) / ((radius.circle^2)*points.tot)
piEst
## [1] 3.148236

Function to simulate pi over many iterations

pi.simulation2 = function(iter, lambda, side.square){
  
  piEst.list = rep(0,iter)
  points.incircle.list = rep(0,iter)
  
  # inscribe a circle within the square
  radius.circle <-  side.square/2
  win.circle <- disc(radius.circle, c(radius.circle,radius.circle))
  
  for(i in 1:iter){
    # HPP with intensity 100 in the unit square
    pp.square <- rpoispp(lambda,win=owin(c(0,side.square),c(0,side.square)))
    points.tot <- npoints(pp.square) 
    
    # find points within the circle
    pp.incircle <- pp.square[win.circle]
    points.incircle <- npoints(pp.incircle)

    # calculate pi using ratio of areas and number of points
    pihat <- ((side.square^2)*points.incircle) / ((radius.circle^2)*points.tot)
    piEst.list[i] = pihat
  }
  return( piEst.list)
}

iter = 1000
lambda = 1000
side.square = .5
est = pi.simulation2(iter,lambda, side.square)
mean(est)
## [1] 3.146699
var(est)
## [1] 0.01027661
# we plot the histogram and see it appears to be binomial. why?
hist(est,nclass = 15)

Tests for CSR

Testing for CSR with Ripley’s K

# create polygon (square) area for finpines data
minfx <- min(finpines$x)
maxfx <- max(finpines$x)
minfy <- min(finpines$y)
maxfy <- max(finpines$y)

polyfin <- as.points(c(minfx,maxfx,maxfx,minfx),c(minfy,minfy,maxfy,maxfy))
UL.khat.fin <- Kenv.csr(length(finpines$x), polyfin, nsim=99, seq(0.05,3.5,0.15))
## Doing simulation  1 
## Doing simulation  2 
## Doing simulation  3 
## Doing simulation  4 
## Doing simulation  5 
## Doing simulation  6 
## Doing simulation  7 
## Doing simulation  8 
## Doing simulation  9 
## Doing simulation  10 
## Doing simulation  11 
## Doing simulation  12 
## Doing simulation  13 
## Doing simulation  14 
## Doing simulation  15 
## Doing simulation  16 
## Doing simulation  17 
## Doing simulation  18 
## Doing simulation  19 
## Doing simulation  20 
## Doing simulation  21 
## Doing simulation  22 
## Doing simulation  23 
## Doing simulation  24 
## Doing simulation  25 
## Doing simulation  26 
## Doing simulation  27 
## Doing simulation  28 
## Doing simulation  29 
## Doing simulation  30 
## Doing simulation  31 
## Doing simulation  32 
## Doing simulation  33 
## Doing simulation  34 
## Doing simulation  35 
## Doing simulation  36 
## Doing simulation  37 
## Doing simulation  38 
## Doing simulation  39 
## Doing simulation  40 
## Doing simulation  41 
## Doing simulation  42 
## Doing simulation  43 
## Doing simulation  44 
## Doing simulation  45 
## Doing simulation  46 
## Doing simulation  47 
## Doing simulation  48 
## Doing simulation  49 
## Doing simulation  50 
## Doing simulation  51 
## Doing simulation  52 
## Doing simulation  53 
## Doing simulation  54 
## Doing simulation  55 
## Doing simulation  56 
## Doing simulation  57 
## Doing simulation  58 
## Doing simulation  59 
## Doing simulation  60 
## Doing simulation  61 
## Doing simulation  62 
## Doing simulation  63 
## Doing simulation  64 
## Doing simulation  65 
## Doing simulation  66 
## Doing simulation  67 
## Doing simulation  68 
## Doing simulation  69 
## Doing simulation  70 
## Doing simulation  71 
## Doing simulation  72 
## Doing simulation  73 
## Doing simulation  74 
## Doing simulation  75 
## Doing simulation  76 
## Doing simulation  77 
## Doing simulation  78 
## Doing simulation  79 
## Doing simulation  80 
## Doing simulation  81 
## Doing simulation  82 
## Doing simulation  83 
## Doing simulation  84 
## Doing simulation  85 
## Doing simulation  86 
## Doing simulation  87 
## Doing simulation  88 
## Doing simulation  89 
## Doing simulation  90 
## Doing simulation  91 
## Doing simulation  92 
## Doing simulation  93 
## Doing simulation  94 
## Doing simulation  95 
## Doing simulation  96 
## Doing simulation  97 
## Doing simulation  98 
## Doing simulation  99
khat.fin<-khat(finpines.pts, polyfin, seq(0.05,3.5,0.15))
# plot of Khat-pi t^2 from data
plot(seq(0.05,3.5,0.15), khat.fin-pi*seq(0.05,3.5,0.15)^2, type="l", xlab="Distance", ylab="Estimated K-pi h^2")

# plot upper bound
lines(seq(0.05,3.5,0.15), UL.khat.fin$upper-pi*seq(0.05,3.5,0.15)^2, lty=2)

# plot lower bound
lines(seq(0.05,3.5,0.15), UL.khat.fin$lower-pi*seq(0.05,3.5,0.15)^2, lty=2)

Other plots with Ripley’s K (Variance stabilizing)

# L functions
l<-function(k,h)
{sqrt(k/pi)-h}

# plot Lhat from data
plot(seq(0.05,3.5,0.15), l(khat.fin, seq(0.05,3.5,0.15)), type="l", xlab="Distance", ylab="Estimated L",ylim=c(-.1,.2))

# plot upper bound of Lhat
lines(seq(0.05,3.5,0.15), l(UL.khat.fin$upper,seq(0.05,3.5,0.15)), lty=2)

# plot lower bound of Lhat
lines(seq(0.05,3.5,0.15), l(UL.khat.fin$lower,seq(0.05,3.5,0.15)), lty=2)

Different R function for Ripley’s K

K<-Kest(finpines)
plot(K)

plot(K, cbind(r, sqrt(iso/pi)) ~ r)

plot(K, cbind(trans,iso,border) - theo ~ r)

plot(envelope(finpines,Kest))
## 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.

distance-based functions testing for CSR

# G test: Nearest neighbor distribution test
Gtest<-Gest(finpines)
plot(envelope(finpines,Gest))
## 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.

# F test: Point to nearest event distribution test
Ftest<-Fest(finpines)
plot(envelope(finpines,Fest))
## 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.

# J test: J(r) = (1-G(r)/(1-F(r)))
Jtest<-Jest(finpines)
plot(envelope(finpines,Jest))
## 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.