#!/usr/bin/env python3 """M1-LONG selection diagnostics — RESEARCH ONLY. Uses published v2 authorized pairs + micro extract. Compares usable-pair origins vs eligible origins that fail retention on observed covariates. """ from __future__ import annotations import hashlib import json import math from datetime import datetime, timezone from pathlib import Path import numpy as np import pandas as pd ROOT = Path(__file__).resolve().parents[2] MICRO = ROOT / "data/raw/michigan_sca/microdata/sca_micro_2015_2023_pilot.csv" PAIRS = ROOT / "data/m1/v2/m1long_pairs_origin_le2019.csv" OUT = ROOT / "data/m1/v2" DOCS = ROOT / "docs/m1" SAMPLE_FRESH = 1 # verify against builder ORIGIN_MAX = 201912 LEFT_CENSOR = 201810 VALID_ATT = {1, 3, 5} def sha256_file(p: Path) -> str: h = hashlib.sha256() with p.open("rb") as f: for chunk in iter(lambda: f.read(1 << 20), b""): h.update(chunk) return h.hexdigest() def std_diff(m1, m0, s1, s0): pooled = math.sqrt(0.5 * (s1 ** 2 + s0 ** 2)) if pooled == 0 or np.isnan(pooled): return None return float((m1 - m0) / pooled) def main(): # confirm SAMPLE_FRESH from protocol proto = json.loads((ROOT / "data/m1/v2/m1long_v2_origin_le2019_protocol.json").read_text()) # fall back: check builder constants via import import sys sys.path.insert(0, str(Path(__file__).resolve().parent)) from build_m1_long_pilot import SAMPLE_FRESH as SF, VALID_ATT as VA sample_fresh = SF valid_att = set(VA) micro = pd.read_csv(MICRO) pairs = pd.read_csv(PAIRS) pairs["origin_yyyymm"] = pairs["origin_yyyymm"].astype(int) pairs["origin_id"] = pairs["origin_id"].astype(int) usable_keys = set(zip(pairs["origin_id"], pairs["origin_yyyymm"])) m = micro.copy() m["YYYYMM"] = m["YYYYMM"].astype(int) m["ID"] = m["ID"].astype(int) elig = m[ (m["SAMPLE"].isin(sample_fresh if isinstance(sample_fresh, (set, list, tuple, frozenset)) else [sample_fresh])) & (m["PEXP"].isin(valid_att)) & (m["YYYYMM"] >= LEFT_CENSOR) & (m["YYYYMM"] <= ORIGIN_MAX) ].copy() elig["retained"] = [((int(i), int(y)) in usable_keys) for i, y in zip(elig["ID"], elig["YYYYMM"])] for col in ["AGE", "SEX", "EDUC", "REGION", "PAGO", "PEXP", "YYYYMM", "METHOD"]: elig[col] = pd.to_numeric(elig[col], errors="coerce") elig["origin_year"] = elig["YYYYMM"] // 100 elig["origin_month"] = elig["YYYYMM"] % 100 retained = elig[elig["retained"]] dropped = elig[~elig["retained"]] def summarize_cont(col): def stats(df): s = df[col].dropna() return { "n": int(s.shape[0]), "missing": int(df[col].isna().sum()), "mean": float(s.mean()) if len(s) else None, "std": float(s.std(ddof=1)) if len(s) > 1 else None, } r, d = stats(retained), stats(dropped) sd = None if all(x is not None for x in (r["mean"], d["mean"], r["std"], d["std"])): sd = std_diff(r["mean"], d["mean"], r["std"], d["std"]) return {"retained": r, "dropped": d, "std_diff_retained_minus_dropped": sd} def summarize_cat(col): def props(df): s = df[col].dropna() vc = s.value_counts(normalize=True).sort_index() counts = s.value_counts().sort_index() out_p, out_c = {}, {} for k, v in vc.items(): key = str(int(k)) if float(k).is_integer() else str(k) out_p[key] = float(v) for k, v in counts.items(): key = str(int(k)) if float(k).is_integer() else str(k) out_c[key] = int(v) return {"n": int(s.shape[0]), "missing": int(df[col].isna().sum()), "proportions": out_p, "counts": out_c} r, d = props(retained), props(dropped) keys = sorted(set(r["proportions"]) | set(d["proportions"]), key=lambda x: (len(x), x)) diffs = {} for k in keys: pr = r["proportions"].get(k, 0.0) pd_ = d["proportions"].get(k, 0.0) pooled = math.sqrt(0.5 * (pr * (1 - pr) + pd_ * (1 - pd_))) sd = ((pr - pd_) / pooled) if pooled > 0 else None diffs[k] = {"prop_retained": pr, "prop_dropped": pd_, "diff": pr - pd_, "std_diff": sd} return {"retained": r, "dropped": d, "level_diffs": diffs} cont = {c: summarize_cont(c) for c in ["AGE", "YYYYMM"]} cats = {c: summarize_cat(c) for c in ["SEX", "EDUC", "REGION", "PAGO", "PEXP", "origin_year", "origin_month", "METHOD"]} # logistic diagnostic feature_cols = [] X_parts = [] age = elig["AGE"].astype(float) age_z = (age - age.mean()) / (age.std(ddof=0) or 1) X_parts.append(age_z.fillna(0).to_numpy().reshape(-1, 1)) feature_cols.append("age_z") for col, prefix in [("SEX", "sex"), ("EDUC", "educ"), ("REGION", "region"), ("PAGO", "pago"), ("PEXP", "pexp")]: dummies = pd.get_dummies(elig[col].astype("Int64"), prefix=prefix, dummy_na=True) if dummies.shape[1] > 1: dummies = dummies.iloc[:, 1:] X_parts.append(dummies.to_numpy(dtype=float)) feature_cols.extend(list(dummies.columns)) y_z = (elig["origin_year"] - elig["origin_year"].mean()) / (elig["origin_year"].std(ddof=0) or 1) X_parts.append(y_z.to_numpy().reshape(-1, 1)) feature_cols.append("origin_year_z") X = np.hstack(X_parts) y = elig["retained"].astype(int).to_numpy() Xd = np.hstack([np.ones((X.shape[0], 1)), X]) feature_cols_i = ["intercept"] + feature_cols beta = np.zeros(Xd.shape[1]) converged = False for _ in range(50): eta = Xd @ beta p = 1 / (1 + np.exp(-np.clip(eta, -30, 30))) W = np.clip(p * (1 - p), 1e-6, None) z = eta + (y - p) / W try: XtW = Xd.T * W beta_new = np.linalg.solve(XtW @ Xd, XtW @ z) except np.linalg.LinAlgError: beta_new = beta + 0.05 * (Xd.T @ (y - p)) if np.max(np.abs(beta_new - beta)) < 1e-6: beta = beta_new converged = True break beta = beta_new eta = Xd @ beta p = 1 / (1 + np.exp(-np.clip(eta, -30, 30))) ll = float(np.sum(y * np.log(np.clip(p, 1e-12, 1)) + (1 - y) * np.log(np.clip(1 - p, 1e-12, 1)))) p0 = float(y.mean()) ll0 = float(np.sum(y * np.log(p0) + (1 - y) * np.log(1 - p0))) mcfadden = 1 - (ll / ll0) if ll0 != 0 else None coefs = {feature_cols_i[i]: float(beta[i]) for i in range(len(beta))} # 3x3 on pairs u = pairs.copy() # column names may be pexp/pago lower pexp_col = "pexp" if "pexp" in u.columns else "PEXP" pago_col = "pago" if "pago" in u.columns else "PAGO" u[pexp_col] = pd.to_numeric(u[pexp_col], errors="coerce") u[pago_col] = pd.to_numeric(u[pago_col], errors="coerce") ct = pd.crosstab(u[pexp_col], u[pago_col]) crosstab = { str(int(r)): {str(int(c)): int(ct.loc[r, c]) for c in ct.columns} for r in ct.index } pack = { "title": "M1-LONG selection diagnostics (origin≤2019)", "generated_at": datetime.now(timezone.utc).isoformat(), "research_only": True, "score_authorized_M1": False, "L1_untouched": True, "response_to": [ "54396241-9636-4773-9f13-ee4ff8a06ee8", "astra:7b24493b-7368-451d-9d66-71d9211659fb", ], "definition": { "eligible_origins": f"SAMPLE={sample_fresh} (fresh), valid PEXP, YYYYMM in [{LEFT_CENSOR},{ORIGIN_MAX}]", "usable_pairs": "published v2 m1long_pairs_origin_le2019.csv (gap 11–13)", "N_eligible": int(len(elig)), "N_retained_usable": int(len(retained)), "N_dropped": int(len(dropped)), "retention_rate": float(len(retained) / len(elig)) if len(elig) else None, "N_usable_pairs_file": int(len(pairs)), "sample_fresh_codes": list(sample_fresh) if isinstance(sample_fresh,(set,list,tuple,frozenset)) else [sample_fresh], }, "unavailable_covariates": { "income": "NOT in pilot extract", "employment": "NOT in pilot extract", "party_id": "NOT in pilot extract", "note": "Re-extract from SDA required; no proxies invented.", }, "continuous": cont, "categorical": cats, "retention_logit_diagnostic": { "converged": converged, "n": int(len(y)), "base_rate": p0, "loglik": ll, "loglik_null": ll0, "mcfadden_pseudo_r2": mcfadden, "coefficients": coefs, "caveat": "Observed-covariate retention model cannot establish ignorable attrition for unobserved factors. Diagnostic only — not production weights.", }, "unweighted_3x3_pexp_by_pago": crosstab, "michigan_longitudinal_weighting": { "status": "public_docs_and_extract_check", "findings": [ "Extract has WT and WT_HH cross-section weights only.", "No dedicated longitudinal/reinterview attrition weight column in this extract.", "IPW from this diagnostic = proposed sensitivity only; needs independent defensibility review.", ], "codebook": "https://sda.umsurvey.org/sca/Doc/sca0001.htm", "sda": "https://sda.umsurvey.org/", }, "replication_artifacts": { "raw_extract": "data/raw/michigan_sca/microdata/sca_micro_2015_2023_pilot.csv", "raw_sha256": "8f74769178d2918818906ac41a99d795131af4592b17beefa4b1d2ef05771dcf", "pairs_v2": "data/m1/v2/m1long_pairs_origin_le2019.csv", "source_urls": [ "https://sda.umsurvey.org/", "https://sda.umsurvey.org/sda-public/cgi-bin/hsda?setupfile=harcsda&datasetname=sca&ui=1", "https://sda.umsurvey.org/sca/Doc/sca0001.htm", ], "scripts": [ "srp/m1/build_m1_long_pilot.py", "srp/m1/build_m1_long_v2.py", "srp/m1/build_m1_long_selection_diag.py", ], "note": "Do not republish raw microdata on public Hub. Hashes + scripts enable private-box replication.", }, "pre2020_note": "Origins ≤2019 include destinations in 2020 (e.g. 2019→2020). Pre-2020 origins ≠ pre-2020 outcomes.", } OUT.mkdir(parents=True, exist_ok=True) json_path = OUT / "m1long_v2_selection_diagnostics.json" json_path.write_text(json.dumps(pack, indent=2) + "\n") rows = [] for col, block in cont.items(): rows.append({ "covariate": col, "type": "continuous", "retained_n": block["retained"]["n"], "dropped_n": block["dropped"]["n"], "retained_mean": block["retained"]["mean"], "dropped_mean": block["dropped"]["mean"], "std_diff": block["std_diff_retained_minus_dropped"], "retained_missing": block["retained"]["missing"], "dropped_missing": block["dropped"]["missing"], }) csv_path = OUT / "m1long_v2_selection_balance.csv" pd.DataFrame(rows).to_csv(csv_path, index=False) age_sd = cont["AGE"]["std_diff_retained_minus_dropped"] memo = f"""# M1-LONG selection diagnostics (origin ≤ 2019) **Status:** RESEARCH ONLY · `score_authorized(M1)=false` · L1 untouched **Response to:** Astra/ChatGPT `54396241` / `7b24493b` (CONDITIONAL_ACCEPT_FEASIBILITY) ## Cohort | Group | N | |------|---| | Eligible origins (fresh, valid PEXP, 201810–201912) | {len(elig)} | | Retained as usable 11–13m pairs | {len(retained)} | | Dropped (eligible but not usable pair) | {len(dropped)} | | Retention rate | {len(retained)/len(elig):.3f} | | Usable pairs file (v2) | {len(pairs)} | ## Unavailable in this extract Income, employment, and party ID are **not** in `sca_micro_2015_2023_pilot.csv`. Re-extract from SDA required. No proxies invented. ## Continuous balance (retained − dropped) | Covariate | Retained mean | Dropped mean | Std diff | Retained N | Dropped N | |-----------|---------------|--------------|----------|------------|-----------| | AGE | {cont['AGE']['retained']['mean']:.3f} | {cont['AGE']['dropped']['mean']:.3f} | {age_sd:.3f} | {cont['AGE']['retained']['n']} | {cont['AGE']['dropped']['n']} | | YYYYMM | {cont['YYYYMM']['retained']['mean']:.1f} | {cont['YYYYMM']['dropped']['mean']:.1f} | {cont['YYYYMM']['std_diff_retained_minus_dropped']:.3f} | {cont['YYYYMM']['retained']['n']} | {cont['YYYYMM']['dropped']['n']} | ## Categorical (see JSON for full tables) Key std diffs (retained − dropped prop, bernoulli pooled): SEX/EDUC/REGION/PAGO/PEXP levels in `m1long_v2_selection_diagnostics.json`. ## Retention logit (diagnostic only) - Converged: `{converged}` - Base retention rate: `{p0:.3f}` - McFadden pseudo-R²: `{mcfadden:.4f}` - **Caveat (Astra/Codex):** observed-covariate models cannot establish ignorable attrition for unobserved factors. Not production weights. ## Michigan longitudinal weights Cross-section `WT`/`WT_HH` present; no dedicated longitudinal/reinterview attrition weight in this extract. IPW from this model = sensitivity proposal only. ## Pre-2020 origins vs outcomes Origins ≤2019 can have 2020 destinations. Pre-2020 origins ≠ pre-2020 outcomes. ## Replication pack (private) - Raw extract sha256: `8f74769178d2918818906ac41a99d795131af4592b17beefa4b1d2ef05771dcf` - Scripts: `build_m1_long_pilot.py`, `build_m1_long_v2.py`, `build_m1_long_selection_diag.py` - Artifacts: `data/m1/v2/m1long_v2_selection_diagnostics.json`, `m1long_v2_selection_balance.csv` - Do **not** publish raw microdata on public Hub. ## Unweighted 3×3 (PEXP × PAGO) See JSON. No IPW sensitivity applied this tick. """ DOCS.mkdir(parents=True, exist_ok=True) memo_path = DOCS / "WAVE2B_M1_LONG_SELECTION_DIAGNOSTICS.md" memo_path.write_text(memo) sums = [] for p in sorted(OUT.iterdir()): if p.is_file() and p.name != "SHA256SUMS.txt": sums.append(f"{sha256_file(p)} {p.name}") (OUT / "SHA256SUMS.txt").write_text("\n".join(sums) + "\n") print(json.dumps({ "N_eligible": int(len(elig)), "N_retained": int(len(retained)), "N_dropped": int(len(dropped)), "retention": float(len(retained) / len(elig)), "pairs_file": int(len(pairs)), "mcfadden": mcfadden, "age_std_diff": age_sd, "json_sha": sha256_file(json_path), "memo_sha": sha256_file(memo_path), "balance_sha": sha256_file(csv_path), "converged": converged, }, indent=2)) if __name__ == "__main__": main()