1.2 First-order Analysis: Spatial Intensity

Part: 1: Point Pattern Analysis (Chapter 1.2 of 5)

Dataset: Buffalo crime incidents, 2018 (data/crime_data_2018.csv, data/buffalo_city.shp)

Focus: Estimating and mapping how event density varies across a study region: quadrat counts, kernel density estimation, bandwidth selection, and relative-risk surfaces

Library: spatialstats.pointpattern (PySpatialStats)

---

Table of Contents

1. Overview

2. What is Spatial Intensity?

3. Estimating Intensity: Two Approaches

   • 3.1 The quadrat-count test of CSR

   • 3.2 Kernel density estimation

   • 3.3 Bandwidth selection

   • 3.4 Edge correction

4. Python Setup

5. Rebuild the Point Pattern

6. Quadrat Analysis of Crime Intensity

7. Kernel Density Estimation of Crime Intensity

8. Bandwidth Selection in Practice

9. Relative Risk: Comparing Crime Types

10. Summary and Conclusions

11. Resources

---

1. Overview

Chapter 1.1 established that Buffalo's 2018 crime events are not uniformly spread across the city — a first, informal impression from comparing the real data to a simulated CSR pattern. This chapter makes that impression precise: it estimates the intensity function $\lambda(\mathbf{s})$, the local density of events per unit area, two different ways, and formally tests whether a constant intensity (CSR) is consistent with the data at all.

Step	What we do
Theory	Define intensity; contrast binned (quadrat) and smooth (kernel) estimators
Quadrat test	Partition the city into cells; test observed counts against the Poisson counts CSR predicts
Kernel density	Estimate a smooth intensity surface with edge correction for the irregular city boundary
Bandwidth	Compare a rule-of-thumb bandwidth (Scott) against one chosen by cross-validation
Relative risk	Map where one crime type is over- or under-represented relative to crime overall

