# Sample size for the Chapter 8 trial. Unit of randomisation and analysis: the
# 5-person team. Each team does m = 10 estimation tasks with known answers.
# Per task, the anchoring index c = |G - R| / (|G - F| + |G - R|) on log
# estimates (G agreed group estimate, F mean private estimate of the first two
# speakers, R mean of the other three); team score = mean c over its tasks.
# Assumptions (to be confirmed): SD of c across tasks sigma = 0.25; share of
# that variance lying between teams (ICC) rho = 0.30; smallest difference of
# interest 0.10; two-sided alpha 0.05; power 0.90; normal approximation + 1.
import math, random
from statistics import NormalDist
Z = NormalDist().inv_cdf(0.975) + NormalDist().inv_cdf(0.90)
SIGMA, M = 0.25, 10

def teams_per_arm(delta, rho, m=M, sigma=SIGMA):
    var_team = sigma ** 2 * (rho + (1 - rho) / m)      # variance of a team mean
    return math.ceil(2 * Z ** 2 * var_team / delta ** 2) + 1

print("teams per arm for 90 % power (10 tasks per team)")
print("  difference | ICC rho = 0.2 | 0.3 | 0.5")
for delta in (0.08, 0.10, 0.12):
    print(f"     {delta:.2f}    | {teams_per_arm(delta, 0.2):13d} |"
          f" {teams_per_arm(delta, 0.3):3d} | {teams_per_arm(delta, 0.5):3d}")
n = teams_per_arm(0.10, 0.30); rec = math.ceil(n / 0.85)
print(f"chosen: 0.10, rho = 0.30 -> {n} analysable teams per arm;"
      " allowing 15 % loss")
print(f"  recruit {rec} teams per arm: {2 * rec} teams, {2 * rec * 5} people")
naive = math.ceil(2 * Z ** 2 * SIGMA ** 2 / 0.10 ** 2 / M) + 1
print(f"ignoring clustering (tasks as independent) would suggest only {naive}"
      " teams per arm;")
print(f"  design effect 1 + (m - 1) rho = {1 + (M - 1) * 0.30:.2f}")
rng = random.Random(15)
def sim_power(teams, delta, rho, sims=2000):
    hits = 0
    for _ in range(sims):
        means = []
        for arm_shift in (0.0, delta):
            sd_b, sd_w = SIGMA * math.sqrt(rho), SIGMA * math.sqrt(1 - rho)
            means.append([arm_shift + rng.gauss(0, sd_b) +
                          sum(rng.gauss(0, sd_w) for _ in range(M)) / M
                          for _ in range(teams)])
        a, b = means
        ma, mb = sum(a) / teams, sum(b) / teams
        va = sum((x - ma) ** 2 for x in a) / (teams - 1)
        vb = sum((x - mb) ** 2 for x in b) / (teams - 1)
        hits += abs(mb - ma) / math.sqrt((va + vb) / teams) > 1.96
    return hits / sims
p1, p2 = sim_power(n, 0.10, 0.30), sim_power(naive, 0.10, 0.30)
print(f"simulation check: power with {n} teams per arm = {p1:.2f};"
      f" with {naive} = {p2:.2f}")
# teams per arm for 90 % power (10 tasks per team)
#   difference | ICC rho = 0.2 | 0.3 | 0.5
#      0.08    |            59 |  77 | 114
#      0.10    |            38 |  50 |  74
#      0.12    |            27 |  35 |  52
# chosen: 0.10, rho = 0.30 -> 50 analysable teams per arm; allowing 15 % loss
#   recruit 59 teams per arm: 118 teams, 590 people
# ignoring clustering (tasks as independent) would suggest only 15 teams per arm;
#   design effect 1 + (m - 1) rho = 3.70
# simulation check: power with 50 teams per arm = 0.91; with 15 = 0.46

# Decision rule for the Chapter 8 trial, fixed before any data are collected.
# diff = arm A (private written judgments first, anonymous summary, then the
# usual discussion) minus arm B (usual discussion only) in the team-level
# anchoring index; NEGATIVE diff = less anchoring on the first two speakers.
# 95 % CI from the pre-specified linear model (arm + organisation stratum).
# Smallest difference of interest 0.10 on the 0-1 index.
SESOI = 0.10

def verdict(diff, lo, hi, fidelity, contamination):
    if fidelity < 0.80:
        return (f"diff {diff:+.2f}: NOT INTERPRETABLE - private judgments first"
                f" in only {fidelity:.0%} of A sessions")
    note = ("\n    caution: arm B teams used written rounds in"
            f" {contamination:.0%} of sessions") if contamination > 0.20 else ""
    if hi <= -SESOI:   v = "A reduces anchoring by at least 0.10"
    elif hi < 0:       v = "A reduces anchoring; not shown to reach 0.10"
    elif lo > 0:       v = "A INCREASES anchoring"
    elif lo >= -SESOI and hi <= SESOI: v = "no difference worth acting on"
    else:              v = "unclear: the interval is too wide to decide"
    return f"diff {diff:+.2f} ({lo:+.2f} to {hi:+.2f}): {v}{note}"

for case in ((-0.18, -0.25, -0.11, 0.95, 0.05), (-0.07, -0.14, 0.00, 0.92, 0.10),
             (-0.12, -0.19, -0.05, 0.90, 0.05), (0.01, -0.06, 0.08, 0.93, 0.08),
             (-0.20, -0.27, -0.13, 0.70, 0.05), (-0.15, -0.22, -0.08, 0.90, 0.30)):
    print(verdict(*case))
print("group error |log agreed estimate - log truth| is secondary: reported with")
print("  intervals, never used to override the primary verdict")
# diff -0.18 (-0.25 to -0.11): A reduces anchoring by at least 0.10
# diff -0.07 (-0.14 to +0.00): unclear: the interval is too wide to decide
# diff -0.12 (-0.19 to -0.05): A reduces anchoring; not shown to reach 0.10
# diff +0.01 (-0.06 to +0.08): no difference worth acting on
# diff -0.20: NOT INTERPRETABLE - private judgments first in only 70% of A sessions
# diff -0.15 (-0.22 to -0.08): A reduces anchoring; not shown to reach 0.10
#     caution: arm B teams used written rounds in 30% of sessions
# group error |log agreed estimate - log truth| is secondary: reported with
#   intervals, never used to override the primary verdict