8.3 Estimation and Variance

Part: 8 — Spatial Sampling and Design (Chapter 8.3 of 5)

Dataset: diabetes_usa.csv (all 3,107 US counties — repeated sampling against this known population)

Focus: estimate_mean's design-matched variance; evaluate_design / compare_designs; the full reference table, reproduced

Library: spatialstats.sampling (PySpatialStats)

---

Table of Contents

1. Overview

2. Python Setup

3. Variance Depends on the Design

4. Simulating Against a Known Population

5. The Full Comparison, 800 Replicates

6. Does Spatial Balance Only Help When the Variable Is Spatially Structured?

7. PPS-GRTS: Unbiased Totals with Unequal Probabilities

8. Summary and Conclusions

9. Resources

---

1. Overview

Chapter 8.2 ranked designs by a purely geometric spatial-balance score. This chapter asks the question that actually matters: does a better balance score translate into a smaller true sampling error? Because diabetes_usa.csv is a genuine census, this is not a theoretical question — it can be answered directly, by drawing hundreds of samples under each design and comparing the resulting estimates to the one number that is never in doubt: the true population mean.

Step	What we do
Design-matched variance	estimate_mean's variance='auto' picks SRS, stratified, local-mean, or ultimate-cluster variance automatically
Simulate	evaluate_design / compare_designs: hundreds of replicate samples, each design, against the known truth
The full table	Monte Carlo s.d., design effect, spatial balance, CI coverage — this Part's own reference table, reproduced live
A necessary condition	Spatial balance only pays off when the outcome variable is itself spatially structured — shown directly, not just asserted

---

2. Python Setup

import sys, warnings
from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

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.sampling import (
    simple_random, systematic, stratified, two_stage, grts, estimate_mean, evaluate_design, compare_designs
)

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,
})
COLORS = {"blue": "#2b6cb0", "red": "#c0392b", "green": "#2f855a", "orange": "#dd6b20",
          "grey": "#718096", "purple": "#6b46c1", "teal": "#2c7a7b"}

usa = pd.read_csv(DATA / "diabetes_usa.csv")
frame = usa[["X", "Y"]]
y = usa["Diabetes_per"].to_numpy()
true_mean = y.mean()
print(f"spatialstats : version {spatialstats.__version__}")
print(f"true population mean: {true_mean:.4f}")

spatialstats : version 0.2.0
true population mean: 8.5966

---

3. Variance Depends on the Design

estimate_mean(y, sample, variance='auto') picks the correct variance formula from the sample's own design metadata — an SRS-style formula for simple random samples, a stratified formula for stratified samples, local-mean variance (Section 4) for spatially balanced or systematic samples (using the sample's own geometry: neighbouring selected units are similar when $y$ is spatially structured, which lets the variance formula exploit that structure directly), and an ultimate-cluster/ratio-linearised formula for two-stage samples.

srs = simple_random(frame, n=100, seed=1)
sys_hilbert = systematic(frame, n=100, order="hilbert", seed=1)
two_stage_sample = two_stage(frame, clusters=usa["State"].to_numpy(), n_clusters=12, n_per_cluster=8, seed=1)

for name, s in [("SRS", srs), ("Systematic (Hilbert)", sys_hilbert), ("Two-stage 12x8", two_stage_sample)]:
    est = estimate_mean(y, s, variance="auto")
    print(f"{name:22s}: mean={est.mean:.3f}  se={est.se:.3f}  variance_method={est.variance_method}")

SRS                   : mean=8.720  se=0.171  variance_method=srs
Systematic (Hilbert)  : mean=8.537  se=0.118  variance_method=local
Two-stage 12x8        : mean=8.220  se=0.328  variance_method=cluster

Each design automatically receives the variance formula appropriate to how it was drawn — estimate_mean never silently applies an SRS-style formula to a systematically or cluster-sampled dataset, which would misstate the true uncertainty in either direction depending on the design.

---

4. Simulating Against a Known Population

evaluate_design(values, design_fn, n_reps, seed) draws n_reps independent replicate samples under a design (supplied as a function of a random generator), computes estimate_mean for each, and returns the realized distribution of estimates — from which the true Monte Carlo standard deviation, coverage, and design effect can all be computed directly, rather than relying on any single sample's own reported standard error.

srs_eval = evaluate_design(y, lambda rng: simple_random(frame, n=100, seed=rng), n_reps=200, seed=42)
print(srs_eval)

