Exercises day 4

Exercises day 4

Try fitting homogeneous Poisson, inhomogeneous Poisson and clustered point process models to the hardcore pattern simulated on day 3.

library(spatstat)

## Loading required package: spatstat.data

## Loading required package: spatstat.univar

## spatstat.univar 3.1-1

## Loading required package: spatstat.geom

## spatstat.geom 3.3-5

## Loading required package: spatstat.random

## spatstat.random 3.3-2

## Loading required package: spatstat.explore

## Loading required package: nlme

## spatstat.explore 3.3-4

## Loading required package: spatstat.model

## Loading required package: rpart

## spatstat.model 3.3-4

## Loading required package: spatstat.linnet

## spatstat.linnet 3.2-5

## 
## spatstat 3.3-1 
## For an introduction to spatstat, type 'beginner'

# no seed set, so expect different results:
set.seed(1357)
ppi = rpoispp(function(x,y,...) 500 * x)
cl <- rThomas(100, .02, 5)
hc <- rHardcore(0.05, 1.5, square(50)) 
ppm(hc, ~1) # homogeneous Poisson

## Stationary Poisson process
## Fitted to point pattern dataset 'hc'
## Intensity: 0.0348
##              Estimate      S.E.   CI95.lo   CI95.hi Ztest      Zval
## log(lambda) -3.358138 0.1072113 -3.568268 -3.148008   *** -31.32263

ppm(hc, ~x) # inhomogeneous Poisson

## Nonstationary Poisson process
## Fitted to point pattern dataset 'hc'
## 
## Log intensity:  ~x
## 
## Fitted trend coefficients:
##  (Intercept)            x 
## -3.209065155 -0.006118179 
## 
##                 Estimate        S.E.     CI95.lo    CI95.hi Ztest        Zval
## (Intercept) -3.209065155 0.206680317 -3.61415113 -2.8039792   *** -15.5267091
## x           -0.006118179 0.007446249 -0.02071256  0.0084762        -0.8216458

kppm(hc)

## Stationary cluster point process model
## Fitted to point pattern dataset 'hc'
## Fitted by minimum contrast
## 	Summary statistic: K-function
## 
## Uniform intensity:	0.0348
## 
## Cluster model: Thomas process
## Fitted cluster parameters:
##    kappa    scale 
## 114.0657 376.2674 
## Mean cluster size:  0.0003050873 points
## 
## Cluster strength: phi =  4.928e-09
## Sibling probability: psib =  4.928e-09

and for the long output:

kppm(hc) |> summary()

## Stationary cluster point process model
## Fitted to point pattern dataset 'hc'
## Fitted by minimum contrast
## 	Summary statistic: K-function
## Minimum contrast fit (object of class "minconfit")
## Model: Thomas process
## Fitted by matching theoretical K function to hc
## 
## Internal parameters fitted by minimum contrast ($par):
##       kappa      sigma2 
##    114.0657 141577.1448 
## 
## Fitted cluster parameters:
##    kappa    scale 
## 114.0657 376.2674 
## Mean cluster size:  0.0003050873 points
## 
## Converged successfully after 61 function evaluations
## 
## Starting values of parameters:
##    kappa   sigma2 
##  0.03480 43.19339 
## Domain of integration: [ 0 , 12.5 ]
## Exponents: p= 2, q= 0.25
## 
## ----------- TREND  -----
## Point process model
## Fitted to data: X
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = X, trend = trend, rename.intercept = FALSE, covariates = covariates, 
##     covfunargs = covfunargs, use.gam = use.gam, forcefit = TRUE, 
##     improve.type = ppm.improve.type, improve.args = ppm.improve.args, 
##     nd = nd, eps = eps)
## Edge correction: "border"
## 	[border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
## 
## Data pattern:
## Planar point pattern:  87 points
## Average intensity 0.0348 points per square unit
## Window: rectangle = [0, 50] x [0, 50] units
## Window area = 2500 square units
## 
## Dummy quadrature points:
##      32 x 32 grid of dummy points, plus 4 corner points
##      dummy spacing: 1.5625 units
## 
## Original dummy parameters: =
## Planar point pattern:  1028 points
## Average intensity 0.411 points per square unit
## Window: rectangle = [0, 50] x [0, 50] units
## Window area = 2500 square units
## Quadrature weights:
##      (counting weights based on 32 x 32 array of rectangular tiles)
## All weights:
## 	range: [1.22, 2.44]	total: 2500
## Weights on data points:
## 	range: [1.22, 1.22]	total: 106
## Weights on dummy points:
## 	range: [1.22, 2.44]	total: 2390
## --------------------------------------------------------------------------------
## FITTED :
## 
## Stationary Poisson process
## 
## ---- Intensity: ----
## 
## 
## Uniform intensity:
## [1] 0.0348
## 
##              Estimate      S.E.   CI95.lo   CI95.hi Ztest      Zval
## (Intercept) -3.358138 0.1072113 -3.568268 -3.148008   *** -31.32263
## 
## ----------- gory details -----
## 
## Fitted regular parameters (theta):
## (Intercept) 
##   -3.358138 
## 
## Fitted exp(theta):
## (Intercept) 
##      0.0348 
## 
## ----------- CLUSTER  -----------
## Model: Thomas process
## 
## Fitted cluster parameters:
##    kappa    scale 
## 114.0657 376.2674 
## Mean cluster size:  0.0003050873 points
## 
## Final standard error and CI
## (allowing for correlation of cluster process):
##              Estimate      S.E.   CI95.lo   CI95.hi Ztest      Zval
## (Intercept) -3.358138 0.1072113 -3.568268 -3.148008   *** -31.32262
## 
## ----------- cluster strength indices ----------
## Sibling probability 4.927673e-09
## Count overdispersion index (on original window): 1
## Cluster strength: 4.927673e-09
## 
## Spatial persistence index (over window): 0.9912098
## 
## Bound on distance from Poisson process (over window): 3.729809e-05
## 	 = min (1, 0.0530771, 0.0265799, 3.729809e-05, 0.006107173)
## 
## Bound on distance from MIXED Poisson process (over window): 2.015028e-07
## 
## Intensity of parents of nonempty clusters: 0.03479469
## Mean number of offspring in a nonempty cluster: 1.000153
## Intensity of parents of clusters of more than one offspring point: 5.30744e-06
## Ratio of parents to parents-plus-offspring: 0.999695 (where 1 = Poisson 
## process)
## Probability that a typical point belongs to a nontrivial cluster: 0.0003050408
## 
## >>> The Poisson process is a better fit <<<

