# Balancing a crossover: simulated OSPAN in SD units, period correlation 0.7, practice +0.35 SD
# at the second session, true effect +0.20 SD, optional carry-over +0.10 SD to a desk session
# after a phone-free one. Estimate = mean within-person (B - A). Values illustrative.
import random, math
random.seed(5)
RHO, PRACTICE, EFFECT, N, REPS = 0.7, 0.35, 0.20, 160, 400
NOISE = math.sqrt((1 - RHO) / RHO)                # person SD 1 -> period correlation rho

def estimate(order, carry=0.0):
    diffs = []
    for i in range(N):
        person = random.gauss(0, 1)
        seq = order(i)                             # "AB" or "BA"
        score = {}
        for period, cond in enumerate(seq):
            y = person + random.gauss(0, NOISE) + PRACTICE * period
            y += EFFECT if cond == "B" else 0.0
            if period == 1 and seq == "BA":
                y += carry                         # desk session after a phone-free session
            score[cond] = y
        diffs.append(score["B"] - score["A"])
    return sum(diffs) / N

def mean_est(order, carry=0.0):
    return sum(estimate(order, carry) for _ in range(REPS)) / REPS

print(f"true effect {EFFECT:+.2f} SD; practice effect {PRACTICE:+.2f} SD at the second session")
print(f"  everyone AB (desk first)        : estimate {mean_est(lambda i: 'AB'):+.2f}")
print(f"  everyone BA (other room first)  : estimate {mean_est(lambda i: 'BA'):+.2f}")
bal = lambda i: "AB" if i % 2 else "BA"
print(f"  balanced AB/BA                  : estimate {mean_est(bal):+.2f}")
print(f"  balanced, with carry-over +0.10 : estimate "
      f"{mean_est(bal, carry=0.10):+.2f}")
print("balancing cancels the practice effect; carry-over still biases the estimate, so it must be")
print("prevented by design (wash-out of at least 7 days, parallel task forms)")
# true effect +0.20 SD; practice effect +0.35 SD at the second session
#   everyone AB (desk first)        : estimate +0.56
#   everyone BA (other room first)  : estimate -0.15
#   balanced AB/BA                  : estimate +0.20
#   balanced, with carry-over +0.10 : estimate +0.14
# balancing cancels the practice effect; carry-over still biases the estimate, so it must be
# prevented by design (wash-out of at least 7 days, parallel task forms)

# Sample size for the AB/BA crossover: A = own phone face down and silenced on the desk,
# B = phone locked in another room. Primary outcome OSPAN absolute score; effect in units of the
# between-person SD; SD of the within-person difference = SD * sqrt(2 (1 - rho)).
# Two-sided alpha 0.05, power 0.90, paired t approximation.
import math
from statistics import NormalDist
N = NormalDist()

def t_quantile(p, df):
    # Cornish-Fisher style approximation to the t quantile (good enough for df > 20)
    z = N.inv_cdf(p)
    return z + (z ** 3 + z) / (4 * df) + (5 * z ** 5 + 16 * z ** 3 + 3 * z) / (96 * df ** 2)

def n_pairs(d, rho, alpha=0.05, power=0.90):
    dz = d / math.sqrt(2 * (1 - rho))
    n = 10
    for _ in range(50):
        n_new = math.ceil(((t_quantile(1 - alpha / 2, n - 1) + t_quantile(power, n - 1)) / dz) ** 2)
        if n_new == n: break
        n = n_new
    return n

print("participants completing both periods, 90 % power")
print(" effect (SD) | rho = 0.6 | rho = 0.7 | rho = 0.8")
for d in (0.15, 0.20, 0.25, 0.30):
    print(f"    {d:.2f}     | {n_pairs(d, 0.6):9d} | {n_pairs(d, 0.7):9d} | {n_pairs(d, 0.8):9d}")
n = n_pairs(0.20, 0.7)
n_even = n + n % 2
recruit = math.ceil(n_even / 0.85); recruit += recruit % 2
print(f"chosen: 0.20 SD, rho = 0.7 -> {n} completers (rounded to {n_even} for equal sequences);")
print(f"allowing 15 % who miss the second session: randomise {recruit}"
      f" ({recruit // 2} per sequence)")
# participants completing both periods, 90 % power
#  effect (SD) | rho = 0.6 | rho = 0.7 | rho = 0.8
#     0.15     |       376 |       283 |       189
#     0.20     |       213 |       160 |       108
#     0.25     |       137 |       103 |        70
#     0.30     |        96 |        73 |        49
# chosen: 0.20 SD, rho = 0.7 -> 160 completers (rounded to 160 for equal sequences);
# allowing 15 % who miss the second session: randomise 190 (95 per sequence)

# Decision rule, fixed in advance. d = B (phone in another room) minus A (phone on desk) in OSPAN
# absolute score, in between-person SD units; positive = better without the phone. 95 % CI from
# the mixed model (condition, period, sequence; participant random). SESOI 0.20 SD. Carry-over
# is reported but never used to switch analyses.
SESOI = 0.20

def verdict(d, lo, hi, adherence_ok, harms):
    if harms:
        return f"d = {d:+.2f}: PAUSE and report - {harms} serious adverse event(s) judged related"
    if not adherence_ok:
        return (f"d = {d:+.2f}: NOT INTERPRETABLE - phone placement verified in fewer than"
                f" 90 % of sessions")
    if lo >= SESOI:    v = "phone absence improves OSPAN by at least the SESOI"
    elif lo > 0:       v = "phone absence improves OSPAN; size not shown to reach 0.20 SD"
    elif hi < 0:       v = "phone absence LOWERS OSPAN"
    elif lo >= -SESOI and hi <= SESOI: v = "no difference worth acting on (equivalence)"
    else:              v = "unclear: the interval is too wide to decide"
    return f"d = {d:+.2f} ({lo:+.2f} to {hi:+.2f}): {v}"

for case in ((0.34, 0.21, 0.47, True, 0), (0.15, 0.02, 0.28, True, 0),
             (0.03, -0.10, 0.16, True, 0), (-0.05, -0.27, 0.17, True, 0),
             (0.30, 0.17, 0.43, False, 0), (0.25, 0.12, 0.38, True, 1)):
    print(verdict(*case))
print("SART commission errors, RT variability and self-rated phone urges are secondary:")
print("reported with intervals, never decisive")
# d = +0.34 (+0.21 to +0.47): phone absence improves OSPAN by at least the SESOI
# d = +0.15 (+0.02 to +0.28): phone absence improves OSPAN; size not shown to reach 0.20 SD
# d = +0.03 (-0.10 to +0.16): no difference worth acting on (equivalence)
# d = -0.05 (-0.27 to +0.17): unclear: the interval is too wide to decide
# d = +0.30: NOT INTERPRETABLE - phone placement verified in fewer than 90 % of sessions
# d = +0.25: PAUSE and report - 1 serious adverse event(s) judged related
# SART commission errors, RT variability and self-rated phone urges are secondary:
# reported with intervals, never decisive