import sys, warnings
from pathlib import Path

import numpy as np
import pandas as pd
import geopandas as gpd
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.weights import Queen
from spatialstats.bayes import expected_counts, smr, BYM

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"}

counties = gpd.read_file(DATA / "diabetes_northeast.shp")
w = Queen(counties, transform="binary")
O = counties["Dia_count"].to_numpy()
POP = counties["POP_Total"].to_numpy()
E = expected_counts(POP, observed=O)
print(f"spatialstats : version {spatialstats.__version__}")

spatialstats : version 0.2.0

fit = BYM(O, E, w, model="bym", backend="gibbs", draws=2000, tune=2000, chains=4, thin=10, seed=1)

fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))
counties.assign(rr=fit.rr_mean).plot(column="rr", cmap="RdBu_r", ax=axes[0], edgecolor="k", linewidth=0.2,
                                     legend=True, legend_kwds={"label": "posterior mean RR", "shrink": 0.7})
axes[0].set_axis_off()
axes[0].set_title("Point Estimate Alone")

width = fit.rr_hi - fit.rr_lo
counties.assign(width=width).plot(column="width", cmap="viridis", ax=axes[1], edgecolor="k", linewidth=0.2,
                                  legend=True, legend_kwds={"label": "95% CI width", "shrink": 0.7})
axes[1].set_axis_off()
axes[1].set_title("Credible-interval Width (the missing half of the story)")
plt.tight_layout()
plt.show()

print(f"CI width range: [{width.min():.3f}, {width.max():.3f}]  "
      f"(narrowest-to-widest ratio: {width.max() / width.min():.1f}x)")

CI width range: [0.009, 0.228]  (narrowest-to-widest ratio: 26.2x)

rng = np.random.default_rng(1)
O_thin = rng.binomial(O.astype(int), 0.02)

suppress_threshold = 10
to_suppress = O_thin < suppress_threshold
print(f"counties with a thinned count below {suppress_threshold}: {int(to_suppress.sum())} of {len(O_thin)}")
counties.loc[to_suppress, ["COUNTY", "STATE"]].assign(thinned_count=O_thin[to_suppress])

counties with a thinned count below 10: 3 of 217

rng = np.random.default_rng(0)
# Two synthetic "counties." Within each, x and y are INDEPENDENT (zero individual-level correlation).
county_A_x = rng.normal(30, 5, 2000)    # e.g. individual obesity measure
county_A_y = rng.normal(5, 1, 2000)     # e.g. individual diabetes indicator, unrelated to x within the county
county_B_x = rng.normal(45, 5, 2000)    # a county with higher average obesity...
county_B_y = rng.normal(9, 1, 2000)     # ...and higher average diabetes, but STILL unrelated to x within it

within_corr_A = np.corrcoef(county_A_x, county_A_y)[0, 1]
within_corr_B = np.corrcoef(county_B_x, county_B_y)[0, 1]
print(f"within-county A correlation (individual level): {within_corr_A:.3f}")
print(f"within-county B correlation (individual level): {within_corr_B:.3f}")

county_means_x = [county_A_x.mean(), county_B_x.mean()]
county_means_y = [county_A_y.mean(), county_B_y.mean()]
between_corr = np.corrcoef(county_means_x, county_means_y)[0, 1]
print(f"between-county correlation (what an aggregate map would show): {between_corr:.3f}")

within-county A correlation (individual level): -0.017
within-county B correlation (individual level): 0.021
between-county correlation (what an aggregate map would show): 1.000

Step	What we do
Uncertainty	Why a smoothed point estimate alone, without its interval, misrepresents what the model actually knows
Small-cell suppression	A concrete demonstration of when and why raw counts should not be published at all
Ecological fallacy	A worked numerical example of why a county-level association cannot be read as an individual-level one
Communicating risk	Concrete do's and don'ts for presenting smoothed disease-mapping results

	COUNTY	STATE	thinned_count
89	Hamilton County	New York	9
142	Cameron County	Pennsylvania	6
209	Grand Isle County	Vermont	8

Reference	Notes
Lawson, A. B. (2018). Bayesian Disease Mapping (3rd ed.). CRC Press.	Chapter 9 covers reporting and communication standards for disease-mapping results
Robinson, W. S. (1950). Ecological correlations and the behavior of individuals. American Sociological Review, 15, 351–357.	The original ecological fallacy paper, and origin of the numerical structure used in Section 5

Paper	Contribution
Elliott, P. & Wartenberg, D. (2004). Spatial epidemiology: current approaches and future challenges. Environmental Health Perspectives, 112, 998–1006.	Communication and ethical standards for spatial epidemiology
Rushton, G. et al. (2006). Geocoding in cancer research: a review. American Journal of Preventive Medicine, 30, S16–S24.	Small-cell disclosure risk in geocoded health data

Object	Role
`BYM.rr_lo`, `.rr_hi`	The credible-interval bounds used in Section 3
`BYM.exceedance()`	The alternative, uncorrected uncertainty summary from Chapter 5.6
Chapter 5.11	Closing hand-offs and exercises

5.10 Reporting, Interpretation and Ethics¶

Table of Contents¶

1. Overview¶

2. Python Setup¶

3. Report Uncertainty, Not Just a Point Estimate¶

4. Small-cell Suppression¶

5. The Ecological Fallacy¶

6. Communicating Risk Responsibly¶

7. Summary and Conclusions¶

What we did¶

Practical takeaways¶

Next steps¶

8. Resources¶

Textbooks¶

Key papers¶

In this library¶