abfactor

Bayesian A/B testing with Bayes factors, built on the aggregated data experimentation platforms already produce, and checked by simulation.

• Bayes factors, not just P(B > A). Evidence for and against "no effect", directional tests, and exact pairwise comparisons in multi-arm tests.
• Default priors that suit online experiments. The fixed-effects construction of Rouder et al. (2012), with the prior scale on the model scale (log odds, log rate) so tiny real effects aren't swamped by Lindley's paradox.
• Experiments that learn from each other. A shared effect scale and a spike-and-slab across a batch, with each experiment's Bayes factor computed given the rest.
• Aggregated data done properly. Per-user squares and cross products give dispersion for counts, real variances for per-user means and delta-method SEs for ratios.
• Time and pre-periods. One-off conversions as a discrete-time hazard on the risk set, A/A periods that catch broken buckets, CUPED and difference-in-differences.
• Sequential monitoring that keeps its false-positive guarantee.

Install

pip install -e ".[dev]"
pytest              # fast checks, ~25 s
pytest -m slow      # also checks the quadrature engine against NUTS and SMC

How it fits together

An analysis is four choices (see docs/design.md):

Data come in as Stats: sufficient statistics with dimensions experiment × arm × period (n, x, x__sq, x__x__y: plain column names any warehouse accepts), built from user rows, event times, a long frame of query output, or arrays.

Quick start

from abfactor import Analysis, Stats, Proportion, DefaultFixedEffects

stats = Stats.from_arrays(arms=["control", "b", "c"], n=[[20000, 20000, 20000]],
                          converted=[[1061, 1102, 1131]])
result = Analysis(Proportion("converted"), baseline=0.05, prior=DefaultFixedEffects(r=0.1)).fit(stats)

result.omnibus      # log BF10: does any arm differ?
result.pairwise     # each variant vs control: BF10, BF+0, P(better), contrast
result.summary()    # plus the relative lift with a 90% interval

From user-level data, with a ratio metric:

from abfactor import Ratio, NormalEffects

stats = Stats.from_users(users, metrics={"revenue": "revenue", "orders": "boxes"}, control="control")
aov = Analysis(Ratio("revenue", "orders"), baseline=50, prior=NormalEffects(5)).fit(stats)

A batch of experiments that learn from each other:

from abfactor import SpikeAndSlab

result = Analysis(Proportion("converted"), 0.05, prior=DefaultFixedEffects(0.5),
                  pooling=SpikeAndSlab()).fit(batch)
result.omnibus["P_effect"]       # each experiment's probability of having any effect
result.pooling["w_mean"]         # share of experiments with effects

One-off conversions with an A/A period before the treatment, guarding against broken buckets:

from abfactor import TimeToEvent

stats = Stats.from_events(users, event_period="converted_period", n_periods=2, control="control")
result = Analysis(TimeToEvent("conversion"), 0.05, treatment_start=1, offsets="free").fit(stats)

Straight from a warehouse query, one row per experiment and arm:

SELECT experiment_id, variation,
       COUNT(*)               AS n,
       SUM(revenue)           AS revenue,  SUM(revenue * revenue) AS revenue__sq,
       SUM(orders)            AS orders,   SUM(orders * orders)   AS orders__sq,
       SUM(revenue * orders)  AS revenue__x__orders
FROM user_level GROUP BY 1, 2

stats = Stats.from_frame(query_result, control="control", experiment="experiment_id", arm="variation")

Every experiment in a Stats has the same arms. When variants differ between experiments (two arms here, three there, or different names), Stats.batches_from_frame returns one Stats per set of arms; Analysis.fit takes the list, pools across all of it, and returns one Result per part.

Sequential monitoring, pooling across completed experiments only:

from abfactor import Monitor

prior = Analysis(Proportion("converted"), 0.05, prior=DefaultFixedEffects(0.5)).learn_prior(completed)
monitor = Monitor(Analysis(Proportion("converted"), 0.05, prior=prior), threshold=20)
monitor.update(stats_today)      # cumulative counts; stops the first time the evidence reaches 20

For conversion rates the monitor uses an exact e-process (method="exact"): the false-positive rate stays below 1/threshold however often you look and whatever the baseline. Other metrics use the Bayes factor (method="bayes"), checked by simulation. See docs/sequential.md.

Noise between identical arms

In real experiments identical arms can differ by more than sampling noise allows. In the Upworthy headline archive the excess was about 13% of the click-through rate, and ignoring it made one A/A test in 16–40 show strong evidence of an effect; in ASOS's e-commerce tests it was at most 0.6% (docs/real_data.md). Measure it from your own A/A tests:

from abfactor.aa import estimate_arm_sd, aa_check

est = estimate_arm_sd(Proportion("converted"), aa_stats)     # SD on the model scale, with a 95% interval
analysis = Analysis(Proportion("converted"), 0.05, arm_sd=est["arm_sd"])
aa_check(analysis, aa_stats)                                 # share of A/A tests with BF10 > 3, > 20, < 1/3

Many variants and metrics

abfactor.multiple adjusts for asking many questions at once (docs/multiple.md):

from abfactor.multiple import adjust_direction, joint, ebh

claims = adjust_direction(results, p0=0.5)   # sign claims with prior odds adjusted for the family
joint(results, claims=claims)                # P(all those claims hold), from the draws
ebh(monitor.e_values().E, alpha=0.05)        # false discovery rate with exact e-values

Several metrics on the same users (conversion, revenue, AOV) can be analysed jointly, so that "revenue improves and conversion doesn't drop" gets one calibrated probability: MultiMetric({...}).fit(stats).decide("revenue", guardrails={"conversion": 0.01}).

