6  Spatial and spatiotemporal log-Gaussian Cox processes

library(stelfi)
library(ggplot2)

6.1 The model and your data

A log-Gaussian Cox process (LGCP) is a model for the locations of events in a study region. Unlike a Hawkes process, events do not cause other events. Instead, the intensity (the expected number of events per unit area) varies smoothly over space because of some unobserved environmental “surface”: events are more common where the surface is high and rarer where it is low. In the simplest case the intensity at location \(s\) is

\[\lambda(s) = \text{exp}(\beta_0 + G(s)), \tag{6.1}\]

where \(\beta_0\) is the intercept, the log intensity at an “average” location, and \(G(s)\) is a Gaussian random field, a smooth random surface with mean zero. Because \(G(s)\) is estimated from the data, the model can capture spatial clustering of events without you having to explain where it comes from. Covariates (e.g., population density) can be added to the log intensity, \(\beta_0 + \beta_1 x_1(s) + ... + G(s)\), and the random field then captures the spatial variation that the covariates do not explain.

Plotted below is a realisation of an LGCP within a disc shaped region, overlain on the random field \(G(s)\): events are more common where the field is high (yellow).

sf <- sf::st_as_sf(spatstat.geom::disc())
smesh <- fmesher::fm_mesh_2d(loc.domain = sf::st_coordinates(sf)[, 1:2], max.edge = 0.2, cutoff = 0.1)
set.seed(1)
sim <- sim_lgcp(parameters = c(beta = 3, log_tau = 1.6, log_kappa = 1.95), sf = sf, smesh = smesh)
pts <- smesh$loc[rep(1:nrow(smesh$loc), times = sim$y), 1:2]
set.seed(23)
pts <- data.frame(x = pts[, 1] + runif(nrow(pts), -0.05, 0.05), y = pts[, 2] + runif(nrow(pts), -0.05, 0.05))
pts <- pts[spatstat.geom::inside.owin(pts, w = spatstat.geom::disc()), ]
show_field(sim$x, smesh, sf = sf, clip = TRUE) + theme_void() + theme(legend.position = "none") +
    geom_point(data = pts, aes(x = x, y = y), size = 2)

The random field is described by two parameters, which stelfi reports on an interpretable scale:

• the range \(r\), the distance at which the correlation between values of the field drops to about 0.1. A small range means the intensity varies over short distances (small “hot spots”), a large range means it varies smoothly over large distances; and
• the standard deviation \(\sigma\) of the field, how much the log intensity varies from place to place. For example \(\sigma = 1\) means that, on the log scale, the intensity typically varies by \(\pm 1\) from \(\beta_0\), i.e., by a factor of about \(e^{1} \approx 2.7\).

Internally the field is parameterised by \(\tau\) and \(\kappa\), so get_coefs() reports log_tau and log_kappa (the parameters that are actually estimated) along with range \(= \sqrt{8}/\kappa\) and stdev \(= 1/\sqrt{4\pi\kappa^2\tau^2}\). Following Lindgren, Rue, and Lindström (2011), the field has a Matérn covariance function with smoothness \(\nu = 1\), and it is represented using the stochastic partial differential equation (SPDE) approach on a triangulation (a mesh) of the study region.

The mesh controls how finely the field can vary in space and how much computation is needed, so it is an important part of fitting any spatial model in stelfi. How to build a mesh, and how to check whether it is adequate, is covered in Chapter 5.

6.1.1 What your data need to look like

fit_lgcp() needs three things:

• locs, a data frame of event locations with columns named x and y (and t for a spatiotemporal model, Section 6.5);
• sf, an sf polygon of the study region; and
• smesh, a mesh (Chapter 5).

Things to be aware of:

• Use a projected coordinate system. The field is defined in terms of distances, so longitude and latitude (in degrees) should be projected to a planar coordinate system (e.g., UTM) using sf::st_transform(). Working in kilometres rather than metres keeps the numbers well scaled and makes the parameters easy to interpret: the intensity is then per km\(^2\) and the range is in km. The events, the study region and the mesh must all be in the same coordinates (it is often easiest to drop the CRS once everything is projected, see Section 6.3.2).
• The study region should be where events could have been observed. The model assumes that every part of the study region was observed with equal effort. If, for example, sightings can only be made where there are people, then the intensity estimates the rate of reported events, and the random field will pick up variation in observation effort as well as variation in the underlying process. Including a covariate that describes effort (e.g., population density) can help. Events outside the study region are ignored.
• Events are summarised by mesh node. stelfi counts the number of events closest to each mesh node (i.e., in each cell of the “dual” mesh) and models these counts. Detail at scales smaller than the mesh triangles is therefore lost, which is another reason to check that results are not sensitive to the mesh.
• Duplicate locations are allowed for fitting (unlike tied times in a Hawkes process), but they need to be jittered for some of the goodness-of-fit checks (Section 6.4).
• Covariates must be given at the mesh nodes, as a matrix with one row per mesh node (smesh$n rows) and one column per covariate. If your covariate is a raster or is only known at some locations, it needs to be extracted or interpolated at the mesh nodes, see Section 6.3.2.
• Events do not interact. An LGCP describes clustering caused by a smooth environment, not clustering caused by events triggering each other. If the latter is plausible, consider a spatiotemporal self-exciting model (Chapter 8).
• Number of events. The field parameters (in particular the range) are hard to estimate with few events or when the range is large relative to the study region. A few hundred events is a reasonable minimum.

6.2 A simulated example

Data were simulated from an LGCP (using spatstat, Baddeley, Rubak, and Turner (2015)) within a disc of radius 1, with intercept \(\beta_0 = 4\), field standard deviation \(\sigma = 1\), and range \(r = \sqrt{8} \times 0.2 \approx\) 0.566. In spatstat.random::rLGCP() the field has a Matérn covariance with smoothness nu \(= 1\) (the same as in stelfi), mu is the intercept, var is the variance of the field (\(\sigma^2\)), and scale is \(1/\kappa\), so that the range is \(\sqrt{8} \times\) scale.

win <- spatstat.geom::disc()
set.seed(1234)
sim <- spatstat.random::rLGCP("matern", mu = 4, var = 1, scale = 0.2, nu = 1, win = win)
locs <- data.frame(x = sim$x, y = sim$y)
truth <- c(beta = 4, range = sqrt(8) * 0.2, sigma = 1)

The data are a data frame of 272 event locations:

