Comparing two spatial patterns with sppt()

The spatial point pattern test asks whether two sets of events, for example residential burglaries in two different years, are distributed the same way across a set of areal units. sppt() answers that unit by unit, and then summarises the whole map.

Getting started

The package includes residential (BNER) and commercial (BNEC) burglary counts for Vancouver from 2003 to 2025, in 127 census tracts (vancouver_ct) and 1,016 dissemination areas (vancouver_da).

To compare residential burglary in 2003 (the base) with 2025 (the test), name the identifier column and the two count columns. alpha = 0.05 is the default and is shown here for clarity.

res <- sppt(
  data      = vancouver_ct,
  group_col = "id",
  count_col = c("BNER2003", "BNER2025"),
  B         = 999,
  alpha     = 0.05,
  seed      = 42,
  create_maps = FALSE
)
#> 
#> ==================================================================
#> Spatial Point Pattern Test
#> ==================================================================
#>   Base variable : BNER2003  (N = 6875)
#>   Test variable : BNER2025  (N = 885)
#>   Spatial units : 127
#>   lambda_min    : 6.969
#>   Bootstrap B   : 999
#>   P-value alpha : 0.050
#>   Decision tier : Medium-high scale (6 <= lambda_min < 17): PV alpha = 0.05 recommended.
#> ------------------------------------------------------------------
#>   Robust S-index  : 0.6378  (46 of 127 units dissimilar at alpha = 0.050)
#>     BNER2025 > BNER2003 : 13 units
#>     BNER2003 > BNER2025 : 33 units
#>   Sameness test   : S = 0.6378 vs 0.9500 expected, p < 0.0001 (exact binomial)
#>     Threshold     : S >= 0.9182 at tau = 0.05
#>     Verdict       : more different than two realisations of one process
#>   Arrangement     : Moran's I = +0.3451 on signed evidence, p = 0.0010 (999 perms)
#>     Local clusters: 10 units where BNER2025 > BNER2003, 9 units where BNER2003 > BNER2025
#> 
#>   Two-way verdict
#>                        | arrangement random    | arrangement clustered
#>     -----------------------------------------------------------------
#>     count at chance    | same process          | few differences, but concentrated
#>     count above chance | many, diffuse          | differences real and localised <<<
#>     This result        : count above chance, arrangement clustered
#> ==================================================================

The console report gives:

lambda_min, the mean number of events per unit in the rarer pattern, which governs power, and the decision tier that follows from it;
the robust S index, the share of units (not empty in both years) whose share of events did not change significantly;
the sameness test, which asks whether that S index is lower than two realisations of one process would give;
the arrangement test, Moran’s I on the signed evidence surface, which asks whether the differences cluster in space.

summary() repeats the key results at any time.

summary(res)
#> 
#> ==================================================================
#> SPPT Summary
#> ==================================================================
#>   Base             : BNER2003
#>   Test             : BNER2025
#>   Spatial units    : 127  (127 eligible, 0 zero-zero excluded)
#>   lambda_min       : 6.969
#>   Bootstrap B      : 999
#>   Robust S-index   : 0.6378  (alpha = 0.050; 46 dissimilar)
#>   Sameness test    : p < 0.0001  (expected S = 0.9500, SE = 0.0193, exact binomial)
#>   Sameness cut-off : S >= 0.9182 at tau = 0.05  [different]
#>   Arrangement      : Moran's I = +0.3451 on signed evidence, p = 0.0010
#>   Local clusters   : 10 (BNER2025 > BNER2003), 9 (BNER2003 > BNER2025)
#>   Decision tier    : Medium-high scale (6 <= lambda_min < 17): PV alpha = 0.05 recommended.
#> ==================================================================

Each unit gets a p-value and a direction: 1 where the test pattern has the larger share, -1 where the base has, and 0 where the difference is not significant.

table(res$sig_pvalue_dir, useNA = "ifany")
#> 
#> -1  0  1 
#> 33 81 13

The p-value is the bounded Monte Carlo form \(p = (1 + r)/(B + 1)\), so it is never exactly zero; with B = 999 the smallest possible value is 0.001. B = 4999 (the default) is recommended for final analyses.

Maps

With create_maps = TRUE (the default), sppt() prints a difference map and a significance map and stores both on the result. sppt_map() redraws them, with the legend, north arrow and scale bar wherever suits your study area.

sppt_map(res, "pvalue")

sppt_map(res, "pvalue", legend_position = "right",
         north_arrow = "tl", scale_bar = "br", title = FALSE)

The maps are ordinary ggplot objects, so anything else can be changed by adding to them, for example sppt_map(res, "difference") + ggplot2::scale_fill_viridis_c().

Saving results

Nothing is written to disk unless you ask. To save the maps and the result, name a directory:

res <- sppt(vancouver_ct, group_col = "id",
            count_col = c("BNER2003", "BNER2025"),
            export_maps = TRUE, export_results = TRUE,
            export_format = "gpkg", export_dir = "sppt_output")

export_format = "rds" keeps the summary and the maps with the data; readRDS() restores everything, and summary() and sppt_map() work on the reloaded object. The result is also an ordinary sf object, so sf::st_write() works on it directly.

All options

The call below lists every argument with its default value.

res <- sppt(
  data               = vancouver_ct,
  group_col          = "id",
  count_col          = c("BNER2003", "BNER2025"), # base, then test
  B                  = 4999,    # bootstrap replicates; odd avoids ties
  seed               = NULL,    # set for reproducible results
  alpha              = 0.05,    # unit-level significance level
  calculate_pvalues  = TRUE,    # bootstrap p-values
  sameness_tau       = 0.05,    # significance level of the sameness test
  arrangement        = TRUE,    # Moran's I on the signed evidence surface
  nsim_moran         = 999,     # permutations for global Moran's I
  calculate_ci       = FALSE,   # also compute marginal CI overlap (2009)
  conf_level         = 0.683,   # CI level; matches alpha = 0.05
  fix_base           = FALSE,   # TRUE = do not bootstrap the base pattern
  use_percentages    = TRUE,    # percentages (0-100) rather than proportions
  verbose            = TRUE,    # print the console report
  create_maps        = TRUE,    # build, print and store the maps
  map_crs            = NULL,    # CRS for the maps; NULL keeps that of data
  export_maps        = FALSE,   # save maps as PNG in export_dir
  export_dir         = NULL,    # required for any export
  map_dpi            = 300,     # resolution of exported maps
  export_results     = FALSE,   # save the result
  export_format      = "gpkg",  # gpkg, shp, csv, txt or rds
  export_results_dir = NULL     # defaults to export_dir
)

Confidence interval overlap

The original test (Andresen 2009) compared bootstrap confidence intervals. calculate_ci = TRUE adds that method alongside the p-values, using the same bootstrap replicates. The default conf_level = 0.683 (about one standard deviation per tail) matches the performance of alpha = 0.05 p-values; 0.95 reproduces the original, more conservative analysis.

Fixed base

fix_base = TRUE treats the base pattern as known rather than sampled. That suits a comparison against a theoretical or population distribution, but the test is then asymmetric, and the sameness test is not reported.