1"""F4-conf-01: extract printed values from two source extracts and derive two quantities.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/; writes out/results.json. No randomness."""4import json, math, re, os56DCR_ID, RI_ID = "arxiv-2609.27732v1", "arxiv-2606.14621v1"7DCR = f"data/sources/{DCR_ID}.txt"8RI = f"data/sources/{RI_ID}.txt"9Z = 1.96 # normal critical value for two-sided 95%101112def read(p):13 with open(p, encoding="utf-8") as f:14 return f.read()151617def nlines(p): # equivalent of wc -l18 return read(p).count("\n")192021def cells(line):22 return [c.strip() for c in line.split("|")]232425def num(tok):26 """Printed number token -> float/int, keeping printed value (strip stars, $, commas, %)."""27 t = re.sub(r"[*$,%()]", "", tok).strip()28 return int(t) if re.fullmatch(r"-?\d+", t) else float(t)293031def decimals(tok):32 t = re.sub(r"[*$,%()]", "", tok).strip()33 return len(t.split(".")[1]) if "." in t else 0343536def table_block(text, start_marker, end_marker):37 i = text.index(start_marker)38 j = text.index(end_marker, i)39 return [l for l in text[i:j].splitlines() if l.strip()]404142def regression_table(lines, header_prefix):43 """Parse a regression table rendered as row | cells, SE on the next line '| (se) | ...'."""44 hdr = next(l for l in lines if l.startswith(header_prefix))45 cols = [c for c in cells(hdr)[1:]]46 rows, last = {}, None47 for l in lines:48 c = cells(l)49 if l.startswith("|") and c[1].startswith("("):50 rows[last]["se"] = c[1:]51 rows[last]["se_line"] = l52 elif len(c) > 1 and l is not hdr and not l.startswith("Note"):53 last = c[0]54 rows[last] = {"est": c[1:], "line": l}55 return cols, rows565758dcr = read(DCR)59ri = read(RI)60matched = []6162# ---- DCR Table 663t6 = table_block(dcr, "Table 6: Overall Average", "Note: The table reports the overall ATT")64cols6, rows6 = regression_table(t6, "| 1 km")65cols6 = cols6 # ['1 km','3 km','5 km']66assert cols6 == ["1 km", "3 km", "5 km"], cols667# ---- DCR Table 768t7 = table_block(dcr, "Table 7: Overall ATT by Outcome", "Note: Each cell reports")69cols7, rows7 = regression_table(t7, "Outcome |")70assert cols7 == ["1 km", "3 km", "5 km"], cols771# ---- DCR Table A19 (consistency check only)72ta19 = table_block(dcr, "Table A19:", "Note: All specifications")73colsA, rowsA = regression_table(ta19, "| CS-DR")7475def t_est(rows, cols, r, c):76 k = cols.index(c)77 matched.append(rows[r]["line"])78 return rows[r]["est"][k]7980def t_se(rows, cols, r, c):81 k = cols.index(c)82 matched.append(rows[r]["se_line"])83 return rows[r]["se"][k]8485R = {}86L = {}87def loc(src, tab, row, col, extra=""):88 return f'{src} | Table {tab} | row "{row}" | column "{col}"' + extra8990spec = [91 ("dcr_att_5km", rows6, cols6, "ATT", "5 km", "6", "est"),92 ("dcr_att_5km_se", rows6, cols6, "ATT", "5 km", "6", "se"),93 ("dcr_att_3km", rows6, cols6, "ATT", "3 km", "6", "est"),94 ("dcr_att_3km_se", rows6, cols6, "ATT", "3 km", "6", "se"),95 ("dcr_obs_5km", rows6, cols6, "Observations", "5 km", "6", "est"),96 ("dcr_narcotics_att_3km", rows7, cols7, "Narcotics", "3 km", "7", "est"),97 ("dcr_narcotics_att_3km_se", rows7, cols7, "Narcotics", "3 km", "7", "se"),98 ("dcr_fatal_att_1km", rows7, cols7, "Drug-related (fatal)", "1 km", "7", "est"),99 ("dcr_fatal_att_1km_se", rows7, cols7, "Drug-related (fatal)", "1 km", "7", "se"),100 ("dcr_hiv_att_5km", rows7, cols7, "HIV", "5 km", "7", "est"),101 ("dcr_hiv_att_5km_se", rows7, cols7, "HIV", "5 km", "7", "se"),102]103tokens = {}104for key, rows, cols, r, c, tab, kind in spec:105 tok = t_est(rows, cols, r, c) if kind == "est" else t_se(rows, cols, r, c)106 tokens[key] = tok107 R[key] = num(tok)108 L[key] = loc(DCR_ID, tab, r, c, " (standard error in parentheses on the line below the coefficient)" if kind == "se" else "")109110hiv_tok = tokens["dcr_hiv_att_5km"]111R["dcr_hiv_att_5km_p_value"] = None # Table 7 prints only significance stars, no p-values112L["dcr_hiv_att_5km_p_value"] = (f"not reported in {DCR_ID} (Table 7 prints only significance stars; "113 f'cell row "HIV" column "5 km" = "{hiv_tok}" carries no star)')114R["dcr_hiv_att_5km_significant_10pct"] = "*" in hiv_tok115116# derived CI117b, se = R["dcr_att_5km"], R["dcr_att_5km_se"]118R["dcr_att_5km_ci_low"] = round(b - Z * se, 4)119R["dcr_att_5km_ci_high"] = round(b + Z * se, 4)120L["dcr_att_5km_ci_low"] = L["dcr_att_5km_ci_high"] = (121 f'{DCR_ID} | Table 6 | row "ATT" | column "5 km" | derived: estimate -/+ 1.96 x printed SE')122123# ---- RI Table 2124i2 = ri.index("\nOutcome | Treatment | Control |") + 1125t2 = [l for l in ri[i2:].splitlines()[:7]]126hdr2 = cells(t2[0])127rows2 = {cells(l)[0]: dict(zip(hdr2, cells(l)), line=l) for l in t2[1:]}128ATE = "Average treatment effect estimate"129130def parse_effect(cell):131 """'$1,153; p=0.022; 95% CI: [166, 2,141]' -> ('1,153', '=0.022', '166', '2,141')."""132 m = re.match(r"^\s*(-?\$?-?[\d,\.]+)\s*(?:pp|weeks|wks)?;\s*p\s*([=<])\s*([\d\.]+);\s*95% CI:\s*\[(.+)\]\s*$", cell)133 if not m:134 raise ValueError(cell)135 lo, hi = [p.strip() for p in m.group(4).split(", ")] # bounds separated by comma+space; thousands commas have no space136 return m.group(1).replace("$", ""), m.group(2) + m.group(3), lo, hi137138wage = rows2["Wage"]; matched.append(wage["line"])139w_est, w_p, w_lo, w_hi = parse_effect(wage[ATE])140R["ri_wage_ate"] = num(w_est)141R["ri_wage_ate_ci_low"] = num(w_lo)142R["ri_wage_ate_ci_high"] = num(w_hi)143for k in ["ri_wage_ate", "ri_wage_ate_ci_low", "ri_wage_ate_ci_high"]:144 L[k] = loc(RI_ID, "2", "Wage", ATE) + {"ri_wage_ate": " (point estimate)",145 "ri_wage_ate_ci_low": " (95% CI lower bound)",146 "ri_wage_ate_ci_high": " (95% CI upper bound)"}[k]147R["ri_wage_ate_se_reported"] = ("SE" in wage[ATE]) or ("(" in wage[ATE])148R["ri_wage_ate_se_derived"] = round((R["ri_wage_ate_ci_high"] - R["ri_wage_ate_ci_low"]) / (2 * Z), 2)149L["ri_wage_ate_se_derived"] = (f'{RI_ID} | Table 2 | row "Wage" | column "{ATE}" | derived: '150 f'(CI upper - CI lower)/(2 x 1.96); SE not reported in {RI_ID}')151ctrl_tok = re.match(r"^\$([\d,]+)", wage["Control"]).group(1)152R["ri_control_mean_wage"] = num(ctrl_tok)153L["ri_control_mean_wage"] = loc(RI_ID, "2", "Wage", "Control") + " (mean; SD in parentheses not used)"154155re_row = rows2["Reemployment"]; matched.append(re_row["line"])156r_est, r_p, r_lo, r_hi = parse_effect(re_row[ATE])157R["ri_reemployment_ate_pp"] = num(r_est)158L["ri_reemployment_ate_pp"] = loc(RI_ID, "2", "Reemployment", ATE) + " (point estimate, pp)"159160ui = rows2["Weeks on UI"]; matched.append(ui["line"])161u_est, u_p, u_lo, u_hi = parse_effect(ui[ATE])162R["ri_weeks_ui_ate"] = num(u_est)163R["ri_weeks_ui_ate_ci_low"] = num(u_lo)164R["ri_weeks_ui_ate_ci_high"] = num(u_hi)165L["ri_weeks_ui_ate"] = loc(RI_ID, "2", "Weeks on UI", ATE) + " (point estimate, weeks)"166L["ri_weeks_ui_ate_ci_low"] = loc(RI_ID, "2", "Weeks on UI", ATE) + " (95% CI lower bound)"167L["ri_weeks_ui_ate_ci_high"] = loc(RI_ID, "2", "Weeks on UI", ATE) + " (95% CI upper bound)"168169# ---- consistency checks on the sources (feed data_issues)170issues = []171a19_se = rowsA["ATT"]["se"][colsA.index("CS-DR")]172a19_est = rowsA["ATT"]["est"][colsA.index("CS-DR")]173if num(a19_se) != R["dcr_att_5km_se"]:174 issues.append({"issue": f"Internal inconsistency in {DCR_ID}: the 5 km CS-DR ATT is {a19_est} in both Table 6 and Table A19, "175 f"but its SE is printed as {tokens['dcr_att_5km_se']} in Table 6 and {a19_se} in Table A19 (same estimator, buffer and window as described).",176 "handling": "Used Table 6 values as the key definitions require; A19 value not used; the derived CI depends on this choice (with 0.218 it would be wider)."})177obs = [num(x) for x in rows6["Observations"]["est"]]178tr = [num(x) for x in rows6["Treated regions"]["est"]]179ctrl = [num(x) for x in rows6["Control regions"]["est"]]180per_region = [o / (t + c) for o, t, c in zip(obs, tr, ctrl)]181implied_ctrl_25 = [o / 25 - t for o, t in zip(obs, tr)]182if "640 regions" in dcr and "608" in rows6["Control regions"]["line"]:183 issues.append({"issue": f"Internal inconsistency in {DCR_ID}: Tables 6/7/A19 and figure notes state 608 never-treated control regions, but the Table A19 note says 640 never-treated comparison regions. "184 f"Observations/(treated+608 control) is non-integer ({', '.join(f'{x:.3f}' for x in per_region)}), whereas Observations/25 - treated = {implied_ctrl_25} for 1/3/5 km, i.e. the printed Observations are arithmetically consistent with 640 control regions x 25 years (1998-2022).",185 "handling": "Reported Observations exactly as printed (17250); the discrepancy is flagged, not resolved; no extracted key depends on the control-region count."})186fatal5 = num(rows7["Drug-related (fatal)"]["est"][2])187pct_fatal5 = 100 * (math.exp(fatal5) - 1)188issues.append({"issue": f"Text/table inconsistency in {DCR_ID}: the abstract says fatal outcomes 'fall by 35%', but Table 7 5 km fatal ATT {fatal5} implies {pct_fatal5:.1f}% (exp(b)-1), and Section 6.1 says 'approximately 30% at 5 km'; 35% appears to read the log-point coefficient as a percent.",189 "handling": "Not used for any key; flagged only."})190issues.append({"issue": f"Interpretive inconsistency in {DCR_ID}: the text calls HIV estimates 'close to zero', whereas Table 7 prints {rows7['HIV']['est'][1]} (3 km) and {rows7['HIV']['est'][2]} (5 km), insignificant but not small in magnitude.",191 "handling": "Table values reported as printed; text characterisation not used."})192pc = rowsA["% change"]["est"]193bad = []194for col, e, p in zip(colsA, rowsA["ATT"]["est"], pc):195 calc = 100 * (math.exp(num(e)) - 1)196 if abs(round(calc, 1) - num(p)) > 0.05:197 bad.append(f"{col}: printed {p}, exp({num(e)})-1 = {calc:.2f}%")198if bad:199 issues.append({"issue": f"Minor rounding discrepancy in {DCR_ID} Table A19 '% change' row vs exp(ATT)-1 of the printed (rounded) ATT: " + "; ".join(bad) + ". Likely computed from unrounded coefficients.",200 "handling": "Not used for any key; flagged only."})201issues.append({"issue": f"Rendering defect in {RI_ID} Section 3 text: 'p ¡ 0.001' (an inverted exclamation mark in place of '<').",202 "handling": "Text only; no key read from that sentence."})203issues.append({"issue": f"Sign convention in {RI_ID}: the text reports the weeks-on-UI CI as magnitudes of a reduction '[1.75 weeks, 2.24 weeks]' whereas Table 2 prints [-2.24, -1.75].",204 "handling": "Used Table 2 values as printed, lower bound first (-2.24, -1.75)."})205issues.append({"issue": f"Table re-setting in {RI_ID}: Tables 2-3 were re-set by the extract editor into one row per outcome (per extract header).",206 "handling": "Accepted the re-set layout; columns identified by the editor-provided header; Table 3 'Difference in means' column equals Table 2 ATE column (checked)."})207t3_line = next(l for l in ri.splitlines() if l.startswith("Wage | $1,153"))208assert cells(t3_line)[1] == wage[ATE]209210# symmetry check of wage CI (for the SE derivation)211half_lo = R["ri_wage_ate"] - R["ri_wage_ate_ci_low"]212half_hi = R["ri_wage_ate_ci_high"] - R["ri_wage_ate"]213z_implied = R["ri_wage_ate"] / R["ri_wage_ate_se_derived"]214p_implied = math.erfc(z_implied / math.sqrt(2))215216derivation_notes = {217 "dcr_att_5km_ci": (f"Normal-approximation Wald interval, critical value 1.96: {b} -/+ 1.96 x {se} "218 f"= [{R['dcr_att_5km_ci_low']}, {R['dcr_att_5km_ci_high']}] (rounded to 4 decimals). Inputs from Table 6 (SE clustered at Medstat region level). "219 f"Assumes asymptotic normality; a t critical value with few clusters (50 treated regions) would widen it slightly. Table A19 prints SE 0.218 for the same spec."),220 "ri_wage_ate_se": (f"SE = (upper - lower)/(2 x 1.96) = ({R['ri_wage_ate_ci_high']} - {R['ri_wage_ate_ci_low']})/3.92 = {R['ri_wage_ate_se_derived']} (rounded to 2 decimals). "221 f"Assumes a symmetric normal-based 95% CI (half-widths {half_lo} and {half_hi}, symmetric up to rounding). "222 f"Check: implied z = {z_implied:.3f}, two-sided p = {p_implied:.4f} vs printed p={w_p.lstrip('=')}."),223}224225n_input = {"README.md": nlines("data/README.md"),226 f"sources/{DCR_ID}.txt": nlines(DCR),227 f"sources/{RI_ID}.txt": nlines(RI)}228229results = {"task_id": "F4-conf-01"}230order = ["dcr_att_5km", "dcr_att_5km_se", "dcr_att_3km", "dcr_att_3km_se", "dcr_obs_5km",231 "dcr_narcotics_att_3km", "dcr_narcotics_att_3km_se", "dcr_fatal_att_1km", "dcr_fatal_att_1km_se",232 "dcr_hiv_att_5km", "dcr_hiv_att_5km_se", "dcr_hiv_att_5km_p_value", "dcr_hiv_att_5km_significant_10pct",233 "dcr_att_5km_ci_low", "dcr_att_5km_ci_high", "ri_wage_ate", "ri_wage_ate_ci_low", "ri_wage_ate_ci_high",234 "ri_wage_ate_se_reported", "ri_wage_ate_se_derived", "ri_control_mean_wage", "ri_reemployment_ate_pp",235 "ri_weeks_ui_ate", "ri_weeks_ui_ate_ci_low", "ri_weeks_ui_ate_ci_high"]236for k in order:237 results[k] = R[k]238results["locators"] = {k: L[k] for k in order if k not in ("dcr_hiv_att_5km_significant_10pct", "ri_wage_ate_se_reported")}239results["derivation_notes"] = derivation_notes240results["feasibility"] = "complete"241results["feasibility_reason"] = None242results["data_issues"] = issues243results["n_input_rows"] = n_input244results["n_analysed"] = len(order)245results["random_seed"] = None246results["printed_tokens"] = {k: tokens[k] for k in tokens}247results["printed_tokens"].update({"ri_wage_ate": wage[ATE], "ri_reemployment_ate_pp": re_row[ATE],248 "ri_weeks_ui_ate": ui[ATE], "ri_control_mean_wage": wage["Control"]})249250os.makedirs("out", exist_ok=True)251with open("out/results.json", "w", encoding="utf-8") as f:252 json.dump(results, f, indent=2, ensure_ascii=False)253254print("Matched source lines:")255for l in dict.fromkeys(matched):256 print(" ", l)257print(json.dumps({k: results[k] for k in order}, indent=1))258print(f"{len(issues)} data issues recorded; wrote out/results.json")259