head(locs)

           x          y
1 -0.1800667  0.8802885
2  0.5430666 -0.5046758
3  0.2112133  0.8551978
4  0.1208994  0.9408004
5  0.2303162 -0.8321345
6 -0.4731538  0.5615277

The study region is a disc of radius 1, which is converted to an sf polygon. The plot below shows the events and the true (simulated) random field.

sf <- sf::st_sf(geometry = sf::st_geometry(sf::st_as_sf(spatstat.geom::disc())))

xs <- attr(sim, "Lambda")$xcol; ys <- attr(sim, "Lambda")$yrow
pxl <- expand.grid(x = xs, y = ys)
pxl$rf <- log(as.vector(t(attr(sim, "Lambda")$v))) - 4 ## minus the intercept
ggplot() +
    geom_tile(data = pxl, aes(x = x, y = y, fill = rf)) +
    labs(fill = "G(s)") + xlab("") + ylab("") +
    scale_fill_viridis_c(option = "D", na.value = NA) +
    coord_equal() + theme_void() +
    geom_sf(data = sf, fill = NA, linewidth = 1) +
    geom_point(data = locs, aes(x = x, y = y), size = 0.8)

6.2.1 Fitting the model using stelfi

First a mesh is built around the events (Chapter 5). Here max.edge is 0.1 inside the region, about a fifth of the true range.

smesh <- fmesher::fm_mesh_2d(loc = locs[, 1:2], max.edge = c(0.1, 0.3), cutoff = 0.03)
smesh$n ## number of mesh nodes

[1] 1124

The model is fitted using fit_lgcp() with starting values for \(\beta_0\) (beta), \(\log(\tau)\) (log_tau) and \(\log(\kappa)\) (log_kappa). A sensible starting value for the intercept is the log of the average number of events per unit area.

area <- as.numeric(sf::st_area(sf))
fit <- fit_lgcp(locs = locs, sf = sf, smesh = smesh,
                parameters = c(beta = log(nrow(locs) / area), log_tau = log(1), log_kappa = log(1)))
stelfi <- get_coefs(fit)
stelfi

            Estimate Std. Error
beta       4.1456922  0.1798657
log_tau   -3.2271510  0.2588465
log_kappa  2.2036142  0.2577273
range      0.3122680  0.0804800
stdev      0.7850755  0.1059713

The estimated random field is extracted using get_fields() and plotted using show_field(). Compare this with the true field above.

get_fields(fit, smesh) |>
    show_field(smesh = smesh, sf = sf, clip = TRUE) + theme_void() +
    geom_point(data = locs, aes(x = x, y = y), size = 0.8)

6.2.2 Comparison with inlabru and spatstat

The same model is fitted below using inlabru (Bachl et al. (2019)) and spatstat (Baddeley, Rubak, and Turner (2015)). When comparing the results bear in mind that:

• stelfi uses maximum likelihood (the field is integrated out using the Laplace approximation) so there are no priors and the estimates have standard errors.
• inlabru is Bayesian (using INLA), so priors are needed for the field and it reports posterior means and standard deviations. Penalised complexity priors are used via inla.spde2.pcmatern(): prior.range = c(0.05, 0.01) means \(P(r < 0.05) = 0.01\) and prior.sigma = c(5, 0.01) means \(P(\sigma > 5) = 0.01\), both deliberately weak. To make inlabru as close as possible to a maximum likelihood fit, the hyperparameters are not integrated over (int.strategy = "eb"). The same mesh is used for both.
• spatstat fits the model by minimum contrast, matching the theoretical \(K\) function of the LGCP to the empirical one. This does not use a mesh and does not give standard errors. Its Matérn scale parameter (reported as alpha in modelpar) is \(1/\kappa\), so the range is \(\sqrt{8} \times\) scale, and it reports the variance sigma2 rather than the standard deviation.

Note that inlabru requires the geometry column of the study region to be called geometry, which is why sf was created using sf::st_sf(geometry = ...) above.

library(inlabru)
library(INLA)
sf_locs <- sf::st_as_sf(locs, coords = c("x", "y"))
matern <- INLA::inla.spde2.pcmatern(smesh, prior.range = c(0.05, 0.01), prior.sigma = c(5, 0.01))
cmp <- geometry ~ field(geometry, model = matern) + Intercept(1)
fit_bru <- lgcp(cmp, sf_locs, samplers = sf, domain = list(geometry = smesh),
                options = bru_options(control.inla = list(int.strategy = "eb")))
bru <- rbind(fit_bru$summary.fixed[, 1:2], fit_bru$summary.hyperpar[, 1:2])
bru

                     mean        sd
Intercept       4.1897729 0.1942438
Range for field 0.4284088 0.1123945
Stdev for field 0.7812107 0.1223857

ss <- spatstat.model::lgcp.estK(sim, covmodel = list(model = "matern", nu = 1))
ss$modelpar

   sigma2     alpha        mu 
0.6509685 0.1252826 4.1359896 

fmt <- function(x) paste0(round(x[1], 3), " (", round(x[2], 3), ")")
params <- data.frame(Truth = round(truth, 3),
                     stelfi = apply(stelfi[c("beta", "range", "stdev"), ], 1, fmt),
                     inlabru = apply(bru, 1, fmt),
                     spatstat = round(c(ss$modelpar["mu"], sqrt(8) * ss$modelpar["alpha"],
                                        sqrt(ss$modelpar["sigma2"])), 3))
knitr::kable(t(params), col.names = c("$\\beta_0$", "$r$", "$\\sigma$"), escape = FALSE)

	\(\beta_0\)	\(r\)	\(\sigma\)
Truth	4.000	0.566	1.000
stelfi	4.146 (0.18)	0.312 (0.08)	0.785 (0.106)
inlabru	4.19 (0.194)	0.428 (0.112)	0.781 (0.122)
spatstat	4.136	0.354	0.807

stelfi and inlabru give very similar estimates, as they fit the same model on the same mesh. spatstat uses a very different method but tells a similar story. The intercept is well estimated by all three, while the range and the standard deviation of the field are much harder to pin down from a single dataset (note the larger standard errors).

6.3 Applied examples

6.3.1 Sasquatch sightings in the USA

The sasquatch data (see Chapter 10) are the locations of 972 reported sasquatch (bigfoot) sightings in the contiguous USA between 2000 and 2005, collated by the Bigfoot Field Researchers Organization. This example is also used in Jones-Todd and van Helsdingen (2024).

