Q1: Read in event observations and plot a point process dataframe

library(spatstat)
## Warning: package 'spatstat' was built under R version 3.2.3
## Loading required package: nlme
## 
## spatstat 1.44-1       (nickname: 'Gift Horse') 
## For an introduction to spatstat, type 'beginner'
setwd("C:\\Users\\Yang\\Desktop\\Spatial Analysis\\Lab Work")
sppdat <- read.table("Lab_01.txt", header=T)
as.ppp(sppdat, c(min(sppdat$x), max(sppdat$x), min(sppdat$y), max(sppdat$y)))
## Planar point pattern: 19 points
## window: rectangle = [3, 9] x [2, 9] units
plot(sppdat)

Q2: Generate poison, random, SSI and regular matern process and plot them Sequential Spatial Inhibition and Matern type I are both regular ones

set.seed(4887)
spois <- rpoispp(50)
srand <- rpoint(50)
sssi <- rSSI(.1, 50, win = owin(c(0,1), c(0,1)))
smaternI <- rMaternI(50, .02, win=owin(c(0,1), c(0,1)))

par(mfrow=c(2,2), mex=.5)
plot(spois)
plot(srand)
plot(sssi)
plot(smaternI)

By changing the parameter to increase the inhibition distance, thus making the points look more regular

scluster <- rThomas(17, .08, 3)
scluster2 <- rMatClust(17, .08, 3)

par(mfrow=c(1,2), mex=.4)
plot(scluster)
plot(scluster2)

generate and plot my own processes (This function generates n independent, identically distributed random points with common probability density proportional to f; f could be numeric constant or a function)

sownProp <- rpoint(50, function(x,y) {x^2 + y^2})
par(mfrow=c(1,1), mex=.5)
plot(sownProp)

Q3: Perform the quadrat test on two Matern Cluster processes

Generate and cluster process and plot them

clus1 <- rMatClust(30, .04, 4)
clus2 <- rMatClust (30, .08, 4)
par(mfrow=c(1,2), mex = .5)
plot(clus1)
plot(clus2)

Counting events on quadrats

count1 <- quadratcount(clus1, nx=2, ny=2)
count2 <- quadratcount(clus2, nx=2, ny=2)
exp1 <- sum(count1)/4
exp2 <- sum(count2)/4

Quadrat test

chi1 <- sum((count1 - exp1)^2)/(exp1)
chi2 <- sum((count2 - exp2)^2)/(exp2)

2 * (1-pchisq(chi1, 3))
## [1] 0.0001710985
2 * (1-pchisq(chi2, 3))
## [1] 0.6512354
ID1 <- sum((count1 - exp1)^2)/3/(exp1)
ID2 <- sum((count2 - exp2)^2)/3/(exp2)

ID1
## [1] 7.144578
ID2
## [1] 1.15427
2 * (1-pchisq(3*ID1, 3))
## [1] 0.0001710985
2 * (1-pchisq(3*ID2, 3))
## [1] 0.6512354

According the class note, ID1 is >> 1, so we think that the piont process is cluster; ID2 is close to 1, the point process is CSR. And P value for cluster1 is small, we can not accept H0 that is it is not CSR; P value for cluster2 is high, we accept the H0 that it is CSR

Q4: Peform Poisson lack-of-fit test to test the simdat (Jananese black pine data) to figure out if these data are CSR and at what scale

First, just look at the data

data(simdat)
plot(c(0,10), c(0,10), xlim=c(0,10), ylim=c(0,10), type="n", xlab = "x", ylab = "y")
plot(simdat, add=T)
title(main="simdat")

By looking at the data pattern, it seems not CSR

Using lack-of-fit test to Poisson Distribution and we have 7 categories here

cut <- c(0:5)/.5
xcut <- cut(simdat$x, cut)
class(xcut)
## [1] "factor"
ycut <- cut(simdat$y, cut)
mytab <- table(xcut, ycut)
mytab
##         ycut
## xcut     (0,2] (2,4] (4,6] (6,8] (8,10]
##   (0,2]      5    10     2     3      8
##   (2,4]      9     2    12    14     11
##   (4,6]     12     8     7    10     11
##   (6,8]      9     9     6    10     10
##   (8,10]     7     7     3     8      4
uhat <- sum(mytab) / length(mytab)
expected <- rep(0,7)
observed <- rep(0,7)

for(k in 0:5){
  print(k)
       expected[k+1]<- 25*dpois(k,uhat)
       observed[k+1]<- sum(mytab==k)
      }
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
round(expected,4)
## [1] 0.0095 0.0745 0.2936 0.7711 1.5191 2.3941 0.0000
observed
## [1] 0 0 2 2 1 1 0
expected[7] <- 25-sum(expected[1:6])
observed[7]<- sum(mytab>=6)
round(expected,4)
## [1]  0.0095  0.0745  0.2936  0.7711  1.5191  2.3941 19.9381
observed
## [1]  0  0  2  2  1  1 19
chisq<-sum(((expected-observed)^2)/expected)
chisq
## [1] 12.99428
pvalue<-pchisq(chisq,5,lower.tail=F)
pvalue
## [1] 0.02343241

The P value is 0.2 so we may reject H0 and this is not CSR.

Q5: Use simulation to roughly approximate the power of the d1 test

Create a function to calculates the d1 statistic

d.one <- function (a)
{
  nnd <- sort(nndist(a))
  dmin <- nnd[1]
  d.one <- a$n*(a$n-1)*pi*dmin^2
}

Then create a loop that simulates an SSI processes and returns the d1 statistic for each simulation; run 100 sumilations with an inhibition distance of .05 and calculate how many times the null hypothesis of CSR is rejected

dsim.ssi <- function(m,n){
  set.seed(4887)
  temp <- rep(0, m)
  for(i in 1:m)
    {
    mypatt <- rSSI(n, 80)
    d1 <- d.one(mypatt)
    temp[i] <- (1-pchisq(d1, 2))*2
    }
  return(temp)
}

simResults <- dsim.ssi(100, .05)
sum(simResults < .05)
## [1] 100

It seems that all 100 times the null hypothesis of CSR is rejected

Then if we change the inhibition distance t0 0.02, 0.015 and 0.01

simResults <- dsim.ssi(100, .02)
sum(simResults < .05)
## [1] 100
simResults <- dsim.ssi(100, .015)
sum(simResults < .05)
## [1] 23
simResults <- dsim.ssi(100, .01)
sum(simResults < .05)
## [1] 7

If we change to 0.02, still 100 rejected; but if we change inhibition distance to 0.015, only 23 rejection and if the inhibition distance is 0.01, only 7 rejection. It suggested that the power of the d1 test is low, the inhibition distance is very important for conclusion

Class Resources: Dr. Cristine Morgan’s Lab/Applied Spatial Analysis