# Reviewer's sensitivity check, run on the bundle's own results/associations.json (no new NHANES analysis).
# 1) Recomputes the plan-6.5 "difference from published" test: z = (b_rep - b_pub)/sqrt(se_pub^2 + se_rep^2),
#    se_pub from the published 95% CI / 3.92 on the analysis scale, BH over 40; checks it matches the bundle.
# 2) Repeats it with se_pub floored at the bundle's own design-based SE on the paper's cycles (original.se),
#    since the paper itself reports several published intervals narrower than the design supports.
import json, math
A = json.load(open('bundle/results/associations.json'))
def tr(m, x): return math.log(x) if m == 'OR' else x
def p2(z): return math.erfc(abs(z) / math.sqrt(2))
def bh(ps):
    n = len(ps); idx = sorted(range(n), key=lambda i: ps[i]); q = [0]*n; m = 1.0
    for rank in range(n, 0, -1):
        i = idx[rank-1]; m = min(m, ps[i]*n/rank); q[i] = m
    return q
ids, pa, pb, info = [], [], [], []
for k, a in A.items():
    m = a['measure']; P = a['published']
    bpub = tr(m, P['estimate']); sepub = (tr(m, P['high']) - tr(m, P['low'])) / 3.919928
    br, ser = a['replication']['b'], a['replication']['se']; seo = a['original']['se']
    ids.append(k); pa.append(p2((br-bpub)/math.hypot(sepub, ser))); pb.append(p2((br-bpub)/math.hypot(max(sepub, seo), ser)))
    info.append((a.get('published_difference_p'), a.get('published_difference_q'), sepub, seo))
qa, qb = bh(pa), bh(pb)
print('id | bundle p, q | reviewer p, q (same method) | p, q with published SE floored at design SE | published SE | design SE')
for i, k in enumerate(ids):
    if qa[i] < 0.2 or qb[i] < 0.2 or (info[i][1] or 1) < 0.2:
        print(k, '|', round(info[i][0], 5), round(info[i][1], 4), '|', round(pa[i], 5), round(qa[i], 4), '|', round(pb[i], 5), round(qb[i], 4), '|', round(info[i][2], 4), round(info[i][3], 4))
print('significant (q<0.05): bundle', sum((x[1] or 1) < 0.05 for x in info), '| reviewer same method', sum(q < 0.05 for q in qa), '| design-SE floor', sum(q < 0.05 for q in qb))