Saving and reproducing results

result.to_netcdf("checkout_test.nc")           # draws, Bayes factors, settings, package versions
result = Result.from_netcdf("checkout_test.nc")
result.to_frame()                               # flat table: one row per experiment and variant

Analysis.from_settings(result.settings).fit(stats)   # recomputes the same draws exactly

Small samples

The default engine approximates each metric's likelihood by a normal on its model scale. Measured against exact integration (docs/accuracy.md), it never overstated the evidence for an effect, but:

• with fewer than 5 conversions (or non-conversions) in an arm, intervals can be far too narrow: fit raises a SmallSampleWarning, and engine="smc" or "nuts" should be used;
• with 5 to 49, Bayes factors can understate the evidence by up to a factor of about 2.3 near the decision thresholds: the experiment is listed in result.diagnostics["small_sample"]["conservative_bf"], so borderline calls can be confirmed with engine="smc".

Choosing the prior scale

Bayes factors depend on how big the alternative expects effects to be. A 10% lift on a 5% conversion rate is a log odds ratio of about 0.1. With r = 0.5 (Rouder et al.'s standardised default) the alternative bets on effects around 0.7 and the null wins almost regardless of the data. Use r around 0.1 for conversion rates, or pool across experiments and let the batch learn the scale. Real effect sizes vary a lot: pairwise effects were about 0.02 on the log scale in ASOS's e-commerce tests and 0.4 in Upworthy's headline tests (docs/real_data.md), so learn_prior on your own completed experiments beats any default.

What the simulations show

From python harness/harness.py --quick (60 experiments per scenario; full results in harness/results/):

• Headline: the recommended analyses' calibration ratio (how big the actual error is, as a multiple of the stated uncertainty; 1 is ideal) is 0.87–1.02, never overconfident beyond noise. The five deliberately misspecified baselines score 1.11–2.12, every one overconfident (scorecard).
• Sequential monitoring: stopping at BF10 ≥ 20 gave 0.8% false positives over 30 looks (the guarantee is 5%), against 33% for peeking at a z-test.
• Overdispersed counts: a Poisson model gave false positives of 22% in one tail; the quasi-Poisson correction and the per-user mean brought it back to 3–12%.
• AOV: treating orders as independent gave intervals that held the true lift 62–68% of the time; the delta method 85–92%.
• Broken buckets: with an A/A period, modelling arm offsets cut false positives from 62% to 2–5%.
• CUPED: at a pre/post correlation of 0.62 it cut the SD of the lift by 20% (theory: 21%), with level intervals now calibrated too (89–91% coverage at 200 experiments).

To see what each analysis assumed and how far the data moved it, add --report harness_report.html: for a couple of experiments per scenario and model it plots the prior and the posterior against the truth and the raw estimates, arm against arm on the natural scale and each contrast on the model scale (harness/report.py). The same draws are available for any analysis from Analysis.sample_prior(stats).

Known limitations

What 1.0 deliberately does not do is listed in docs/api.md: experiments that share users, learned arm offsets with the default engine, the exact sequential guarantee for metrics other than conversion rates, and difference in differences for ratios. Beyond those:

• Bayes factors from the default engine can understate the evidence when an arm has fewer than about 50 events (never overstate it); with fewer than 5, use engine="smc" or "nuts" (docs/accuracy.md). Posterior intervals are calibrated down to 5 events per arm.
• Real-data checks cover conversion rates (Upworthy, with A/A arms) and counts and revenue (ASOS, without A/A arms); see docs/real_data.md. Arm noise and typical effect sizes differed between the two by more than an order of magnitude, so measure both on your own experiments.

Also included

• abfactor.brease: the BREASE framework (Irons & Cinelli 2023) in full; also usable as prior=Brease(...) for a two-arm binary experiment.
• abfactor.anova: Rouder et al.'s standardised-effects ANOVA Bayes factors, for reference.
• abfactor.baselines: independent-beta and abtest-style logit-transform Bayes factors.
• examples/: a tour of the API (quickstart.py) and reproductions of the BREASE and abtest papers.
• harness/sbc.py: simulation-based calibration of every engine and metric (docs/calibration.md); harness/upworthy.py: the real-data checks.
• docs/methods/: the methods explained step by step, from the maths to the abfactor code: Bayes factors, metrics beyond conversion and many experiments at once.

Licence, stability and roadmap

MIT licensed (LICENSE). What counts as public API, and what 1.0 will promise about it, is in docs/api.md. The documentation site, with these guides, the methods, the harness scorecards and an API reference built from the docstrings, builds with pip install -e ".[docs]" then mkdocs serve (preview) or mkdocs build. The remaining work before 1.0 is in docs/roadmap.md.

References

• Rouder, Morey, Speckman & Province (2012). Default Bayes factors for ANOVA designs. Journal of Mathematical Psychology 56(5).
• Rouder, Speckman, Sun, Morey & Iverson (2009). Bayesian t tests for accepting and rejecting the null hypothesis. Psychonomic Bulletin & Review 16(2).
• Irons & Cinelli (2023). Causally sound priors for binary experiments. arXiv.
• Gronau, Raj K. N. & Wagenmakers (2021). Informed Bayesian inference for the A/B test. Journal of Statistical Software 100(17).
• Deng, Knoblich & Lu (2018). Applying the delta method in metric analytics. KDD.
• Wedderburn (1974). Quasi-likelihood functions. Biometrika 61(3).
• Dickey & Lientz (1970); Wagenmakers et al. (2010): the Savage–Dickey density ratio.

Full bibliography: docs/references.bib.