By the end you will have quantified how far Buffalo crime departs from CSR in intensity alone (before any question of event-to-event interaction, which is Chapter 1.3's job), and produced maps a crime analyst could actually use to prioritise areas.

---

2. What is Spatial Intensity?

Recall from Chapter 1.1 the intensity function

$$ \lambda(\mathbf{s}) = \lim_{|d\mathbf{s}| \to 0} \frac{\mathbb{E}\,N(d\mathbf{s})}{|d\mathbf{s}|}, $$

the expected number of events per unit area at location $\mathbf{s}$. Under CSR, $\lambda(\mathbf{s}) = \lambda$ is a single constant, estimated by the crude global rate $\hat\lambda = n / |W|$ (Chapter 1.1's pp.intensity). This chapter asks a sharper question: does $\lambda$ vary over $\mathbf{s}$, and if so, how?

Two estimation strategies, with different strengths:

Approach	Idea	Good for
Quadrat counts	Partition $W$ into cells; compare observed counts to what a constant $\lambda$ predicts	A formal test of CSR (a single p-value); simple, transparent
Kernel density	Smooth each event into a small "bump" of probability mass and sum the bumps	A continuous map of $\hat\lambda(\mathbf{s})$; better for visual and comparative work

They answer the same underlying question at different resolutions: quadrat counts summarise intensity variation into one test statistic, kernel density spreads it into a full surface. Both are used below, and both point to the same conclusion.

---

3. Estimating Intensity: Two Approaches

3.1 The quadrat-count test of CSR

Divide the bounding box of $W$ into an $n_x \times n_y$ grid of cells. Under CSR, the count in cell $c$ is

$$ O_c \sim \text{Poisson}\big(\lambda \cdot |c \cap W|\big), $$

so the expected count is $E_c = n \cdot |c \cap W| / |W|$ — proportional to the area of the cell that actually lies inside the window, which matters for a city boundary as irregular as Buffalo's. The Pearson chi-square statistic

$$ X^2 = \sum_{c} \frac{(O_c - E_c)^2}{E_c} \;\overset{\text{approx.}}{\sim}\; \chi^2_{k-1} \quad \text{under CSR} $$

(with $k$ the number of cells with $E_c > 0$) tests departure from a constant intensity. A companion statistic, the index of dispersion $\text{IoD} = s^2_O / \bar{O}$ (variance over mean of the cell counts), is 1 under CSR, greater than 1 for clustered/overdispersed patterns, and less than 1 for regular ones — informal, but intuitive, and does not depend on the $\chi^2$ approximation holding.

The chi-square approximation degrades when expected counts are small (a common rule of thumb: at least 80% of cells should have $E_c \geq 5$). quadrat_test also offers a Monte Carlo p-value (permutations=), drawing count vectors from the multinomial distribution implied by CSR — valid regardless of cell size.

3.2 Kernel density estimation

The kernel estimator replaces each event with a small bump of mass and sums them:

$$ \hat\lambda(\mathbf{s}) = \sum_{i=1}^n \frac{1}{h^2}\, k\!\left(\frac{\mathbf{s} - \mathbf{s}_i}{h}\right), $$

where $k$ is a kernel function (Gaussian, or the compact Epanechnikov kernel) and $h$ the bandwidth. This package implements $k$ as a product (separable) kernel, $k(\mathbf{u}) = k_1(u_x)\,k_1(u_y)$, which makes the computation an efficient matrix product rather than an $O(n \times n_{\text{grid}}^2)$ loop.

3.3 Bandwidth selection

$h$ controls a bias–variance trade-off: small $h$ tracks individual events closely (low bias, high variance — a speckled, noisy surface); large $h$ averages over a wide area (low variance, high bias — a smooth surface that can wash out genuine local structure). Two choices are compared below:

• Scott's rule of thumb: $h = \hat\sigma \cdot n^{-1/6}$, with $\hat\sigma$ the mean coordinate standard deviation — fast, data-driven only through the pattern's overall spread, and typically oversmooths clustered spatial data (it was derived for roughly Gaussian-shaped densities, not point-process clusters).
• Likelihood cross-validation: choose $h$ to maximise the leave-one-out log-likelihood $\sum_i \log \hat\lambda_{-i}(\mathbf{s}_i)$ — adapts to the data, but is well known to undersmooth clustered patterns (it keeps shrinking $h$ to fit the tightest local clusters more closely).

Section 8 shows both are visibly biased in opposite directions here, which is itself the practical lesson: neither is a substitute for looking at the map.

3.4 Edge correction

Near the window boundary, a naive kernel estimate is biased low: part of the kernel's mass would fall outside $W$ (where there could not possibly have been an event), so the sum under-counts. kde(..., edge_correction=True) divides by the fraction of kernel mass that actually lands inside $W$ (the "uniform" or Jones correction), which is exact for any window shape, not just rectangles — important for Buffalo's non-convex boundary (Chapter 1.1, Section 10).

---

4. Python Setup

import sys, warnings
from pathlib import Path

import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from IPython.display import Markdown, display

warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=DeprecationWarning)


def find_repo_root(start: Path = Path.cwd()) -> Path:
    for p in [start, *start.parents]:
        if (p / "pyproject.toml").exists() and (p / "data").is_dir():
            return p
    raise FileNotFoundError("Run this notebook from inside the PySpatialStats repository.")


ROOT = find_repo_root()
DATA = ROOT / "data"
try:
    import spatialstats
except ImportError:
    sys.path.insert(0, str(ROOT))
    import spatialstats

from spatialstats.pointpattern import Window, PointPattern, quadrat_test, kde, bandwidth_scott, bandwidth_cv, relative_risk

plt.rcParams.update({
    "figure.dpi": 100, "savefig.dpi": 100, "font.size": 10, "axes.titlesize": 11,
    "axes.titleweight": "bold", "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25, "legend.frameon": False,
})
C = {"blue": "#2b6cb0", "red": "#c0392b", "green": "#2f855a", "orange": "#dd6b20",
     "grey": "#718096", "purple": "#6b46c1", "teal": "#2c7a7b"}

print(f"spatialstats : version {spatialstats.__version__}")
print(f"data folder  : {DATA.relative_to(ROOT)}/")

spatialstats : version 0.2.0
data folder  : data/

---

5. Rebuild the Point Pattern

Chapter 1.1 diagnosed a coordinate reference system error in the source shapefiles (declared UTM zone 11N; Buffalo is in zone 17N) and fixed it by reprojecting through longitude/latitude to EPSG:26917. We repeat that short, corrected build here so this notebook is self-contained; see Chapter 1.1 for the full diagnosis.

from pyproj import Transformer

tbl = pd.read_csv(DATA / "crime_data_2018.csv")
city = gpd.read_file(DATA / "buffalo_city.shp")

CRS_CORRECT = 26917  # NAD83 / UTM zone 17N -- the correct zone for Buffalo (Chapter 1.1)
city_fixed = city.to_crs(4326).to_crs(CRS_CORRECT)
transformer = Transformer.from_crs(4326, CRS_CORRECT, always_xy=True)
x_m, y_m = transformer.transform(tbl["long"].to_numpy(), tbl["lat"].to_numpy())
coords_all = np.column_stack([x_m, y_m])

win_city = Window.from_gdf(city_fixed)
pp_all = PointPattern(coords_all, window=win_city, marks=tbl["incident_type"].to_numpy(), clip=True)

print(pp_all)
print(f"global (crude) intensity : {pp_all.intensity * 1e6:.2f} events / km²")

PointPattern(n=13948, intensity=0.0001318, window=Window(polygon, area=1.0583e+08, bounds=(670377.394, 4743812.67, 680296.782, 4759265.808)))
global (crude) intensity : 131.80 events / km²

---

6. Quadrat Analysis of Crime Intensity

We start with an $8 \times 8$ grid over the city and a Monte Carlo p-value (199 multinomial simulations) alongside the analytical $\chi^2$ approximation.

qd = quadrat_test(pp_all, nx=8, ny=8, permutations=199, seed=1)

print(qd.summary())
p_analytic = "< 0.0001" if qd["p_value"] < 1e-4 else f"= {qd['p_value']:.4f}"
display(Markdown(
    f"**Index of dispersion = {qd['index_of_dispersion']:.1f}** — the variance of cell counts is "
    f"{qd['index_of_dispersion']:.0f} times the mean. Under CSR this ratio would be 1. Both the analytical "
    f"$\\chi^2$ test (p {p_analytic}) and the Monte Carlo test (p = {qd['p_mc']:.3f}, the smallest value 199 "
    "permutations can report) reject CSR decisively."
))

Quadrat test of CSR
===================
  chi2                 5572.08
  df                   57
  p_value              0
  p_clustered          0
  p_regular            1
  index_of_dispersion  166.867
  n_cells              58
  min_expected         6.91743
  p_mc                 0.005

Index of dispersion = 166.9 — the variance of cell counts is 167 times the mean. Under CSR this ratio would be 1. Both the analytical $\chi^2$ test (p < 0.0001) and the Monte Carlo test (p = 0.005, the smallest value 199 permutations can report) reject CSR decisively.

fig, ax = plt.subplots(figsize=(7, 7))
qd.plot(ax=ax, cmap="YlOrRd")
city_fixed.boundary.plot(ax=ax, color="k", lw=0.4)
ax.set_title(f"Observed counts, {8}×{8} quadrats\nIoD={qd['index_of_dispersion']:.0f}, "
             f"$\\chi^2$={qd['chi2']:.0f}, df={qd['df']}")
plt.tight_layout(); plt.show()

Sensitivity to cell size

The chi-square approximation needs expected counts that are not too small. We check a range of grid resolutions, reporting the minimum expected count per cell (min_expected) alongside the index of dispersion, so the trade-off is visible rather than assumed.

rows = []
for nx in [4, 6, 8, 10, 12, 16]:
    q = quadrat_test(pp_all, nx, nx)
    rows.append({"grid": f"{nx}×{nx}", "n_cells": q["n_cells"], "min_expected": q["min_expected"],
                 "chi2": q["chi2"], "index_of_dispersion": q["index_of_dispersion"]})
sens = pd.DataFrame(rows).set_index("grid")
sens.round(2)

	n_cells	min_expected	chi2	index_of_dispersion
grid
4×4	15	210.23	2289.57	273.40
6×6	33	38.06	4087.72	210.25
8×8	58	6.92	5572.08	166.87
10×10	86	0.31	6388.13	113.51
12×12	123	1.28	6794.74	85.74
16×16	217	0.03	9136.28	64.00

The index of dispersion falls as cells get smaller — expected, since a very fine grid starts to bump into genuinely unpopulated cells (parks, industrial land, the river) that shrink the variance contributed by "real" crime hot spots relative to zero-count cells everywhere. min_expected also falls well below the conventional rule-of-thumb of 5 once the grid is finer than about $8\times 8$, so the $\chi^2$ approximation (not the Monte Carlo p-value, which stays valid) becomes questionable there. Every resolution still shows overdispersion far beyond anything CSR would produce — the qualitative conclusion is robust even where the exact chi-square p-value is not.

---

7. Kernel Density Estimation of Crime Intensity

kde returns a KDEResult holding the surface on a regular grid, scaled so it integrates to (approximately) the event count $n$ — a useful sanity check on any kernel implementation.

h_scott = bandwidth_scott(pp_all)
kd = kde(pp_all, bandwidth=h_scott, grid_size=200)

print(f"Scott's-rule bandwidth : {h_scott:.0f} m")
print(f"grid                   : {kd.density.shape[1]} x {kd.density.shape[0]} cells")
print(f"integral of density    : {kd.integral:.0f}  (n = {pp_all.n:,}; "
      f"{100 * (kd.integral / pp_all.n - 1):+.1f}% — discretisation/edge-correction error)")

fig, ax = plt.subplots(figsize=(7, 7))
kd.plot(ax=ax, points=pp_all, cmap="inferno")
city_fixed.boundary.plot(ax=ax, color="w", lw=0.5, alpha=0.6)
ax.set_title(f"Kernel intensity of crime, 2018 (h = {h_scott:.0f} m)")
plt.tight_layout(); plt.show()

Scott's-rule bandwidth : 580 m
grid                   : 129 x 200 cells
integral of density    : 14160  (n = 13,948; +1.5% — discretisation/edge-correction error)

---

8. Bandwidth Selection in Practice

bandwidth_cv maximises leave-one-out likelihood over a geometric grid of candidate bandwidths (subsampling large patterns for speed — the default cap is 5,000 events, comfortably above what is needed to locate the optimum). We compare it against Scott's rule.

h_cv, candidates, scores = bandwidth_cv(pp_all, seed=1)

print(f"Scott's rule            : h = {h_scott:.0f} m")
print(f"Likelihood CV           : h = {h_cv:.0f} m")
print(f"ratio (Scott / CV)      : {h_scott / h_cv:.2f}×")

fig, ax = plt.subplots(figsize=(6.5, 4.2))
ax.plot(candidates, scores, "o-", color=C["blue"], ms=4)
ax.axvline(h_cv, color=C["red"], ls="--", label=f"CV optimum: {h_cv:.0f} m")
ax.axvline(h_scott, color=C["grey"], ls=":", label=f"Scott's rule: {h_scott:.0f} m")
ax.set_xlabel("bandwidth h (m)"); ax.set_ylabel("leave-one-out log-likelihood")
ax.set_title("Cross-validation score by bandwidth")
ax.legend(fontsize=9)
plt.tight_layout(); plt.show()

Scott's rule            : h = 580 m
Likelihood CV           : h = 136 m
ratio (Scott / CV)      : 4.27×

h_mid = float(np.sqrt(h_cv * h_scott))   # a visual midpoint on a log scale, for contrast only
bandwidths = {"Cross-validation (undersmooths)": h_cv, "Geometric midpoint": h_mid,
              "Scott's rule (oversmooths)": h_scott}

fig, axes = plt.subplots(1, 3, figsize=(16, 5.6))
for ax, (label, h) in zip(axes, bandwidths.items()):
    k = kde(pp_all, bandwidth=h, grid_size=180)
    im = ax.imshow(k.density * 1e6, origin="lower", extent=k.extent, cmap="inferno")
    plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label="events / km²")
    ax.set_aspect("equal")
    ax.set_title(f"{label}\nh = {h:.0f} m")