Preparing the data

The sightings are given as longitude/latitude. They, and the boundary of the contiguous USA (from the maps package), are projected to the Conus Albers projection (EPSG:5070) and converted from metres to kilometres.

data("sasquatch", package = "stelfi")
sf::sf_use_s2(FALSE)
## boundary of the contiguous USA: projected and in km
us <- maps::map("usa", fill = TRUE, plot = FALSE) |>
    sf::st_as_sf() |>
    sf::st_make_valid() |>
    sf::st_transform(crs = 5070)
us_km <- sf::st_sf(geometry = sf::st_union(sf::st_geometry(us)) / 1000)
## sighting locations: projected and in km
locs_sq <- sf::st_transform(sasquatch, crs = 5070) |>
    sf::st_coordinates() |>
    as.data.frame()
names(locs_sq) <- c("x", "y")
locs_sq <- locs_sq / 1000
head(locs_sq)

          x        y
1 -515.6556 2280.102
2 1532.8509 2407.822
3 -398.2722 2817.326
4 -505.9007 2766.177
5 1734.8802 2163.947
6 -832.7920 1547.099

Fitting the model

The mesh has a maximum edge length of 200 km (the contiguous USA is about 4500 km across).

smesh_km <- fmesher::fm_mesh_2d(loc = locs_sq[, 1:2], max.edge = c(200, 500), cutoff = 100)
area_sq <- as.numeric(sf::st_area(us_km))
fit_sq <- fit_lgcp(locs = locs_sq, sf = us_km, smesh = smesh_km,
                   parameters = c(beta = log(nrow(locs_sq) / area_sq),
                                  log_tau = log(1), log_kappa = log(1)))
coefs <- get_coefs(fit_sq)
coefs

            Estimate  Std. Error
beta       -9.937128   0.3865470
log_tau     4.046748   0.1069927
log_kappa  -5.662232   0.1818064
range     813.993537 147.9892483
stdev       1.419028   0.1880643

From these estimates, the range of the field is around 814 km, so sightings are clustered at a regional scale, and the standard deviation of 1.42 means the intensity of sightings varies by a large factor between regions. The estimated random field and the estimated intensity (expected number of sightings per km\(^2\)) are plotted using show_field() and show_lambda().

get_fields(fit_sq, smesh_km) |>
    show_field(smesh = smesh_km, sf = us_km, clip = TRUE) + theme_void()

show_lambda(fit_sq, smesh = smesh_km, sf = us_km, clip = TRUE) + theme_void()

The hot spots (e.g., the Pacific Northwest) are clear. Remember, though, that these are reported sightings: the model cannot distinguish places where sasquatch are common from places where many people are out looking for them (Section 6.1.1).

6.3.2 Terrorism events in Iraq

The iraq_terrorism data (see Chapter 10) are the locations and dates of terrorism events in Iraq carried out by ISIL between 2013 and 2017. This example shows the full stelfi workflow for a region with an awkward boundary, compares stelfi with inlabru on several meshes, explores how sensitive the results are to the mesh, and includes covariates in the intensity. A spatiotemporal version is in Section 6.5.3.

Preparing the data

The event locations are projected to UTM zone 38N (EPSG:32638), which covers most of Iraq, and converted to kilometres. The boundary of Iraq is obtained from the rnaturalearth package and transformed in the same way.

data("iraq_terrorism", package = "stelfi")
## event locations (lon/lat) -> UTM zone 38N -> km
locs_km <- as.data.frame(sf::st_coordinates(sf::st_transform(iraq_terrorism, crs = 32638)) / 1000)
names(locs_km) <- c("x", "y")
## the boundary of Iraq in the same coordinates
iraq_sf <- rnaturalearth::ne_countries(country = "Iraq", scale = "medium", returnclass = "sf") |>
    sf::st_transform(crs = 32638) |>
    sf::st_union()                     ## collapse to a single polygon
iraq_sf_km <- sf::st_sf(geometry = sf::st_geometry(iraq_sf) / 1000)
## once in km the coordinates are just planar numbers, so drop the CRS
## so that the events, the boundary and the mesh are all treated consistently
sf::st_crs(iraq_sf_km) <- NA
points_km <- sf::st_as_sf(locs_km, coords = c("x", "y")) ## used by inlabru

Building the meshes

To show how much the mesh matters, three very different meshes are built. Each follows the boundary of Iraq (boundary), is placed around the events, and has a finer inner and coarser outer part (max.edge = c(inner, outer)).

• Too coarse (max.edge = c(200, 400), cutoff = 100): triangles a couple of hundred km across, larger than many of the clusters of attacks;
• Adequate (max.edge = c(40, 150), cutoff = 20): triangles a few tens of km across. The relatively large cutoff means that events in the same town share a node, so there are no very small triangles; and
• Too fine (max.edge = c(10, 100), cutoff = 1): there is a node at almost every distinct event location, so there are many tiny triangles where attacks are clustered.

mesh_args <- list(too_coarse = list(max.edge = c(200, 400), cutoff = 100),
                  adequate   = list(max.edge = c(40, 150), cutoff = 20),
                  too_fine   = list(max.edge = c(10, 100), cutoff = 1))
meshes <- lapply(mesh_args, function(a) {
    fmesher::fm_mesh_2d(loc = locs_km[, 1:2],
                        boundary = fmesher::fm_as_segm(iraq_sf_km),
                        max.edge = a$max.edge, cutoff = a$cutoff)
})
## number of nodes, and the size of the triangles (km^2)
t(sapply(meshes, function(m) {
    area <- as.numeric(sf::st_area(mesh_2_sf(m)))
    c(nodes = m$n, triangles = length(area), median_area = round(median(area)),
      smallest_area = round(min(area), 2), below_10km2 = sum(area < 10))
}))

           nodes triangles median_area smallest_area below_10km2
too_coarse    65       106        8431       4128.56           0
adequate     851      1662         338        164.21           0
too_fine   12186     24309          19          0.00         970

(The triangle areas are calculated directly from mesh_2_sf(), since meshmetrics() (Section 5.3) is slow for a mesh with tens of thousands of triangles.) The three meshes, the events, and the boundary of Iraq (blue) are shown below.