Using krige.cv:

carry out a leave-one-out cross validation using the four fitted variogram models of day 3, and compute the root mean square prediction error for all four models. Which one is favourable in this respect?

# create no2.sf:
library(sf)

## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE

no2 <- read.csv(system.file("external/no2.csv",
    package = "gstat"))
crs <- st_crs("EPSG:32632") # a csv doesn't carry a CRS!
st_as_sf(no2, crs = "OGC:CRS84", coords =
    c("station_longitude_deg", "station_latitude_deg")) |>
    st_transform(crs) -> no2.sf
# sample variogram:
library(gstat)

## 
## Attaching package: 'gstat'

## The following object is masked from 'package:spatstat.explore':
## 
##     idw

v = variogram(NO2~1, no2.sf)
fit.e = fit.variogram(v, vgm(1, "Exp", 50000))
fit.s = fit.variogram(v, vgm(1, "Sph", 50000))
fit.m1 = fit.variogram(v, vgm(1, "Mat", 50000, kappa = .3))
fit.m2 = fit.variogram(v, vgm(1, "Mat", 50000, kappa = 4))
rmse = function(x) sqrt(mean(x$residual^2))
krige.cv(NO2~1, no2.sf, fit.e) |> rmse()

## [1] 3.870342

krige.cv(NO2~1, no2.sf, fit.s) |> rmse()

## [1] 4.203151

krige.cv(NO2~1, no2.sf, fit.m1) |> rmse()

## [1] 3.745997

krige.cv(NO2~1, no2.sf, fit.m2) |> rmse()

## [1] 4.079682

What causes the differences between the mean and the variance of the conditional simulations (left) and the mean and variance obtained by kriging (right)?

The sample size: for a larger sample size (larger number of simulations) they become closer.

"https://github.com/edzer/sdsr/raw/main/data/de_nuts1.gpkg" |>
  read_sf() |>
  st_transform(crs) -> de
library(stars)

## Loading required package: abind

g2 = st_as_stars(st_bbox(de))
g3 = st_crop(g2, de)
v.fit = fit.e
cs = krige(NO2~1, no2.sf, g3, v.fit, nsim = 500, nmax = 30)

## drawing 500 GLS realisations of beta...
## [using conditional Gaussian simulation]

csm = st_apply(cs, c("x", "y"), mean)
csm$kriging = krige(NO2~1, no2.sf, g3, v.fit)[1]

## [using ordinary kriging]

plot(merge(csm), breaks = "equal")

csv = st_apply(cs, c("x", "y"), var)
csv$kr_var = krige(NO2~1, no2.sf, g3, v.fit)[2]

## [using ordinary kriging]

plot(merge(csv), breaks = "equal")

When comparing (above) the sample variogram …

of simulated fields with the variogram model used to simulate them, where do you see differences? Can you explain why these are not identical, or hypothesize under which circumstances these would become more similar (or identical)?

• differences are larger for larger lag values
• if we would make the extent (area) larger, we get more independent information for any given distance (lag), and obtain spatial correlations values that are closer to the ones used for generating the Gaussian process.

Under which practical data analysis problem would conditional simulations be more useful than the kriging prediction and kriging variance maps?

Whenever the maps are input to a model that is sensitive to errors in predictions and/or to correlations of prediction errors, e.g. a model that does a non-linear aggregation over space and time, such as a rainfall-runoff model.