plt.tight_layout(); plt.show()

The contrast is exactly what the theory in Section 3.3 predicts: the cross-validation bandwidth resolves many small, distinct hot spots — plausibly true local structure, but also plausibly some noise, since likelihood CV is known to undersmooth clustered spatial patterns; Scott's rule washes almost everything into two or three broad regions, useful for a citywide overview but too coarse to act on locally. Neither is automatically "right": the geometric midpoint is shown only as a visual compromise, not a principled recommendation. In practice, report the bandwidth used, and show at least two to let the reader judge — exactly as this section does.

---

9. Relative Risk: Comparing Crime Types

A map of where crime happens is dominated by where people and property are. A sharper operational question is: where is one crime type over-represented relative to crime overall? relative_risk(cases, controls) answers this with the ratio of two kernel density estimates — cases (e.g. Assault) against controls (every other incident type) — so that areas with generally high crime no longer automatically dominate the map.

assault = pp_all.by_mark("Assault")
other = pp_all.subset(pp_all.marks != "Assault")
print(f"Assault (cases)      : n = {assault.n:,}")
print(f"other crime (control): n = {other.n:,}")

rr = relative_risk(assault, other, bandwidth="scott", grid_size=180, log=True)

fig, ax = plt.subplots(figsize=(7.5, 7))
vmax = np.nanmax(np.abs(rr.density))
im = ax.imshow(rr.density, origin="lower", extent=rr.extent, cmap="RdBu_r", vmin=-vmax, vmax=vmax)
plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label="log relative risk (Assault vs. rest)")
city_fixed.boundary.plot(ax=ax, color="k", lw=0.4)
ax.set_aspect("equal")
ax.set_title("Where is Assault over- or under-represented\nrelative to crime overall?")
plt.tight_layout(); plt.show()

