Chapter 40: Experimental Design & A/B Testing

Why Correlation Is Not Enough

The central problem in data-driven decision making is attribution: did our change cause the outcome, or were both the change and the outcome caused by something else? Consider a product team that ships a new onboarding flow on a Tuesday and sees a 5 % lift in day-7 retention. Was it the onboarding flow? Or did a marketing campaign also go live that day, attracting higher-intent users who would have retained anyway? Without experimental control, you cannot tell.

The Rubin potential-outcomes framework makes this precise. For each unit $i$ (a user, a session, a request), define $Y_i(1)$ as the outcome if treated and $Y_i(0)$ as the outcome if not treated. The individual treatment effect is $\tau_i = Y_i(1) - Y_i(0)$. The fundamental problem of causal inference is that you can never observe both potential outcomes for the same unit at the same time — you only ever see one.

The Average Treatment Effect (ATE) is:

In observational data, $\mathbb{E}[Y \mid T=1] - \mathbb{E}[Y \mid T=0]$ is a biased estimator of the ATE whenever the assignment $T$ is correlated with potential outcomes — i.e., whenever there is confounding.

Randomization as a Confounder Killer

Randomization breaks the dependence between treatment assignment and all covariates, both observed and unobserved. When assignment $T$ is independent of $(Y_i(0), Y_i(1))$, the simple difference-in-means becomes an unbiased estimator of the ATE:

This is the magic of the RCT: you do not need to measure or model the confounders. They balance in expectation across arms. This is also why "we tried to control for everything" in observational analysis is never fully satisfying — you can only control for what you measured.

Randomization Units

Choosing the right randomization unit is a critical design decision:

• User-level randomization is the most common. A user always gets the same variant. This prevents cross-contamination and gives a stable user experience.
• Session-level randomization increases statistical power (more units) but is invalid if there is any spillover between sessions of the same user (e.g., a user sees both variants and changes behavior).
• Cluster-level randomization (randomize by city, school, device type) is necessary when treatment has network effects — if treating one user affects another, user-level randomization violates the SUTVA (Stable Unit Treatment Value Assumption) because control units are contaminated by nearby treatment units.

Checking Balance

Even with randomization, always run a balance check (also called an A/A test or covariate check) before interpreting results. Regress treatment assignment on pre-experiment covariates. In a well-randomized experiment, none of them should be predictive:

import pandas as pd
import statsmodels.formula.api as smf

# df has columns: variant (0/1), age, country_us, device_mobile, prior_purchases
model = smf.logit(
    'variant ~ age + country_us + device_mobile + prior_purchases',
    data=df
).fit(disp=False)

print(model.summary2().tables[1][['Coef.', 'P>|z|']])
# All p-values should be large (no covariate predictive of assignment)

A failed balance check signals a Sample Ratio Mismatch or a buggy randomization implementation — both of which invalidate the experiment. We return to SRM in Section 5.

Internal vs External Validity

An experiment has internal validity if its estimate of the ATE for the tested population is unbiased. It has external validity if the ATE generalizes to other populations or contexts. Online A/B tests typically have strong internal validity (randomization is clean) but more limited external validity: results on your current user base may not hold when you scale, enter new markets, or test on a qualitatively different cohort.

The Metric Hierarchy

Before running a single user through an experiment, you need to answer: what are you trying to move, and what must you not break? A mature experimentation culture distinguishes three metric types:

• Primary metric (OEC — Overall Evaluation Criterion): the single number that decides the experiment. Typically a rate or a mean — conversion rate, revenue per user, session length. It must be sensitive enough to detect your expected effect and causally linked to the long-run business goal.
• Secondary metrics: additional signals that help interpret the primary. If conversion rises but average order value falls, you want to know.
• Guardrail metrics: metrics that must not degrade. Page load time, error rate, customer support contact rate. A treatment that wins on conversion but breaks latency is not shippable.

Goodhart's Law applies here: once a metric becomes a target it ceases to be a good measure. Design your OEC to be hard to game and aligned with user value.

Statistical Framework

The classical two-sided hypothesis test for an experiment:

Four key quantities are interrelated — fix any three and the fourth is determined:

• <!--MATHBLOCK5-->: Type I error rate (false positive rate). Typically 0.05.
• <!--MATHBLOCK6-->: Type II error rate (false negative rate). Typically 0.20, giving power <!--MATHBLOCK7-->.
• <!--MATHBLOCK8-->: Minimum Detectable Effect (MDE) — the smallest effect size worth detecting.
• <!--MATHBLOCK9-->: required sample size per variant.

For a two-sample z-test on proportions $p_C$ and $p_T = p_C + \delta$:

For means with standard deviation $\sigma$ (assuming equal variance):

Notice the $\delta^2$ in the denominator: halving your MDE quadruples the required sample size. This is the most important relationship in experiment design.

Choosing the MDE

The MDE is not a statistical choice — it is a business decision. Ask: "What is the smallest effect that would change our product decision?" If an effect smaller than 0.5 % absolute conversion lift is too small to act on (cost of engineering, risk of regression, etc.), then setting MDE = 0.5 % is correct. You are not looking for a 0.1 % effect you would never ship anyway.

A common mistake is setting MDE to whatever the experiment will be powered to detect given a fixed duration. That inverts the logic: MDE drives sample size, not the other way around.

Sample Size Calculation in Python

from scipy.stats import norm
import numpy as np

def sample_size_proportions(
    p_control: float,
    mde: float,          # absolute lift, e.g. 0.02 for +2pp
    alpha: float = 0.05,
    power: float = 0.80,
    two_sided: bool = True,
) -> int:
    """Required sample size per variant for a two-proportion z-test."""
    p_treatment = p_control + mde
    z_alpha = norm.ppf(1 - alpha / (2 if two_sided else 1))
    z_beta  = norm.ppf(power)
    pooled_var = p_control * (1 - p_control) + p_treatment * (1 - p_treatment)
    n = (z_alpha + z_beta)**2 * pooled_var / mde**2
    return int(np.ceil(n))

# Example: checkout conversion rate 5%, want to detect +1pp lift
n = sample_size_proportions(p_control=0.05, mde=0.01)
print(f'Required per variant: {n:,}')
print(f'Total (2 variants):   {2*n:,}')

Duration Estimation

Given required sample size $n$ per variant and a daily traffic count $D$ split evenly:

where $k$ is the number of variants (including control). Practical considerations:

• Run for at least one full business cycle (usually a full week) to avoid day-of-week bias.
• Do not run so long that the experiment population changes substantially (holiday season, major product launch).
• If the computed duration exceeds a few months, revisit your MDE — perhaps you are trying to detect an effect that is not worth detecting.

Ratio Metrics and the Delta Method

Many business metrics are ratios: revenue per user = total revenue / total users, or click-through rate = clicks / impressions when impressions vary per user. The variance of a ratio $\hat{r} = \bar{X} / \bar{Y}$ is not simply $\text{Var}(X)/\text{Var}(Y)$. Use the delta method:

Using the wrong variance formula for ratio metrics is one of the most common sources of inflated false positive rates in online experimentation.

The Analysis Pipeline

Once the experiment has run for the pre-specified duration and collected the pre-specified sample, analysis follows a disciplined pipeline:

1. Check for Sample Ratio Mismatch (SRM) — see Section 5.
2. Compute the test statistic and p-value for the primary metric.
3. Compute and report the confidence interval for the effect, not just the p-value.
4. Check secondary and guardrail metrics.
5. Make a ship/no-ship decision based on the pre-registered criteria.

Never add more data after "peeking" at results and finding borderline p-values. The stopping rule must be fixed before the experiment starts.

Two-Sample Tests for Common Metric Types

Proportions (Conversion Rate, CTR)

For a binary outcome, use the two-proportion z-test. The test statistic under $H_0: p_C = p_T$ is:

where $\hat{p}$ is the pooled conversion rate. In Python:

from statsmodels.stats.proportion import proportions_ztest
import numpy as np

# Control: 1200 conversions out of 24000 visitors (5.0%)
# Treatment: 1380 conversions out of 24000 visitors (5.75%)
count = np.array([1380, 1200])
nobs  = np.array([24000, 24000])

stat, p_value = proportions_ztest(count, nobs)
print(f'z-statistic: {stat:.4f}')
print(f'p-value:     {p_value:.4f}')

