"""Locked analysis per run3/prereg_side.md (sha d2585e4a...). One pass -> side_results.json.""" import json, math, random, statistics, sys sys.path.insert(0, ".") from collections import defaultdict random.seed(20260823) R = lambda tag: random.Random(20260823 + hash(tag) % 1000) rows = json.load(open("run3/review_rows.json")) tlen = {x["review_id"]: x["text_len"] for x in rows} rev = {x["review_id"]: x for x in json.load(open("run3/rev_table.json"))} out = {} # ---------- build analysis frame ---------- byp = defaultdict(list) for r in rev.values(): if r["rigour"] is not None: byp[r["paper"]].append(r["review_id"]) frame = [] for rid, r in rev.items(): if r["rigour"] is None: continue peers = [x for x in byp[r["paper"]] if x != rid] if not peers: continue loo = statistics.mean(rev[x]["rigour"] for x in peers) frame.append({ "rid": rid, "paper": r["paper"], "reviewer": r["reviewer"], "rig": r["rigour"], "nov": r["novelty"], "cla": r["clarity"], "sig": r["significance"], "absdev": abs(r["rigour"] - loo), "absdev_nov": abs(r["novelty"] - statistics.mean(rev[x]["novelty"] for x in peers)), "absdev_cla": abs(r["clarity"] - statistics.mean(rev[x]["clarity"] for x in peers)), "absdev_sig": abs(r["significance"] - statistics.mean(rev[x]["significance"] for x in peers)), "q": r["quality_score"], "rated": (r.get("quality_score") is not None and r.get("ratings_count", 0) >= 1), "text_len": tlen.get(rid, 0), "same_op": bool(r.get("same_operator")), "n_other": len(peers), "created": r["created_at"], "ratings_count": r.get("ratings_count", 0), }) frame.sort(key=lambda x: (x["paper"], x["created"])) for p in set(f["paper"] for f in frame): idx = [i for i, f in enumerate(frame) if f["paper"] == p] for pos, i in enumerate(idx): frame[i]["position"] = pos t0 = min(f["created"] for f in frame) HOURS = lambda ts: (__import__("datetime").datetime.fromisoformat(ts.replace("Z", "+00:00")) - __import__("datetime").datetime.fromisoformat(t0.replace("Z", "+00:00"))).total_seconds() / 3600 for f in frame: f["hours"] = HOURS(f["created"]) rated = [f for f in frame if f["rated"]] print("frame:", len(frame), "rated:", len(rated)) out["n_frame"] = len(frame); out["n_rated"] = len(rated) def pearson(xs, ys): n = len(xs); mx = sum(xs)/n; my = sum(ys)/n cov = sum((x-mx)*(y-my) for x, y in zip(xs, ys)) vx = sum((x-mx)**2 for x in xs); vy = sum((y-my)**2 for y in ys) return cov/(vx*vy)**.5 def zscore(xs): m = sum(xs)/len(xs); v = sum((x-m)**2 for x in xs)/len(xs) s = v**.5 if v > 0 else 1.0 return [(x-m)/s for x in xs] # ---------- A1a replication anchor ---------- r_a1a = pearson([f["absdev"] for f in rated], [f["q"] for f in rated]) out["A1a_r_absdev_quality"] = round(r_a1a, 4) print("A1a r(absdev,q) =", round(r_a1a, 4)) # ---------- A1b PRIMARY: OLS with reviewer FE + text-len control, cluster-robust ---------- def ols_fe(rows_): """q ~ b1*z_absdev + b2*z_log_textlen + reviewer FE (within demeaning on all vars).""" ys = [f["q"] for f in rows_] xa = zscore([f["absdev"] for f in rows_]) xb = zscore([math.log1p(f["text_len"]) for f in rows_]) groups = defaultdict(list) for i, f in enumerate(rows_): groups[f["reviewer"]].append(i) yd = ys[:]; ad = xa[:]; bd = xb[:] for g, idxs in groups.items(): if len(idxs) < 2: continue my = sum(ys[i] for i in idxs)/len(idxs) ma = sum(xa[i] for i in idxs)/len(idxs) mb = sum(xb[i] for i in idxs)/len(idxs) for i in idxs: yd[i] = ys[i]-my; ad[i] = xa[i]-ma; bd[i] = xb[i]-mb # regress yd on [ad, bd] (+ intercept) n = len(yd) X = [[1.0, ad[i], bd[i]] for i in range(n)] # normal equations 3x3 XtX = [[sum(X[i][a]*X[i][b] for i in range(n)) for b in range(3)] for a in range(3)] Xty = [sum(X[i][a]*yd[i] for i in range(n)) for a in range(3)] # gaussian elimination import copy M = [row[:] + [Xty[j]] for j, row in enumerate(XtX)] for col in range(3): piv = max(range(col, 3), key=lambda rr: abs(M[rr][col])) M[col], M[piv] = M[piv], M[col] pv = M[col][col] if abs(pv) < 1e-12: return None, None, groups M[col] = [v/pv for v in M[col]] for rr in range(3): if rr != col and M[rr][col] != 0: fac = M[rr][col] M[rr] = [v - fac*w for v, w in zip(M[rr], M[col])] beta = [M[j][3] for j in range(3)] resid = [yd[i] - (beta[0] + beta[1]*ad[i] + beta[2]*bd[i]) for i in range(n)] dof = n - 3 - (len(groups)-1) # cluster-robust vcov for b1 (CR1) meat = [0.0]*9 for g, idxs in groups.items(): s = [0.0, 0.0, 0.0] for i in idxs: wgt = math.sqrt(len(groups)/(len(groups)-1)) * (n-1)/dof if len(groups) > 1 else 1 ui = resid[i] s[0] += X[i][0]*ui; s[1] += X[i][1]*ui; s[2] += X[i][2]*ui for a in range(3): for bb in range(3): meat[a*3+bb] += s[a]*s[bb] # (XtX)^-1 via elimination on identity inv = [[1.0 if i == j else 0.0 for j in range(3)] for i in range(3)] A = [row[:] for row in XtX] for col in range(3): piv = max(range(col, 3), key=lambda rr: abs(A[rr][col])) A[col], A[piv] = A[piv], A[col] inv[col], inv[piv] = inv[piv], inv[col] pv = A[col][col] A[col] = [v/pv for v in A[col]]; inv[col] = [v/pv for v in inv[col]] for rr in range(3): if rr != col and A[rr][col] != 0: fac = A[rr][col] A[rr] = [v - fac*w for v, w in zip(A[rr], A[col])] inv[rr] = [v - fac*w for v, w in zip(inv[rr], inv[col])] var_b1 = sum(inv[1][a]*meat[a*3+bb]*inv[bb][1] for a in range(3) for bb in range(3)) se = var_b1**.5 t = beta[1]/se if se > 0 else float('nan') return beta[1], {"se_cluster": se, "t": t, "n": n, "n_reviewers": len(groups)}, groups b1, info, _ = ols_fe(rated) out["A1b_b1"] = round(b1, 4); out["A1b_info"] = info print("A1b b1 =", round(b1, 4), "| cluster SE:", round(info["se_cluster"], 4), "| t:", round(info["t"], 2)) # bootstrap CI over reviewers for b1 (cluster bootstrap, B=10000) rng = R("boot") revs_list = sorted(set(f["reviewer"] for f in rated)) B = 10000 bs = [] by_rev = defaultdict(list) for f in rated: by_rev[f["reviewer"]].append(f) import warnings warnings.filterwarnings("ignore") for _ in range(B): sample = [] for _ in revs_list: sample.extend(by_rev[rng.choice(revs_list)]) try: bb, _, _ = ols_fe(sample) if bb is not None: bs.append(bb) except Exception: pass bs.sort() ci = [bs[int(0.025*len(bs))], bs[int(0.975*len(bs))-1]] out["A1b_boot_ci95"] = [round(ci[0], 4), round(ci[1], 4)] print("A1b cluster-bootstrap CI95:", [round(c, 4) for c in ci]) # ---------- H2 robustness variants ---------- rob = {} def variant(name, rows_, add=None): rr = rows_ if add == "drop_sameop": rr = [f for f in rows_ if not f["same_op"]] b, inf, _ = ols_fe(rr) rob[name] = round(b, 4) variant("base", rated) variant("drop_same_operator", rated, "drop_sameop") def ols_fe_extra(rows_, extra): qs = [f["q"] for f in rows_]; es = [float(f[extra]) for f in rows_] mq = sum(qs)/len(qs); me = sum(es)/len(es) qr = [v-mq for v in qs]; er = [v-me for v in es] bee = sum(a*b for a, b in zip(er, qr))/sum(v*v for v in er) qres = [qr[i]-bee*er[i]+mq for i in range(len(qs))] rows2 = [dict(f) for f in rows_] for i, f in enumerate(rows2): f["q"] = qres[i] b, inf, _ = ols_fe(rows2) return b rob["add_same_operator_dummy"] = round(ols_fe_extra(rated, "same_op"), 4) rob["add_z_log_n_other"] = round(ols_fe_extra(rated, "n_other"), 4) rob["add_time_control"] = round(ols_fe_extra(rated, "hours"), 4) # EB shrinkage variant: shrink within-reviewer demeaning toward grand mean by k/(k+5) ys = [f["q"] for f in rated]; gm = sum(ys)/len(ys) groups = defaultdict(list) for i, f in enumerate(rated): groups[f["reviewer"]].append(i) shr = [] for g, idxs in groups.items(): k = len(idxs) w = k/(k+5) mg = sum(ys[i] for i in idxs)/k for i in idxs: shr.append(gm + w*(ys[i]-mg)) rows_shr = [dict(f) for f in rated] for i, f in enumerate(rows_shr): f["q"] = shr[i] b_eb, _, _ = ols_fe(rows_shr) rob["empirical_bayes_shrinkage"] = round(b_eb, 4) out["H2_robustness_b1"] = rob print("H2 robustness:", json.dumps(rob)) json.dump(out, open("run3/side_results.json", "w"), indent=1) print("saved stage 1")