truth         8.596611
bias         -0.011802
sd            0.155794
rmse          0.155851
mean_se       0.152920
coverage      0.940000
mean_width    0.599437
balance       0.343245
dtype: float64

---

5. The Full Comparison, 800 Replicates

Every design from Chapter 8.2, compared on identical terms — compare_designs runs evaluate_design for each one and reports Monte Carlo standard deviation, design effect (relative to whichever design serves as the baseline in this same call), coverage, and spatial balance in one table. 800 replicates matches this Part's own reference validation.

strata_urban_rural = usa["Urban_Rural"].to_numpy()
designs = {
    "SRS": lambda rng: simple_random(frame, n=100, seed=rng),
    "Systematic (Hilbert)": lambda rng: systematic(frame, n=100, order="hilbert", seed=rng),
    "GRTS": lambda rng: grts(frame, n=100, seed=rng),
    "Stratified+GRTS": lambda rng: stratified(frame, strata=strata_urban_rural, n=100, allocation="proportional",
                                              within="grts", seed=rng),
    "Two-stage 12x8": lambda rng: two_stage(frame, clusters=usa["State"].to_numpy(), n_clusters=12,
                                            n_per_cluster=8, seed=rng),
    "Two-stage 25x4": lambda rng: two_stage(frame, clusters=usa["State"].to_numpy(), n_clusters=25,
                                            n_per_cluster=4, seed=rng),
}
results = compare_designs(y, designs, n_reps=800, seed=42)
results.round(3)

	truth	bias	sd	rmse	mean_se	coverage	mean_width	balance	deff
SRS	8.597	-0.005	0.153	0.153	0.154	0.949	0.604	0.334	1.000
Systematic (Hilbert)	8.597	0.001	0.118	0.118	0.128	0.972	0.501	0.081	0.600
GRTS	8.597	-0.000	0.125	0.124	0.126	0.945	0.492	0.114	0.664
Stratified+GRTS	8.597	-0.002	0.124	0.124	0.151	0.975	0.594	0.145	0.659
Two-stage 12x8	8.597	0.012	0.310	0.310	0.324	0.938	1.271	1.920	4.122
Two-stage 25x4	8.597	0.012	0.208	0.208	0.257	0.979	1.008	0.657	1.848

fig, axes = plt.subplots(1, 2, figsize=(13, 5))
order = results.sort_values("sd").index
axes[0].barh(order, results.loc[order, "sd"], color=COLORS["blue"])
axes[0].set_xlabel("Monte Carlo s.d.")
axes[0].set_title("Realized sampling variability, 800 replicates")
axes[0].invert_yaxis()

axes[1].barh(order, results.loc[order, "coverage"] * 100, color=COLORS["green"])
axes[1].axvline(95, color=COLORS["red"], linestyle="--", label="nominal 95%")
axes[1].set_xlabel("95% CI coverage (%)")
axes[1].set_title("Confidence-interval coverage")
axes[1].invert_yaxis()
axes[1].legend()
plt.tight_layout()
plt.savefig("design_comparison_800reps.png", dpi=100, bbox_inches="tight")
plt.show()

Systematic (Hilbert) and GRTS both cut Monte Carlo standard deviation substantially below SRS, with genuinely valid (at or above nominal) confidence-interval coverage — not just a narrower interval that happens to be wrong more often. Both two-stage designs show design effects above 1: fewer independent regions visited means less genuinely new spatial information per county sampled, even though the raw count of counties sampled ($n=91$–100) looks similar to every other design in this table. This is the full, live reproduction of this Part's own reference validation table.

---

6. Does Spatial Balance Only Help When the Variable Is Spatially Structured?