patchwork::wrap_plots(lapply(names(meshes), function(m) {
    ggplot() +
        geom_sf(data = mesh_2_sf(meshes[[m]]), fill = NA, linewidth = 0.1) +
        geom_sf(data = iraq_sf_km, fill = NA, col = "blue", linewidth = 0.8) +
        geom_sf(data = points_km, size = 0.2, col = "red") +
        theme_void() + ggtitle(m)
}), nrow = 1)

Fitting with stelfi

The model is fitted using each mesh, and the time taken (in seconds) is recorded using system.time().

area_km2 <- as.numeric(sf::st_area(iraq_sf_km))
beta_start <- log(nrow(locs_km) / area_km2)
start <- c(beta = beta_start, log_tau = log(1), log_kappa = log(1))
fits_stelfi <- list(); time_stelfi <- c()
for (m in names(meshes)) {
    time_stelfi[m] <- system.time(
        fits_stelfi[[m]] <- fit_lgcp(locs = locs_km, sf = iraq_sf_km, smesh = meshes[[m]],
                                     parameters = start)
    )["elapsed"]
}
get_coefs(fits_stelfi$adequate)

            Estimate Std. Error
beta      -9.9442475 0.52076890
log_tau    0.9717464 0.09948255
log_kappa -3.3800065 0.12009682
range     83.0736230 9.97687820
stdev      3.1353723 0.24396259

The estimated random field and intensity surface (events per km\(^2\)) for the adequate mesh:

get_fields(fits_stelfi$adequate, meshes$adequate) |>
    show_field(smesh = meshes$adequate, sf = iraq_sf_km, clip = TRUE) + theme_void()

show_lambda(fits_stelfi$adequate, smesh = meshes$adequate, sf = iraq_sf_km, clip = TRUE) +
    theme_void()

Comparing with inlabru

The same model is fitted using inlabru on the same three meshes, with weak penalised complexity priors (\(P(r < 5\text{ km}) = 0.01\) and \(P(\sigma > 10) = 0.01\)) and empirical Bayes (see Section 6.2), again recording the time taken.

fit_bru <- function(mesh) {
    matern <- INLA::inla.spde2.pcmatern(mesh, prior.sigma = c(10, 0.01), prior.range = c(5, 0.01))
    cmp <- geometry ~ mySmooth(geometry, model = matern) + Intercept(1)
    inlabru::lgcp(cmp, data = points_km, samplers = iraq_sf_km,
                  domain = list(geometry = mesh),
                  options = inlabru::bru_options(
                      control.inla = list(int.strategy = "eb", strategy = "gaussian",
                                          diagonal = 1e-4, control.vb = list(enable = FALSE))))
}
fits_inlabru <- list(); time_inlabru <- c()
for (m in names(meshes)) {
    time_inlabru[m] <- system.time(fits_inlabru[[m]] <- fit_bru(meshes[[m]]))["elapsed"]
}

Sensitivity to the mesh

The table gives, for each mesh and package, the number of mesh nodes, the time taken to fit the model (in seconds; this will vary from computer to computer), and the estimates (standard error for stelfi, posterior standard deviation for inlabru, in brackets). Because the coordinates are in km, \(\beta_0\) is the log intensity per km\(^2\) and \(r\) is in km.

rows <- lapply(names(meshes), function(m) {
    s <- get_coefs(fits_stelfi[[m]])[c("beta", "range", "stdev"), 1:2]
    b <- rbind(fits_inlabru[[m]]$summary.fixed[, 1:2], fits_inlabru[[m]]$summary.hyperpar[, 1:2])
    data.frame(mesh = m, nodes = meshes[[m]]$n, package = c("stelfi", "inlabru"),
               time = round(c(time_stelfi[m], time_inlabru[m]), 1),
               beta0 = c(fmt(s[1, ]), fmt(b[1, ])),
               range = c(fmt(s[2, ]), fmt(b[2, ])),
               sigma = c(fmt(s[3, ]), fmt(b[3, ])))
})
knitr::kable(do.call(rbind, rows), row.names = FALSE,
             col.names = c("Mesh", "Nodes", "Package", "Time (s)", "$\\beta_0$", "$r$", "$\\sigma$"),
             escape = FALSE)

Mesh	Nodes	Package	Time (s)	\(\beta_0\)	\(r\)	\(\sigma\)
too_coarse	65	stelfi	0.1	-8.081 (1.161)	362.387 (112.719)	2.381 (0.462)
too_coarse	65	inlabru	1.1	-7.238 (1.348)	616.884 (238.843)	2.275 (0.591)
adequate	851	stelfi	0.7	-9.944 (0.521)	83.074 (9.977)	3.135 (0.244)
adequate	851	inlabru	2.1	-8.073 (0.604)	135.057 (18.16)	3.01 (0.323)
too_fine	12186	stelfi	30.9	-11.493 (0.2)	8.462 (0.509)	5.442 (0.244)
too_fine	12186	inlabru	17.2	-8.657 (0.216)	22.435 (1.27)	4.21 (0.19)

The intensity surfaces estimated by stelfi using the three meshes are shown below.

lambdas <- lapply(names(meshes), function(m) {
    show_lambda(fits_stelfi[[m]], smesh = meshes[[m]], sf = iraq_sf_km, clip = TRUE) +
        theme_void() + ggtitle(m) + theme(legend.position = "none")
})
patchwork::wrap_plots(lambdas, nrow = 1)

• Too coarse. With only a few dozen nodes both fits are almost instant, but the triangles are larger than the clusters of attacks. The field can only vary smoothly over hundreds of km, so the estimated range is several hundred km, the hot spots are smeared out, and the intensity surface is blocky.
• Too fine. With over ten thousand nodes both fits take far longer (here stelfi was the slower of the two for this mesh: TMB runs on a single thread by default, whereas INLA uses several). The field now resolves individual towns: the estimated range drops to a few km and the standard deviation increases, and the two packages disagree much more, since their approximations of the integral of the intensity over Iraq (stelfi assigns each event and each piece of area to its nearest mesh node, inlabru uses integration points) differ most when the mesh is this fine. Very little is gained for the extra computation: the large scale pattern of the intensity is the same as for the adequate mesh.
• Adequate. The fits take a second or two, the intensity surface shows the hot spots around Baghdad and along the Tigris and Euphrates, and the two packages tell a similar story: a field with a large standard deviation (around 3, so the intensity varies by orders of magnitude across Iraq) and a range of around 100 km (inlabru’s is somewhat larger). The intercepts differ more between the packages, partly because of the different integration schemes and partly because the intercept is only weakly identified when the field is this variable.

