# Reviewer's sensitivity: the bundle's estimator (log RR, DerSimonian-Laird tau^2, modified
# Hartung-Knapp interval, t with k-1 df), re-implemented from the paper's Methods, on its three
# 28-day contrasts, then with IVERCOR-COVID19 (Vallejos et al. 2021: hospitalization 14/250 vs 21/251) added.
import math
from scipy.stats import t as tdist
def pool(rows):
    y = [math.log((a/n)/(c/m)) for a, n, c, m in rows]
    v = [1/a - 1/n + 1/c - 1/m for a, n, c, m in rows]
    w = [1/x for x in v]; k = len(rows)
    yb = sum(wi*yi for wi, yi in zip(w, y)) / sum(w)
    Q = sum(wi*(yi-yb)**2 for wi, yi in zip(w, y))
    tau2 = max(0, (Q-k+1) / (sum(w) - sum(wi*wi for wi in w)/sum(w)))
    ws = [1/(vi+tau2) for vi in v]
    mu = sum(wi*yi for wi, yi in zip(ws, y)) / sum(ws)
    se = math.sqrt(max(1, sum(wi*(yi-mu)**2 for wi, yi in zip(ws, y))/(k-1)) / sum(ws))
    q = tdist.ppf(0.975, k-1)
    fe = math.sqrt(1/sum(w))
    return dict(k=k, rr=math.exp(mu), lo=math.exp(mu-q*se), hi=math.exp(mu+q*se), tau2=tau2,
                fixed=math.exp(yb), flo=math.exp(yb-1.96*fe), fhi=math.exp(yb+1.96*fe))
main = [(34, 2157, 27, 1806), (10, 817, 9, 774), (7, 708, 2, 722)]
fmt = lambda d: {k: round(v, 3) if isinstance(v, float) else v for k, v in d.items()}
print("bundle's three contrasts:", fmt(pool(main)))
print("plus IVERCOR-COVID19:   ", fmt(pool(main + [(14, 250, 21, 251)])))