1"""H2 — Regional assortment diversity vs regional scale (n <= 17; preregistered: Spearman rho with CI; OLS HC3 secondary),2plus H7-V3 (incremental information of the ECI beyond population, business count and income; leave-one-region-out).3External covariates: ISQ population 2025p, ISQ business locations with employees Dec 2025, ISQ disposable income per capita 2024p,4Census 2021 land area (density) — see data/external/README.md. Output experiments/h2_scale.json"""5from __future__ import annotations6import json7import numpy as np, pandas as pd8from scipy import stats9from src.common import ROOT, SEED, get_logger, ledger, write_json, PROC_DIR, EXP_DIR, P1011log = get_logger("models.h2")12PROC = PROC_DIR13COV = {"log_population": ("population", "ISQ population estimate, 1 July 2025 (provisional)"), "log_business_locations": ("business_locations_with_employees", "ISQ business locations with employees, Dec 2025"),14 "disposable_income": ("disposable_income_per_capita_cad", "ISQ disposable income per capita, 2024p"), "log_density": ("population_density_per_km2", "Census 2021 population density (2021 ER boundaries)")}151617def spearman_ci(x, y, rng, B=5000, boot_y=None):18 """Region-level bootstrap combined with outcome resampling draws (boot_y: array B x n of outcome draws, optional)."""19 n = len(x); rhos = []20 for b in range(B):21 idx = rng.integers(0, n, n); yy = y if boot_y is None else boot_y[b % boot_y.shape[0]]22 if len(set(idx)) < 4: continue23 rhos.append(stats.spearmanr(x[idx], yy[idx]).correlation)24 rhos = np.array(rhos); return float(stats.spearmanr(x, y).correlation), [float(np.nanpercentile(rhos, 2.5)), float(np.nanpercentile(rhos, 97.5))]252627def ols_hc3(x, y):28 X = np.column_stack([np.ones_like(x), x]); b = np.linalg.lstsq(X, y, rcond=None)[0]; e = y - X @ b; h = np.diag(X @ np.linalg.pinv(X.T @ X) @ X.T)29 V = np.linalg.pinv(X.T @ X) @ (X.T * (e ** 2 / (1 - h) ** 2)) @ X @ np.linalg.pinv(X.T @ X); return float(b[1]), float(np.sqrt(V[1, 1]))303132def main():33 rng = np.random.default_rng(SEED)34 ext = pd.read_csv(ROOT / "data" / "external" / "regions_external.csv").set_index("region_name")35 div = pd.read_csv(PROC / "diversity_region.csv").set_index("region"); ind = pd.read_csv(PROC / "h7_region_indicators.csv").set_index("region")36 sp = pd.read_parquet(PROC / "store_panel.parquet"); cnt = sp[sp.region_final != "Inconnue"].groupby("region_final").size()37 df = div.join(ext, how="inner"); df["n_stores"] = cnt.reindex(df.index)38 for k, (col, _) in COV.items():39 df[k] = np.log(df[col]) if k.startswith("log") else df[col]40 out = {"regions_in_external": int(len(ext)), "regions_matched": int(len(df)), "correlations": {}, "covariate_sources": {k: v[1] for k, v in COV.items()},41 "note": "17 regions; rarefied diversity needs >= 15 stores (Nord-du-Québec, 3 stores, excluded). Associations only; CIs from 5 000 region-level bootstrap draws combined with the 2 000 rarefaction draws of the outcome; OLS slope with HC3 SE. Power: with n = 16 the smallest |rho| detectable with 80 % power at alpha = 0.05 is about 0.66.", "power_n16_rho_80pct": 0.66}42 drp = PROC / "diversity_rarefied_draws.parquet"; draws = pd.read_parquet(drp) if drp.exists() else None43 for outc in ("H_rarefied15", "H_norm_M", "eff_cat_M"):44 for k in COV:45 d = df[[outc, k]].dropna(); x = d[k].values.astype(float); y = d[outc].values.astype(float)46 if len(d) < 6: continue47 by = draws[list(d.index)].values if (outc == "H_rarefied15" and draws is not None and all(r in draws.columns for r in d.index)) else None48 rho, ci = spearman_ci(x, y, rng, boot_y=by); b, se = ols_hc3(x, y)49 out["correlations"][f"{outc} ~ {k}"] = {"n": int(len(d)), "rho": round(rho, 3), "ci": [round(ci[0], 3), round(ci[1], 3)], "ols_b": round(b, 4), "ols_se_hc3": round(se, 4)}50 # stores per capita vs scale (coverage gradient)51 df["log_stores_per_100k"] = np.log(df.n_stores / df.population * 1e5)52 for k in ("log_population", "log_density", "disposable_income"):53 d = df[["log_stores_per_100k", k]].dropna(); rho, ci = spearman_ci(d[k].values.astype(float), d.log_stores_per_100k.values, rng)54 out["correlations"][f"log_stores_per_100k ~ {k}"] = {"n": int(len(d)), "rho": round(rho, 3), "ci": [round(ci[0], 3), round(ci[1], 3)], "ols_b": None, "ols_se_hc3": None}55 # H7-V3: ECI on covariates, leave-one-region-out predictive correlation56 e = ind.join(df[["log_population", "log_business_locations", "disposable_income"]], how="inner").dropna(subset=["eci"])57 X = np.column_stack([np.ones(len(e)), e.log_population, e.log_business_locations, e.disposable_income]); y = e.eci.values58 b = np.linalg.lstsq(X, y, rcond=None)[0]; r2 = 1 - ((y - X @ b) ** 2).sum() / ((y - y.mean()) ** 2).sum()59 pred = np.array([(X[i] @ np.linalg.lstsq(np.delete(X, i, 0), np.delete(y, i), rcond=None)[0]) for i in range(len(y))])60 v3 = {"n": int(len(y)), "R2_eci_on_covariates": round(float(r2), 3), "loo_pred_spearman": round(float(stats.spearmanr(pred, y).correlation), 3), "criterion_rho_ge_0.5": bool(stats.spearmanr(pred, y).correlation >= 0.5)}61 xd = e.join(ind[["diversity"]], rsuffix="_i"); v3["partial_corr_eci_diversity_given_covariates"] = round(float(stats.spearmanr(y - X @ b, xd.diversity - X @ np.linalg.lstsq(X, xd.diversity.values, rcond=None)[0]).correlation), 3)62 out["H7_V3"] = v363 df.reset_index().rename(columns={"index": "region"}).to_csv(PROC / "h2_region_table.csv", index=False)64 write_json(out, EXP_DIR / "h2_scale.json", "h2"); ledger("models.h2_scale", **{k: v for k, v in out.items() if k != "correlations"}, primary=out["correlations"].get("H_rarefied15 ~ log_population"))65 log.info(json.dumps({k: v for k, v in out["correlations"].items() if "H_rarefied15" in k or "per_100k" in k})); log.info(json.dumps(v3))666768if __name__ == "__main__":69 main()70