1"""F4-conf-01: extract printed treatment effects / uncertainty from two arXiv extracts.23Reads only data/; writes out/results.json deterministically. No randomness.4Run from workspace root: ~/wp29/venv/bin/python out/analysis.py5"""6import json, math, os, re, hashlib78DATA = "data"9DCR_ID, RI_ID = "arxiv-2609.27732v1", "arxiv-2606.14621v1"10DCR_F = os.path.join(DATA, "sources", DCR_ID + ".txt")11RI_F = os.path.join(DATA, "sources", RI_ID + ".txt")12Z = 1.96 # normal critical value for two-sided 95%131415def lines(path):16 with open(path, encoding="utf-8") as f:17 return f.read().split("\n")181920def cells(line):21 return [c.strip() for c in line.split("|")]222324def num(s):25 """Parse a printed number: strip stars, $, %, thousands commas, parentheses."""26 s = s.replace("$", "").replace(",", "").replace("*", "").strip().strip("()")27 return int(s) if re.fullmatch(r"-?\d+", s) else float(s)282930def stars(s):31 m = re.search(r"(\*+)\s*$", s.strip())32 return len(m.group(1)) if m else 0333435def table_block(L, title_prefix, end_prefix="Note:"):36 i = next(k for k, l in enumerate(L) if l.startswith(title_prefix))37 j = next(k for k in range(i + 1, len(L)) if L[k].startswith(end_prefix))38 return i, [(k, L[k]) for k in range(i + 1, j) if L[k].strip()]394041def reg_table(L, title_prefix):42 """Parse a regression table: header, coefficient rows followed by SE rows."""43 i, blk = table_block(L, title_prefix)44 hdr = cells(blk[0][1])45 cols = [c for c in hdr if c and c != "Outcome"]46 rows, last = {}, None47 for k, l in blk[1:]:48 c = cells(l)49 if c[0] == "" and last is not None and c[1].startswith("("):50 rows[last]["se"] = dict(zip(cols, c[1:]))51 rows[last]["se_line"] = (k + 1, l)52 else:53 last = c[0]54 rows[last] = {"est": dict(zip(cols, c[1:])), "line": (k + 1, l)}55 return cols, rows565758def show(tag, ln):59 print(f"[{tag}] line {ln[0]}: {ln[1]}")606162def main():63 DL, RL = lines(DCR_F), lines(RI_F)64 res, loc, issues = {}, {}, []6566 def L(src, t, row, col, extra=""):67 return f'{src} | Table {t} | row "{row}" | column "{col}"' + extra6869 # ---------------- DCR Table 6 ----------------70 c6, t6 = reg_table(DL, "Table 6:")71 show("T6 ATT", t6["ATT"]["line"]); show("T6 SE", t6["ATT"]["se_line"])72 show("T6 Obs", t6["Observations"]["line"])73 res["dcr_att_5km"] = num(t6["ATT"]["est"]["5 km"])74 res["dcr_att_5km_se"] = num(t6["ATT"]["se"]["5 km"])75 res["dcr_att_3km"] = num(t6["ATT"]["est"]["3 km"])76 res["dcr_att_3km_se"] = num(t6["ATT"]["se"]["3 km"])77 res["dcr_obs_5km"] = num(t6["Observations"]["est"]["5 km"])78 loc["dcr_att_5km"] = L(DCR_ID, 6, "ATT", "5 km")79 loc["dcr_att_5km_se"] = L(DCR_ID, 6, "ATT", "5 km", " (standard error in parentheses, line below)")80 loc["dcr_att_3km"] = L(DCR_ID, 6, "ATT", "3 km")81 loc["dcr_att_3km_se"] = L(DCR_ID, 6, "ATT", "3 km", " (standard error in parentheses, line below)")82 loc["dcr_obs_5km"] = L(DCR_ID, 6, "Observations", "5 km")8384 # ---------------- DCR Table 7 ----------------85 c7, t7 = reg_table(DL, "Table 7:")86 for r in ("Narcotics", "Drug-related (fatal)", "HIV"):87 show("T7 " + r, t7[r]["line"]); show("T7 SE " + r, t7[r]["se_line"])88 spec = [("dcr_narcotics_att_3km", "Narcotics", "3 km"),89 ("dcr_fatal_att_1km", "Drug-related (fatal)", "1 km"),90 ("dcr_hiv_att_5km", "HIV", "5 km")]91 for key, r, c in spec:92 res[key] = num(t7[r]["est"][c]); loc[key] = L(DCR_ID, 7, r, c)93 res[key + "_se"] = num(t7[r]["se"][c])94 loc[key + "_se"] = L(DCR_ID, 7, r, c, " (standard error in parentheses, line below)")95 hiv_cell = t7["HIV"]["est"]["5 km"]96 # exact p-value: search extract for any printed p-value97 pv_hits = [l for l in DL if re.search(r"\bp\s*=\s*0?\.\d", l)]98 print(f"[DCR] lines with an exact 'p=' value in extract: {len(pv_hits)}")99 res["dcr_hiv_att_5km_p_value"] = None100 loc["dcr_hiv_att_5km_p_value"] = ("not reported in " + DCR_ID + " (Table 7 prints only the estimate, "101 "its SE in parentheses and significance stars; no exact p-value is printed)")102 res["dcr_hiv_att_5km_significant_10pct"] = stars(hiv_cell) >= 1103104 # ---------------- derived CI (DCR) ----------------105 b, s = res["dcr_att_5km"], res["dcr_att_5km_se"]106 res["dcr_att_5km_ci_low"] = round(b - Z * s, 4)107 res["dcr_att_5km_ci_high"] = round(b + Z * s, 4)108 for k in ("dcr_att_5km_ci_low", "dcr_att_5km_ci_high"):109 loc[k] = (f'{DCR_ID} | Table 6 | row "ATT" | column "5 km" | derived from printed ATT '110 f'and SE: ATT {"-" if k.endswith("low") else "+"} 1.96 x SE')111112 # ---------------- RESEA Table 2 ----------------113 i2 = next(k for k, l in enumerate(RL) if l.startswith("Outcome | Treatment | Control"))114 hdr2 = cells(RL[i2])115 t2 = {}116 for k in range(i2 + 1, len(RL)):117 if not RL[k].strip(): break118 c = cells(RL[k]); t2[c[0]] = (dict(zip(hdr2[1:], c[1:])), (k + 1, RL[k]))119 for r in ("Wage", "Reemployment", "Weeks on UI"):120 show("T2 " + r, t2[r][1])121 ATE = "Average treatment effect estimate"122123 def ate_parse(cell):124 m = re.match(r"\s*(-?\$?-?[\d,\.]+)\s*(pp|weeks|wks)?;\s*p\s*([=<])\s*([\d\.]+);\s*95% CI:\s*\[([^\]]+)\]", cell)125 est = num(m.group(1)); lo, hi = [num(x) for x in m.group(5).split(", ")]126 return est, m.group(3) + m.group(4), lo, hi127128 w_est, w_p, w_lo, w_hi = ate_parse(t2["Wage"][0][ATE])129 r_est, r_p, r_lo, r_hi = ate_parse(t2["Reemployment"][0][ATE])130 u_est, u_p, u_lo, u_hi = ate_parse(t2["Weeks on UI"][0][ATE])131 res["ri_wage_ate"], res["ri_wage_ate_ci_low"], res["ri_wage_ate_ci_high"] = w_est, w_lo, w_hi132 res["ri_wage_ate_se_reported"] = bool(re.search(r"\bSE\b|standard error|\(\d", t2["Wage"][0][ATE]))133 res["ri_wage_ate_se_derived"] = round((w_hi - w_lo) / (2 * Z), 2)134 ctrl = re.match(r"\$([\d,]+)", t2["Wage"][0]["Control"]).group(1)135 res["ri_control_mean_wage"] = num(ctrl)136 res["ri_reemployment_ate_pp"] = r_est137 res["ri_weeks_ui_ate"], res["ri_weeks_ui_ate_ci_low"], res["ri_weeks_ui_ate_ci_high"] = u_est, u_lo, u_hi138 loc["ri_wage_ate"] = L(RI_ID, 2, "Wage", ATE, " (point estimate)")139 loc["ri_wage_ate_ci_low"] = L(RI_ID, 2, "Wage", ATE, " (95% CI, lower bound)")140 loc["ri_wage_ate_ci_high"] = L(RI_ID, 2, "Wage", ATE, " (95% CI, upper bound)")141 loc["ri_wage_ate_se_derived"] = L(RI_ID, 2, "Wage", ATE,142 " | derived from printed 95% CI: (upper - lower) / (2 x 1.96)")143 loc["ri_control_mean_wage"] = L(RI_ID, 2, "Wage", "Control", " (mean)")144 loc["ri_reemployment_ate_pp"] = L(RI_ID, 2, "Reemployment", ATE, " (point estimate)")145 loc["ri_weeks_ui_ate"] = L(RI_ID, 2, "Weeks on UI", ATE, " (point estimate)")146 loc["ri_weeks_ui_ate_ci_low"] = L(RI_ID, 2, "Weeks on UI", ATE, " (95% CI, lower bound)")147 loc["ri_weeks_ui_ate_ci_high"] = L(RI_ID, 2, "Weeks on UI", ATE, " (95% CI, upper bound)")148149 # ---------------- consistency checks -> data_issues ----------------150 # (a) Table 6 vs Table A19 SE for the same 5 km CS-DR ATT151 cA, tA = reg_table(DL, "Table A19:")152 a19_att, a19_se = num(tA["ATT"]["est"]["CS-DR"]), num(tA["ATT"]["se"]["CS-DR"])153 show("A19 ATT", tA["ATT"]["line"]); show("A19 SE", tA["ATT"]["se_line"])154 if a19_att == res["dcr_att_5km"] and a19_se != res["dcr_att_5km_se"]:155 issues.append({"issue": f"Same 5 km CS-DR ATT ({a19_att}) printed with SE {t6['ATT']['se']['5 km']} in Table 6 "156 f"but SE {tA['ATT']['se']['CS-DR']} in Table A19 (CS-DR column) of {DCR_ID}.",157 "handling": "Used Table 6 (the table named in the key definitions) for dcr_att_5km_se and the derived CI; "158 "the significance star (**) is the same under either SE."})159 # (b) control-region count: 608 (Tables 6/7, Fig. 6) vs 640 (A19 note); observations arithmetic160 a19_note = next(l for l in DL if "never-treated comparison regions (" in l)161 m = re.search(r"\((\d+) regions\)", a19_note)162 n_ctrl6 = num(t6["Control regions"]["est"]["5 km"])163 obs = {c: num(t6["Observations"]["est"][c]) for c in c6}164 trt = {c: num(t6["Treated regions"]["est"][c]) for c in c6}165 per_reg = {c: obs[c] / (trt[c] + n_ctrl6) for c in c6}166 per_reg640 = {c: obs[c] / (trt[c] + 640) for c in c6}167 print(f"[check] obs/(treated+608) = {per_reg}; obs/(treated+640) = {per_reg640}")168 if m and int(m.group(1)) != n_ctrl6:169 issues.append({"issue": f"Control-region count inconsistent within {DCR_ID}: Tables 6/7 and the Figure 6 note say "170 f"{n_ctrl6} never-treated regions; the Table A19 note says {m.group(1)}. Table 6 observations "171 f"divided by (treated + {n_ctrl6}) are non-integers ({', '.join(f'{v:.2f}' for v in per_reg.values())}), "172 f"whereas divided by (treated + 640) they equal {', '.join(f'{v:g}' for v in per_reg640.values())}.",173 "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)."})174 # (c) abstract: fatal outcomes 'fall by 35%' vs Table 7 5 km implied change175 fat5 = num(t7["Drug-related (fatal)"]["est"]["5 km"]); narc5 = num(t7["Narcotics"]["est"]["5 km"])176 pct = lambda x: 100 * (math.exp(x) - 1)177 print(f"[check] exp(ATT)-1: fatal 5km {pct(fat5):.1f}%, narcotics 5km {pct(narc5):.1f}%, overall 5km {pct(res['dcr_att_5km']):.1f}%")178 if any("fall by 35%" in l for l in DL):179 issues.append({"issue": f"Abstract of {DCR_ID} says fatal outcomes 'fall by 35%', but Table 7 fatal 5 km ATT {fat5} implies "180 f"{pct(fat5):.1f}% (Section 6.1 says ~30% at 5 km); 35% matches the overall 5 km ATT. The text also calls HIV "181 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).",182 "handling": "Narrative statements not used; all values taken from the tables named in the keys."})183 # (d) RESEA text gives weeks-on-UI CI as positive magnitudes184 if any("[1.75 weeks, 2.24 weeks]" in l for l in RL):185 issues.append({"issue": f"{RI_ID} text reports the weeks-on-UI CI as magnitudes '[1.75 weeks, 2.24 weeks]' (a reduction), "186 f"while Table 2 prints [{u_lo}, {u_hi}].",187 "handling": "Used Table 2 signs and order (lower bound first)."})188 # (e) rendering defect189 if any("p ¡ 0.001" in l for l in RL):190 issues.append({"issue": f"Rendering defect in {RI_ID} text: 'p ¡ 0.001' (should read p < 0.001).",191 "handling": "Not used for any key; Table 2 values used."})192 # star / SE consistency (not an issue if consistent)193 from statistics import NormalDist194 nd = NormalDist()195 star_ok = []196 for tab, rows, cols in (("6", t6, c6), ("7", t7, c7)):197 for r, d in rows.items():198 if "se" not in d: continue199 for c in cols:200 e, se_ = num(d["est"][c]), num(d["se"][c])201 p = 2 * (1 - nd.cdf(abs(e / se_)))202 exp_st = 3 if p < .01 else 2 if p < .05 else 1 if p < .10 else 0203 star_ok.append((tab, r, c, round(p, 4), stars(d["est"][c]), exp_st))204 mism = [x for x in star_ok if x[4] != x[5]]205 print(f"[check] star vs normal-approx p from est/SE: {len(star_ok)} cells, {len(mism)} mismatches {mism}")206 hiv_p_approx = [x[3] for x in star_ok if x[1] == "HIV" and x[2] == "5 km"][0]207 # RESEA CI symmetry208 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}")209210 # ---------------- baseline / sensitivity ----------------211 res["baseline_dcr_att_5km_ci_rule2"] = [round(b - 2 * s, 4), round(b + 2 * s, 4)]212 res["baseline_ri_wage_ate_se_rule2"] = round((w_hi - w_lo) / 4, 2)213 res["check_dcr_hiv_att_5km_normal_approx_p_not_printed"] = hiv_p_approx214215 out = {}216 order = ["dcr_att_5km", "dcr_att_5km_se", "dcr_att_3km", "dcr_att_3km_se", "dcr_obs_5km",217 "dcr_narcotics_att_3km", "dcr_narcotics_att_3km_se", "dcr_fatal_att_1km", "dcr_fatal_att_1km_se",218 "dcr_hiv_att_5km", "dcr_hiv_att_5km_se", "dcr_hiv_att_5km_p_value", "dcr_hiv_att_5km_significant_10pct",219 "dcr_att_5km_ci_low", "dcr_att_5km_ci_high", "ri_wage_ate", "ri_wage_ate_ci_low", "ri_wage_ate_ci_high",220 "ri_wage_ate_se_reported", "ri_wage_ate_se_derived", "ri_control_mean_wage", "ri_reemployment_ate_pp",221 "ri_weeks_ui_ate", "ri_weeks_ui_ate_ci_low", "ri_weeks_ui_ate_ci_high"]222 for k in order: out[k] = res[k]223 out["locators"] = {k: loc[k] for k in order if k in loc}224 out["derivation_notes"] = {225 "dcr_att_5km_ci": "Normal-approximation (Wald) 95% CI from Table 6 printed values: ATT +/- 1.96 x SE = "226 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 "227 "CS-DR ATT with the printed clustered SE; a 2.0 rule gives "228 f"{res['baseline_dcr_att_5km_ci_rule2']}.",229 "ri_wage_ate_se": "SE implied by the Table 2 printed 95% CI assuming a symmetric normal (Wald) interval: "230 f"(upper - lower) / (2 x 1.96) = ({w_hi} - ({w_lo})) / 3.92; rounded to 2 decimals. With N = 23,549 a "231 f"t critical value is practically 1.96; a 2.0 rule gives {res['baseline_ri_wage_ate_se_rule2']}. "232 "Printed bounds are rounded to whole dollars, so the derived SE is accurate to about +/-0.2."}233 out["task_id"] = "F4-conf-01"234 out["feasibility"] = "complete"235 out["feasibility_reason"] = None236 out["data_issues"] = issues237 nrows = {}238 for p in (os.path.join(DATA, "README.md"), DCR_F, RI_F):239 with open(p, "rb") as f: nrows[os.path.relpath(p, DATA)] = f.read().count(b"\n")240 out["n_input_rows"] = nrows241 out["n_analysed"] = len(order)242 out["random_seed"] = None243 out["baseline_dcr_att_5km_ci_rule2"] = res["baseline_dcr_att_5km_ci_rule2"]244 out["baseline_ri_wage_ate_se_rule2"] = res["baseline_ri_wage_ate_se_rule2"]245 out["check_dcr_hiv_att_5km_normal_approx_p_not_printed"] = res["check_dcr_hiv_att_5km_normal_approx_p_not_printed"]246 out["input_sha256"] = {os.path.relpath(p, DATA): hashlib.sha256(open(p, "rb").read()).hexdigest()247 for p in (os.path.join(DATA, "README.md"), DCR_F, RI_F)}248 os.makedirs("out", exist_ok=True)249 with open(os.path.join("out", "results.json"), "w") as f:250 json.dump(out, f, indent=2, ensure_ascii=False); f.write("\n")251 print("Wrote out/results.json")252253254if __name__ == "__main__":255 main()256