CUPED and Variance Reduction

25.1 Why the pre-period is a gift (40 min)

Concept. Most of the variation in a customer's revenue during a test has nothing to do with the test: heavy shoppers spend a lot whichever group they're in. That variation is noise for the comparison. But it's predictable: a customer's revenue last month predicts their revenue this month.

CUPED (Controlled-experiment Using Pre-Experiment Data, Deng, Xu, Kohavi and Walker at Microsoft, 2013) subtracts the predictable part:

Y_cuped = Y − θ × (X − mean(X))        θ = cov(X, Y) / var(X)

where X is the same metric measured before the test. Because X is pre-treatment, both groups have the same expected X, so the expected difference between groups is unchanged. The variance falls by a factor of (1 − ρ²), where ρ is the correlation between X and Y. CUPED is the regression of Day 24.3 written as a recipe, and it's now standard at large experimentation platforms.

25.2 CUPED by hand (50 min)

Create day25/cuped_by_hand.py:

# day25/cuped_by_hand.py   (python -m day25.cuped_by_hand)
"""
CUPED (Deng et al., Microsoft, 2013): subtract the part of the outcome that the pre-period predicts.
    Y_cuped = Y - theta * (X - mean(X)),   theta = cov(X, Y) / var(X)
Same expected difference between groups (X is pre-treatment), smaller variance: by a factor (1 - rho^2).
"""
import numpy as np
import polars as pl
import statsmodels.formula.api as smf

from causalkit.experiment_data import email_v2_results

df = email_v2_results()
y, x, t = (df[c].to_numpy().astype(float) for c in ["revenue", "pre_revenue", "treated"])

theta = np.cov(x, y)[0, 1] / np.var(x, ddof=1)
rho = np.corrcoef(x, y)[0, 1]
y_cuped = y - theta * (x - x.mean())


def diff_and_se(v: np.ndarray) -> tuple[float, float]:
    a, b = v[t == 1], v[t == 0]
    return a.mean() - b.mean(), np.sqrt(a.var(ddof=1) / len(a) + b.var(ddof=1) / len(b))


d0, se0 = diff_and_se(y)
d1, se1 = diff_and_se(y_cuped)
print(f"theta = {theta:.3f}   corr(pre_revenue, revenue) = {rho:.3f}   1 - rho^2 = {1 - rho**2:.3f}")
print(f"Plain:  effect {d0:+.3f}  SE {se0:.3f}")
print(f"CUPED:  effect {d1:+.3f}  SE {se1:.3f}   (variance x {(se1 / se0) ** 2:.3f})")

reg = smf.ols("revenue ~ treated + pre_revenue", df.to_pandas()).fit(cov_type="HC2")
print(f"Regression on treated + pre_revenue: {reg.params['treated']:+.3f}  SE {reg.bse['treated']:.3f}"
      "   (CUPED is this regression, done in two steps)")
print(f"Equivalent sample size: CUPED is worth {1 / (1 - rho**2):.2f}x as many customers")

python -m day25.cuped_by_hand

theta = 0.299   corr(pre_revenue, revenue) = 0.493   1 - rho^2 = 0.757
Plain:  effect +2.446  SE 1.011
CUPED:  effect +2.553  SE 0.880   (variance x 0.757)
Regression on treated + pre_revenue: +2.553  SE 0.880   (CUPED is this regression, done in two steps)
Equivalent sample size: CUPED is worth 1.32x as many customers

With a correlation of 0.49, CUPED keeps 76% of the variance: the same precision as running the plain test on 32% more customers. Free, from data the company already has.

25.3 1,000 replays: CUPED and capping (55 min)

Concept. Revenue has a long right tail: a few very large orders add enormous noise. Capping (winsorising) replaces values above a high percentile with that percentile. It reduces variance further, but it changes the estimand: you now measure the effect on capped revenue. Whether that's acceptable is a business choice, made in the plan, before the data arrive.

Create day25/cuped_replays.py:

# day25/cuped_replays.py   (python -m day25.cuped_replays)
"""
1,000 replays with fresh randomisation: plain, CUPED, and CUPED on revenue capped at the 99th percentile.
Capping (winsorising) cuts variance further, but changes WHAT is estimated: be explicit about it.
"""
import numpy as np
import polars as pl

from causalkit.lab import population, run, truth

pop = population()
x = pop["pre_revenue"].to_numpy()
rng = np.random.default_rng(25)


def cuped(y: np.ndarray, x: np.ndarray) -> np.ndarray:
    theta = np.cov(x, y)[0, 1] / np.var(x, ddof=1)
    return y - theta * (x - x.mean())