Note that the estimated range does not settle down as the mesh is refined, for either package. The reason is the data: attacks are concentrated in particular towns and cities, so there is spatial structure at every scale down to a few kilometres, and whatever mesh is used the field picks up variation down to the smallest scale the mesh can represent. This also means that the rule of thumb of Section 5.2.1 (inner max.edge no more than about a fifth of the range) cannot be satisfied. In this situation:

• the range should be interpreted as “the scale of variation resolved at this mesh”, and the mesh reported alongside the results;
• choose the mesh to match the scale of the question being asked (e.g., differences between regions versus between neighbourhoods), avoiding many small triangles that add computation without changing the large scale pattern; and
• consider whether the fine scale clustering is better explained by covariates (e.g., the locations of towns, population density) or by events triggering further events (Chapter 8).

Including covariates

Covariates are included in the log intensity, \(\lambda(s) = \text{exp}(\beta_0 + \beta_1 x_1(s) + ... + G(s))\), using the covariates argument of fit_lgcp(). The covariates must be a matrix with one column per covariate and one row per mesh node, and the starting value beta becomes a vector: the intercept followed by one coefficient per covariate. It is a good idea to standardise covariates (using scale()), so that the coefficients are the change in log intensity per standard deviation increase in the covariate and the optimiser works with parameters on a similar scale. The adequate mesh is used below.

smesh <- meshes$adequate

Distance to Baghdad. Here the covariate is the distance (in km) from each mesh node to Baghdad (33.30°N, 44.40°E).

baghdad <- sf::st_as_sf(data.frame(lon = 44.40, lat = 33.30), coords = c("lon", "lat"), crs = 4326) |>
    sf::st_transform(crs = 32638)
baghdad_km <- sf::st_coordinates(baghdad) / 1000
mesh_coords <- smesh$loc[, 1:2]
baghdad_dist <- sqrt((mesh_coords[, 1] - baghdad_km[1])^2 + (mesh_coords[, 2] - baghdad_km[2])^2)
baghdad_scaled <- as.matrix(scale(baghdad_dist))
sd(baghdad_dist) ## km per unit of the scaled covariate

[1] 134.3086

fit_baghdad <- fit_lgcp(locs = locs_km, sf = iraq_sf_km, smesh = smesh,
                        parameters = list(beta = c(beta_start, 0), log_tau = log(1), log_kappa = log(1)),
                        covariates = baghdad_scaled)
get_coefs(fit_baghdad)

            Estimate Std. Error
beta      -9.9469657  0.4618613
beta      -1.2816661  0.4223855
log_tau    0.9017762  0.1091130
log_kappa -3.2594699  0.1217000
range     73.6401720  8.9620054
stdev      2.9807699  0.2164733

The second beta is the change in the log intensity per standard deviation (see above, in km) increase in distance from Baghdad, so \(\text{exp}(\hat{\beta}_1)\) is the multiplicative change in the intensity. A negative coefficient means events are more intense closer to Baghdad. Note that the mesh extends beyond the boundary of Iraq, so the standardisation uses nodes outside the country too.

Population density. In these data population density is only recorded at the event locations, not at the mesh nodes, so it has to be interpolated to the mesh nodes. Here inverse distance weighting (gstat::idw()) is used.

points_pop <- points_km
points_pop$popdensity <- iraq_terrorism$popdensity
mesh_sf <- sf::st_as_sf(as.data.frame(smesh$loc[, 1:2]) |> setNames(c("x", "y")), coords = c("x", "y"))
idw_fit <- gstat::idw(formula = popdensity ~ 1, locations = points_pop,
                      newdata = mesh_sf, debug.level = 0)
scaled_popdensity <- as.matrix(scale(idw_fit$var1.pred))
fit_pop <- fit_lgcp(locs = locs_km, sf = iraq_sf_km, smesh = smesh,
                    parameters = list(beta = c(beta_start, 0), log_tau = log(1), log_kappa = log(1)),
                    covariates = scaled_popdensity)
get_coefs(fit_pop)

             Estimate Std. Error
beta      -10.0145140  0.4314508
beta       -1.7120502  0.1886892
log_tau     0.9140417  0.1043827
log_kappa  -3.2184006  0.1143297
range      70.6770795  8.0804916
stdev       2.8259560  0.2037467

Note

Care is needed when interpreting this coefficient. Because population density is only known at the events, the interpolated surface is informed only by where events happened, not by where people live in general, so it is not an independent description of the population. A population density raster covering the whole of Iraq would be preferable. It is also worth remembering that the random field can absorb spatial structure shared with a covariate (spatial confounding), so check the standard error of the coefficient, and compare the model with and without the covariate (e.g., by looking at how much the estimated field changes).

6.4 Goodness-of-fit

A natural question is whether point patterns simulated from the fitted model look like the observed pattern. This can be checked by comparing summary functions of the observed pattern with those of patterns simulated from the fitted model, using spatstat (Baddeley, Rubak, and Turner (2015)). Two nearest-neighbour summaries are used here:

• the \(G\) function, the distribution function of the distance from a typical event to its nearest other event. For a clustered pattern \(G(r)\) rises faster than for a completely random (Poisson) pattern since events tend to have a close neighbour; and
• the \(F\) function (empty space function), the distribution function of the distance from an arbitrary location to the nearest event. For a clustered pattern \(F(r)\) rises more slowly than for a Poisson pattern since there are large empty gaps between clusters.

If the fitted LGCP is adequate, the observed curves should lie within the envelope formed by the curves of the simulated patterns. Distances must be measured in planar coordinates, so this is done for the sasquatch model (Section 6.3.1), which was fitted in km.

6.4.1 Preparing the observed pattern

The observed events are converted to a spatstat point pattern (ppp) within the study region (owin). Some sightings share exactly the same location, which gives nearest-neighbour distances of zero, so the pattern is jittered slightly using rjitter().

library(spatstat.geom)
library(spatstat.explore)
us_window <- as.owin(us_km)
set.seed(1)
obs_ppp <- ppp(x = locs_sq$x, y = locs_sq$y, window = us_window)
obs_ppp <- rjitter(obs_ppp, retry = TRUE)

6.4.2 Simulating from the fitted model