print(f"log relative risk range: [{np.nanmin(rr.density):.2f}, {np.nanmax(rr.density):.2f}]")
print("positive (red) = Assault more common here than the citywide mix would predict;"
      " negative (blue) = less common")

Assault (cases)      : n = 2,765
other crime (control): n = 11,183

log relative risk range: [-1.42, 0.81]
positive (red) = Assault more common here than the citywide mix would predict; negative (blue) = less common

Repeating this for a second type shows the surfaces genuinely differ by crime type, not just by population density (which would make every relative-risk map look the same).

tov = pp_all.by_mark("Theft of Vehicle")
other_tov = pp_all.subset(pp_all.marks != "Theft of Vehicle")
rr_tov = relative_risk(tov, other_tov, bandwidth="scott", grid_size=180, log=True)

fig, axes = plt.subplots(1, 2, figsize=(14, 6.4))
for ax, r, title in zip(axes, [rr, rr_tov], ["Assault vs. rest", "Theft of Vehicle vs. rest"]):
    vmax = np.nanmax(np.abs(r.density))
    im = ax.imshow(r.density, origin="lower", extent=r.extent, cmap="RdBu_r", vmin=-vmax, vmax=vmax)
    plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label="log relative risk")
    city_fixed.boundary.plot(ax=ax, color="k", lw=0.4)
    ax.set_aspect("equal"); ax.set_title(title)
