"""Replication script for F4-conf-01. Parses the two source extracts in data/sources/, prints the matched source lines, and writes out/results.json. Deterministic; no randomness.""" import json, os, re, math from scipy.stats import norm ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) DCR_ID, RI_ID = "arxiv-2609.27732v1", "arxiv-2606.14621v1" P_DCR = os.path.join(ROOT, "data", "sources", DCR_ID + ".txt") P_RI = os.path.join(ROOT, "data", "sources", RI_ID + ".txt") P_README = os.path.join(ROOT, "data", "README.md") def read(p): with open(p, encoding="utf-8") as f: return f.read() def wc_l(p): return read(p).count("\n") dcr = read(P_DCR).splitlines() ri = read(P_RI).splitlines() def block(lines, start_pat, end_pat): i0 = next(i for i, l in enumerate(lines) if l.startswith(start_pat)) i1 = next(i for i in range(i0 + 1, len(lines)) if lines[i].startswith(end_pat)) return lines[i0:i1] def cells(line): return [c.strip() for c in line.split("|")] def parse_reg_table(blk, header_line_idx_pred): """Return header cols and dict row -> (coef cells, se cells or None).""" nonempty = [l for l in blk if l.strip()] hdr = next(l for l in nonempty if header_line_idx_pred(l)) cols = cells(hdr)[1:] rows = {} for k, l in enumerate(nonempty): c = cells(l) if l is hdr or "|" not in l or c[0] == "": continue se = None if k + 1 < len(nonempty) and cells(nonempty[k + 1])[0] == "" and nonempty[k + 1].strip().startswith("| ("): se = cells(nonempty[k + 1])[1:] print(" matched:", l, "\n ", nonempty[k + 1]) else: print(" matched:", l) rows[c[0]] = (c[1:], se) return cols, rows def num(s): m = re.match(r"^\(?(-?[\d,]*\.?\d+)\)?(\**)", s.replace("$", "").strip()) return m.group(1).replace(",", ""), m.group(2) def as_printed(s): v, _ = num(s) return int(v) if re.fullmatch(r"-?\d+", v) else float(v) print("== DCR Table 6 ==") t6 = block(dcr, "Table 6:", "Note:") c6, r6 = parse_reg_table(t6, lambda l: l.strip().startswith("| 1 km")) print("== DCR Table 7 ==") t7 = block(dcr, "Table 7:", "Note:") c7, r7 = parse_reg_table(t7, lambda l: l.startswith("Outcome |")) print("== DCR Table A19 ==") ta = block(dcr, "Table A19:", "Note:") ca, ra = parse_reg_table(ta, lambda l: l.strip().startswith("| CS-DR")) def g(cols, rows, row, col, se=False): j = cols.index(col) c, s = rows[row] return (s if se else c)[j] def loc(src, tab, row, col, extra=""): return f'{src} | Table {tab} | row "{row}" | column "{col}"' + extra R, L, issues = {}, {}, [] def put(key, cols, rows, tab, row, col, se=False): raw = g(cols, rows, row, col, se) R[key] = as_printed(raw) L[key] = loc(DCR_ID, tab, row, col, " (standard error, line below coefficient)" if se else "") put("dcr_att_5km", c6, r6, 6, "ATT", "5 km") put("dcr_att_5km_se", c6, r6, 6, "ATT", "5 km", se=True) put("dcr_att_3km", c6, r6, 6, "ATT", "3 km") put("dcr_att_3km_se", c6, r6, 6, "ATT", "3 km", se=True) put("dcr_obs_5km", c6, r6, 6, "Observations", "5 km") put("dcr_narcotics_att_3km", c7, r7, 7, "Narcotics", "3 km") put("dcr_narcotics_att_3km_se", c7, r7, 7, "Narcotics", "3 km", se=True) put("dcr_fatal_att_1km", c7, r7, 7, "Drug-related (fatal)", "1 km") put("dcr_fatal_att_1km_se", c7, r7, 7, "Drug-related (fatal)", "1 km", se=True) put("dcr_hiv_att_5km", c7, r7, 7, "HIV", "5 km") put("dcr_hiv_att_5km_se", c7, r7, 7, "HIV", "5 km", se=True) # HIV p-value: search the whole extract for any printed p-value near HIV hiv_p = [l for l in dcr if "HIV" in l and re.search(r"p\s*[=<]\s*0?\.\d", l)] print("Lines mentioning HIV with an exact p-value:", hiv_p) R["dcr_hiv_att_5km_p_value"] = None L["dcr_hiv_att_5km_p_value"] = (f"not reported in {DCR_ID} (Table 7 prints only the coefficient, SE in parentheses " "and significance stars; no exact p-value is printed anywhere in the extract)") hiv_stars = num(g(c7, r7, "HIV", "5 km"))[1] R["dcr_hiv_att_5km_significant_10pct"] = len(hiv_stars) >= 1 # derived 95% CI for Table 6 5km ATT Z = 1.96 b, se = R["dcr_att_5km"], R["dcr_att_5km_se"] R["dcr_att_5km_ci_low"] = round(b - Z * se, 4) R["dcr_att_5km_ci_high"] = round(b + Z * se, 4) for k in ("dcr_att_5km_ci_low", "dcr_att_5km_ci_high"): L[k] = f'{DCR_ID} | Table 6 | row "ATT" | column "5 km" | derived from printed ATT (-0.431) and SE (0.210)' # ---------- RI ---------- print("== RI Table 2 ==") t2 = block(ri, "Outcome | Treatment | Control", "Table 2:") hdr2 = cells(t2[0]) rows2 = {} for l in t2[1:]: if l.strip(): print(" matched:", l) c = cells(l); rows2[c[0]] = c[1:] cols2 = hdr2[1:] ATE = "Average treatment effect estimate" def parse_effect(s): m = re.match(r"^(-?)\$?(-?[\d,]*\.?\d+)\s*(pp|weeks|wks)?;\s*p\s*([=<])\s*([\d.]+);\s*95% CI:\s*\[([^,\]]+(?:,\d{3})*),\s*([^\]]+)\]$", s) # split CI robustly: bounds may contain thousands separators ("2,141") est = m.group(1) + m.group(2).replace(",", "") ci = re.search(r"\[(.*)\]", s).group(1) parts = re.split(r",\s+", ci) lo, hi = [p.replace(",", "") for p in parts] return est, m.group(4) + m.group(5), lo, hi def tonum(v): return int(v) if re.fullmatch(r"-?\d+", v) else float(v) w_est, w_p, w_lo, w_hi = parse_effect(rows2["Wage"][cols2.index(ATE)]) R["ri_wage_ate"] = tonum(w_est) R["ri_wage_ate_ci_low"], R["ri_wage_ate_ci_high"] = tonum(w_lo), tonum(w_hi) for k, what in (("ri_wage_ate", "point estimate"), ("ri_wage_ate_ci_low", "95% CI lower bound"), ("ri_wage_ate_ci_high", "95% CI upper bound")): L[k] = f'{RI_ID} | Table 2 | row "Wage" | column "{ATE}" ({what})' se_printed = any(re.search(r"\bSE\b|standard error", l, re.I) for l in t2) R["ri_wage_ate_se_reported"] = bool(se_printed) R["ri_wage_ate_se_derived"] = round((R["ri_wage_ate_ci_high"] - R["ri_wage_ate_ci_low"]) / (2 * Z), 2) L["ri_wage_ate_se_derived"] = (f'{RI_ID} | Table 2 | row "Wage" | column "{ATE}" | derived from printed 95% CI ' f'[166, 2,141]; SE not reported in {RI_ID}') cm = rows2["Wage"][cols2.index("Control")] R["ri_control_mean_wage"] = tonum(re.match(r"\$([\d,]+)", cm).group(1).replace(",", "")) L["ri_control_mean_wage"] = f'{RI_ID} | Table 2 | row "Wage" | column "Control" (mean; SD printed separately)' r_est, r_p, r_lo, r_hi = parse_effect(rows2["Reemployment"][cols2.index(ATE)]) R["ri_reemployment_ate_pp"] = tonum(r_est) L["ri_reemployment_ate_pp"] = f'{RI_ID} | Table 2 | row "Reemployment" | column "{ATE}" (point estimate)' u_est, u_p, u_lo, u_hi = parse_effect(rows2["Weeks on UI"][cols2.index(ATE)]) R["ri_weeks_ui_ate"] = tonum(u_est) R["ri_weeks_ui_ate_ci_low"], R["ri_weeks_ui_ate_ci_high"] = tonum(u_lo), tonum(u_hi) for k, what in (("ri_weeks_ui_ate", "point estimate"), ("ri_weeks_ui_ate_ci_low", "95% CI lower bound"), ("ri_weeks_ui_ate_ci_high", "95% CI upper bound")): L[k] = f'{RI_ID} | Table 2 | row "Weeks on UI" | column "{ATE}" ({what})' # ---------- consistency checks ---------- def stars_expected(b, s): p = 2 * (1 - norm.cdf(abs(b / s))) return ("***" if p < .01 else "**" if p < .05 else "*" if p < .10 else ""), p star_mismatch = [] for tab, cols, rows in (("6", c6, r6), ("7", c7, r7), ("A19", ca, ra)): for row, (c, s) in rows.items(): if s is None: continue for j, col in enumerate(cols): if j >= len(c) or not c[j]: continue bv, st = num(c[j]); sv, _ = num(s[j]) exp, p = stars_expected(float(bv), float(sv)) if exp != st: star_mismatch.append(f"Table {tab} {row}/{col}: printed '{st}' vs normal-approx '{exp}' (p={p:.3f})") print("Star/SE mismatches:", star_mismatch) a19_se = as_printed(g(ca, ra, "ATT", "CS-DR", se=True)) if a19_se != R["dcr_att_5km_se"]: issues.append({"issue": f"Internal inconsistency in {DCR_ID}: the 5 km CS-DR overall ATT (-0.431**) has SE " f"{R['dcr_att_5km_se']} in Table 6 but SE {a19_se} in Appendix Table A19 (same estimator, same 5 km buffer).", "handling": "Used Table 6 as the key definition requires (SE 0.210); derived CI uses 0.210. " "With 0.218 the CI would be [" + f"{b - Z * a19_se:.4f}, {b + Z * a19_se:.4f}" + "]."}) n_ctrl = [l for l in dcr if "640 regions" in l] obs = {col: as_printed(g(c6, r6, "Observations", col)) for col in c6} tr = {col: as_printed(g(c6, r6, "Treated regions", col)) for col in c6} implied = {col: obs[col] / 25 - tr[col] for col in c6} print("Observations/25 years (1998-2022) minus treated regions:", implied) if n_ctrl: issues.append({"issue": f"Internal inconsistency in {DCR_ID}: Tables 6/7/A19 state 608 never-treated control regions, but the " "Table A19 note says 640 never-treated regions. The printed observation counts (16,175 / 16,875 / 17,250) " "equal 25 years x (640 + 7 / 35 / 50 treated regions), i.e. are consistent with 640 rather than 608 controls " f"(implied controls: {implied}). The 25-year panel length (1998-2022) is our assumption.", "handling": "Reported dcr_obs_5km = 17250 exactly as printed in Table 6; flagged, not corrected."}) issues.append({"issue": f"Text-vs-table wording in {DCR_ID}: the abstract says fatal outcomes 'fall by 35%', while Table 7 " "5 km fatal ATT is -0.354 log points (exp(-0.354)-1 = -29.8%, which Section 6.1 reports as ~30%); the text also " "calls HIV estimates 'close to zero' although the 5 km HIV ATT is -0.340 (insignificant).", "handling": "No effect on extracted keys; values taken from tables only."}) issues.append({"issue": f"In {RI_ID} the results paragraph prints the Weeks-on-UI CI as magnitudes '[1.75 weeks, 2.24 weeks]' " "for a 'reduction', whereas Table 2 prints the signed interval [-2.24, -1.75].", "handling": "Used Table 2 values with sign as printed (lower bound -2.24, upper bound -1.75)."}) issues.append({"issue": f"In {RI_ID} the printed wage-ATE CI [166, 2,141] is very slightly asymmetric around $1,153 " "(987 below, 988 above), consistent with rounding of printed values; the CI method ('described in the main text') is not " "in the extract.", "handling": "Derived SE from the full CI width/(2*1.96); implied z = 1153/SE gives p close to the printed p=0.022 (checked)."}) if star_mismatch: issues.append({"issue": "Significance stars not matching a normal approximation of coef/SE: " + "; ".join(star_mismatch), "handling": "Stars reported as printed; flagged only."}) z_w = R["ri_wage_ate"] / R["ri_wage_ate_se_derived"] p_w = 2 * (1 - norm.cdf(z_w)) print(f"Check: wage ATE z = {z_w:.3f}, implied two-sided p = {p_w:.4f} (printed {w_p})") print(f"Check: HIV 5km z = {R['dcr_hiv_att_5km']/R['dcr_hiv_att_5km_se']:.3f}; stars printed: '{hiv_stars}'") 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"] out = {k: R[k] for k in order} out["locators"] = {k: L[k] for k in order if k in L} out["derivation_notes"] = { "dcr_att_5km_ci": "Normal approximation, critical value 1.96: CI = ATT +/- 1.96*SE = -0.431 +/- 1.96*0.210 " "(Table 6 printed values; rounded to 4 decimals). Assumes asymptotic normality of the CS-DR ATT " "and that the printed clustered SE is the relevant SE.", "ri_wage_ate_se": "Normal approximation, critical value 1.96: SE = (upper - lower)/(2*1.96) = (2141 - 166)/3.92 " "(Table 2 printed CI; rounded to 2 decimals). Assumes the printed CI is symmetric and normal-based; with " "N = 23,549 a t critical value would be indistinguishable from 1.96.", } out["task_id"] = "F4-conf-01" out["feasibility"] = "complete" out["feasibility_reason"] = None out["data_issues"] = issues out["n_input_rows"] = {"README.md": wc_l(P_README), f"{DCR_ID}.txt": wc_l(P_DCR), f"{RI_ID}.txt": wc_l(P_RI)} out["n_analysed"] = len(order) out["random_seed"] = None os.makedirs(os.path.join(ROOT, "out"), exist_ok=True) with open(os.path.join(ROOT, "out", "results.json"), "w") as f: json.dump(out, f, indent=2) print(json.dumps({k: out[k] for k in order}, indent=1))