1"""H1 — Regional specialization: RAA matrix, Cramér's V, store-level region-reshuffle permutation test, store bootstrap CIs.2Preregistered: scheme M primary, decision requires p < 0.05 under M and T (Holm within the confirmatory family, applied in src/models/summary.py).3Outputs experiments/h1_raa.json, data/processed/raa_matrix.csv (point estimates + CIs per sample/scheme).4"""5from __future__ import annotations6import json7import numpy as np, pandas as pd8from src.common import ROOT, P, REGIONS, SEED, get_logger, ledger, register_output, write_json910log = get_logger("models.h1")11PROC = ROOT / "data" / "processed"12N_PERM, N_BOOT = 2000, 200013CATS = None141516def weights(sc: pd.DataFrame, scheme: str) -> np.ndarray:17 sizes = sc.groupby("store_id").n.transform("sum")18 if scheme == "L": return sc.n.values.astype(float)19 if scheme == "M": return (sc.n / sizes).values20 cut = sc.groupby("store_id").n.sum().quantile(0.99)21 return np.where(sizes > cut, 0.0, sc.n.values.astype(float))222324def table(region_codes, cat_codes, w, R, C):25 t = np.zeros((R, C)); np.add.at(t, (region_codes, cat_codes), w); return t262728def cramers_v(t):29 n = t.sum(); r = t.sum(1, keepdims=True); c = t.sum(0, keepdims=True)30 e = r @ c / n; m = e > 031 chi2 = ((t[m] - e[m]) ** 2 / e[m]).sum()32 return float(np.sqrt(chi2 / (n * (min(t.shape) - 1))))333435def raa(t):36 n = t.sum(); return (t / t.sum(1, keepdims=True)) / (t.sum(0, keepdims=True) / n)373839def run(sc, sample, scheme, rng):40 d = sc[(sc.region_final != "Inconnue")] if sample == "S_R" else sc[sc.region_db.notna()].assign(region_final=lambda x: x.region_dbfirst)41 d = d[d.region_final != "Inconnue"].reset_index(drop=True)42 regions = sorted(d.region_final.unique()); cats = sorted(d.category.unique())43 rcode = d.region_final.map({r: i for i, r in enumerate(regions)}).values; ccode = d.category.map({c: i for i, c in enumerate(cats)}).values44 w = weights(d, scheme); R, C = len(regions), len(cats)45 t = table(rcode, ccode, w, R, C); v = cramers_v(t); A = raa(t)46 # permutation: reshuffle region labels across STORES47 stores = d.store_id.unique(); sreg = d.drop_duplicates("store_id").set_index("store_id").region_final48 sidx = d.store_id.map({s: i for i, s in enumerate(stores)}).values49 srcode = np.array([ {r: i for i, r in enumerate(regions)}[sreg[s]] for s in stores])50 vperm = np.empty(N_PERM)51 for k in range(N_PERM):52 perm = rng.permutation(srcode); vperm[k] = cramers_v(table(perm[sidx], ccode, w, R, C))53 p = float((np.sum(vperm >= v) + 1) / (N_PERM + 1))54 # store bootstrap55 vb = np.empty(N_BOOT); Ab = np.empty((N_BOOT, R, C))56 rows_by_store = [np.where(sidx == i)[0] for i in range(len(stores))]57 for b in range(N_BOOT):58 draw = rng.integers(0, len(stores), len(stores))59 idx = np.concatenate([rows_by_store[i] for i in draw])60 tb = table(rcode[idx], ccode[idx], w[idx], R, C)61 vb[b] = cramers_v(tb); Ab[b] = raa(np.where(tb == 0, 1e-12, tb))62 lo, hi = np.nanpercentile(Ab, 2.5, axis=0), np.nanpercentile(Ab, 97.5, axis=0)63 mat = pd.DataFrame({"sample": sample, "scheme": scheme, "region": np.repeat(regions, C), "category": np.tile(cats, R),64 "n": t.ravel(), "raa": A.ravel(), "raa_lo": lo.ravel(), "raa_hi": hi.ravel()})65 mat["sig_gt1"] = mat.raa_lo > 1; mat["sig_lt1"] = mat.raa_hi < 166 res = {"sample": sample, "scheme": scheme, "stores": int(len(stores)), "regions": R, "categories": C, "cramers_v": round(v, 4),67 "v_boot_percentile_note": "store resampling with replacement duplicates stores and biases V upward; not used for inference",68 "v_boot_p2.5_97.5": [round(float(np.percentile(vb, 2.5)), 4), round(float(np.percentile(vb, 97.5)), 4)],69 "perm_null_mean": round(float(vperm.mean()), 4), "perm_null_p5": round(float(np.percentile(vperm, 5)), 4), "excess_v": round(float(v - vperm.mean()), 4),70 "perm_null_p95": round(float(np.percentile(vperm, 95)), 4), "p_perm": p, "cells_sig_gt1": int(mat.sig_gt1.sum()), "cells_sig_lt1": int(mat.sig_lt1.sum()), "cells": R * C}71 return res, mat727374def main():75 sc = pd.read_parquet(PROC / "store_category.parquet")76 rng = np.random.default_rng(SEED)77 results, mats = [], []78 for sample in ("S_R", "S_R0"):79 for scheme in ("M", "L", "T"):80 res, mat = run(sc, sample, scheme, rng); results.append(res); mats.append(mat); log.info(json.dumps(res))81 # merchant-count RAA (stores instead of listings) for comparison82 d = sc[sc.region_final != "Inconnue"].drop_duplicates(["store_id", "category"]).assign(n=1)83 res_m, mat_m = run(d, "S_R", "L", rng); res_m["scheme"] = "stores"; mat_m["scheme"] = "stores"84 results.append(res_m); mats.append(mat_m)85 pd.concat(mats).to_csv(PROC / "raa_matrix.csv", index=False); register_output(PROC / "raa_matrix.csv", "h1")86 write_json({"n_perm": N_PERM, "n_boot": N_BOOT, "results": results}, ROOT / "experiments" / "h1_raa.json", "h1"); ledger("models.h1_raa", results=results)878889if __name__ == "__main__":90 main()91