import math
from statistics import NormalDist

N = NormalDist()

def n_per_group(sd, delta, rho, alpha=0.05, power=0.80):
    """Two parallel groups, outcome analysed with baseline as covariate (ANCOVA): the residual SD is
    sd * sqrt(1 - rho^2), where rho is the baseline-follow-up correlation. Two-sided alpha."""
    s = sd * math.sqrt(1 - rho ** 2)
    za, zb = N.inv_cdf(1 - alpha / 2), N.inv_cdf(power)
    return math.ceil(2 * ((za + zb) * s / delta) ** 2 + za ** 2 / 4)

# Primary outcome: 7-day mean of waking lnRMSSD (natural log of ms), spontaneous breathing, at week 8.
# Between-person SD about 0.5 ln units; smallest worthwhile difference 0.15 (about 16 % in RMSSD).
for rho in (0.6, 0.7, 0.8):
    print(f"baseline correlation {rho}: {n_per_group(0.5, 0.15, rho)} per group for 80 % power, "
          f"{n_per_group(0.5, 0.15, rho, power=0.90)} for 90 %")
n = n_per_group(0.5, 0.15, 0.7, power=0.90)
print(f"plan: {n} completers per group ({2 * n} in all); recruit {math.ceil(2 * n / 0.8)} allowing 20 % loss")
# Secondary outcome, CAR as AUCi (nmol/L x h) from samples at waking, +30 and +45 min on 2 days.
print(f"CAR AUCi, SD 2.5, difference 1.0, rho 0.4: {n_per_group(2.5, 1.0, 0.4)} per group for 80 % power")
# baseline correlation 0.6: 113 per group for 80 % power, 151 for 90 %
# baseline correlation 0.7: 90 per group for 80 % power, 121 for 90 %
# baseline correlation 0.8: 64 per group for 80 % power, 86 for 90 %
# plan: 121 completers per group (242 in all); recruit 303 allowing 20 % loss
# CAR AUCi, SD 2.5, difference 1.0, rho 0.4: 84 per group for 80 % power

SESOI = 0.15     # ln units of waking RMSSD: the smallest difference worth acting on (about 16 %)

def verdict(diff, lo, hi):
    """diff = intervention minus usual routine at week 8 (ANCOVA, adjusted for baseline), with its 95 %
    confidence interval. Fixed before any data are seen."""
    if lo >= SESOI:
        v = "increase of at least the smallest worthwhile size"
    elif lo > 0:
        v = "increase; not shown to reach the smallest worthwhile size"
    elif hi < 0:
        v = "decrease compared with usual routine"
    elif lo >= -SESOI and hi <= SESOI:
        v = "no difference worth acting on"
    else:
        v = "unclear"
    return f"lnRMSSD {diff:+.2f} ({lo:+.2f} to {hi:+.2f}): {v}"

# Invented results, to show how the trial will be reported.
print(verdict(0.06, -0.04, 0.16))
print(verdict(0.21, 0.16, 0.26))
print(verdict(0.02, -0.08, 0.12))
print("CAR (AUCi), perceived stress and sleep are reported with intervals but cannot make the trial")
print("positive. A rise or a fall in CAR is reported as found: the direction is not fixed in advance.")
# lnRMSSD +0.06 (-0.04 to +0.16): unclear
# lnRMSSD +0.21 (+0.16 to +0.26): increase of at least the smallest worthwhile size
# lnRMSSD +0.02 (-0.08 to +0.12): no difference worth acting on
# CAR (AUCi), perceived stress and sleep are reported with intervals but cannot make the trial
# positive. A rise or a fall in CAR is reported as found: the direction is not fixed in advance.