Online A/B testing is causal inference at industrial scale: Microsoft, Netflix, and DoorDash each run tens of thousands of experiments per year, and the methods they have invented to handle that volume — CUPED, anytime-valid sequential tests, switchback designs, and Bayesian decision rules — are now the field’s modern toolkit. This session shows that S04’s blocking is the same thing as tech-industry stratified A/B testing, introduces CUPED and its ML extensions (CUPAC, MLRATE, CURE) as Lin-(2013) adjustment with the pre-period outcome, and then resolves the two open promises from S04: Fisher’s permutation logic becomes A/B-test p-values, and the Bayesian credible interval becomes the proper cure for the frequentist CI misreadings. The session closes with what breaks when textbook A/B assumptions fail — sample ratio mismatch, SUTVA violations in marketplaces, and the switchback designs that DoorDash and Lyft use to handle them — and a surrogate-index slide showing the limits of “long-horizon decisions on short-horizon proxies.” The PBL session works through five progressive exercises spanning all five stacks; Exercise 5 (SRM diagnostic) directly prepares the Quiz II material that gates Assignment II.
Bing once attributed a $100M per year annual revenue lift to a single ad-layout experiment that a junior engineer had built on a slow afternoon (Kohavi, Tang, and Xu 2020, ch. 1). Across Microsoft, Booking, Netflix, LinkedIn, and Airbnb the modal experiment outcome looks nothing like that: Kohavi’s “1/3, 1/3, 1/3” rule of thumb says that roughly one third of features ship with a positive effect, one third are flat, and one third hurt the OEC and must be reverted. Stated more bluntly: about 85% of product-team hypotheses fail to ship as proposed. That asymmetry — the rare blockbuster, the routine bust — is precisely why platform-scale experimentation pays off. The expected value of a single feature is small and noisy; the expected value of a pipeline that screens hundreds of features per quarter through a trustworthy A/B platform is enormous.
The implication for this course is operational: the techniques in S04 (random assignment, OLS with HC-robust SEs, Lin adjustment) are necessary but not sufficient. Once experiments run continuously, in parallel, on millions of users, new failure modes appear that textbook RCTs never had to confront — and new tools exist to handle them. 4 put it directly: online controlled experimentation creates statistical challenges that are not encountered in traditional applications, and the unprecedented scale at which experimentation is now conducted in industry has driven much of the recent innovation in this area. The rest of this session catalogues those tools.
Every A/B platform forces one decision before the first user is randomized: which single number defines winning? Kohavi calls this the Overall Evaluation Criterion (OEC) (Kohavi, Tang, and Xu 2020, ch. 7). A good OEC is
Around the OEC sit two supporting tiers. Driver metrics are the secondary numbers an analyst inspects when the OEC moves (clicks per session, time on page, scroll depth). Guardrail metrics are the constraint metrics that block a launch even when the OEC improves: page-load time, crash rate, error rate, complaint volume. Data-quality metrics are the diagnostic layer — sample-ratio mismatch percentage, missing-data fraction, instrumentation-error counts — that decide whether the experiment is even interpretable.
The cleanest example from Kohavi ch. 7 makes the OEC-design point viscerally: a search engine could plausibly optimise for queries per user, on the logic that more queries = more user engagement = more ad inventory. But a user typing five queries in a row before clicking a result is failing to find what they want. The metric that aligns with the long-term business is sessions per user — the user came back tomorrow because today’s session resolved their need. Optimising queries-per-user can directly destroy sessions-per-user; the wrong OEC choice ships features that look good on the launch dashboard and erode user retention over the following months. Once you have lived through one of those mistakes, you stop arguing about the OEC and start fighting for it.
The OEC framework would not exist if not for a corpus of painful experiences. 15 documented five real Bing experiments whose headline results were ultimately driven not by the feature change itself but by operational artefacts: click-tracking failures inflating treatment metrics, primacy effects that decayed over the first week, carryover from a prior experiment running on the same users, and an analysis that ignored the experiment’s actual sample-ratio mismatch. The lesson is uncomfortable: at platform scale, most “interesting” results are wrong before they are interesting. The five puzzles motivate the diagnostic discipline of when A/B breaks (Twyman’s law, SRM checks, SUTVA, surrogate-index caveats); they also explain why modern A/B platforms ship with extensive automated guardrails rather than treating the analyst’s judgment as the last line of defence.
A direct continuation lives in 17, who replicate the 55% click-through-rate lift reported in 16 Study 2 (rounded vs. sharp CTA buttons on an e-commerce homepage; \(n = 919\) visits, \(p = 0.037\)). Three pre-registered replications at SeaWorld and Coop Norway — with 2.8M, 2.2M, and 1.9M users (each over 2,000× the original sample) — returned estimated lifts of 0.16% (\(p = 0.20\)), 0.29% (\(p = 0.60\)), 0.73% (\(p = 0.09\)): two orders of magnitude smaller, all 95% confidence intervals containing zero. BAC’s Study 2 was simply under-powered: no pre-declared MDE, no power calculation, and a 919-visit sample in a domain where real effects sit well under 1%. (Their companion Study 1 had separate problems — a 44:56 sample-ratio mismatch with \(z = 46\), plus Google Ads’ non-random assignment — but the +55% headline came from Study 2, not Study 1.) The winner’s-curse mechanism explains the 55% perfectly: any underpowered study that produces a “significant” result must have over-estimated the true effect (the rejected hypotheses include the smaller, true effects). When tens of thousands of experiments run per year, even tiny SRM rates produce a steady stream of plausible-looking false positives. The operational moral is that good A/B platforms automate the diagnostics, and good experimenters declare an MDE and power before the test starts.
This stack cashes the S04 forward pointer: stratified A/B testing is blocked randomization with tech-industry vocabulary. Microsoft’s ExP, Netflix’s experimentation platform, Booking, LinkedIn’s XLNT, and Uber’s Morpheus all stratify on the same kinds of features S04’s blocked-RA estimator handled — activity decile × market × device, country × plan × pre-period-minutes, platform × device × country, and so on. The estimator is unchanged from S04:
\[\hat\tau_{\text{strat}} \;=\; \sum_k \pi_k \big(\bar Y^{(1)}_k - \bar Y^{(0)}_k\big), \qquad \pi_k = n_k / n,\]
with \(S_k\) indexing the strata and \(\pi_k\) their relative size. Exercise 1 reapplies the S04 Exercise 3 simulation under this jargon, confirming the variance ratio is identical.
CUPED — Controlled-experiment Using Pre-Experiment Data (Deng et al. 2013) — is the analysis-stage cousin of stratification. Instead of designing the randomization to balance on a pre-period covariate, CUPED leaves the randomization alone and adjusts the outcome at analysis time using the pre-experiment outcome \(Y^{\text{pre}}\). The estimator is one line:
\[Y^{\text{cuped}}_i \;=\; Y_i \;-\; \theta\big(Y^{\text{pre}}_i - \bar Y^{\text{pre}}\big), \qquad \theta \;=\; \frac{\text{Cov}(Y, Y^{\text{pre}})}{\text{Var}(Y^{\text{pre}})}.\]
The ATE is then estimated by the difference of group means on \(Y^{\text{cuped}}\). The variance of the adjusted outcome shrinks by a factor of \((1 - \rho^2)\) where \(\rho = \text{Cor}(Y, Y^{\text{pre}})\):
\[\text{Var}(Y^{\text{cuped}}) \;=\; \text{Var}(Y)\,\bigl(1 - \rho^2\bigr).\]
A small algebra check: \(\text{Var}(Y^{\text{cuped}}) = \text{Var}(Y - \theta Y^{\text{pre}}) = \text{Var}(Y) - 2\theta\,\text{Cov}(Y, Y^{\text{pre}}) + \theta^2 \text{Var}(Y^{\text{pre}})\). Plug in the optimal \(\theta = \text{Cov}(Y, Y^{\text{pre}})/\text{Var}(Y^{\text{pre}})\) and the second and third terms combine to \(-\text{Cov}(Y, Y^{\text{pre}})^2 / \text{Var}(Y^{\text{pre}}) = -\rho^2 \text{Var}(Y)\), giving the \((1 - \rho^2)\) result.
The connection to S04 is immediate: setting \(X^{\text{pre}} = Y^{\text{pre}}\) in 2’s interaction estimator (with mean-centred covariates) reproduces CUPED with a single regression. CUPED is Lin (2013), specialised to the single-covariate case where the covariate is the pre-experiment outcome. The reason CUPED is presented as a separate method is operational: production platforms apply the \(\theta\)-adjustment in a single pass over a pre-computed \(Y^{\text{pre}}\) table, which is much cheaper than fitting per-experiment interaction OLS at petabyte scale.
The Microsoft headline number: \(\rho \approx 0.7\) for typical engagement metrics, giving \(\text{Var}(Y^{\text{cuped}}) / \text{Var}(Y) \approx 1 - 0.49 \approx 0.51\) — roughly a 50% variance reduction, which under the standard \(n \propto 1/\text{Var}\) scaling translates into 2× experiment throughput per quarter at fixed power (Deng et al. 2013). Microsoft attributes the doubling of Bing’s experimentation velocity post-2013 to the rollout of this single technique.
CUPED’s pre-period outcome is only one of many predictors of \(Y\). Three production extensions push further.
CUPAC — DoorDash’s “Control Using Predictions As Covariates” — replaces \(Y^{\text{pre}}\) with a gradient-boosted regression prediction \(\hat f(X^{\text{pre}})\) trained on a rich set of pre-period features (user activity decile, market, device, hour-of-day patterns, prior-week order counts). When the predictor explains substantially more variance than the raw \(Y^{\text{pre}}\) alone, the residual variance shrinks further and the SE drops. DoorDash reports approximately 40% additional test-length reduction on top of CUPED in their production rollout (DoorDash Engineering Blog, “Improving Experimental Power Through Control Using Predictions As Covariates”). The estimator is still unbiased under random assignment for the same reason CUPED is: \(\hat f(X^{\text{pre}})\) is a function of pre-treatment variables and is therefore exogenous to \(D\).
MLRATE — Machine Learning Regression-Adjusted Treatment Effect (Guo et al. 2021) — is the academically-grounded version. It uses cross-fitting and a doubly-robust adjustment to remain unbiased even when the ML predictor is misspecified. On 48 Facebook metrics, MLRATE reported variance reductions of up to 70%; the cross-fitting machinery directly previews S08’s Double Machine Learning treatment of nuisance functions. The cost is operational complexity — one extra training pass per experiment per cross-fitting fold.
CURE — Statsig’s April 2025 release — builds on CUPED with a non-parametric residual model that captures non-linear relationships between \(Y^{\text{pre}}\) and \(Y\). Statsig reports 10–20% additional variance reduction on top of CUPED on its production traffic (“Announcing CURE”, Statsig Engineering Blog, 2025). The method is recent and not yet peer-reviewed; flag accordingly when citing in coursework.
(Kohavi, Tang, and Xu 2020, ch. 22) gives the practitioner discipline. Three items must hold for CUPED to be valid and effective.
The closing punchline: variance reduction is experimental throughput. Every method in this stack converts variance into more experiments per quarter; CUPED’s ~50% conversion is the well-trodden floor, and CUPAC / MLRATE / CURE push further when the operational cost is justified.
S04 left two threads dangling. The first was Fisher’s randomization test from the previous session, which we said would return next session as the engine of A/B-test p-values. The second was the previous session’s table of CI misinterpretations, which we said would be cured by the Bayesian credible interval introduced here. Both promises close in this stack.
The textbook A/B test of S04 is the fixed-horizon mode: declare a target sample size \(n^*\) before the experiment starts, run until \(n = n^*\), run a single two-sided t-test (or HC-robust OLS) on the OEC, and decide. The t-test mechanics are unchanged from the previous session’s hypothesis-testing procedure; the new question is how to choose \(n^*\) in the first place.
The minimum detectable effect (MDE) framework gives the answer. Suppose we want to detect a true effect of size \(\delta\) at significance level \(\alpha\) (typically 0.05) with power \(1 - \beta\) (typically 0.80). The required per-arm sample size is
\[n \;\approx\; \frac{2\,(z_{\alpha/2} + z_\beta)^2 \,\sigma^2}{\delta^2},\]
where \(\sigma^2\) is the within-arm variance of the OEC. With \(\alpha = 0.05\) (so \(z_{\alpha/2} = 1.96\)) and \(\beta = 0.20\) (so \(z_\beta = 0.84\)), the numerator’s \(z\)-sum is \(2.80\), and \(n \approx 15.7 \,\sigma^2 / \delta^2\) per arm. For a metric with coefficient of variation \(\sigma/\mu = 2\) and a target relative lift of 1% (\(\delta = 0.01\mu\)), plugging in \(\sigma^2/\delta^2 = 4\mu^2/(0.01\mu)^2 = 40{,}000\) gives roughly 627,000 users per arm (~1.25M total). This is why platform-scale \(n\) in the high six figures or low millions is normal: tiny relative effects on noisy engagement metrics demand exactly this much data — and the variance-reduction methods are how teams compress the requirement. R’s pwr::pwr.t.test() and Python’s statsmodels.stats.power.tt_ind_solve_power() invert the formula for you.
Two practical consequences. First, power is a planning tool, not a result: report it before the experiment starts, not after. Second, do not pre-emptively halt the experiment before reaching \(n^*\) just because an early-looking effect “looks good.” That move is the peeking pathology, and the fixed-horizon machinery is mathematically defenceless against it.
The Fisher permutation test — the permutation / Fisher randomization test from the previous session — is the second move here. The permutation test you saw last session is exactly the engine of Optimizely’s exact A/B p-values. When the classical t-test assumptions fail — small \(n\), heavy-tailed outcomes, contaminated metrics — modern A/B platforms (Optimizely, VWO) fall back to a Monte Carlo Fisher randomization test that produces an exact finite-sample p-value: re-shuffle the treatment labels \(B\) times, compute \(\hat\tau_\pi\) on each shuffle, report the fraction with \(|\hat\tau_\pi| \geq |\hat\tau_{\text{obs}}|\). The S04 Fisher slide and this fixed-horizon block are the same procedure, packaged for two different audiences.
Suppose you set up a fixed-horizon A/B test with \(\alpha = 5\%\) but, because the dashboard updates daily, you give yourself permission to “stop early if the result looks definitive.” Each day’s t-test, viewed in isolation, has Type-I rate 5%. The probability that at least one of 30 daily tests rejects under the null is much higher — in the dependent-test setting of an A/B platform, simulations show the cumulative false-positive rate at 30 daily peeks is roughly 30%. The five-, ten-, and thirty-day curves climb steeply: about 14% at five peeks, 22% at ten, 30% at thirty.
The peeking pathology is not academic. 7 ran a field study on 4,964 A/B tests at a major online marketplace and found that 28–37% of nominally significant results were false discoveries — a direct consequence of peeking and selective stopping interacting with the marketplace’s continuous-monitoring dashboard.
The simulation below reproduces the inflated rate empirically. Run a moderate number of null-hypothesis A/B tests, check daily for 30 days, count the fraction where any peek rejects at \(p < 0.05\).
set.seed(2026)
peeks <- 30
n_per_arm_per_peek <- 50
trials <- 500 # bump to 10,000 for the headline ~30% on a fast machine
any_reject <- vapply(seq_len(trials), function(i) {
reject <- FALSE
Yt <- numeric(0); Yc <- numeric(0)
for (k in seq_len(peeks)) {
Yt <- c(Yt, rnorm(n_per_arm_per_peek))
Yc <- c(Yc, rnorm(n_per_arm_per_peek))
if (t.test(Yt, Yc)$p.value < 0.05) { reject <- TRUE; break }
}
reject
}, logical(1))
mean(any_reject)[1] 0.294import numpy as np
from scipy import stats
rng = np.random.default_rng(2026)
peeks, n_step, trials = 30, 50, 500 # bump trials to 10,000 for headline ~30%
any_reject = np.zeros(trials, dtype=bool)
for i in range(trials):
Yt, Yc = np.empty(0), np.empty(0)
for k in range(peeks):
Yt = np.concatenate([Yt, rng.normal(size=n_step)])
Yc = np.concatenate([Yc, rng.normal(size=n_step)])
if stats.ttest_ind(Yt, Yc).pvalue < 0.05:
any_reject[i] = True
break
print(round(any_reject.mean(), 3))0.256The reported fraction will land in the 0.20–0.30 range at trials = 500; bumping to trials = 10000 tightens it on the headline ~0.30 number. One nominal 5% test, run 30 times with early stopping, becomes a 30% test in practice.
The peeking pathology has three modern cures: alpha-spending boundaries, sequential / anytime-valid confidence sequences, and Bayesian decision rules. The latter two are this stack’s main material.
A confidence sequence \(\{C_n\}_{n \geq 1}\) is a sequence of intervals with the time-uniform coverage property
\[\Pr\big(\theta \in C_n \;\forall n \geq 1\big) \;\geq\; 1 - \alpha.\]
Compare this to a fixed-horizon CI, which has \(\Pr(\theta \in C_{n^*}) \geq 1 - \alpha\) for one specific \(n^*\) chosen in advance. The confidence sequence permits stopping at any \(n\) while preserving the family-wise error rate — it is, by construction, peeking-safe (Howard et al. 2021).
9 derive an especially tight family of confidence sequences by re-framing inference as a betting game: at each new observation, an imaginary gambler places a fractional bet against the null. Their accumulated wealth is a non-negative martingale, and Ville’s inequality translates the wealth into an anytime-valid \(p\)-value. Inverting the procedure — collecting the set of candidate parameter values whose wealth process has not yet exceeded \(1/\alpha\) — gives a confidence sequence that is typically 20–40% tighter than the mixture-supermartingale bounds of 8 on bounded outcomes. Both constructions are production-ready in the confseq R package and betting-cs Python implementations.
Sequential is now the industry default. As of 2026, Eppo, Statsig, and Optimizely all ship sequential-by-default; fixed-horizon is the opt-in mode. The reversal is recent — as recently as 2020, fixed-horizon was the default and sequential was an “advanced” option. The Berman & Van den Bulte 2022 false-discovery evidence drove the shift. Microsoft ExP is the lone holdout among major vendors, still defaulting to fixed-horizon with planned-look alpha-spending (a legacy of its 2007 internal design).
Sequential CSs are not free. At the planned sample size \(n^*\), a confidence sequence is typically 15–25% wider than the fixed-horizon CI at the same \(n^*\) — the law-of-iterated-logarithm correction that any anytime-valid procedure must pay (the half-width scales as \(\sqrt{2 \log \log n / n}\) instead of \(1/\sqrt{n}\)). The honest trade-off: if you genuinely run a one-shot experiment with no monitoring, fixed-horizon delivers a tighter interval. If you’re going to monitor and may stop early, the peeking-freedom more than pays for the extra width, since you typically stop well before \(n^*\).
Netflix Sequential is the production capstone. 10 unify the two threads of this deck — variance reduction (Stack 2) and anytime-valid inference (this stack) — into a single estimator: the CUPED-adjusted residual \(Y^{\text{cuped}}\) feeds into a betting confidence sequence that the experimenter can monitor without peeking penalty. The Netflix Tech Blog “Sequential Testing at Netflix” Pt. 1 (continuous outcomes) and Pt. 2 (counting processes for clicks/conversions) walk through the production code, including the bias-correction term that the CUPED residual introduces into the martingale construction. The operational payoff: variance reduction and sequential testing compound multiplicatively. CUPED halves variance; sequential testing eliminates the peeking penalty. Together they reduce time-to-decision by ~4× versus naïve fixed-horizon analysis.
The previous session’s table of CI misinterpretations gave you the four-row CI-misreadings table: students consistently want frequentist CIs to deliver direct probability statements about the parameter, and the frequentist machinery refuses. The Bayesian credible interval is the cure. Here is the three-column version that completes the S04 contract.
The mechanics are old and clean. For a binary conversion-rate experiment with \(\theta_A\) and \(\theta_B\) the per-arm true conversion rates, the Beta-Binomial conjugate gives a closed-form posterior. With a uniform \(\text{Beta}(1, 1)\) prior and \(s_j\) successes from \(n_j\) trials in arm \(j\), the posterior is \(\theta_j \mid \text{data} \sim \text{Beta}(1 + s_j, 1 + n_j - s_j)\) — where \(s_A\) = observed conversions in arm A, \(n_A\) = users in arm A (likewise for B). Draw \(M = 100{,}000\) samples from each marginal posterior and compute \(\Delta^{(m)} = \theta_B^{(m)} - \theta_A^{(m)}\). The posterior probability that B beats A is the empirical fraction of \(\Delta^{(m)} > 0\); the 95% credible interval is the empirical 2.5%–97.5% quantile range of \(\Delta^{(m)}\).
The expected-loss decision rule ships only when the cost of being wrong is small enough. Define the loss of shipping B when A is actually better as \(\text{loss}(B \text{ over } A) = \max(\theta_A - \theta_B, 0)\). The posterior expected loss is \(\mathbb{E}[\max(\theta_A - \theta_B, 0) \mid \text{data}] \approx \frac{1}{M} \sum_m \max(\theta_A^{(m)} - \theta_B^{(m)}, 0)\). The decision rule is: ship B if expected loss < \(\varepsilon\), where \(\varepsilon\) is a business-set threshold calibrated to the metric scale (for absolute conversion rates around 6%, \(\varepsilon = 10^{-4}\) is the Eppo/Statsig default; for relative-lift metrics, \(\varepsilon \approx 0.005\) of baseline is more common). Both vendors implement this as the default Bayesian decision rule.
Exercise 4 walks through the full Beta-Binomial / expected-loss / decision-rule pipeline on the same conversion-rate scenario (Kohavi, Tang, and Xu 2020, ch. 17). The three modes side by side. Fixed-horizon \(p\)-values, confidence sequences, and Bayesian credible intervals answer different questions: fixed-horizon delivers \(P(\text{data} \mid H_0)\) at a single pre-committed \(n\) and is not peeking-safe; confidence sequences deliver \(P(\theta \in C_n \;\forall n) \geq 1 - \alpha\) and are peeking-safe by construction, at the cost of 15–25% wider intervals; credible intervals deliver the direct \(P(\theta \in [L,U] \mid \text{data}) = 0.95\) statement and are coherent under any stopping rule, with widths roughly matching fixed-horizon under a flat prior. A decision heuristic: pre-committed one-shot launches (regulatory submission, single launch gate) → fixed-horizon for tightness; continuously-monitored experiments → sequential for peeking-safety; multi-armed asymmetric-cost decisions → Bayesian expected loss. In 2026 tech, sequential is the safe default ~80% of the time; Bayesian is gaining ground in conversion-rate testing; fixed-horizon survives in regulated industries.
Even with a perfect OEC, a tuned CUPED, and a sequential-valid inference mode, the textbook A/B contract fails in four distinct ways. Each is a separate diagnostic.
A sample ratio mismatch is a discrepancy between the intended assignment ratio (typically 50/50) and the observed count of users in each arm. SRMs occur in roughly 6% of Microsoft A/B tests and 10% of LinkedIn’s zoomed-in / opt-in tests (Fabijan et al. 2019). Causes: bot filtering applied asymmetrically across arms, redirect failures that drop one arm’s users at the load balancer, late-binding instrumentation that misses early-treated users, opt-out flows that disproportionately exit one arm. The diagnostic is a chi-squared test on the observed counts against the intended ratio:
# Intended 50/50, observed: 5640 control vs 6360 treatment
chisq.test(c(5640, 6360), p = c(0.5, 0.5))
Chi-squared test for given probabilities
data: c(5640, 6360)
X-squared = 43.2, df = 1, p-value = 4.942e-11from scipy.stats import chisquare
print(chisquare(f_obs=[5640, 6360], f_exp=[6000, 6000]))Power_divergenceResult(statistic=np.float64(43.2), pvalue=np.float64(4.942068131098458e-11))The Kohavi threshold for action is p < 0.001 (Kohavi, Tang, and Xu 2020, ch. 21). When SRM is detected, the only correct response is to stop and debug the pipeline. Do not interpret the experiment’s headline result — the bias from differential dropout almost always exceeds the apparent effect size, and reporting a “winner” from an SRM’d experiment is the single most common form of A/B malpractice. Exercise 5 trains exactly this diagnostic on three datasets with intended 50/50 splits but different observed deviations.
The Stable Unit Treatment Value Assumption (SUTVA) says unit \(i\)’s outcome depends only on unit \(i\)’s treatment — \(Y_i = Y_i(D_i)\), not \(Y_i(D_1, \ldots, D_n)\). Three canonical violation patterns in tech:
In each case, \(Y_i\) depends on \(D_j\) for some \(j \neq i\) — SUTVA fails. The DAG has cross-arrows \(D_i \to Y_j\) and \(D_j \to Y_i\) for any pair \(i, j\) sharing a market, a feed, or a supply pool. User-level randomization estimates an effect diluted by spillover — usually toward zero in two-sided marketplaces (untreated users absorb part of the treatment’s positive externality), occasionally away from zero (e.g., when the treatment cannibalises a shared resource).
12 report a meta-experiment at Airbnb: they ran user-level randomization on half the markets and market-level cluster randomization on the other half, then compared the estimated treatment effects across designs. Cluster-randomized estimates were 32% larger than user-level estimates — evidence that the user-level design was attenuating the true effect through interference. The design-stage fix is cluster randomization: randomize at the level of the spillover unit (market, city, geo-region) rather than at the user level. The analysis-stage cost is reduced effective sample size (the variance scales with the number of clusters, not the number of users), but the bias reduction is large enough that two-sided marketplaces routinely accept the trade.
Switchback is the canonical SUTVA fix when all users share a single global treatment (e.g., a dispatch algorithm, a marketplace pricing rule). Instead of splitting users, split time: alternate the global treatment between \(D = 1\) and \(D = 0\) in short windows. DoorDash runs dispatch-algorithm tests in 30-minute windows and reports approximately 30% faster iteration versus user-level designs (DoorDash Engineering Blog, “Switchback Tests and Randomized Experimentation under Network Effects at DoorDash”). The canonical paper is 11, who derive design and analysis for switchback experiments with autocorrelated outcomes.
The cost is temporal autocorrelation: outcomes within a window are not independent, and outcomes near a window boundary carry over from the previous window’s treatment. 13 show that the optimal window length is a function of the carryover decay rate: if the treatment effect on the relevant outcome decays in roughly 5 minutes, 30-minute windows leave about 80% of the post-boundary variance “clean” and the bias is small; if decay is slower (e.g., 20 minutes), 30-minute windows leave only the last 10 minutes clean and the bias dominates. Picking the window length is therefore not a free parameter — it requires a pre-experiment estimate of the carryover function.
Some decisions need to be made on metrics that take months to mature: 90-day retention, lifetime value, churn at year one. 14’s surrogate index method predicts the long-horizon metric from a vector of short-horizon proxies (sessions per user, day-7 retention, week-2 engagement) using a regression trained on historical data. Across 1,098 test arms from 200 Netflix A/B tests, 14-day surrogate predictions agreed with 63-day ground truth on ~95% of launch decisions — a useful decision-time accelerator. But the recall on actual launch winners was only 65–79%: 21–35% of features that would have won at 63 days were rejected at 14 days based on the surrogate.
The right reading is operational: surrogate methods save decision-time, not bias. A team that ships surrogate-positive features and reverts surrogate-negative features will be roughly right on 95% of decisions but will miss one in five real wins. Whether that trade is acceptable depends on the cost of a missed launch versus the cost of a delayed launch — a business decision, not a statistical one. The slide deck flags this as use with caution.
Heterogeneous treatment effects –> S09. Today’s A/B framework assumes you ship the feature globally if the average effect is positive. In practice, Netflix’s HTE Tech Blog catalogues persistent across-segment heterogeneity: firmware-dependent UI features that win on iOS but lose on Android; country-level homepage variants where the same change wins in Latin America and loses in Northern Europe; push-vs-email notification policies that segment by tenure. S09 introduces formal HTE estimators (causal forests, R-learner, meta-learners) that turn these segment-level patterns into ship-by-segment policies.
Bandits as an alternative. Multi-armed bandits trade off exploration and exploitation continuously and can outperform A/B at very-high-traffic homepage-personalisation problems — Spotify’s 2025 contextual-bandit deployment is the current state-of-the-art; Stitch Fix’s algorithms blog has an accessible architecture write-up. Bandits are out of scope for this course (A/B remains the default for one-shot product launches and any inference task where you need a clean estimate of a causal effect), but the one-slide flag on the deck prevents the perennial “why not bandits” question.
Geo-experiments –> S11. When user-level randomization is impossible (advertising attribution, brand campaigns, regulated industries) tech firms run geo-experiments: ship the feature in a subset of cities or countries and use the un-shipped geographies as comparison units. The analysis machinery is synthetic control — a weighted combination of un-shipped units that matches the treated unit’s pre-period trajectory. S11 develops this estimator in full.
Assignment II hook. All of S04 and S05 are examined in Assignment II, due 2026-05-31. Expect to (a) run an SRM check on a provided dataset; (b) compute CUPED-adjusted ATE with HC-robust SE; (c) pick one inference mode and justify the choice in writing; (d) report a brief HTE exploration as a closing forward-pointer to S09.
Five exercises, increasing in difficulty and total budget ≈ 145 minutes:
| # | Topic | Slide stack | Time |
|---|---|---|---|
| 1 | Stratified A/B in design and analysis (two stratifiers × two designs) | Variance Reduction | ~25 min |
| 2 | CUPED –> CUPAC on the blocked-RA dataset | Variance Reduction | ~30 min |
| 3 | Sample size, peeking, and confidence sequences | Inference Modes | ~40 min |
| 4 | Bayesian A/B with expected-loss decision rule | Inference Modes | ~30 min |
| 5 | SRM diagnostic on three deliberately-flawed datasets | When A/B Breaks | ~30 min |
All exercises are bilingual: every R chunk has a Python equivalent in a panel-tabset.
Goal / LO. Industry uses “stratified A/B” for two procedures with the same asymptotic variance but different mechanics:
The exercise builds one DGP with potential outcomes \((Y(0), Y(1))\) available, then constructs three observed datasets from it: one under simple RA, two under blocked RA on different stratification keys. The same four estimators are applied to each combination, producing three punchlines from the resulting tables:
Y ~ D + factor(S) regression. The two diverge whenever \(p_k\) varies across strata.os_type (binary; \(R^2 \approx 0.06\) on \(Y\)) and activity_decile (10 deciles of n_sess_pre; \(R^2 \approx 0.32\) on \(Y\)). The strong-stratifier table’s worst estimator beats the weak-stratifier table’s best estimator — variance reduction is governed by the predictive power of \(S\), not the cleverness of the formula applied afterwards.The blocked-on-os_type dataset built here flows directly into Exercise 2, where CUPED and CUPAC layer on top of the same blocked-RA design.
The DGP mirrors the S04 in-app feature scenario but exposes both potential outcomes \((Y_i(0), Y_i(1))\) so we can apply three assignment mechanisms to the same population. Treatment effects are heterogeneous by os_type (\(\tau = 2\) min for Android, \(\tau = 4\) min for iOS) — this is what makes the precision-weighted vs sample-share-weighted estimators numerically diverge under simple RA. Two pre-period stratification keys are constructed: the binary os_type (weak predictor of \(Y\)) and a 10-level activity_decile derived from n_sess_pre (strong predictor).
suppressMessages({library(dplyr); library(estimatr); library(randomizr)})
set.seed(20260519L)
n <- 2000
os_type <- rbinom(n, 1, 0.60)
n_sess_pre <- rgamma(n, shape = 2, scale = 5)
activity_decile <- ntile(n_sess_pre, 10) # 10 deciles of pre-period sessions
Y_pre <- 18 + 6 * os_type + 0.7 * n_sess_pre + rnorm(n, sd = 5)
Y0 <- 8 + 0.6 * Y_pre + 0.2 * n_sess_pre + rnorm(n, sd = 5)
tau_k_dgp <- ifelse(os_type == 0, 2, 4) # heterogeneous treatment effect
Y1 <- Y0 + tau_k_dgp
set.seed(20260520L)
D_simple <- rbinom(n, 1, 0.5) # design 1: simple RA
D_blocked_os <- block_ra(blocks = os_type, prob = 0.5) # design 2: blocked on os_type
D_blocked_dec <- block_ra(blocks = activity_decile, prob = 0.5) # design 3: blocked on decile
mk_df <- function(D) tibble::tibble(user_id = 1:n, os_type, activity_decile,
n_sess_pre, Y_pre,
D = D, Y = D * Y1 + (1 - D) * Y0)
inapp_simple <- mk_df(D_simple)
inapp_blocked_os <- mk_df(D_blocked_os)
inapp_blocked_dec <- mk_df(D_blocked_dec)
# For Exercise 2 continuity, the canonical "blocked" dataset uses os_type
inapp_blocked <- inapp_blocked_os
cat(sprintf("Sample ATE = %.3f\n", mean(Y1 - Y0)))Sample ATE = 3.200cat(sprintf("R^2 of os_type on Y0 = %.3f (weak)\n",
summary(lm(Y0 ~ factor(os_type)))$r.squared))R^2 of os_type on Y0 = 0.058 (weak)cat(sprintf("R^2 of activity_decile on Y0 = %.3f (strong)\n",
summary(lm(Y0 ~ factor(activity_decile)))$r.squared))R^2 of activity_decile on Y0 = 0.323 (strong)cat(sprintf("cor(Y, Y_pre) = %.3f (for Ex 2 CUPED)\n",
cor(inapp_simple$Y, inapp_simple$Y_pre)))cor(Y, Y_pre) = 0.737 (for Ex 2 CUPED)import numpy as np, pandas as pd
import statsmodels.formula.api as smf
rng = np.random.default_rng(20260519)
n = 2000
os_type = rng.binomial(1, 0.60, n)
n_sess_pre = rng.gamma(2, 5, n)
activity_decile = pd.qcut(n_sess_pre, 10, labels=False) + 1
Y_pre = 18 + 6 * os_type + 0.7 * n_sess_pre + rng.normal(0, 5, n)
Y0 = 8 + 0.6 * Y_pre + 0.2 * n_sess_pre + rng.normal(0, 5, n)
tau_k_dgp = np.where(os_type == 0, 2, 4)
Y1 = Y0 + tau_k_dgp
def block_ra(blocks, prob, rng):
D = np.zeros(len(blocks), dtype=int)
for b in np.unique(blocks):
idx = np.where(blocks == b)[0]
k = int(round(len(idx) * prob))
D[rng.choice(idx, k, replace=False)] = 1
return D
rng2 = np.random.default_rng(20260520)
D_simple = rng2.binomial(1, 0.5, n)
D_blocked_os = block_ra(os_type, 0.5, rng2)
D_blocked_dec = block_ra(activity_decile, 0.5, rng2)
def mk_df(D):
return pd.DataFrame({"user_id": np.arange(n), "os_type": os_type,
"activity_decile": activity_decile,
"n_sess_pre": n_sess_pre, "Y_pre": Y_pre,
"D": D, "Y": D * Y1 + (1 - D) * Y0})
inapp_simple = mk_df(D_simple)
inapp_blocked_os = mk_df(D_blocked_os)
inapp_blocked_dec = mk_df(D_blocked_dec)
inapp_blocked = inapp_blocked_os # canonical blocked dataset for Ex 2
r2_os = smf.ols("Y0 ~ C(os_type)", data=pd.DataFrame({"Y0":Y0,"os_type":os_type})).fit().rsquared
r2_dec = smf.ols("Y0 ~ C(activity_decile)", data=pd.DataFrame({"Y0":Y0,"activity_decile":activity_decile})).fit().rsquared
print(f"Sample ATE = {(Y1-Y0).mean():.3f}")Sample ATE = 3.153print(f"R^2 of os_type on Y0 = {r2_os:.3f} (weak)")R^2 of os_type on Y0 = 0.058 (weak)print(f"R^2 of activity_decile on Y0 = {r2_dec:.3f} (strong)")R^2 of activity_decile on Y0 = 0.298 (strong)print(f"cor(Y, Y_pre) = {np.corrcoef(inapp_simple.Y, inapp_simple.Y_pre)[0,1]:.3f} (for Ex 2 CUPED)")cor(Y, Y_pre) = 0.723 (for Ex 2 CUPED)# Helper that applies the four estimators to a dataset using `key` as the strata column.
# Works for any K-level stratifier; uses estimatr::lm_lin for the saturated (d) form.
four_estimators <- function(d, key) {
k_sym <- rlang::sym(key)
# (a) Unstratified
fit_a <- lm_robust(Y ~ D, data = d, se_type = "HC1")
# (b) Manual weighted: sum_k pi_k * tau_hat_k
strat <- d |> group_by(!!k_sym) |>
summarise(tau_k = mean(Y[D == 1]) - mean(Y[D == 0]),
n_k = n(), p_k = mean(D), .groups = "drop") |>
mutate(pi_k = n_k / sum(n_k))
tau_b <- sum(strat$pi_k * strat$tau_k)
# (c) OLS additive with strata dummies -- YOUR CODE
# (d) Lin (2013) via estimatr::lm_lin with factor(key) -- YOUR CODE
# Return a 4-row table with tau_hat and HC1 SE
}
# Apply to (inapp_simple, inapp_blocked_os, inapp_blocked_dec) under each key.
# Compare especially: does (b) = (d) for the strong stratifier too?
# Under blocked RA on a key, do all four estimators coincide?def four_estimators(d, key):
fit_a = smf.ols("Y ~ D", data=d).fit(cov_type="HC1")
g = d.groupby(key)
tau_k = g.apply(lambda x: x.loc[x.D == 1, "Y"].mean() - x.loc[x.D == 0, "Y"].mean())
pi_k = g.size() / len(d)
tau_b = float((pi_k * tau_k).sum())
# (c) Additive: smf.ols("Y ~ D + C(key)", ...) -- YOUR CODE
# (d) Lin saturated: mean-center K dummies, then `Y ~ D + dummies_c + D:dummies_c` -- YOUR CODE
# Return a 4-row DataFrameTo see what design-stage blocking buys under each stratifier, re-randomize 2,000 times under each of three designs — simple RA, blocked on os_type, blocked on activity_decile — holding the potential outcomes fixed. Histograms of the unstratified diff-in-means \(\hat\tau\) across draws expose the sampling-distribution argument the single-shot HC1 SEs cannot.
library(ggplot2); library(tidyr)
set.seed(20260521L)
reps <- 2000
po <- tibble::tibble(os_type, activity_decile, Y0, Y1)
draws <- replicate(reps, {
Ds <- rbinom(n, 1, 0.5)
Dbo <- block_ra(blocks = po$os_type, prob = 0.5)
Dbd <- block_ra(blocks = po$activity_decile, prob = 0.5)
Ys <- Ds * po$Y1 + (1 - Ds) * po$Y0
Ybo <- Dbo * po$Y1 + (1 - Dbo) * po$Y0
Ybd <- Dbd * po$Y1 + (1 - Dbd) * po$Y0
c(simple = mean(Ys [Ds == 1]) - mean(Ys [Ds == 0]),
`blocked / os` = mean(Ybo[Dbo == 1]) - mean(Ybo[Dbo == 0]),
`blocked / decile` = mean(Ybd[Dbd == 1]) - mean(Ybd[Dbd == 0]))
})
hist_data <- as.data.frame(t(draws)) |>
setNames(c("simple", "blocked / os_type", "blocked / activity_decile")) |>
pivot_longer(everything(), names_to = "design", values_to = "tau_hat") |>
mutate(design = factor(design,
levels = c("simple", "blocked / os_type", "blocked / activity_decile")))
ggplot(hist_data, aes(x = tau_hat, fill = design)) +
geom_histogram(position = "identity", alpha = 0.5, bins = 40, color = "white") +
geom_vline(xintercept = mean(Y1 - Y0), linetype = "dashed") +
scale_fill_manual(values = c("simple" = "#FF7E15",
"blocked / os_type" = "#7BAFD4",
"blocked / activity_decile" = "#0072B2")) +
labs(x = expression(hat(tau)~"(unstratified diff-in-means)"),
y = "count", fill = "Design",
caption = "Dashed line = sample ATE. 2,000 re-randomizations of fixed potential outcomes.") +
theme_minimal()
mc_se <- function(x) sd(x) / sqrt(2 * (length(x) - 1))
for (k in c("simple", "blocked / os", "blocked / decile"))
cat(sprintf("Empirical SD, %-18s = %.4f (MC-SE = %.4f)\n",
k, sd(draws[k, ]), mc_se(draws[k, ])))Empirical SD, simple = 0.3325 (MC-SE = 0.0053)
Empirical SD, blocked / os = 0.3283 (MC-SE = 0.0052)
Empirical SD, blocked / decile = 0.2809 (MC-SE = 0.0044)import matplotlib.pyplot as plt
rng3 = np.random.default_rng(20260521)
reps = 2000
sim_draws = np.empty(reps); blk_os = np.empty(reps); blk_dec = np.empty(reps)
for i in range(reps):
Ds = rng3.binomial(1, 0.5, n)
Dbo = block_ra(os_type, 0.5, rng3)
Dbd = block_ra(activity_decile, 0.5, rng3)
Ys = Ds * Y1 + (1 - Ds) * Y0
Ybo = Dbo * Y1 + (1 - Dbo) * Y0
Ybd = Dbd * Y1 + (1 - Dbd) * Y0
sim_draws[i] = Ys [Ds == 1].mean() - Ys [Ds == 0].mean()
blk_os[i] = Ybo[Dbo == 1].mean() - Ybo[Dbo == 0].mean()
blk_dec[i] = Ybd[Dbd == 1].mean() - Ybd[Dbd == 0].mean()
fig, ax = plt.subplots(figsize=(8.5, 3.8))
all_draws = np.concatenate([sim_draws, blk_os, blk_dec])
bins = np.linspace(all_draws.min(), all_draws.max(), 41)
ax.hist(sim_draws, bins=bins, alpha=0.5, color="#FF7E15",
label="simple", edgecolor="white")
ax.hist(blk_os, bins=bins, alpha=0.5, color="#7BAFD4",
label="blocked / os_type", edgecolor="white")
ax.hist(blk_dec, bins=bins, alpha=0.5, color="#0072B2",
label="blocked / activity_decile", edgecolor="white")
ax.axvline((Y1 - Y0).mean(), linestyle="--", color="black")
ax.set_xlabel(r"$\hat\tau$ (unstratified diff-in-means)")
ax.set_ylabel("count")
ax.legend(title="Design", loc="upper right")
plt.tight_layout(); plt.show()
for name, x in [("simple", sim_draws), ("blocked / os", blk_os),
("blocked / decile", blk_dec)]:
print(f"Empirical SD, {name:18s} = {x.std(ddof=1):.4f} "
f"(MC-SE = {x.std(ddof=1) / np.sqrt(2*(len(x)-1)):.4f})")Empirical SD, simple = 0.3323 (MC-SE = 0.0053)
Empirical SD, blocked / os = 0.3200 (MC-SE = 0.0051)
Empirical SD, blocked / decile = 0.2838 (MC-SE = 0.0045)The three-design histogram makes the lesson tactile:
os_type overlap almost perfectly — their empirical SDs differ by less than one MC standard error (~0.33 vs ~0.33). Blocking on a weak stratifier buys essentially zero design-stage variance reduction, exactly as the single-shot SE table predicted.activity_decile is visibly tighter — its empirical SD is ~0.28, a ~15% SD reduction over simple RA. The strong stratifier’s predictive power compounds into a real design-stage gain even for the simplest (unstratified) estimator.Combined message: the histogram and the SE table tell the same story by two different routes. Picking the stratifier well is what moves the needle; whether you stratify at design time vs analysis time is a second-order operational choice (with the four practical reasons in the Reflection above still favouring blocked RA).
The empirical SDs above (sampling distribution of the unstratified estimator across designs) and the HC1 SEs from the Solution tables (single-shot SE of the stratified estimators on simple-RA data) fit together as the cells of a design × analysis matrix:
| Naive estimator (a), \(\text{Var}(\hat\tau)\) | Stratified estimator (c), HC1 SE | |
|---|---|---|
| Simple RA | empirical SD ≈ 0.333 | (with os_type) ≈ 0.328; (with activity_decile) ≈ 0.286 |
Blocked RA on os_type |
empirical SD ≈ 0.328 | ≈ 0.327 |
Blocked RA on activity_decile |
empirical SD ≈ 0.281 | ≈ 0.285 |
The diagonal in bold is the Neyman–Cochran result made empirical: the naive estimator’s sampling SE under blocked-RA-on-\(S\) equals the analysis-stratified estimator’s SE under simple-RA-with-\(S\)-as-strata. Design-stage stratification and analysis-stage stratification are two roads to the same variance; on the matching diagonal, they coincide.
This is what your intuition predicts and what the single-shot (a) HC1 SE column visibly hides. Under blocked-RA-on-activity_decile, the true sampling SE of \(\hat\tau_{\text{naive}}\) is ≈ 0.28 — but the (a) row of the Solution table reports HC1 SE = 0.343, identical to the simple-RA row. The HC1 estimator computes residuals from lm(Y ~ D) regardless of how \(D\) was assigned; those residuals include between-stratum variability that the design has already balanced. So HC1 over-estimates the sampling variance of the naive estimator under blocked RA.
Goal / LO. Stack analysis-stage variance reduction on top of the design-stage blocking from Exercise 1. The same inapp_blocked dataset feeds four estimators:
Three things to verify:
os_type is already implicitly absorbed by blocked RA, but n_sess_pre carries unique signal beyond \(Y^{\text{pre}}\) — that’s CUPAC’s edge.Data. Reuses inapp_blocked built in Exercise 1 (no fresh DGP needed). This is the data-continuity payoff of the unified setup: same n = 2000 users, same potential outcomes, blocked-randomized on os_type, with \(Y^{\text{pre}}\) already available as a column.
Fill in the three estimators on inapp_blocked; assemble a 3-row results table with tau_hat, HC1 SE, and the variance ratio relative to the naive estimator.
# (a) Naive (blocked RA only) --- diff-in-means with HC1 SE
fit_naive <- # ___ YOUR CODE (use lm_robust(... , se_type = "HC1"))
# (b) Lin (2013) regression adjustment: mean-centred Y_pre + treatment interaction
inapp_blocked <- inapp_blocked |> mutate(Ypre_c = Y_pre - mean(Y_pre))
fit_lin <- # ___ YOUR CODE (lm_robust with Y ~ D * Ypre_c)
# (c) CUPED with theta = Cov(Y, Y_pre) / Var(Y_pre)
theta_cuped <- # ___ YOUR CODE
inapp_blocked <- inapp_blocked |>
mutate(Y_cuped = Y - theta_cuped * (Y_pre - mean(Y_pre)))
fit_cuped <- # ___ YOUR CODE (lm_robust of Y_cuped on D)
# Assemble a 3-row data.frame with: estimator, tau_hat, se_hc1, var_ratio
# var_ratio = (se_i / se_naive)^2# (a) Naive
fit_naive = # ___ YOUR CODE (smf.ols(..., data=inapp_blocked).fit(cov_type="HC1"))
# (b) Lin (2013): mean-centred Y_pre + treatment interaction
inapp_blocked["Ypre_c"] = inapp_blocked["Y_pre"] - inapp_blocked["Y_pre"].mean()
fit_lin = # ___ YOUR CODE (Y ~ D * Ypre_c)
# (c) CUPED
theta_cuped = # ___ YOUR CODE
inapp_blocked["Y_cuped"] = (inapp_blocked.Y
- theta_cuped * (inapp_blocked.Y_pre - inapp_blocked.Y_pre.mean()))
fit_cuped = # ___ YOUR CODE
# Build a 3-row DataFrame with: estimator, tau_hat, se_hc1, var_ratio.Goal / LO. Ground sample-size intuition (a), reproduce the inflated false-positive rate under peeking (b), and compare two cures: a fixed Pocock alpha-spending boundary (c) and an anytime-valid confidence sequence (d). The exercise demonstrates why sequential CSs are the modern industry default — not because Pocock is wrong, but because it is K-specific and you rarely know K in advance.
Setup / Task.
pwr::pwr.t.test() (R) or statsmodels.stats.power.tt_ind_solve_power() (Python), compute the per-arm sample size \(n^*\) needed to detect a \(\delta = 1\%\) relative lift on a continuous metric with \(\sigma / \mu = 2\), at \(\alpha = 0.05\) and power \(= 0.80\). Confirm the slide’s headline figure \(n^* \approx 6.3 \times 10^5\) per arm.n_per_arm_per_peek = 50, peek 30 times, stop at first \(|t| > 1.96\), and report the cumulative false-positive rate at each peek.library(ggplot2); library(tibble); library(pwr)
# (a) Sample size for delta = 1% relative lift, sigma/mu = 2, alpha = 0.05, power = 0.80
d <- 0.01 / 2 # Cohen's d = relative lift / coefficient of variation
mde_res <- # ___ YOUR CODE: pwr.t.test(d, sig.level, power, type = "two.sample")
# Print n* per arm (should be ≈ 6.3e5)
# (b)-(d) Peeking simulation. Build a function that takes a critical_value
# (constant for naive/Pocock; function of n for CS) and returns the cumulative
# FPR curve across 30 peeks of 50 obs/arm.
run_simulation <- function(critical_value, trials = 2000, peeks = 30, n_step = 50) {
# ___ YOUR CODE: loop trials × peeks, accumulate Yt/Yc by n_step per peek,
# compute t.test, compare |t| to (numeric or function-of-n) critical value,
# mark the FIRST rejection per trial, return cumulative FPR per peek.
}
# Confidence-sequence critical value (Howard et al. 2021 LIL form):
# c(n) = sqrt(2 * log(log(n) + 2) + log(2 / alpha))
cs_bound <- # ___ YOUR CODE (function of n)
# Run three simulations (same seed for paired comparison)
set.seed(2026); cum_fp_naive <- run_simulation(1.96)
set.seed(2026); cum_fp_pocock <- run_simulation(2.413)
set.seed(2026); cum_fp_cs <- run_simulation(cs_bound)
# Plot the three cumulative-FPR curves with a horizontal nominal-alpha line at 0.05.
# Produce a 3-row table with the final-peek FPR for each procedure.import numpy as np, pandas as pd
from scipy import stats
from statsmodels.stats.power import tt_ind_solve_power
import matplotlib.pyplot as plt
# (a) Sample size
d = 0.01 / 2
n_star = # ___ YOUR CODE: tt_ind_solve_power(effect_size=d, alpha, power, alternative)
# (b)-(d) Peeking simulation
def run_simulation(critical_value, trials=2000, peeks=30, n_step=50):
rng = np.random.default_rng(2026)
# ___ YOUR CODE: same logic as R version
pass
def cs_bound(n, alpha=0.05):
return # ___ YOUR CODE
cum_fp_naive = run_simulation(1.96)
cum_fp_pocock = run_simulation(2.413)
cum_fp_cs = run_simulation(cs_bound)
# Plot + summary table.library(ggplot2); library(tibble); library(pwr)
# (a) Sample size --- Cohen's d = (relative lift) / (coefficient of variation)
# Spec: delta = 1% relative lift, sigma/mu = 2, alpha = 0.05, power = 0.80
delta_rel <- 0.01 # 1% relative lift
cv <- 2 # coefficient of variation: sigma / mu
d <- delta_rel / cv # = 0.005 (mu cancels because both are proportional to mu)
mde_res <- pwr.t.test(d = d, sig.level = 0.05, power = 0.80,
type = "two.sample")
cat(sprintf("Cohen's d = %.4f\nn* per arm = %s\n",
d, format(ceiling(mde_res$n), big.mark = ",")))Cohen's d = 0.0050
n* per arm = 627,910# (b)-(d) Peeking simulation; critical_value is either a constant
# (naive / Pocock) or a function of the cumulative per-arm sample size (CS).
set.seed(2026)
run_simulation <- function(critical_value, trials = 2000, peeks = 30,
n_step = 50) {
ever_reject <- matrix(FALSE, nrow = trials, ncol = peeks)
for (i in seq_len(trials)) {
Yt <- numeric(0); Yc <- numeric(0); rejected <- FALSE
for (k in seq_len(peeks)) {
Yt <- c(Yt, rnorm(n_step))
Yc <- c(Yc, rnorm(n_step))
tt <- t.test(Yt, Yc)
cv <- if (is.function(critical_value))
critical_value(length(Yt)) else critical_value
if (!rejected && abs(tt$statistic) > cv) rejected <- TRUE
if (rejected) { ever_reject[i, k:peeks] <- TRUE; break }
}
}
colMeans(ever_reject)
}
cs_bound <- function(n, alpha = 0.05)
sqrt(2 * log(log(n) + 2) + log(2 / alpha))
set.seed(2026); cum_fp_naive <- run_simulation(1.96)
set.seed(2026); cum_fp_pocock <- run_simulation(2.413)
set.seed(2026); cum_fp_cs <- run_simulation(cs_bound)
df_plot <- tibble(
peeks = rep(seq_len(30), 3),
cum_fp = c(cum_fp_naive, cum_fp_pocock, cum_fp_cs),
procedure = factor(rep(c("Naive (z=1.96)",
"Pocock K=5 (z=2.413)",
"Confidence sequence (LIL)"), each = 30),
levels = c("Naive (z=1.96)",
"Pocock K=5 (z=2.413)",
"Confidence sequence (LIL)"))
)
ggplot(df_plot, aes(peeks, cum_fp, colour = procedure)) +
geom_line(linewidth = 0.9) +
geom_hline(yintercept = 0.05, linetype = "dashed", colour = "grey30") +
annotate("text", x = 30, y = 0.05, label = " nominal α = 0.05",
hjust = 1, vjust = -0.5, size = 3.2, colour = "grey30") +
scale_colour_manual(values = c("Naive (z=1.96)" = "#FF7E15",
"Pocock K=5 (z=2.413)" = "#005e73",
"Confidence sequence (LIL)" = "#00C1D4")) +
labs(x = "Number of peeks (days)", y = "Cumulative Type-I error rate",
colour = NULL) +
theme(legend.position = "bottom")
tibble(procedure = c("Naive (z=1.96)", "Pocock K=5 (z=2.413)",
"Confidence sequence (LIL)"),
final_fp = c(tail(cum_fp_naive, 1),
tail(cum_fp_pocock, 1),
tail(cum_fp_cs, 1)))# A tibble: 3 × 2
procedure final_fp
<chr> <dbl>
1 Naive (z=1.96) 0.297
2 Pocock K=5 (z=2.413) 0.116
3 Confidence sequence (LIL) 0.0545import numpy as np, pandas as pd
from scipy import stats
from statsmodels.stats.power import tt_ind_solve_power
import matplotlib.pyplot as plt
# (a) Sample size --- Cohen's d = (relative lift) / (coefficient of variation)
# Spec: delta = 1% relative lift, sigma/mu = 2, alpha = 0.05, power = 0.80
delta_rel = 0.01 # 1% relative lift
cv = 2 # coefficient of variation: sigma / mu
d = delta_rel / cv # = 0.005 (mu cancels because both are proportional to mu)
n_star = tt_ind_solve_power(effect_size=d, alpha=0.05, power=0.80,
alternative="two-sided")
print(f"Cohen's d = {d:.4f}")Cohen's d = 0.0050print(f"n* per arm = {int(np.ceil(n_star)):,}")n* per arm = 627,910# (b)-(d) Peeking simulation
def run_simulation(critical_value, trials=2000, peeks=30, n_step=50):
rng = np.random.default_rng(2026)
ever_reject = np.zeros((trials, peeks), dtype=bool)
for i in range(trials):
Yt, Yc = np.empty(0), np.empty(0); rejected = False
for k in range(peeks):
Yt = np.concatenate([Yt, rng.normal(size=n_step)])
Yc = np.concatenate([Yc, rng.normal(size=n_step)])
tt = stats.ttest_ind(Yt, Yc)
cv = (critical_value(len(Yt)) if callable(critical_value)
else critical_value)
if not rejected and abs(tt.statistic) > cv:
rejected = True
if rejected:
ever_reject[i, k:] = True; break
return ever_reject.mean(axis=0)
def cs_bound(n, alpha=0.05):
return np.sqrt(2 * np.log(np.log(n) + 2) + np.log(2 / alpha))
cum_fp_naive = run_simulation(1.96)
cum_fp_pocock = run_simulation(2.413)
cum_fp_cs = run_simulation(cs_bound)
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.plot(range(1, 31), cum_fp_naive, color="#FF7E15",
linewidth=1.8, label="Naive (z=1.96)")
ax.plot(range(1, 31), cum_fp_pocock, color="#005e73",
linewidth=1.8, label="Pocock K=5 (z=2.413)")
ax.plot(range(1, 31), cum_fp_cs, color="#00C1D4",
linewidth=1.8, label="Confidence sequence (LIL)")
ax.axhline(0.05, linestyle="--", color="grey", linewidth=0.9)
ax.annotate(r" nominal $\alpha = 0.05$", xy=(30, 0.05), xytext=(30, 0.052),
ha="right", va="bottom", color="grey", fontsize=9)
ax.set_xlabel("Number of peeks (days)")
ax.set_ylabel("Cumulative Type-I error rate")
ax.legend(loc="best")
plt.tight_layout(); plt.show()
print(pd.DataFrame({
"procedure": ["Naive (z=1.96)", "Pocock K=5 (z=2.413)",
"Confidence sequence (LIL)"],
"final_fp": [cum_fp_naive[-1], cum_fp_pocock[-1], cum_fp_cs[-1]]
})) procedure final_fp
0 Naive (z=1.96) 0.2850
1 Pocock K=5 (z=2.413) 0.1195
2 Confidence sequence (LIL) 0.0435Expected output.
(b)–(d) Three-line plot + summary table:
| Procedure | Critical value | Cumulative FPR at peek 30 |
|---|---|---|
| Naive | \(z = 1.96\) (constant) | \(\approx 0.29\) |
| Pocock \(K = 5\) | \(z = 2.41\) (constant) | \(\approx 0.13\) |
| Confidence sequence (LIL) | \(c(n_k)\) (adaptive) | \(\approx 0.05\) |
The plot shows the naive procedure climbing steeply, the Pocock procedure attenuating to ~13% (not 5%), and the CS bound oscillating near the dotted 5% reference line.
Reflection. Two operational lessons.
confseq R / Python packages expose the rigorous versions. The cost of anytime-validity is a slightly wider rejection threshold at the planned endpoint (~15–25% wider than fixed-horizon, per slide L1484) — but you typically stop well before that endpoint anyway.By 2026, Eppo, Statsig, and Optimizely all ship CS-based sequential testing as the default mode. Microsoft ExP is the lone holdout among major vendors, still on alpha-spending with planned interim looks.
Goal / LO. Construct a Bayesian credible interval and an expected-loss decision rule for a conversion-rate A/B test, and contrast all three quantities (frequentist CI, credible interval, expected loss) on the same data.
Setup / Task. Control arm: \(n_A = 12{,}000\), \(s_A = 720\) conversions (\(\hat p_A = 0.060\)). Treatment arm: \(n_B = 12{,}000\), \(s_B = 762\) conversions (\(\hat p_B = 0.0635\)). (a) Compute the frequentist two-proportion z-test CI for \(p_B - p_A\). (b) Place a uniform \(\text{Beta}(1,1)\) prior on each rate, draw \(M = 100{,}000\) samples from the posteriors \(\text{Beta}(1 + s_j, 1 + n_j - s_j)\), and compute the 95% credible interval for \(p_B - p_A\). (c) Compute the posterior probability \(\Pr(p_B > p_A \mid \text{data})\). (d) Compute the expected loss of shipping B, \(\mathbb{E}[\max(p_A - p_B, 0) \mid \text{data}]\), and apply the ship-if-loss-< 0.0001 decision rule. Produce a comparison table and write a one-sentence decision recommendation.
set.seed(2026)
n_A <- 12000; s_A <- 720
n_B <- 12000; s_B <- 762
# (a) Frequentist 95% CI for p_B - p_A (two-proportion z-test)
p_A <- s_A / n_A; p_B <- s_B / n_B
diff_hat <- p_B - p_A
se_diff <- # ___ YOUR CODE: sqrt(p_A(1-p_A)/n_A + p_B(1-p_B)/n_B)
ci_freq <- # ___ YOUR CODE: diff_hat ± 1.96 * se_diff
# (b) Beta(1,1) prior + Beta-Binomial posterior, M = 100k draws per arm
M <- 100000
theta_A <- # ___ YOUR CODE: rbeta(M, 1 + s_A, 1 + n_A - s_A)
theta_B <- # ___ YOUR CODE
diff_post <- theta_B - theta_A
ci_cred <- # ___ YOUR CODE: 2.5%/97.5% quantiles of diff_post
# (c) Posterior probability that B beats A
p_B_better <- # ___ YOUR CODE
# (d) Expected loss of shipping B; ship-if-loss < 1e-4
expected_loss <- # ___ YOUR CODE: mean(pmax(theta_A - theta_B, 0))
decision <- # ___ YOUR CODE: "SHIP B" or "DO NOT SHIP"
# Assemble a 5-row comparison table (frequentist CI, credible interval,
# P(B>A), expected loss, decision) and write a one-sentence recommendation.import numpy as np, pandas as pd
rng = np.random.default_rng(2026)
n_A, s_A = 12000, 720
n_B, s_B = 12000, 762
# (a) Frequentist CI
p_A, p_B = s_A / n_A, s_B / n_B
diff_hat = p_B - p_A
se_diff = # ___ YOUR CODE
ci_freq = # ___ YOUR CODE
# (b)-(d) Beta-Binomial posterior, M = 100k
M = 100000
theta_A = # ___ YOUR CODE: rng.beta(1 + s_A, 1 + n_A - s_A, size=M)
theta_B = # ___ YOUR CODE
diff_post = theta_B - theta_A
ci_cred = # ___ YOUR CODE
p_B_better = # ___ YOUR CODE
expected_loss = # ___ YOUR CODE
decision = # ___ YOUR CODE
# Build the 5-row comparison DataFrame and print.Goal / LO. Operational discipline: students leave able to flag a sample-ratio mismatch in their own data, refuse to interpret an SRM’d experiment, and (bonus) localise the cause of the mismatch with a stratified diagnostic.
Setup / Task. Three datasets, each intending a 50/50 user split between control and treatment but with observed counts as below:
| Dataset | \(n_A\) (control) | \(n_B\) (treatment) | Intent |
|---|---|---|---|
srm_clean.csv |
5,990 | 6,010 | Clean — small random fluctuation around 50/50 |
srm_buggy.csv |
5,640 | 6,360 | Buggy — 6% drift, classic redirect-failure pattern |
srm_borderline.csv |
5,910 | 6,090 | Borderline — 1.5% drift; stratify before deciding |
Each file also contains a device_type column with three levels (mobile, desktop, tablet). (a) Run a chi-squared test on each dataset against the 50/50 expectation. (b) Apply the Kohavi threshold \(p < 0.001\): for each, decide whether the experiment is safe to interpret. (c) For the buggy dataset, write two sentences hypothesising the cause (bot filtering, redirect failure, opt-out flow). Bonus: stratify the buggy dataset by device_type and identify which device is driving the mismatch.
library(dplyr); library(here); library(readr)
# (a) Build a helper that runs chi-square against the 50/50 expectation and
# returns a 1-row tibble with dataset name, n_A, n_B, p-value, safe-flag.
srm_check <- function(path, name) {
df <- read_csv(here("materials/data", path), show_col_types = FALSE)
counts <- table(df$D)
chi <- # ___ YOUR CODE: chisq.test(counts, p = c(0.5, 0.5))
tibble::tibble(dataset = name,
n_A = counts[["0"]], n_B = counts[["1"]],
chi_p = chi$p.value,
safe = chi$p.value >= 0.001)
}
# Apply to all three datasets and bind into a single table.
# ___ YOUR CODE: bind_rows(srm_check(...) for each file)
# (Bonus) Stratify the buggy dataset by device_type and run chi-square per device.
buggy <- read_csv(here("materials/data/srm_buggy.csv"), show_col_types = FALSE)
# ___ YOUR CODE: group_by(device_type), within-group chi-square,
# identify which device(s) drive the mismatch.import pandas as pd
from scipy.stats import chisquare
def srm_check(path, name):
df = pd.read_csv(f"../../materials/data/{path}")
counts = df["D"].value_counts().sort_index()
n_A, n_B = int(counts.loc[0]), int(counts.loc[1])
chi = # ___ YOUR CODE: chisquare against [n/2, n/2]
return {"dataset": name, "n_A": n_A, "n_B": n_B,
"chi_p": chi.pvalue, "safe": chi.pvalue >= 0.001}
# Apply to the three datasets and assemble.
# ___ YOUR CODE
# Bonus: stratify the buggy dataset by device_type.
buggy = pd.read_csv("../../materials/data/srm_buggy.csv")
# ___ YOUR CODE: per-device chi-square; print or table the offending device.pwr::pwr.t.test() for the MDE / power calculations in fixed-horizon inference. confseq for anytime-valid confidence sequences (betting and mixture variants) implementing 8 and 9. bayestestR for ready-made Bayesian summaries and credible intervals. xgboost for the CUPAC predictor in Exercise 2. randomizr::block_ra() for blocked / stratified random assignment from the design side.statsmodels.stats.power.tt_ind_solve_power for power calculations. scipy.stats.beta and scipy.stats.chisquare for the Bayesian and SRM machinery in Exercises 4 and 5. xgboost for the CUPAC predictor in Exercise 2. statsmodels.formula.api.ols(...).fit(cov_type="HC1") continues from S04 as the HC-robust regression workhorse.