est = {"plain": [], "CUPED": [], "CUPED + cap at p99": []}
for r in range(1000):
    z = rng.permutation(np.repeat([0, 1], pop.height // 2))
    y = run(pop, z, "email_v2", seed=2000 + r)["revenue"].to_numpy()
    cap = np.percentile(y, 99)                       # pooled cap: same rule for both groups
    for name, v in [("plain", y), ("CUPED", cuped(y, x)), ("CUPED + cap at p99", cuped(np.minimum(y, cap), x))]:
        est[name].append(v[z == 1].mean() - v[z == 0].mean())

target = truth("email_v2")["revenue"]
rows = [{"estimator": k, "mean": round(np.mean(v), 3), "SD": round(np.std(v), 3),
         "SD vs plain": f"{np.std(v) / np.std(est['plain']):.0%}"} for k, v in est.items()]
print(pl.DataFrame(rows))
print(f"True effect on (uncapped) revenue: {target:+.3f}")

python -m day25.cuped_replays

│ estimator          ┆ mean  ┆ SD    ┆ SD vs plain │
│ plain              ┆ 3.494 ┆ 1.002 ┆ 100%        │
│ CUPED              ┆ 3.505 ┆ 0.887 ┆ 88%         │
│ CUPED + cap at p99 ┆ 3.287 ┆ 0.826 ┆ 82%         │
True effect on (uncapped) revenue: +3.480

• CUPED is unbiased and 12% less noisy.
• Capping cuts noise another 6 points, but its average (3.29) is 6% below the true effect on revenue: it estimates the effect on capped revenue, which is smaller because the email also raises some big orders.

Use capping when the decision tolerates it ("we care about typical customers, not whales"), state it in the plan, and apply the same cap to both groups.

25.4 Choosing the covariate (40 min)

The best CUPED covariate is the one most correlated with the outcome, usually the same metric in the pre-period. More covariates can help a little; use regression (Day 24) to combine them.

Try it. For orders, revenue and conversion, which pre-period covariate gives the largest variance reduction: pre_orders or pre_revenue?

# day25/which_covariate.py   (python -m day25.which_covariate)
import numpy as np

from causalkit.experiment_data import email_v2_results

df = email_v2_results()
for metric in ["orders", "revenue", "converted"]:
    y = df[metric].to_numpy().astype(float)
    for cov in ["pre_orders", "pre_revenue"]:
        rho = np.corrcoef(df[cov].to_numpy(), y)[0, 1]
        print(f"{metric:<10} with {cov:<12} rho = {rho:.3f}  variance kept = {1 - rho**2:.2f}")

orders     with pre_orders   rho = 0.447  variance kept = 0.80
orders     with pre_revenue  rho = 0.359  variance kept = 0.87
revenue    with pre_orders   rho = 0.369  variance kept = 0.86
revenue    with pre_revenue  rho = 0.493  variance kept = 0.76
converted  with pre_orders   rho = 0.259  variance kept = 0.93
converted  with pre_revenue  rho = 0.206  variance kept = 0.96

25.5 Day 25 mini-project: can CUPED rescue the guardrail? (70 min)

Goal. Day 22 found that detecting a $0.50 per-customer margin loss would take 110,000 customers. How far does CUPED close the gap?

Create day25/guardrail_plan.py:

# day25/guardrail_plan.py   (python -m day25.guardrail_plan)
"""
Mini-project: Day 22 showed the free-shipping margin guardrail needs ~110,000 customers.
How much does CUPED (with one or two pre-period covariates) bring that down?
"""
import numpy as np
import polars as pl
import statsmodels.formula.api as smf
from scipy.stats import norm

from causalkit.lab import population, run

pop = population()
last_month = run(pop, np.zeros(pop.height), "aa_test", seed=99).join(
    pop.select("customer_id", "pre_revenue", "pre_orders", "loyalty"), on="customer_id").to_pandas()

Z = norm.ppf(0.975) + norm.ppf(0.80)
LOSS = 0.50                                         # the margin loss per customer we must be able to see
rows = []
for label, covs in [("no CUPED", []), ("CUPED: pre_revenue", ["pre_revenue"]),
                    ("CUPED: pre_revenue + pre_orders + loyalty", ["pre_revenue", "pre_orders", "loyalty"])]:
    if covs:
        r2 = smf.ols("margin ~ " + " + ".join(covs), last_month).fit().rsquared
    else:
        r2 = 0.0
    sd = last_month["margin"].std() * np.sqrt(1 - r2)
    n = 2 * Z**2 * sd**2 / LOSS**2
    rows.append({"design": label, "R2 from pre-period": round(r2, 3),
                 "customers per group": int(round(n, -2)), "total": int(round(2 * n, -2))})
with pl.Config(fmt_str_lengths=50, tbl_width_chars=110):
    print(pl.DataFrame(rows))
print("Available: 40,000 customers. Day 30 will need another idea as well: see 'non-inferiority'.")

python -m day25.guardrail_plan

│ design                                    ┆ R2 from pre-period ┆ customers per group ┆ total  │
│ no CUPED                                  ┆ 0.0                ┆ 55400               ┆ 110800 │
│ CUPED: pre_revenue                        ┆ 0.255              ┆ 41300               ┆ 82600  │
│ CUPED: pre_revenue + pre_orders + loyalty ┆ 0.267              ┆ 40600               ┆ 81300  │
Available: 40,000 customers. Day 30 will need another idea as well: see 'non-inferiority'.

CUPED saves about 28,000 customers, but 81,000 is still twice what NorthMart has. The remaining gap is a question of framing, not statistics: a guardrail doesn't need to detect a small loss; it needs to rule out a loss too large to accept. That's a non-inferiority test, which Day 30 uses.

25.6 Save your work with Git (15 min)

git add .
git commit -m "Day 25: CUPED, capping, covariate choice, guardrail planning"
git push

Day 25 checklist

• I can compute CUPED by hand and explain why it doesn't bias the estimate
• I can predict CUPED's variance reduction from the pre-period correlation
• I can explain what capping buys and what it changes about the estimand
• I can choose the best pre-period covariate for a metric