The function below simulates one point pattern from the fitted LGCP using sim_lgcp() and the estimated parameters. This simulates a new random field each time (the model is checked as a whole, not just the estimated field). sim_lgcp() returns the number of events per mesh node, so each node is repeated as many times as it has events and the resulting pattern is jittered (by half the mesh cutoff) so that events are spread around the node.

generate_simulated_ppp <- function() {
    sim <- sim_lgcp(sf = us_km, smesh = smesh_km,
                    parameters = c(beta = unname(coefs["beta", 1]),
                                   log_tau = unname(coefs["log_tau", 1]),
                                   log_kappa = unname(coefs["log_kappa", 1])))
    if (sum(sim$y) < 2) return(NULL) ## need two events for nearest-neighbour distances
    idx <- rep(1:nrow(smesh_km$loc), times = sim$y)
    sim_ppp <- ppp(x = smesh_km$loc[idx, 1], y = smesh_km$loc[idx, 2],
                   window = us_window, check = FALSE)
    rjitter(sim_ppp, radius = 50, retry = TRUE)
}

6.4.3 Simulation envelopes

A simulation envelope gives, for each distance \(r\), the range (here the minimum and maximum) of the summary function over many simulated patterns. With 39 simulations a pointwise envelope based on the minimum and maximum corresponds to a test with significance level \(2/(39 + 1) = 0.05\) at each \(r\) (more simulations, e.g., 99, are preferable but take longer). Passing the list of LGCP simulations to the simulate argument of envelope() is important, otherwise the envelope is formed from completely random patterns and instead tests whether there is any spatial structure at all.

set.seed(3)
sim_list <- lapply(1:39, function(i) generate_simulated_ppp())
sim_list <- Filter(Negate(is.null), sim_list)
g_envelope <- envelope(Y = obs_ppp, fun = Gest, simulate = sim_list,
                       nsim = length(sim_list), savefuns = TRUE, verbose = FALSE)
plot(g_envelope, main = "Sasquatch LGCP: G function", xlab = "Distance (km)", legend = FALSE)

f_envelope <- envelope(Y = obs_ppp, fun = Fest, simulate = sim_list,
                       nsim = length(sim_list), savefuns = TRUE, verbose = FALSE)
plot(f_envelope, main = "Sasquatch LGCP: F function", xlab = "Distance (km)", legend = FALSE)

In each plot the solid black line is the observed function, the dashed red line is the mean of the simulated patterns, and the grey band is the envelope. Where the observed curve leaves the band the fitted model is not reproducing that aspect of the data. Points to bear in mind when reading these plots:

• at distances smaller than the mesh resolution (here about 100 to 200 km) the simulated patterns are only as realistic as the jittering used to spread events around their mesh node, so differences at very short distances say more about the mesh than about the model; and
• the \(G\) function is about how close events are to each other. If events are much closer together than the model expects (more short nearest-neighbour distances), this suggests clustering at a finer scale than the field can capture, e.g., because events trigger one another (Chapter 8) or because of duplicated reports.

6.5 Spatiotemporal LGCP

6.5.1 The model and your data

If the event locations are recorded over a period of time, the spatial pattern may change over time. A spatiotemporal LGCP allows this by having a separate random field at each of a set of time knots \(t_1, t_2, ..., t_K\):

\[\lambda_k(s) = \text{exp}(\beta_0 + G_k(s)),\]

where the fields are linked over time by a first order autoregressive (AR(1)) process,

\[G_k(s) = \rho\, G_{k-1}(s) + \epsilon_k(s),\]

and the \(\epsilon_k(s)\) are independent random fields. All of the fields share the same range and standard deviation. The new parameter is \(\rho\) (\(-1 < \rho < 1\)), the correlation between the fields at consecutive time knots: \(\rho\) close to 1 means the spatial pattern of events is very persistent over time, \(\rho\) close to 0 means the pattern at one knot is almost unrelated to the pattern at the previous one. get_coefs() reports rho along with the parameters of the spatial model, and atanh_rho is the transformed parameter used when giving starting values.

To fit the model you need, in addition to a spatial model (Section 6.1.1):

• a column t in locs giving the time of each event; and
• a temporal mesh tmesh, built using fmesher::fm_mesh_1d() with the knots as its argument, e.g., fmesher::fm_mesh_1d(0:5) for knots at \(t = 0, 1, ..., 5\).

Each event is assigned to its nearest time knot, so knot \(t_k\) covers the time period from halfway between \(t_{k-1}\) and \(t_k\) to halfway between \(t_k\) and \(t_{k+1}\), and the first and last knots cover only half a period. Things to be aware of:

• The times and knots must be on the same scale. For example, if t is in years since the start of the study, the knots should be too. Check tmesh$loc to see the knots that were created.
• All times must lie between the first and last knots (fit_lgcp() gives an error otherwise). For example, with yearly data starting at \(t = 0\) and finishing within year 5, use continuous (decimal) times and knots 0:5, so the last knot covers the final half year. Knots whose period contains no events are allowed, but they carry little information.
• The cost grows with the number of knots, since there is a field at every knot: roughly the number of mesh nodes times the number of knots. A coarser spatial mesh than for a purely spatial model is often needed.

6.5.2 A simulated example

Data were simulated from a spatiotemporal LGCP within a disc of radius 2 over the time period \([0, 5)\), with \(\beta_0 = 2\), range \(r = 1\), standard deviation \(\sigma = 1\), and \(\rho = 0.7\), with one field per unit of time. The simulation has three steps:

1. Independent fields. At each of the \(K = 6\) time knots (\(t = 0, 1, ..., 5\)) an independent random field \(\epsilon_k(s)\) is simulated at the mesh nodes using sim_lgcp() (only the simulated field, $x, is used). sim_lgcp() takes the field parameters on the internal scale, so the range and standard deviation are first converted to \(\kappa = \sqrt{8}/r\) and \(\tau = 1/(\sigma\sqrt{4\pi}\kappa)\).
2. AR(1) fields. The fields are linked over time, \(G_k(s) = \rho\, G_{k-1}(s) + \sqrt{1 - \rho^2}\, \epsilon_k(s)\). The factor \(\sqrt{1 - \rho^2}\) keeps the standard deviation of every \(G_k(s)\) equal to \(\sigma\).
3. Events. The events in the time period covered by each knot (Section 6.5) are simulated by thinning: candidate events are simulated uniformly over the disc at the maximum of the intensity \(\lambda_{\max}\), and each candidate at location \(s\) is kept with probability \(\lambda_k(s)/\lambda_{\max}\), with the field at \(s\) interpolated from the mesh nodes using fmesher::fm_basis(). The kept events are given times drawn uniformly within the period.

