1"""Harmonise five printed effect sizes (F4-conf-03). Reads only data/, writes out/results.json.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, math, os, re, csv45ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))6D = os.path.join(ROOT, "data")7SRC = os.path.join(D, "sources")8OUT = os.path.join(ROOT, "out")910def read_lines(path):11 with open(path, encoding="utf-8") as f:12 return f.read().split("\n")1314def wc_l(path):15 with open(path, "rb") as f:16 return f.read().count(b"\n")1718def find_row(lines, pattern, fname):19 """Return (line_number, cells) of the unique line matching regex pattern; print it."""20 hits = [(i + 1, l) for i, l in enumerate(lines) if re.match(pattern, l)]21 assert len(hits) == 1, f"pattern {pattern!r} matched {len(hits)} lines in {fname}"22 n, l = hits[0]23 print(f"[{fname}:{n}] {l}")24 return n, [c.strip() for c in l.split("|")]2526def num(cell):27 """Parse a printed number, stripping significance stars, %, 'pp' and thousands separators."""28 m = re.match(r"^\s*([-+]?\d[\d,]*\.?\d*)", cell.replace("*", ""))29 assert m, cell30 return float(m.group(1).replace(",", ""))3132effects = json.load(open(os.path.join(D, "effects.json"), encoding="utf-8"))33files = {34 "arxiv-2609.27732v1": "arxiv-2609.27732v1.txt",35 "arxiv-2606.14621v1": "arxiv-2606.14621v1.txt",36 "arxiv-2505.04555v3": "arxiv-2505.04555v3.txt",37 "arxiv-2510.03658v1": "arxiv-2510.03658v1.txt",38}39L = {k: read_lines(os.path.join(SRC, v)) for k, v in files.items()}40H = {}4142# ---- E1: DCR overall ATT, Table 6, 5 km -> relative % change ------------------------43s = "arxiv-2609.27732v1"44_, hdr = find_row(L[s], r"^\| 1 km \| 3 km \| 5 km$", files[s])45n_att, att = find_row(L[s], r"^ATT \| -0\.399", files[s])46n_pct, pct = find_row(L[s], r"^% change \| -32\.9%", files[s])47col = hdr.index("5 km") # header has a leading empty cell, same layout as data rows48b1 = num(att[col])49v1 = 100 * (math.exp(b1) - 1)50H["E1"] = dict(value=round(v1, 4), unit="percent (relative change vs. counterfactual hospitalisation rate)",51 transformation=f"100*(exp(b)-1) with b = {att[col].replace('*','')} (printed ATT, log points) -> 100*(exp({b1})-1) = {v1:.4f}%",52 harmonisable=True, missing_quantity=None,53 locator=f"arxiv-2609.27732v1 | Table 6 | row 'ATT' | column '5 km' (line {n_att}); cross-check: row '% change' | column '5 km' = {pct[col]} (line {n_pct})")54print(f"E1: b={b1} -> {v1:.4f}% (printed % change {pct[col]})")5556# ---- E4: DCR fatal, Table 7, 5 km -> absolute change per 100,000 --------------------57_, hdr7 = find_row(L[s], r"^Outcome \| 1 km \| 3 km \| 5 km$", files[s])58n_fat, fat = find_row(L[s], r"^Drug-related \(fatal\) \|", files[s])59b4 = num(fat[hdr7.index("5 km")])60rel4 = 100 * (math.exp(b4) - 1)61# Search the extract for any printed baseline/mean level of the fatal rate.62lvl = [l for l in L[s] if re.search(r"(mean|baseline|pre-treatment level|average rate)", l, re.I)63 and re.search(r"fatal", l, re.I) and re.search(r"\d+(\.\d+)? (per|hospitali[sz]ations per)", l, re.I)]64print(f"E4: b={b4}; relative change would be {rel4:.4f}%; lines printing a fatal-rate level: {len(lvl)}")65H["E4"] = dict(value=None, unit="hospitalisations per 100,000 inhabitants per year",66 transformation=(f"absolute change = R0*(exp(b)-1) with b = {b4} (printed ATT, log points); "67 f"exp(b)-1 = {rel4/100:.4f} is computable, but R0 (counterfactual/baseline fatal rate) is not printed"),68 harmonisable=False,69 missing_quantity=("baseline (counterfactual, e.g. pre-opening treated-region mean) level of the fatal drug-related "70 "hospitalisation rate per 100,000 inhabitants per year in the 5 km treated regions; not printed in the extract"),71 locator=f"arxiv-2609.27732v1 | Table 7 | row 'Drug-related (fatal)' | column '5 km' (line {n_fat})")7273# ---- E2: RESEA reemployment ATE (pp) -> % of control mean ---------------------------74s = "arxiv-2606.14621v1"75_, hdr2 = find_row(L[s], r"^Outcome \| Treatment \| Control \| Treatment-control difference", files[s])76n_re, re_row = find_row(L[s], r"^Reemployment \| \d+\.\d%", files[s])77d2 = num(re_row[hdr2.index("Average treatment effect estimate")].split(";")[0])78p2 = num(re_row[hdr2.index("Control")])79v2 = 100 * d2 / p280H["E2"] = dict(value=round(v2, 4), unit="percent (relative change vs. control-group mean reemployment rate)",81 transformation=f"100*Delta/p with Delta = {d2} pp (printed ATE) and p = {p2}% (printed control mean) -> 100*{d2}/{p2} = {v2:.4f}%",82 harmonisable=True, missing_quantity=None,83 locator=(f"arxiv-2606.14621v1 | Table 2 | row 'Reemployment' | column 'Average treatment effect estimate' (1.47 pp) "84 f"and column 'Control' (75.0%) (line {n_re})"))85print(f"E2: {d2}/{p2} -> {v2:.4f}%")8687# ---- E3: MW employment elasticity * printed % MW change ----------------------------88s = "arxiv-2505.04555v3"89n_el, el = find_row(L[s], r"^Employment elasticity w\.r\.t MW \|", files[s])90n_g, g = find_row(L[s], r"^% MW changes \|", files[s])91n_note, _ = find_row(L[s], r"^Note: The table reports the effects of a minimum wage increase", files[s])92eta, gg = num(el[1]), num(g[1])93v3 = 100 * eta * gg94H["E3"] = dict(value=round(v3, 4), unit="percent (relative change in employment)",95 transformation=f"100*eta*g with eta = {eta} (printed elasticity) and g = {gg} (printed '% MW changes', proportional) -> 100*{eta}*{gg} = {v3:.4f}%",96 harmonisable=True, missing_quantity=None,97 locator=(f"arxiv-2505.04555v3 | Table 2 | row 'Employment elasticity w.r.t MW' | column 'Estimate (SE)' (line {n_el}); "98 f"row '% MW changes' | column 'Estimate (SE)' (line {n_g})"))99print(f"E3: {eta}*{gg} -> {v3:.4f}%")100101# ---- E5: LPM coefficient -> percentage points --------------------------------------102s = "arxiv-2510.03658v1"103_, hdr5 = find_row(L[s], r"^\| \(1\) \| \(2\) \| \(3\) \| \(4\)$", files[s])104n_t, tr = find_row(L[s], r"^Treatment \| ", files[s])105n_m, mean = find_row(L[s], r"^Mean Dep\. Var\. \|", files[s])106n_pc, pc = find_row(L[s], r"^% Change \|", files[s])107c4 = hdr5.index("(4)")108b5 = num(tr[c4])109v5 = 100 * b5110H["E5"] = dict(value=round(v5, 4), unit="percentage points",111 transformation=f"100*beta with beta = {b5} (printed LPM coefficient, 0-1 scale) -> {v5:.4f} pp",112 harmonisable=True, missing_quantity=None,113 locator=f"arxiv-2510.03658v1 | Table 3 | row 'Treatment' | column '(4)' (line {n_t})")114print(f"E5: {b5} -> {v5:.4f} pp (check: beta/mean = {b5/num(mean[c4]):.5f}, printed % Change {pc[c4]})")115116# ---- consistency diagnostics used for data_issues ----------------------------------117s = "arxiv-2609.27732v1"118_, a19 = find_row(L[s], r"^ATT \| -0\.431\*\* \| -0\.438", files[s])119_, a19se = find_row(L[s], r"^\| \(0\.218\)", files[s])120_, t6se = find_row(L[s], r"^\| \(0\.412\)", files[s])121_, a19n = find_row(L[s], r"^Note: All specifications use the 5 km buffer", files[s])122abs_line = [l for l in L[s] if "which fall by 35%" in l]123mech_line = [l for l in L[s] if "by approximately 30% at 5 km" in l]124print(f"check E4 text: abstract '35%' present={bool(abs_line)}, Sec 6.1 '30%' present={bool(mech_line)}, exp(-0.354)-1={rel4:.2f}%")125126data_issues = [127 {"issue": "arxiv-2609.27732v1: the abstract says fatal drug-related hospitalisations 'fall by 35%%', but Table 7 (5 km) prints -0.354 log points, i.e. 100*(exp(-0.354)-1) = %.1f%%, matching Section 6.1 ('approximately 30%% at 5 km'); 35%% is the overall (Table 6) effect." % rel4,128 "handling": "Table 7 cell used as the printed estimate for E4; abstract figure not used. E4 is not harmonisable for a separate reason (missing baseline level)."},129 {"issue": "arxiv-2609.27732v1: the same CS-DR 5 km ATT (-0.431) is printed with SE (0.210) in Table 6 but (0.218) in Appendix Table A19; A19 note says 640 never-treated regions while its table and Table 6 say 608.",130 "handling": "Only the point estimate (identical in both tables) is used; SEs and region counts are not needed and were not used."},131 {"issue": "arxiv-2505.04555v3: Table 2 prints '% MW changes' = 0.047 (4.7%), while the table note describes the 2023 increases as 'about 4%'.",132 "handling": "Per the target definition, the table value 0.047 is used for E3; with 'about 4%%' the implied change would be 100*%s*0.04 = %.4f%% (sensitivity only, not the reported value)." % (eta, 100 * eta * 0.04)},133 {"issue": "arxiv-2505.04555v3: Table 2 prints the same SE (0.091) for 'Affected employment' and 'Emp. elasticity w.r.t affected wage' (possible duplication); point estimates are internally consistent (-0.268/0.366 = -0.732).",134 "handling": "Not used; no uncertainty is propagated."},135 {"issue": "arxiv-2510.03658v1: Table 3 prints '%% Change' = -0.3685 for column (4), whereas -0.0962/0.261 = %.5f (rounds to -0.3686); the text also attributes the 0.0855 estimate to 'table 2' and calls an 18.79%% relative change a 'percentage change' while interpreting it as lives saved per 100 deaths." % (b5 / num(mean[c4])),136 "handling": "Not material for E5, whose target (percentage points) needs only the printed coefficient -0.0962; noted for transparency."},137 {"issue": "arxiv-2606.14621v1: the ATE (1.47 pp) is block(week)-adjusted while the 75.0% control mean is the raw aggregate control-group mean; the abstract rounds the ATE to 1.5 pp.",138 "handling": "Used the printed control mean 75.0% as the counterfactual level p, as the target definition specifies (percent of the control-group mean); 1.47 (table) used rather than 1.5 (abstract)."},139]140141order = ["E1", "E2", "E3", "E4", "E5"]142assert [e["id"] for e in effects] == order143harm = {k: H[k] for k in order}144n_h = sum(1 for k in order if harm[k]["value"] is not None)145n_nh = sum(1 for k in order if harm[k]["harmonisable"] is False)146n_input_rows = {"effects.json": len(effects), "README.md": wc_l(os.path.join(D, "README.md"))}147for v in files.values():148 n_input_rows[v] = wc_l(os.path.join(SRC, v))149150res = {151 "task_id": "F4-conf-03",152 "harmonised": harm,153 "n_harmonised": n_h,154 "n_not_harmonisable": n_nh,155 "feasibility": "complete",156 "feasibility_reason": None,157 "data_issues": data_issues,158 "n_input_rows": n_input_rows,159 "n_analysed": len(effects),160 "random_seed": None,161 # P2 baseline: naive readings without the proper transformation (for comparison only)162 "baseline_naive": {163 "E1_logpoints_times_100_pct": round(100 * b1, 4),164 "E2_printed_pp": d2,165 "E3_printed_elasticity": eta,166 "E4_logpoints_times_100_pct": round(100 * b4, 4),167 "E5_printed_coefficient": b5,168 "note": "Naive values = printed numbers read at face value (log points x100 treated as percent; no rescaling). Not harmonised estimates.",169 },170}171os.makedirs(OUT, exist_ok=True)172with open(os.path.join(OUT, "results.json"), "w", encoding="utf-8") as f:173 json.dump(res, f, indent=2, ensure_ascii=False)174with open(os.path.join(OUT, "exclusions.csv"), "w", newline="", encoding="utf-8") as f:175 w = csv.writer(f)176 w.writerow(["step", "rule", "n_before", "n_after", "n_removed", "rationale"])177 w.writerow([1, "keep all items listed in data/effects.json", len(effects), len(effects), 0,178 "No item excluded; E4 is analysed and reported as not harmonisable rather than removed."])179print(f"n_harmonised={n_h} n_not_harmonisable={n_nh} feasibility=complete -> out/results.json written")180