1"""Phase 6 — QDACI components by region with store-bootstrap CIs, and composite instability across weightings.2Components: (1) rarefied diversity; (2) specialization = mean |log RAA| (scheme M); (3) complexity — omitted when H7 V1 fails;3(4a) concentration = within-region Gini of catalogue sizes; (4b) coverage = stores per 100 000 population; (5) price dispersion = H5 region4effect (regions with >= 10 matched stores only). Composites: equal weights on rank-normalised components, PCA first component, Borda.5Publication rule (preregistered): no composite if the median pairwise Spearman rho across composites < 0.7."""6from __future__ import annotations7import json8import numpy as np, pandas as pd9from scipy import stats10from src.common import ROOT, SEED, get_logger, ledger, register_output, write_json, PROC_DIR, EXP_DIR, P11from src.models.h1_raa import weights, table, raa12from src.models.h3_longtail import gini1314log = get_logger("models.qdaci")15PROC = PROC_DIR161718def spec_index(d, regions, cats):19 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)}).values20 t = table(rcode, ccode, weights(d, "M"), len(regions), len(cats)); A = raa(t); A[t == 0] = np.nan21 return np.nanmean(np.abs(np.log(A)), axis=1)222324def main():25 rng = np.random.default_rng(SEED)26 sc = pd.read_parquet(PROC / "store_category.parquet"); d = sc[sc.region_final != "Inconnue"].reset_index(drop=True)27 regions = sorted(d.region_final.unique()); cats = sorted(d.category.unique())28 div = pd.read_csv(PROC / "diversity_region.csv").set_index("region"); h7 = json.load(open(EXP_DIR / "h7_complexity.json"))29 ext = pd.read_csv(ROOT / "data" / "external" / "regions_external.csv").set_index("region_name"); sp = pd.read_parquet(PROC / "store_panel.parquet")30 spec = pd.Series(spec_index(d, regions, cats), index=regions)31 # bootstrap specialization32 stores = d.store_id.unique(); sidx = d.store_id.map({s: i for i, s in enumerate(stores)}).values; rows_by_store = [np.where(sidx == i)[0] for i in range(len(stores))]33 sb = np.array([spec_index(d.iloc[np.concatenate([rows_by_store[i] for i in rng.integers(0, len(stores), len(stores))])], regions, cats) for _ in range(500)])34 comp = pd.DataFrame(index=regions); comp["diversity_rarefied"] = div.H_rarefied15.reindex(regions); comp["diversity_lo"] = div.H_rarefied15_lo.reindex(regions); comp["diversity_hi"] = div.H_rarefied15_hi.reindex(regions)35 comp["specialization"] = spec; comp["specialization_lo"] = np.nanpercentile(sb, 2.5, axis=0); comp["specialization_hi"] = np.nanpercentile(sb, 97.5, axis=0)36 g = sp[sp.region_final != "Inconnue"].groupby("region_final"); comp["concentration_gini"] = g.n_listings.apply(lambda v: gini(v.values) if len(v) >= 5 else np.nan).reindex(regions)37 comp["coverage_stores_per_100k"] = (g.size() / ext.population.reindex(g.size().index) * 1e5).reindex(regions)38 eff = pd.read_csv(PROC / "h5_region_effects.csv"); eff = eff[(eff.label == "primary_AB") & (eff.n_stores >= 10)].set_index("region"); comp["price_effect"] = eff.effect.reindex(regions); comp["price_effect_se"] = eff.se.reindex(regions)39 include_eci = h7["V1_bootstrap_rank_stability"]["pass_eci"] and h7["V2_reshuffle_null"]["pass"]40 if include_eci: comp["complexity_eci"] = pd.read_csv(PROC / "h7_region_indicators.csv").set_index("region").eci.reindex(regions)41 core = [c for c in ("diversity_rarefied", "specialization", "concentration_gini", "coverage_stores_per_100k", "complexity_eci") if c in comp.columns]42 X = comp[core].dropna(); R = X.rank() / len(X)43 # orientation: higher = more diverse, more specialised, less concentrated (invert Gini), more covered44 R["concentration_gini"] = 1 - R["concentration_gini"]45 eq = R.mean(axis=1); borda = R.rank().sum(axis=1)46 Z = (X - X.mean()) / X.std(); Z["concentration_gini"] = -Z["concentration_gini"]; U, S, Vt = np.linalg.svd(Z.values - Z.values.mean(0), full_matrices=False); pc = pd.Series(Z.values @ Vt[0], index=Z.index)47 if np.corrcoef(pc, eq)[0, 1] < 0: pc = -pc48 rhos = [stats.spearmanr(a, b).correlation for a, b in ((eq, pc), (eq, borda), (pc, borda))]49 # bootstrap instability of the equal-weight ranking: resample stores, recompute diversity (merchant-equal, all stores) & specialization & gini50 rank_rho = []51 for b in range(200):52 idx = np.concatenate([rows_by_store[i] for i in rng.integers(0, len(stores), len(stores))]); dd = d.iloc[idx]53 spb = spec_index(dd, regions, cats); divb = {}54 for r_, gg in dd.groupby("region_final"):55 w = gg.n / gg.groupby("store_id").n.transform("sum"); p = (w.groupby(gg.category).sum() / w.sum()).values; p = p[p > 0]; divb[r_] = -(p * np.log(p)).sum() / np.log(20)56 Rb = pd.DataFrame({"d": pd.Series(divb).reindex(X.index), "s": pd.Series(spb, index=regions).reindex(X.index)}).rank() / len(X)57 rank_rho.append(stats.spearmanr(Rb.mean(axis=1), R[["diversity_rarefied", "specialization"]].mean(axis=1)).correlation)58 res = {"components": core, "price_dispersion_regions": int(comp.price_effect.notna().sum()), "complexity_included": bool(include_eci),59 "composite_agreement_spearman": {"equal_vs_pca": round(rhos[0], 3), "equal_vs_borda": round(rhos[1], 3), "pca_vs_borda": round(rhos[2], 3), "median": round(float(np.median(rhos)), 3)},60 "pca_first_component_share": round(float(S[0] ** 2 / (S ** 2).sum()), 3), "pca_loadings": dict(zip(core, np.round(Vt[0], 3).tolist())),61 "bootstrap_rank_stability_equal_weights_median_rho": round(float(np.nanmedian(rank_rho)), 3),62 "publish_composite": bool(np.median(rhos) >= 0.7), "note": "Components describe observed online catalogues of Quebec-based merchants, not production or output. Complexity omitted because the ECI failed the preregistered stability criterion." if not include_eci else "Components describe observed online catalogues, not production."}63 comp["composite_equal"] = eq.reindex(regions); comp["composite_pca"] = pc.reindex(regions); comp["composite_borda"] = borda.reindex(regions)64 comp.reset_index().rename(columns={"index": "region"}).to_csv(PROC / "qdaci_components.csv", index=False); register_output(PROC / "qdaci_components.csv", "qdaci")65 write_json(res, EXP_DIR / "qdaci.json", "qdaci"); ledger("models.qdaci", **res); log.info(json.dumps(res))666768if __name__ == "__main__":69 main()70