#!/usr/bin/env python3 """Plans and reads lifecycle holdouts: minimum detectable effect, duration, sample ratio mismatch and the interval for a lift. Usage: python3 scripts/holdout_math.py plan --baseline 0.21 --eligible-per-period 1200 --control-share 0.10 --periods 3 [--alpha 0.05] [--power 0.8] [--target-mde-pp 3] python3 scripts/holdout_math.py result --treat-n 3310 --treat-conv 745 --control-n 290 --control-conv 58 [--expected-control-share 0.10] [--alpha 0.05] plan: minimum detectable effect (MDE) in percentage points for the split and duration, using MDE = (z(1 - alpha/2) + z(power)) * sqrt(p * (1 - p) * (1/n_t + 1/n_c)). With --target-mde-pp it also prints the periods needed at this split and at 50/50. result: rates, absolute and relative difference, an unpooled 95% interval, and with --expected-control-share a sample ratio mismatch (SRM) chi-square test with 1 degree of freedom, p = erfc(sqrt(chi2 / 2)), FAIL below 0.0005. Small counts: when an arm has fewer than 5 conversions or 5 non-conversions, the unpooled interval is too narrow. result then prints NOT RELIABLE, a Newcombe score interval built from two Wilson intervals and a two-sided Fisher exact test, and reports not proven unless both exclude no difference. When an expected SRM count is below 5, it uses a two-sided exact binomial test, not chi-square. Assumes a binary outcome per unit, independent units and no spillover between arms. Clustered units, repeated outcomes and revenue need a different analysis. Python 3.9+, standard library only. Exit status: 1 on SRM FAIL, 2 on bad input, 0 otherwise. """ import argparse import math import sys from statistics import NormalDist SRM_THRESHOLD = 0.0005 SMALL_COUNT = 5 # below this many conversions or non-conversions in an arm, the normal interval is unreliable TIE = 1 + 1e-7 # tolerance when comparing exact probabilities ASSUMPTIONS = ("Assumptions: binary outcome per unit, independent units, no spillover between arms, two-sided test." " Clustered units, repeated outcomes and revenue need a different analysis.") def bad_input(message): print("holdout_math: " + message, file=sys.stderr) sys.exit(2) def check_share(name, value, low=0.0, high=1.0): if not low < value < high: bad_input("%s must be between %s and %s (exclusive)" % (name, low, high)) def mde(p, n_t, n_c, alpha, power): z = NormalDist().inv_cdf(1 - alpha / 2) + NormalDist().inv_cdf(power) return z * math.sqrt(p * (1 - p) * (1.0 / n_t + 1.0 / n_c)) def wilson(x, n, z): """Wilson score interval for one proportion.""" p = x / n centre = (p + z * z / (2 * n)) / (1 + z * z / n) half = z * math.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / (1 + z * z / n) return max(0.0, centre - half), min(1.0, centre + half) def newcombe(x_t, n_t, x_c, n_c, z): """Newcombe's hybrid score interval for p_t - p_c, built from two Wilson intervals.""" p_t, p_c = x_t / n_t, x_c / n_c l_t, u_t = wilson(x_t, n_t, z) l_c, u_c = wilson(x_c, n_c, z) diff = p_t - p_c return (diff - math.sqrt((p_t - l_t) ** 2 + (u_c - p_c) ** 2), diff + math.sqrt((u_t - p_t) ** 2 + (p_c - l_c) ** 2)) def log_comb(n, k): return math.lgamma(n + 1) - math.lgamma(k + 1) - math.lgamma(n - k + 1) def fisher_two_sided(x_t, n_t, x_c, n_c): """Two-sided Fisher exact test on the 2x2 table: sums every table with the same margins that is no more likely than the observed one.""" total = x_t + x_c low, high = max(0, total - n_t), min(total, n_c) base = log_comb(n_t + n_c, total) probs = [math.exp(log_comb(n_c, k) + log_comb(n_t, total - k) - base) for k in range(low, high + 1)] observed = probs[x_c - low] return min(1.0, sum(p for p in probs if p <= observed * TIE)) def binomial_two_sided(k, n, share): """Two-sided exact binomial test of k successes in n trials against share.""" log_p, log_q = math.log(share), math.log(1 - share) probs = [math.exp(log_comb(n, i) + i * log_p + (n - i) * log_q) for i in range(n + 1)] observed = probs[k] return min(1.0, sum(p for p in probs if p <= observed * TIE)) def periods_needed(p, per_period, share, target, alpha, power): """Periods at which the MDE first reaches the target (MDE falls with 1 / sqrt(periods)).""" one = mde(p, per_period * (1 - share), per_period * share, alpha, power) return (one / target) ** 2 def plan(args): check_share("--baseline", args.baseline) check_share("--control-share", args.control_share) check_share("--alpha", args.alpha) check_share("--power", args.power, 0.5, 1.0) if args.eligible_per_period <= 0 or args.periods <= 0: bad_input("--eligible-per-period and --periods must be above 0") total = args.eligible_per_period * args.periods n_c = total * args.control_share n_t = total - n_c effect = mde(args.baseline, n_t, n_c, args.alpha, args.power) print("holdout_math.py plan") print(ASSUMPTIONS) print("Inputs: baseline %.1f%%, %s eligible units per period, control share %s, %s periods, alpha %s," " power %s" % (args.baseline * 100, format(args.eligible_per_period, ",.0f"), "%g%%" % (args.control_share * 100), "%g" % args.periods, "%g" % args.alpha, "%.2f" % args.power)) print("Sample: %s units (treatment %s, control %s)" % ( format(total, ",.0f"), format(n_t, ",.0f"), format(n_c, ",.0f"))) print("Minimum detectable effect: %.2f percentage points (%.1f%% to %.1f%%, relative %.1f%%)" % ( effect * 100, args.baseline * 100, (args.baseline + effect) * 100, effect / args.baseline * 100)) expected_small = min(n_c * args.baseline, n_c * (1 - args.baseline)) if expected_small < SMALL_COUNT: print("Note: the control arm expects %.1f conversions or non-conversions, fewer than %d, so this normal" " approximation is rough and the readout will need the small-count method." % ( expected_small, SMALL_COUNT)) if args.target_mde_pp is not None: target = args.target_mde_pp / 100.0 if target <= 0: bad_input("--target-mde-pp must be above 0") here = periods_needed(args.baseline, args.eligible_per_period, args.control_share, target, args.alpha, args.power) even = periods_needed(args.baseline, args.eligible_per_period, 0.5, target, args.alpha, args.power) print("Target %.2f pp: needs %.1f periods at a %g%% control share (control %s units), or %.1f periods at" " 50/50 (%s units per arm)" % ( args.target_mde_pp, here, args.control_share * 100, format(here * args.eligible_per_period * args.control_share, ",.0f"), even, format(even * args.eligible_per_period * 0.5, ",.0f"))) if effect > target: print("WARN MH09 power: the minimum detectable effect (%.2f pp) is larger than the target (%.2f pp);" " this design cannot detect the smallest effect that would change the decision in %g periods." % (effect * 100, args.target_mde_pp, args.periods)) print("Summary: 0 FAIL, 1 WARN. Nothing was changed.") else: print("PASS MH09 power: the minimum detectable effect (%.2f pp) is within the target (%.2f pp)." % (effect * 100, args.target_mde_pp)) print("Summary: 0 FAIL, 0 WARN. Nothing was changed.") else: print("NOT CHECKED MH09: pass --target-mde-pp, the smallest effect that would change the decision.") print("Summary: 0 FAIL, 0 WARN. Nothing was changed.") return 0 def result(args): if args.treat_n <= 0 or args.control_n <= 0: bad_input("--treat-n and --control-n must be above 0") if not (0 <= args.treat_conv <= args.treat_n and 0 <= args.control_conv <= args.control_n): bad_input("conversions must be between 0 and the arm size") check_share("--alpha", args.alpha) p_t = args.treat_conv / args.treat_n p_c = args.control_conv / args.control_n diff = p_t - p_c se = math.sqrt(p_t * (1 - p_t) / args.treat_n + p_c * (1 - p_c) / args.control_n) z = NormalDist().inv_cdf(1 - args.alpha / 2) low, high = diff - z * se, diff + z * se level = (1 - args.alpha) * 100 print("holdout_math.py result") print(ASSUMPTIONS) print("Treatment: %s of %s = %.2f%%. Control: %s of %s = %.2f%%." % ( format(args.treat_conv, ","), format(args.treat_n, ","), p_t * 100, format(args.control_conv, ","), format(args.control_n, ","), p_c * 100)) relative = ("%+.1f%%" % (diff / p_c * 100)) if p_c > 0 else "undefined (control rate is 0)" print("Difference: %+.2f percentage points (relative lift %s)" % (diff * 100, relative)) print("%g%% interval for the difference (unpooled): %+.2f to %+.2f percentage points" % ( level, low * 100, high * 100)) smallest = min(args.treat_conv, args.treat_n - args.treat_conv, args.control_conv, args.control_n - args.control_conv) if smallest < SMALL_COUNT: s_low, s_high = newcombe(args.treat_conv, args.treat_n, args.control_conv, args.control_n, z) p_exact = fisher_two_sided(args.treat_conv, args.treat_n, args.control_conv, args.control_n) print("NOT RELIABLE: an arm has fewer than %d conversions or %d non-conversions, so the normal interval" " above is too narrow." % (SMALL_COUNT, SMALL_COUNT)) print("%g%% score interval (Newcombe, from two Wilson intervals): %+.2f to %+.2f percentage points." " Fisher exact test, two-sided: p %s." % ( level, s_low * 100, s_high * 100, "< 0.0001" if p_exact < 0.0001 else "= %.4f" % p_exact)) if p_exact < args.alpha and not s_low <= 0 <= s_high: print("The exact test and the score interval exclude 0. Check MH03 to MH08 before reporting it as" " incremental impact.") else: print("Not proven: the exact test or the score interval does not exclude 0.") elif low <= 0 <= high: print("Not proven: the interval includes 0.") else: print("The interval excludes 0. Check MH03 to MH08 before reporting it as incremental impact.") fails = 0 if args.expected_control_share is not None: check_share("--expected-control-share", args.expected_control_share) total = args.treat_n + args.control_n exp_c = total * args.expected_control_share exp_t = total - exp_c observed = args.control_n / total if min(exp_c, exp_t) < SMALL_COUNT: p_value = binomial_two_sided(args.control_n, total, args.expected_control_share) detail = ("observed control share %.2f%% against the intended %.2f%%; an expected count is %.1f, below" " %d, so a two-sided exact binomial test was used: p = %.6f" % ( observed * 100, args.expected_control_share * 100, min(exp_c, exp_t), SMALL_COUNT, p_value)) else: chi2 = (args.control_n - exp_c) ** 2 / exp_c + (args.treat_n - exp_t) ** 2 / exp_t p_value = math.erfc(math.sqrt(chi2 / 2)) detail = "observed control share %.2f%% against the intended %.2f%%; chi-square %.2f, p = %.6f" % ( observed * 100, args.expected_control_share * 100, chi2, p_value) if p_value < SRM_THRESHOLD: fails += 1 print("FAIL MH05 sample-ratio: %s (below %s). Find the cause (assignment, execution, log processing or" " analysis) before reading the lift." % (detail, SRM_THRESHOLD)) else: print("PASS MH05 sample-ratio: %s (not below %s)." % (detail, SRM_THRESHOLD)) else: print("NOT CHECKED MH05: pass --expected-control-share to test the split.") print("Summary: %d FAIL. Nothing was changed." % fails) return 1 if fails else 0 def main(): parser = argparse.ArgumentParser(description="Holdout planning and readout maths.") sub = parser.add_subparsers(dest="mode") p = sub.add_parser("plan") p.add_argument("--baseline", type=float, required=True, help="control conversion rate, for example 0.21") p.add_argument("--eligible-per-period", type=float, required=True) p.add_argument("--control-share", type=float, required=True, help="for example 0.10") p.add_argument("--periods", type=float, required=True) p.add_argument("--alpha", type=float, default=0.05) p.add_argument("--power", type=float, default=0.8) p.add_argument("--target-mde-pp", type=float, help="smallest effect that would change the decision, in pp") r = sub.add_parser("result") r.add_argument("--treat-n", type=int, required=True) r.add_argument("--treat-conv", type=int, required=True) r.add_argument("--control-n", type=int, required=True) r.add_argument("--control-conv", type=int, required=True) r.add_argument("--expected-control-share", type=float) r.add_argument("--alpha", type=float, default=0.05) args = parser.parse_args() if args.mode == "plan": sys.exit(plan(args)) if args.mode == "result": sys.exit(result(args)) parser.print_help() sys.exit(2) if __name__ == "__main__": main()