plt.tight_layout(); plt.show()

corr = np.corrcoef(np.nan_to_num(rr.density.ravel()), np.nan_to_num(rr_tov.density.ravel()))[0, 1]
print(f"correlation between the two log-relative-risk surfaces: {corr:.3f}  "
      "(well short of 1: these are genuinely different spatial patterns, not both just 'more crime')")

correlation between the two log-relative-risk surfaces: 0.298  (well short of 1: these are genuinely different spatial patterns, not both just 'more crime')

---

10. Summary and Conclusions

What we did

1. Defined the intensity function $\lambda(\mathbf{s})$ and contrasted quadrat (binned, test-oriented) and kernel (smooth, mapping-oriented) estimation.
2. Ran the quadrat-count test at $8\times 8$ resolution: index of dispersion far above 1, and both the analytical and Monte Carlo tests reject CSR decisively — confirming Chapter 1.1's visual impression.
3. Checked sensitivity to grid resolution and flagged where the chi-square approximation becomes unreliable (small expected counts), while noting the qualitative conclusion (overdispersion) is robust throughout.
4. Estimated a kernel intensity surface with the package's uniform edge correction, and checked its integral against the known event count.
5. Compared Scott's-rule and cross-validation bandwidths, and showed they disagree by roughly a factor of four here, in the direction theory predicts (CV undersmooths, Scott's rule oversmooths).
6. Mapped relative risk for two crime types against the rest of crime, showing that "risk relative to the local baseline" is a different — and more actionable — map than raw density.

