1"""Replication script for F4-conf-01. Parses the two source extracts in data/sources/,2prints the matched source lines, and writes out/results.json. Deterministic; no randomness."""3import json, os, re, math4from scipy.stats import norm56ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))7DCR_ID, RI_ID = "arxiv-2609.27732v1", "arxiv-2606.14621v1"8P_DCR = os.path.join(ROOT, "data", "sources", DCR_ID + ".txt")9P_RI = os.path.join(ROOT, "data", "sources", RI_ID + ".txt")10P_README = os.path.join(ROOT, "data", "README.md")1112def read(p):13 with open(p, encoding="utf-8") as f:14 return f.read()1516def wc_l(p):17 return read(p).count("\n")1819dcr = read(P_DCR).splitlines()20ri = read(P_RI).splitlines()2122def block(lines, start_pat, end_pat):23 i0 = next(i for i, l in enumerate(lines) if l.startswith(start_pat))24 i1 = next(i for i in range(i0 + 1, len(lines)) if lines[i].startswith(end_pat))25 return lines[i0:i1]2627def cells(line):28 return [c.strip() for c in line.split("|")]2930def parse_reg_table(blk, header_line_idx_pred):31 """Return header cols and dict row -> (coef cells, se cells or None)."""32 nonempty = [l for l in blk if l.strip()]33 hdr = next(l for l in nonempty if header_line_idx_pred(l))34 cols = cells(hdr)[1:]35 rows = {}36 for k, l in enumerate(nonempty):37 c = cells(l)38 if l is hdr or "|" not in l or c[0] == "":39 continue40 se = None41 if k + 1 < len(nonempty) and cells(nonempty[k + 1])[0] == "" and nonempty[k + 1].strip().startswith("| ("):42 se = cells(nonempty[k + 1])[1:]43 print(" matched:", l, "\n ", nonempty[k + 1])44 else:45 print(" matched:", l)46 rows[c[0]] = (c[1:], se)47 return cols, rows4849def num(s):50 m = re.match(r"^\(?(-?[\d,]*\.?\d+)\)?(\**)", s.replace("$", "").strip())51 return m.group(1).replace(",", ""), m.group(2)5253def as_printed(s):54 v, _ = num(s)55 return int(v) if re.fullmatch(r"-?\d+", v) else float(v)5657print("== DCR Table 6 ==")58t6 = block(dcr, "Table 6:", "Note:")59c6, r6 = parse_reg_table(t6, lambda l: l.strip().startswith("| 1 km"))60print("== DCR Table 7 ==")61t7 = block(dcr, "Table 7:", "Note:")62c7, r7 = parse_reg_table(t7, lambda l: l.startswith("Outcome |"))63print("== DCR Table A19 ==")64ta = block(dcr, "Table A19:", "Note:")65ca, ra = parse_reg_table(ta, lambda l: l.strip().startswith("| CS-DR"))6667def g(cols, rows, row, col, se=False):68 j = cols.index(col)69 c, s = rows[row]70 return (s if se else c)[j]7172def loc(src, tab, row, col, extra=""):73 return f'{src} | Table {tab} | row "{row}" | column "{col}"' + extra7475R, L, issues = {}, {}, []7677def put(key, cols, rows, tab, row, col, se=False):78 raw = g(cols, rows, row, col, se)79 R[key] = as_printed(raw)80 L[key] = loc(DCR_ID, tab, row, col, " (standard error, line below coefficient)" if se else "")8182put("dcr_att_5km", c6, r6, 6, "ATT", "5 km")83put("dcr_att_5km_se", c6, r6, 6, "ATT", "5 km", se=True)84put("dcr_att_3km", c6, r6, 6, "ATT", "3 km")85put("dcr_att_3km_se", c6, r6, 6, "ATT", "3 km", se=True)86put("dcr_obs_5km", c6, r6, 6, "Observations", "5 km")87put("dcr_narcotics_att_3km", c7, r7, 7, "Narcotics", "3 km")88put("dcr_narcotics_att_3km_se", c7, r7, 7, "Narcotics", "3 km", se=True)89put("dcr_fatal_att_1km", c7, r7, 7, "Drug-related (fatal)", "1 km")90put("dcr_fatal_att_1km_se", c7, r7, 7, "Drug-related (fatal)", "1 km", se=True)91put("dcr_hiv_att_5km", c7, r7, 7, "HIV", "5 km")92put("dcr_hiv_att_5km_se", c7, r7, 7, "HIV", "5 km", se=True)9394# HIV p-value: search the whole extract for any printed p-value near HIV95hiv_p = [l for l in dcr if "HIV" in l and re.search(r"p\s*[=<]\s*0?\.\d", l)]96print("Lines mentioning HIV with an exact p-value:", hiv_p)97R["dcr_hiv_att_5km_p_value"] = None98L["dcr_hiv_att_5km_p_value"] = (f"not reported in {DCR_ID} (Table 7 prints only the coefficient, SE in parentheses "99 "and significance stars; no exact p-value is printed anywhere in the extract)")100hiv_stars = num(g(c7, r7, "HIV", "5 km"))[1]101R["dcr_hiv_att_5km_significant_10pct"] = len(hiv_stars) >= 1102103# derived 95% CI for Table 6 5km ATT104Z = 1.96105b, se = R["dcr_att_5km"], R["dcr_att_5km_se"]106R["dcr_att_5km_ci_low"] = round(b - Z * se, 4)107R["dcr_att_5km_ci_high"] = round(b + Z * se, 4)108for k in ("dcr_att_5km_ci_low", "dcr_att_5km_ci_high"):109 L[k] = f'{DCR_ID} | Table 6 | row "ATT" | column "5 km" | derived from printed ATT (-0.431) and SE (0.210)'110111# ---------- RI ----------112print("== RI Table 2 ==")113t2 = block(ri, "Outcome | Treatment | Control", "Table 2:")114hdr2 = cells(t2[0])115rows2 = {}116for l in t2[1:]:117 if l.strip():118 print(" matched:", l)119 c = cells(l); rows2[c[0]] = c[1:]120cols2 = hdr2[1:]121ATE = "Average treatment effect estimate"122123def parse_effect(s):124 m = re.match(r"^(-?)\$?(-?[\d,]*\.?\d+)\s*(pp|weeks|wks)?;\s*p\s*([=<])\s*([\d.]+);\s*95% CI:\s*\[([^,\]]+(?:,\d{3})*),\s*([^\]]+)\]$", s)125 # split CI robustly: bounds may contain thousands separators ("2,141")126 est = m.group(1) + m.group(2).replace(",", "")127 ci = re.search(r"\[(.*)\]", s).group(1)128 parts = re.split(r",\s+", ci)129 lo, hi = [p.replace(",", "") for p in parts]130 return est, m.group(4) + m.group(5), lo, hi131132def tonum(v):133 return int(v) if re.fullmatch(r"-?\d+", v) else float(v)134135w_est, w_p, w_lo, w_hi = parse_effect(rows2["Wage"][cols2.index(ATE)])136R["ri_wage_ate"] = tonum(w_est)137R["ri_wage_ate_ci_low"], R["ri_wage_ate_ci_high"] = tonum(w_lo), tonum(w_hi)138for k, what in (("ri_wage_ate", "point estimate"), ("ri_wage_ate_ci_low", "95% CI lower bound"),139 ("ri_wage_ate_ci_high", "95% CI upper bound")):140 L[k] = f'{RI_ID} | Table 2 | row "Wage" | column "{ATE}" ({what})'141se_printed = any(re.search(r"\bSE\b|standard error", l, re.I) for l in t2)142R["ri_wage_ate_se_reported"] = bool(se_printed)143R["ri_wage_ate_se_derived"] = round((R["ri_wage_ate_ci_high"] - R["ri_wage_ate_ci_low"]) / (2 * Z), 2)144L["ri_wage_ate_se_derived"] = (f'{RI_ID} | Table 2 | row "Wage" | column "{ATE}" | derived from printed 95% CI '145 f'[166, 2,141]; SE not reported in {RI_ID}')146cm = rows2["Wage"][cols2.index("Control")]147R["ri_control_mean_wage"] = tonum(re.match(r"\$([\d,]+)", cm).group(1).replace(",", ""))148L["ri_control_mean_wage"] = f'{RI_ID} | Table 2 | row "Wage" | column "Control" (mean; SD printed separately)'149r_est, r_p, r_lo, r_hi = parse_effect(rows2["Reemployment"][cols2.index(ATE)])150R["ri_reemployment_ate_pp"] = tonum(r_est)151L["ri_reemployment_ate_pp"] = f'{RI_ID} | Table 2 | row "Reemployment" | column "{ATE}" (point estimate)'152u_est, u_p, u_lo, u_hi = parse_effect(rows2["Weeks on UI"][cols2.index(ATE)])153R["ri_weeks_ui_ate"] = tonum(u_est)154R["ri_weeks_ui_ate_ci_low"], R["ri_weeks_ui_ate_ci_high"] = tonum(u_lo), tonum(u_hi)155for k, what in (("ri_weeks_ui_ate", "point estimate"), ("ri_weeks_ui_ate_ci_low", "95% CI lower bound"),156 ("ri_weeks_ui_ate_ci_high", "95% CI upper bound")):157 L[k] = f'{RI_ID} | Table 2 | row "Weeks on UI" | column "{ATE}" ({what})'158159# ---------- consistency checks ----------160def stars_expected(b, s):161 p = 2 * (1 - norm.cdf(abs(b / s)))162 return ("***" if p < .01 else "**" if p < .05 else "*" if p < .10 else ""), p163164star_mismatch = []165for tab, cols, rows in (("6", c6, r6), ("7", c7, r7), ("A19", ca, ra)):166 for row, (c, s) in rows.items():167 if s is None:168 continue169 for j, col in enumerate(cols):170 if j >= len(c) or not c[j]:171 continue172 bv, st = num(c[j]); sv, _ = num(s[j])173 exp, p = stars_expected(float(bv), float(sv))174 if exp != st:175 star_mismatch.append(f"Table {tab} {row}/{col}: printed '{st}' vs normal-approx '{exp}' (p={p:.3f})")176print("Star/SE mismatches:", star_mismatch)177178a19_se = as_printed(g(ca, ra, "ATT", "CS-DR", se=True))179if a19_se != R["dcr_att_5km_se"]:180 issues.append({"issue": f"Internal inconsistency in {DCR_ID}: the 5 km CS-DR overall ATT (-0.431**) has SE "181 f"{R['dcr_att_5km_se']} in Table 6 but SE {a19_se} in Appendix Table A19 (same estimator, same 5 km buffer).",182 "handling": "Used Table 6 as the key definition requires (SE 0.210); derived CI uses 0.210. "183 "With 0.218 the CI would be [" + f"{b - Z * a19_se:.4f}, {b + Z * a19_se:.4f}" + "]."})184n_ctrl = [l for l in dcr if "640 regions" in l]185obs = {col: as_printed(g(c6, r6, "Observations", col)) for col in c6}186tr = {col: as_printed(g(c6, r6, "Treated regions", col)) for col in c6}187implied = {col: obs[col] / 25 - tr[col] for col in c6}188print("Observations/25 years (1998-2022) minus treated regions:", implied)189if n_ctrl:190 issues.append({"issue": f"Internal inconsistency in {DCR_ID}: Tables 6/7/A19 state 608 never-treated control regions, but the "191 "Table A19 note says 640 never-treated regions. The printed observation counts (16,175 / 16,875 / 17,250) "192 "equal 25 years x (640 + 7 / 35 / 50 treated regions), i.e. are consistent with 640 rather than 608 controls "193 f"(implied controls: {implied}). The 25-year panel length (1998-2022) is our assumption.",194 "handling": "Reported dcr_obs_5km = 17250 exactly as printed in Table 6; flagged, not corrected."})195issues.append({"issue": f"Text-vs-table wording in {DCR_ID}: the abstract says fatal outcomes 'fall by 35%', while Table 7 "196 "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 "197 "calls HIV estimates 'close to zero' although the 5 km HIV ATT is -0.340 (insignificant).",198 "handling": "No effect on extracted keys; values taken from tables only."})199issues.append({"issue": f"In {RI_ID} the results paragraph prints the Weeks-on-UI CI as magnitudes '[1.75 weeks, 2.24 weeks]' "200 "for a 'reduction', whereas Table 2 prints the signed interval [-2.24, -1.75].",201 "handling": "Used Table 2 values with sign as printed (lower bound -2.24, upper bound -1.75)."})202issues.append({"issue": f"In {RI_ID} the printed wage-ATE CI [166, 2,141] is very slightly asymmetric around $1,153 "203 "(987 below, 988 above), consistent with rounding of printed values; the CI method ('described in the main text') is not "204 "in the extract.",205 "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)."})206if star_mismatch:207 issues.append({"issue": "Significance stars not matching a normal approximation of coef/SE: " + "; ".join(star_mismatch),208 "handling": "Stars reported as printed; flagged only."})209210z_w = R["ri_wage_ate"] / R["ri_wage_ate_se_derived"]211p_w = 2 * (1 - norm.cdf(z_w))212print(f"Check: wage ATE z = {z_w:.3f}, implied two-sided p = {p_w:.4f} (printed {w_p})")213print(f"Check: HIV 5km z = {R['dcr_hiv_att_5km']/R['dcr_hiv_att_5km_se']:.3f}; stars printed: '{hiv_stars}'")214215order = ["dcr_att_5km", "dcr_att_5km_se", "dcr_att_3km", "dcr_att_3km_se", "dcr_obs_5km",216 "dcr_narcotics_att_3km", "dcr_narcotics_att_3km_se", "dcr_fatal_att_1km", "dcr_fatal_att_1km_se",217 "dcr_hiv_att_5km", "dcr_hiv_att_5km_se", "dcr_hiv_att_5km_p_value", "dcr_hiv_att_5km_significant_10pct",218 "dcr_att_5km_ci_low", "dcr_att_5km_ci_high", "ri_wage_ate", "ri_wage_ate_ci_low", "ri_wage_ate_ci_high",219 "ri_wage_ate_se_reported", "ri_wage_ate_se_derived", "ri_control_mean_wage", "ri_reemployment_ate_pp",220 "ri_weeks_ui_ate", "ri_weeks_ui_ate_ci_low", "ri_weeks_ui_ate_ci_high"]221out = {k: R[k] for k in order}222out["locators"] = {k: L[k] for k in order if k in L}223out["derivation_notes"] = {224 "dcr_att_5km_ci": "Normal approximation, critical value 1.96: CI = ATT +/- 1.96*SE = -0.431 +/- 1.96*0.210 "225 "(Table 6 printed values; rounded to 4 decimals). Assumes asymptotic normality of the CS-DR ATT "226 "and that the printed clustered SE is the relevant SE.",227 "ri_wage_ate_se": "Normal approximation, critical value 1.96: SE = (upper - lower)/(2*1.96) = (2141 - 166)/3.92 "228 "(Table 2 printed CI; rounded to 2 decimals). Assumes the printed CI is symmetric and normal-based; with "229 "N = 23,549 a t critical value would be indistinguishable from 1.96.",230}231out["task_id"] = "F4-conf-01"232out["feasibility"] = "complete"233out["feasibility_reason"] = None234out["data_issues"] = issues235out["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)}236out["n_analysed"] = len(order)237out["random_seed"] = None238os.makedirs(os.path.join(ROOT, "out"), exist_ok=True)239with open(os.path.join(ROOT, "out", "results.json"), "w") as f:240 json.dump(out, f, indent=2)241print(json.dumps({k: out[k] for k in order}, indent=1))242