# Confidence interval for the absolute difference
from statsmodels.stats.proportion import confint_proportions_2indep
ci = confint_proportions_2indep(1380, 24000, 1200, 24000, method='newcomb')
print(f'95% CI for lift: ({ci[0]*100:.3f}pp, {ci[1]*100:.3f}pp)')

Continuous Metrics (Revenue, Session Length)

For continuous outcomes, use Welch's t-test (does not assume equal variance):

from scipy import stats
import numpy as np

rng = np.random.default_rng(42)
# Simulate skewed revenue data (log-normal)
control   = rng.lognormal(mean=2.0, sigma=1.5, size=10000)
treatment = rng.lognormal(mean=2.1, sigma=1.5, size=10000)  # ~10% lift

t_stat, p_val = stats.ttest_ind(treatment, control, equal_var=False)
diff = treatment.mean() - control.mean()
se   = np.sqrt(treatment.var(ddof=1)/len(treatment) +
               control.var(ddof=1)/len(control))
ci_low, ci_high = diff - 1.96*se, diff + 1.96*se
print(f'Mean diff:   ${diff:.2f}')
print(f'95% CI:      (${ci_low:.2f}, ${ci_high:.2f})')
print(f'p-value:     {p_val:.4f}')

Important: Welch's t-test is robust even for non-normal distributions when sample sizes are large (CLT applies), but for very heavy-tailed distributions (e.g., revenue with whales), consider the Mann-Whitney U test or a bootstrap CI.

Interpreting the Confidence Interval, Not Just the p-value

A p-value below 0.05 only tells you the sign of the effect is unlikely to be noise. The confidence interval tells you the magnitude of the effect and whether it is practically significant. A lift of +0.02 % with a p-value of 0.001 is statistically significant but almost certainly not worth shipping.

Always plot the effect size and its CI. A forest plot across multiple metrics (primary + guardrails) provides an at-a-glance summary. Report relative lift for business audiences and absolute effect for technical audiences.

Multiple Comparisons

If you test $m$ metrics simultaneously at $\alpha = 0.05$, the probability of at least one false positive is $1 - (1-0.05)^m$. For $m=14$ metrics this is nearly 50 %.

Common corrections:

• Bonferroni: use <!--MATHBLOCK8--> per test. Conservative; controls family-wise error rate (FWER).
• Benjamini-Hochberg (BH): controls the False Discovery Rate (FDR) — the expected fraction of discoveries that are false. Less conservative than Bonferroni, preferred when you have many secondary metrics.
• Pre-registration: commit to exactly one primary metric before running. No correction needed for that single test.

from statsmodels.stats.multitest import multipletests
import numpy as np

# 12 secondary metric p-values from an experiment
pvalues = np.array([0.001, 0.043, 0.032, 0.15, 0.62, 0.041,
                    0.29, 0.002, 0.95, 0.08, 0.021, 0.53])

reject_bf, p_adj_bf, _, _ = multipletests(pvalues, method='bonferroni')
reject_bh, p_adj_bh, _, _ = multipletests(pvalues, method='fdr_bh')

for i, (p, rbf, rbh) in enumerate(zip(pvalues, reject_bf, reject_bh)):
    print(f'Metric {i+1:2d}: raw p={p:.3f}  Bonferroni={"sig" if rbf else "ns "}  BH={"sig" if rbh else "ns "}')

The BH procedure is generally recommended for secondary metric analysis in online experiments. For guardrail metrics, consider FWER control (Bonferroni or Holm) since any false negative (missing a regression) is costly.

The Variance Problem

Experiment duration is proportional to the variance of your metric. If you can reduce variance by 50 %, you halve the required sample size — or equivalently, run the same experiment twice as fast. CUPED (Controlled-experiment Using Pre-Experiment Data), introduced by Deng et al. at Microsoft in 2013, is the most practical and widely deployed technique for this.

The insight is that pre-experiment observations of a user's behavior (last week's revenue, last month's sessions) are strong predictors of the user's in-experiment outcome. Including them as covariates absorbs noise without affecting unbiasedness (because they are independent of treatment assignment by construction).

The CUPED Estimator

Let $Y_i$ be the in-experiment outcome and $X_i$ be a pre-experiment covariate (e.g., same metric in a prior period). Define the CUPED-adjusted outcome:

