library(spatstat)
The command rpoispp(100)
generates realisations of the
Poisson process with intensity \(\lambda =
100\) in the unit square.
Repeat the command plot(rpoispp(100))
several times
to build your intuition about the appearance of a completely random
pattern of points.
Let’s plot it three times:
replicate(3, plot(rpoispp(lambda = 100), main = ""))
As can be seen, the points (unsurprisingly) are much more random that want one might think. “Randomly” drawing points on a piece of paper one would usually draw a point pattern that is more regular (i.e. the points are repulsive).
Try the same thing with intensity \(\lambda = 1.5\).
For brevity we only do it once here:
plot(rpoispp(lambda = 1.5), main = "")
Here we expect 1.5 points in the plot each time.
Returning to the Japanese Pines data,
Fit the uniform Poisson point process model to the Japanese Pines data
ppm(japanesepines~1)
We fit the Poisson process model with the given command and print the output:
m.jp <- ppm(japanesepines ~ 1)
print(m.jp)
## Stationary Poisson process
## Intensity: 65
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## log(lambda) 4.174387 0.1240347 3.931284 4.417491 *** 33.65499
Read off the fitted intensity. Check that this is the correct value of the maximum likelihood estimate of the intensity.
We extract the coeficient with the coef
function, and
compare to the straightforward estimate obtained by `intensity``:
unname(exp(coef(m.jp)))
## [1] 65
intensity(japanesepines)
## [1] 65
As seen, they agree exactly.
The japanesepines
dataset is believed to exhibit spatial
inhomogeneity.
Plot a kernel smoothed intensity estimate.
Plot the kernel smoothed intensity estimate selecting the bandwidth
with bw.scott
:
jp.dens <- density(japanesepines, sigma = bw.scott)
plot(jp.dens)
plot(japanesepines, col = "white", cex = .4, pch = 16, add = TRUE)
Fit the Poisson point process models with loglinear intensity
(trend formula ~x+y
) and log-quadratic intensity (trend
formula ~polynom(x,y,2)
) to the Japanese Pines data.
We fit the two models with ppm
:
jp.m <- ppm(japanesepines ~ x + y)
jp.m2 <- ppm(japanesepines ~ polynom(x, y, 2) )
extract the fitted coefficients for these models using
coef
.
coef(jp.m)
## (Intercept) x y
## 4.0670790 -0.2349641 0.4296171
coef(jp.m2)
## (Intercept) x y I(x^2) I(x * y) I(y^2)
## 4.0645501 1.1436854 -1.5613621 -0.7490094 -1.2009245 2.5061569
Plot the fitted model intensity (using
plot(predict(fit))
)
par(mar=rep(0,4))
plot(predict(jp.m), main = "")
plot(predict(jp.m, se=TRUE)$se, main = "")
plot(predict(jp.m2), main = "")
plot(predict(jp.m2, se=TRUE)$se, main = "")
perform the Likelihood Ratio Test for the null hypothesis of a
loglinear intensity against the alternative of a log-quadratic
intensity, using anova
.
anova(jp.m, jp.m2)
## Analysis of Deviance Table
##
## Model 1: ~x + y Poisson
## Model 2: ~x + y + I(x^2) + I(x * y) + I(y^2) Poisson
## Npar Df Deviance
## 1 3
## 2 6 3 3.3851
Generate 10 simulated realisations of the fitted log-quadratic
model, and plot them, using plot(simulate(fit, nsim=10))
where fit
is the fitted model.
par(mar=rep(0.5,4))
plot(simulate(jp.m2, nsim=10), main = "")
## Generating 10 simulated patterns ...1, 2, 3, 4, 5, 6, 7, 8, 9, 10.