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_json, PROC_DIR, EXP_DIR910log = get_logger("models.h1")11PROC = PROC_DIR12N_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 if sample == "S_R": d = sc[sc.region_final != "Inconnue"]41 elif sample == "S_R0": d = sc[sc.region_db.notna()].assign(region_final=lambda x: x.region_db)42 else: d = sc.assign(region_final=lambda x: x.region_pipfirst)43 d = d[d.region_final != "Inconnue"].reset_index(drop=True)44 regions = sorted(d.region_final.unique()); cats = sorted(d.category.unique())45 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)}).values46 w = weights(d, scheme); R, C = len(regions), len(cats)47 t = table(rcode, ccode, w, R, C); v = cramers_v(t); A = raa(t)48 # permutation: reshuffle region labels across STORES49 stores = d.store_id.unique(); sreg = d.drop_duplicates("store_id").set_index("store_id").region_final50 sidx = d.store_id.map({s: i for i, s in enumerate(stores)}).values51 srcode = np.array([ {r: i for i, r in enumerate(regions)}[sreg[s]] for s in stores])52 vperm = np.empty(N_PERM)53 for k in range(N_PERM):54 perm = rng.permutation(srcode); vperm[k] = cramers_v(table(perm[sidx], ccode, w, R, C))55 p = float((np.sum(vperm >= v) + 1) / (N_PERM + 1))56 # store bootstrap57 vb = np.empty(N_BOOT); Ab = np.empty((N_BOOT, R, C))58 rows_by_store = [np.where(sidx == i)[0] for i in range(len(stores))]59 for b in range(N_BOOT):60 draw = rng.integers(0, len(stores), len(stores))61 idx = np.concatenate([rows_by_store[i] for i in draw])62 tb = table(rcode[idx], ccode[idx], w[idx], R, C)63 vb[b] = cramers_v(tb); Ab[b] = raa(np.where(tb == 0, 1e-12, tb))64 lo, hi = np.nanpercentile(Ab, 2.5, axis=0), np.nanpercentile(Ab, 97.5, axis=0)65 # calibration of the cell-level intervals: expected number of cells with CI > 1 / < 1 under the merchant-reshuffle null66 # (20 permutations x 300 store bootstraps; no multiplicity control is applied to the 340 cells, so the raw counts67 # must be read against these null counts — referee 1)68 null_gt, null_lt = [], []69 for k in range(20):70 perm = rng.permutation(srcode); rp = perm[sidx]; Ab_ = np.empty((300, R, C))71 for b in range(300):72 draw = rng.integers(0, len(stores), len(stores)); idx = np.concatenate([rows_by_store[i] for i in draw])73 tb_ = table(rp[idx], ccode[idx], w[idx], R, C); Ab_[b] = raa(np.where(tb_ == 0, 1e-12, tb_))74 lo_, hi_ = np.nanpercentile(Ab_, 2.5, axis=0), np.nanpercentile(Ab_, 97.5, axis=0)75 null_gt.append(int((lo_ > 1).sum())); null_lt.append(int((hi_ < 1).sum()))76 mat = pd.DataFrame({"sample": sample, "scheme": scheme, "region": np.repeat(regions, C), "category": np.tile(cats, R),77 "n": t.ravel(), "raa": A.ravel(), "raa_lo": lo.ravel(), "raa_hi": hi.ravel()})78 mat["sig_gt1"] = mat.raa_lo > 1; mat["sig_lt1"] = mat.raa_hi < 179 res = {"sample": sample, "scheme": scheme, "stores": int(len(stores)), "regions": R, "categories": C, "cramers_v": round(v, 4),80 "v_boot_percentile_note": "store resampling with replacement duplicates stores and biases V upward; not used for inference",81 "v_boot_p2.5_97.5": [round(float(np.percentile(vb, 2.5)), 4), round(float(np.percentile(vb, 97.5)), 4)],82 "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),83 "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,84 "null_cells_sig_gt1_mean_max": [round(float(np.mean(null_gt)), 1), int(max(null_gt))], "null_cells_sig_lt1_mean_max": [round(float(np.mean(null_lt)), 1), int(max(null_lt))]}85 return res, mat868788def main():89 sc = pd.read_parquet(PROC / "store_category.parquet")90 rng = np.random.default_rng(SEED)91 results, mats = [], []92 for sample in ("S_R", "S_R0", "S_Rpip"):93 for scheme in ("M", "L", "T"):94 res, mat = run(sc, sample, scheme, rng); results.append(res); mats.append(mat); log.info(json.dumps(res))95 # merchant-count RAA (stores instead of listings) for comparison96 d = sc[sc.region_final != "Inconnue"].drop_duplicates(["store_id", "category"]).assign(n=1)97 res_m, mat_m = run(d, "S_R", "L", rng); res_m["scheme"] = "stores"; mat_m["scheme"] = "stores"98 results.append(res_m); mats.append(mat_m)99 pd.concat(mats).to_csv(PROC / "raa_matrix.csv", index=False); register_output(PROC / "raa_matrix.csv", "h1")100 write_json({"n_perm": N_PERM, "n_boot": N_BOOT, "results": results}, EXP_DIR / "h1_raa.json", "h1"); ledger("models.h1_raa", results=results)101102103if __name__ == "__main__":104 main()105