Every design so far has been evaluated on Diabetes_per, a variable with real spatial structure (Chapter 8.1, Section 3's map). If GRTS's advantage comes specifically from exploiting that structure, it should disappear on a variable with no spatial structure at all.

rng_noise = np.random.default_rng(0)
y_random = rng_noise.normal(loc=y.mean(), scale=y.std(), size=len(y))   # same mean/sd, but spatially shuffled

designs_noise = {
    "SRS": lambda rng: simple_random(frame, n=100, seed=rng),
    "GRTS": lambda rng: grts(frame, n=100, seed=rng),
}
results_noise = compare_designs(y_random, designs_noise, n_reps=800, seed=42)
print("On a spatially UNSTRUCTURED variable (same mean/sd as Diabetes_per, but randomly shuffled):")
print(results_noise[["sd", "deff", "coverage"]].round(3))

On a spatially UNSTRUCTURED variable (same mean/sd as Diabetes_per, but randomly shuffled):
         sd   deff  coverage
SRS   0.153  1.000     0.955
GRTS  0.148  0.931     0.958

Confirmed directly: GRTS's design effect collapses to essentially 1.0 on the spatially unstructured version of the same variable — no efficiency gain at all, though coverage remains correctly calibrated either way. This is not a limitation of GRTS; it is the mechanism working exactly as advertised. Spatial balance improves estimation precisely by exploiting spatial structure in $y$ — with none to exploit, there is nothing for a spatially balanced sample to gain over a random one.

---

7. PPS-GRTS: Unbiased Totals with Unequal Probabilities

PPS sampling (size-weighted inclusion probabilities) deliberately oversamples large-population counties. The Horvitz-Thompson total should still be unbiased — a claim worth checking directly, by simulation, rather than trusting the algebra alone.

pop = usa["POP_Total"].to_numpy()
true_total = y.sum()

def pps_grts(rng):
    return grts(frame, n=100, size=pop, seed=rng)

totals = []
for _ in range(300):
    s = pps_grts(np.random.default_rng())
    est = estimate_mean(y, s)
    totals.append(est.total)
totals = np.array(totals)

rel_bias = (totals.mean() - true_total) / true_total * 100
print(f"true total (sum of Diabetes_per across all counties): {true_total:,.1f}")
print(f"mean Horvitz-Thompson total, 300 PPS-GRTS replicates: {totals.mean():,.1f}")
print(f"relative bias: {rel_bias:+.2f}%")

true total (sum of Diabetes_per across all counties): 26,709.7
mean Horvitz-Thompson total, 300 PPS-GRTS replicates: 26,552.9
relative bias: -0.59%

The mean Horvitz-Thompson total across replicates lands within a couple of percent of the true total (printed above) — small relative to the sample-to-sample spread (Section 5's design effects), confirming this library's PPS-GRTS implementation delivers on the unbiasedness guarantee unequal-probability sampling promises, not merely a plausible-looking approximation. The remaining gap reflects Monte Carlo noise from a moderate replicate count, not a systematic flaw — it should shrink further with more replicates, exactly as an unbiased estimator should.

---

8. Summary and Conclusions

This chapter turned Chapter 8.2's geometric spatial-balance ranking into measured statistical efficiency, against a population where the truth is always known.

What we did

1. Confirmed estimate_mean automatically matches its variance formula to each sample's own design.
2. Ran evaluate_design/compare_designs with 800 replicates — this Part's own reference validation, reproduced live: GRTS and systematic sampling substantially reduce Monte Carlo standard deviation versus SRS with valid CI coverage; two-stage designs show design effects above 1.
3. Confirmed directly that GRTS's efficiency gain vanishes on a spatially unstructured variable — design effect collapses to $\approx 1$ once the same variable is randomly shuffled across counties, isolating exactly why spatial balance helps.
4. Confirmed PPS-GRTS's Horvitz-Thompson total is unbiased to within a couple of percent by direct simulation, not just algebraic assertion.

Practical takeaways

• Always let estimate_mean pick the variance formula (variance="auto") rather than assuming an SRS-style formula applies to any design — using the wrong one can badly misstate a confidence interval's actual coverage.
• Before expecting a spatially balanced design to outperform SRS, check whether the target variable is actually spatially structured; if it is not, the extra design complexity buys nothing.
• A design effect measured once, honestly, by simulation (as in this chapter) is far more trustworthy than one assumed from a textbook table — real designs on real populations can differ from idealized theory.

Next steps

Chapter 8.4 (Monitoring-Network Design) applies these same tools to a genuinely unsolved, practical problem: where should new PM2.5 monitors go, given the gaps in California's existing 159-site network from Part 7?

---

9. Resources

Textbooks

Reference	Notes
Särndal, C.-E., Swensson, B. & Wretman, J. (1992). op. cit. (Chapter 8.1)	The design-matched variance formulas behind Section 3

Key papers

Paper	Contribution
Stevens, D. L. & Olsen, A. R. (2003). Variance estimation for spatially balanced samples of environmental resources. Environmetrics, 14, 593–610.	The local-mean variance estimator used for GRTS/systematic samples in Section 3

In this library

Object	Role
spatialstats.sampling.estimate_mean	Section 3
spatialstats.sampling.evaluate_design, compare_designs	Sections 4–7, the core of this chapter
Chapter 8.2	The spatial-balance ranking this chapter turns into measured efficiency