disc_sf <- sf::st_sf(geometry = sf::st_geometry(sf::st_as_sf(spatstat.geom::disc(radius = 2))))
smesh <- fmesher::fm_mesh_2d(boundary = fmesher::fm_as_segm(disc_sf), max.edge = 0.2, cutoff = 0.05)
tmesh <- fmesher::fm_mesh_1d(0:5)
k <- tmesh$n
truth_st <- c(beta = 2, range = 1, sigma = 1, rho = 0.7)
kap <- unname(sqrt(8) / truth_st["range"]); tau <- unname(1 / (truth_st["sigma"] * sqrt(4 * pi) * kap))
set.seed(3)
## 1. independent fields, one column per time knot
e <- sapply(1:k, function(i)
    as.vector(sim_lgcp(c(beta = 0, log_tau = log(tau), log_kappa = log(kap)), sf = disc_sf, smesh = smesh)$x))
## 2. AR(1) fields
x <- e
for (i in 2:k) x[, i] <- truth_st["rho"] * x[, i - 1] + sqrt(1 - truth_st["rho"]^2) * e[, i]
## 3. events, by thinning, for the time period covered by each knot
locs <- do.call(rbind, lapply(1:k, function(i) {
    lo <- max(0, tmesh$loc[i] - 0.5); hi <- min(5, tmesh$loc[i] + 0.5)
    lmax <- exp(truth_st["beta"] + max(x[, i]))
    n0 <- rpois(1, lmax * 16 * (hi - lo)) ## candidates over the 4 by 4 square around the disc
    cand <- data.frame(x = runif(n0, -2, 2), y = runif(n0, -2, 2))
    cand <- cand[sqrt(cand$x^2 + cand$y^2) < 2, ]
    A <- fmesher::fm_basis(smesh, loc = as.matrix(cand))
    keep <- runif(nrow(cand)) < exp(truth_st["beta"] + as.vector(A %*% x[, i])) / lmax
    data.frame(cand[keep, ], t = runif(sum(keep), lo, hi))
}))
locs <- locs[locs$t < 5, ]
locs <- locs[order(locs$t), ]
rownames(locs) <- NULL

The data are a data frame of 891 events with columns x, y, and t:

head(locs)

           x           y           t
1  1.9298853  0.31397382 0.008389111
2 -0.1897221  1.06251820 0.009234515
3  1.7794098  0.41180309 0.010330225
4  1.6102736  0.02557874 0.011357767
5  0.2621431 -0.15527212 0.018003215
6  1.7512718  0.49230152 0.034940886

The spatial mesh follows the boundary of the study region and the temporal mesh has knots at \(0, 1, ..., 5\). These are the meshes used to simulate the data, repeated here so that the code needed to fit the model is all in one place.

disc_sf <- sf::st_sf(geometry = sf::st_geometry(sf::st_as_sf(spatstat.geom::disc(radius = 2))))
smesh <- fmesher::fm_mesh_2d(boundary = fmesher::fm_as_segm(disc_sf), max.edge = 0.2, cutoff = 0.05)
tmesh <- fmesher::fm_mesh_1d(0:5)
tmesh$loc

[1] 0 1 2 3 4 5

fit_st <- fit_lgcp(locs = locs, sf = disc_sf, smesh = smesh, tmesh = tmesh,
                   parameters = c(beta = 1, log_tau = log(1), log_kappa = log(1), atanh_rho = 0))
stelfi_st <- get_coefs(fit_st)
stelfi_st

            Estimate Std. Error
rho        0.5998841 0.07272397
beta       1.8393273 0.32568893
log_tau   -2.1073093 0.10586796
log_kappa  0.7947107 0.13262724
range      1.2776342 0.16944910
stdev      1.0482127 0.09213091

The same model can be fitted using inlabru, by giving the field a group (the time) with an AR(1) control.group, and including time in the domain of the point process.

matern <- INLA::inla.spde2.pcmatern(smesh, prior.range = c(0.1, 0.01), prior.sigma = c(5, 0.01))
pts <- sf::st_as_sf(locs, coords = c("x", "y"))
cmp <- ~ field(geometry, model = matern, group = t, group_mapper = bm_fmesher(tmesh),
               control.group = list(model = "ar1")) + Intercept(1)
fit_bru_st <- lgcp(cmp, data = pts, samplers = disc_sf,
                   domain = list(geometry = smesh, t = tmesh),
                   formula = geometry + t ~ field + Intercept,
                   options = bru_options(control.inla = list(int.strategy = "eb")))
bru_st <- rbind(fit_bru_st$summary.fixed[, 1:2], fit_bru_st$summary.hyperpar[, 1:2])
bru_st

                        mean         sd
Intercept          1.9663568 0.43998045
Range for field    1.7118848 0.24379471
Stdev for field    1.1097764 0.11783784
GroupRho for field 0.6048416 0.07667948

params <- data.frame(Truth = truth_st,
                     stelfi = apply(stelfi_st[c("beta", "range", "stdev", "rho"), ], 1, fmt),
                     inlabru = apply(bru_st, 1, fmt))
knitr::kable(t(params), col.names = c("$\\beta_0$", "$r$", "$\\sigma$", "$\\rho$"), escape = FALSE)

	\(\beta_0\)	\(r\)	\(\sigma\)	\(\rho\)
Truth	2.0	1.0	1.0	0.7
stelfi	1.839 (0.326)	1.278 (0.169)	1.048 (0.092)	0.6 (0.073)
inlabru	1.966 (0.44)	1.712 (0.244)	1.11 (0.118)	0.605 (0.077)

Both packages give similar estimates and recover the true values reasonably well. The estimated field at any time knot can be extracted using get_fields() (which, with tmesh, returns a matrix with one column per knot) and the intensity plotted using show_lambda() with the timestamp argument (the index of the knot).

plots <- lapply(c(1, 3, 6), function(k) {
    show_lambda(fit_st, smesh = smesh, sf = disc_sf, tmesh = tmesh, timestamp = k, clip = TRUE) +
        theme_void() + ggtitle(paste("t =", tmesh$loc[k]))
})
patchwork::wrap_plots(plots, nrow = 1)

6.5.3 An applied example: terrorism events in Iraq by year