where

This is exactly the coefficient you would get from regressing $Y$ on $X$. The variance reduction is:

where $\rho$ is the Pearson correlation between $Y$ and $X$. If $\rho = 0.7$ (typical for revenue or DAU metrics with one-week lag), variance is cut by $1 - 0.49 = 51\%$ — a factor-of-two speedup.

Crucially, because $X_i$ is independent of treatment assignment, $\mathbb{E}[\tilde{Y}_{\text{treatment}} - \tilde{Y}_{\text{control}}] = \mathbb{E}[Y_{\text{treatment}} - Y_{\text{control}}]$. CUPED does not introduce bias.

Key Implementation Details

Theta is estimated on the full data (both variants combined), not per-variant. Using per-variant theta would absorb part of the treatment effect into the adjustment.

The covariate must be pre-experiment. Using a covariate that could itself be affected by the treatment (e.g., clicks during the experiment) would introduce bias via collider conditioning.

Multiple covariates: you can use multiple pre-experiment features. The optimal linear adjustment is:

where $\boldsymbol{\theta} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{Y}$ (OLS coefficients). This is equivalent to ANCOVA.

import numpy as np
import pandas as pd
from scipy import stats

rng = np.random.default_rng(42)
n = 5000

# Pre-experiment metric (e.g., sessions last week)
X = rng.normal(10, 3, n)

# In-experiment outcome with true treatment effect of 1.0
true_effect = 1.0
variant = rng.integers(0, 2, n)
Y = 5.0 + true_effect * variant + 0.7 * X + rng.normal(0, 2.5, n)

# Naive estimate
t_naive = stats.ttest_ind(Y[variant==1], Y[variant==0])
diff_naive = Y[variant==1].mean() - Y[variant==0].mean()
se_naive = np.sqrt(Y[variant==1].var(ddof=1)/sum(variant==1) +
                   Y[variant==0].var(ddof=1)/sum(variant==0))

# CUPED adjustment (theta from full data)
theta = np.cov(Y, X, ddof=1)[0,1] / np.var(X, ddof=1)
Y_cuped = Y - theta * (X - X.mean())

t_cuped = stats.ttest_ind(Y_cuped[variant==1], Y_cuped[variant==0])
diff_cuped = Y_cuped[variant==1].mean() - Y_cuped[variant==0].mean()
se_cuped = np.sqrt(Y_cuped[variant==1].var(ddof=1)/sum(variant==1) +
                   Y_cuped[variant==0].var(ddof=1)/sum(variant==0))

var_reduction = 1 - Y_cuped.var(ddof=1) / Y.var(ddof=1)

print(f'True effect:              {true_effect:.4f}')
print(f'Naive estimate:           {diff_naive:.4f}  SE={se_naive:.4f}  p={t_naive.pvalue:.4f}')
print(f'CUPED estimate:           {diff_cuped:.4f}  SE={se_cuped:.4f}  p={t_cuped.pvalue:.4f}')
print(f'Variance reduction:       {var_reduction*100:.1f}%')
print(f'Equivalent sample saving: {var_reduction*100:.1f}% fewer users needed')

CUPED vs ANCOVA

CUPED and ANCOVA (Analysis of Covariance) are mathematically equivalent when $\theta$ is estimated by OLS on the combined data. ANCOVA regresses the outcome directly:

The coefficient $\hat{\tau}$ is the treatment effect estimate, and its standard error incorporates the variance reduction automatically. In practice, ANCOVA via statsmodels or pingouin is often cleaner for multi-covariate adjustments.

import statsmodels.formula.api as smf

df = pd.DataFrame({'Y': Y, 'T': variant, 'X': X})
model = smf.ols('Y ~ T + X', data=df).fit()
print(model.summary().tables[1])
# Coefficient on T is identical to CUPED difference-in-means estimate

When CUPED Fails

• New users have no pre-experiment history. CUPED only helps for returning users. Stratified analysis (CUPED for returning, standard test for new) is common.
• Low correlation: if <!--MATHBLOCK19-->, the variance reduction is below 4 % — not worth the implementation complexity.
• Short experiments: if the experiment window is only 1-2 days but you use a 30-day pre-experiment window, the covariate may not be predictive of the narrow in-experiment window.