Practical takeaways

• A significant quadrat test does not say where intensity is elevated, only that it is — kernel density is the complementary tool for that.
• Bandwidth choice is not a detail: the same data can look like two dozen hot spots or two broad regions depending on $h$ alone. Report the bandwidth, and show more than one.
• A crime-type map that just reproduces the overall crime density map is not telling you anything new about that crime type — relative risk strips out the shared background.
• None of this chapter's methods use information about interaction between events (do assaults cluster near other assaults, beyond what the intensity surface predicts?) — that is Chapter 1.3's second-order question, logically separate from everything estimated here.

Next steps

Chapter 1.3 (Second-order and Distance-based Summary Functions) takes the interaction question head on: nearest-neighbour statistics, Ripley's K and L, and the pair correlation function, all of which ask whether events are closer together (or farther apart) than their intensity alone would predict.

---

11. Resources

Textbooks

Reference	Notes
Baddeley, A., Rubak, E. & Turner, R. (2015). Spatial Point Patterns: Methodology and Applications with R. CRC Press.	Chapters 6–7: quadrat counts and kernel intensity estimation
Silverman, B. W. (1986). Density Estimation for Statistics and Data Analysis. Chapman & Hall.	The classical reference on kernel density estimation and bandwidth selection
Diggle, P. J. (2013). Statistical Analysis of Spatial and Spatio-Temporal Point Patterns (3rd ed.). CRC Press.	Intensity estimation and relative risk (Ch. 4, 9)
Wand, M. P. & Jones, M. C. (1995). Kernel Smoothing. Chapman & Hall.	Edge correction and bandwidth selection in depth

Key papers

Paper	Contribution
Diggle, P. J. (1985). A kernel method for smoothing point process data. Applied Statistics, 34, 138–147.	Origin of kernel intensity estimation for point patterns
Scott, D. W. (1992). Multivariate Density Estimation. Wiley.	Scott's rule-of-thumb bandwidth
Kelsall, J. E. & Diggle, P. J. (1995). Non-parametric estimation of spatial variation in relative risk. Statistics in Medicine, 14, 2335–2342.	The relative-risk (case/control density ratio) method used in Section 9

In this library

Object	Role
spatialstats.pointpattern.quadrat_test	Chi-square and Monte Carlo test of CSR from binned counts
spatialstats.pointpattern.kde	Edge-corrected kernel intensity surface; KDEResult.integral, .at(), .plot()
spatialstats.pointpattern.bandwidth_scott, bandwidth_cv	Two bandwidth-selection rules
spatialstats.pointpattern.relative_risk	Case/control kernel density ratio
Chapter 1.1	Window, PointPattern, the CRS correction reused here
Chapter 1.3	Second-order (interaction) functions: ripley_k, ripley_l, pair_correlation