17. Chapter D. Spatial Statistics#
17.1. Learning Objectives#
After completing this chapter, you should be able to:
Distinguish CSR（完全空间随机）from conditional random labeling.
Simulate homogeneous Poisson event patterns and aggregate occupancy in disks.
Derive positive covariance caused by overlapping spatial supports.
Show why shared events invalidate unrestricted transaction-label permutations.
Contrast equal disks with uneven event intensity and point-process references.
Connection to STARM: The examples mirror the binary outcome, fixed antecedent indicators, lagged transactions, and overlapping spatial supports in the STARM manuscript. They are teaching examples, not reproductions of manuscript results.
Teaching philosophy: start with the general statistical problem, build mathematical intuition, and use STARM as one illustrative case study. The new sections precede the original computational labs.
17.2. Why Study This? — Motivation, Scope, and Real-World Problems#
Ordinary statistics often assumes independent observations. Geographical data violate this because nearby locations share environments, interactions, and even the same events. Spatial statistics (空间统计学) studies point locations, geographic surfaces, spatial dependence, and uncertainty while respecting where observations occur. Applications include epidemiology, ecology, crime analysis, geosciences, and public-health surveillance.
Questions it answers: Is an observed cluster unusual? How should we model event intensity? What does proximity do to variance and inference? Does the spatial support or aggregation unit change a result?
17.3. Conceptual Roadmap#
Scientific problem → statistical representation → assumptions → mathematics → computational experiment → diagnostics → interpretation.
The chapters are standalone: STARM is a later application, not the reason these mathematical ideas exist.
17.4. Foundations and Mathematical Derivations#
17.4.1. Three representations, three statistical questions#
Point patterns (点模式) model the event locations \(\{s_j\}\); e.g. disease cases or ignitions.
Geostatistical fields (地统计场) model a continuous \(Z(s)\); e.g. temperature, pollution. Spatial covariance is \(C(h)=\mathrm{Cov}[Z(s),Z(s+h)]\) for a stationary field.
Areal/lattice data (区域数据) summarize attributes in polygons; spatial neighbors can be encoded through \(W=[w_{ij}]\).
17.4.2. CSR, intensity, and Ripley’s K#
Under homogeneous complete spatial randomness (CSR, 完全空间随机性), events follow a Poisson process with intensity \(\lambda\): for region \(A\), \(N(A)\sim\mathrm{Poisson}(\lambda|A|)\), and counts on disjoint regions are independent. The void probability is \(P[N(A)=0]=e^{-\lambda|A|}\). Under CSR (ignoring edge issues), the theoretical \(K\)-function equals \(K(r)=\pi r^2\). Observed clustering can also reflect inhomogeneous intensity, not genuine interaction.
17.4.3. Why overlapping supports induce covariance#
Let \(Y_A=1\{N(A)>0\}\) and \(Y_B=1\{N(B)>0\}\) under homogeneous Poisson CSR. Since
and \(P(Y_A=0)=e^{-\lambda|A|}\), \(P(Y_B=0)=e^{-\lambda|B|}\), it follows by expanding \(\mathrm{Cov}(Y_A,Y_B)\) that
Because \(|A\cup B|=|A|+|B|-|A\cap B|\), this is positive for overlapping supports of positive area. This is dependence even when points follow CSR.
17.4.4. Moran’s I and why aggregation matters#
For values \(z_i\) on areas and spatial weights \(w_{ij}\), Moran’s \(I\) is
Interpretation depends on \(W\) and the null used to calibrate it. The MAUP (可变空间单元问题) means different zoning or aggregation scales can alter statistical summaries and relationships.
17.4.5. Inference needs a spatial reference#
Random-label models, spatial point-process models, Gaussian random fields, block resampling and spatial regression each condition on different features. No one reference is universal. Use a model matched to the scientific question and the observed dependence structure.
17.4.6. Critical thinking checkpoint#
Before running code, write down which observations are assumed independent or exchangeable and explain why. Revisit this statement after inspecting the visualizations.
17.4.7. Why spatial hotspot maps need multiplicity control#
A spatial analyst may test every county, every disk center, several radii, outcomes, and time windows. Even if each valid local test uses \(\alpha=0.05\), displaying every unadjusted “significant” location creates a multiple-comparisons problem. Overlapping polygons or disks also reuse events, so their p-values are dependent even under complete spatial randomness. This dependence is a consequence of the observation geometry, not proof of a causal spatial process.
For a map-wide family of \(M\) tests, \(V\) counts false hotspot declarations, \(R\) all declarations, and \(S=R-V\) true ones. The FDP is the realized share \(V/\max(R,1)\); its repeated-study expectation is FDR; FWER is the probability \(P(V\ge1)\) of any false hotspot. Under the complete null, FDR=FWER; when some hotspots are real, FDR can be much smaller. Ask what a false declaration costs before choosing which risk to control.
17.4.8. How the methods arose and what their guarantees mean here#
Bonferroni: many nominal tests make at least one false alarm likely. The union bound allocates \(\alpha/M\) per valid test, ensuring FWER \(\le\alpha\) under arbitrary dependence. Use it for a small prespecified set or when any false hotspot could trigger a costly intervention. It can be conservative for thousands of map cells.
BH: Benjamini and Hochberg (1995) shifted the target from “no false result anywhere” to the expected false share of the reported list, seeking more power when true spatial signals exist. Sort p-values and reject through the largest rank \(k\) satisfying \(p_{(k)}\le kq/M\). The guarantee needs independent valid p-values or the appropriate positive-dependence condition; overlapping supports do not automatically satisfy that condition.
BY: Benjamini and Yekutieli (2001) made the BH thresholds more stringent by dividing \(q\) by \(H_M=\sum_{j=1}^M1/j\). BY controls FDR under arbitrary dependence among valid p-values, at a power cost. It does not rescue p-values produced by an invalid spatial null.
17.4.9. A practical spatial workflow#
Define the family: all locations, radii, outcomes, and scales inspected—not only the final map layer.
Choose a spatial null that preserves the features required by the question: CSR for homogeneous point locations, an inhomogeneous point-process model for varying intensity, or a dependence-preserving permutation/block design when justified.
Check marginal p-value calibration under that null before adjustment. For overlapping disks, compare row shuffling with full point-pattern simulation; the two answer different questions.
Choose FWER for confirmatory, high-consequence hotspot declarations; choose FDR for a clearly labeled exploratory map. Use BH only when its dependence conditions are credible, BY for arbitrary dependence of valid p-values, or a spatial max-statistic procedure when a design-valid joint null is available.
Map adjusted q-values or corrected decisions alongside effect sizes and uncertainty. Do not interpret a corrected p-value as a measure of practical importance or causation.
Historical anchors: Benjamini & Hochberg (1995), FDR as a practical alternative to FWER; Benjamini & Yekutieli (2001), FDR under dependency; NIST, Bonferroni inequality and simultaneous comparisons.
17.5. Additional Derivation Visualization — General-Purpose Example#
The following self-contained experiment illustrates a general statistical principle, not STARM-specific results. Check how the plot corresponds to the formulas above.
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist
rng=np.random.default_rng(77)
pts=rng.uniform(size=(150,2))
fig,axs=plt.subplots(1,2,figsize=(10,4))
axs[0].scatter(*pts.T,s=13);axs[0].set(xlim=(0,1),ylim=(0,1),title='Homogeneous Poisson points conditional on N')
xy=rng.uniform(size=(20,2))
Z=np.sin(2*np.pi*xy[:,0])+np.sin(2*np.pi*xy[:,1])
D=cdist(xy,xy);W=((D>0)&(D<.3)).astype(float);S=W.sum()
center=Z-Z.mean();I=len(Z)/S*(center@W@center)/(center@center)
axs[1].scatter(*xy.T,c=Z,s=60);axs[1].set(xlim=(0,1),ylim=(0,1),title=f'A spatial field | Moran I={I:.2f}')
for a in axs:a.set_aspect('equal')
plt.tight_layout();plt.show()
print('Illustrative Moran I depends on chosen neighborhood distance and null.')
Illustrative Moran I depends on chosen neighborhood distance and null.
flowchart LR
A[Simulate spatial point pattern] --> B[Place sample anchors]
B --> C[Build disk supports]
C --> D[Count events per disk]
D --> E[Compare overlap dependence]
E --> F[Select justified spatial null]
This chapter follows the same structure as the supplied MCLP tutorial: conceptual workflow, formal definition, assumptions, a hand-worked example, executable Python experiments, sensitivity checks, common pitfalls, and a quiz.
17.6. Formal Definition and Core Equations#
17.6.1. Notation#
Symbol |
Meaning |
|---|---|
\(N(A)\) |
point count in spatial set \(A\) |
\(\lambda\) |
homogeneous Poisson intensity |
\(Y_A=\mathbf{1}\{N(A)>0\}\) |
occupied-disk indicator |
$ |
A\cap B |
For disjoint sets under a homogeneous Poisson process, event counts are independent. For overlapping sets, common events induce positive dependence even under CSR.
17.7. Core Mechanisms and Concepts#
17.7.1. Spatial CSR does not mean independent disk occupancies#
A Poisson pattern is generated in continuous geography. When the same event lies inside two sampling disks, those two transaction outcomes are logically linked.
17.7.2. Three distinct null constructions#
CSR: simulate events over the spatial domain, then aggregate to supports.
Point-mark random labeling: hold pooled event locations fixed while permuting their marks (when appropriate).
Transaction-row shuffling: hold disk-level covariates fixed and reallocate occupancy indicators; this imposes a separate label exchangeability hypothesis.
17.7.3. Geography of support#
Radius affects area and overlap, and hence occupancy prevalence and dependence. A sample of many disks does not automatically provide that many independent observations.
17.8. Assumptions, Strengths, Limitations, and Trade-offs#
Assume homogeneous Poisson intensity and correctly specified disk geometry for the analytic covariance. Strength: exact overlap-based covariance. Limitation: wildfire has seasonal variation and spatially varying intensities, and event processes may be clustered beyond Poisson.
Trade-off: larger disks may collect more events but increase overlap and reduce the effective amount of independent information.
17.9. Worked Toy Example#
Take two equal disks with radius \(r\). They are disjoint when center separation \(d\ge2r\). If they intersect, the intersection area is positive and so is covariance of occupancy. This dependence appears even if events follow perfect homogeneous CSR and the antecedent does not affect intensity.
17.10. Python Lab: Set Up the Dataset#
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import hypergeom, binom
SEED = 41020001
rng = np.random.default_rng(SEED)
print("Reproducible random seed:", SEED)
Reproducible random seed: 41020001
17.11. Simulate points and disks#
Inspect event points and supporting disks in a 100×100 square.
from matplotlib.patches import Circle
L=100;r=9;lam=.006
anchors=np.array([[25,50],[39,50],[65,50],[85,50]])
npnts=rng.poisson(lam*L*L);points=rng.uniform(0,L,size=(npnts,2))
fig,ax=plt.subplots(figsize=(8,4));ax.scatter(points[:,0],points[:,1],s=8,alpha=.6)
for a in anchors:ax.add_patch(Circle(a,r,fill=False,linewidth=2))
ax.scatter(anchors[:,0],anchors[:,1],marker='x',s=80,color='black');ax.set(xlim=(0,L),ylim=(0,L),aspect='equal',title='Poisson points and sampled disk supports');plt.show()
Interpretation. The first two supports overlap; points inside their shared area contribute to both outcomes.
17.12. Derive and plot overlap covariance#
Use the exact equal-disk intersection formula.
def intersection_area(r,d):
d=np.asarray(d,float);a=np.zeros_like(d);z=d<2*r
a[z]=2*r*r*np.arccos(d[z]/(2*r))-.5*d[z]*np.sqrt(np.maximum(0,4*r*r-d[z]**2))
return a
def occupancy_cov(r,d,lam):
area=np.pi*r*r;v=intersection_area(r,d)
return np.exp(-2*lam*area)*np.expm1(lam*v)
d=np.linspace(0,2.8*r,160)
fig,ax=plt.subplots(figsize=(7,4));ax.plot(d,occupancy_cov(r,d,lam));ax.axvline(2*r,linestyle='--',color='black');ax.set(xlabel='Disk-center separation',ylabel='Covariance of occupancy',title='Poisson overlap covariance');plt.show()
Interpretation. For separated disks the covariance goes to zero; overlapping occupancy indicators remain correlated even with CSR.
17.13. Estimate covariance from repeated Poisson realizations#
Repeatedly generate event points and compute two disk outcomes.
reps=6000;pair=anchors[:2];draws=np.empty((reps,2))
for b in range(reps):
p=rng.uniform(0,L,size=(rng.poisson(lam*L*L),2))
draws[b]=[((p-a)**2).sum(axis=1).min(initial=np.inf)<=r*r for a in pair]
emp=np.cov(draws,rowvar=False,ddof=0)[0,1]
analytic=occupancy_cov(r,np.linalg.norm(pair[0]-pair[1]),lam).item()
print({'empirical_covariance':emp,'analytic_covariance':analytic})
{'empirical_covariance': np.float64(0.009701333333333343), 'analytic_covariance': 0.009615884832970982}
Interpretation. A shared-event covariance is a property of the data-generating mechanism, not of an association between X and Y.
17.14. Compare expected occupancy under different intensities#
Uneven baseline intensity creates different Bernoulli margins even without shared events.
area=np.pi*r*r
print(pd.DataFrame({'intensity':[.001,.003,.006,.01], 'P_disk_occupied':[1-np.exp(-l*area) for l in [.001,.003,.006,.01]]}).round(3))
intensity P_disk_occupied
0 0.001 0.225
1 0.003 0.534
2 0.006 0.783
3 0.010 0.922
Interpretation. Unequal margins violate uniform allocation of a fixed total unless the null design accounts for them.
17.15. Common Pitfalls and Misunderstandings#
CSR is a statement about spatial event generation, not automatic exchangeability of disk-level binary outcomes.
Nonoverlapping disks can still be nonexchangeable under spatial heterogeneity.
Point-level label permutation and row-level outcome permutation test different nulls.
A raster neighborhood summary uses geometric support; this differs from rule support (transaction frequency).
Increasing Monte Carlo iterations does not remove overlap-induced dependence.
17.17. Suggested Reading and STARM Crosswalk#
Agrawal & Srikant (1994): frequent-itemset search (background).
Hope (1968); Phipson & Smyth (2010): randomization/Monte Carlo p-values (Chapter B).
Benjamini & Hochberg (1995); Benjamini & Yekutieli (2001): multiple testing and dependence (Chapter C).
Baddeley, Rubak & Turner (2015): spatial point patterns (Chapter D).
Winkler et al. (2014): exchangeability blocks and permutation designs (Chapters B, D, E).
STARM manuscript: Sections 3.3–3.6 for reference assumptions; Section 4 for simulations; Sections 5–6 for empirical limitations.
Practice task: Write a short reviewer note distinguishing a descriptive association from an inferential discovery, explicitly naming the null model and its assumptions.