Peeking and Alpha Inflation

The single most common and damaging mistake in online experimentation is peeking: checking results daily, stopping when p < 0.05, and calling it a win. This is not a minor issue — it catastrophically inflates the false positive rate.

To understand why, consider repeatedly testing as data arrives. If you check at every new sample until $n = 10{,}000$ using a fixed threshold $\alpha = 0.05$, the true probability of ever seeing p < 0.05 under $H_0$ (no effect) can exceed 30 %. You are no longer running a 5 % test.

import numpy as np
from scipy import stats

rng = np.random.default_rng(123)
N_SIMS = 10_000
N_MAX  = 5_000
CHECK_EVERY = 50

fp_peek = 0   # false positive with peeking
fp_fixed = 0  # false positive with fixed horizon

for _ in range(N_SIMS):
    # Under H0: no true effect
    ctrl = rng.normal(0, 1, N_MAX)
    trt  = rng.normal(0, 1, N_MAX)
    ever_significant = False
    for n in range(CHECK_EVERY, N_MAX+1, CHECK_EVERY):
        _, p = stats.ttest_ind(trt[:n], ctrl[:n])
        if p < 0.05:
            ever_significant = True
            break
    if ever_significant:
        fp_peek += 1
    _, p_final = stats.ttest_ind(trt, ctrl)
    if p_final < 0.05:
        fp_fixed += 1

print(f'Fixed-horizon false positive rate: {fp_fixed/N_SIMS:.3f} (target 0.05)')
print(f'Peeking false positive rate:       {fp_peek/N_SIMS:.3f} (inflated!)')

Sequential Testing: Valid Early Stopping

If business requirements demand the ability to stop early, use a sequential testing procedure designed for it. Two modern approaches are widely deployed:

Alpha-Spending Functions (O'Brien-Fleming)

Spend the alpha budget across $K$ pre-planned interim looks. The O'Brien-Fleming boundary is very conservative early and approximately $z_{\alpha/2}$ at the final look:

This gives valid Type I error control with up to $K$ peeks.

E-values and the mSPRT

The mixture Sequential Probability Ratio Test (mSPRT) gives an e-value — a nonnegative random variable with $\mathbb{E}[e] \leq 1$ under $H_0$ — that can be monitored continuously without inflating false positives. You reject when $e \geq 1/\alpha$. E-values compose multiplicatively and are robust to optional stopping:

where $\pi(\theta)$ is a mixing prior over effect sizes.

In practice, Optimizely, Netflix, and many others have adopted variants of mSPRT or Bayesian sequential methods. The key point: you cannot use a fixed-horizon p-value for early stopping without a correction.

Sample Ratio Mismatch (SRM)

A Sample Ratio Mismatch occurs when the observed split between variants differs significantly from the intended split. For a 50/50 experiment receiving $N$ total users, a chi-squared test checks:

from scipy.stats import chisquare

# Expected: 50/50 split over 100,000 users
observed = [51_200, 48_800]  # control, treatment
expected = [50_000, 50_000]

chi2, p = chisquare(f_obs=observed, f_exp=expected)
print(f'chi2={chi2:.2f}, p={p:.6f}')
if p < 0.001:
    print('SRM DETECTED — results should not be trusted.')

SRM is caused by:

• Triggering bugs: users assigned to treatment but not exposed to the experiment (e.g., bot filtering applied only to control).
• Hash collisions: poor randomization hash functions that cluster users.
• Redirect latency: treatment involves a redirect; some users abort before the redirect completes and are not logged.
• Survivorship bias: users who experience errors in the treatment group churn and are excluded from analysis.

Always run SRM checks before interpreting any results. An SRM p-value below 0.001 is grounds for invalidating the experiment.

Novelty and Primacy Effects

Novelty effect: users engage more with a new UI simply because it is new, not because it is better. This inflates short-term metrics. Mitigation: run the experiment long enough (typically 2-4 weeks) for novelty to decay, or analyze by user tenure (new vs returning users respond differently).

Primacy effect (the inverse): users initially resist a change, but adoption grows over time. A short experiment would miss the long-run benefit.

Both effects mean that very short experiments — even if statistically powered — can produce misleading estimates. The holdback experiment (keeping a small percentage on control permanently post-launch) is the gold standard for measuring long-run effects.

Network Effects and SUTVA Violations

In social products, treating one user can affect their connections (network effects). Standard A/B testing assumes SUTVA: the outcome for user $i$ depends only on $i$'s own assignment, not on others'. When this is violated:

• Cluster randomization: randomize at the level of the network cluster (graph community, geographic region) to minimize spillover.
• Ego-network randomization: randomize based on whether the user's social graph is predominantly in treatment or control.
• Switchback experiments (temporal holdouts): alternate treatment and control over time periods in the same market — useful for marketplace experiments (Uber, Lyft) where supply and demand interact.

The Full Experiment Checklist

A disciplined experiment passes through four phases. Skipping any phase undermines the validity of conclusions.

Phase 1: Pre-Experiment Design

• Define the primary metric (OEC) and commit to it. Changing it post-hoc is HARKing (Hypothesizing After Results are Known).
• Define guardrail metrics and the degradation thresholds that trigger a no-ship decision.
• Choose the randomization unit and verify SUTVA plausibility.
• Compute the MDE from a business decision perspective.
• Run a power analysis to determine sample size and duration.
• Decide on the analysis method: fixed-horizon, sequential (mSPRT), or Bayesian.
• Specify the allocation (50/50 for maximum power; unequal splits if risk aversion favors protecting control).
• Document everything in an experiment brief before the experiment starts.

Phase 2: Instrumentation and Launch

• Implement the feature flag with a deterministic, well-tested hash function.
• Log both assignment events (when a user is bucketed) and outcome events (conversion, revenue).
• Run an A/A test (both arms get control) for 24-48 hours to validate that your logging pipeline produces a 50/50 split and baseline metrics match.
• Check SRM on Day 1.

Phase 3: Monitoring (Without Peeking)

• If using a fixed-horizon design, do not analyze the primary metric until the pre-specified end date.
• You may monitor data-quality signals: SRM, logging rate, crash rate, guardrail metrics. These are operational checks, not significance tests of the primary metric.
• If using sequential testing, ensure your monitoring infrastructure applies the correct boundary.

Phase 4: Analysis and Decision

import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.proportion import proportions_ztest, confint_proportions_2indep
from statsmodels.stats.multitest import multipletests

def analyze_experiment(df: pd.DataFrame,
                       primary_metric: str,
                       secondary_metrics: list,
                       guardrail_metrics: list,
                       alpha: float = 0.05):
    """
    Analyze a binary-outcome A/B experiment.
    df must have columns: variant (0=ctrl, 1=trt), + metric columns.
    """
    ctrl = df[df.variant == 0]
    trt  = df[df.variant == 1]

    # 1. SRM check
    from scipy.stats import chisquare
    chi2, srm_p = chisquare([len(ctrl), len(trt)])
    srm_flag = srm_p < 0.001
    print(f'=== SRM Check ===' )
    print(f'  n_ctrl={len(ctrl):,}, n_trt={len(trt):,}')
    print(f'  chi2={chi2:.2f}, p={srm_p:.4f}  {"*** SRM DETECTED ***" if srm_flag else "OK"}')
    if srm_flag:
        print('  WARNING: Results should not be trusted.')
        return

    # 2. Primary metric
    print(f'\n=== Primary Metric: {primary_metric} ===')
    c_conv = ctrl[primary_metric].sum()
    t_conv = trt[primary_metric].sum()
    z, p = proportions_ztest([t_conv, c_conv], [len(trt), len(ctrl)])
    ci = confint_proportions_2indep(t_conv, len(trt), c_conv, len(ctrl))
    lift = trt[primary_metric].mean() - ctrl[primary_metric].mean()
    print(f'  Ctrl rate:  {ctrl[primary_metric].mean()*100:.3f}%')
    print(f'  Trt rate:   {trt[primary_metric].mean()*100:.3f}%')
    print(f'  Lift:       {lift*100:+.3f}pp')
    print(f'  95% CI:     ({ci[0]*100:+.3f}pp, {ci[1]*100:+.3f}pp)')
    print(f'  p-value:    {p:.4f}  {"SIGNIFICANT" if p < alpha else "not significant"}')

    # 3. Guardrail metrics (Bonferroni corrected)
    print(f'\n=== Guardrail Metrics (Bonferroni alpha={alpha/len(guardrail_metrics):.4f}) ===')
    for m in guardrail_metrics:
        _, p_g = stats.ttest_ind(trt[m], ctrl[m])
        adj_alpha = alpha / len(guardrail_metrics)
        print(f'  {m}: diff={trt[m].mean()-ctrl[m].mean():+.4f}  '
              f'p={p_g:.4f}  {"GUARDRAIL VIOLATED" if p_g < adj_alpha and trt[m].mean() < ctrl[m].mean() else "OK"}')

    # 4. Secondary metrics (BH FDR)
    print(f'\n=== Secondary Metrics (BH FDR {alpha}) ===')
    pvals = []
    diffs = []
    for m in secondary_metrics:
        _, p_s = stats.ttest_ind(trt[m], ctrl[m])
        pvals.append(p_s)
        diffs.append(trt[m].mean() - ctrl[m].mean())
    reject, pvals_adj, _, _ = multipletests(pvals, method='fdr_bh')
    for m, d, p_s, p_adj, sig in zip(secondary_metrics, diffs, pvals, pvals_adj, reject):
        print(f'  {m}: diff={d:+.4f}  raw_p={p_s:.4f}  adj_p={p_adj:.4f}  {"sig" if sig else "ns"}')

# --- Demo usage ---
rng = np.random.default_rng(77)
N = 30_000
df_exp = pd.DataFrame({
    'variant': rng.integers(0, 2, N),
})
df_exp['conversion'] = rng.binomial(
    1, np.where(df_exp.variant==1, 0.058, 0.050), N)
df_exp['revenue']    = rng.exponential(
    np.where(df_exp.variant==1, 12, 12), N) * df_exp['conversion']
df_exp['latency_ms'] = rng.normal(
    np.where(df_exp.variant==1, 202, 200), 40, N)
df_exp['support_contact'] = rng.binomial(
    1, 0.02, N)

analyze_experiment(
    df_exp,
    primary_metric='conversion',
    secondary_metrics=['revenue'],
    guardrail_metrics=['latency_ms', 'support_contact'],
)

Making the Ship Decision

A ship decision is not automatic on p < 0.05. Ask:

• Is the CI fully above the MDE threshold (practical significance), or just barely nudging zero?
• Are all guardrail metrics passing?
• Is there a plausible causal story for the effect (mechanism), or does it look like noise in one segment?
• Does the effect hold across major segments (mobile/desktop, new/returning) or is it driven entirely by one noisy slice?

Heterogeneous Treatment Effects (HTE): segment-level analyses (subgroup analysis) are hypothesis-generating, not confirmatory. If an overall non-significant experiment shows a significant effect in one subgroup, you must rerun a targeted experiment for that subgroup before acting.

Bayesian Experiment Analysis

A Bayesian approach replaces p-values with posterior probabilities and credible intervals, which are more intuitively interpretable and handle optional stopping more gracefully. For a conversion rate experiment with a Beta prior:

After observing $c_C$ conversions from $n_C$ users (and similarly for treatment), the posteriors are:

The probability that treatment is better is $P(p_T > p_C | \text{data})$, computed by Monte Carlo:

from scipy.stats import beta
import numpy as np

alpha0, beta0 = 1, 1  # uniform prior
n_ctrl, c_ctrl = 15000, 750   # 5.0% conversion
n_trt,  c_trt  = 15000, 870   # 5.8% conversion

samples = 1_000_000
post_ctrl = beta(alpha0 + c_ctrl, beta0 + n_ctrl - c_ctrl).rvs(samples)
post_trt  = beta(alpha0 + c_trt,  beta0 + n_trt  - c_trt ).rvs(samples)

prob_better = (post_trt > post_ctrl).mean()
expected_lift = (post_trt - post_ctrl).mean()
ci_95 = np.percentile(post_trt - post_ctrl, [2.5, 97.5])

print(f'P(treatment > control): {prob_better:.4f}')
print(f'Expected lift:          {expected_lift*100:.3f}pp')
print(f'95% credible interval:  ({ci_95[0]*100:.3f}pp, {ci_95[1]*100:.3f}pp)')