The iraq_terrorism data (Section 6.3.2) also record the year, month and day of each event, so the intensity can be allowed to change over time. Here time is measured in years since the start of 2013 (as a decimal, using the month and day) and there is one knot per year.

locs_st <- as.data.frame(sf::st_coordinates(sf::st_transform(iraq_terrorism, crs = 32638)) / 1000)
names(locs_st) <- c("x", "y")
locs_st$t <- (iraq_terrorism$iyear - 2013) + (iraq_terrorism$imonth - 1) / 12 +
    (iraq_terrorism$iday - 1) / 365
locs_st <- locs_st[order(locs_st$t), ]
range(locs_st$t)

[1] 0.2965753 4.9988584

The spatial mesh follows the boundary of Iraq and is coarser than for the spatial model because there are six fields to estimate. The temporal knots are at \(t = 0, 1, ..., 5\) (the start of each year), so the final knot covers the second half of 2017.

smesh_st <- fmesher::fm_mesh_2d(boundary = fmesher::fm_as_segm(iraq_sf_km), max.edge = 50, cutoff = 20)
smesh_st$n

[1] 496

tmesh_st <- fmesher::fm_mesh_1d(0:5)

Starting values matter for this model. For the intercept, the log of the average number of events per km\(^2\) per year is a sensible start. For the range a starting value of around 300 km is used (\(\kappa = \sqrt{8}/300\)).

n_years <- 5
start <- c(beta = log(nrow(locs_st) / (area_km2 * n_years)),
           log_tau = log(0.1), log_kappa = log(sqrt(8) / 300), atanh_rho = 0.2)
fit_iraq_st <- fit_lgcp(locs = locs_st, sf = iraq_sf_km, smesh = smesh_st, tmesh = tmesh_st,
                        parameters = start)
coefs_st <- get_coefs(fit_iraq_st)
coefs_st

             Estimate Std. Error
rho         0.9220853 0.01286857
beta      -11.5955606 0.56488818
log_tau     1.2540996 0.08731190
log_kappa  -3.5544586 0.09573922
range      98.9069365 9.46927254
stdev       2.8146752 0.19977800

The estimated correlation between the fields of consecutive years is \(\hat{\rho} =\) 0.92, so the spatial pattern of events is very persistent from year to year. The estimated intensity in each year is plotted below.

plots <- lapply(1:6, function(k) {
    show_lambda(fit_iraq_st, smesh = smesh_st, sf = iraq_sf_km, tmesh = tmesh_st,
                timestamp = k, clip = TRUE) +
        theme_void() + ggtitle(paste("t =", tmesh_st$loc[k]))
})
patchwork::wrap_plots(plots, nrow = 2)

If the model does not converge, things to try are a coarser spatial mesh, fewer knots, different starting values (in particular beta, log_tau and atanh_rho), or fitting a spatial LGCP (Section 6.3.2) to each year first to get an idea of sensible values for \(\beta_0\), \(\tau\), and \(\kappa\).

6.5.4 Goodness-of-fit

The simulation envelopes of Section 6.4 can be used for a single time period. Below a slightly different check is used: rather than simulating new random fields, events are simulated from the estimated intensity for one knot (an inhomogeneous Poisson process), and the \(G\) function of these is compared with that of the observed events in the corresponding period. This asks whether, given the estimated intensity surface, the events are more clustered than expected. Here the knot at \(t = 2\) (mid 2014 to mid 2015) is checked.

knot <- 3 ## the index of the knot at t = 2
in_period <- abs(locs_st$t - tmesh_st$loc[knot]) < 0.5
obs <- locs_st[in_period, ]
iraq_window <- as.owin(iraq_sf_km)
set.seed(1)
obs_ppp_st <- rjitter(ppp(obs$x, obs$y, window = iraq_window), retry = TRUE)
## estimated log intensity at the mesh nodes for this knot
eta <- coefs_st["beta", 1] + get_fields(fit_iraq_st, smesh_st, tmesh_st)[, knot]
## simulate from the estimated intensity (per km^2 per year) by thinning
sim_from_fitted <- function() {
    bb <- sf::st_bbox(iraq_sf_km)
    lmax <- exp(max(eta))
    n0 <- rpois(1, lmax * (bb["xmax"] - bb["xmin"]) * (bb["ymax"] - bb["ymin"]))
    cand <- ppp(runif(n0, bb["xmin"], bb["xmax"]), runif(n0, bb["ymin"], bb["ymax"]),
                window = owin(bb[c("xmin", "xmax")], bb[c("ymin", "ymax")]))
    cand <- cand[iraq_window]
    lambda <- exp(as.vector(fmesher::fm_basis(smesh_st, loc = cbind(cand$x, cand$y)) %*% eta))
    cand[runif(cand$n) < lambda / lmax]
}
sims <- lapply(1:39, function(i) sim_from_fitted())
g_env_st <- envelope(obs_ppp_st, fun = Gest, simulate = sims, nsim = 39, verbose = FALSE)

plot(g_env_st, main = "Iraq, t = 2: G function", xlab = "Distance (km)", legend = FALSE)

The observed \(G\) function lies well above the envelope at short distances: many events are much closer to another event than would be expected given the estimated intensity. This is not surprising: the mesh (and the knots) can only describe variation in intensity at scales of tens of kilometres and a year, whereas attacks are often concentrated in particular towns and occur in bursts. A finer mesh would help, but this kind of small scale clustering in space and time is also exactly what a self-exciting model is designed to capture, see Chapter 8.

Bachl, Fabian E., Finn Lindgren, David L. Borchers, and Janine B. Illian. 2019. “inlabru: An R Package for Bayesian Spatial Modelling from Ecological Survey Data.” Methods in Ecology and Evolution 10: 760–66. https://doi.org/10.1111/2041-210X.13168.

Baddeley, Adrian, Ege Rubak, and Rolf Turner. 2015. Spatial Point Patterns: Methodology and Applications with R. London: Chapman; Hall/CRC Press.

Jones-Todd, Charlotte M., and Alec B. M. van Helsdingen. 2024. “stelfi: An R Package for Fitting Hawkes and Log-Gaussian Cox Point Process Models.” Ecology and Evolution 14: e11005. https://doi.org/10.1002/ece3.11005.

Lindgren, Finn, Håvard Rue, and Johan Lindström. 2011. “An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach (with Discussion).” Journal of the Royal Statistical Society B 73 (4): 423–98.