## Code by Georgi Z. Georgiev (2023). Developed with the help of Ronny Kohavi & Alex Deng ## Simulation with a true effects distribution mimicking Bing's data from Figure 3a of Azevedo E.M., Deng A., Olea J.L.M., Rao J., and Weyl E.G. (2019) "A/B Testing with Fat Tails", Journal of Political Economy (128), 12. doi: 10.1086/710607 ## Key parameters follow the paper: 20 million users per test group in a one-sided test at 0.025 (statsig is z-score > 1.96), equal to a two-sided test at 0.05, and the base rate is 0.25. These can be adjusted at will to explore other scenarios. ## The adjusted FPR formula works since it uses the empirical p-values instead of alpha (sig_threshold). The FPR computed using equation (4) does not work, as expected, even if the target power of 0.8 is replaced with the mean power against the true effect sizes library(distr) ## Basic simulation parameters sim_runs = 100000; sample_size = 20000000; #sample size per test group base_rate = 0.25; # baseline conversion rate sig_threshold = 0.025; # 0.025 one-sided = 0.05 two-sided. ## Define the distribution of true effects so it has heavy tails. Not necessary and can be replaced with a simple normal distribution, but produces true effects more aligned with those observed in practice. myMix = UnivarMixingDistribution(Norm(mean=0, sd=0.0002), Norm(mean=0, sd=0.0007), mixCoeff=c(0.8, 0.2)); rmyMix <- r(myMix); ## Draw random true effects from the distribution true_effects = rmyMix(sim_runs); #true_effects <- rep(0,sim_runs) #uncomment for nil effect sims (a.k.a. A/A tests, should result in FPR = 1.0, Fr(SS|H0) = sig_threshold) ## Assign random outcomes for the test groups obs_control = rbinom(n=sim_runs, size=sample_size, prob=base_rate); obs_variant = rbinom(n=sim_runs, size=sample_size, prob=base_rate + true_effects); ## Calculate power = 1-beta against the true (otherwise unknown) effect size, replaces P(SS|H1) in the FPR equation power = power.prop.test(n = sample_size, p1 = base_rate, p2 = base_rate + true_effects[true_effects > 0], sig.level = sig_threshold, alternative='one.sided')\$power ## p-value for absolute difference of proportions calcPval = function(i) { return(prop.test(x = c(obs_variant[i], obs_control[i]), n = c(sample_size, sample_size), alternative='greater')\$p.value); } pvals = sapply(1:sim_runs, calcPval); rejected = pvals < sig_threshold; # statistically significant outcome rejecting H0, or not? falsepos = sum(rejected[true_effects <= 0]); # count of false positive tests (statistically significant with true H0) truepos = sum(rejected[true_effects > 0]); # count tests as true positive (statistically significant with false H0) fpf = falsepos / (falsepos + truepos); # frequency of false positives among all positives trueH0 = length(true_effects[true_effects <= 0]); # count of true H0s FrSSH0 = falsepos / trueH0; # Fr(ss|H0) # frequency of rejected H0s mean_power = mean(power, na.rm = T); # mean power against the true effect size determined analytically #mean_power = truepos / (sim_runs - trueH0) # empirical mean power against the true effect size (can be used to check pre-test vs empirical, but pre-test is the required input as per the FPR formula) ## False positive risk equation calculation true_odds = 1/1 # with a median of zero of the true effects distribution and simple superiority test the true odds are 1/1 fpr = FrSSH0 / mean_power * true_odds / (1 + FrSSH0 / mean_power * true_odds) fpreq4 = sig_threshold / 0.8 * true_odds / (1 + sig_threshold / 0.8 * true_odds) # FPR computed following equation (4) ## Output of FPR and other statistics and sanity checks cat(paste('False positive freq. (empirical): ', fpf, '(', round(fpf * 100, 2), '%)', '\nFalse positive risk (formula) : ', fpr, '(', round(fpr * 100, 2), '%)', '\nFalse positive risk (equation 4): ', fpreq4, '(', round(fpreq4 * 100, 2), '%)', '\nMean true effect: ', mean(true_effects), '\nMedian true effect: ', median(true_effects), '(', round(median(true_effects),6), ')', '\nTrue effect sigma: ', sd(true_effects), '\nFr(true_effect <= 0) = ', trueH0 / length(true_effects), '\nFr(ss|H0): ', FrSSH0, '(', round(FrSSH0 * 100, 2), '%)', '\nMean power vs true effects under H1 (Fr(ss|H1)): ', mean_power, '(', round(mean_power * 100, 2), '%)'), sep='\n') ## Examine the input for the FPR equation (comment out the unnecessary ones) and sanity checks hist(true_effects, breaks=40); hist(pvals, breaks=20); #hist(power, breaks=20);