1"""F4-conf-03: harmonise five printed effect sizes. Reads data/ only; writes out/results.json.2No randomness. Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, math, re, os45ROOT = os.getcwd()6D = os.path.join(ROOT, "data")7SRC = os.path.join(D, "sources")8910def read_lines(path):11 with open(path, encoding="utf-8") as f:12 return f.read().split("\n")131415def find(lines, prefix, fname):16 """Return (line_no, line) of the unique line starting with prefix; print it."""17 hits = [(i + 1, l) for i, l in enumerate(lines) if l.startswith(prefix)]18 assert len(hits) == 1, (fname, prefix, hits)19 n, l = hits[0]20 print(f"[{fname}:{n}] {l}")21 return n, l222324def cells(line):25 return [c.strip() for c in line.split("|")]262728def num(s):29 m = re.search(r"-?\d[\d,]*\.?\d*", s.replace("−", "-"))30 return float(m.group(0).replace(",", ""))313233effects = json.load(open(os.path.join(D, "effects.json"), encoding="utf-8"))34files = {35 "dcr": "arxiv-2609.27732v1.txt",36 "resea": "arxiv-2606.14621v1.txt",37 "mw": "arxiv-2505.04555v3.txt",38 "cs": "arxiv-2510.03658v1.txt",39}40L = {k: read_lines(os.path.join(SRC, v)) for k, v in files.items()}4142# n_input_rows: wc -l semantics (newline count) for text, top-level items for JSON43n_input_rows = {"effects.json": len(effects)}44for name in ["README.md"] + ["sources/" + v for v in files.values()]:45 with open(os.path.join(D, name), "rb") as f:46 n_input_rows[name] = f.read().count(b"\n")4748H = {}49issues = []5051# ---------------- E1: log-point ATT -> relative percent change ----------------52print("\nE1 source lines:")53_, hdr = find(L["dcr"], "| 1 km | 3 km | 5 km", files["dcr"])54n_att, att = find(L["dcr"], "ATT | -0.399", files["dcr"])55_, pct = find(L["dcr"], "% change | -32.9%", files["dcr"])56b1 = num(cells(att)[3]) # 5 km column57assert cells(hdr)[3] == "5 km"58v1 = 100 * (math.exp(b1) - 1)59printed_pct_5km = num(cells(pct)[3])60H["E1"] = {61 "value": round(v1, 2),62 "unit": "percent change relative to the counterfactual hospitalisation rate",63 "transformation": f"100*(exp(b)-1) with b = {cells(att)[3]} (log points) -> 100*(exp({b1})-1) = {v1:.4f} %; consistent with printed '% change' {cells(pct)[3]}",64 "harmonisable": True,65 "missing_quantity": None,66 "locator": f"arxiv-2609.27732v1 | Table 6 | row 'ATT' | column '5 km' (line {n_att}); cross-check row '% change', column '5 km'",67}68assert abs(round(v1, 1) - printed_pct_5km) < 1e-96970# ---------------- E2: pp effect -> percent of control mean ----------------71print("\nE2 source lines:")72_, hdr2 = find(L["resea"], "Outcome | Treatment | Control |", files["resea"])73n_re, re_row = find(L["resea"], "Reemployment | 75.2%", files["resea"])74c = cells(re_row)75assert cells(hdr2)[2] == "Control" and cells(hdr2)[4] == "Average treatment effect estimate"76p_ctrl = num(c[2]) # 75.0 (%)77delta = num(c[4].split(";")[0]) # 1.47 pp78v2 = 100 * delta / p_ctrl79H["E2"] = {80 "value": round(v2, 2),81 "unit": "percent of the control-group mean reemployment rate",82 "transformation": f"100*Delta/p with Delta = {c[4].split(';')[0]} and p = {c[2]} (control mean) -> 100*{delta}/{p_ctrl} = {v2:.4f} %",83 "harmonisable": True,84 "missing_quantity": None,85 "locator": f"arxiv-2606.14621v1 | Table 2 | row 'Reemployment' | columns 'Average treatment effect estimate' (Delta) and 'Control' (p) (line {n_re})",86}8788# ---------------- E3: elasticity x printed MW change ----------------89print("\nE3 source lines:")90n_el, el = find(L["mw"], "Employment elasticity w.r.t MW |", files["mw"])91n_g, gline = find(L["mw"], "% MW changes |", files["mw"])92_, note = find(L["mw"], "Note: The table reports the effects of a minimum wage increase", files["mw"])93eta = num(cells(el)[1])94g = num(cells(gline)[1])95v3 = 100 * eta * g96H["E3"] = {97 "value": round(v3, 2),98 "unit": "percent change in employment",99 "transformation": f"100*eta*g with eta = {cells(el)[1]} and g = {cells(gline)[1]} (proportional MW change, read as 4.7%; the table note says 'about 4%') -> 100*{eta}*{g} = {v3:.4f} %",100 "harmonisable": True,101 "missing_quantity": None,102 "locator": f"arxiv-2505.04555v3 | Table 2 | row 'Employment elasticity w.r.t MW' | column 'Estimate (SE)' (line {n_el}); row '% MW changes' | column 'Estimate (SE)' (line {n_g})",103}104issues.append({105 "issue": "arxiv-2505.04555v3 Table 2: the row '% MW changes' prints 0.047, a label suggesting percent units while the value is only consistent with a proportion (4.7%); the table note describes the change as 'about 4%'.",106 "handling": "Read 0.047 as a proportional change g = 0.047 (4.7%), the only reading compatible with the note ('about 4%'); the printed table value (not the rounded note value) is used."})107108# ---------------- E4: log-point ATT on fatal hosp. -> absolute change per 100k ----------------109print("\nE4 source lines:")110_, hdr4 = find(L["dcr"], "Outcome | 1 km | 3 km | 5 km", files["dcr"])111n_f, fat = find(L["dcr"], "Drug-related (fatal) |", files["dcr"])112b4 = num(cells(fat)[3])113rel4 = 100 * (math.exp(b4) - 1)114# search the extract for any printed level of the outcome (mean/baseline rate); none expected115level_hits = [l for l in L["dcr"] if re.search(r"(mean|baseline rate|pre-treatment (mean|level)|average rate)", l, re.I)116 and re.search(r"fatal", l, re.I)]117print(f"Lines printing a level of the fatal-hospitalisation rate: {len(level_hits)}")118H["E4"] = {119 "value": None,120 "unit": "hospitalisations per 100,000 inhabitants per year",121 "transformation": f"Absolute change = R0*(exp(b)-1) with b = {cells(fat)[3]} (log points), i.e. relative change 100*(exp({b4})-1) = {rel4:.4f} %; R0 (counterfactual/baseline fatal rate per 100,000 in treated regions) is not printed in the extract, so the absolute change cannot be computed",122 "harmonisable": False,123 "missing_quantity": "baseline (counterfactual, e.g. pre-opening or control-group mean) level of the fatal drug-related hospitalisation rate per 100,000 inhabitants per year in the 5 km treated regions",124 "locator": f"arxiv-2609.27732v1 | Table 7 | row 'Drug-related (fatal)' | column '5 km' (line {n_f})",125}126assert not level_hits127128# ---------------- E5: LPM coefficient -> percentage points ----------------129print("\nE5 source lines:")130_, hdr5 = find(L["cs"], "| (1) | (2) | (3) | (4)", files["cs"])131n_t, tr = find(L["cs"], "Treatment |", files["cs"])132_, mean5 = find(L["cs"], "Mean Dep. Var. |", files["cs"])133_, pc5 = find(L["cs"], "% Change |", files["cs"])134b5 = num(cells(tr)[4])135v5 = 100 * b5136H["E5"] = {137 "value": round(v5, 2),138 "unit": "percentage points",139 "transformation": f"100*beta with beta = {cells(tr)[4]} (change in probability, 0-1 scale) -> 100*({b5}) = {v5:.4f} pp",140 "harmonisable": True,141 "missing_quantity": None,142 "locator": f"arxiv-2510.03658v1 | Table 3 | row 'Treatment' | column '(4)' (line {n_t})",143}144145# ---------------- consistency checks on the sources (recorded as data issues) ----------------146print("\nSource-consistency checks:")147_, abs_dcr = find(L["dcr"], "This paper estimates the causal effect of drug consumption room", files["dcr"])148issues.append({149 "issue": f"arxiv-2609.27732v1: the abstract states that fatal outcomes 'fall by 35%', but Table 7 (5 km) prints {cells(fat)[3]} log points, i.e. 100*(exp({b4})-1) = {rel4:.1f}%; Section 6.1 says 'approximately 30%', consistent with Table 7. The 35% matches the overall ATT (Table 6), not the fatal outcome.",150 "handling": "Table 7 is used as the authoritative printed estimate for E4; the abstract figure is not used."})151n_a19, a19 = find(L["dcr"], "ATT | -0.431** | -0.438**", files["dcr"])152a19_se = L["dcr"][n_a19 + 1]153t6_se = L["dcr"][n_att + 1]154print(f"[{files['dcr']}:{n_att+2}] {t6_se}")155print(f"[{files['dcr']}:{n_a19+2}] {a19_se}")156issues.append({157 "issue": f"arxiv-2609.27732v1: the 5 km CS-DR ATT -0.431 has SE {cells(t6_se)[3]} in Table 6 but {cells(a19_se)[1]} in Table A19; the Table A19 note mentions 640 never-treated comparison regions while the table body and Tables 6-7 print 608.",158 "handling": "Does not affect the point estimate used in E1 (identical -0.431 in both tables); Table 6 is used as the cited location. No uncertainty is propagated."})159_, a19pct = find(L["dcr"], "% change | -35.0% | -35.5%", files["dcr"])160mism = []161for b, p in zip(cells(a19)[1:], cells(a19pct)[1:]):162 bb, pp = num(b), num(p)163 implied = 100 * (math.exp(bb) - 1)164 if abs(round(implied, 1) - pp) > 1e-9:165 mism.append(f"{b} -> {implied:.2f}% vs printed {p}")166print("A19 % change mismatches:", mism)167if mism:168 issues.append({"issue": "arxiv-2609.27732v1 Table A19: printed '% change' differs in the last decimal from 100*(exp(ATT)-1) for: " + "; ".join(mism) + ".",169 "handling": "Minor rounding inconsistency in a robustness table not used for any item; recorded only."})170ratio = b5 / num(cells(mean5)[4])171print(f"Tossou col (4): beta/mean = {ratio:.4f} vs printed % Change {cells(pc5)[4]}")172issues.append({173 "issue": f"arxiv-2510.03658v1: Table 3 column (4) prints '% Change' {cells(pc5)[4]} while {b5}/{num(cells(mean5)[4])} = {ratio:.4f} (last-digit difference, presumably from an unrounded mean). The abstract's '0.0855 (a 18.79 percentage change)' refers to a 'table 2' (mean 0.455) that is not in the extract.",174 "handling": "E5 uses the Table 3 column (4) coefficient as printed and converts only the scale (x100 -> pp); the abstract/Table 2 numbers and the '% Change' row are not used."})175issues.append({176 "issue": "Rendering defects: arxiv-2510.03658v1 footnote markers are rendered inline (e.g. 'percent^9^9 9') and footnote 11 is interleaved with main text; arxiv-2606.14621v1 renders 'p < 0.001' as 'p ¡ 0.001'.",177 "handling": "No numbers used in E1-E5 are affected; recorded only."})178issues.append({179 "issue": "arxiv-2606.14621v1 Table 2: the control-group mean (75.0%) is a raw aggregate mean, whereas the ATE (1.47 pp) adjusts for randomisation week (blocks); the abstract rounds the ATE to 1.5 pp.",180 "handling": "Following the target definition ('percent of the control-group mean'), the printed control mean 75.0% is used as p and the Table 2 ATE 1.47 pp (not the rounded 1.5) as Delta; the denominator mismatch is noted as a limitation."})181182n_h = sum(1 for v in H.values() if v["value"] is not None)183n_nh = sum(1 for v in H.values() if v["harmonisable"] is False)184results = {185 "task_id": "F4-conf-03",186 "harmonised": H,187 "n_harmonised": n_h,188 "n_not_harmonisable": n_nh,189 "feasibility": "complete",190 "feasibility_reason": None,191 "data_issues": issues,192 "n_input_rows": n_input_rows,193 "n_analysed": len(effects),194 "random_seed": None,195 # P2 baseline: the printed estimates on their reported definitions, untransformed196 "baseline_printed_values": {"E1": b1, "E2": delta, "E3": eta, "E4": b4, "E5": b5},197 "baseline_printed_units": {"E1": "log points", "E2": "percentage points", "E3": "elasticity (unitless)",198 "E4": "log points", "E5": "probability change (0-1 scale)"},199}200assert [e["id"] for e in effects] == list(H.keys())201assert n_h + n_nh == len(effects)202with open(os.path.join(ROOT, "out", "results.json"), "w", encoding="utf-8") as f:203 json.dump(results, f, indent=2, ensure_ascii=False)204print("\nSummary:")205for k, v in H.items():206 print(f" {k}: value={v['value']} [{v['unit']}] harmonisable={v['harmonisable']}")207print(f"n_harmonised={n_h} n_not_harmonisable={n_nh} feasibility=complete")208