import math import random import statistics import unittest from statistics import NormalDist from power import plan, calculate, design_effect, rollout_probability, accounts_for_estimate_rule BASE = {'p0': .3, 'p1': .35, 'practical_threshold': .05, 'decision_rule': 'lower_bound_above_zero', 'family_alpha': .05, 'primary_hypotheses': 3, 'power': .8, 'users_per_account': 20, 'icc_values': [0, .05, .1, .2, .3], 'base_icc': .05, 'size_cv_values': [1.0], 'allocation': '1:1 accounts', 'attrition': 0, 'activation_completion_days': 7, 'reopen_observation_days': 7, 'followup_days': 14, 'operating_characteristic_gains': [0, .025, .05, .075], 'estimate_rule_planning_gains': [.06, .07]} class PowerTests(unittest.TestCase): def test_independent_reference(self): # R power.prop.test documentation's n=50, p1=.50, p2=.75 case # inverted using a normal-tail root, independently of plan's closed form. a, b, alpha, target = .5, .75, .05, .8 normal = NormalDist() def achieved(n): cutoff = normal.inv_cdf(1-alpha/2)*math.sqrt(2*((a+b)/2)*(1-(a+b)/2)/n) return normal.cdf(((b-a)-cutoff)/math.sqrt((a*(1-a)+b*(1-b))/n)) lo, hi = 1., 10000. for _ in range(100): mid = (lo+hi)/2 if achieved(mid) < target: lo = mid else: hi = mid self.assertAlmostEqual(plan(a,b,alpha,target,1,0)['independent_users_per_arm_unrounded'], hi, places=9) # The primary reference reports three decimal places, not extra precision. self.assertEqual(round(achieved(50), 3), .740) def test_known_z_reference(self): # External standard quantiles 1.95996398454 and 0.841621233573. expected = (1.95996398454*math.sqrt(.495)+.841621233573*math.sqrt(.49))**2/.01 self.assertAlmostEqual(plan(.5,.6,.05,.8,1,0)['independent_users_per_arm_unrounded'],expected,places=7) def test_design_effect_from_covariance(self): m, rho, p = 20, .05, .3 variance = p*(1-p) full_sum = sum(variance if i==j else rho*variance for i in range(m) for j in range(m)) measured = full_sum/(m*variance) self.assertAlmostEqual(plan(.3,.35,.05,.8,m,rho)['design_effect'],measured) def test_zero_icc(self): self.assertEqual(plan(.3,.35,.05,.8,20,0)['design_effect'],1) def test_singleton_cluster(self): self.assertEqual(plan(.3,.35,.05,.8,1,.9)['design_effect'],1) self.assertEqual(plan(.3,.35,.05,.8,1,.9)['total_accounts'],plan(.3,.35,.05,.8,1,0)['total_accounts']) def test_rounding_is_up_to_whole_accounts(self): for m in [1,3,20,101]: for rho in [0,.05,.5,1]: row=plan(.3,.35,.05,.8,m,rho) need=row['independent_users_per_arm_unrounded']*row['design_effect'] self.assertGreaterEqual(row['accounts_per_arm']*m,need) self.assertLess((row['accounts_per_arm']-1)*m,need) self.assertEqual(row['total_users'],row['total_accounts']*m) def test_increasing_icc_increases_sample(self): counts=[plan(.3,.35,.05,.8,20,r)['accounts_per_arm'] for r in [0,.05,.1,.2,.3]] self.assertEqual(counts,sorted(set(counts))) def test_tighter_alpha_and_more_power_cost_sample(self): baseline=plan(.3,.35,.05,.8,20,.05)['total_accounts'] self.assertGreater(plan(.3,.35,.05/3,.8,20,.05)['total_accounts'],baseline) self.assertGreater(plan(.3,.35,.05,.9,20,.05)['total_accounts'],baseline) def test_smaller_effect_costs_sample(self): self.assertGreater(plan(.3,.32,.05,.8,20,.05)['total_accounts'],plan(.3,.35,.05,.8,20,.05)['total_accounts']) def test_arm_swap_and_outcome_complement(self): base=plan(.3,.35,.05,.8,20,.05)['independent_users_per_arm_unrounded'] for a,b in [(.35,.3),(.7,.65)]: self.assertAlmostEqual(plan(a,b,.05,.8,20,.05)['independent_users_per_arm_unrounded'],base) def test_invalid_inputs_fail(self): base=[.3,.35,.05,.8,20,.05] for index,values in [(0,[0,1,float('nan'),True]),(1,[.3,-1,float('inf')]),(2,[0,1]),(3,[.5,1]),(4,[0,1.5,True]),(5,[-.1,1.1])]: for v in values: args=base.copy();args[index]=v with self.assertRaises(ValueError):plan(*args) for cv in [-.1,float('nan'),True]: with self.assertRaises(ValueError):plan(*base,cv) def test_unimplemented_designs_fail(self): a=dict(BASE) self.assertEqual(calculate(a)['rows'][0]['alpha'],.05/3) for key,value in [('allocation','2:1 accounts'),('attrition',.1),('primary_hypotheses',0),('icc_values',[0,0]), ('decision_rule','estimate_threshold'),('practical_threshold',.04),('followup_days',21), ('size_cv_values',[0]),('operating_characteristic_gains',[.025,.05]),('estimate_rule_planning_gains',[.05])]: with self.assertRaises(ValueError):calculate(a|{key:value}) def test_unequal_cluster_design_effect_from_covariance(self): # User-weighted mean over accounts of unequal size: enumerate the covariance of every user pair. sizes, rho, p = [4, 8, 14, 20, 26, 48], .05, .3 users = sum(sizes) mean = users/len(sizes) cv = math.sqrt(sum((s-mean)**2 for s in sizes)/len(sizes))/mean variance = p*(1-p) pair_sum = sum(sum(variance if i==j else rho*variance for i in range(s) for j in range(s)) for s in sizes) measured = pair_sum/(users*variance) self.assertEqual(mean, 20) self.assertAlmostEqual(design_effect(20, rho, cv), measured, places=12) self.assertEqual(design_effect(20, rho, 0), plan(.3,.35,.05,.8,20,rho)['design_effect']) effects=[design_effect(20,rho,c) for c in [0,.5,1,1.5]] self.assertEqual(effects,sorted(set(effects))) def test_sized_design_reaches_power_and_one_fewer_account_does_not(self): for rho in [0,.05,.1,.2,.3]: for cv in [0,1.0]: row=plan(.3,.35,.05/3,.8,20,rho,cv) self.assertGreaterEqual(rollout_probability(.3,.05,.05/3,row['accounts_per_arm'],20,rho,cv),.8) self.assertLess(rollout_probability(.3,.05,.05/3,row['accounts_per_arm']-1,20,rho,cv),.8) def test_operating_characteristics(self): result=calculate(dict(BASE)) oc={round(r['true_gain'],3):r for r in result['operating_characteristics']['rows']} self.assertEqual(result['operating_characteristics']['accounts_per_arm'],180) # No true gain: a rollout needs a false positive in the favourable tail only. self.assertAlmostEqual(oc[0]['lower_bound_rule'],.05/3/2,places=12) # The earlier rule (estimate must also reach 5 points) is a coin flip at a true 5-point gain. self.assertAlmostEqual(oc[.05]['earlier_estimate_threshold_rule'],.5,places=12) self.assertGreaterEqual(oc[.05]['lower_bound_rule'],.8) probabilities=[oc[g]['lower_bound_rule'] for g in sorted(oc)] self.assertEqual(probabilities,sorted(probabilities)) for row in oc.values(): self.assertLessEqual(row['earlier_estimate_threshold_rule'],row['lower_bound_rule']) def test_estimate_rule_sizing_matches_search(self): for gain in [.055,.06,.07,.09]: row=accounts_for_estimate_rule(.3,gain,.05,.05/3,.8,20,.05) k=1 while rollout_probability(.3,gain,.05/3,k,20,.05,estimate_threshold=.05)<.8:k+=1 self.assertEqual(row['accounts_per_arm'],k) with self.assertRaises(ValueError):accounts_for_estimate_rule(.3,.05,.05,.05/3,.8,20,.05) def test_simulated_decision_power(self): # Seed, replications and tolerance were fixed in the declaration before this test first ran. rng=random.Random(20261010) p0,p1,rho,m,k,reps,alpha=.3,.35,.05,20,180,2000,.05/3 z=NormalDist().inv_cdf(1-alpha/2) def arm(p): a=p*(1-rho)/rho;b=(1-p)*(1-rho)/rho return [rng.binomialvariate(m,rng.betavariate(a,b))/m for _ in range(k)] new=old=0 for _ in range(reps): control=arm(p0);treated=arm(p1) estimate=statistics.fmean(treated)-statistics.fmean(control) se=math.sqrt(statistics.variance(treated)/k+statistics.variance(control)/k) above=estimate-z*se>0 new+=above;old+=above and estimate>=.05 self.assertLessEqual(abs(new/reps-rollout_probability(p0,.05,alpha,k,m,rho)),.04) self.assertLessEqual(abs(old/reps-rollout_probability(p0,.05,alpha,k,m,rho,estimate_threshold=.05)),.04) if __name__=='__main__':unittest.main()