1"""Robustness variants (preregistered §5): H1 Cramér's V with store-reshuffle permutation p, H2-type diversity range and H3 Gini2recomputed on (i) S_R0, (ii) certification/association sources only, (iii) excluding Montréal, (iv) excluding the 1 % largest stores,3(v) 8 macro-categories, (vi) 'autre' dropped, (vii) categories reassigned by the audit confusion matrix (100 draws, if the audit exists),4(viii) excluding resellers/collectives, (ix) region conflicts dropped. Output experiments/robustness.json"""5from __future__ import annotations6import json7import numpy as np, pandas as pd8from src.common import ROOT, SEED, get_logger, ledger, write_json, PROC_DIR, EXP_DIR, P, BASE_PROC, BASE_EXP9from src.models.h1_raa import weights, table, cramers_v10from src.models.h3_longtail import gini1112log = get_logger("validation.robustness")13PROC = PROC_DIR14MACRO = {"erable": "food", "epicerie": "food", "cafe_the": "food", "chocolat": "food", "boulangerie": "food", "viande_poisson": "food", "fromage_laitier": "food", "boissons": "food",15 "alcool": "alcohol", "sante_beaute": "beauty", "maison": "home", "mode": "fashion", "bijoux": "fashion", "art": "art", "enfants": "kids_pets", "animaux": "kids_pets",16 "plein_air": "outdoor_garden", "jardin": "outdoor_garden", "industriel": "industrial", "autre": "other"}17CERT_SOURCES = {"alimentsduquebec_ent", "lesproduitsduquebec_ent", "lesproduitsduquebec_det", "cartv_bio", "ctaq", "erabledici", "fromagesdici", "cidreduquebec", "vinsduquebec",18 "distilleriesduquebec", "ambq_profils", "artisansaloeuvre", "afmq", "apiculteursduquebec", "boeufdici", "zoneboreale", "createursdesaveurs", "croquezoutaouais",19 "arretsgourmands", "goutezy", "laurentidesjenmange", "goutcotenord", "gaspesiegourmande", "charlevoix", "mauriciemiam", "gardemangerduquebec", "saveursbsl", "goutezat"}202122def h1(d, rng, scheme="M", nperm=1000):23 d = d.reset_index(drop=True); regions = sorted(d.region_final.unique()); cats = sorted(d.category.unique())24 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)}).values25 w = weights(d, scheme); R, C = len(regions), len(cats); v = cramers_v(table(rcode, ccode, w, R, C))26 stores = d.store_id.unique(); sidx = d.store_id.map({s: i for i, s in enumerate(stores)}).values27 sreg = d.drop_duplicates("store_id").set_index("store_id").region_final.reindex(stores).map({r: i for i, r in enumerate(regions)}).values28 null = np.array([cramers_v(table(rng.permutation(sreg)[sidx], ccode, w, R, C)) for _ in range(nperm)])29 return {"V": round(v, 4), "null_mean": round(float(null.mean()), 4), "p": float((np.sum(null >= v) + 1) / (nperm + 1)), "stores": int(len(stores)), "regions": R}303132def ent_range(d):33 out = {}34 for r, g in d.groupby("region_final"):35 w = g.n / g.groupby("store_id").n.transform("sum"); p = (w.groupby(g.category).sum() / w.sum()).values; p = p[p > 0]36 out[r] = float(-(p * np.log(p)).sum() / np.log(d.category.nunique()))37 s = pd.Series(out); return {"min": round(float(s.min()), 3), "max": round(float(s.max()), 3), "top": s.idxmax(), "bottom": s.idxmin()}383940def main():41 rng = np.random.default_rng(SEED)42 sc = pd.read_parquet(PROC / "store_category.parquet"); sp = pd.read_parquet(PROC / "store_panel.parquet").set_index("id")43 base = sc[sc.region_final != "Inconnue"]44 variants = {"baseline_S_R": base, "S_R0_dbfirst_regions": sc[sc.region_db.notna()].assign(region_final=lambda x: x.region_dbfirst).query("region_final != 'Inconnue'"),45 "cert_assoc_sources_only": base[base.first_source.isin(CERT_SOURCES)], "excl_Montreal": base[base.region_final != "Montréal"],46 "excl_top1pct_stores": base[~base.store_id.isin(sc.groupby("store_id").n.sum().pipe(lambda s: s[s > s.quantile(0.99)]).index)],47 "macro8_categories": base.assign(category=base.category.map(MACRO)).groupby(["store_id", "category", "region_final", "region_source", "first_source", "store_kind"], as_index=False).n.sum(),48 "drop_autre": base[base.category != "autre"], "excl_resellers_collectives": base[~base.store_kind.isin(["revendeur", "collectif"])],49 "excl_region_conflicts": base[~base.store_id.map(sp.region_conflict).fillna(False)]}50 res = {}51 for k, d in variants.items():52 res[k] = {"h1_M": h1(d, rng, "M"), "h1_T": h1(d, rng, "T"), "entropy_M": ent_range(d), "gini_listings": round(gini(d.groupby("store_id").n.sum().values), 4)}53 log.info("%s %s", k, json.dumps(res[k]))54 audit = BASE_EXP / "category_audit.json"55 if audit.exists():56 conf = json.load(open(audit)).get("confusion_rule_vs_consensus", {})57 if conf:58 cats = sorted(sc.category.unique()); Vs = []59 P_ = {r: {c: v for c, v in row.items()} for r, row in conf.items()}60 for k in range(100):61 d = base.copy()62 def draw(c):63 row = P_.get(c); 64 if not row: return c65 ks, vs = list(row.keys()), np.array(list(row.values()), float); return rng.choice(ks, p=vs / vs.sum())66 d["category"] = [draw(c) for c in d.category]67 d = d.groupby(["store_id", "category", "region_final", "region_source", "first_source", "store_kind"], as_index=False).n.sum()68 Vs.append(h1(d, rng, "M", nperm=200))69 res["confusion_reassignment_100draws"] = {"V_mean": round(float(np.mean([v["V"] for v in Vs])), 4), "V_min": round(float(min(v["V"] for v in Vs)), 4), "share_p_lt_05": round(float(np.mean([v["p"] < 0.05 for v in Vs])), 3)}70 write_json(res, EXP_DIR / "robustness.json", "robustness"); ledger("validation.robustness", variants=list(res))717273if __name__ == "__main__":74 main()75