"""F4-conf-01: extract printed treatment effects / uncertainty from two arXiv extracts. Reads only data/; writes out/results.json deterministically. No randomness. Run from workspace root: ~/wp29/venv/bin/python out/analysis.py """ import json, math, os, re, hashlib DATA = "data" DCR_ID, RI_ID = "arxiv-2609.27732v1", "arxiv-2606.14621v1" DCR_F = os.path.join(DATA, "sources", DCR_ID + ".txt") RI_F = os.path.join(DATA, "sources", RI_ID + ".txt") Z = 1.96 # normal critical value for two-sided 95% def lines(path): with open(path, encoding="utf-8") as f: return f.read().split("\n") def cells(line): return [c.strip() for c in line.split("|")] def num(s): """Parse a printed number: strip stars, $, %, thousands commas, parentheses.""" s = s.replace("$", "").replace(",", "").replace("*", "").strip().strip("()") return int(s) if re.fullmatch(r"-?\d+", s) else float(s) def stars(s): m = re.search(r"(\*+)\s*$", s.strip()) return len(m.group(1)) if m else 0 def table_block(L, title_prefix, end_prefix="Note:"): i = next(k for k, l in enumerate(L) if l.startswith(title_prefix)) j = next(k for k in range(i + 1, len(L)) if L[k].startswith(end_prefix)) return i, [(k, L[k]) for k in range(i + 1, j) if L[k].strip()] def reg_table(L, title_prefix): """Parse a regression table: header, coefficient rows followed by SE rows.""" i, blk = table_block(L, title_prefix) hdr = cells(blk[0][1]) cols = [c for c in hdr if c and c != "Outcome"] rows, last = {}, None for k, l in blk[1:]: c = cells(l) if c[0] == "" and last is not None and c[1].startswith("("): rows[last]["se"] = dict(zip(cols, c[1:])) rows[last]["se_line"] = (k + 1, l) else: last = c[0] rows[last] = {"est": dict(zip(cols, c[1:])), "line": (k + 1, l)} return cols, rows def show(tag, ln): print(f"[{tag}] line {ln[0]}: {ln[1]}") def main(): DL, RL = lines(DCR_F), lines(RI_F) res, loc, issues = {}, {}, [] def L(src, t, row, col, extra=""): return f'{src} | Table {t} | row "{row}" | column "{col}"' + extra # ---------------- DCR Table 6 ---------------- c6, t6 = reg_table(DL, "Table 6:") show("T6 ATT", t6["ATT"]["line"]); show("T6 SE", t6["ATT"]["se_line"]) show("T6 Obs", t6["Observations"]["line"]) res["dcr_att_5km"] = num(t6["ATT"]["est"]["5 km"]) res["dcr_att_5km_se"] = num(t6["ATT"]["se"]["5 km"]) res["dcr_att_3km"] = num(t6["ATT"]["est"]["3 km"]) res["dcr_att_3km_se"] = num(t6["ATT"]["se"]["3 km"]) res["dcr_obs_5km"] = num(t6["Observations"]["est"]["5 km"]) loc["dcr_att_5km"] = L(DCR_ID, 6, "ATT", "5 km") loc["dcr_att_5km_se"] = L(DCR_ID, 6, "ATT", "5 km", " (standard error in parentheses, line below)") loc["dcr_att_3km"] = L(DCR_ID, 6, "ATT", "3 km") loc["dcr_att_3km_se"] = L(DCR_ID, 6, "ATT", "3 km", " (standard error in parentheses, line below)") loc["dcr_obs_5km"] = L(DCR_ID, 6, "Observations", "5 km") # ---------------- DCR Table 7 ---------------- c7, t7 = reg_table(DL, "Table 7:") for r in ("Narcotics", "Drug-related (fatal)", "HIV"): show("T7 " + r, t7[r]["line"]); show("T7 SE " + r, t7[r]["se_line"]) spec = [("dcr_narcotics_att_3km", "Narcotics", "3 km"), ("dcr_fatal_att_1km", "Drug-related (fatal)", "1 km"), ("dcr_hiv_att_5km", "HIV", "5 km")] for key, r, c in spec: res[key] = num(t7[r]["est"][c]); loc[key] = L(DCR_ID, 7, r, c) res[key + "_se"] = num(t7[r]["se"][c]) loc[key + "_se"] = L(DCR_ID, 7, r, c, " (standard error in parentheses, line below)") hiv_cell = t7["HIV"]["est"]["5 km"] # exact p-value: search extract for any printed p-value pv_hits = [l for l in DL if re.search(r"\bp\s*=\s*0?\.\d", l)] print(f"[DCR] lines with an exact 'p=' value in extract: {len(pv_hits)}") res["dcr_hiv_att_5km_p_value"] = None loc["dcr_hiv_att_5km_p_value"] = ("not reported in " + DCR_ID + " (Table 7 prints only the estimate, " "its SE in parentheses and significance stars; no exact p-value is printed)") res["dcr_hiv_att_5km_significant_10pct"] = stars(hiv_cell) >= 1 # ---------------- derived CI (DCR) ---------------- b, s = res["dcr_att_5km"], res["dcr_att_5km_se"] res["dcr_att_5km_ci_low"] = round(b - Z * s, 4) res["dcr_att_5km_ci_high"] = round(b + Z * s, 4) for k in ("dcr_att_5km_ci_low", "dcr_att_5km_ci_high"): loc[k] = (f'{DCR_ID} | Table 6 | row "ATT" | column "5 km" | derived from printed ATT ' f'and SE: ATT {"-" if k.endswith("low") else "+"} 1.96 x SE') # ---------------- RESEA Table 2 ---------------- i2 = next(k for k, l in enumerate(RL) if l.startswith("Outcome | Treatment | Control")) hdr2 = cells(RL[i2]) t2 = {} for k in range(i2 + 1, len(RL)): if not RL[k].strip(): break c = cells(RL[k]); t2[c[0]] = (dict(zip(hdr2[1:], c[1:])), (k + 1, RL[k])) for r in ("Wage", "Reemployment", "Weeks on UI"): show("T2 " + r, t2[r][1]) ATE = "Average treatment effect estimate" def ate_parse(cell): m = re.match(r"\s*(-?\$?-?[\d,\.]+)\s*(pp|weeks|wks)?;\s*p\s*([=<])\s*([\d\.]+);\s*95% CI:\s*\[([^\]]+)\]", cell) est = num(m.group(1)); lo, hi = [num(x) for x in m.group(5).split(", ")] return est, m.group(3) + m.group(4), lo, hi w_est, w_p, w_lo, w_hi = ate_parse(t2["Wage"][0][ATE]) r_est, r_p, r_lo, r_hi = ate_parse(t2["Reemployment"][0][ATE]) u_est, u_p, u_lo, u_hi = ate_parse(t2["Weeks on UI"][0][ATE]) res["ri_wage_ate"], res["ri_wage_ate_ci_low"], res["ri_wage_ate_ci_high"] = w_est, w_lo, w_hi res["ri_wage_ate_se_reported"] = bool(re.search(r"\bSE\b|standard error|\(\d", t2["Wage"][0][ATE])) res["ri_wage_ate_se_derived"] = round((w_hi - w_lo) / (2 * Z), 2) ctrl = re.match(r"\$([\d,]+)", t2["Wage"][0]["Control"]).group(1) res["ri_control_mean_wage"] = num(ctrl) res["ri_reemployment_ate_pp"] = r_est res["ri_weeks_ui_ate"], res["ri_weeks_ui_ate_ci_low"], res["ri_weeks_ui_ate_ci_high"] = u_est, u_lo, u_hi loc["ri_wage_ate"] = L(RI_ID, 2, "Wage", ATE, " (point estimate)") loc["ri_wage_ate_ci_low"] = L(RI_ID, 2, "Wage", ATE, " (95% CI, lower bound)") loc["ri_wage_ate_ci_high"] = L(RI_ID, 2, "Wage", ATE, " (95% CI, upper bound)") loc["ri_wage_ate_se_derived"] = L(RI_ID, 2, "Wage", ATE, " | derived from printed 95% CI: (upper - lower) / (2 x 1.96)") loc["ri_control_mean_wage"] = L(RI_ID, 2, "Wage", "Control", " (mean)") loc["ri_reemployment_ate_pp"] = L(RI_ID, 2, "Reemployment", ATE, " (point estimate)") loc["ri_weeks_ui_ate"] = L(RI_ID, 2, "Weeks on UI", ATE, " (point estimate)") loc["ri_weeks_ui_ate_ci_low"] = L(RI_ID, 2, "Weeks on UI", ATE, " (95% CI, lower bound)") loc["ri_weeks_ui_ate_ci_high"] = L(RI_ID, 2, "Weeks on UI", ATE, " (95% CI, upper bound)") # ---------------- consistency checks -> data_issues ---------------- # (a) Table 6 vs Table A19 SE for the same 5 km CS-DR ATT cA, tA = reg_table(DL, "Table A19:") a19_att, a19_se = num(tA["ATT"]["est"]["CS-DR"]), num(tA["ATT"]["se"]["CS-DR"]) show("A19 ATT", tA["ATT"]["line"]); show("A19 SE", tA["ATT"]["se_line"]) if a19_att == res["dcr_att_5km"] and a19_se != res["dcr_att_5km_se"]: issues.append({"issue": f"Same 5 km CS-DR ATT ({a19_att}) printed with SE {t6['ATT']['se']['5 km']} in Table 6 " f"but SE {tA['ATT']['se']['CS-DR']} in Table A19 (CS-DR column) of {DCR_ID}.", "handling": "Used Table 6 (the table named in the key definitions) for dcr_att_5km_se and the derived CI; " "the significance star (**) is the same under either SE."}) # (b) control-region count: 608 (Tables 6/7, Fig. 6) vs 640 (A19 note); observations arithmetic a19_note = next(l for l in DL if "never-treated comparison regions (" in l) m = re.search(r"\((\d+) regions\)", a19_note) n_ctrl6 = num(t6["Control regions"]["est"]["5 km"]) obs = {c: num(t6["Observations"]["est"][c]) for c in c6} trt = {c: num(t6["Treated regions"]["est"][c]) for c in c6} per_reg = {c: obs[c] / (trt[c] + n_ctrl6) for c in c6} per_reg640 = {c: obs[c] / (trt[c] + 640) for c in c6} print(f"[check] obs/(treated+608) = {per_reg}; obs/(treated+640) = {per_reg640}") if m and int(m.group(1)) != n_ctrl6: issues.append({"issue": f"Control-region count inconsistent within {DCR_ID}: Tables 6/7 and the Figure 6 note say " f"{n_ctrl6} never-treated regions; the Table A19 note says {m.group(1)}. Table 6 observations " f"divided by (treated + {n_ctrl6}) are non-integers ({', '.join(f'{v:.2f}' for v in per_reg.values())}), " f"whereas divided by (treated + 640) they equal {', '.join(f'{v:g}' for v in per_reg640.values())}.", "handling": "Reported dcr_obs_5km exactly as printed in Table 6 (17,250); flagged, not resolved (the panel structure is not stated in the extract)."}) # (c) abstract: fatal outcomes 'fall by 35%' vs Table 7 5 km implied change fat5 = num(t7["Drug-related (fatal)"]["est"]["5 km"]); narc5 = num(t7["Narcotics"]["est"]["5 km"]) pct = lambda x: 100 * (math.exp(x) - 1) print(f"[check] exp(ATT)-1: fatal 5km {pct(fat5):.1f}%, narcotics 5km {pct(narc5):.1f}%, overall 5km {pct(res['dcr_att_5km']):.1f}%") if any("fall by 35%" in l for l in DL): issues.append({"issue": f"Abstract of {DCR_ID} says fatal outcomes 'fall by 35%', but Table 7 fatal 5 km ATT {fat5} implies " f"{pct(fat5):.1f}% (Section 6.1 says ~30% at 5 km); 35% matches the overall 5 km ATT. The text also calls HIV " f"estimates 'close to zero' although the 5 km HIV ATT is {t7['HIV']['est']['5 km']} (SE {t7['HIV']['se']['5 km']}, no star).", "handling": "Narrative statements not used; all values taken from the tables named in the keys."}) # (d) RESEA text gives weeks-on-UI CI as positive magnitudes if any("[1.75 weeks, 2.24 weeks]" in l for l in RL): issues.append({"issue": f"{RI_ID} text reports the weeks-on-UI CI as magnitudes '[1.75 weeks, 2.24 weeks]' (a reduction), " f"while Table 2 prints [{u_lo}, {u_hi}].", "handling": "Used Table 2 signs and order (lower bound first)."}) # (e) rendering defect if any("p ¡ 0.001" in l for l in RL): issues.append({"issue": f"Rendering defect in {RI_ID} text: 'p ¡ 0.001' (should read p < 0.001).", "handling": "Not used for any key; Table 2 values used."}) # star / SE consistency (not an issue if consistent) from statistics import NormalDist nd = NormalDist() star_ok = [] for tab, rows, cols in (("6", t6, c6), ("7", t7, c7)): for r, d in rows.items(): if "se" not in d: continue for c in cols: e, se_ = num(d["est"][c]), num(d["se"][c]) p = 2 * (1 - nd.cdf(abs(e / se_))) exp_st = 3 if p < .01 else 2 if p < .05 else 1 if p < .10 else 0 star_ok.append((tab, r, c, round(p, 4), stars(d["est"][c]), exp_st)) mism = [x for x in star_ok if x[4] != x[5]] print(f"[check] star vs normal-approx p from est/SE: {len(star_ok)} cells, {len(mism)} mismatches {mism}") hiv_p_approx = [x[3] for x in star_ok if x[1] == "HIV" and x[2] == "5 km"][0] # RESEA CI symmetry print(f"[check] wage CI half-widths: {w_est - w_lo}, {w_hi - w_est}; weeks: {u_est - u_lo:.2f}, {u_hi - u_est:.2f}") # ---------------- baseline / sensitivity ---------------- res["baseline_dcr_att_5km_ci_rule2"] = [round(b - 2 * s, 4), round(b + 2 * s, 4)] res["baseline_ri_wage_ate_se_rule2"] = round((w_hi - w_lo) / 4, 2) res["check_dcr_hiv_att_5km_normal_approx_p_not_printed"] = hiv_p_approx out = {} order = ["dcr_att_5km", "dcr_att_5km_se", "dcr_att_3km", "dcr_att_3km_se", "dcr_obs_5km", "dcr_narcotics_att_3km", "dcr_narcotics_att_3km_se", "dcr_fatal_att_1km", "dcr_fatal_att_1km_se", "dcr_hiv_att_5km", "dcr_hiv_att_5km_se", "dcr_hiv_att_5km_p_value", "dcr_hiv_att_5km_significant_10pct", "dcr_att_5km_ci_low", "dcr_att_5km_ci_high", "ri_wage_ate", "ri_wage_ate_ci_low", "ri_wage_ate_ci_high", "ri_wage_ate_se_reported", "ri_wage_ate_se_derived", "ri_control_mean_wage", "ri_reemployment_ate_pp", "ri_weeks_ui_ate", "ri_weeks_ui_ate_ci_low", "ri_weeks_ui_ate_ci_high"] for k in order: out[k] = res[k] out["locators"] = {k: loc[k] for k in order if k in loc} out["derivation_notes"] = { "dcr_att_5km_ci": "Normal-approximation (Wald) 95% CI from Table 6 printed values: ATT +/- 1.96 x SE = " f"{t6['ATT']['est']['5 km'].rstrip('*')} +/- 1.96 x {t6['ATT']['se']['5 km'].strip('()')}; rounded to 4 decimals. Assumes approx. normal sampling distribution of the " "CS-DR ATT with the printed clustered SE; a 2.0 rule gives " f"{res['baseline_dcr_att_5km_ci_rule2']}.", "ri_wage_ate_se": "SE implied by the Table 2 printed 95% CI assuming a symmetric normal (Wald) interval: " f"(upper - lower) / (2 x 1.96) = ({w_hi} - ({w_lo})) / 3.92; rounded to 2 decimals. With N = 23,549 a " f"t critical value is practically 1.96; a 2.0 rule gives {res['baseline_ri_wage_ate_se_rule2']}. " "Printed bounds are rounded to whole dollars, so the derived SE is accurate to about +/-0.2."} out["task_id"] = "F4-conf-01" out["feasibility"] = "complete" out["feasibility_reason"] = None out["data_issues"] = issues nrows = {} for p in (os.path.join(DATA, "README.md"), DCR_F, RI_F): with open(p, "rb") as f: nrows[os.path.relpath(p, DATA)] = f.read().count(b"\n") out["n_input_rows"] = nrows out["n_analysed"] = len(order) out["random_seed"] = None out["baseline_dcr_att_5km_ci_rule2"] = res["baseline_dcr_att_5km_ci_rule2"] out["baseline_ri_wage_ate_se_rule2"] = res["baseline_ri_wage_ate_se_rule2"] out["check_dcr_hiv_att_5km_normal_approx_p_not_printed"] = res["check_dcr_hiv_att_5km_normal_approx_p_not_printed"] out["input_sha256"] = {os.path.relpath(p, DATA): hashlib.sha256(open(p, "rb").read()).hexdigest() for p in (os.path.join(DATA, "README.md"), DCR_F, RI_F)} os.makedirs("out", exist_ok=True) with open(os.path.join("out", "results.json"), "w") as f: json.dump(out, f, indent=2, ensure_ascii=False); f.write("\n") print("Wrote out/results.json") if __name__ == "__main__": main()